Estimation of extended mixed models using latent classes and latent processes.

Description

Functions for the estimation of latent class mixed models (LCMM), joint latent class mixed models for longitudinal and survival data (JLCM) and latent process mixed models (with or without latent classes of trajectory) for univariate and multivariate longitudinal outcomes of different types including curvilinear and ordinal outcomes. All the models are estimated in a maximum likelihood framework using an iterative algorithm. The package also provides various post fit functions.

Details

The package includes for the moment the estimation of :

• latent class mixed models for Gaussian longitudinal outcomes using hlme function,

• latent class mixed models for other quantitative, bounded quantitative (curvilinear) and discrete (ordinal/binary) longitudinal outcomes using lcmm function,

• mixed models (with and without latent classes) for multivariate longitudinal outcomes of different nature using multlcmm function (this includes a longitudinal IRT model for homogeneous and heterogeneous data),

• joint latent class mixed models for a Gaussian (or curvilinear) longitudinal outcome and a right-censored (potentially left-truncated and of multiple causes) time-to-event using Jointlcmm function,

• joint latent class mixed models for multivariate longitudinal outcomes and a right-censored (potentially left-truncated and of multiple causes) time-to-event using mpjlcmm function.

Please report any bug or comment regarding the package for future updates VIA GITHUB ONLY.

Author(s)

Cecile Proust-Lima, Viviane Philipps, Amadou Diakite and Benoit Liquet

cecile.proust-lima@inserm.fr

References

Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lcmm. Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02

Lin, Turnbull, McCulloch and Slate (2002). Latent class models for joint analysis of longitudinal biomarker and event process data: application to longitudinal prostate-specific antigen readings and prostate cancer. Journal of the American Statistical Association 97, 53-65.

Muthen and Shedden (1999). Finite mixture modeling with mixture outcomes using the EM algorithm. Biometrics 55, 463-9

Proust and Jacqmin-Gadda (2005). Estimation of linear mixed models with a mixture of distribution for the random-effects. Comput Methods Programs Biomed 78:165-73

Proust, Jacqmin-Gadda, Taylor, Ganiayre, and Commenges (2006). A nonlinear model with latent process for cognitive evolution using multivariate longitudinal data. Biometrics 62, 1014-24.

Proust-Lima, Dartigues and Jacqmin-Gadda (2011). Misuse of the linear mixed model when evaluating risk factors of cognitive decline. Amer J Epidemiol 174(9), 1077-88

Proust-Lima and Taylor (2009). Development and validation of a dynamic prognostic tool for prostate cancer recurrence using repeated measures of post-treatment PSA: a joint modelling approach. Biostatistics 10, 535-49.

Proust-Lima, Sene, Taylor, Jacqmin-Gadda (2014). Joint latent class models for longitudinal and time-to-event data: a review. Statistical Methods in Medical Research 23, 74-90.

Proust-Lima, Amieva, Jacqmin-Gadda (2013). Analysis of multivariate mixed longitudinal data: A flexible latent process approach. Br J Math Stat Psychol 66(3), 470-87.

Proust-Lima, Philipps, Perrot, Blanchin, Sebille (2021). Modeling repeated self-reported outcome data: a continuous-time longitudinal Item Response Theory model. arXiv:210913064. http://arxiv.org/abs/2109.13064

Proust-Lima, Dartigues, Jacqmin-Gadda (2016). Joint modeling of repeated multivariate cognitive measures and competing risks of dementia and death: a latent process and latent class approach. Stat Med;35(3):382-98

Proust-Lima, Philipps, Dartigues, Bennett, Glymour, Jacqmin-Gadda, et al (2019). Are latent variable models preferable to composite score approaches when assessing risk factors of change? Evaluation of type-I error and statistical power in longitudinal cognitive studies. Stat Methods Med Res;28(7):1942-57

Verbeke and Lesaffre (1996). A linear mixed-effects model with heterogeneity in the random-effects population. Journal of the American Statistical Association 91, 217-21

See Also

Useful links:

• https://cecileproust-lima.github.io/lcmm/

• Report bugs at https://github.com/CecileProust-Lima/lcmm/issues

Difference of expected prognostic cross-entropy (EPOCE) estimators and its 95% tracking interval between two joint latent class models estimated with Jointlcmm

Description

This function computes the difference of 2 EPOCE estimates (CVPOL or MPOL) and its 95% tracking interval between two joint latent class models estimated using Jointlcmm and evaluated using epoce function. Difference in CVPOL is computed when the EPOCE was previously estimated on the same dataset as used for estimation (using an approximated cross-validation), and difference in MPOL is computed when the EPOCE was previously estimated on an external dataset.

Usage

Diffepoce(epoceM1, epoceM2)

Arguments

Details

This function does not apply for the moment with multiple causes of event (competing risks).

From the EPOCE estimates and the individual contributions to the prognostic observed log-likelihood obtained with epoce function on the same dataset from two different estimated joint latent class models, the difference of CVPOL (or MPOL) and its 95% tracking interval is computed. The 95% tracking interval is:

Delta(MPOL) +/- qnorm(0.975)*sqrt(VARIANCE) for an external dataset

Delta(CVPOL) +/- qnorm(0.975)*sqrt(VARIANCE) for the dataset used in Jointlcmm

where Delta(CVPOL) (or Delta(MPOL)) is the difference of CVPOL (or MPOL) of the two joint latent class models, and VARIANCE is the empirical variance of the difference of individual contributions to the prognostic observed log-likelihoods of the two joint latent class models.

See Commenges et al. (2012) and Proust-Lima et al. (2012) for further details.

Value

Author(s)

Cecile Proust-Lima and Amadou Diakite

References

Commenges, Liquet and Proust-Lima (2012). Choice of prognostic estimators in joint models by estimating differences of expected conditional Kullback-Leibler risks. Biometrics 68(2), 380-7.

Proust-Lima, Sene, Taylor, Jacqmin-Gadda (2014). Joint latent class models for longitudinal and time-to-event data: a review. Statistical Methods in Medical Research 23, 74-90.

See Also

Jointlcmm, epoce, summary.Diffepoce

Examples

## Not run: 
#### estimation with 2 latent classes (ng=2)
m2 <- Jointlcmm(fixed= Ydep1~Time*X1,random=~Time,mixture=~Time,subject='ID'
,survival = Surv(Tevent,Event)~ X1+X2 ,hazard="Weibull"
,hazardtype="PH",ng=2,data=data_lcmm,
B=c( 0.7608, -9.4974,  1.0242,  1.4331,  0.1063 , 0.6714, 10.4679, 11.3178,
 -2.5671, -0.5386,  1.4616, -0.0605,  0.9489,  0.1020,  0.2079,  1.5045),logscale=TRUE)
m1 <- Jointlcmm(fixed= Ydep1~Time*X1,random=~Time,subject='ID'
,survival = Surv(Tevent,Event)~ X1+X2 ,hazard="Weibull"
,hazardtype="PH",ng=1,data=data_lcmm,
B=c(-7.6634,  0.9136,  0.1002,  0.6641, 10.5675, -1.6589,  1.4767, -0.0806,
  0.9240,0.5643,  1.2277,  1.5004))

## EPOCE computation for predictions times from 1 to 6 on the dataset used
## for estimation of m.
VecTime <- c(1,3,5,7,9,11,13,15)
cvpol1 <- epoce(m1,var.time="Time",pred.times=VecTime)
cvpol1
cvpol2 <- epoce(m2,var.time="Time",pred.times=VecTime)
cvpol2
DeltaEPOCE <- Diffepoce(cvpol1,cvpol2)
summary(DeltaEPOCE)
plot(DeltaEPOCE,bty="l")

## End(Not run)

For internal use only ...

Description

For internal use only ...

Conditional probabilities and item information given specified latent process values for lcmm or multlcmm object with ordinal outcomes.

Description

The function computes the conditional probability and information function of each level of each ordinal outcome and the information function at the item level. Confidence bands (and median) can be computed by a Monte Carlo approximation.

Usage

ItemInfo(
  x,
  lprocess,
  condRE_Y = FALSE,
  nsim = 200,
  draws = FALSE,
  ndraws = 2000,
  ...
)

Arguments

Value

An object of class ItemInfo with values :

- ItemInfo: If draws=FALSE, returns a matrix with 3 columns: the first column indicates the name of the outcome, the second indicates the latent process value and the last is the computed Fisher information. If draws=TRUE, returns a matrix with 5 columns: the name of the outcome, the latent process value and the 50%, 2.5% and 97.5% percentiles of the approximated posterior distribution of information.

- LevelInfo: If draws=FALSE, returns a matrix with 5 columns: the first column indicates the name of the outcome, the second indicates the outcome's level, the third indicates the latent process value and the two last contain the probability and Fisher information. If draws=TRUE, returns a matrix with 5 columns: the name of the outcome, the outcome's level, the latent process value and the 50%, 2.5% and 97.5% percentiles of the approximated posterior distribution of the probability and information.

- object: the model from which the computations are done.

- IC: indicator specifying if confidence intervals are computed.

Author(s)

Cecile Proust-Lima, Viviane Philipps

Examples

## Not run: 
## This is a toy example to illustrate the information functions.
## The binary outcomes are arbitrarily created, please do not
## consider them as relevent indicators.
data_lcmm$Yord1 <- as.numeric(data_lcmm$Ydep1>10)
data_lcmm$Yord2 <- as.numeric(data_lcmm$Ydep2>25)
m <- multlcmm(Yord1+Yord2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="thresholds")
info <- ItemInfo(m,lprocess=seq(-4,4,length.out=100),draws=TRUE)
plot(info)
par(mfrow=c(1,2))
plot(info, which="LevelInfo", outcome="Yord1")
plot(info, which="LevelInfo", outcome="Yord2")
plot(info, which="LevelProb", outcome="Yord1")
plot(info, which="LevelProb", outcome="Yord2")

## End(Not run)

Estimation of joint latent class models for longitudinal and time-to-event data

Description

This function fits joint latent class mixed models for a longitudinal outcome and a right-censored (possibly left-truncated) time-to-event. The function handles competing risks and Gaussian or non Gaussian (curvilinear) longitudinal outcomes. For curvilinear longitudinal outcomes, normalizing continuous functions (splines or Beta CDF) can be specified as in lcmm.

Usage

Jointlcmm(
  fixed,
  mixture,
  random,
  subject,
  classmb,
  ng = 1,
  idiag = FALSE,
  nwg = FALSE,
  survival,
  hazard = "Weibull",
  hazardtype = "Specific",
  hazardnodes = NULL,
  hazardrange = NULL,
  TimeDepVar = NULL,
  link = NULL,
  intnodes = NULL,
  epsY = 0.5,
  range = NULL,
  cor = NULL,
  data,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  maxiter = 100,
  nsim = 100,
  prior = NULL,
  pprior = NULL,
  logscale = FALSE,
  subset = NULL,
  na.action = 1,
  posfix = NULL,
  partialH = FALSE,
  verbose = FALSE,
  returndata = FALSE,
  var.time = NULL,
  nproc = 1,
  clustertype = NULL
)

jlcmm(
  fixed,
  mixture,
  random,
  subject,
  classmb,
  ng = 1,
  idiag = FALSE,
  nwg = FALSE,
  survival,
  hazard = "Weibull",
  hazardtype = "Specific",
  hazardnodes = NULL,
  hazardrange = NULL,
  TimeDepVar = NULL,
  link = NULL,
  intnodes = NULL,
  epsY = 0.5,
  range = NULL,
  cor = NULL,
  data,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  maxiter = 100,
  nsim = 100,
  prior = NULL,
  pprior = NULL,
  logscale = FALSE,
  subset = NULL,
  na.action = 1,
  posfix = NULL,
  partialH = FALSE,
  verbose = FALSE,
  returndata = FALSE,
  var.time = NULL,
  nproc = 1,
  clustertype = NULL
)

Arguments

Details

A. BASELINE RISK FUNCTIONS

For the baseline risk functions, the following parameterizations were considered. Be careful, parametrisations changed in lcmm_V1.5:

1. With the "Weibull" function: 2 parameters are necessary w_1 and w_2 so that the baseline risk function a_0(t) = w_1^2 * w_2^2 * (w_1^2 * t)^(w_2^2 - 1) if logscale=FALSE and a_0(t) = exp(w_1) * exp(w_2) * (t)^(exp(w_2) - 1) if logscale=TRUE.

2. with the "piecewise" step function and nz nodes (y_1,...y_nz), nz-1 parameters are necesssary p_1,...p_nz-1 so that the baseline risk function a_0(t) = p_j^2 for y_j < t =< y_j+1 if logscale=FALSE and a_0(t) = exp(p_j) for y_j < t =< y_j+1 if logscale=TRUE.

3. with the "splines" function and nz nodes (y_1,...y_nz), nz+2 parameters are necessary s_1,...s_nz+2 so that the baseline risk function a_0(t) = sum_j s_j^2 M_j(t) if logscale=FALSE and a_0(t) = sum_j exp(s_j) M_j(t) if logscale=TRUE where M_j is the basis of cubic M-splines.

Two parametrizations of the baseline risk function are proposed (logscale=TRUE or FALSE) because in some cases, especially when the instantaneous risks are very close to 0, some convergence problems may appear with one parameterization or the other. As a consequence, we recommend to try the alternative parameterization (changing logscale option) when a joint latent class model does not converge (maximum number of iterations reached) where as convergence criteria based on the parameters and likelihood are small.

B. THE VECTOR OF PARAMETERS B

The parameters in the vector of initial values B or in the vector of maximum likelihood estimates best are included in the following order: (1) ng-1 parameters are required for intercepts in the latent class membership model, and if covariates are included in classmb, ng-1 parameters should be entered for each one; (2) parameters for the baseline risk function: 2 parameters for each Weibull, nz-1 for each piecewise constant risk and nz+2 for each splines risk; this number should be multiplied by ng if specific hazard is specified; otherwise, ng-1 additional proportional effects are expected if PH hazard is specified; otherwise nothing is added if common hazard is specified. In the presence of competing events, the number of parameters should be adapted to the number of causes of event; (3) for all covariates in survival, ng parameters are required if the covariate is inside a mixture(), otherwise 1 parameter is required. Covariates parameters should be included in the same order as in survival. In the presence of cause-specific effects, the number of parameters should be multiplied by the number of causes; (4) for all covariates in fixed, one parameter is required if the covariate is not in mixture, ng parameters are required if the covariate is also in mixture. Parameters should be included in the same order as in fixed; (5) the variance of each random-effect specified in random (including the intercept) if idiag=TRUE and the inferior triangular variance-covariance matrix of all the random-effects if idiag=FALSE; (6) only if nwg=TRUE, ng-1 parameters for class-specific proportional coefficients for the variance covariance matrix of the random-effects; (7) the variance of the residual error.

C. CAUTION

Some caution should be made when using the program:

(1) As the log-likelihood of a latent class model can have multiple maxima, a careful choice of the initial values is crucial for ensuring convergence toward the global maximum. The program can be run without entering the vector of initial values (see point 2). However, we recommend to systematically enter initial values in B and try different sets of initial values.

(2) The automatic choice of initial values that we provide requires the estimation of a preliminary linear mixed model. The user should be aware that first, this preliminary analysis can take time for large datatsets and second, that the generated initial values can be very not likely and even may converge slowly to a local maximum. This is a reason why several alternatives exist. The vector of initial values can be directly specified in B the initial values can be generated (automatically or randomly) from a model with ng=. Finally, function gridsearch performs an automatic grid search.

(3) Convergence criteria are very strict as they are based on derivatives of the log-likelihood in addition to the parameter and log-likelihood stability. In some cases, the program may not converge and reach the maximum number of iterations fixed at 150. In this case, the user should check that parameter estimates at the last iteration are not on the boundaries of the parameter space. If the parameters are on the boundaries of the parameter space, the identifiability of the model is critical. This may happen especially when baseline risk functions involve splines (value close to the lower boundary - 0 with logscale=F -infinity with logscale=F) or classmb parameters that are too high or low (perfect classification) or linkfunction parameters. When identifiability of some parameters is suspected, the program can be run again from the former estimates by fixing the suspected parameters to their value with option posfix. This usually solves the problem. An alternative is to remove the parameters of the Beta of Splines link function from the inverse of the Hessian with option partialH. If not, the program should be run again with other initial values. Some problems of convergence may happen when the instantaneous risks of event are very low and "piecewise" or "splines" baseline risk functions are specified. In this case, changing the parameterization of the baseline risk functions with option logscale is recommended (see paragraph A for details).

Value

The list returned is:

Author(s)

Cecile Proust Lima, Amadou Diakite and Viviane Philipps

cecile.proust-lima@inserm.fr

References

Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lcmm. Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02

Lin, H., Turnbull, B. W., McCulloch, C. E. and Slate, E. H. (2002). Latent class models for joint analysis of longitudinal biomarker and event process data: application to longitudinal prostate-specific antigen readings and prostate cancer. Journal of the American Statistical Association 97, 53-65.

Proust-Lima, C. and Taylor, J. (2009). Development and validation of a dynamic prognostic tool for prostate cancer recurrence using repeated measures of post-treatment PSA: a joint modelling approach. Biostatistics 10, 535-49.

Jacqmin-Gadda, H. and Proust-Lima, C. (2010). Score test for conditional independence between longitudinal outcome and time-to-event given the classes in the joint latent class model. Biometrics 66(1), 11-9

Proust-Lima, Sene, Taylor and Jacqmin-Gadda (2014). Joint latent class models of longitudinal and time-to-event data: a review. Statistical Methods in Medical Research 23, 74-90.

See Also

postprob, plot.Jointlcmm, plot.predict, epoce

Examples

## Not run: 
#### Example of a joint latent class model estimated for a varying number
# of latent classes: 
# The linear mixed model includes a subject- (ID) and class-specific 
# linear trend (intercept and Time in fixed, random and mixture components)
# and a common effect of X1 and its interaction with time over classes 
# (in fixed).
# The variance of the random intercept and slopes are assumed to be equal 
# over classes (nwg=F).
# The covariate X3 predicts the class membership (in classmb). 
# The baseline hazard function is modelled with cubic M-splines -3 
# nodes at the quantiles- (in hazard) and a proportional hazard over 
# classes is assumed (in hazardtype). Covariates X1 and X2 predict the 
# risk of event (in survival) with a common effect over classes for X1
# and a class-specific effect of X2.
# !CAUTION: for illustration, only default initial values where used but 
# other sets of initial values should be tried to ensure convergence
# towards the global maximum.


#### estimation with 1 latent class (ng=1): independent models for the 
# longitudinal outcome and the time of event
m1 <- Jointlcmm(fixed= Ydep1~X1*Time,random=~Time,subject='ID',
survival = Surv(Tevent,Event)~ X1+X2 ,hazard="3-quant-splines",
hazardtype="PH",ng=1,data=data_lcmm)
summary(m1)
#Goodness-of-fit statistics for m1:
#    maximum log-likelihood: -3944.77 ; AIC: 7919.54  ;  BIC: 7975.09  

## End(Not run)

#### estimation with 2 latent classes (ng=2)
m2 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",ng=2,data=data_lcmm,
B=c(0.64,-0.62,0,0,0.52,0.81,0.41,0.78,0.1,0.77,-0.05,10.43,11.3,-2.6,
-0.52,1.41,-0.05,0.91,0.05,0.21,1.5))
summary(m2)
#Goodness-of-fit statistics for m2:
#       maximum log-likelihood: -3921.27; AIC: 7884.54; BIC: 7962.32  

## Not run: 
#### estimation with 3 latent classes (ng=3)
m3 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~ X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",ng=3,data=data_lcmm,
B=c(0.77,0.4,-0.82,-0.27,0,0,0,0.3,0.62,2.62,5.31,-0.03,1.36,0.82,
-13.5,10.17,10.24,11.51,-2.62,-0.43,-0.61,1.47,-0.04,0.85,0.04,0.26,1.5))
summary(m3)
#Goodness-of-fit statistics for m3:
#       maximum log-likelihood: -3890.26 ; AIC: 7834.53;  BIC: 7934.53  

#### estimation with 4 latent classes (ng=4)
m4 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~ X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",ng=4,data=data_lcmm,
B=c(0.54,-0.42,0.36,-0.94,-0.64,-0.28,0,0,0,0.34,0.59,2.6,2.56,5.26,
-0.1,1.27,1.34,0.7,-5.72,10.54,9.02,10.2,11.58,-2.47,-2.78,-0.28,-0.57,
1.48,-0.06,0.61,-0.07,0.31,1.5))
summary(m4)
#Goodness-of-fit statistics for m4:
#   maximum log-likelihood: -3886.93 ; AIC: 7839.86;  BIC: 7962.09  


##### The model with 3 latent classes is retained according to the BIC  
##### and the conditional independence assumption is not rejected at
##### the 5% level. 
# posterior classification
plot(m3,which="postprob")
# Class-specific predicted baseline risk & survival functions in the 
# 3-class model retained (for the reference value of the covariates) 
plot(m3,which="baselinerisk",bty="l")
plot(m3,which="baselinerisk",ylim=c(0,5),bty="l")
plot(m3,which="survival",bty="l")
# class-specific predicted trajectories in the 3-class model retained 
# (with characteristics of subject ID=193)
data <- data_lcmm[data_lcmm$ID==193,]
plot(predictY(m3,var.time="Time",newdata=data,bty="l"))
# predictive accuracy of the model evaluated with EPOCE
vect <- 1:15
cvpl <- epoce(m3,var.time="Time",pred.times=vect)
summary(cvpl)
plot(cvpl,bty="l",ylim=c(0,2))
############## end of example ##############

## End(Not run)

Standard methods for estimated models

Description

coef and vcov for hlme, lcmm, mutlcmm, Jointlcmm, mpjlcmm, externSurv, and externX models, fixef, ranef, fitted and residuals methods for estimated hlme, lcmm, mutlcmm and Jointlcmm models.

Usage

## S3 method for class 'hlme'
coef(object, ...)

Arguments

Value

For coef, the vector of the estimates.

For vcov, the variance-covariance matrix of the estimates.

For fixef : - for hlme, lcmm and multlcmm objects, a list containing the fixed effects estimates in the class-membership model and in the longitudinal model. - for Jointlcmm objects, a list containing the fixed effects estimates in the class-membership model, the survival model and in the longitudinal model.

For ranef, a matrix (nrow=number of subjects, ncol=number of covariates with random effect) containing the individual random effects.

For fitted, a vector containing the subject-specific predictions extracted from object.

For residuals, a vector containing the subject-specific residuals extracted from object.

Author(s)

Cecile Proust-Lima, Viviane Philipps

Variance-covariance of the estimates

Description

This function provides the variance-covariance matrix of the estimates. vcov is an alias for it.

Usage

VarCov(x)

Arguments

Value

a matrix containing the variance-covariance of the estimates. For the parameters of the matrix of variance-covariance of the random effects, the Cholesky transformed parameters are considered so that VarCov provides the covariance matrix of function estimates with cholesky=TRUE.

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

estimates

Estimates, standard errors and Wald test for the parameters of the variance-covariance matrix of the random effects.

Description

Fromm the Cholesky transformed parameters, this function provides estimates, standard errors and Wald test for the parameters of the variance-covariance matrix of the random effects.

Usage

VarCovRE(Mod)

Arguments

Value

a matrix containing the estimates of the parameters of the variance-covariance matrix of the random effects, their standard errors, and, for the covariance parameters, the Wald statistic and the associated p-value.

Author(s)

Cecile Proust-Lima, Lionelle Nkam and Viviane Philipps

Percentage of variance explained by the (latent class) linear mixed model regression

Description

The function provides the percentage of variance explained by the (latent class) linear mixed regression in a model estimated with hlme, lcmm, multlcmm or Jointlcmm.

Usage

VarExpl(x, values)

Arguments

Value

For hlme, lcmm, and Jointlcmm objects, the function returns a matrix with 1 row and ng (ie the number of latent classes) columns containing (the class specific) percentages of variance explained by the linear mixed regression.

For multlcmm objects, the function returns a matrix containing (the class specific) percentages of variance explained by the linear mixed regression for each outcome. The resulting matrix is composed of as many rows as outcomes and as many columns as latent classes.

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

hlme, lcmm, multlcmm, Jointlcmm

Examples

## Not run: 
m1 <- multlcmm(Ydep1+Ydep2~1+Time*X2+contrast(X2),random=~1+Time,
subject="ID",randomY=TRUE,link=c("4-manual-splines","3-manual-splines"),
intnodes=c(8,12,25),data=data_lcmm, 
B=c(-1.071, -0.192,  0.106, -0.005, -0.193,  1.012,  0.870,  0.881,
  0.000,  0.000, -7.520,  1.401,  1.607 , 1.908,  1.431,  1.082,
 -7.528,  1.135 , 1.454 , 2.328, 1.052))

# variation percentages explained by linear mixed regression
VarExpl(m1,data.frame(Time=0))

## End(Not run)

Multivariate Wald Test

Description

This function provides multivariate and univariate Wald tests for combinations of parameters from hlme, lcmm, multlcmm, Jointlcmm or mpjlcmm models.

Usage

WaldMult(Mod, pos = NULL, contrasts = NULL, name = NULL, value = NULL)

Arguments

Value

If contrasts is NULL, the function returns a matrix with 1 row and 2 columns containing the value of the Wald test's statistic and the associated p-value.

If contrasts is not NULL, the function returns a matrix with 1 row and 4 columns containing the value of the coefficient (dot product of pos and contrasts), his standard deviation, the value of the Wald test's statistic and the associated p-value.

Author(s)

Cecile Proust-Lima, Lionelle Nkam and Viviane Philipps

Predicted cumulative incidence of event according to a profile of covariates

Description

This function computes the predicted cumulative incidence of each cause of event according to a profile of covariates from a joint latent class model. Confidence bands can be computed by a Monte-Carlo method.

Usage

cuminc(x, time, draws = FALSE, ndraws = 2000, integrateOptions = NULL, ...)

Arguments

Value

An object of class cuminc containing as many matrices as profiles defined by the covariates values. Each of these matrices contains the event-specific cumulative incidences in each latent class at the different times specified.

Author(s)

Viviane Philipps and Cecile Proust-Lima

See Also

Jointlcmm, plot.Jointlcmm, plot.cuminc

Examples

m2 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",ng=2,data=data_lcmm,
B=c(0.64,-0.62,0,0,0.52,0.81,0.41,0.78,0.1,0.77,-0.05,10.43,11.3,-2.6,
-0.52,1.41,-0.05,0.91,0.05,0.21,1.5))

par(mfrow=c(1,2))
plot(cuminc(m2,time=seq(0,20),X1=0,X2=0), ylim=c(0,1))
plot(cuminc(m2,time=seq(0,20),X1=0,X2=1), ylim=c(0,1))

Simulated dataset for hlme function

Description

The data were simulated from a 3-latent class linear mixed model. Repeated data for 100 subjects were simulated. The three latent classes are predicted by X2 and X3. In each latent class, Y follows a linear mixed model including intercept and time both with correlated random-effects and class-specific fixed effects. In addition, X1 and X1*time have a common impact over classes on the Y trajectory.

Format

A data frame with 326 observations on the following 9 variables.

ID

subject identification number

Y

longitudinal outcome

Time

time of measurement

X1

binary covariate

X2

binary covariate

X3

binary covariate

See Also

hlme, postprob, summary.lcmm, plot.predict

Simulated dataset for lcmm and Jointlcmm functions

Description

The data were simulated from a joint latent class mixed model with 3 classes. Repeated data of 3 longitudinal outcomes (Ydep1, Ydep2, Ydep3) and censored time of event (Tevent, Event) with delayed entry (Tentry) were simulated for a total of 300 subjects. The three latent classes were predicted by the continuous covariate X3. In each latent class, the longitudinal outcome Ydep1 followed a linear mixed model including intercept, time and squared time both with correlated random-effects and class-specific fixed effects. In addition, the binary covariate X1 and its interaction with time X1:Time had a common impact (over classes) on the Ydep1 trajectory. The longitudinal ordinal outcomes Ydep2 and Ydep3 were generated from Ydep1 using threshold models with respectively 30 and 10 thresholds. In each latent class, the time of event followed a class-specific Weibull hazard with a common proportional effect of the binary covariate X2. Both time of entry Tentry and time of censoring had a uniform distribution

Format

A data frame with 1678 observations over 300 different subjects and 22 variables.

ID

subject identification number

Ydep1

longitudinal continuous outcome

Ydep2

longitudinal ordinal outcome with 31 levels

Ydep3

longitudinal ordinal outcome with 11 levels

Tentry

delayed entry for the time-to-event

Tevent

observed time-to-event: either censoring time or time of event

Event

indicator that Tevent is the time of event

Time

time of measurement

X1

binary covariate

X2

binary covariate

X3

continuous covariate

X4

categorical covariate

See Also

Jointlcmm, lcmm, hlme

Individual dynamic predictions from a joint latent class model

Description

This function computes individual dynamic predictions and 95% confidence bands. Given a joint latent class model, a landmark time s, a horizon time t and measurements until time s, the predicted probability of event in the window [s,s+t] is calculated. Confidence bands can be provided using a Monte Carlo method.

Usage

dynpred(
  model,
  newdata,
  event = 1,
  landmark,
  horizon,
  var.time,
  fun.time = identity,
  na.action = 1,
  draws = FALSE,
  ndraws = 2000
)

Arguments

Value

A list containing :

Author(s)

Cecile Proust-Lima, Viviane Philipps

References

Proust-Lima, Sene, Taylor and Jacqmin-Gadda (2014). Joint latent class models of longitudinal and time-to-event data: a review. Statistical Methods in Medical Research 23, 74-90.

See Also

plot.dynpred, Jointlcmm, predictY, plot.predict

Examples

## Joint latent class model with 2 classes :
m32 <- Jointlcmm(Ydep1~Time*X1,mixture=~Time,random=~Time,subject="ID",
classmb=~X3,ng=2,survival=Surv(Tevent,Event)~X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",data=data_lcmm,
B = c(0.641, -0.6217, 0, 0, 0.5045, 0.8115, -0.4316, 0.7798, 0.1027, 
0.7704, -0.0479, 10.4257, 11.2972, -2.5955, -0.5234, 1.4147, 
-0.05, 0.9124, 0.0501, 0.2138, 1.5027))

## Predictions at landmark 10 and 12 for horizon 3, 5 and 10 for two subjects :

dynpred(m32,landmark=c(10,12),horizon=c(3,5,10),var.time="Time",
fun.time=function(x){10*x},newdata=data_lcmm[1:8,])
## Not run: 
dynpred(m32,landmark=c(10,12),horizon=c(3,5,10),var.time="Time",
fun.time=function(x){10*x},newdata=data_lcmm[1:8,],draws=TRUE,ndraws=2000)

## End(Not run)

Estimators of the Expected Prognostic Observed Cross-Entropy (EPOCE) for evaluating predictive accuracy of joint latent class models estimated using Jointlcmm

Description

This function computes estimators of the Expected Prognostic Observed Cross-Entropy (EPOCE) for evaluating the predictive accuracy of joint latent class models estimated using Jointlcmm. On the same data as used for estimation of the Jointlcmm object, this function computes both the Mean Prognostic Observed Log-Likelihood (MPOL) and the Cross-Validated Observed Log-Likelihood (CVPOL), two estimators of EPOCE. The latter corrects the MPOL estimate for over-optimism by approximated cross-validation. On external data, this function only computes the Mean Prognostic Observed Log-Likelihood (MPOL).

Usage

epoce(
  model,
  pred.times,
  var.time,
  fun.time = identity,
  newdata = NULL,
  subset = NULL,
  na.action = 1
)

Arguments

Details

This function does not apply for the moment with multiple causes of event (competing risks).

EPOCE assesses the prognostic information of a joint latent class model. It relies on information theory.

MPOL computed at time s equals minus the mean individual contribution to the conditional log-likelihood of the time to event given the longitudinal data up to the time of prediction s and given the subject is still at risk of event in s.

CVPOL computed at time s equals MPOL at time s plus a penalty term that corrects for over-optimism when computing predictive accuracy measures on the same dataset as used for estimation. This penalty term is computed from the inverse of the Hessian of the joint log-likelihood and the product of the gradients of the contributions to respectively the joint log-likelihood and the conditional log-likelihood.

The theory of EPOCE and its estimators MPOL and CVPOL is given in Commenges et al. (2012), and further detailed and illustrated for joint models in Proust-Lima et al. (2013).

Value

Author(s)

Cecile Proust-Lima and Amadou Diakite

References

Commenges, Liquet and Proust-Lima (2012). Choice of prognostic estimators in joint models by estimating differences of expected conditional Kullback-Leibler risks. Biometrics 68(2), 380-7.

Proust-Lima, Sene, Taylor and Jacqmin-Gadda (2014). Joint latent class models of longitudinal and time-to-event data: a review. Statistical Methods in Medical Research 23, 74-90.

See Also

Jointlcmm, print.epoce, summary.epoce, plot.epoce

Examples

## Not run: 
## estimation of a joint latent class model with 2 latent classes (ng=2)
# (see the example section of Jointlcmm for details about
#  the model specification)

m <- Jointlcmm(fixed= Ydep1~Time*X1,random=~Time,mixture=~Time,subject='ID'
,survival = Surv(Tevent,Event)~ X1+X2 ,hazard="Weibull"
,hazardtype="PH",ng=2,data=data_lcmm,logscale=TRUE,
B=c(0.7608, -9.4974 , 1.0242,  1.4331 , 0.1063 , 0.6714, 10.4679, 11.3178,
 -2.5671, -0.5386,  1.4616, -0.0605,  0.9489,  0.1020 , 0.2079,  1.5045))
summary(m)

## Computation of the EPOCE on the same dataset as used for
# estimation of m with times at predictions from 1 to 15 
VecTime <- c(1,3,5,7,9,11,13,15)
cvpl <- epoce(m,var.time="Time",pred.times=VecTime)
summary(cvpl)
plot(cvpl,bty="l",ylim=c(0,2))

## End(Not run)

Maximum likelihood estimates

Description

This function provides the vector of maximum likelihood estimates of a model estimated with hlme, lcmm, multlcmm, Jointlcmm, mpjlcmm, externSurv, or externX.

Usage

estimates(x, cholesky = TRUE)

Arguments

Value

a vector with all estimates of the model.

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

VarCov, hlme, lcmm, multlcmm, Jointlcmm

Estimation of a secondary regression model after the estimation of a primary latent class model

Description

This function fits regression models to relate a latent class structure (stemmed from a latent class model estimated within lcmm package) with either an external outcome or external class predictors. Two inference techniques are implemented. They both account for the classification error in the posterior class assignment:

- a 2-stage estimation using the joint likelihood of the primary latent class model and of the secondary/ external regression;

- a conditional regression of the external outcome given the underlying latent class structure, or of the underlying class structure given external covariates.

It returns an object of one of the lcmm package classes.

Usage

externVar(
  model,
  fixed,
  mixture,
  random,
  subject,
  classmb,
  survival,
  hazard = "Weibull",
  hazardtype = "Specific",
  hazardnodes = NULL,
  TimeDepVar = NULL,
  logscale = FALSE,
  idiag = FALSE,
  nwg = FALSE,
  randomY = NULL,
  link = NULL,
  intnodes = NULL,
  epsY = NULL,
  cor = NULL,
  nsim = NULL,
  range = NULL,
  data,
  longitudinal,
  method,
  varest,
  M = 200,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  maxiter = 100,
  posfix,
  partialH = FALSE,
  verbose = FALSE,
  nproc = 1
)

Arguments

Details

A. DATA STRUCTURE

The data argument must follow specific structure. It must include all the data necessary to compute the posterior classification probabilities (so a longitudinal format usually) as well as the information for the secondary analysis. For time-invariant variables in the secondary analyses: - if used as an external outcome: the information should not be duplicated at each row of the subject. It should appear once for each individual. - if used as an external covariate: the information can be duplicated at each row of the subject (as usual)

B. VARIANCE ESTIMATION

The two techniques rely on a sequential analysis (two-stage analysis) so the variance calculation should account for both the uncertainty in the first and the second stage. Not taking into account the first-stage uncertainty by specifying varest="none" may lead to the underestimation of the final variance. When possible, Method varest="Hessian" which relies on the combination of Hessians from the primary and secondary models is recommended. However, it may become numerically intensive when the primary latent class model includes a high number of parameters. As an alternative, especially when the primary model is complex and the second model includes a limited number of parameters, the parametric Bootstrap method varest="paramBoot" can be favored.

Value

an object of class externVar and externSurv for external survival outcomes, externX for external class predictors, and hlme, lcmm, or multlcmm for external longitudinal or cross-sectional outcomes.

Author(s)

Maris Dussartre, Cecile Proust-Lima and Viviane Philipps

Examples

## Not run: 


###### 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)


## End(Not run)

Marginal predictions of the longitudinal outcome(s) in their natural scale from lcmm, Jointlcmm or multlcmm objects

Description

The function computes the marginal predictions of the longitudinal outcome(s) in their natural scale on the individual data used for the estimation from lcmm, Jointlcmm or multlcmm objects.

Usage

fitY(x)

Arguments

Value

For lcmm and Jointlcmm objects, returns a matrix with ng+1 columns containing the subject identifier and the ng class-specific marginal predicted values.

For multlcmm objects, returns a matrix with ng+2 columns containing the subject identifier, the outcome indicator and the ng class-specific predicted values.

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

predictY, plot.lcmm

Automatic grid search

Description

This function provides an automatic grid search for latent class mixed models estimated with hlme, lcmm, multlcmm and Jointlcmm functions.

Usage

gridsearch(m, rep, maxiter, minit, cl = NULL)

Arguments

Details

The function permits the estimation of a model from a grid of random initial values to reduce the odds of a convergence towards a local maximum.

The function was inspired by the emEM technique described in Biernacki et al. (2003). It consists in:

1. randomly generating rep sets of initial values for m from the estimates of minit (this is done internally using option B=random(minit) rep times)

2. running the optimization algorithm for the model specified in m from the rep sets of initial values with a maximum number of iterations of maxit each time.

3. retaining the estimates of the random initialization that provides the best log-likelihood after maxiter iterations.

4. running the optimization algorithm from these estimates for the final estimation.

Value

an object of class hlme, lcmm, multlcmm, Jointlcmm or mpjlcmm corresponding to the call specified in m.

Author(s)

Cecile Proust-Lima and Viviane Philipps

References

Biernacki C, Celeux G, Govaert G (2003). Choosing Starting Values for the EM Algorithm for Getting the Highest Likelihood in Multivariate Gaussian Mixture models. Computational Statistics and Data Analysis, 41(3-4), 561-575.

Examples

## Not run: 
# initial model with ng=1 for the random initial values
m1 <- hlme(Y ~ Time * X1, random =~ Time, subject = 'ID', ng = 1, 
      data = data_hlme)

# gridsearch with 10 iterations from 50 random departures
m2d <- gridsearch(rep = 50, maxiter = 10, minit = m1, hlme(Y ~ Time * X1,
      mixture =~ Time, random =~ Time, classmb =~ X2 + X3, subject = 'ID',
          ng = 2, data = data_hlme))
        
## End(Not run)

Estimation of latent class linear mixed models

Description

This function fits linear mixed models and latent class linear mixed models (LCLMM) also known as growth mixture models or heterogeneous linear mixed models. The LCLMM consists in assuming that the population is divided in a finite number of latent classes. Each latent class is characterised by a specific trajectory modelled by a class-specific linear mixed model. Both the latent class membership and the trajectory can be explained according to covariates. This function is limited to a mixture of Gaussian outcomes. For other types of outcomes, please see function lcmm. For multivariate longitudinal outcomes, please see multlcmm.

Usage

hlme(
  fixed,
  mixture,
  random,
  subject,
  classmb,
  ng = 1,
  idiag = FALSE,
  nwg = FALSE,
  cor = NULL,
  data,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  prior,
  pprior = NULL,
  maxiter = 500,
  subset = NULL,
  na.action = 1,
  posfix = NULL,
  verbose = FALSE,
  returndata = FALSE,
  var.time = NULL,
  partialH = FALSE,
  nproc = 1,
  clustertype = NULL
)

Arguments

Details

A. THE VECTOR OF PARAMETERS B

The parameters in the vector of initial values B or equivalently in the vector of maximum likelihood estimates best are included in the following order:

(1) ng-1 parameters are required for intercepts in the latent class membership model, and when covariates are included in classmb, ng-1 paramaters should be entered for each covariate;

(2) for all covariates in fixed, one parameter is required if the covariate is not in mixture, ng paramaters are required if the covariate is also in mixture;

(3) the variance of each random-effect specified in random (including the intercept) when idiag=TRUE, or the inferior triangular variance-covariance matrix of all the random-effects when idiag=FALSE;

(4) only when nwg=TRUE, ng-1 parameters are required for the ng-1 class-specific proportional coefficients in the variance covariance matrix of the random-effects;

(5) when cor is specified, 1 parameter corresponding to the variance of the Brownian motion should be entered with cor=BM and 2 parameters corresponding to the correlation and the variance parameters of the autoregressive process should be entered

(6) the standard error of the residual error.

B. CAUTIONS

Some caution should be made when using the program:

(1) As the log-likelihood of a latent class model can have multiple maxima, a careful choice of the initial values is crucial for ensuring convergence toward the global maximum. The program can be run without entering the vector of initial values (see point 2). However, we recommend to systematically enter initial values in B and try different sets of initial values.

(2) The automatic choice of initial values we provide requires the estimation of a preliminary linear mixed model. The user should be aware that first, this preliminary analysis can take time for large datatsets and second, that the generated initial values can be very not likely and even may converge slowly to a local maximum. This is the reason why several alternatives exist. The vector of initial values can be directly specified in B the initial values can be generated (automatically or randomly) from a model with ng=. Finally, function gridsearch performs an automatic grid search.

(3) Convergence criteria are very strict as they are based on the derivatives of the log-likelihood in addition to the parameter stability and log-likelihood stability. In some cases, the program may not converge and reach the maximum number of iterations fixed at 100. In this case, the user should check that parameter estimates at the last iteration are not on the boundaries of the parameter space. If the parameters are on the boundaries of the parameter space, the identifiability of the model is critical. This may happen especially with splines parameters that may be too close to 0 (lower boundary) or classmb parameters that are too high or low (perfect classification). When identifiability of some parameters is suspected, the program can be run again from the former estimates by fixing the suspected parameters to their value with option posfix. This usually solves the problem. An alternative is to remove the parameters of the Beta of Splines link function from the inverse of the Hessian with option partialH. If not, the program should be run again with other initial values, with a higher maximum number of iterations or less strict convergence tolerances.

Value

The list returned is:

Author(s)

Cecile Proust-Lima, Benoit Liquet and Viviane Philipps

cecile.proust-lima@inserm.fr

References

Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lcmm. Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02

Verbeke G and Lesaffre E (1996). A linear mixed-effects model with heterogeneity in the random-effects population. Journal of the American Statistical Association 91, 217-21

Muthen B and Shedden K (1999). Finite mixture modeling with mixture outcomes using the EM algorithm. Biometrics 55, 463-9

Proust C and Jacqmin-Gadda H (2005). Estimation of linear mixed models with a mixture of distribution for the random-effects. Computer Methods Programs Biomedicine 78, 165-73

See Also

postprob, plot.hlme, summary, predictY

Examples

##### Example of a latent class model estimated for a varying number
# of latent classes: 
# The model includes a subject- (ID) and class-specific linear 
# trend (intercept and Time in fixed, random and mixture components)
# and a common effect of X1 and its interaction with time over classes 
# (in fixed). 
# The variance of the random intercept and slope are assumed to be equal 
# over classes (nwg=F).
# The covariate X3 predicts the class membership (in classmb).
#
# !CAUTION: initialization of mixed models with latent classes is 
# of most importance because of the problem of multimodality of the likelihood.
# Calls m2a-m2d illustrate the different implementations for the 
# initial values.

### homogeneous linear mixed model (standard linear mixed model) 
### with correlated random-effects
m1<-hlme(Y~Time*X1,random=~Time,subject='ID',ng=1,data=data_hlme)
summary(m1)

### latent class linear mixed model with 2 classes

# a. automatic specification from G=1 model estimates:
m2a<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
         ng=2,data=data_hlme,B=m1)
         
# b. vector of initial values provided by the user:
m2b<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
         ng=2,data=data_hlme,B=c(0.11,-0.74,-0.07,20.71,
                                 29.39,-1,0.13,2.45,-0.29,4.5,0.36,0.79,0.97))
 
# c. random draws from G = 1 model estimates:
m2c<-hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
          ng=2,data=data_hlme,B=random(m1))

# d. gridsearch with 50 departures and 10 iterations of the algorithm 
#     (see function gridsearch for details)
## Not run: 
m2d <- gridsearch(rep = 50, maxiter = 10, minit = m1, hlme(Y ~ Time * X1, 
mixture =~ Time, random =~ Time, classmb =~ X2 + X3, subject = 'ID', ng = 2, 
data = data_hlme))


## End(Not run)  
          


# summary of the estimation process
summarytable(m1, m2a, m2b, m2c)

# summary of m2a
summary(m2a)

# posterior classification
postprob(m2a)

# plot of predicted trajectories using some newdata
newdata<-data.frame(Time=seq(0,5,length=100),
X1=rep(0,100),X2=rep(0,100),X3=rep(0,100))
plot(predictY(m2a,newdata,var.time="Time"),legend.loc="right",bty="l")

Estimation of mixed-effect models and latent class mixed-effect models for different types of outcomes (continuous Gaussian, continuous non-Gaussian or ordinal)

Description

This function fits mixed models and latent class mixed models for different types of outcomes. It handles continuous longitudinal outcomes (Gaussian or non-Gaussian) as well as bounded quantitative, discrete and ordinal longitudinal outcomes. The different types of outcomes are taken into account using parameterized nonlinear link functions between the observed outcome and the underlying latent process of interest it measures. At the latent process level, the model estimates a standard linear mixed model or a latent class linear mixed model when heterogeneity in the population is investigated (in the same way as in function hlme). It should be noted that the program also works when no random-effect is included. Parameters of the nonlinear link function and of the latent process mixed model are estimated simultaneously using a maximum likelihood method.

Usage

lcmm(
  fixed,
  mixture,
  random,
  subject,
  classmb,
  ng = 1,
  idiag = FALSE,
  nwg = FALSE,
  link = "linear",
  intnodes = NULL,
  epsY = 0.5,
  cor = NULL,
  data,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  maxiter = 100,
  nsim = 100,
  prior,
  pprior = NULL,
  range = NULL,
  subset = NULL,
  na.action = 1,
  posfix = NULL,
  partialH = FALSE,
  verbose = FALSE,
  returndata = FALSE,
  var.time = NULL,
  nproc = 1,
  clustertype = NULL,
  computeDiscrete = NULL
)

Arguments

Details

A. THE PARAMETERIZED LINK FUNCTIONS

lcmm function estimates mixed models and latent class mixed models for different types of outcomes by assuming a parameterized link function for linking the outcome Y(t) with the underlying latent process L(t) it measures. To fix the latent process dimension, we chose to constrain the (first) intercept of the latent class mixed model at the latent process level at 0 and the standard error of the gaussian error of measurement at 1. These two parameters are replaced by additional parameters in the parameterized link function :

1. With the "linear" link function, 2 parameters are required that correspond directly to the intercept and the standard error: (Y - b1)/b2 = L(t).

2. With the "beta" link function, 4 parameters are required for the following transformation: [ h(Y(t)',b1,b2) - b3]/b4 where h is the Beta CDF with canonical parameters c1 and c2 that can be derived from b1 and b2 as c1=exp(b1)/[exp(b2)*(1+exp(b1))] and c2=1/[exp(b2)*(1+exp(b1))], and Y(t)' is the rescaled outcome i.e. Y(t)'= [ Y(t) - min(Y(t)) + epsY ] / [ max(Y(t)) - min(Y(t)) +2*epsY ].

3. With the "splines" link function, n+2 parameters are required for the following transformation b_1 + b_2*I_1(Y(t)) + ... + b_(n+2) I_(n+1)(Y(t)), where I_1,...,I_(n+1) is the basis of quadratic I-splines. To constraint the parameters to be positive, except for b_1, the program estimates b_k^* (for k=2,...,n+2) so that b_k=(b_k^*)^2.

4. With the "thresholds" link function for an ordinal outcome in levels 0,...,C. A maximumn of C parameters are required for the following transformation: Y(t)=c <=> b_c < L(t) <= b_(c+1) with b_0 = - infinity and b_(C+1)=+infinity. The number of parameters is reduced if some levels do not have any information. For example, if a level c is not observed in the dataset, the corresponding threshold b_(c+1) is constrained to be the same as the previous one b_c. The number of parameters in the link function is reduced by 1.

To constraint the parameters to be increasing, except for the first parameter b_1, the program estimates b_k^* (for k=2,...C) so that b_k=b_(k-1)+(b_k^*)^2.

Details of these parameterized link functions can be found in the referred papers.

B. THE VECTOR OF PARAMETERS B

The parameters in the vector of initial values B or in the vector of maximum likelihood estimates best are included in the following order: (1) ng-1 parameters are required for intercepts in the latent class membership model, and if covariates are included in classmb, ng-1 paramaters should be entered for each one; (2) for all covariates in fixed, one parameter is required if the covariate is not in mixture, ng paramaters are required if the covariate is also in mixture; When ng=1, the intercept is not estimated and no parameter should be specified in B. When ng>1, the first intercept is not estimated and only ng-1 parameters should be specified in B; (3) the variance of each random-effect specified in random (including the intercept) if idiag=TRUE and the inferior triangular variance-covariance matrix of all the random-effects if idiag=FALSE; (4) only if nwg=TRUE, ng-1 parameters for class-specific proportional coefficients for the variance covariance matrix of the random-effects; (5) In contrast with hlme, due to identifiability purposes, the standard error of the Gaussian error is not estimated (fixed at 1), and should not be specified in B; (6) The parameters of the link function: 2 for "linear", 4 for "beta", n+2 for "splines" with n nodes and the number of levels minus one for "thresholds".

C. CAUTIONS REGARDING THE USE OF THE PROGRAM

Some caution should be made when using the program. convergence criteria are very strict as they are based on derivatives of the log-likelihood in addition to the parameter and log-likelihood stability. In some cases, the program may not converge and reach the maximum number of iterations fixed at 100. In this case, the user should check that parameter estimates at the last iteration are not on the boundaries of the parameter space. If the parameters are on the boundaries of the parameter space, the identifiability of the model is critical. This may happen especially with splines parameters that may be too close to 0 (lower boundary) or classmb parameters that are too high or low (perfect classification). When identifiability of some parameters is suspected, the program can be run again from the former estimates by fixing the suspected parameters to their value with option posfix. This usually solves the problem. An alternative is to remove the parameters of the Beta of Splines link function from the inverse of the Hessian with option partialH. If not, the program should be run again with other initial values, with a higher maximum number of iterations or less strict convergence tolerances.

Specifically when investigating heterogeneity (that is with ng>1): (1) As the log-likelihood of a latent class model can have multiple maxima, a careful choice of the initial values is crucial for ensuring convergence toward the global maximum. The program can be run without entering the vector of initial values (see point 2). However, we recommend to systematically enter initial values in B and try different sets of initial values. (2) The automatic choice of initial values we provide requires the estimation of a preliminary linear mixed model. The user should be aware that first, this preliminary analysis can take time for large datatsets and second, that the generated initial values can be very not likely and even may converge slowly to a local maximum. This is the reason why several alternatives exist. The vector of initial values can be directly specified in B the initial values can be generated (automatically or randomly) from a model with ng=. Finally, function gridsearch performs an automatic grid search.

D. NUMERICAL INTEGRATION WITH THE THRESHOLD LINK FUNCTION

With exception for the threshold link function, maximum likelihood estimation implemented in lcmm does not require any numerical integration over the random-effects so that the estimation procedure is relatively fast. See Proust et al. (2006) for more details on the estimation procedure.

However, with the threshold link function and when at least one random-effect is specified, a numerical integration over the random-effects distribution is required in each computation of the individual contribution to the likelihood which complicates greatly the estimation procedure. For the moment, we do not allow any option regarding the numerical integration technics used. 1. When a single random-effect is specified, we use a standard non-adaptive Gaussian quadrature with 30 points. 2. When at least two random-effects are specified, we use a multivariate non-adaptive Gaussian quadrature implemented by Genz (1996) in HRMSYM Fortran subroutine.

Further developments should allow for adaptive technics and more options regarding the numerical integration technic.

E. POSTERIOR DISCRETE LIKELIHOOD

Models involving nonlinear continuous link functions assume the continuous data while the model with a threshold model assumes discrete data. As a consequence, comparing likelihoods or criteria based on the likelihood (as AIC) for these models is not possible as the former are based on a Lebesgue measure and the latter on a counting measure. To make the comparison possible, we compute the posterior discrete likelihood for all the models with a nonlinear continuous link function. This posterior likelihood considers the data as discrete; it is computed at the MLE (maximum likelihood estimates) using the counting measure so that models with threshold or continuous link functions become comparable. Further details can be found in Proust-Lima, Amieva, Jacqmin-Gadda (2012).

In addition to the Akaike information criterion based on the discrete posterior likelihood, we also compute a universal approximate cross-validation criterion to compare models based on a different measure. See Commenges, Proust-Lima, Samieri, Liquet (2015) for further details.

Value

The list returned is:

Author(s)

Cecile Proust-Lima, Amadou Diakite, Benoit Liquet and Viviane Philipps

cecile.proust-lima@inserm.fr

References

Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lcmm. Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02

Genz and Keister (1996). Fully symmetric interpolatory rules for multiple integrals over infinite regions with gaussian weight. Journal of Computational and Applied Mathematics 71: 299-309.

Proust and Jacqmin-Gadda (2005). Estimation of linear mixed models with a mixture of distribution for the random-effects. Comput Methods Programs Biomed 78: 165-73.

Proust, Jacqmin-Gadda, Taylor, Ganiayre, and Commenges (2006). A nonlinear model with latent process for cognitive evolution using multivariate longitudinal data. Biometrics 62: 1014-24.

Proust-Lima, Dartigues and Jacqmin-Gadda (2011). Misuse of the linear mixed model when evaluating risk factors of cognitive decline. Amer J Epidemiol 174(9): 1077-88.

Proust-Lima, Amieva and Jacqmin-Gadda (2013). Analysis of multivariate mixed longitudinal data : a flexible latent process approach, British Journal of Mathematical and Statistical Psychology 66(3): 470-87.

Commenges, Proust-Lima, Samieri, Liquet (2015). A universal approximate cross-validation criterion for regular risk functions. Int J Biostat. 2015 May;11(1):51-67

See Also

postprob, plot.lcmm, plot.predict, hlme

Examples

## Not run: 
#### Estimation of homogeneous mixed models with different assumed link
#### functions, a quadratic mean trajectory for the latent process and 
#### correlated random intercept and slope (the random quadratic slope 
#### was removed as it did not improve the fit of the data).
#### -- comparison of linear, Beta and 3 different splines link functions --
# linear link function
m10<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="linear")
summary(m10)
# Beta link function
m11<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="beta")
summary(m11)
plot(m11,which="linkfunction",bty="l")
# I-splines with 3 equidistant nodes
m12<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="3-equi-splines")
summary(m12)
# I-splines with 5 nodes at quantiles
m13<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="5-quant-splines")
summary(m13)
# I-splines with 5 nodes, and interior nodes entered manually
m14<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="5-manual-splines",intnodes=c(10,20,25))
summary(m14)
plot(m14,which="linkfunction",bty="l")


# Thresholds
# Especially for the threshold link function, we recommend to estimate 
# models with increasing complexity and use estimates of previous ones 
# to specify plausible initial values (we remind that estimation of
# models with threshold link function involves a computationally demanding 
# numerical integration -here of size 3)
m15<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1
,data=data_lcmm,link="thresholds",maxiter=100,
B=c(-0.8379, -0.1103,  0.3832,  0.3788 , 0.4524, -7.3180,  0.5917,  0.7364,
 0.6530, 0.4038,  0.4290,  0.6099,  0.6014 , 0.5354 , 0.5029 , 0.5463,
 0.5310 , 0.5352, 0.6498,  0.6653,  0.5851,  0.6525,  0.6701 , 0.6670 ,
 0.6767 , 0.7394 , 0.7426, 0.7153,  0.7702,  0.6421))
summary(m15)
plot(m15,which="linkfunction",bty="l")

#### Plot of estimated different link functions:
#### (applicable for models that only differ in the "link function" used. 
####  Otherwise, the latent process scale is different and a rescaling
####  is necessary)
plot(m10,which="linkfunction",col=1,xlab="latent process",ylab="marker",
bty="l",xlim=c(-10,5),legend=NULL)
plot(m11,which="linkfunction",add=TRUE,col=2,legend=NULL)
plot(m12,which="linkfunction",add=TRUE,col=3,legend=NULL)
plot(m13,which="linkfunction",add=TRUE,col=4,legend=NULL)
plot(m14,which="linkfunction",add=TRUE,col=5,legend=NULL)
plot(m15,which="linkfunction",add=TRUE,col=6,legend=NULL)
legend(x="bottomright",legend=c("linear","beta","spl_3e","spl_5q","spl_5m","thresholds"),
col=1:6,lty=1,inset=.02,box.lty=0)

#### Estimation of 2-latent class mixed models with different assumed link 
#### functions with individual and class specific linear trend
#### for illustration, only default initial values where used but other
#### sets of initial values should also be tried to ensure convergence 
#### towards the golbal maximum
# Linear link function
m20<-lcmm(Ydep2~Time,random=~Time,subject='ID',mixture=~Time,ng=2,
idiag=TRUE,data=data_lcmm,link="linear",B=c(-0.98,0.79,-2.09,
-0.81,0.19,0.55,24.49,2.24))
summary(m20)
postprob(m20)
# Beta link function
m21<-lcmm(Ydep2~Time,random=~Time,subject='ID',mixture=~Time,ng=2,
idiag=TRUE,data=data_lcmm,link="beta",B=c(-0.1,-0.56,-0.4,-1.77,
0.53,0.14,0.6,-0.83,0.73,0.09))
summary(m21)
postprob(m21)
# I-splines link function (and 5 nodes at quantiles)
m22<-lcmm(Ydep2~Time,random=~Time,subject='ID',mixture=~Time,ng=2,
idiag=TRUE,data=data_lcmm,link="5-quant-splines",B=c(0.12,0.63,
-1.76,-0.39,0.51,0.13,-7.37,1.05,1.28,1.96,1.3,0.93,1.05))
summary(m22)
postprob(m22)

data <- data_lcmm[data_lcmm$ID==193,]
plot(predictL(m22,var.time="Time",newdata=data,bty="l")


## End(Not run)    

Wrapper to the Fortran subroutines computing the log-likelihood

Description

Log-likelihood of hlme, lcmm, multlcmm, Jointlcmm and mpjlcmm models. The argument's specification is not straightforward, so these functions are usually not directly used.

Usage

loglikhlme(
  b,
  Y0,
  X0,
  prior0,
  pprior0,
  idprob0,
  idea0,
  idg0,
  idcor0,
  ns0,
  ng0,
  nv0,
  nobs0,
  nea0,
  nmes0,
  idiag0,
  nwg0,
  ncor0,
  npm0,
  fix0,
  nfix0,
  bfix0
)

logliklcmm(
  b,
  Y0,
  X0,
  prior0,
  pprior0,
  idprob0,
  idea0,
  idg0,
  idcor0,
  ns0,
  ng0,
  nv0,
  nobs0,
  nea0,
  nmes0,
  idiag0,
  nwg0,
  ncor0,
  npm0,
  epsY0,
  idlink0,
  nbzitr0,
  zitr0,
  minY0,
  maxY0,
  ide0,
  fix0,
  nfix0,
  bfix0
)

loglikmultlcmm(
  b,
  Y0,
  X0,
  prior0,
  pprior0,
  idprob0,
  idea0,
  idg0,
  idcor0,
  idcontr0,
  ny0,
  ns0,
  ng0,
  nv0,
  nobs0,
  nea0,
  nmes0,
  idiag0,
  nwg0,
  ncor0,
  nalea0,
  npm0,
  epsY0,
  idlink0,
  nbzitr0,
  zitr0,
  uniqueY0,
  indiceY0,
  nvalSPLORD0,
  fix0,
  nfix0,
  bfix0,
  methInteg0,
  nMC0,
  dimMC0,
  seqMC0,
  chol0
)

loglikJointlcmm(
  b,
  Y0,
  X0,
  prior0,
  pprior0,
  tentr0,
  tevt0,
  devt0,
  ind_survint0,
  idprob0,
  idea0,
  idg0,
  idcor0,
  idcom0,
  idspecif0,
  idtdv0,
  idlink0,
  epsY0,
  nbzitr0,
  zitr0,
  uniqueY0,
  nvalSPL0,
  indiceY0,
  typrisq0,
  risqcom0,
  nz0,
  zi0,
  ns0,
  ng0,
  nv0,
  nobs0,
  nmes0,
  nbevt0,
  nea0,
  nwg0,
  ncor0,
  idiag0,
  idtrunc0,
  logspecif0,
  npm0,
  fix0,
  nfix0,
  bfix0
)

loglikmpjlcmm(
  b,
  K0,
  ny0,
  nbevt0,
  ng0,
  ns0,
  Y0,
  nobs0,
  X0,
  nv0,
  Xns0,
  nv20,
  prior0,
  Tentr0,
  Tevt0,
  Devt0,
  ind_survint0,
  idnv0,
  idnv20,
  idspecif0,
  idlink0,
  epsY0,
  nbzitr0,
  zitr0,
  uniqueY0,
  nvalSPL0,
  indiceY0,
  typrisq0,
  risqcom0,
  nz0,
  zi0,
  nmes0,
  nea0,
  nw0,
  ncor0,
  nalea0,
  idiag0,
  idtrunc0,
  logspecif0,
  npm0,
  fix0,
  contrainte0,
  nfix0,
  bfix0
)

Arguments

Value

the log-likelihood

Author(s)

Cecile Proust-Lima, Viviane Philipps

Estimation of multivariate joint latent class mixed models

Description

This function fits joint latent class models for multivariate longitudinal markers and competing causes of event. It is a multivariate extension of the Jointlcmm function. It defines each longitudinal dimension as a latent process (mp in mpjlcmm is for multivariate processes), possibly measured by sereval continuous markers (Gaussian or curvilinear). For the moment, theses processes are assumed independent given the latent classes. The (optional) survival part handles competing risks, right censoring and left truncation. The specification of the function is similar to other estimating functions of the package.

Usage

mpjlcmm(
  longitudinal,
  subject,
  classmb,
  ng,
  survival,
  hazard = "Weibull",
  hazardtype = "Specific",
  hazardnodes = NULL,
  TimeDepVar = NULL,
  data,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  maxiter = 100,
  nsim = 100,
  prior,
  logscale = FALSE,
  subset = NULL,
  na.action = 1,
  posfix = NULL,
  partialH = FALSE,
  verbose = FALSE,
  nproc = 1,
  clustertype = NULL
)

Arguments

Author(s)

Cecile Proust Lima and Viviane Philipps

Examples

## Not run: 
paquid$age65 <- (paquid$age-65)/10

##############################################################################
###                          EXAMPLE 1 :                                   ###
###two outcomes measuring the same latent process along with dementia onset###
##############################################################################

## multlcmm model for MMSE and BVRT for 1 class
mult1 <- multlcmm(MMSE+BVRT~age65+I(age65^2)+CEP+male,random=~age65+I(age65^2),
subject="ID",link=c("5-quant-splines","4-quant-splines"),data=paquid)
summary(mult1)

## joint model for 1 class
m1S <- mpjlcmm(longitudinal=list(mult1),subject="ID",ng=1,data=paquid,
survival=Surv(age_init,agedem,dem)~1)
summary(m1S)


##### joint model for 2 classes #####

## specify longitudinal model for 2 classes, without estimation
mult2 <- multlcmm(MMSE+BVRT~age65+I(age65^2)+CEP+male,random=~age65+I(age65^2),
subject="ID",link=c("5-quant-splines","4-quant-splines"),ng=2,
mixture=~age65+I(age65^2),data=paquid,B=random(mult1),maxiter=0)

## estimation of the associated joint model 
m2S <- mpjlcmm(longitudinal=list(mult2),subject="ID",ng=2,data=paquid,
survival=Surv(age_init,agedem,dem)~1)

## estimation by a grid search with 50 replicates, initial values
## randomly generated from m1S
m2S_b <- gridsearch(mpjlcmm(longitudinal=list(mult2),subject="ID",ng=2,
data=paquid,survival=Surv(age_init,agedem,dem)~1), minit=m1S, rep=50, maxiter=30)


##### joint model for 3 classes #####
mult3 <- multlcmm(MMSE+BVRT~age65+I(age65^2)+CEP+male,random=~age65+I(age65^2),
subject="ID",link=c("5-quant-splines","4-quant-splines"),ng=3,
mixture=~age65+I(age65^2),data=paquid,B=random(mult1),maxiter=0)

m3S <- mpjlcmm(longitudinal=list(mult3),subject="ID",ng=3,data=paquid,
survival=Surv(age_init,agedem,dem)~1)

m3S_b <- gridsearch(mpjlcmm(longitudinal=list(mult3),subject="ID",ng=3,
data=paquid,survival=Surv(age_init,agedem,dem)~1), minit=m1S, rep=50, maxiter=30)


##### summary of the models #####

summarytable(m1S,m2S,m2S_b,m3S,m3S_b)



##### post-fit #####

## update longitudinal models :
mod2 <- update(m2S)

mult2_post <- mod2[[1]]
## -> use the available functions for multlcmm on the mult2_post object

## fit of the longitudinal trajectories
par(mfrow=c(2,2))
plot(mult2_post,"fit","age65",marg=TRUE,shades=TRUE,outcome=1)
plot(mult2_post,"fit","age65",marg=TRUE,shades=TRUE,outcome=2)

plot(mult2_post,"fit","age65",marg=FALSE,shades=TRUE,outcome=1)
plot(mult2_post,"fit","age65",marg=FALSE,shades=TRUE,outcome=2)


## predicted trajectories
dpred <- data.frame(age65=seq(0,3,0.1),male=0,CEP=0)

predL <- predictL(mult2_post,newdata=dpred,var.time="age65",confint=TRUE)
plot(predL,shades=TRUE) # in the latent process scale


predY <- predictY(mult2_post,newdata=dpred,var.time="age65",draws=TRUE)

plot(predY,shades=TRUE,ylim=c(0,30),main="MMSE") #in the 0-30 scale for MMSE
plot(predY,shades=TRUE,ylim=c(0,15),outcome=2,main="BVRT") #in 0-15 for BVRT

## baseline hazard and survival curves :
plot(m2S,"hazard")
plot(m2S,"survival")

## posteriori probabilities and classification :
postprob(m2S)



####################################################################################
###                              EXAMPLE 2 :                                     ###
### two latent processes measured each by one outcome along with dementia onset  ###
####################################################################################

## define the two longitudinal models

mMMSE1 <- lcmm(MMSE~age65+I(age65^2)+CEP,random=~age65+I(age65^2),subject="ID",
link="5-quant-splines",data=paquid)

mCESD1 <- lcmm(CESD~age65+I(age65^2)+male,random=~age65+I(age65^2),subject="ID",
link="5-quant-splines",data=paquid)


## joint estimation

mm1S <- mpjlcmm(longitudinal=list(mMMSE1,mCESD1),subject="ID",ng=1,data=paquid,
survival=Surv(age_init,agedem,dem)~CEP+male)


## with 2 latent classes

mMMSE2 <- lcmm(MMSE~age65+I(age65^2)+CEP,random=~age65+I(age65^2),subject="ID",
link="5-quant-splines",data=paquid,ng=2,mixture=~age65+I(age65^2),
B=random(mMMSE1),maxiter=0)

mCESD2 <- lcmm(CESD~age65+I(age65^2)+male,random=~age65+I(age65^2),subject="ID",
link="5-quant-splines",data=paquid,ng=2,mixture=~age65+I(age65^2),
B=random(mCESD1),maxiter=0)

mm2S <- mpjlcmm(longitudinal=list(mMMSE2,mCESD2),subject="ID",ng=2,data=paquid,
survival=Surv(age_init,agedem,dem)~CEP+mixture(male),classmb=~CEP+male)

mm2S_b <- gridsearch(mpjlcmm(longitudinal=list(mMMSE2,mCESD2),subject="ID",ng=2,
data=paquid,survival=Surv(age_init,agedem,dem)~CEP+mixture(male),
classmb=~CEP+male),minit=mm1S,rep=50,maxiter=50)

summarytable(mm1S,mm2S,mm2S_b)


mod1_biv <- update(mm1S)
mod2_biv <- update(mm2S)

## -> use post-fit functions as in exemple 1

## End(Not run)

Estimation of multivariate mixed-effect models and multivariate latent class mixed-effect models for multivariate longitudinal outcomes of possibly multiple types (continuous Gaussian, continuous non-Gaussian/curvilinear, ordinal) that measure the same underlying latent process.

Description

This function constitutes a multivariate extension of function lcmm. It fits multivariate mixed models and multivariate latent class mixed models for multivariate longitudinal outcomes of different types. It handles continuous longitudinal outcomes (Gaussian or non-Gaussian, curvilinear) as well as ordinal longitudinal outcomes (with cumulative probit measurement model). The model assumes that all the outcomes measure the same underlying latent process defined as their common factor, and each outcome is related to this latent common factor by a specific parameterized link function. At the latent process level, the model estimates a standard linear mixed model or a latent class linear mixed model when heterogeneity in the population is investigated (in the same way as in functions hlme and lcmm). Parameters of the nonlinear link functions and of the latent process mixed model are estimated simultaneously using a maximum likelihood method.

Usage

multlcmm(
  fixed,
  mixture,
  random,
  subject,
  classmb,
  ng = 1,
  idiag = FALSE,
  nwg = FALSE,
  randomY = FALSE,
  link = "linear",
  intnodes = NULL,
  epsY = 0.5,
  cor = NULL,
  data,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  maxiter = 100,
  nsim = 100,
  prior,
  pprior = NULL,
  range = NULL,
  subset = NULL,
  na.action = 1,
  posfix = NULL,
  partialH = FALSE,
  verbose = FALSE,
  returndata = FALSE,
  methInteg = "QMC",
  nMC = NULL,
  var.time = NULL,
  nproc = 1,
  clustertype = NULL
)

mlcmm(
  fixed,
  mixture,
  random,
  subject,
  classmb,
  ng = 1,
  idiag = FALSE,
  nwg = FALSE,
  randomY = FALSE,
  link = "linear",
  intnodes = NULL,
  epsY = 0.5,
  cor = NULL,
  data,
  B,
  convB = 1e-04,
  convL = 1e-04,
  convG = 1e-04,
  maxiter = 100,
  nsim = 100,
  prior,
  pprior = NULL,
  range = NULL,
  subset = NULL,
  na.action = 1,
  posfix = NULL,
  partialH = FALSE,
  verbose = FALSE,
  returndata = FALSE,
  methInteg = "QMC",
  nMC = NULL,
  var.time = NULL,
  nproc = 1,
  clustertype = NULL
)

Arguments

Details

A. THE PARAMETERIZED LINK FUNCTIONS

multlcmm function estimates multivariate latent class mixed models for different types of outcomes by assuming a parameterized link function for linking each outcome Y_k(t) with the underlying latent common factor L(t) they measure. To fix the latent process dimension, we chose to constrain at the latent process level the (first) intercept of the latent class mixed model at 0 and the standard error of the first random effect at 1.

1. With the "linear" link function, 2 parameters are required for the following transformation (Y(t) - b1)/b2

2. With the "beta" link function, 4 parameters are required for the following transformation: [ h(Y(t)',b1,b2) - b3]/b4 where h is the Beta CDF with canonical parameters c1 and c2 that can be derived from b1 and b2 as c1=exp(b1)/[exp(b2)*(1+exp(b1))] and c2=1/[exp(b2)*(1+exp(b1))], and Y(t)' is the rescaled outcome i.e. Y(t)'= [ Y(t) - min(Y(t)) + epsY ] / [ max(Y(t)) - min(Y(t)) +2*epsY ].

3. With the "splines" link function, n+2 parameters are required for the following transformation b_1 + b_2*I_1(Y(t)) + ... + b_(n+2) I_(n+1)(Y(t)), where I_1,...,I_(n+1) is the basis of quadratic I-splines. To constraint the parameters to be positive, except for b_1, the program estimates b_k^* (for k=2,...,n+2) so that b_k=(b_k^*)^2. This parameterization may lead in some cases to problems of convergence that we are currently addressing.

4. With the "thresholds" link function for an ordinal outcome with levels 0,...,C, C-1 parameters are required for the following transformation: Y(t)=c <=> b_c < L(t) <= b_(c+1) with b_0 = - infinity and b_(C+1)=+infinity. To constraint the parameters to be increasing, except for the first parameter b_1, the program estimates b_k^* (for k=2,...C-1) so that b_k=b_(k-1)+(b_k^*)^2.

Details of these parameterized link functions can be found in the papers: Proust-Lima et al. (Biometrics 2006), Proust-Lima et al. (BJMSP 2013), Proust-Lima et al. (arxiv 2021 - https://arxiv.org/abs/2109.13064)

B. THE VECTOR OF PARAMETERS B

The parameters in the vector of initial values B or in the vector of maximum likelihood estimates best are included in the following order: (1) ng-1 parameters are required for intercepts in the latent class membership model, and if covariates are included in classmb, ng-1 paramaters should be entered for each one; (2) for all covariates in fixed, one parameter is required if the covariate is not in mixture, ng paramaters are required if the covariate is also in mixture; When ng=1, the intercept is not estimated and no intercept parameter should be specified in B. When ng>1, the first intercept is not estimated and only ng-1 intercept parameters should be specified in B; (3) for all covariates included with contrast() in fixed, one supplementary parameter per outcome is required excepted for the last outcome for which the parameter is not estimated but deduced from the others; (4) if idiag=TRUE, the variance of each random-effect specified in random is required excepted the first one (usually the intercept) which is constrained to 1. (5) if idiag=FALSE, the inferior triangular variance-covariance matrix of all the random-effects is required excepted the first variance (usually the intercept) which is constrained to 1. (6) only if nwg=TRUE and ng>1, ng-1 parameters for class-specific proportional coefficients for the variance covariance matrix of the random-effects; (7) if cor is specified, the standard error of the Brownian motion or the standard error and the correlation parameter of the autoregressive process; (8) the standard error of the outcome-specific Gaussian errors (one per outcome); (9) if randomY=TRUE, the standard error of the outcome-specific random intercept (one per outcome); (10) the parameters of each parameterized link function: 2 for "linear", 4 for "beta", n+2 for "splines" with n nodes.

C. CAUTIONS REGARDING THE USE OF THE PROGRAM

Some caution should be made when using the program. Convergence criteria are very strict as they are based on the derivatives of the log-likelihood in addition to the parameter and log-likelihood stability. In some cases, the program may not converge and reach the maximum number of iterations fixed at 100. In this case, the user should check that parameter estimates at the last iteration are not on the boundaries of the parameter space.

If the parameters are on the boundaries of the parameter space, the identifiability of the model is critical. This may happen especially with splines parameters that may be too close to 0 (lower boundary) or classmb parameters that are too high or low (perfect classification). When identifiability of some parameters is suspected, the program can be run again from the former estimates by fixing the suspected parameters to their value with option posfix. This usually solves the problem. An alternative is to remove the parameters of the Beta or Splines link function from the inverse of the Hessian with option partialH.

If not, the program should be run again with other initial values, with a higher maximum number of iterations or less strict convergence tolerances.

Specifically when investigating heterogeneity (that is with ng>1): (1) As the log-likelihood of a latent class model can have multiple maxima, a careful choice of the initial values is crucial for ensuring convergence toward the global maximum. The program can be run without entering the vector of initial values (see point 2). However, we recommend to systematically enter initial values in B and try different sets of initial values. (2) The automatic choice of initial values we provide requires the estimation of a preliminary linear mixed model. The user should be aware that first, this preliminary analysis can take time for large datatsets and second, that the generated initial values can be very not likely and even may converge slowly to a local maximum. This is the reason why several alternatives exist. The vector of initial values can be directly specified in B the initial values can be generated (automatically or randomly) from a model with ng=. Finally, function gridsearch performs an automatic grid search.

D. NUMERICAL INTEGRATION WITH THE THRESHOLD LINK FUNCTION

When dealing only with continuous outcomes, the computation of the likelihood does not require any numerical integration over the random-effects, so that the estimation procedure is relatively fast. When at least one ordinal outcome is modeled, a numerical integration over the random-effects is required in each computation of the individual contribution to the likelihood. This achieved using a Monte-Carlo procedure. We allow three options: the standard Monte-Carlo simulations, as well as antithetic Monte-Carlo and quasi Monte-Carlo methods as proposed in Philipson et al (2020).

Value

The list returned is:

Author(s)

Cecile Proust-Lima and Viviane Philipps

cecile.proust-lima@inserm.fr

References

Proust-Lima C, Philipps V, Liquet B (2017). Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lcmm. Journal of Statistical Software, 78(2), 1-56. doi:10.18637/jss.v078.i02

Proust and Jacqmin-Gadda (2005). Estimation of linear mixed models with a mixture of distribution for the random-effects. Comput Methods Programs Biomed 78: 165-73.

Proust, Jacqmin-Gadda, Taylor, Ganiayre, and Commenges (2006). A nonlinear model with latent process for cognitive evolution using multivariate longitudinal data. Biometrics 62, 1014-24.

Proust-Lima, Dartigues and Jacqmin-Gadda (2011). Misuse of the linear mixed model when evaluating risk factors of cognitive decline. Amer J Epidemiol 174(9): 1077-88.

Proust-Lima, Amieva, Jacqmin-Gadda (2013). Analysis of multivariate mixed longitudinal data: A flexible latent process approach. Br J Math Stat Psychol 66(3): 470-87.

Commenges, Proust-Lima, Samieri, Liquet (2012). A universal approximate cross-validation criterion and its asymptotic distribution, Arxiv.

Philipson, Hickey, Crowther, Kolamunnage-Dona (2020). Faster Monte Carlo estimation of semiparametric joint models of time-to-event and multivariate longitudinal data. Computational Statistics & Data Analysis 151.

Proust-Lima, Philipps, Perrot, Blanchin, Sebille (2021). Modeling repeated self-reported outcome data: a continuous-time longitudinal Item Response Theory model. https://arxiv.org/abs/2109.13064

See Also

postprob, plot.multlcmm, predictL, predictY lcmm

Examples

## Not run: 
# Latent process mixed model for two curvilinear outcomes. Link functions are 
# aproximated by I-splines, the first one has 3 nodes (i.e. 1 internal node 8),
# the second one has 4 nodes (i.e. 2 internal nodes 12,25)

m1 <- multlcmm(Ydep1+Ydep2~1+Time*X2+contrast(X2),random=~1+Time,
subject="ID",randomY=TRUE,link=c("4-manual-splines","3-manual-splines"),
intnodes=c(8,12,25),data=data_lcmm)

# to reduce the computation time, the same model is estimated using 
# a vector of initial values
m1 <- multlcmm(Ydep1+Ydep2~1+Time*X2+contrast(X2),random=~1+Time,
subject="ID",randomY=TRUE,link=c("4-manual-splines","3-manual-splines"),
intnodes=c(8,12,25),data=data_lcmm, 
B=c(-1.071, -0.192,  0.106, -0.005, -0.193,  1.012,  0.870,  0.881,
  0.000,  0.000, -7.520,  1.401,  1.607 , 1.908,  1.431,  1.082,
 -7.528,  1.135 , 1.454 , 2.328, 1.052))


# output of the model
summary(m1)
# estimated link functions
plot(m1,which="linkfunction")
# variation percentages explained by linear mixed regression
VarExpl(m1,data.frame(Time=0))

#### Heterogeneous latent process mixed model with linear link functions 
#### and 2 latent classes of trajectory 
m2 <- multlcmm(Ydep1+Ydep2~1+Time*X2,random=~1+Time,subject="ID",
link="linear",ng=2,mixture=~1+Time,classmb=~1+X1,data=data_lcmm,
B=c( 18,-20.77,1.16,-1.41,-1.39,-0.32,0.16,-0.26,1.69,1.12,1.1,10.8,
1.24,24.88,1.89))
# summary of the estimation
summary(m2)
# posterior classification
postprob(m2)
# longitudinal predictions in the outcomes scales for a given profile of covariates 
newdata <- data.frame(Time=seq(0,5,length=100),X1=0,X2=0,X3=0)
predGH <- predictY(m2,newdata,var.time="Time",methInteg=0,nsim=20) 
head(predGH)

## End(Not run)

Longitudinal data on cognitive and physical aging in the elderly

Description

The dataset consists in a subsample of the Paquid prospective cohort study. Repeated measures cognitive measures (MMSE, IST, BVRT psychometric tests), physical dependency (HIER) and depression sympatomatology (CESD) were collected over a maximum period of 20 years along with dementia information (age at dementia diagnosis, dementia diagnosis information). Time-independent socio-demographic information is also provided (CEP, male, age_init).

Format

A data frame with 2250 observations over 500 subjects and 12 variables:

ID

subject identification number

MMSE

score at the Mini-Mental State Examination (MMSE), a psychometric test of global cognitive functioning (integer in range 0-30)

BVRT

score at the Benton Visual Retention Test (BVRT), a psychometric test of spatial memory (integer in range 0-15)

IST

score at the Isaacs Set Test (IST) truncated at 15 seconds, a test of verbal memory (integer in range 0-40)

HIER

score of physical dependency (0=no dependency, 1=mild dependency, 2=moderate dependency, 3=severe dependency)

CESD

score of a short self-report scale CES-D designed to measure depressive symptomatology in the general population (integer in range 0-52)

age

age at the follow-up visit

dem

indicator of positive diagnosis of dementia

agedem

age at dementia diagnosis for dem=1 and at last contact for dem=0

age_init

age at entry in the cohort

CEP

binary indicator of educational level (CEP=1 for subjects who graduated from primary school; CEP=0 otherwise)

male

binary indicator for gender (male=1 for men; male=0 for women)

References

Letenneur, L., Commenges, D., Dartigues, J. F., & Barberger-Gateau, P. (1994). Incidence of dementia and Alzheimer's disease in elderly community residents of southwestern France. International Journal of Epidemiology, 23 (6), 1256-61.

Examples

summary(paquid)

Permutation of the latent classes

Description

This function allows a re-ordering of the latent classes of an estimated model.

Usage

permut(m, order, estim = TRUE)

Arguments

Value

An object of the same class as m, with reordered classes, or the initial object with new coefficients if estim is FALSE.

Author(s)

Viviane Philipps and Cecile Proust-Lima

Examples

## Estimation of a hlme model with 2 classes
m2 <- hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
         ng=2,data=data_hlme,B=c(0.11,-0.74,-0.07,20.71,
                                 29.39,-1,0.13,2.45,-0.29,4.5,0.36,0.79,0.97))

## Exchange class 2 and class 1
m2b <- permut(m2,order=c(2,1))

Plot of a fitted model

Description

This function produces different plots (residuals, goodness-of-fit, estimated link functions, estimated baseline risk/survival and posterior probabilities distributions) of a fitted object of class hlme, lcmm, multlcmm or Jointlcmm.

Usage

## S3 method for class 'hlme'
plot(x, which = "residuals", var.time, break.times, marg, subset, shades, ...)

## S3 method for class 'lcmm'
plot(x, which = "residuals", var.time, break.times, marg, subset, shades, ...)

## S3 method for class 'multlcmm'
plot(
  x,
  which = "residuals",
  var.time,
  break.times,
  marg,
  outcome,
  subset,
  shades,
  ...
)

## S3 method for class 'Jointlcmm'
plot(
  x,
  which = "residuals",
  var.time,
  break.times,
  marg,
  event,
  subset,
  shades,
  ...
)

## S3 method for class 'mpjlcmm'
plot(x, which, event, ...)

## S3 method for class 'externSurv'
plot(x, which = "hazard", event, ...)

## S3 method for class 'externX'
plot(x, which = "postprob", event, ...)

Arguments

Details

With which="residuals", this function provides the marginal residuals against the marginal predictions, the subject-specific residuals against the subject-specific predictions, a normal QQ-plot with confidence bands for the marginal residuals and a normal QQ-plot with confidence bands for the subject-specific residuals.

With which="postprob", the function provides the histograms of the posterior class-membership probabilities stemmed from a Jointlcmm, lcmm, hlme or multlcmm object.

With which="link" or which="linkfunction", the function displays the estimated transformation(s) specified in the option link of lcmm and multlcmm functions. It corresponds to the (non)linear parameterized link estimated between the oberved longitudinal outcome and the underlying latent process.

With which="fit", the function provides the class-specific weighted marginal and subject-specific mean predicted trajectories with time and the class-specific weighted mean observed trajectories and their 95% confidence bounds. The predicted and observed class-specific values are weighted means within each time interval; For each observation or prediction (in the transformed scale if appropriate), the weights are the class-specific (posterior with subject-specific or marginal otherwise) probabilities to belong to the latent class.

With which="baselinerisk" or which="hazard", the function displays the estimated baseline risk functions for the time-to-event of interest in each latent class.

With which="survival", the function displays the estimated event-free probabilities (survival functions) for the time-to-event of interest in each latent class.

Author(s)

Cecile Proust-Lima, Viviane Philipps and Benoit Liquet

See Also

hlme, lcmm, multlcmm, Jointlcmm

Examples

###################### fit, residuals and postprob 

# estimation of the model
m<-lcmm(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,
subject='ID',ng=2,data=data_hlme,B=c(0.41,0.55,-0.18,-0.41,
-14.26,-0.34,1.33,13.51,24.65,2.98,1.18,26.26,0.97))  

# fit
plot(m,which="fit",marg=FALSE,var.time="Time",bty="n")
# residuals plot
plot(m)
# postprob plot
plot(m,which="postprob") 


###################### fit, linkfunctions

#### Estimation of homogeneous mixed models with different assumed link
#### functions, a quadratic mean trajectory for the latent process with 
#### independent random intercept, slope and quadratic slope
#### (comparison of linear, Beta and 3 and 5 splines link functions)
## Not run: 

# linear link function
m10<-lcmm(Ydep2~Time+I(Time^2),random=~Time+I(Time^2),subject='ID',ng=1,
          data=data_lcmm,link="linear",
          B=c(-0.7454, -0.2031,  0.2715,  0.2916 , 0.6114, -0.0064,  0.0545,
              0.0128, 25.3795, 2.2371))
            
# Beta link function
m11<-lcmm(Ydep2~Time+I(Time^2),random=~Time+I(Time^2),subject='ID',ng=1,
          data=data_lcmm,link="beta",B=c(-0.9109, -0.0831,  0.5194,  0.1910 ,
          0.8984, -0.0179, -0.0636,  0.0045,  0.5514, -0.7692,  0.7037,  0.0899))
          
# fit 
par(mfrow=c(2,1),mar=c(4,4,1,1))
plot(m11,which="fit",var.time="Time",bty="l",ylim=c(-3,0))
plot(m11,which="fit",var.time="Time",marg=FALSE,bty="l",ylim=c(-3,0))

# I-splines with 3 equidistant nodes
m12<-lcmm(Ydep2~Time+I(Time^2),random=~Time+I(Time^2),subject='ID',ng=1,
          data=data_lcmm,link="3-equi-splines",B=c(-0.9272, -0.0753 , 0.5304, 
          0.1950,  0.9260, -0.0204, -0.0739 , 0.0059, -7.8369,  0.9228 ,-1.4689,
          2.0396,  1.8102))

# I-splines with 5 nodes, and interior nodes entered manually
m13<-lcmm(Ydep2~Time+I(Time^2),random=~Time+I(Time^2),subject='ID',ng=1,
          data=data_lcmm,link="5-manual-splines",intnodes=c(10,20,25),
          B=c(-0.9315, -0.0739 , 0.5254 , 0.1933,  0.9418, -0.0206, -0.0776,
          0.0064, -7.8645, 0.7470,  1.2080,  1.5537 , 1.7558 , 1.3386 , 1.0982))

# Plot of estimated different link functions:
# (applicable for models that only differ in the "link function" used. 
# Otherwise, the latent process scale is different and a rescaling
# is necessary)
plot(m10,which="linkfunction",bty="l")
plot(m11,which="linkfunction",bty="l",add=TRUE,col=2)
plot(m12,which="linkfunction",bty="l",add=TRUE,col=3)
plot(m13,which="linkfunction",bty="l",add=TRUE,col=4)
legend("topleft",legend=c("linear","beta","3-Isplines","5-Isplines"),
col=1:4,lty=1,bty='n')

## End(Not run)


###################### fit, baselinerisk and survival
## Not run: 
#### estimation with 3 latent classes (ng=3) - see Jointlcmm 
#### help for details on the model
m3 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~ X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",ng=3,data=data_lcmm,
B=c(0.7576, 0.4095, -0.8232, -0.2737, 0, 0, 0, 0.2838, -0.6338, 
2.6324, 5.3963, -0.0273, 1.3979, 0.8168, -15.041, 10.164, 10.2394, 
11.5109, -2.6219, -0.4553, -0.6055, 1.473, -0.0383, 0.8512, 0.0389, 
0.2624, 1.4982))

# fit
plot(m3,which="fit",var.time="Time",bty="l")
plot(m3,which="fit",var.time="Time",marg=FALSE,bty="l",ylim=c(0,15))


# Class-specific predicted baseline risk & survival functions in the 
# 3-class model retained (for the reference value of the covariates) 
plot(m3,which="baselinerisk",bty="l")
plot(m3,which="baselinerisk",ylim=c(0,5),bty="l")
plot(m3,which="survival",bty="l")

## End(Not run)

Plots

Description

This function displays plots related to predictive accuracy functions: epoce and Diffepoce.

Usage

## S3 method for class 'Diffepoce'
plot(x, ...)

## S3 method for class 'epoce'
plot(x, ...)

Arguments

Details

These functions do not apply for the moment with multiple causes of event (competing risks).

For epoce objects, the function displays the EPOCE estimate (either MPOL or CVPOL) according to the time of prediction. For Diffepoce objects, plot displays the difference in EPOCE estimates (either MPOL or CVPOL) and its 95% tracking interval between two joint latent class models

Value

Returns plots related to epoce and Diffepoce

Author(s)

Cecile Proust-Lima and Viviane Philipps

See Also

epoce,Diffepoce

Examples

## Not run: 
# estimation of the joint latent class model
m3 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",ng=3,data=data_lcmm,
B=c(0.7667, 0.4020, -0.8243, -0.2726, 0.0000, 0.0000, 0.0000, 0.3020,
-0.6212, 2.6247, 5.3139, -0.0255, 1.3595, 0.8172, -11.6867, 10.1668,
10.2355, 11.5137, -2.6209, -0.4328, -0.6062, 1.4718, -0.0378, 0.8505,
0.0366, 0.2634, 1.4981))
# predictive accuracy of the model evaluated with EPOCE
VecTime <- c(1,3,5,7,9,11,13,15)
cvpl <- epoce(m3,var.time="Time",pred.times=VecTime)
summary(cvpl)
plot(cvpl,bty="l",ylim=c(0,2))

## End(Not run)

Plot of information functions

Description

This function plots the information functions stemmed from a lcmm or multlcmm object with ordinal outcomes modeled via threshold links.

Usage

## S3 method for class 'ItemInfo'
plot(
  x,
  which = "ItemInfo",
  outcome = "all",
  legend.loc = "topright",
  legend = NULL,
  add = FALSE,
  shades = TRUE,
  ...
)

Arguments

Author(s)

Viviane Philipps and Cecile Proust-Lima

Plot of predicted cumulative incidences according to a profile of covariates

Description

This function displays the predicted cause-specific cumulative incidences derived from a joint latent class model according to a profile of covariates. does. ~~

Usage

## S3 method for class 'cuminc'
plot(
  x,
  profil = 1,
  event = 1,
  add = FALSE,
  legend,
  legend.loc = "topleft",
  ...
)

Arguments

Value

returns NULL

Author(s)

Viviane Philipps and Cecile Proust-Lima

See Also

Jointlcmm, plot.Jointlcmm, cuminc

Plot of individual dynamic predictions

Description

This function provides a graphical representation of individual dynamic predictions obtained from a joint latent class model and plots simultaneously the observed outcome.

Usage

## S3 method for class 'dynpred'
plot(x, subject = NULL, landmark = NULL, horizon = NULL, add = FALSE, ...)

Arguments

Details

Two types of plot are provided for the moment :

- if one horizon is selected (and one or several landmarks), each prediction is represented by a point at the landmark time. If available, the predictions are surrounded by confidence intervals.

- if several horizons (t1, t2, etc) and only one landmark (s) is selected, a line linking the predictions (placed at abscissa s+t1, s+t2, etc) is drawn. Confidence bands (if available) are represented as dotted lines.

Value

returns NULL

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

dynpred

Examples

## Not run: 

## Joint latent class model with 2 classes :
m32 <- Jointlcmm(Ydep1~Time*X1,mixture=~Time,random=~Time,subject="ID",
classmb=~X3,ng=2,survival=Surv(Tevent,Event)~X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",data=data_lcmm,B = c(0.64, -0.62, 
0, 0, 0.52, 0.81, 0.41, 0.78, 0.1, 0.77, -0.05, 10.43, 11.3, -2.6, -0.52, 1.41, 
-0.05, 0.91, 0.05, 0.21, 1.5))

## Predictions at landmark 10 and 12 for horizon 3, 5 and 10 for two subjects :
dynpred.m32 <- dynpred(m32,landmark=c(10,12),horizon=c(3,5,10),var.time="Time",
fun.time=function(x){10*x},newdata=data_lcmm[4:8,],draws=TRUE,ndraws=2000)

## Plot of the predictions at landmark 10 for horizon 3,5,10 :
plot(dynpred.m32,landmark=10)

## Plot of the predictions at landmark 10 and 12 for horizon 3 :
plot(dynpred.m32,horizon=3)

## End(Not run)

Plot of predicted trajectories and link functions

Description

This function provides the class-specific predicted trajectories stemmed from a hlme, lcmm, multlcmm or Jointlcmm object.

Usage

## S3 method for class 'predictL'
plot(x, legend.loc = "topright", legend, add = FALSE, shades = FALSE, ...)

## S3 method for class 'predictY'
plot(
  x,
  outcome = 1,
  legend.loc = "topright",
  legend,
  add = FALSE,
  shades = FALSE,
  ...
)

## S3 method for class 'predictYcond'
plot(x, legend.loc = "topleft", legend, add = FALSE, shades = TRUE, ...)

Arguments

Author(s)

Cecile Proust-Lima, Benoit Liquet and Viviane Philipps

See Also

hlme, lcmm, Jointlcmm, multlcmm

Examples

################# Prediction from linear latent class model
## fitted model
m<-lcmm(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,
subject='ID',ng=2,data=data_hlme,B=c(0.41,0.55,-0.18,-0.41,
-14.26,-0.34,1.33,13.51,24.65,2.98,1.18,26.26,0.97))
## newdata for predictions plot
newdata<-data.frame(Time=seq(0,5,length=100),
X1=rep(0,100),X2=rep(0,100),X3=rep(0,100))
plot(predictL(m,newdata,var.time="Time"),legend.loc="right",bty="l")
## data from the first subject for predictions plot
firstdata<-data_hlme[1:3,]
plot(predictL(m,firstdata,var.time="Time"),legend.loc="right",bty="l")

 ## Not run: 
################# Prediction from a joint latent class model
## fitted model - see help of Jointlcmm function for details on the model
m3 <- Jointlcmm(fixed= Ydep1~Time*X1,mixture=~Time,random=~Time,
classmb=~X3,subject='ID',survival = Surv(Tevent,Event)~X1+mixture(X2),
hazard="3-quant-splines",hazardtype="PH",ng=3,data=data_lcmm,
B=c(0.7576, 0.4095, -0.8232, -0.2737, 0, 0, 0, 0.2838, -0.6338, 
2.6324, 5.3963, -0.0273, 1.398, 0.8168, -15.041, 10.164, 10.2394, 
11.5109, -2.6219, -0.4553, -0.6055, 1.473, -0.0383, 0.8512, 0.0389, 
0.2624, 1.4982))
# class-specific predicted trajectories 
#(with characteristics of subject ID=193)
data <- data_lcmm[data_lcmm$ID==193,]
plot(predictY(m3,newdata=data,var.time="Time"),bty="l")

## End(Not run)

Posterior classification stemmed from a hlme, lcmm, multlcmm or Jointlcmm estimation

Description

This function provides informations about the posterior classification stemmed from a hlme, lcmm, multlcmm, Jointlcmm, mpjlcmm, externSurv or externX object.

Usage

postprob(x, threshold = c(0.7, 0.8, 0.9), ...)

Arguments

Details

This function provides the number of subjects classified a posteriori in each latent class, the percentage of subjects classified with a posterior probability above a certain threshold, and the classification table that contains the mean of the posterior probability of belonging to each latent class over the subjects classified in each of the latent classes. This table aims at evaluating the quality of the posterior classification. For hlme, lcmm objects, the posterior classification and the classification table are derived from the posterior class-membership probabilities given the vector of repeated measures that are contained in pprob output matrix. For a Jointlcmm object, the first posterior classification and the classification table are derived from the posterior class-membership probabilities given the vector of repeated measures and the time-to-event information (that are contained in columns probYT1, probYT2, etc in pprob output matrix). The second posterior classification is derived from the posterior class-membership probabilities given only the vector of repeated measures (that are contained in columns probY1, probY2, etc in pprob output matrix).

Value

A list containing the posterior classification, the posterior classification table and the percentage of subjects classified with a posterior probability above the given thresholds.

Note

This function can only be used with latent class mixed models and joint latent class mixed models that include at least 2 latent classes

Author(s)

Cecile Proust-Lima, Benoit Liquet and Viviane Philipps

See Also

Jointlcmm, lcmm, hlme, plot.lcmm

Examples

m<-lcmm(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,
subject='ID',ng=2,data=data_hlme,B=c(0.41,0.55,-0.18,-0.41,
-14.26,-0.34,1.33,13.51,24.65,2.98,1.18,26.26,0.97))
postprob(m)

Posterior classification and class-membership probabilities

Description

This function provides the posterior classification and posterior individual class-membership probabilities for external data.

Usage

predictClass(model, newdata, subject = NULL)

Arguments

Value

a matrix with 2+ng columns: the grouping structure, the predicted class and the ng posterior class-membership probabilities.

Author(s)

Sasha Cuau, Viviane Philipps, Cecile Proust-Lima

Examples

## Not run: 
library(NormPsy)
paquid$normMMSE <- normMMSE(paquid$MMSE)
paquid$age65 <- (paquid$age - 65)/10
m2b <- hlme(normMMSE ~ age65+I(age65^2)+CEP, random =~ age65+I(age65^2), subject = 'ID',
data = paquid, ng = 2, mixture =~ age65+I(age65^2), B = c(0, 60, 40, 0, -4, 0, -10, 10,
212.869397, -216.421323,456.229910, 55.713775, -145.715516, 59.351000, 10.072221))
predictClass(m2b, newdata=paquid[1:6,])

## End(Not run)

Class-specific marginal predictions in the latent process scale for lcmm, Jointlcmm and multlcmm objects

Description

This function provides a matrix containing the class-specific predicted trajectories computed in the latent process scale, that is the latent process underlying the curvilinear outcome(s), for a profile of covariates specified by the user. This function applies only to lcmm and multlcmm objects. The function plot.predict provides directly the plot of these class-specific predicted trajectories. The function predictY provides the class-specific predicted trajectories computed in the natural scale of the outcome(s).

Usage

predictL(x, newdata, var.time, na.action = 1, confint = FALSE, ...)

Arguments

Value

An object of class predictL with values :

- pred : a matrix containing the class-specific predicted values in the latent process scale, the lower and the upper limits of the confidence intervals (if calculated).

- times : the var.time variable from newdata

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

plot.predict, predictY, lcmm

Examples

#### Prediction from a 2-class model with a Splines link function
## Not run: 
## fitted model
m<-lcmm(Ydep2~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,
subject='ID',ng=2,data=data_lcmm,link="splines",B=c(
-0.175,      -0.191,       0.654,      -0.443, 
-0.345,      -1.780,       0.913,       0.016, 
 0.389,       0.028,       0.083,      -7.349, 
 0.722,       0.770,       1.376,       1.653, 
 1.640,       1.285))
summary(m)
## predictions for times from 0 to 5 for X1=0
newdata<-data.frame(Time=seq(0,5,length=100),
X1=rep(0,100),X2=rep(0,100),X3=rep(0,100))
predictL(m,newdata,var.time="Time")
## predictions for times from 0 to 5 for X1=1
newdata$X1 <- 1
predictY(m,newdata,var.time="Time")

## End(Not run)

Predictions of the random-effects

Description

The function computes the predicted values of the random effects given observed data provided in input. With multiple latent classes, these predictions are averaged over classes using the posterior class-membership probabilities.

Usage

predictRE(model, newdata, subject = NULL)

Arguments

Value

a matrix containing the grouping structure and the predicted random-effects.

Author(s)

Sasha Cuau, Viviane Philipps, Cecile Proust-Lima

Examples

## Not run: 
library(NormPsy)
paquid$normMMSE <- normMMSE(paquid$MMSE)
paquid$age65 <- (paquid$age - 65)/10
m2b <- hlme(normMMSE ~ age65+I(age65^2)+CEP, random =~ age65+I(age65^2), subject = 'ID',
data = paquid, ng = 2, mixture =~ age65+I(age65^2), B = c(0, 60, 40, 0, -4, 0, -10, 10,
212.869397, -216.421323,456.229910, 55.713775, -145.715516, 59.351000, 10.072221))
predictRE(m2b,newdata=paquid[1:6,])

## End(Not run)

Predictions (marginal and possibly subject-specific in some cases) of a hlme, lcmm, multlcmm or Jointlcmm object in the natural scale of the longitudinal outcome(s) computed from a profile of covariates (marginal) or individual data (subject specific in case of hlme).

Description

For hlme and Jointlcmm objects, the function computes the predicted values of the longitudinal marker (in each latent class of ng>1) for a specified profile of covariates. For lcmm and multlcmm objects, the function computes predicted values in the natural scale of the outcomes for a specified profile of covariates. For linear and threshold links, the predicted values are computed analytically. For splines and Beta links, a Gauss-Hermite or Monte-Carlo integration are used to numerically compute the predictions. In addition, for any type of link function, confidence bands (and median) can be computed by a Monte Carlo approximation of the posterior distribution of the predicted values.

Usage

## S3 method for class 'Jointlcmm'
predictY(
  x,
  newdata,
  var.time,
  methInteg = 0,
  nsim = 20,
  draws = FALSE,
  ndraws = 2000,
  na.action = 1,
  ...
)

## S3 method for class 'hlme'
predictY(
  x,
  newdata,
  var.time,
  draws = FALSE,
  na.action = 1,
  marg = TRUE,
  subject = NULL,
  ...
)

## S3 method for class 'lcmm'
predictY(
  x,
  newdata,
  var.time,
  methInteg = 0,
  nsim = 20,
  draws = FALSE,
  ndraws = 2000,
  na.action = 1,
  ...
)

predictY(x, newdata, var.time, ...)

## S3 method for class 'multlcmm'
predictY(
  x,
  newdata,
  var.time,
  methInteg = 0,
  nsim = 20,
  draws = FALSE,
  ndraws = 2000,
  na.action = 1,
  ...
)

Arguments

Value

An object of class predictY with values :

- pred : a matrix with the same rows (number and order) as in newdata.

For hlme objects and lcmm or Jointlcmm with draws=FALSE, returns a matrix with ng columns corresponding to the ng class-specific vectors of predicted values computed at the point estimate

For objects of class lcmm or Jointlcmm with draws=TRUE, returns a matrix with ng*3 columns representing the ng class-specific 50%, 2.5% and 97.5% percentiles of the approximated posterior distribution of the class-specific predicted values.

For objects of class multlcmm with draws=FALSE, returns a matrix with ng+1 columns: the first column indicates the name of the outcome which is predicted and the ng subsequent columns correspond to the ng class-specific vectors of predicted values computed at the point estimate

For objects of class multlcmm with draws=TRUE, returns a matrix with ng*3+1 columns: the first column indicates the name of the outcome which is predicted and the ng*3 subsequent columns correspond to the ng class-specific 50%, 2.5% and 97.5% percentiles of the approximated posterior distribution of the class-specific predicted values.

For objects of class hlme with marg=FALSE, returns a matrix with 2+ng columns : the grouping structure, subject-specific predictions (pred_ss) averaged over classes and the class-specific subject-specific predictions (with the number of the latent class: pred_ss_1,pred_ss_2,...)

- times : the var.time variable from newdata

Author(s)

Cecile Proust-Lima, Viviane Philipps, Sasha Cuau

See Also

lcmm, multlcmm, hlme, Jointlcmm

Examples

#### Prediction from a 2-class model with a Splines link function
## Not run: 
## fitted model
m<-lcmm(Ydep2~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,
subject='ID',ng=2,data=data_lcmm,link="splines",B=c(
-0.175,      -0.191,       0.654,      -0.443, 
-0.345,      -1.780,       0.913,       0.016, 
 0.389,       0.028,       0.083,      -7.349, 
 0.722,       0.770,       1.376,       1.653, 
 1.640,       1.285))
summary(m)
## predictions for times from 0 to 5 for X1=0
newdata<-data.frame(Time=seq(0,5,length=100),
X1=rep(0,100),X2=rep(0,100),X3=rep(0,100))
pred0 <- predictY(m,newdata,var.time="Time")
head(pred0)
## Option draws=TRUE to compute a MonteCarlo 
# approximation of the predicted value distribution 
# (quite long with ndraws=2000 by default)
\dontrun{
pred0MC <- predictY(m,newdata,draws=TRUE,var.time="Time")
}
## predictions for times from 0 to 5 for X1=1
newdata$X1 <- 1
pred1 <- predictY(m,newdata,var.time="Time")
## Option draws=TRUE to compute a MonteCarlo 
# approximation of the predicted value distribution 
# (quite long with ndraws=2000 by default)
\dontrun{
pred1MC <- predictY(m,newdata,draws=TRUE,var.time="Time")
}

## End(Not run)

Marginal predictions in the natural scale of a pre-transformed outcome

Description

The function computes the predicted values of the longitudinal marker (in each latent class if ng>1) for a specified profile of covariates, when a non-parameterized pre-transformation was applied (e.g., log, square root). A Gauss-Hermite or Monte-Carlo integration is used to numerically compute the back-transformed predictions.

Usage

predictYback(
  x,
  newdata,
  var.time,
  methInteg = 0,
  nsim = 20,
  draws = FALSE,
  ndraws = 2000,
  na.action = 1,
  back,
  ...
)

Arguments

Value

An object of class predictY.

Examples

data_lcmm$transfYdep2 <- sqrt(30 - data_lcmm$Ydep2)
m1 <- hlme(transfYdep2 ~ Time, random=~ Time, subject="ID", data = data_lcmm)
pred1 <- predictYback(m1, newdata = data.frame(Time = seq(0, 3, 0.1)), 
var.time = "Time", back = function(x) {30 - x^2})
plot(pred1)

Conditional predictions of a lcmm, multlcmm or Jointlcmm object in the natural scale of the longitudinal outcome(s) for specified latent process values.

Description

The function computes the predicted values of the longitudinal markers in their natural scale for specified values of the latent process. For splines and Beta links, a Gauss-Hermite integration is used to numerically compute the predictions. In addition, for any type of link function, confidence bands (and median) can be computed by a Monte Carlo approximation of the posterior distribution of the predicted values.

Usage

predictYcond(
  x,
  lprocess,
  condRE_Y = FALSE,
  nsim = 200,
  draws = FALSE,
  ndraws = 2000,
  ...
)

Arguments

Value

An object of class predictYcond with values :

- pred : If draws=FALSE, returns a matrix with 3 columns : the first column indicates the name of the outcome, the second indicates the latent process value and the last is the computed prediction. If draws=TRUE, returns a matrix with 5 columns : the name of the outcome, the latent process value and the 50%, 2.5% and 97.5% percentiles of the approximated posterior distribution of predicted values.

- object : the model from which the predictions are computed.

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

predictY, predictlink

Examples

## Not run: 
m12 <- lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="3-equi-splines")
predm12 <- predictYcond(m12,lprocess=seq(-8,2,length.out=100),draws=TRUE)
plot(predm12)

## End(Not run)

Confidence intervals for the estimated link functions from lcmm, Jointlcmm and multlcmm

Description

This function provides 95% confidence intervals around the estimated transformation given in estimlink attribute of lcmm, Jointlcmm and multlcmm objects. It can also be used to evaluate the link functions at other values than those given in attribute estimlink of lcmm, Jointlcmm or multlcmm object.

Usage

predictlink(x, ndraws, Yvalues, ...)

Arguments

Value

An object of class predictlink with values :

- pred :

For a lcmm or Jointlcmm object, a data frame containing the values at which the transformation is evaluated, the transformed values and the lower and the upper limits of the confidence intervals (if ndraws>0).

For a multlcmm object, a data frame containing the indicator of the outcome, the values at which the transformations are evaluated,the transformed values and the lower and the upper limits of the confidence intervals (if ndraws>0).

- object : the object from which the link function is predicted

Author(s)

Cecile Proust-Lima and Viviane Philipps

See Also

lcmm, multlcmm, plot.lcmm, plot.predictlink

Examples

 ## Not run: 

## Univariate mixed model with splines link funciton
m14<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_lcmm,link="5-manual-splines",intnodes=c(10,20,25),
B=c(-0.89255, -0.09715, 0.56335, 0.21967, 0.61937, -7.90261, 0.75149, 
-1.22357, 1.55832, 1.75324, 1.33834, 1.0968))

##Transformed values of several scores and their confidence intervals
transf.m14 <- predictlink(m14,ndraws=2000,Yvalues=c(0,1,7:30))
plot(transf.m14)


## Multivariate mixed model with splines link functions
m1 <- multlcmm(Ydep1+Ydep2~1+Time*X2+contrast(X2),random=~1+Time,
subject="ID",randomY=TRUE,link=c("4-manual-splines","3-manual-splines"),
intnodes=c(8,12,25),data=data_lcmm,
B=c(-1.071, -0.192,  0.106, -0.005, -0.193,  1.012,  0.870,  0.881,
  0.000,  0.000, -7.520,  1.401,  1.607 , 1.908,  1.431,  1.082,
 -7.528,  1.135 , 1.454 , 2.328, 1.052))
##Confidence intervals for the transformed values (given in m1$estimlink)
transf.m1 <- predictlink(m1,ndraws=200)
plot(transf.m1)

## End(Not run)

Brief summary of a hlme, lcmm, Jointlcmm,multlcmm, epoce or Diffepoce objects

Description

The function provides a brief summary of hlme, lcmm,multlcmm or Jointlcmm estimations, and epoce or Diffepoce computations.

Usage

## S3 method for class 'lcmm'
print(x, ...)

Arguments

Author(s)

Cecile Proust-Lima, Viviane Philipps, Amadou Diakite and Benoit Liquet

See Also

hlme, lcmm, Jointlcmm, epoce, Diffepoce

Simulated dataset simdataHADS

Description

The data mimic the PREDIALA study described and analyzed in Proust-Lima et al (2021 - https://arxiv.org/abs/2109.13064). The study aims to describe the trajectories of depressive symptomatology of patients suffering end-stage renal disease and registered on the renal transplant waiting list. Repeated measures of anxiety and depression (HADS) were simulated at different times of measurement for 561 subjects. Four time-independent covariates were also generated: group (dialyzed or pre-emptive), sex and age at entry in the cohort and time on the waiting list at entry in the cohort.

Format

A data frame with 1140 observations on the following 13 variables.

grp

group with 0=dialyzed and 1=preemptive

sex

sex with 0=woman and 1=man

age

age at entry in the cohort

hads_2

item 2 of HADS measuring depression

hads_4

item 4 of HADS measuring depression

hads_6

item 6 of HADS measuring depression

hads_8

item 8 of HADS measuring depression

hads_10

item 10 of HADS measuring depression

hads_12

item 12 of HADS measuring depression

hads_14

item 14 of HADS measuring depression

ID

subject identification number

time

time of measurement

time_entry

time on the waiting list at entry in the cohort

Data simulation according to models from lcmm package

Description

This function simulates a sample according to a model estimated with hlme, lcmm, multlcmm or Jointlcmm functions.

Usage

## S3 method for class 'lcmm'
simulate(
  object,
  nsim,
  seed,
  times,
  tname = NULL,
  n,
  Xbin = NULL,
  Xcont = NULL,
  entry = 0,
  dropout = NULL,
  pMCAR = 0,
  ...
)

Arguments

Value

a data frame with one line per observation and one column per variable. Variables appears in the following order : subject id, measurement time, entry time, binary covariates, continuous covariates, longitudinal outcomes, latent class, entry time, survival time, event indicator.

Author(s)

Viviane Philipps and Cecile Proust-Lima

Examples

## estimation of a 2 classes mixed model
m2 <- hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',
         ng=2,data=data_hlme,B=c(0.11,-0.74,-0.07,20.71,
                                 29.39,-1,0.13,2.45,-0.29,4.5,0.36,0.79,0.97))

## simulate according to model m2 with same number of subjects and
## same measurement times as in data_lcmm. Binary covariates X1 and X2 are simulated
## according to a Bernoulli distribution with probability p=0.5, continuous covariate
## X3 is simulated according to a Gaussian distribution with mean=1 and sd=1 :
dsim1 <- simulate(m2, times=data_hlme[,c("ID","Time")],
                  Xbin=list(X1=0.5, X2=0.5), Xcont=list(X3=c(1,1)))

## simulate a dataset of 300 subjects according to the same model
## with new observation times, equally spaced and ranging from 0 to 3 :
dsim2 <- simulate(m2, times=c(0,3,0.5,0), n=300, tname="Time",
                  Xbin=list(X1=0.5, X2=0.5), Xcont=list(X3=c(1,1)))

Summary of a hlme, lcmm, Jointlcmm, multlcmm, mpjlcmm, externSurv, externX epoce or Diffepoce objects

Description

The function provides a summary of hlme, lcmm, multlcmm and Jointlcmm estimations, or epoce and Diffepoce computations.

Usage

## S3 method for class 'lcmm'
summary(object, ...)

Arguments

Value

For epoce or Diffepoce objects, returns NULL. For hlme, lcmm, Jointlcmm or multlcmm returns also a matrix containing the fixed effect estimates in the longitudinal model, their standard errors, Wald statistics and p-values

Author(s)

Cecile Proust-Lima, Viviane Philipps, Amadou Diakite and Benoit Liquet

See Also

hlme, lcmm, multlcmm, Jointlcmm, epoce, Diffepoce

Summary of models

Description

This function provides a plot summarizing the results of different models fitted by hlme, lcmm, multlcmm, Jointlcmm, mpjlcmm or externVar.

Usage

summaryplot(
  m1,
  ...,
  which = c("BIC", "entropy", "ICL"),
  mfrow = c(1, length(which)),
  xaxis = "G"
)

Arguments

Details

Can be reported the usual criteria used to assess the fit and the clustering of the data: - maximum log-likelihood L (the higher the better) - number of parameters P, number of classes G, convergence criterion (1 = converged) - AIC (the lower the better) computed as -2L+2P - BIC (the lower the better) computed as -2L+ P log(N) where N is the number of subjects - SABIC (the lower the better) computed as -2L+ P log((N+2)/24) - Entropy (the closer to one the better) computed as 1+sum[pi_ig*log(pi_ig)]/(N*log(G)) where pi_ig is the posterior probability that subject i belongs to class g - ICL (the lower the better) computed in two ways : ICL1 = BIC - sum[pi_ig*log(pi_ig)] or ICL2 = BIC - 2*sum(log(max(pi_ig)), where the max is taken over the classes for each subject. - %class computed as the proportion of each class based on c_ig

Author(s)

Sasha Cuau, Viviane Philipps, Cecile Proust-Lima

See Also

summary, summarytable

Examples

## Not run: 
library(NormPsy)
paquid$normMMSE <- normMMSE(paquid$MMSE)
paquid$age65 <- (paquid$age - 65)/10
m1 <- hlme(normMMSE~age65+I(age65^2)+CEP, random=~age65+I(age65^2), subject='ID', data=paquid)
m2 <- hlme(normMMSE~age65+I(age65^2)+CEP, random=~age65+I(age65^2), subject='ID', data=paquid,
ng = 2, mixture=~age65+I(age65^2), B=m1)
m3g <- gridsearch(hlme(normMMSE~age65+I(age65^2)+CEP, random=~age65+I(age65^2), subject='ID',
data=paquid, ng=3, mixture=~age65+I(age65^2)), rep=100, maxiter=30, minit=m1)
summaryplot(m1, m2, m3g, which=c("BIC","entropy","ICL"),bty="l",pch=20,col=2)

## End(Not run)

Summary of models

Description

This function provides a table summarizing the results of different models fitted by hlme, lcmm, multlcmm, Jointlcmm, mpjlcmm or externVar.

Usage

summarytable(
  m1,
  ...,
  which = c("G", "loglik", "npm", "BIC", "%class"),
  display = TRUE
)

Arguments

Details

Can be reported the usual criteria used to assess the fit and the clustering of the data: - maximum log-likelihood L (the higher the better) - number of parameters P, number of classes G, convergence criterion (1 = converged) - AIC (the lower the better) computed as -2L+2P - BIC (the lower the better) computed as -2L+ P log(N) where N is the number of subjects - SABIC (the lower the better) computed as -2L+ P log((N+2)/24) - Entropy (the closer to one the better) computed as 1+sum[pi_ig*log(pi_ig)]/(N*log(G)) where pi_ig is the posterior probability that subject i belongs to class g - ICL (the lower the better) computed in two ways : ICL1 = BIC - sum[pi_ig*log(pi_ig)] or ICL2 = BIC - 2*sum(log(max(pi_ig)), where the max is taken over the classes for each subject. - %class computed as the proportion of each class based on c_ig

Value

a matrix giving for each model the values of the requested indexes. By default, the number a latent classes, the log-likelihood, the number of parameters, the BIC and the posterior probability of the latent classes.

Author(s)

Cecile Proust-Lima, Viviane Philipps

See Also

summary, hlme, lcmm, multlcmm, Jointlcmm

Update the longitudinal submodels

Description

This function updates the longitudinal submodels of a mpjlcmm object by injecting the estimated parameters and their variances in each hlme/lcmm/multlcmm model used to define the multi-process joint model. The same (uni-process) models as specified in the mpjlcmm call are returned, with updated outputs for best, V, conv, cholesky, pred, predRE, predRE_Y, pprob. All postfit functions (plots, predictions, ...) can then be called on the updated hlme/lcmm/multlcmm models. See mpjlcmm's help page for examples.

Usage

## S3 method for class 'mpjlcmm'
update(object, ...)

Arguments

Value

A list of hlme/lcmm/multlcmm models. The models appear in the same order as specified in the call to the mpjlcmm function.

Cross classifications

Description

This function crosses the posterior classifications of two estimated models

Usage

xclass(m1, m2)

Arguments

Value

the contingency table of the two classifications

Author(s)

Viviane Philipps and Cecile Proust-Lima

Examples

## Estimation of the models
m2 <- hlme(Y~Time*X1,mixture=~Time,random=~Time,classmb=~X2+X3,subject='ID',ng=2,
data=data_hlme,B=c(0.11,-0.74,-0.07,20.71,29.39,-1,0.13,2.45,-0.29,4.5,0.36,0.79,0.97))
m3 <- hlme(fixed = Y ~ Time * X1, mixture = ~Time, random = ~Time,subject = "ID",
classmb = ~X2 + X3, ng = 3, data = data_hlme,B=c(-0.21, 0.31, -2.11, -0.81, -0.24,
-0.18, 25.4, 20.09, 30.18, -0.43, -1.1, 0.25, 2.37, -0.29, 2.34, 0.03, 0.74, 0.97))

## Compare the classifications
xclass(m2,m3)
# The 39 subjects in class 2 of m3 come from class 1 of m2.
# In the same way, all the subjects in class 3 come from class 2 of m2.
# Class 1 of m3 mixes subject from class 1 and class 2 of m2.