Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Common and unique process model
##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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.