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:

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

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

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]"))

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

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.