if (FALSE) {
###### Estimation of the primary latent class model ######
# this is a linear latent class mixed model for Ydep1
# with 2 classes and a linear trajectory
set.seed(1234)
PrimMod <- hlme(Ydep1~Time,random=~Time,subject='ID',ng=1,data=data_lcmm)
PrimMod2 <- hlme(Ydep1~Time,mixture=~Time,random=~Time,subject='ID',
ng=2,data=data_lcmm,B=random(PrimMod))
###### Example 1: Relationship between the latent class structure and #
# external class predictors ######
# We consider here 4 external predictors X1-X4.
# estimation of the secondary multinomial logistic model with total variance
# computed with the Hessian
XextHess <- externVar(PrimMod2,
classmb = ~X1 + X2 + X3 + X4,
subject = "ID",
data = data_lcmm,
method = "twoStageJoint")
summary(XextHess)
# estimation of a secondary multinomial logistic model with total variance
# computed with parametric Bootstrap (much longer). When planning to use
# the bootstrap estimator, we recommend running first the analysis
# with option varest = "none" which is faster but which underestimates
# the variance. And then use these values as plausible initial values when
# running the estimation with varest = "paramBoot" to obtain a valid
# variance of the parameters.
XextNone <- externVar(PrimMod2,
classmb = ~X1 + X2 + X3 + X4,
subject = "ID",
data = data_lcmm,
varest = "none",
method = "twoStageJoint")
XextBoot <- externVar(PrimMod2,
classmb = ~X1 + X2 + X3 + X4,
subject = "ID",
data = data_lcmm,
varest = "paramBoot",
method = "twoStageJoint",
B = XextNone$best)
summary(XextBoot)
###### Example 2: Relationship between a latent class structure and #
# external outcome (repeatedly measured over time) ######
# We want to estimate a linear mixed model for Ydep2 with a linear trajectory
# adjusted on X1.
# estimation of the secondary linear mixed model with total variance
# computed with the Hessian
YextHess = externVar(PrimMod2, #primary model
fixed = Ydep2 ~ Time*X1, #secondary model
random = ~Time, #secondary model
mixture = ~Time, #secondary model
subject="ID",
data=data_lcmm,
method = "twoStageJoint")
# estimation of a secondary linear mixed model with total variance
# computed with parametric Bootstrap (much longer). When planning to use
# the bootstrap estimator, we recommend running first the analysis
# with option varest = "none" which is faster but which underestimates
# the variance. And then use these values as plausible initial values when
# running the estimation with varest = "paramBoot" to obtain a valid
# variance of the parameters.
YextNone = externVar(PrimMod2, #primary model
fixed = Ydep2 ~ Time*X1, #secondary model
random = ~Time, #secondary model
mixture = ~Time, #secondary model
subject="ID",
data=data_lcmm,
varest = "none",
method = "twoStageJoint")
YextBoot = externVar(PrimMod2, #primary model
fixed = Ydep2 ~ Time*X1, #secondary model
random = ~Time, #secondary model
mixture = ~Time, #secondary model
subject="ID",
data=data_lcmm,
method = "twoStageJoint",
B = YextNone$best,
varest= "paramBoot")
summary(YextBoot)
###### Example 3: Relationship between a latent class structure and #
# external outcome (survival) ######
# We want to estimate a proportional hazard model (with proportional hazard
# across classes) for time to event Tevent (indicator Event) and assuming
# a splines baseline risk with 3 knots.
# estimation of the secondary survival model with total variance
# computed with the Hessian
YextHess = externVar(PrimMod2, #primary model
survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
hazard="3-quant-splines", #secondary model
hazardtype="PH", #secondary model
subject="ID",
data=data_lcmm,
method = "twoStageJoint")
summary(YextHess)
# estimation of a secondary survival model with total variance
# computed with parametric Bootstrap (much longer). When planning to use
# the bootstrap estimator, we recommend running first the analysis
# with option varest = "none" which is faster but which underestimates
# the variance. And then use these values as plausible initial values when
# running the estimation with varest = "paramBoot" to obtain a valid
# variance of the parameters.
YextNone = externVar(PrimMod2, #primary model
survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
hazard="3-quant-splines", #secondary model
hazardtype="PH", #secondary model
subject="ID",
data=data_lcmm,
varest = "none",
method = "twoStageJoint")
YextBoot = externVar(PrimMod2, #primary model
survival = Surv(Tevent,Event)~ X1+mixture(X2), #secondary model
hazard="3-quant-splines", #secondary model
hazardtype="PH", #secondary model
subject="ID",
data=data_lcmm,
method = "twoStageJoint",
B = YextNone$best,
varest= "paramBoot")
summary(YextBoot)
}
Run the code above in your browser using DataLab