lcmm (version 1.6.6)

lcmm: 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,range=NULL,
 subset=NULL, na.action=1)

Arguments

fixed

a two-sided linear formula object for specifying the fixed-effects in the linear mixed model at the latent process level. The response outcome is on the left of ~ and the covariates are separated by + on the right of the ~<

mixture

a one-sided formula object for the class-specific fixed effects in the latent process 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-spec

random

an optional one-sided formula for the random-effects in the latent process 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 inclu

subject

name of the covariate representing the grouping structure.

classmb

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

ng

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

idiag

optional logical for the variance-covariance structure 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 considered.

nwg

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

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 hlme). Other possibilities include "beta" for

intnodes

optional vector of interior nodes. This argument is only required for a I-splines 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.

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.

B

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 and ng>1, a preliminary analysis involving the estimation of a mixed m

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=100.

nsim

number of points used to plot the estimated link function. By default, nsim=100.

prior

name of the covariate containing the prior on the latent class membership. The covariate should be an integer with values in 0,1,...,ng. When there is no prior, the value should be 0. When there is a prior for the subject, the value should be the number o

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 transformation

subset

optional vector giving the line numbers of the data frame data to use. By default, all lines.

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:
• ns
  number of grouping units in the dataset
• ng
  number of latent classes
• loglik
  log-likelihood of the model
• best
  vector of parameter estimates in the same order as specified in B and detailed in section details
• V
  vector 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
• 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
• dataset
  dataset
• N
  internal information used in related functions
• idiag
  internal information used in related functions
• pred
  table of individual predictions and residuals in the underlying latent process scale; 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 transformed observations in the latent process scale (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,...). This output is not available yet when specifying a thresholds transformation.
• pprob
  table of posterior classification and posterior individual class-membership probabilities
• Xnames
  list of covariates included in the model
• predRE
  table containing individual predictions of the random-effects : a column per random-effect, a line per subject. This output is not available yet when specifying a thresholds transformation.
• cholesky
  vector containing the estimates of the Cholesky transformed parameters of the variance-covariance matrix of the random-effects
• estimlink
  table containing the simulated values of the marker and corresponding estimated link function
• epsY
  definite positive real used to rescale the marker in (0,1) when the beta link function is used. By default, epsY=0.5.
• linktype
  indicator of link function type: 0 for linear, 1 for beta, 2 for splines and 3 for thresholds
• linknodes
  vector of nodes useful only for the 'splines' link function

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="" the="" number="" of="" parameters="" is="" reduced="" if="" some="" levels="" do="" not="" have="" any="" information.="" for="" example,="" a="" level="" c="" observed="" in="" dataset,="" corresponding="" threshold="" constrained="" to="" be="" same="" as="" previous="" one="" b_{c}.="" link="" function="" by="" 1.="" <="" p="">

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

We understand that it can be difficult to enter the correct number of parameters in B at the first place. So we recommend to run the program without specifying the initial vector B even if this model does not converge. As the final vector best has exactly the same structure as B (even when the program stops without convergence), it will help defining a satisfying vector of initial values B for next runs.

C. CAUTIONS REGARDING THE USE OF THE PROGRAM

Some caution should be made when using the program. First, 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 should be assessed. 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 a reason why specification of initial values in B should be preferred.

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 treshold 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, Porust-Lima, Samieri, Liquet (2012) for further details.

References

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 (2012). A universal approximate cross-validation criterion and its asymptotic distribution, Arxiv.

See Also

postprob, plot.postprob, plot.linkfunction, plot.predict, hlme

Examples

#### 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 --
data(data_Jointlcmm)
# linear link function
m10<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_Jointlcmm,link="linear")
summary(m10)
# Beta link function
m11<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_Jointlcmm,link="beta")
summary(m11)
plot.linkfunction(m11,bty="l")
# I-splines with 3 equidistant nodes
m12<-lcmm(Ydep2~Time+I(Time^2),random=~Time,subject='ID',ng=1,
data=data_Jointlcmm,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_Jointlcmm,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_Jointlcmm,link="5-manual-splines",intnodes=c(10,20,25))
summary(m14)
plot.linkfunction(m14,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_Jointlcmm,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.linkfunction(m15,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.linkfunction(m10,col=1,xlab="latent process",ylab="marker",bty="l",xlim=c(-10,5),legend=NULL)
plot.linkfunction(m11,add=TRUE,col=2,legend=NULL)
plot.linkfunction(m12,add=TRUE,col=3,legend=NULL)
plot.linkfunction(m13,add=TRUE,col=4,legend=NULL)
plot.linkfunction(m14,add=TRUE,col=5,legend=NULL)
plot.linkfunction(m15,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_Jointlcmm,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_Jointlcmm,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_Jointlcmm,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_Jointlcmm[data_Jointlcmm$ID==193,]
plot.predict(m22,var.time="Time",newdata=data,bty="l")