Robust Growth Mixture Models

1. From profiles to trajectories

robust_lpa() finds latent profiles in variables measured once. When the same variables are measured repeatedly, the question often becomes how people change: do they all follow one average trajectory, or are there subgroups with different courses (stable, slowly declining, rapidly declining)? Growth mixture models (GMM; Verbeke & Lesaffre, 1996; Muthen & Shedden, 1999) answer this question: they are finite mixtures of linear mixed-effects models, in which every latent class has its own mean trajectory, and persons deviate from their class trajectory through random effects. Without random effects the model is a latent class growth analysis (LCGA; Nagin, 1999).

robust_gmm() fits these models for one or several outcomes at once, with

unbalanced data: every person contributes the likelihood of the values actually observed, at his/her own times, so persons with missed visits or who dropped out are kept, and the estimates are valid when drop-out depends on earlier observed values (missing at random);
robust classes: multivariate-t classes (robust_method = "t"), whose heavy tails absorb persons with outlying trajectories or gross errors instead of creating spurious classes;
two estimation engines (EM and MCMC); and
LASSO penalties that identify classes that do not change on an outcome and outcomes that do not differ between classes.

2. Example data

neuro_long contains simulated annual assessments (up to six visits) of 400 persons on three tests, in long format. Three latent classes were simulated: a stable class, a slowly declining class and a fast declining class on Memory and Executive; Speed declines slightly and equally in all classes. Persons drop out more often after a low Memory score, and a few scores are corrupted by gross errors (see ?neuro_long).

data(neuro_long)
head(neuro_long)
#>   ID Year Memory Executive Speed   True_Class  Age Biomarker
#> 1  1    0   53.3      53.7  42.2 Slow decline 60.5       924
#> 2  1    1   45.8      51.1  48.1 Slow decline 60.5       924
#> 3  1    2   49.0      44.1  43.4 Slow decline 60.5       924
#> 4  1    3   46.4      46.0  42.0 Slow decline 60.5       924
#> 5  1    4   41.6      40.9  41.1 Slow decline 60.5       924
#> 6  1    5   39.2      42.4  40.6 Slow decline 60.5       924
table(visits = table(neuro_long$ID))
#> visits
#>   1   2   3   4   5   6 
#>  19  20  26  32  29 274

3. Fitting a robust growth mixture model

The data are in long format: one row per person and visit, with the person identifier (id), the time variable (time, here years since baseline, so that the intercept is the baseline level) and the outcomes. By default every class has a linear trajectory (degree = 1), persons have correlated random intercepts and slopes (random = "slope"), and the random-effect covariance and the residual variances are shared by the classes (re_cov = "equal", resid_var = "equal"), the usual and more stable specification.

fit <- robust_gmm(neuro_long, id = "ID", time = "Year",
                  outcomes = c("Memory", "Executive"), G = 3,
                  robust_method = "t", n_starts = 3)
fit
#> <robust_gmm> EM | G = 3 | persons = 400 | outcomes: Memory, Executive
#> Trajectory: degree 1 | random: slope (full, equal across classes) | robust: multivariate t (nu = 10.27)
#> LogLik = -10536.7 | BIC = 21235.3 | Entropy = 0.801
#> Proportions: C1=0.18, C2=0.50, C3=0.32

The estimated trajectories are on the original scale of the outcomes:

summary(fit)
#> robust_gmm summary -- EM | G = 3 | persons = 400
#> Estimation: robust: multivariate t (nu = 10.27) | degree 1 | random: slope
#> 
#> Class trajectories (original scale):
#>  Class   Outcome (Intercept)   Year
#>      1    Memory      44.966 -5.050
#>      1 Executive      42.888 -3.968
#>      2    Memory      53.393 -0.037
#>      2 Executive      52.273 -0.074
#>      3    Memory      48.556 -2.135
#>      3 Executive      47.854 -1.664
#> 
#> Class sizes (modal assignment):
#>  Class   N Proportion
#>      1  78      0.185
#>      2 198      0.499
#>      3 124      0.316
#> 
#> Random-effect covariance (equal across classes; class 1 shown):
#>                       Memory:(Intercept) Memory:Year Executive:(Intercept)
#> Memory:(Intercept)                18.116       0.375                 9.787
#> Memory:Year                        0.375       0.053                 0.090
#> Executive:(Intercept)              9.787       0.090                22.239
#> Executive:Year                    -0.162       0.083                 0.384
#>                       Executive:Year
#> Memory:(Intercept)            -0.162
#> Memory:Year                    0.083
#> Executive:(Intercept)          0.384
#> Executive:Year                 0.224
#> 
#> Residual variances:
#>         Memory Executive
#> Class_1   6.63     5.639
#> Class_2   6.63     5.639
#> Class_3   6.63     5.639
#> 
#> Fit: LogLik = -10536.74 | Parameters = 27 | AIC = 21127.5 | BIC = 21235.3 | SABIC = 21149.6 | Entropy = 0.801

plot_robust_gmm(fit)

Class mean trajectories over the individual trajectories

fit$probabilities and fit$assignments give the posterior class probabilities and the modal class of every person (in the order of fit$ids); fit$weights gives each person’s robustness weight, which is small for the persons whose trajectory is far from every class (here, the persons with a gross error):

head(sort(fit$weights))
#> [1] 0.1016763 0.1031012 0.1155740 0.1194273 0.1233247 0.1260103

4. How many classes? Robust vs. classical estimation

estimate_gmm_robust() fits several numbers of classes at once. Comparing the classical (Gaussian) and the robust (t) models shows why robustness matters here: the classical model uses an extra class to accommodate a handful of persons with gross errors (note its minimum class size), and BIC then favours too many classes; the t model does not.

classical <- estimate_gmm_robust(neuro_long, id = "ID", time = "Year",
                                 outcomes = c("Memory", "Executive"),
                                 n_classes = 2:4, robust = FALSE, n_starts = 2)
robust <- estimate_gmm_robust(neuro_long, id = "ID", time = "Year",
                              outcomes = c("Memory", "Executive"),
                              n_classes = 2:4, robust_method = "t", n_starts = 2)
classical$fit_table[, c("Model", "LogLik", "BIC", "Entropy", "Min_Size")]
#>      Model    LogLik      BIC   Entropy Min_Size
#> 1 G2_slope -10816.56 21758.93 0.7899989   0.2150
#> 2 G3_slope -10771.35 21698.48 0.8195362   0.1900
#> 3 G4_slope -10743.39 21672.52 0.8552182   0.0025
robust$fit_table[, c("Model", "LogLik", "BIC", "Entropy", "Min_Size")]
#>      Model    LogLik      BIC   Entropy Min_Size
#> 1 G2_slope -10586.20 21304.21 0.7703904   0.2175
#> 2 G3_slope -10536.74 21235.26 0.8009036   0.1950
#> 3 G4_slope -10531.26 21254.25 0.7414121   0.0750

With the t model the log-likelihood is a proper likelihood, so the bootstrapped likelihood ratio test can complement BIC. It simulates data from the fitted G - 1-class model on the observed visit schedule and missingness pattern (slow: use 200 or more samples and several cores):

blrt_gmm_robust(neuro_long, id = "ID", time = "Year",
                outcomes = c("Memory", "Executive"), G = 3,
                robust_method = "t", n_samples = 200, cores = 4)

5. LASSO for trajectories

With several outcomes, two questions are natural: on which outcomes does each class actually change? and which outcomes differentiate the classes at all? Two penalties answer them:

lambda_growth shrinks the growth terms (slope, …) of every class toward zero: a slope set exactly to zero identifies a class that is stable on that outcome;
lambda_diff shrinks the trajectories of the classes toward each other; with group_diff = TRUE all the coefficients of an outcome are penalized together, so an outcome that does not differ between classes is removed from the class separation.

The penalties are adaptive by default: each term is weighted by its unpenalized estimate, so that lambda = z^2 / N sets to zero, roughly, the terms whose Wald statistic is below z. Here z = 3:

N <- length(unique(neuro_long$ID))
fit_l <- robust_gmm(neuro_long, id = "ID", time = "Year",
                    outcomes = c("Memory", "Executive", "Speed"), G = 3,
                    robust_method = "t", n_starts = 2,
                    lambda_growth = 9 / N, lambda_diff = 9 / N, group_diff = TRUE,
                    relax = TRUE)
summary(fit_l)
#> robust_gmm summary -- EM | G = 3 | persons = 400
#> Estimation: robust: multivariate t (nu = 10.30) | degree 1 | random: slope
#> 
#> Class trajectories (original scale):
#>  Class   Outcome (Intercept)   Year
#>      1    Memory      45.136 -5.004
#>      1 Executive      43.115 -3.937
#>      1     Speed      50.105 -0.606
#>      2    Memory      48.822 -2.132
#>      2 Executive      47.988 -1.643
#>      2     Speed      50.105 -0.606
#>      3    Memory      53.208  0.000
#>      3 Executive      52.145  0.000
#>      3     Speed      50.105 -0.606
#> 
#> Class sizes (modal assignment):
#>  Class   N Proportion
#>      1  77      0.183
#>      2 123      0.318
#>      3 200      0.499
#> 
#> Random-effect covariance (equal across classes; class 1 shown):
#>                       Memory:(Intercept) Memory:Year Executive:(Intercept)
#> Memory:(Intercept)                19.256       0.394                10.520
#> Memory:Year                        0.394       0.047                 0.104
#> Executive:(Intercept)             10.520       0.104                23.248
#> Executive:Year                    -0.116       0.073                 0.380
#> Speed:(Intercept)                 10.519      -0.231                11.665
#> Speed:Year                         0.073       0.031                 0.144
#>                       Executive:Year Speed:(Intercept) Speed:Year
#> Memory:(Intercept)            -0.116            10.519      0.073
#> Memory:Year                    0.073            -0.231      0.031
#> Executive:(Intercept)          0.380            11.665      0.144
#> Executive:Year                 0.229            -0.438      0.070
#> Speed:(Intercept)             -0.438            24.004      0.250
#> Speed:Year                     0.070             0.250      0.076
#> 
#> Residual variances:
#>         Memory Executive Speed
#> Class_1  6.875     5.695 6.675
#> Class_2  6.875     5.695 6.675
#> Class_3  6.875     5.695 6.675
#> 
#> Fit: LogLik = -15699.01 | Parameters = 39 | AIC = 31476.0 | BIC = 31631.7 | SABIC = 31507.9 | Entropy = 0.815
#> 
#> LASSO: growth = 0.0225, difference = 0.0225 (group by outcome), adaptive weights, estimates from the relaxed (unpenalized) refit
#> Outcomes differentiating the classes: Memory, Executive (removed: Speed)
#> Growth terms set to zero: 2 of 9

Speed is recognized as an outcome that does not differentiate the classes (the three classes share its trajectory), and the stable class has exactly zero slopes on Memory and Executive. With relax = TRUE the selected model is refitted without penalty (relaxed Lasso), so the reported estimates are not shrunk; its BIC counts only the free coefficients, and is lower than the BIC of the unpenalized three-outcome model:

fit_u <- robust_gmm(neuro_long, id = "ID", time = "Year",
                    outcomes = c("Memory", "Executive", "Speed"), G = 3,
                    robust_method = "t", n_starts = 2)
rbind(unpenalized = fit_u$fit[, c("LogLik", "Parameters", "BIC")],
      lasso_relaxed = fit_l$fit[, c("LogLik", "Parameters", "BIC")])
#>                  LogLik Parameters      BIC
#> unpenalized   -15697.92         45 31665.45
#> lasso_relaxed -15699.01         39 31631.69

Instead of fixing z, estimate_gmm_robust(tune_penalty = "bic") or tune_penalty = "cv" (cross-validation over persons) chooses it from a grid:

estimate_gmm_robust(neuro_long, id = "ID", time = "Year",
                    outcomes = c("Memory", "Executive", "Speed"), n_classes = 3,
                    tune_penalty = "bic", z_grid = c(1.5, 2, 2.5, 3, 4),
                    robust_method = "t", group_diff = TRUE, relax = TRUE)

6. What distinguished the classes at baseline?

bch_robust() relates the trajectory classes to a variable that was not used to estimate them – a baseline characteristic or a distal outcome – correcting for classification error (Bolck, Croon & Hagenaars, 2004). The auxiliary variable must have one value per person, in the order of fit$ids:

baseline <- neuro_long[!duplicated(neuro_long$ID), ]
baseline <- baseline[match(fit$ids, baseline$ID), ]
bch_biomarker <- bch_robust(fit, baseline$Biomarker)
round(bch_biomarker$Profile_Means)
#> Profile_1 Profile_2 Profile_3 
#>       696       895       851
bch_biomarker$ANOVA_Table
#>            Df  Sum_Sq    Mean_Sq  F_value      p_value
#> Class       2 2160216 1080108.18 47.83063 2.449367e-19
#> Residuals 397 8965027   22581.93       NA           NA

For publication, correction = "bootstrap" adds standard errors that account for the uncertainty of the classification (whole persons are resampled and the growth mixture model is refitted every time).

7. Bayesian estimation

The MCMC engine fits the same models by Gibbs sampling (with the LASSO penalties turned into Bayesian-Lasso priors). It starts from the EM solution and reports the WAIC and the Gelman-Rubin diagnostics:

fit_b <- robust_gmm(neuro_long, id = "ID", time = "Year", outcomes = "Memory",
                    G = 3, robust_method = "t", engine = "MCMC",
                    mcmc_iter = 600, n_chains = 2, n_starts = 2)
fit_b
#> <robust_gmm> MCMC | G = 3 | persons = 400 | outcomes: Memory
#> Trajectory: degree 1 | random: slope (full, equal across classes) | robust: multivariate t (nu = 6.58)
#> LogLik = -5450.8 | BIC = 10979.4 | Entropy = 0.738 | WAIC = 10926.0
#> Proportions: C1=0.19, C2=0.31, C3=0.51

plot_mcmc_chains(fit_b, pars = c("beta[1,Memory:Year]", "beta[2,Memory:Year]",
                                 "beta[3,Memory:Year]", "nu"))

MCMC trace plots for the Memory slopes of the three classes

8. Practical recommendations

Express time in meaningful units and center it where the intercept should be interpreted (e.g. years since baseline).
Start with shared variance components (re_cov = "equal", resid_var = "equal") and a random intercept and slope; free them only if the data support it (compare BIC).
Use several starts (n_starts) and check the smallest class size: tiny classes are often spurious.
Prefer robust_method = "t": it keeps a proper likelihood, so BIC, the BLRT and BCH are used as intended.
Use the LASSO to select, and report the relaxed estimates (relax = TRUE).

References

Bolck, A., Croon, M., & Hagenaars, J. (2004). Estimating latent structure models with categorical variables: One-step versus three-step estimators. Political Analysis, 12(1), 3-27.

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

Nagin, D. S. (1999). Analyzing developmental trajectories: A semiparametric, group-based approach. Psychological Methods, 4(2), 139-157.

Pinheiro, J. C., Liu, C., & Wu, Y. N. (2001). Efficient algorithms for robust estimation in linear mixed-effects models using the multivariate t distribution. Journal of Computational and Graphical Statistics, 10(2), 249-276.

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

Xie, B., Pan, W., & Shen, X. (2008). Variable selection in penalized model-based clustering via regularization on grouped parameters. Biometrics, 64(3), 921-930.

Zou, H. (2006). The adaptive lasso and its oracle properties. Journal of the American Statistical Association, 101(476), 1418-1429.