Getting Started with RobustLPA

1. Introduction

Latent Profile Analysis (LPA) groups observations into a small number of unobserved (“latent”) profiles based on a set of continuous indicator variables, by fitting a finite mixture of multivariate normal distributions. It is widely used in psychology, education, and the health sciences to identify subgroups of people who share a similar pattern of scores – without specifying the groups in advance.

Standard (maximum-likelihood) LPA estimation is not robust: a handful of extreme or mismeasured observations can distort the estimated profile means and covariances, sometimes badly enough to change which observations end up in which profile (Garcia-Escudero et al., 2010). RobustLPA provides:

Two estimation engines: classical Expectation-Maximization (EM) and Bayesian MCMC (Gibbs sampling), both with an optional Huber-type robust mode (robust = TRUE, the default) that down-weights outlying observations, or a classical (non-robust) mode (robust = FALSE) for comparison.
Six variance-covariance parameterizations (model = 1:6), from a single shared diagonal covariance to a fully unconstrained covariance per profile, so model complexity can be chosen to fit the data rather than assumed.
LASSO-type regularization of the profile means, with cross-validated penalty selection.
Native handling of missing data via Full Information Maximum Likelihood (FIML) – no need to drop or impute incomplete rows.
Tools for the two questions every LPA analysis has to answer: how many profiles? (estimate_profiles_robust(), the bootstrapped likelihood ratio test blrt_robust()) and how do profiles relate to variables outside the model? (bch_robust(), implementing the Bolck-Croon-Hagenaars three-step method).
Parallel computing (cores =) throughout, for the EM restarts / MCMC chains of a single fit, for grid searches over models and profile counts, and for the bootstrap procedures.

This vignette walks through a complete analysis on the dataset bundled with the package, neuro_data. Progress messages and the log-likelihood-decrease warnings that Huber-weighted robust estimation can occasionally emit (expected behavior, explained in the “Robust estimation” section of ?robust_lpa) are suppressed below to keep the output readable; they do not affect any of the fitted values shown.

2. The example dataset

neuro_data contains simulated neuropsychological test scores and reaction times for 250 people belonging to two true, known groups: “Healthy” (n = 150) and “Pathological” (n = 100). The group label (True_Profile) is included only so that recovered profiles can be checked against ground truth – it is never used for estimation, since LPA is unsupervised.

data(neuro_data)
str(neuro_data)
#> 'data.frame':    250 obs. of  7 variables:
#>  $ ID                 : int  1 2 3 4 5 6 7 8 9 10 ...
#>  $ True_Profile       : chr  "Healthy" "Healthy" "Healthy" "Healthy" ...
#>  $ Memory             : num  81.7 81.6 96.1 80.9 72.3 ...
#>  $ Attention          : num  66.2 76.7 80.3 71.7 77.6 ...
#>  $ Executive_Functions: num  69.3 75.4 78.6 76.9 64.6 ...
#>  $ RT_Stroop          : num  399 446 458 434 361 ...
#>  $ RT_TMT             : num  441 430 375 467 455 ...
table(neuro_data$True_Profile)
#> 
#>      Healthy Pathological 
#>          150          100

Two of the five continuous variables (Attention, Executive_Functions) are, by design, identically distributed in both groups: they carry no group signal and act as “noise” variables. Memory, RT_Stroop, and RT_TMT differ between groups, and the two reaction-time variables additionally differ in variance and in how strongly they correlate with each other – a genuine difference in covariance structure, not just location, between the two groups (see ?neuro_data). A subset of the Pathological group also carries extra, variable-magnitude outlying values on the two reaction-time variables, simulating measurement contamination.

As with any LPA analysis, we standardize the indicators first, since several parts of the package (LASSO shrinkage in particular) are only meaningful on a common scale:

vars <- c("Memory", "Attention", "Executive_Functions", "RT_Stroop", "RT_TMT")
x <- scale(as.matrix(neuro_data[, vars]))
head(x)
#>           Memory  Attention Executive_Functions  RT_Stroop     RT_TMT
#> [1,]  0.41669327 -0.1948912          -0.3123633 -0.6660181 -0.6352654
#> [2,]  0.41163048  1.0144344           0.4423769 -0.2481754 -0.7173859
#> [3,]  1.15527236  1.4246020           0.8401542 -0.1408952 -1.1412761
#> [4,]  0.37585978  0.4405564           0.6296725 -0.3562865 -0.4329801
#> [5,] -0.06607811  1.1194497          -0.8959079 -0.9913752 -0.5272060
#> [6,]  0.03892365 -1.1028652          -1.4977872 -0.5733588 -0.6258035

3. Choosing a variance-covariance model

robust_lpa()‘s model argument selects how the profiles’ covariance matrices are constrained, from most to least parsimonious:

`model`	Variances across profiles	Covariances across profiles
1	Equal (shared)	Zero (diagonal), shared
2	Varying	Zero (diagonal), own
3	Equal (shared)	Equal (shared), full
4	Varying	Shared correlation structure, own variances
5	Equal (shared)	Own correlation structure, shared variances
6	Varying	Varying (fully unconstrained per profile)

More parsimonious models (1-2) are more stable with smaller samples but can under-fit real covariance structure; less parsimonious models (especially 6) can fit better but need more data and are more prone to numerically unstable, near-singular covariance estimates for small or overlapping profiles – robust_lpa() guards against this automatically and will warn if a fitted profile ends up implausibly small (see ?robust_lpa). Section 7 below shows how to let BIC choose among all six objectively, rather than assuming one.

4. Fitting a single model with the EM engine

fit_em <- robust_lpa(x, G = 2, model = 6, n_starts = 5)
fit_em
#> <robust_lpa> EM | model 6 | G = 2 | N = 250
#> LogLik = -1267.0 | AIC = 2616.0 | BIC = 2760.4 | Entropy = 0.864
#> Proportions: P1=0.61, P2=0.39

summary() adds per-profile means and sizes:

summary(fit_em)
#> robust_lpa summary -- EM | model 6 | G = 2 | N = 250
#> 
#> Profile means:
#>                        P1    P2
#> Memory               0.57 -0.89
#> Attention           -0.02  0.04
#> Executive_Functions  0.06 -0.10
#> RT_Stroop           -0.63  0.99
#> RT_TMT              -0.62  0.99
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1 157      0.612
#>       P2  93      0.388
#> 
#> Fit: LogLik = -1267.0 | AIC = 2616.0 | BIC = 2760.4 | Entropy = 0.864

Since neuro_data includes the ground-truth group label, we can check how well the fitted profiles recover it:

table(True_Profile = neuro_data$True_Profile, Assigned = fit_em$assignments)
#>               Assigned
#> True_Profile     1   2
#>   Healthy      145   5
#>   Pathological  12  88

Robust vs. classical estimation

robust = TRUE (the default) down-weights each observation’s contribution to its assigned profile’s mean/covariance once its Mahalanobis distance exceeds a chi-squared cutoff (controlled by alpha). Comparing against robust = FALSE shows the effect of the contamination built into RT_Stroop/RT_TMT:

fit_classical <- robust_lpa(x, G = 2, model = 6, n_starts = 5, robust = FALSE)
rbind(
  robust    = sapply(fit_em$means, `[`, "RT_Stroop"),
  classical = sapply(fit_classical$means, `[`, "RT_Stroop")
)
#>            RT_Stroop  RT_Stroop
#> robust    -0.6295026  0.9918472
#> classical  0.9910631 -0.6269938

5. LASSO regularization

For higher-dimensional indicator sets, lambda applies LASSO-type soft-thresholding shrinkage to the profile means (meaningful only on standardized data, as used throughout this vignette):

fit_lasso <- robust_lpa(x, G = 2, model = 6, n_starts = 3, lambda = 0.15)
summary(fit_lasso)
#> robust_lpa summary -- EM | model 6 | G = 2 | N = 250
#> 
#> Profile means:
#>                        P1    P2
#> Memory               0.34 -0.35
#> Attention            0.00  0.00
#> Executive_Functions  0.00  0.00
#> RT_Stroop           -0.45  0.47
#> RT_TMT              -0.44  0.46
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1 141      0.509
#>       P2 109      0.491
#> 
#> Fit: LogLik = -1305.3 | AIC = 2684.6 | BIC = 2814.9 | Entropy = 0.633

Rather than fixing lambda by hand, estimate_profiles_robust(tune_lasso = TRUE) selects it by k-fold cross-validation (see Section 7).

6. Missing data (FIML)

robust_lpa() handles missing values natively via Full Information Maximum Likelihood – no listwise deletion or imputation needed – for both engines and every variance-covariance model:

x_na <- x
set.seed(1)
na_idx <- cbind(
  sample(nrow(x_na), 15),
  sample(ncol(x_na), 15, replace = TRUE)
)
x_na[na_idx] <- NA
mean(is.na(x_na))
#> [1] 0.012

fit_fiml <- robust_lpa(x_na, G = 2, model = 6, n_starts = 5)
summary(fit_fiml)
#> robust_lpa summary -- EM | model 6 | G = 2 | N = 250
#> 
#> Profile means:
#>                        P1    P2
#> Memory               0.57 -0.89
#> Attention           -0.02  0.02
#> Executive_Functions  0.06 -0.09
#> RT_Stroop           -0.63  0.99
#> RT_TMT              -0.62  0.97
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1 156       0.61
#>       P2  94       0.39
#> 
#> Fit: LogLik = -1256.6 | AIC = 2595.1 | BIC = 2739.5 | Entropy = 0.858

7. Choosing the number of profiles and the covariance model

estimate_profiles_robust() fits every combination of n_profiles and models and collects their fit indices in one table, so models can be compared by AIC/BIC/SABIC rather than assumed in advance:

grid <- estimate_profiles_robust(x, n_profiles = 1:3, models = 1:6, n_starts = 5)
grid$fit_table[order(grid$fit_table$BIC), ]
#>    Model Profiles    LogLik Parameters      AIC      BIC    SABIC    Entropy
#> 17     6        2 -1267.024         41 2616.048 2760.428 2630.455 0.86419639
#> 18     6        3 -1235.911         62 2595.821 2814.152 2617.607 0.91629446
#> 8      3        2 -1339.403         26 2730.806 2822.364 2739.942 0.92443096
#> 9      3        3 -1326.983         32 2717.967 2830.653 2729.211 0.78779567
#> 1      1        1 -1396.510         10 2813.020 2848.235 2816.534 1.00000000
#> 4      2        1 -1396.510         10 2813.020 2848.235 2816.534 1.00000000
#> 15     5        3 -1290.310         52 2684.621 2867.737 2702.892 0.96146470
#> 2      1        2 -1395.965         16 2823.929 2880.273 2829.551 0.01849328
#> 11     4        2 -1360.579         31 2783.158 2892.323 2794.051 0.57536898
#> 12     4        3 -1335.707         42 2755.414 2903.316 2770.172 0.83301661
#> 7      3        1 -1396.510         20 2833.020 2903.450 2840.048 1.00000000
#> 10     4        1 -1396.510         20 2833.020 2903.450 2840.048 1.00000000
#> 13     5        1 -1396.510         20 2833.020 2903.450 2840.048 1.00000000
#> 16     6        1 -1396.510         20 2833.020 2903.450 2840.048 1.00000000
#> 3      1        3 -1392.067         22 2828.133 2905.606 2835.864 0.95448733
#> 5      2        2 -1395.596         21 2833.192 2907.143 2840.571 0.02707980
#> 6      2        3 -1374.018         32 2812.036 2924.723 2823.280 0.86454683
#> 14     5        2 -1383.117         36 2838.233 2965.006 2850.883 0.32571177
#>    Min_Size Max_Size Lambda
#> 17    0.372    0.628      0
#> 18    0.064    0.612      0
#> 8     0.284    0.716      0
#> 9     0.176    0.540      0
#> 1     1.000    1.000      0
#> 4     1.000    1.000      0
#> 15    0.004    0.708      0
#> 2     0.480    0.520      0
#> 11    0.452    0.548      0
#> 12    0.104    0.728      0
#> 7     1.000    1.000      0
#> 10    1.000    1.000      0
#> 13    1.000    1.000      0
#> 16    1.000    1.000      0
#> 3     0.120    0.672      0
#> 5     0.448    0.552      0
#> 6     0.256    0.424      0
#> 14    0.448    0.552      0

neuro_data’s genuine group-level covariance difference (Section 2) was specifically calibrated so that the fully unconstrained model (model = 6) at two profiles fits measurably better than more parsimonious alternatives, despite its larger parameter penalty – if you reproduce this table, Model = 6, Profiles = 2 should be at or very near the top by BIC. Each element of grid$models is a fitted robust_lpa object:

summary(grid$models[["model_6_profiles_2"]])
#> robust_lpa summary -- EM | model 6 | G = 2 | N = 250
#> 
#> Profile means:
#>                        P1    P2
#> Memory              -0.89  0.57
#> Attention            0.04 -0.02
#> Executive_Functions -0.10  0.06
#> RT_Stroop            0.99 -0.63
#> RT_TMT               0.99 -0.62
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1  93      0.388
#>       P2 157      0.612
#> 
#> Fit: LogLik = -1267.0 | AIC = 2616.0 | BIC = 2760.4 | Entropy = 0.864

plot_robust_lpa() accepts either a single fit or a full grid (in which case it plots the lowest-BIC model automatically):

plot_robust_lpa(grid, title = "Best-fitting model (lowest BIC)")

Profile plot of the best-fitting model

Cross-validated LASSO tuning uses the same grid interface:

grid_lasso <- estimate_profiles_robust(
  x, n_profiles = 2, models = 6, n_starts = 3,
  tune_lasso = TRUE, k_folds = 5, lambda_grid = c(0, 0.05, 0.1, 0.2)
)
grid_lasso$fit_table[, c("Model", "Profiles", "BIC", "Lambda")]
#>   Model Profiles      BIC Lambda
#> 1     6        2 2760.428      0

8. Confirming the number of profiles: the bootstrapped likelihood ratio test

BIC alone does not come with a significance test for “is G profiles actually better than G - 1?”. blrt_robust() answers this via parametric bootstrap (Nylund, Asparouhov & Muthen, 2007): it simulates data under the simpler (G - 1)-profile model, refits both models to each simulated dataset, and builds a reference distribution for the observed likelihood ratio. n_samples is kept small below for a fast vignette build; for publication-grade inference use at least 200-500 (and consider cores > 1, see Section 11):

blrt_res <- blrt_robust(x, G = 2, model = 6, n_samples = 20, n_starts = 3)
blrt_res
#> $LRT_Observed
#> [1] 258.9722
#> 
#> $Bootstrap_LRTs
#>  [1] 49.27097 41.25416 35.18181 40.00441 40.98219 44.83117 35.82499 37.53378
#>  [9] 45.01754 37.92280 49.20367 58.41630 37.32596 30.60113 44.64631 38.92990
#> [17] 27.91325 38.86899 33.34014 26.18940
#> 
#> $p_value
#> [1] 0.04761905
#> 
#> $Bootstrap_Failures
#> [1] 0

A small p_value supports keeping the second profile over collapsing to a single one.

9. Bayesian MCMC estimation

The MCMC engine estimates the same variance-covariance models via Gibbs sampling, optionally under a Bayesian Lasso (Laplace) prior on the profile means (prior_laplace), and runs multiple chains by default so convergence can be checked. mcmc_iter is kept small below for a fast vignette build; production analyses should use several thousand iterations:

fit_mcmc <- robust_lpa(x, G = 2, model = 6, engine = "MCMC",
                        mcmc_iter = 500, n_chains = 4, prior_laplace = 0.1)
summary(fit_mcmc)
#> robust_lpa summary -- MCMC | model 6 | G = 2 | N = 250
#> 
#> Profile means:
#>                        P1    P2
#> Memory               0.56 -0.86
#> Attention           -0.03  0.04
#> Executive_Functions  0.06 -0.09
#> RT_Stroop           -0.63  0.96
#> RT_TMT              -0.62  0.95
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1 157      0.605
#>       P2  93      0.395
#> 
#> Fit: LogLik = -1268.0 | AIC = 2618.0 | BIC = 2762.4 | Entropy = 0.845
#> MCMC: 4 chains x 500 iter | Rhat [1.00, 1.02] | ESS [271, 1061]

The summary’s Rhat/ESS range comes from the classic Gelman-Rubin potential scale reduction statistic and effective sample size (fit_mcmc$mcmc_diagnostics has the full per-parameter table); values of Rhat near 1 support convergence. plot_mcmc_chains() draws overlaid per-chain trace plots for visual inspection – pass pars to select a subset of the "mu[...]"/"sigma[...]"/"pi[...]" parameters (see ?plot_mcmc_chains) when there are many:

plot_mcmc_chains(fit_mcmc, pars = c("mu[1,1]", "mu[2,1]", "pi[1]"))

MCMC trace plots for two profile means and one mixing proportion

10. Relating profiles to an outside variable: the BCH method

A common follow-up question is whether the fitted profiles differ on a variable that was not used to estimate them (a distal outcome), while correctly accounting for classification error in the profile assignments (naively comparing group means on the hard-assigned profiles understates this error and biases the comparison). bch_robust() implements the three-step Bolck-Croon-Hagenaars (2004) method for this.

To keep this a genuine “outside variable” rather than one already in the measurement model, this section fits a reduced model that leaves RT_TMT out, so it can legitimately serve as the auxiliary/distal outcome:

x_reduced <- scale(as.matrix(neuro_data[, c("Memory", "Attention",
                                             "Executive_Functions", "RT_Stroop")]))
fit_reduced <- robust_lpa(x_reduced, G = 2, model = 6, n_starts = 5)

bch_res <- bch_robust(fit_reduced, neuro_data$RT_TMT)
bch_res$Profile_Means
#> Profile_1 Profile_2 
#>  440.2881  683.6163
bch_res$ANOVA_Table
#>            Df    Sum_Sq     Mean_Sq  F_value      p_value
#> Class       1 3318369.2 3318369.164 950.9043 7.886382e-87
#> Residuals 248  865445.2    3489.698       NA           NA

$ANOVA_Table’s F-test treats the classification error matrix as fixed, which can understate uncertainty (Vermunt, 2010). correction = "bootstrap" adds a nonparametric approximation to the Bakk, Oberski & Vermunt (2014) sandwich correction – bootstrap standard errors, confidence intervals, and a Wald test – at the cost of refitting the step-1 model n_boot times:

bch_boot <- bch_robust(fit_reduced, neuro_data$RT_TMT,
                        correction = "bootstrap", n_boot = 30)
bch_boot$Bootstrap_Correction
#> $n_boot_used
#> [1] 30
#> 
#> $n_boot_failed
#> [1] 0
#> 
#> $SE
#> Profile_1 Profile_2 
#>  4.632717 14.944461 
#> 
#> $CI_lower
#> Profile_1 Profile_2 
#>  432.0065  655.2415 
#> 
#> $CI_upper
#> Profile_1 Profile_2 
#>  450.1307  706.6348 
#> 
#> $Wald_stat
#> [1] 300.201
#> 
#> $Wald_df
#> [1] 1
#> 
#> $Wald_p_value
#> [1] 2.978374e-67

(As with the BLRT, n_boot is kept small here for a fast vignette build; use several hundred for publication-grade inference.)

11. Parallel computing

Every bootstrap- or restart-based procedure in this package accepts a cores argument: EM random restarts or MCMC chains within a single robust_lpa() call, the model/profile grid (and cross-validation folds) in estimate_profiles_robust(), bootstrap replicates in blrt_robust(), and bootstrap correction replicates in bch_robust(). These are not run in this vignette (CRAN’s check machines cap how many cores a package may use during checks), but the calls are otherwise identical to the sequential versions above:

grid_parallel <- estimate_profiles_robust(x, n_profiles = 1:3, models = 1:6,
                                           n_starts = 5, cores = 4)
fit_mcmc_parallel <- robust_lpa(x, G = 2, model = 6, engine = "MCMC",
                                 mcmc_iter = 2000, n_chains = 4, cores = 4)

If you also parallelize an outer loop (e.g. blrt_robust(cores = )) around calls that themselves use cores, keep the product of the two values at or below your machine’s core count to avoid oversubscription.

12. Summary

Task	Function
Fit one model	`robust_lpa()`
Compare models/profile counts	`estimate_profiles_robust()`, `plot_robust_lpa()`
Test the number of profiles	`blrt_robust()`
Relate profiles to an outside variable	`bch_robust()`
Inspect MCMC convergence	`plot_mcmc_chains()`, `fit$mcmc_diagnostics`
Quick robust centroid (no mixture model)	`robust_mean()`

See the function help pages (?robust_lpa, ?estimate_profiles_robust, ?blrt_robust, ?bch_robust, ?plot_mcmc_chains, ?neuro_data) for full argument documentation, and NEWS.md for what changed in this release.

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.

Bakk, Z., Oberski, D. L., & Vermunt, J. K. (2014). Relating latent class assignments to external variables: Standard errors for correct inference. Political Analysis, 22(4), 520-540.

Garcia-Escudero, L. A., Gordaliza, A., Matran, C., & Mayo-Iscar, A. (2010). A review of robust clustering methods. Advances in Data Analysis and Classification, 4(2-3), 89-109.

Nylund, K. L., Asparouhov, T., & Muthen, B. O. (2007). Deciding on the number of classes in latent class analysis and growth mixture modeling: A Monte Carlo simulation study. Structural Equation Modeling, 14(4), 535-569.

Vermunt, J. K. (2010). Latent class modeling with covariates: Two improved three-step approaches. Political Analysis, 18(4), 450-469.