Common and unique process model

gistfile1.txt

##install ctsem (If after May 2019, CRAN is recommended)
install.packages("devtools")
library(devtools)
install_github("cdriveraus/ctsem")
library(ctsem)


#specify generative model (linear sde only at present for data generation -- on the to do list)
gm <- ctModel(
  type='omx', #omx is older model type still needed for data generation
  Tpoints=600, #number of observation occasions
  n.manifest=2, manifestNames=c('y1','y2'), #observed variables
  n.latent=3, #unobserved processes
  LAMBDA=matrix(c( #effect of unobserved processes on observed variables
    1,0,1,
    0,1,1),byrow=TRUE,ncol=3),
  DRIFT=matrix(c( #gradient of unobserved processes conditional on unobserved processes
    -1,0,0,
    0,-1,0,
    0,0,-.2),byrow=TRUE,ncol=3),
  DIFFUSION=matrix(c( # cholesky factor of uncertainty of unobserved processes (although when standard linear estimation is used, off diagonals are parameterized as unconstrained correlation)
    1,0,0,
    0,1,0,
    0,0,1),byrow=TRUE,ncol=3),
  CINT=matrix(0,3), #continuous time intercept of unobserved processes
  T0MEANS=matrix(0,3), #starting point of unobserved processes
  T0VAR=diag(.01,3), #starting cholesky factor uncertainty of unobserved processes
  MANIFESTVAR=diag(0,2), #diagonal matrix with standard deviation of measurement error
  MANIFESTMEANS=matrix(0,2)) #intercept of observed variables


#generate and plot data
d <- ctGenerate(ctmodelobj = gm,n.subjects = 1,burnin = 0,dtmean = .1,logdtsd = 0,wide = FALSE)

matplot(d[,gm$manifestNames],type='l')


#specify model to fit
m <- gm  #use generative model as basis for model to estimate
 
m$DRIFT <- matrix(c( #specify free parameters using character strings -- here just autoeffects
    'dr11',0,0,
    0,'dr22',0,
    0,0,'dr33'),byrow=TRUE,ncol=3)

m$DIFFUSION <- matrix(c( #free standard deviations, but correlation between processes is captured by process 3
    'diff11',0,0,
    0,'diff22',0,
    0,0,'diff33'),byrow=TRUE,ncol=3)

m$T0MEANS=matrix(c('t0m1','t0m2','t0m3'),3)

m$MANIFESTMEANS <- matrix(c('mm1','mm2'),2)

### uncomment following to include nonlinearity if needed (slower to initialise fit due to compilation, and no data generation available at present)
# m$PARS <- matrix(c('diff11base', 'diff11beta',2,1)) #matrix containing extra parameters
# m$DIFFUSION[1,1] = 'exp(PARS[1,1]) + PARS[2,1] * state[3]^2' #std dev of uncertainty of process 1 is now a function of a baseline value and the square of the state of common process.

sm <- ctStanModel(m) #convert omx to stanct type model (for fitting continuous time systems via stan code)
sm$pars #show fixed and free parameters along with transformations used 
sm$pars$indvarying <- FALSE #no random effects across subjects since only 1 subject (this would have been set automatically)
plot(sm) #plot priors

#fit and summary 
f<- ctStanFit(datalong = d,ctstanmodel = sm, 
  # optimize=TRUE, nopriors=FALSE, optimcontrol=list(finishsamples=1000, isloops=5)) #use this line instead for importance sampling
  optimize=TRUE, nopriors=FALSE) #MAP estimate (on unconstrained params) is set here, get maximum likelihood by setting nopriors=TRUE

#f<-ctStanFit(datalong = d,ctstanmodel = sm, iter=400,chains=4) #disable optimization, enable priors, specify iter, chains, cores, for stan hmc fitting)

summary(f)

#kalman filter output
k<-ctStanKalman(f) #extract estimates from kalman filter
names(k) #objects available
dim(k$ypred) #iterations, time points, variables

matplot(apply(k$ypred,c(2,3),mean),type='l',lwd=2,lty=1) #plot mean predictions
matplot(d[,gm$manifestNames],type='p',add=TRUE,pch=c(3)) #overlay observations