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: Expectation-Maximization (EM) and Bayesian MCMC (Gibbs sampling).
• Two robust estimators: Huber-type down-weighting of outlying observations (robust_method = "huber", the default) and mixtures of multivariate t distributions (robust_method = "t"), a likelihood-based robust model; plus classical Gaussian estimation (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 by full-information maximum likelihood (the exact EM for incomplete data; data augmentation in the MCMC engine) – 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 are suppressed below to keep the output readable.

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 | robust: Huber
#> LogLik = -1267.0 | AIC = 2616.1 | BIC = 2760.4 | Entropy = 0.864
#> Proportions: P1=0.39, P2=0.61

summary() adds per-profile means and sizes:

summary(fit_em)
#> robust_lpa summary -- EM | model 6 | G = 2 | N = 250
#> Estimation: robust: Huber
#> 
#> 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.1 | 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        5 145
#>   Pathological  88  12

fit_t <- robust_lpa(x, G = 2, model = 6, n_starts = 5, robust_method = "t")
fit_classical <- robust_lpa(x, G = 2, model = 6, n_starts = 5, robust = FALSE)
rbind(
  huber     = sort(sapply(fit_em$means, `[`, "RT_Stroop")),
  t         = sort(sapply(fit_t$means, `[`, "RT_Stroop")),
  classical = sort(sapply(fit_classical$means, `[`, "RT_Stroop"))
)
#>            RT_Stroop RT_Stroop
#> huber     -0.6296193 0.9918706
#> t         -0.6264657 0.9927519
#> classical -0.6269693 0.9913560
fit_t$nu
#> [1] 199.9213

set.seed(6)
contaminated <- matrix(rnorm(400 * 3), 400, 3)
idx <- sample(400, 20)
contaminated[idx, ] <- contaminated[idx, ] + matrix(rnorm(60, 0, 15), 20, 3)
sapply(list(
  classical = robust_lpa(contaminated, G = 1, model = 6, robust = FALSE),
  huber     = robust_lpa(contaminated, G = 1, model = 6),
  t         = robust_lpa(contaminated, G = 1, model = 6, robust_method = "t")
), function(f) round(diag(f$covariances[[1]]), 2))
#>      classical huber    t
#> [1,]      7.54  2.35 0.79
#> [2,]     12.45  3.06 0.74
#> [3,]     15.95  3.55 0.83

head(order(fit_t$weights))
#> [1] 133 134  44 159  10 178
round(head(sort(fit_t$weights)), 3)
#> [1] 0.950 0.954 0.957 0.962 0.967 0.968

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
#> Estimation: robust: Huber
#> 
#> Profile means:
#>                        P1    P2
#> Memory               0.32 -0.31
#> Attention            0.00  0.00
#> Executive_Functions  0.00  0.00
#> RT_Stroop           -0.44  0.43
#> RT_TMT              -0.44  0.43
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1 138      0.495
#>       P2 112      0.505
#> 
#> Fit: LogLik = -1307.7 | AIC = 2689.5 | BIC = 2819.8 | Entropy = 0.629

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 by full-information maximum likelihood – no listwise deletion or imputation needed – for both engines and every variance-covariance model. Each incomplete row contributes the likelihood of its observed entries, and the M-step uses the exact EM for incomplete data (missing entries are replaced by their conditional expectations and the corresponding conditional covariance is added), so the estimates are maximum likelihood when data are missing at random:

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
#> Estimation: robust: Huber
#> 
#> Profile means:
#>                        P1    P2
#> Memory              -0.89  0.57
#> Attention            0.03 -0.02
#> Executive_Functions -0.09  0.06
#> RT_Stroop            0.99 -0.63
#> RT_TMT               0.97 -0.63
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1  93       0.39
#>       P2 157       0.61
#> 
#> Fit: LogLik = -1256.3 | AIC = 2594.5 | BIC = 2738.9 | Entropy = 0.859

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.028         41 2616.056 2760.436 2630.462 0.8642464
#> 11     4        2 -1315.465         31 2692.931 2802.096 2703.824 0.9393751
#> 9      3        3 -1316.296         32 2696.592 2809.278 2707.836 0.9300158
#> 12     4        3 -1294.760         42 2673.519 2821.421 2688.277 0.8622138
#> 8      3        2 -1339.403         26 2730.807 2822.365 2739.943 0.9252260
#> 18     6        3 -1245.940         62 2615.879 2834.210 2637.665 0.8617985
#> 14     5        2 -1318.894         36 2709.788 2836.561 2722.438 0.9209662
#> 15     5        3 -1283.323         52 2670.646 2853.762 2688.917 0.9544708
#> 13     5        1 -1396.513         20 2833.026 2903.455 2840.053 1.0000000
#> 7      3        1 -1396.513         20 2833.026 2903.455 2840.053 1.0000000
#> 10     4        1 -1396.513         20 2833.026 2903.455 2840.053 1.0000000
#> 16     6        1 -1396.513         20 2833.026 2903.455 2840.053 1.0000000
#> 3      1        3 -1392.130         22 2828.260 2905.732 2835.990 0.9546128
#> 6      2        3 -1372.870         32 2809.739 2922.426 2820.983 0.8948553
#> 5      2        2 -1428.107         21 2898.214 2972.165 2905.593 0.9570553
#> 2      1        2 -1470.082         16 2972.164 3028.507 2977.786 0.9646582
#> 1      1        1 -1772.511         10 3565.022 3600.237 3568.536 1.0000000
#> 4      2        1 -1772.511         10 3565.022 3600.237 3568.536 1.0000000
#>    Min_Size Max_Size Lambda
#> 17    0.372    0.628      0
#> 11    0.340    0.660      0
#> 9     0.128    0.664      0
#> 12    0.180    0.480      0
#> 8     0.284    0.716      0
#> 18    0.196    0.444      0
#> 14    0.288    0.712      0
#> 15    0.128    0.664      0
#> 13    1.000    1.000      0
#> 7     1.000    1.000      0
#> 10    1.000    1.000      0
#> 16    1.000    1.000      0
#> 3     0.120    0.672      0
#> 6     0.200    0.464      0
#> 5     0.332    0.668      0
#> 2     0.284    0.716      0
#> 1     1.000    1.000      0
#> 4     1.000    1.000      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
#> Estimation: robust: Huber
#> 
#> 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.1 | 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.436      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.9701
#> 
#> $Bootstrap_LRTs
#>  [1] 21.56804 36.72249 34.05213 49.08376 34.09293 39.52475 31.00575 33.87822
#>  [9] 41.46293 17.65204 32.17148 38.85211 30.65625 23.88643 37.75806 35.02918
#> [17] 32.81161 36.66562 20.48867 34.74985
#> 
#> $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, under a Bayesian Lasso (Laplace) prior on the profile means (prior_laplace), and runs multiple chains by default so convergence can be checked. The chains start from a preliminary EM fit with dispersed perturbations, and draws are relabeled to that EM solution to resolve label switching. With robust_method = "t" the sampler is an exact Gibbs sampler for the multivariate-t mixture. 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", robust_method = "t",
                       mcmc_iter = 500, n_chains = 4, prior_laplace = 0.1)
summary(fit_mcmc)
#> robust_lpa summary -- MCMC | model 6 | G = 2 | N = 250
#> Estimation: robust: multivariate t (nu = 48.80)
#> 
#> Profile means:
#>                        P1    P2
#> Memory               0.56 -0.88
#> Attention           -0.02  0.05
#> Executive_Functions  0.07 -0.09
#> RT_Stroop           -0.62  0.97
#> RT_TMT              -0.62  0.97
#> 
#> Profile sizes:
#>  Profile   N Proportion
#>       P1 158      0.611
#>       P2  92      0.389
#> 
#> Fit: LogLik = -1270.8 | AIC = 2625.5 | BIC = 2773.4 | Entropy = 0.854 | WAIC = 2620.2
#> MCMC: 4 chains x 500 iter | Rhat [1.00, 1.02] | ESS [199, 963]

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. WAIC (widely applicable information criterion; lower is better) is the recommended criterion for comparing MCMC fits. 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]", "nu"))

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 
#>  430.2883  664.8225
bch_res$ANOVA_Table
#>            Df  Sum_Sq     Mean_Sq  F_value      p_value
#> Class       1 3285353 3285353.360 906.8481 8.226685e-85
#> Residuals 248  898461    3622.827       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 
#>  7.278576 19.395291 
#> 
#> $CI_lower
#> Profile_1 Profile_2 
#>  414.0777  610.2212 
#> 
#> $CI_upper
#> Profile_1 Profile_2 
#>  443.5781  683.0440 
#> 
#> $Wald_stat
#> [1] 216.0096
#> 
#> $Wald_df
#> [1] 1
#> 
#> $Wald_p_value
#> [1] 6.711578e-49

(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 in estimate_profiles_robust(), bootstrap replicates in blrt_robust(), and bootstrap correction replicates in bch_robust(). One random seed is drawn per unit of work before dispatch, so after set.seed() the results are identical whatever the number of cores (on the same machine; different operating systems or linear-algebra libraries can differ in the last digits). 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.

Peel, D., & McLachlan, G. J. (2000). Robust mixture modelling using the t distribution. Statistics and Computing, 10(4), 339-348.

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