Jointlcmm: 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. In this first version, only the Gaussian longitudinal outcome is handled.

Usage

Jointlcmm(fixed, mixture, random, subject, classmb, ng=1, 
 idiag=FALSE, nwg=FALSE, survival,hazard="Weibull",
 hazardtype="Specific", hazardnodes=NULL,TimeDepVar=NULL,
 cor=NULL,data, B, convB=1e-4, convL=1e-4, convG=1e-4,
 maxiter=150, nsim=100, prior,logscale=FALSE, 
 subset=NULL, na.action=1)

Arguments

fixed

two-sided linear formula object for the fixed-effects in the linear mixed model. The response outcome is on the left of ~ and the covariates are separated by + on the right of the ~. By default, an intercept is incl

mixture

one-sided formula object for the class-specific 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 fixed, the covariates with class-specific regre

random

optional one-sided formula for the random-effects in the linear mixed model. Covariates with a random-effect are separated by +. By default, an intercept is included. If no intercept, -1 should be the first term included.

subject

name of the covariate representing the grouping structure (called subject identifier) specified with ''.

classmb

optional one-sided formula describing the covariates in the class-membership multinomial logistic model. Covariates included are separated by +. No intercept should be included in this formula.

optional number of latent classes considered. If ng=1 (by default) no mixture nor classmb should be specified. If ng>1, mixture is required.

idiag

optional logical for the structure of the variance-covariance matrix of the random-effects. If FALSE, a non structured matrix of variance-covariance is considered (by default). If TRUE a diagonal matrix of variance-covariance is

nwg

optional logical indicating if the variance-covariance of the random-effects is class-specific. If FALSE the variance-covariance matrix is common over latent classes (by default). If TRUE a class-specific proportional parameter m

optional vector containing the initial values for the parameters. The order in which the parameters are included is detailed in details section. If no vector is specified, a preliminary analysis involving the estimation of a standard linear

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 log-likelihood 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 be

survival

two-sided formula object. The left side of the formula corresponds to a surv() object of type "counting" for right-censored and left-truncated data (example: Surv(Time,EntryTime,Indicator)) or of type "right" for right-censored d

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 M-splines baseline ris

hazardtype

optional indicator for the type of baseline risk function when ng>1. By default "Specific" indicates a class-specific baseline risk function. Other possibilities are "PH" for a baseline risk function proportional in each latent class, and "Common" for a b

hazardnodes

optional vector containing interior nodes if splines or piecewise is specified for the baseline hazard function in hazard.

TimeDepVar

optional vector containing an intermediate time corresponding to a change in the risk of event. This time-dependent 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

logscale

optional boolean indicating whether an exponential (logscale=TRUE) or a square (logscale=FALSE -by default) transformation is used to ensure positivity of splines/piecewise parameters in the baseline risk functions. See details section

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 fixed, mixture, random, classmb and subject.

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.

Value

The list returned is:
logliklog-likelihood of the model
bestvector of parameter estimates in the same order as specified in B and detailed in section details
Vvector containing the upper triangle matrix of variance-covariance estimates of Best with exception for variance-covariance parameters of the random-effects for which V contains the variance-covariance estimates of the Cholesky transformed parameters displayed in cholesky
gconvvector of convergence criteria: 1. on the parameters, 2. on the likelihood, 3. on the derivatives
convstatus 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
callthe matched call
niternumber of Marquardt iterations
datasetdataset
predtable of individual predictions and residuals; it includes marginal predictions (pred_m), marginal residuals (resid_m), subject-specific predictions (pred_ss) and subject-specific residuals (resid_ss) averaged over classes, the observation (obs) and finally the class-specific marginal and subject-specific predictions (with the number of the latent class: pred_m_1,pred_m_2,...,pred_ss_1,pred_ss_2,...)
pprobtable of posterior classification and posterior individual class-membership probabilities based on the longitudinal data and the time-to-event data
pprobYtable of posterior classification and posterior individual class-membership probabilities based only on the longitudinal data
predREtable containing individual predictions of the random-effects: a column per random-effect, a line per subject
choleskyvector containing the estimates of the Cholesky transformed parameters of the variance-covariance matrix of the random-effects
CIstatStatistic 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 Chi-square with p degrees of freedom where p indicates the number of random-effects in the longitudinal mixed model. See Jacqmin-Gadda and Proust-Lima (2009) for more details.
predSurvtable 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").
hazardinternal information about the hazard specification used in related functions
specifinternal information used in related functions
Namesinternal information used in related fnctions
Names2internal information used in related functions

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) = exp(w_1)*exp(w_2)(t)^(exp(w_2)-1).

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.

For "piecewise" and "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; (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; (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 specification of initial values in B should be preferred.

(4) 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 should be assessed. 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).

References

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.

Examples

Run this code

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

#### data loading
data(data_Jointlcmm)

#### 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_Jointlcmm)
summary(m1)
#Goodness-of-fit statistics for m1:
#    maximum log-likelihood: -3944.77 ; AIC: 7919.54  ;  BIC: 7975.09  

#### 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_Jointlcmm,
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  


#### 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_Jointlcmm,
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_Jointlcmm,
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
postprob(m3)
# Class-specific predicted baseline risk & survival functions in the 
# 3-class model retained (for the reference value of the covariates) 
plot.baselinerisk(m3,bty="l")
plot.baselinerisk(m3,ylim=c(0,5),bty="l")
plot.survival(m3,bty="l")
# class-specific predicted trajectories in the 3-class model retained 
# (with characteristics of subject ID=193)
data <- data_Jointlcmm[data_Jointlcmm$ID==193,]
plot.predict(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 ##############

Run the code above in your browser using DataLab