--- title: "Getting Started with RobustLPA" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting Started with RobustLPA} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4.2, warning = FALSE, message = FALSE ) ``` ```{r setup} library(RobustLPA) set.seed(2026) # every result below is reproducible, for any number of cores ``` ## 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. ```{r} data(neuro_data) str(neuro_data) table(neuro_data$True_Profile) ``` 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: ```{r} vars <- c("Memory", "Attention", "Executive_Functions", "RT_Stroop", "RT_TMT") x <- scale(as.matrix(neuro_data[, vars])) head(x) ``` ## 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 ```{r} fit_em <- robust_lpa(x, G = 2, model = 6, n_starts = 5) fit_em ``` `summary()` adds per-profile means and sizes: ```{r} summary(fit_em) ``` Since `neuro_data` includes the ground-truth group label, we can check how well the fitted profiles recover it: ```{r} table(True_Profile = neuro_data$True_Profile, Assigned = fit_em$assignments) ``` ### Robust vs. classical estimation `robust = TRUE` (the default) down-weights outlying observations. With `robust_method = "huber"` (the default), each observation's contribution to a profile's mean/covariance is down-weighted once its Mahalanobis distance to that profile's current robust estimates exceeds a chi-squared cutoff (controlled by `alpha`). With `robust_method = "t"`, every profile is a multivariate t distribution whose degrees of freedom `nu` are estimated from the data: heavy tails absorb outliers automatically, and -- unlike Huber weighting -- the model has a proper likelihood, so AIC/BIC, the bootstrapped likelihood ratio test and the BCH method are used exactly as intended. Comparing both against `robust = FALSE` (profile labels are arbitrary in every fit, so the profile means are sorted before comparing): ```{r} 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")) ) fit_t$nu ``` On `neuro_data` the contamination is mild relative to the profile-specific covariance that `model = 6` allows, so the three estimators agree closely and the estimated `nu` sits at its upper bound (the t mixture is then essentially Gaussian): robustness costs little when it is not needed. The differences grow with the size and number of outliers. Below, 5% of the rows of otherwise standard-normal data (true variances 1) receive gross errors; the diagonal of the estimated covariance (for the t model: scale) matrix shows how much each estimator is pulled by them: ```{r} 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)) ``` Every fit also reports a robustness weight per observation (`$weights`, 1 = full weight); the smallest weights point to the most outlying cases: ```{r} head(order(fit_t$weights)) round(head(sort(fit_t$weights)), 3) ``` ## 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): ```{r} fit_lasso <- robust_lpa(x, G = 2, model = 6, n_starts = 3, lambda = 0.15) summary(fit_lasso) ``` 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: ```{r} 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)) fit_fiml <- robust_lpa(x_na, G = 2, model = 6, n_starts = 5) summary(fit_fiml) ``` ## 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: ```{r} grid <- estimate_profiles_robust(x, n_profiles = 1:3, models = 1:6, n_starts = 5) grid$fit_table[order(grid$fit_table$BIC), ] ``` `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: ```{r} summary(grid$models[["model_6_profiles_2"]]) ``` `plot_robust_lpa()` accepts either a single fit or a full grid (in which case it plots the lowest-BIC model automatically): ```{r, fig.alt = "Profile plot of the best-fitting model"} plot_robust_lpa(grid, title = "Best-fitting model (lowest BIC)") ``` Cross-validated LASSO tuning uses the same grid interface: ```{r} 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")] ``` ## 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): ```{r} blrt_res <- blrt_robust(x, G = 2, model = 6, n_samples = 20, n_starts = 3) blrt_res ``` 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: ```{r} 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) ``` 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: ```{r, fig.alt = "MCMC trace plots for two profile means, one mixing proportion and the t degrees of freedom"} 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: ```{r} 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 bch_res$ANOVA_Table ``` `$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: ```{r} bch_boot <- bch_robust(fit_reduced, neuro_data$RT_TMT, correction = "bootstrap", n_boot = 30) bch_boot$Bootstrap_Correction ``` (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: ```{r, eval = FALSE} 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()` | | Latent classes of longitudinal trajectories | `robust_gmm()`, `estimate_gmm_robust()`, `blrt_gmm_robust()`, `plot_robust_gmm()` (see `vignette("robust-growth-mixture")`) | 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.