Description Usage Arguments Details Value Author(s) References See Also Examples
This function fits joint latent class mixed models for a longitudinal
outcome and a rightcensored (possibly lefttruncated) timetoevent. 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
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75  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 = 1e04,
convL = 1e04,
convG = 1e04,
maxiter = 100,
nsim = 100,
prior = NULL,
logscale = FALSE,
subset = NULL,
na.action = 1,
posfix = NULL,
partialH = FALSE,
verbose = TRUE,
returndata = FALSE,
var.time = 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 = 1e04,
convL = 1e04,
convG = 1e04,
maxiter = 100,
nsim = 100,
prior = NULL,
logscale = FALSE,
subset = NULL,
na.action = 1,
posfix = NULL,
partialH = FALSE,
verbose = TRUE,
returndata = FALSE,
var.time = NULL
)

fixed 
twosided linear formula object for the fixedeffects in the
linear mixed model. The response outcome is on the left of 
mixture 
onesided formula object for the classspecific fixed effects
in the linear mixed model (to specify only for a number of latent classes
greater than 1). Among the list of covariates included in 
random 
optional onesided formula for the randomeffects in the
linear mixed model. Covariates with a randomeffect are separated by

subject 
name of the covariate representing the grouping structure (called subject identifier) specified with ”. 
classmb 
optional onesided formula describing the covariates in the
classmembership multinomial logistic model. Covariates included are
separated by 
ng 
optional number of latent classes considered. If 
idiag 
optional logical for the structure of the variancecovariance
matrix of the randomeffects. If 
nwg 
optional logical indicating if the variancecovariance of the
randomeffects is classspecific. If 
survival 
twosided formula object. The left side of the formula
corresponds to a 
hazard 
optional family of hazard function assumed for the survival
model. By default, "Weibull" specifies a Weibull baseline risk function.
Other possibilities are "piecewise" for a piecewise constant risk function
or "splines" for a cubic Msplines baseline risk function. For these two
latter families, the number of nodes and the location of the nodes should be
specified as well, separated by 
hazardtype 
optional indicator for the type of baseline risk function when ng>1. By default "Specific" indicates a classspecific baseline risk function. Other possibilities are "PH" for a baseline risk function proportional in each latent class, and "Common" for a baseline risk function that is common over classes. In the presence of competing events, a vector of hazardtypes should be given. 
hazardnodes 
optional vector containing interior nodes if

hazardrange 
optional vector indicating the range of the survival times (that is the minimum and maximum). By default, the range is defined according to the minimum and maximum observed values of the survival times. The option should be used only for piecewise constant and Splines hazard functions. 
TimeDepVar 
optional vector containing an intermediate time corresponding to a change in the risk of event. This timedependent covariate can only take the form of a time variable with the assumption that there is no effect on the risk before this time and a constant effect on the risk of event after this time (example: initiation of a treatment to account for). 
link 
optional family of link functions to estimate. By default,
"linear" option specifies a linear link function leading to a standard
linear mixed model (homogeneous or heterogeneous as estimated in

intnodes 
optional vector of interior nodes. This argument is only required for a Isplines link function with nodes entered manually. 
epsY 
optional definite positive real used to rescale the marker in (0,1) when the beta link function is used. By default, epsY=0.5. 
range 
optional vector indicating the range of the outcome (that is the minimum and maximum). By default, the range is defined according to the minimum and maximum observed values of the outcome. The option should be used only for Beta and Splines transformations. 
cor 
optional brownian motion or autoregressive process modeling the correlation between the observations. "BM" or "AR" should be specified, followed by the time variable between brackets. By default, no correlation is added. 
data 
optional data frame containing the variables named in

B 
optional specification for the initial values for the parameters.
Three options are allowed: (1) a vector of initial values is entered (the
order in which the parameters are included is detailed in 
convB 
optional threshold for the convergence criterion based on the parameter stability. By default, convB=0.0001. 
convL 
optional threshold for the convergence criterion based on the loglikelihood stability. By default, convL=0.0001. 
convG 
optional threshold for the convergence criterion based on the derivatives. By default, convG=0.0001. 
maxiter 
optional maximum number of iterations for the Marquardt iterative algorithm. By default, maxiter=150. 
nsim 
optional number of points for the predicted survival curves and predicted baseline risk curves. By default, nsim=100. 
prior 
optional name of a covariate containing a prior information about the latent class membership. The covariate should be an integer with values in 0,1,...,ng. Value O indicates no prior for the subject while a value in 1,...,ng indicates that the subject belongs to the corresponding latent class. 
logscale 
optional boolean indicating whether an exponential (logscale=TRUE) or a square (logscale=FALSE by default) transformation is used to ensure positivity of parameters in the baseline risk functions. See details section 
subset 
a specification of the rows to be used: defaults to all rows. This can be any valid indexing vector for the rows of data or if that is not supplied, a data frame made up of the variable used in formula. 
na.action 
Integer indicating how NAs are managed. The default is 1 for 'na.omit'. The alternative is 2 for 'na.fail'. Other options such as 'na.pass' or 'na.exclude' are not implemented in the current version. 
posfix 
Optional vector specifying the indices in vector B of the parameters that should not be estimated. Default to NULL, all parameters are estimated. 
partialH 
optional logical for Piecewise and Splines baseline risk functions and Splines link functions only. Indicates whether the parameters of the baseline risk or link functions can be dropped from the Hessian matrix to define convergence criteria. 
verbose 
logical indicating if information about computation should be reported. Default to TRUE. 
returndata 
logical indicating if data used for computation should be returned. Default to FALSE, data are not returned. 
var.time 
optional character indicating the name of the time variable. 
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^21) 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), nz1 parameters are necesssary p_1,...p_nz1 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 Msplines.
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) ng1 parameters are required for intercepts in the latent class
membership model, and if covariates are included in classmb
, ng1
parameters should be entered for each one; (2) parameters for the baseline
risk function: 2 parameters for each Weibull, nz1 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, ng1 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 causespecific 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 randomeffect specified in
random
(including the intercept) if idiag=TRUE
and the
inferior triangular variancecovariance matrix of all the randomeffects if
idiag=FALSE
; (6) only if nwg=TRUE
, ng1 parameters for
classspecific proportional coefficients for the variance covariance matrix
of the randomeffects; (7) the variance of the residual error.
C. CAUTION
Some caution should be made when using the program:
(1) As the loglikelihood 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 loglikelihood in addition to the parameter and loglikelihood 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).
The list returned is:
loglik 
loglikelihood of the model 
best 
vector of parameter estimates in the same order as specified in

V 
vector containing
the upper triangle matrix of variancecovariance estimates of 
gconv 
vector of convergence criteria: 1. on the parameters, 2. on the likelihood, 3. on the derivatives 
conv 
status of convergence: =1 if the convergence criteria were satisfied, =2 if the maximum number of iterations was reached, =4 or 5 if a problem occured during optimisation 
call 
the matched call 
niter 
number of Marquardt iterations 
pred 
table of
individual predictions and residuals; it includes marginal predictions
(pred_m), marginal residuals (resid_m), subjectspecific predictions
(pred_ss) and subjectspecific residuals (resid_ss) averaged over classes,
the observation (obs) and finally the classspecific marginal and
subjectspecific predictions (with the number of the latent class:
pred_m_1,pred_m_2,...,pred_ss_1,pred_ss_2,...). If 
pprob 
table of posterior classification and posterior individual classmembership probabilities based on the longitudinal data and the timetoevent data 
pprobY 
table of posterior classification and posterior individual classmembership probabilities based only on the longitudinal data 
predRE 
table containing individual predictions of the randomeffects: a column per randomeffect, a line per subject 
cholesky 
vector containing the estimates of the Cholesky transformed parameters of the variancecovariance matrix of the randomeffects 
scoretest 
Statistic of the Score Test for the conditional independence assumption of the longitudinal and survival data given the latent class structure. Under the null hypothesis, the statistics is a Chisquare with p degrees of freedom where p indicates the number of randomeffects in the longitudinal mixed model. See JacqminGadda and ProustLima (2009) for more details. 
predSurv 
table of predictions giving for the window of times to event (called "time"), the predicted baseline risk function in each latent class (called "RiskFct") and the predicted cumulative baseline risk function in each latent class (called "CumRiskFct"). 
hazard 
internal information about the hazard specification used in related functions 
data 
the original data set (if returndata is TRUE) 
Cecile Proust Lima, Amadou Diakite and Viviane Philipps
ProustLima 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), 156. 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 prostatespecific antigen readings and prostate cancer. Journal of the American Statistical Association 97, 5365.
ProustLima, C. and Taylor, J. (2009). Development and validation of a dynamic prognostic tool for prostate cancer recurrence using repeated measures of posttreatment PSA: a joint modelling approach. Biostatistics 10, 53549.
JacqminGadda, H. and ProustLima, C. (2010). Score test for conditional independence between longitudinal outcome and timetoevent given the classes in the joint latent class model. Biometrics 66(1), 119
ProustLima, Sene, Taylor and JacqminGadda (2014). Joint latent class models of longitudinal and timetoevent data: a review. Statistical Methods in Medical Research 23, 7490.
postprob
, plot.Jointlcmm
,
plot.predict
, epoce
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86  #### Example of a joint latent class model estimated for a varying number
# of latent classes:
# The linear mixed model includes a subject (ID) and classspecific
# 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 Msplines 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 classspecific 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.
## Not run:
#### 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="3quantsplines"
,hazardtype="PH",ng=1,data=data_lcmm)
summary(m1)
#Goodnessoffit statistics for m1:
# maximum loglikelihood: 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="3quantsplines",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)
#Goodnessoffit statistics for m2:
# maximum loglikelihood: 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="3quantsplines",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)
#Goodnessoffit statistics for m3:
# maximum loglikelihood: 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="3quantsplines",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)
#Goodnessoffit statistics for m4:
# maximum loglikelihood: 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")
# Classspecific predicted baseline risk & survival functions in the
# 3class 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")
# classspecific predicted trajectories in the 3class 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)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.