Maybe too involved to get into, but here’s a data generation and fitting example from a recent paper I worked on discussing interventions in dynamic systems with a psychological focus… parallelizing over subjects in the Kalman filter would be a big improvement I guess.
####Individual differences in treatment effects in dynamic system
library(ctsem)
library(rstan)
nlatent=2 #number of latent processes
nmanifest=1 #number of manifest variables
tpoints=30 #number of measurement occasions
ntdpred=1 #number of time dependent predictors
nsubjects=30 #number of subjects
TDPREDMEANS=matrix(0,ntdpred*tpoints,1)
TDPREDMEANS[floor(tpoints/2)]=1 #time dependent predictor (intervention) after 50% of observations
genm=ctModel(Tpoints=tpoints,
n.latent=nlatent, n.manifest=nmanifest, n.TDpred=ntdpred,
LAMBDA=matrix(c(1, 0), nrow=nmanifest, ncol=nlatent),
DRIFT=matrix(c(-.2, 0, 1, -.2), nrow=nlatent, ncol=nlatent),
DIFFUSION=matrix(c(.1, 0, 0, 0.00001), nrow=nlatent, ncol=nlatent),
MANIFESTVAR=matrix(c(.1), nrow=nmanifest, ncol=nmanifest),
TDPREDEFFECT=matrix(c(0, .2), nrow=nlatent, ncol=ntdpred),
CINT=matrix(c(0, 0), nrow=nlatent, ncol=1),
TDPREDMEANS=TDPREDMEANS,
MANIFESTMEANS=matrix(c(0), nrow=nmanifest, ncol=1))
library(plyr)
dat=aaply(1:nsubjects, 1, function(x){ #generate data w random parameter in DRIFT
tempm=genm
stdage=rnorm(1)
tempm$DRIFT[2, 2] = -exp(rnorm(1, -2, .5) + stdage * .5)
cbind(ctGenerate(ctmodelobj=tempm, n.subjects=1, burnin=50, logdtsd=.3), stdage)
})
#convert to long format used by Bayesian ctsem
datlong=ctWideToLong(datawide = dat, Tpoints = tpoints, n.manifest = nmanifest,
n.TDpred = ntdpred,n.TIpred = 1,manifestNames = c("Y1"),
TDpredNames = c("TD1"), TIpredNames=c("stdage"))
datlong=ctDeintervalise(datlong) #convert intervals to abs time
#create and fit model
fitm=ctModel(Tpoints=tpoints, type="stanct", n.latent=nlatent, n.manifest=nmanifest,
n.TDpred=ntdpred, n.TIpred=1, TIpredNames = "stdage",
LAMBDA=matrix(c(1, 0), nrow=nmanifest, ncol=nlatent),
DRIFT=matrix(c("drift11", 0, 1, "drift22"), nrow=nlatent, ncol=nlatent),
DIFFUSION=matrix(c("diffusion11", 0, 0, 0.001), nrow=nlatent, ncol=nlatent),
MANIFESTVAR=matrix(c("merror11"), nrow=nmanifest, ncol=nmanifest),
MANIFESTMEANS=matrix(c("mmeans_Y1"), nrow=nmanifest, ncol=nmanifest),
TDPREDEFFECT=matrix(c(0, "tdpredeffect21"), nrow=nlatent, ncol=ntdpred))
#only the persistence and strength of the predictor effect varies across individuals
fitm$pars$indvarying[-c(8,18)] = FALSE
#and thus standardised age can only affect those parameters
fitm$pars$stdage_effect[-c(8,18)] = FALSE
stanmodel=ctStanFit(datlong, fitm, fit=FALSE)
stanfit = stan(model_code = stanmodel$stanmodeltext, data = stanmodel$data, iter=50, chains=3, cores=3)
summary(fit)
ctsemfit = stanmodel
ctsemfit$stanfit = stanfit
class(ctsemfit) = 'ctStanFit'
ctKalman(ctsemfit, subjects = 1:3, timestep=.01, plot=TRUE,
kalmanvec='etasmooth', errorvec='etasmoothcov', legendcontrol = list(x = "topleft"))