cpptsemFromCtsem
cpptsemFromCtsem.Rdtransforms fitted ctsem model to cpptsem model. The implementation of cpptsem closely follows that of ctsemOMX (Driver et al., 2017)
Usage
cpptsemFromCtsem(
ctsemModel,
wideData,
removeD = TRUE,
group = NULL,
groupSpecificParameters = NULL,
silent = FALSE
)Arguments
- ctsemModel
fittet ctsem object
- wideData
Please provide a data set in wide format compatible to ctsemOMX
- removeD
removes the D matrix in the mxObject; This should be set to TRUE
- group
numeric vector indicating which group a person belongs to
- groupSpecificParameters
string vector indicating which parameters should be group specific
- silent
suppress messages
Details
# References
Driver, C. C., Oud, J. H. L., & Voelkle, M. C. (2017). Continuous Time Structural Equation Modelling With R Package ctsem. Journal of Statistical Software, 77(5), 1–36. https://doi.org/10.18637/jss.v077.i05
Examples
if (FALSE) {
library(regCtsem)
addCINT <- FALSE
if(addCINT){
CINT = matrix(c("cint1", "cint2"), nrow = 2, ncol = 1)
}else{
CINT = matrix(c(0,0), nrow = 2, ncol = 1)
}
stationary <- c('T0TRAITEFFECT','T0TIPREDEFFECT')
## ctsem model without trait
AnomAuthmodel1 <- ctModel(LAMBDA = matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2),
Tpoints = 5, n.latent = 2, n.manifest = 2,
MANIFESTVAR=diag(0, 2),
TRAITVAR = NULL,
CINT = CINT)
AnomAuthfit1 <- ctFit(AnomAuth, AnomAuthmodel1, useOptimizer = FALSE, stationary = stationary)
AnomAuthfit1$mxobj$fitfunction$result[[1]]
gradientModel1 <- OpenMx::mxRun(OpenMx::mxModel(AnomAuthfit1$mxobj,
OpenMx::mxComputeSequence(steps=list(
OpenMx::mxComputeNumericDeriv(checkGradient = FALSE,
hessian = FALSE))
)))
centralGrandients <- gradientModel1$compute$steps[[1]]$output[["gradient"]][,"central"]
names(centralGrandients) <- rownames(gradientModel1$compute$steps[[1]]$output[["gradient"]])
centralGrandients
## with cpptsem
cpptsemmodel1 <- cpptsemFromCtsem(ctsemModel = AnomAuthfit1, wideData = AnomAuth)
cpptsemmodel1$computeRAM()
cpptsemmodel1$fitRAM()
cpptsemmodel1$m2LL
cpptsemmodel1$approxRAMGradients((1.1 * 10^(-16))^(1/3))[names(centralGrandients)]
# change parameter values
AnomAuthfit1_1 <- ctFit(AnomAuth, AnomAuthmodel1, useOptimizer = TRUE, stationary = stationary)
newParameters <- omxGetParameters(AnomAuthfit1_1$mxobj)
cpptsemmodel1$setParameterValues(newParameters, names(newParameters))
cpptsemmodel1$computeRAM()
cpptsemmodel1$fitRAM()
cpptsemmodel1$m2LL
# Compute and compare standard errors
ctsemSummary <- summary(AnomAuthfit1_1)
ctsemSummary$ctparameters
computeStandardErrorsDelta(cpptsemObject = cpptsemmodel1,
objective = "ML")
## ctsem model with trait
AnomAuthmodel2 <- ctModel(LAMBDA = matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2),
Tpoints = 5, n.latent = 2, n.manifest = 2,
MANIFESTVAR=diag(0, 2),
TRAITVAR = "auto",
CINT = CINT)
AnomAuthfit2 <- ctFit(AnomAuth, AnomAuthmodel2, useOptimizer = FALSE, stationary = stationary)
AnomAuthfit2$mxobj$fitfunction$result[[1]]
gradientModel2 <- OpenMx::mxRun(OpenMx::mxModel(AnomAuthfit2$mxobj,
OpenMx::mxComputeSequence(steps=list(
OpenMx::mxComputeNumericDeriv(checkGradient = FALSE,
hessian = FALSE))
)))
centralGrandients <- gradientModel2$compute$steps[[1]]$output[["gradient"]][,"central"]
names(centralGrandients) <- rownames(gradientModel2$compute$steps[[1]]$output[["gradient"]])
centralGrandients
## with cpptsem
cpptsemmodel2 <- cpptsemFromCtsem(AnomAuthfit2, wideData = AnomAuth)
cpptsemmodel2$computeRAM()
cpptsemmodel2$fitRAM()
cpptsemmodel2$m2LL
cpptsemmodel2$approxRAMGradients((1.1 * 10^(-16))^(1/3))[names(centralGrandients)]
## Example 3: Kalman Filter
set.seed(175446)
## define the population model:
# set the drift matrix. Note that drift eta_1_eta2 is set to equal 0 in the population.
ct_drift <- matrix(c(-.3,.2,0,-.5), ncol = 2)
generatingModel<-ctsem::ctModel(Tpoints=20,n.latent=2,n.TDpred=0,n.TIpred=0,n.manifest=2,
MANIFESTVAR=diag(0,2),
LAMBDA=diag(1,2),
DRIFT=ct_drift,
DIFFUSION=matrix(c(.5,0,0,.5),2),
CINT=matrix(0,nrow = 2, ncol = 1),
T0MEANS=matrix(0,ncol=1,nrow=2),
T0VAR=diag(1,2), type = "omx")
# simulate a training data and testing data set
traindata <- ctsem::ctGenerate(generatingModel,n.subjects = 20, wide = TRUE)
## Build the analysis model. Note that drift eta1_eta2 is freely estimated
# although it is 0 in the population.
myModel <- ctsem::ctModel(Tpoints=20,n.latent=2,n.TDpred=0,n.TIpred=0,n.manifest=2,
LAMBDA=diag(1,2),
MANIFESTVAR=diag(0,2),
CINT=CINT,
DIFFUSION=matrix(c('eta1_eta1',0,0,'eta2_eta2'),2),
T0MEANS=matrix(0,ncol=1,nrow=2),
T0VAR="auto", type = "omx")
myModel <- ctFit(myModel, dat = traindata, objective = "Kalman", useOptimizer = F, stationary = stationary)
myModel$mxobj$fitfunction$result[[1]]
gradientModel3 <- OpenMx::mxRun(OpenMx::mxModel(myModel$mxobj,
OpenMx::mxComputeSequence(steps=list(
OpenMx::mxComputeNumericDeriv(checkGradient = FALSE,
hessian = FALSE))
)))
centralGrandients <- gradientModel3$compute$steps[[1]]$output[["gradient"]][,"central"]
names(centralGrandients) <- rownames(gradientModel3$compute$steps[[1]]$output[["gradient"]])
centralGrandients
KalmanScores <- mxKalmanScores(myModel$mxobj)
KalmanScores$xUpdated[2:21,]
## with cpptsem
cpptsemmodel3 <- cpptsemFromCtsem(ctsemModel = myModel,wideData = traindata)
cpptsemmodel3$computeAndFitKalman(0) # the 0 tells cpptsem that the fit of the full sample should be computed. Passing 1, 2, ... N to the function will compute the fit for a single sample
cpptsemmodel3$m2LL
cpptsemmodel3$latentScores[1,]
cpptsemmodel3$approxKalmanGradients((1.1 * 10^(-16))^(1/3))[names(centralGrandients)]
## Example 4: Kalman Filter with group or person specific parameter values
set.seed(175446)
## define the population model:
addCINT <- FALSE
if(addCINT){
CINT = matrix(c("cint1", "cint2"), nrow = 2, ncol = 1)
}else{
CINT = matrix(c(0,0), nrow = 2, ncol = 1)
}
stationary <- c('T0TRAITEFFECT','T0TIPREDEFFECT')
# set the drift matrix. Note that drift eta_1_eta2 is set to equal 0 in the population.
ct_drift1 <- matrix(c(-.3,.2,0,-.5), ncol = 2)
ct_drift2 <- matrix(c(-.5,.1,.1,-.2), ncol = 2)
generatingModel1<-ctsem::ctModel(Tpoints=200,n.latent=2,n.TDpred=0,n.TIpred=0,n.manifest=2,
MANIFESTVAR=diag(0,2),
LAMBDA=diag(1,2),
DRIFT=ct_drift1,
DIFFUSION=matrix(c(.5,0,0,.5),2),
CINT=matrix(0,nrow = 2, ncol = 1),
T0MEANS=matrix(0,ncol=1,nrow=2),
T0VAR=diag(1,2), type = "omx")
generatingModel2<-ctsem::ctModel(Tpoints=200,n.latent=2,n.TDpred=0,n.TIpred=0,n.manifest=2,
MANIFESTVAR=diag(0,2),
LAMBDA=diag(1,2),
DRIFT=ct_drift2,
DIFFUSION=matrix(c(.5,0,0,.5),2),
CINT=matrix(0,nrow = 2, ncol = 1),
T0MEANS=matrix(0,ncol=1,nrow=2),
T0VAR=diag(1,2), type = "omx")
traindata1 <- ctsem::ctGenerate(generatingModel1,n.subjects = 10, wide = TRUE)
traindata2 <- ctsem::ctGenerate(generatingModel2,n.subjects = 10, wide = TRUE)
traindata <- rbind(traindata1, traindata2)
## Build the analysis model.
myModel <- ctsem::ctModel(Tpoints=200,n.latent=2,n.TDpred=0,n.TIpred=0,n.manifest=2,
LAMBDA=diag(1,2),
MANIFESTVAR=diag(0,2),
CINT=CINT,
DIFFUSION=matrix(c('eta1_eta1',0,0,'eta2_eta2'),2),
T0MEANS=matrix(0,ncol=1,nrow=2),
T0VAR="auto", type = "omx")
myModel <- ctFit(myModel, dat = traindata, objective = "Kalman", useOptimizer = F, stationary = stationary)
## with cpptsem
cpptsemmodel3 <- cpptsemFromCtsem(ctsemModel = myModel, wideData = traindata, group = c(rep(1,10), rep(2,10)), groupSpecificParameters = c(myModel$mxobj$DRIFT$labels))
startingValues <- cpptsemmodel3$getParameterValues()
m2LLCpptsem <- function(parameters, cpptsemmodel){
cpptsemmodel$setParameterValues(parameters, names(parameters))
# catching all errors from cpptsemmodel
# when parameter values are impossible
invisible(utils::capture.output(KALMAN <- try(cpptsemmodel$computeAndFitKalman(0),
silent = TRUE),
type = "message"))
if(class(KALMAN) == "try-error" | is.na(cpptsemmodel$m2LL)){
return(99999999)
}
return(cpptsemmodel$m2LL)
}
# compute
kalmanCpptsemFit <- Rsolnp::solnp(pars = startingValues,
fun = m2LLCpptsem,
eqfun = NULL, eqB = NULL, ineqfun = NULL, ineqLB = NULL,
ineqUB = NULL, LB = NULL, UB = NULL, control = list(trace = 0),
cpptsemmodel3)
kalmanCpptsemFit$pars
ct_drift1
}