Package {RobustLPA}

Title:	Robust Latent Profile Analysis
Version:	1.1.0
Description:	Provides a comprehensive toolset for estimating Latent Profile Analysis (LPA) models that are robust to multivariate outliers and missing data. By integrating a high-performance 'C++' engine via 'RcppArmadillo', it reliably extracts latent profiles using both Expectation-Maximization (EM) and Markov Chain Monte Carlo (MCMC) Bayesian estimation. Robustness is obtained either by Huber-type down-weighting or by mixtures of multivariate t distributions (a likelihood-based robust model, see Peel and McLachlan (2000) <doi:10.1023/A:1008981510081>). Missing data are handled by full information maximum likelihood with the exact EM treatment of incomplete observations (data augmentation in the MCMC engine). The EM engine also supports LASSO regularization with k-fold cross-validation for penalty tuning; the MCMC engine uses a Bayesian Lasso with Laplace priors, multiple chains, Gelman-Rubin/effective sample size diagnostics and the widely applicable information criterion. It supports six geometric variance-covariance models, along with functions for bootstrapped likelihood ratio tests (BLRT), BCH auxiliary variable analysis, and plotting. For longitudinal data, it fits robust growth mixture models and latent class growth analysis (Muthen and Shedden (1999) <doi:10.1111/j.0006-341X.1999.00463.x>) for one or several outcomes measured on unbalanced occasions, with Gaussian, Huber-weighted or multivariate-t (Pinheiro, Liu and Wu (2001) <doi:10.1198/10618600152628059>) latent classes, by EM and MCMC, and optional adaptive LASSO penalties (Zou (2006) <doi:10.1198/016214506000000735>) that identify stable trajectories and the outcomes that differentiate the classes. For methodological details on the Bootstrapped Likelihood Ratio Test, see Nylund et al. (2007) <doi:10.1080/10705510701575396>. For robust clustering methods, see Garcia-Escudero et al. (2010) <doi:10.1007/s11634-010-0064-5>. For BCH auxiliary variable analysis, see Bolck et al. (2004) <doi:10.1093/pan/mph001>.
License:	GPL (≥ 3)
Encoding:	UTF-8
LinkingTo:	Rcpp, RcppArmadillo
Imports:	Rcpp, ggplot2, stats, utils, bayesplot, coda
Suggests:	parallel, knitr, rmarkdown, testthat (≥ 3.1.5), lme4, nlme
Config/testthat/edition:	3
VignetteBuilder:	knitr
Depends:	R (≥ 3.6)
LazyData:	true
NeedsCompilation:	yes
Config/roxygen2/version:	8.0.0
Packaged:	2026-09-27 20:01:23 UTC; hp
Author:	Valerio Riccardo Aquila [aut, cre]
Maintainer:	Valerio Riccardo Aquila <valerio_aquila@hotmail.it>
Repository:	CRAN
Date/Publication:	2026-09-27 22:50:15 UTC

RobustLPA: Robust Latent Profile Analysis

Description

Provides a comprehensive toolset for estimating Latent Profile Analysis (LPA) models that are robust to multivariate outliers and missing data. By integrating a high-performance 'C++' engine via 'RcppArmadillo', it reliably extracts latent profiles using both Expectation-Maximization (EM) and Markov Chain Monte Carlo (MCMC) Bayesian estimation. Robustness is obtained either by Huber-type down-weighting or by mixtures of multivariate t distributions (a likelihood-based robust model, see Peel and McLachlan (2000) doi:10.1023/A:1008981510081). Missing data are handled by full information maximum likelihood with the exact EM treatment of incomplete observations (data augmentation in the MCMC engine). The EM engine also supports LASSO regularization with k-fold cross-validation for penalty tuning; the MCMC engine uses a Bayesian Lasso with Laplace priors, multiple chains, Gelman-Rubin/effective sample size diagnostics and the widely applicable information criterion. It supports six geometric variance-covariance models, along with functions for bootstrapped likelihood ratio tests (BLRT), BCH auxiliary variable analysis, and plotting. For longitudinal data, it fits robust growth mixture models and latent class growth analysis (Muthen and Shedden (1999) doi:10.1111/j.0006-341X.1999.00463.x) for one or several outcomes measured on unbalanced occasions, with Gaussian, Huber-weighted or multivariate-t (Pinheiro, Liu and Wu (2001) doi:10.1198/10618600152628059) latent classes, by EM and MCMC, and optional adaptive LASSO penalties (Zou (2006) doi:10.1198/016214506000000735) that identify stable trajectories and the outcomes that differentiate the classes. For methodological details on the Bootstrapped Likelihood Ratio Test, see Nylund et al. (2007) doi:10.1080/10705510701575396. For robust clustering methods, see Garcia-Escudero et al. (2010) doi:10.1007/s11634-010-0064-5. For BCH auxiliary variable analysis, see Bolck et al. (2004) doi:10.1093/pan/mph001.

Author(s)

Maintainer: Valerio Riccardo Aquila valerio_aquila@hotmail.it (ORCID)

Authors:

• Valerio Riccardo Aquila valerio_aquila@hotmail.it (ORCID)

BCH Method for Auxiliary Continuous Variables

Description

Applies the 3-step Bolck-Croon-Hagenaars (BCH) method to test the relationship between robust latent profiles and a continuous auxiliary (distal outcome) variable, adjusting for classification error in the profile assignments. Optionally adds a bootstrap correction (correction = "bootstrap") for the F-test's main weakness: treating the classification matrix D as known/fixed (see Details).

Usage

bch_robust(
  model,
  aux_var,
  correction = c("none", "bootstrap"),
  n_boot = 200,
  cores = 1
)

Arguments

Details

Let \hat{p}_{ig} be the posterior probability that observation i belongs to profile g (model$probabilities), and let \hat{C}_i be its modal (hard) assignment (model$assignments). The classification error matrix D is estimated as

D_{t,s} = P(\hat{C} = s \mid C = t) \approx \frac{\sum_{i:\, \hat{C}_i = s} \hat{p}_{it}}{\sum_{i=1}^{n} \hat{p}_{it}}

(Bolck, Croon, & Hagenaars, 2004; Vermunt, 2010), so each row of D sums to 1. The BCH weight matrix is W = D^{-1}, and the classification-error-corrected mean of the auxiliary variable Y for profile t is

\hat{\mu}^{BCH}_t = \frac{\sum_{i=1}^{n} W_{\hat{C}_i,t} Y_i}{\sum_{i=1}^{n} W_{\hat{C}_i,t}}

summed over every observation: each contributes to every profile's mean with a (possibly negative) cross-class weight W_{\hat{C}_i,t}, which is what removes the attenuation bias of a naive comparison of modal-assignment groups. This is unbiased under the BCH assumption that Y is independent of the modal assignment given the true profile.

The main ($ANOVA_Table) significance test is a one-way weighted ANOVA on the equivalent "long" data set (one row per observation per profile, weighted by W_{\hat{C}_i,g}), fit by direct weighted normal equations rather than stats::lm()/stats::aov(), because the BCH weights are frequently negative and base R's weighted-least-squares machinery cannot handle that.

The fixed-D caveat, and the bootstrap correction. $ANOVA_Table's F-test treats D as known/fixed. Bolck et al. (2004) and Vermunt (2010) both note that this understates the true uncertainty, because D is itself estimated from the step-1 model; Vermunt (2010) reports that naive (uncorrected) BCH p-values can be "much too small," particularly with poorly separated profiles or small samples. The literature's analytic fix is a "sandwich" (pseudo-likelihood) variance correction (Bakk, Oberski, & Vermunt, 2014), which requires the Fisher information of the step-1 mixture log-likelihood – intractable to derive analytically here for the Huber-robust EM and Bayesian-Lasso MCMC engines. Setting correction = "bootstrap" instead approximates that correction nonparametrically: it resamples observations with replacement, refits the entire step-1 robust_lpa model (using the exact same specification as model, via its stored $call_args) and recomputes D/W/the profile means on each resample, so the resulting bootstrap variability genuinely reflects step-1 estimation uncertainty (unlike the fixed-D F-test); for robust_gmm() fits whole persons (with all their occasions) are resampled (a cluster bootstrap). This yields bootstrap standard errors and percentile confidence intervals for Profile_Means, plus a Wald chi-square test of "all profile means equal" using the bootstrap covariance – reported in $Bootstrap_Correction, and preferable to $ANOVA_Table's p-value for publication-grade inference. It is not the Bakk et al. (2014) analytic formula; treat it as a practical approximation with the same goal (each refit's arbitrary profile labels are first aligned to model's by an exact optimal assignment of standardized profile means, to avoid mixing different real-world profiles together across resamples). Each refit runs sequentially (cores = 1) inside its bootstrap replicate, whatever cores the original model used, so parallelism happens only across replicates. It is off ("none") by default because it requires n_boot additional full model refits and is therefore substantially slower; use cores > 1 to parallelize it.

Value

A list containing:

Profile_Means

Named numeric vector of BCH bias-corrected profile means of aux_var.

ANOVA_Table

A data.frame with Df, Sum_Sq, Mean_Sq, F_value, and p_value for the "Class" and "Residuals" rows (see the fixed-D caveat in Details).

Classification_Matrix

The G x G classification error matrix D, with D[t, s] = P(assigned s | true profile t); rows sum to 1.

Classification_Weights

The G x G BCH weight matrix W = D^{-1} (rows: assigned profile; columns: true profile).

N_Used

Integer, the number of observations retained after removing missing aux_var values.

Bootstrap_Correction

NULL unless correction = "bootstrap", in which case a list with n_boot_used, n_boot_failed, SE and CI_lower/CI_upper (per profile), and the overall Wald_stat/Wald_df/Wald_p_value test of equal profile means (see Details).

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. doi:10.1093/pan/mph001

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

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. doi:10.1093/pan/mpu003

Examples

data(neuro_data)
# Fit the model on Memory and RT_Stroop only
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
set.seed(1)
fit <- robust_lpa(data = x, G = 2, model = 1, n_starts = 3)
summary(fit)  # profile means for the two fitted variables
# Test RT_TMT (not used to fit the model) as an auxiliary outcome
bch_res <- bch_robust(fit, neuro_data$RT_TMT)
bch_res$Profile_Means
bch_res$ANOVA_Table


# Add the bootstrap classification-uncertainty correction (slower: refits
# the model n_boot times). A small n_boot here is just for a fast demo --
# use several hundred for publication-grade inference.
bch_res_boot <- bch_robust(fit, neuro_data$RT_TMT, correction = "bootstrap", n_boot = 30)
bch_res_boot$Bootstrap_Correction

Bootstrapped Likelihood Ratio Test for Robust Growth Mixture Models

Description

Tests a growth mixture model with G classes against the same model with G - 1 classes by parametric bootstrap (McLachlan, 1987; Nylund, Asparouhov & Muthen, 2007): data are simulated from the fitted G - 1-class model – for the persons, occasions and outcomes actually observed, so the visit schedule and the missingness pattern are reproduced exactly – and both models are refitted to every simulated dataset. With robust_method = "t" the data are simulated from the fitted multivariate-t model and the test compares two proper t-mixture likelihoods; with the Huber estimator the statistic is a heuristic.

Usage

blrt_gmm_robust(
  data,
  id,
  time,
  outcomes,
  G,
  n_samples = 50,
  n_starts = 3,
  cores = 1,
  ...
)

Arguments

Value

A list with LRT_Observed, Bootstrap_LRTs, p_value (with the "+1" correction) and Bootstrap_Failures.

References

McLachlan, G. J. (1987). On bootstrapping the likelihood ratio test statistic for the number of components in a normal mixture. Journal of the Royal Statistical Society: Series C, 36(3), 318-324. doi:10.2307/2347790

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. doi:10.1080/10705510701575396

See Also

robust_gmm, estimate_gmm_robust

Examples

data(neuro_long)
set.seed(1)
blrt_gmm_robust(neuro_long, id = "ID", time = "Year", outcomes = "Memory",
                G = 2, n_samples = 5, n_starts = 2)

Bootstrapped Likelihood Ratio Test for Robust LPA

Description

Compares a robust LPA model with G profiles against a null model with G - 1 profiles using parametric bootstrapping (Nylund et al., 2007): the null model is fit to the observed data, data are simulated from it, and both the null and alternative models are refit to each simulated dataset to build a reference distribution for the likelihood ratio test statistic under H0. Supports FIML simulation conditions (the missingness pattern of the observed data is replicated in every simulated dataset).

Usage

blrt_robust(
  data,
  G,
  model = 6,
  engine = "EM",
  n_samples = 50,
  n_starts = 2,
  cores = 1,
  ...
)

Arguments

Details

Data are simulated from the fitted null model in the same family that is being fitted: a Gaussian mixture for classical and Huber fits, and a multivariate-t mixture (with the fitted nu) for robust_method = "t" fits. With robust_method = "t" the test compares two proper t-mixture likelihoods and is the recommended robust option. With the Huber estimator the statistic is computed from Gaussian log-likelihoods evaluated at robust (non-ML) estimates and the null data are simulated without contamination, so the test is a heuristic whose calibration on contaminated data is not guaranteed.

Value

A list containing:

LRT_Observed

The observed likelihood ratio test statistic (non-negative).

Bootstrap_LRTs

Numeric vector of the successfully-fit bootstrap replicates.

p_value

The empirical p-value, computed with the standard "+1" small-sample correction ((sum(Bootstrap_LRTs >= LRT_Observed) + 1) / (length(Bootstrap_LRTs) + 1)), which avoids reporting an (impossible) exact p-value of 0 from a finite bootstrap.

Bootstrap_Failures

Integer, how many of the n_samples bootstrap replicates failed to converge and were excluded from Bootstrap_LRTs / p_value.

References

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. doi:10.1080/10705510701575396

Examples

# Fast demonstration of the robust BLRT: is a 2nd profile justified over 1?
data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
set.seed(1)
blrt_res <- blrt_robust(x, G = 2, model = 1, n_samples = 5, n_starts = 3)
# Print the summary of the results
blrt_res

Estimate Robust Growth Mixture Models Across Numbers of Classes

Description

Fits robust_gmm for every number of classes in n_classes (and every random-effect specification in random), optionally selecting the LASSO penalties of each model by BIC or by K-fold cross-validation over persons, and collects the fit indices in one table.

Usage

estimate_gmm_robust(
  data,
  id,
  time,
  outcomes,
  n_classes = 1:3,
  random = "slope",
  cores = 1,
  tune_penalty = c("none", "bic", "cv"),
  penalty = c("both", "growth", "diff"),
  z_grid = c(1.5, 2, 2.5, 3, 4),
  k_folds = 5,
  ...
)

Arguments

Details

With tune_penalty = "bic" or "cv", the penalty levels are chosen from z_grid: each value z gives \lambda = z^2/N (N = number of persons). With the default adaptive weights of robust_gmm, a growth term (or class difference) is then set to zero approximately when its Wald statistic is below z in absolute value, which makes the grid directly interpretable (e.g. z = 2 or z = \sqrt{\log N}). "bic" fits every candidate on the full data and keeps the lowest BIC (whose parameter count accounts for the zeros and fusions); "cv" keeps the candidate with the highest held-out log-likelihood over k_folds folds of persons and refits it on the full data.

Value

A list with fit_table (one row per model), models (named "G<g>_<random>") and, when penalties are tuned, tuning (the criterion of every candidate).

See Also

robust_gmm, blrt_gmm_robust

Examples

data(neuro_long)
set.seed(1)
res <- estimate_gmm_robust(neuro_long, id = "ID", time = "Year",
                           outcomes = "Memory", n_classes = 1:2, n_starts = 2)
res$fit_table

Estimate Robust Latent Profile Models Across Profiles and Models

Description

Fits every combination of n_profiles and models with robust_lpa and collects their fit indices in one table.

Usage

estimate_profiles_robust(
  data,
  n_profiles = 1:3,
  models = c(1, 2, 3, 4, 5, 6),
  engine = "EM",
  cores = 1,
  n_starts = 5,
  lambda = 0,
  tune_lasso = FALSE,
  k_folds = 5,
  lambda_grid = c(0.01, 0.05, 0.1, 0.2),
  ...
)

Arguments

Value

A list containing fit_table (one row per successfully fitted model, with a Lambda column) and models (the fitted robust_lpa objects, named "model_<m>_profiles_<G>").

Examples

data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
set.seed(1)
res <- estimate_profiles_robust(x, n_profiles = 1:2, models = 1, n_starts = 3)
res$fit_table
summary(res$models[[1]])

Simulated Neuropsychological Dataset for Robust LPA

Description

A synthetic dataset of neuropsychological test scores and reaction times for two latent groups, "Healthy" and "Pathological", designed as a worked example for every estimation path in this package: the Expectation-Maximization and MCMC engines, all six variance-covariance parameterizations, robust vs. classical estimation, LASSO regularization, model/profile selection (estimate_profiles_robust), the bootstrapped likelihood ratio test (blrt_robust), and the BCH auxiliary-variable method (bch_robust).

Usage

neuro_data

Format

A data frame with 250 rows and 7 variables:

ID

Unique identifier for each participant.

True_Profile

The true latent group, "Healthy" (n = 150) or "Pathological" (n = 100). Not used for estimation (LPA is unsupervised); included so that recovered profiles can be checked against ground truth, e.g. table(neuro_data$True_Profile, fit$assignments).

Memory

Simulated memory test score. Differs in mean between groups.

Attention

Simulated attention test score. Identical distribution in both groups (no group signal); a noise variable.

Executive_Functions

Simulated executive functions score. Identical distribution in both groups (no group signal); a noise variable.

RT_Stroop

Reaction time in milliseconds. Differs in mean, variance, and correlation with RT_TMT between groups; a subset of Pathological observations carry an additional outlying shift.

RT_TMT

Reaction time in milliseconds. Differs in mean, variance, and correlation with RT_Stroop between groups; a subset of Pathological observations carry an additional outlying shift.

Details

The two groups differ in more than location: Attention and Executive_Functions have identical means and variances in both groups (deliberate noise variables, carrying no group signal), while Memory, RT_Stroop, and RT_TMT differ in mean between groups, and RT_Stroop/RT_TMT additionally differ in variance and in their correlation with each other (0.35 in Healthy vs. 0.90 in Pathological). This last feature is intentional: it is a genuine, whole-group difference in covariance structure (not merely in means), so that model = 6 (a fully unconstrained covariance matrix per profile) is the best-fitting parameterization for this dataset by BIC at G = 2 – run estimate_profiles_robust(scale(neuro_data[, 3:7]), models = 1:6, n_profiles = 1:3) and inspect $fit_table to see this directly. A more parsimonious model (e.g. model = 3, a single covariance matrix shared across profiles) fits these data measurably worse, illustrating why the six parameterizations exist and how to choose among them.

On top of this, 6% of the Pathological observations (chosen at random) receive an additional, positive, randomly-sized shift on RT_Stroop and RT_TMT (drawn from a Gamma distribution, so the contamination varies in severity rather than landing on a single fixed value) – measurement-error-like outliers on top of the two groups' otherwise multivariate-normal structure. These are what robust = TRUE (the default of robust_lpa) down-weights via Huber-type estimation; compare robust = TRUE vs. robust = FALSE fits to see their effect on the estimated Pathological-profile covariance.

The contamination magnitude and rate were calibrated (by direct grid search across the six variance-covariance models, replicated over multiple random seeds) so that fitting model = 6 with G = 2 reliably wins by BIC over both more-parsimonious models at G = 2 and less-parsimonious models at G = 3, and recovers True_Profile with about 93% accuracy (classical, Huber or t estimation, model = 6, G = 2).

Source

Simulated data for testing and documentation purposes.

Simulated Longitudinal Neuropsychological Dataset for Robust Growth Mixture Models

Description

A synthetic longitudinal dataset (long format: one row per person and annual visit) designed as a worked example for robust_gmm: three latent classes of cognitive change, three outcomes, unbalanced follow-up with missed visits and drop-out, and a few gross data-entry errors.

Usage

neuro_long

Format

A data frame with one row per person and visit, and 8 variables:

ID

Person identifier.

Year

Years since baseline (0 to 5).

Memory, Executive, Speed

Simulated test scores (NA when not observed).

True_Class

The true latent class (not used for estimation).

Age

Age at baseline (constant within person).

Biomarker

A baseline biomarker level (constant within person).

Details

400 persons are assigned to three latent classes ("Stable", "Slow decline", "Fast decline"; about 50/30/20%) and followed for up to six annual visits (Year 0 to 5). Memory and Executive (T-score-like metric) decline at class-specific rates (Memory: 0, -2 and -5 points per year; Executive: 0, -1.5 and -4 points per year), whereas Speed declines by 0.5 points per year in every class (an outcome that does not differentiate the classes, useful to illustrate the group LASSO of robust_gmm(lambda_diff = , group_diff = TRUE); the "Stable" class has exactly zero slopes on Memory and Executive, useful to illustrate lambda_growth). Within classes, persons have correlated random intercepts (SD 5) and slopes (SD 0.4) on every outcome, and residual errors with SD 2.5. After each visit a person drops out with a probability that increases as the last observed Memory score decreases (missing at random); 8% of the follow-up visits are missed and 5% of the single test scores are missing. For 4% of the persons one score is corrupted by a gross error of 25 to 40 points. Age and Biomarker are baseline characteristics that differ between the classes, for illustrating bch_robust on a robust_gmm() fit. The generating code (with its fixed seed) is installed with the package: run source(system.file("scripts", "generate_neuro_long.R", package = "RobustLPA")) to rebuild the data set.

Source

Simulated; see Details.

Examples

data(neuro_long)
head(neuro_long)
table(table(neuro_long$ID))   # number of visits per person

Plot MCMC Trace for Robust LPA Models

Description

Draws multi-chain trace plots for the MCMC engine of robust_lpa using bayesplot. By default, the profile means, profile variances (the diagonal of each covariance matrix), and mixing proportions are shown; use pars to select a subset.

Usage

plot_mcmc_chains(model, pars = NULL)

Arguments

Details

mcmc_chain_cpp (called internally by robust_lpa(engine = "MCMC"), once per chain) returns its draws as a nested list (mu_chain, sigma_chain, pi_chain), not the array/matrix format bayesplot::mcmc_trace() expects. This function reshapes the raw per-chain draws stored in model$mcmc_draws$chains into an [iterations, chains, parameters] array before calling bayesplot::mcmc_trace(), so every chain is shown overlaid as a separate colored trace – the standard visual convergence check (well-mixed, overlapping chains suggest convergence; chains that stay visually separated suggest they have not converged, consistent with a high Gelman-Rubin \hat{R}; see model$mcmc_diagnostics).

Parameter names follow the pattern "mu[g,j]" (mean of variable j in profile g), "sigma[g,j]" (variance of variable j in profile g), "pi[g]" (mixing proportion of profile g), and, for robust_method = "t" fits with an estimated nu, "nu" (the t degrees of freedom). Off-diagonal covariance terms are not included by default to keep the default plot readable; inspect model$mcmc_draws$chains[[1]]$sigma_chain directly if you need those.

Value

A ggplot object generated by bayesplot::mcmc_trace().

Examples

data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
fit <- suppressWarnings(robust_lpa(x, G = 2, model = 2, engine = "MCMC",
                                    mcmc_iter = 200, n_chains = 2))
summary(fit)  # profile means and MCMC convergence diagnostics
plot_mcmc_chains(fit)

Plot the Class Trajectories of a Robust Growth Mixture Model

Description

Draws, for every outcome, the estimated mean trajectory of each latent class over the observed time range, optionally over the observed individual trajectories coloured by modal class.

Usage

plot_robust_gmm(
  model,
  outcomes = NULL,
  individual = TRUE,
  max_individuals = 300,
  title = "Latent class trajectories"
)

Arguments

Value

A ggplot object.

See Also

robust_gmm

Examples

data(neuro_long)
set.seed(1)
fit <- robust_gmm(neuro_long, id = "ID", time = "Year", outcomes = "Memory",
                  G = 2, n_starts = 2)
plot_robust_gmm(fit)

Plot Robust Latent Profiles

Description

Automatically generates a professional profile plot using ggplot2 from an estimated robust LPA model.

Usage

plot_robust_lpa(
  model,
  which_model = NULL,
  title = "Robust Latent Profiles",
  xlab = "Variables",
  ylab = "Value",
  var_labels = NULL,
  legend_title = "Class"
)

Arguments

Value

A ggplot object.

Examples

data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
fit <- suppressWarnings(robust_lpa(data = x, G = 2, model = 1, n_starts = 3))
print(fit)  # concise overview (print.robust_lpa())
plot_robust_lpa(fit)

Print a Fitted Robust Growth Mixture Model

Description

Compact overview of a model fitted by robust_gmm: engine, number of classes and persons, outcomes, random-effect and robustness settings, headline fit indices, class proportions and, if used, the LASSO penalties.

Usage

## S3 method for class 'robust_gmm'
print(x, ...)

Arguments

Value

x, invisibly.

See Also

robust_gmm, summary.robust_gmm

Print a Fitted Robust LPA Model

Description

A short, four-line-or-fewer overview of a model fitted by robust_lpa: engine/model/profiles/N, the headline fit indices (log-likelihood, AIC, BIC, entropy), the mixing proportions, and, for the MCMC engine, the chain configuration. It deliberately omits profile means and full diagnostics – use summary.robust_lpa for those.

Usage

## S3 method for class 'robust_lpa'
print(x, ...)

Arguments

Value

x, invisibly.

See Also

robust_lpa, summary.robust_lpa

Examples

data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
fit <- robust_lpa(x, G = 2, model = 2, n_starts = 3, max_iter = 20)
print(fit)

Print a Summarized Robust Growth Mixture Model

Description

Print a Summarized Robust Growth Mixture Model

Usage

## S3 method for class 'summary.robust_gmm'
print(x, digits = 3, ...)

Arguments

Value

x, invisibly.

Print a Summarized Robust LPA Model

Description

Prints the object returned by summary.robust_lpa: profile means, profile sizes/mixing proportions, the headline fit indices, and, for the MCMC engine, the chain configuration and Gelman-Rubin \hat{R} / effective sample size ranges.

Usage

## S3 method for class 'summary.robust_lpa'
print(x, digits = 2, ...)

Arguments

Value

x, invisibly.

Robust Growth Mixture Model (Latent Class Growth Analysis)

Description

Fits a growth mixture model (GMM) – a finite mixture of linear mixed-effects models for longitudinal data (Verbeke & Lesaffre, 1996; Muthen & Shedden, 1999) – or, with random = "none", a latent class growth analysis (LCGA; Nagin, 1999), with one or several outcomes measured repeatedly on the same persons. Each latent class ("profile of change") has its own polynomial mean trajectory for every outcome; within a class, persons deviate from it through correlated random effects and residual errors. The model can be estimated by maximum likelihood (EM engine) or by Gibbs sampling (MCMC engine), in a classical (Gaussian), Huber-weighted or multivariate-t (robust) version, optionally with LASSO penalties on the class trajectories.

Usage

robust_gmm(
  data,
  id,
  time,
  outcomes,
  G,
  degree = 1,
  random = c("slope", "intercept", "none"),
  re_structure = c("full", "block", "diagonal"),
  re_cov = c("equal", "varying"),
  resid_var = c("equal", "varying"),
  engine = "EM",
  robust = TRUE,
  robust_method = c("huber", "t"),
  nu = NULL,
  alpha = 0.05,
  lambda_growth = 0,
  lambda_diff = 0,
  group_diff = FALSE,
  adaptive = TRUE,
  relax = FALSE,
  n_starts = 5,
  max_iter = 500,
  tol = 1e-08,
  init = c("kmeans", "random"),
  mcmc_iter = 2000,
  n_chains = 4,
  cores = 1
)

Arguments

Value

An object of class "robust_gmm", a list with:

engine, robust_method, nu

Estimation settings (nu: the estimated or fixed t degrees of freedom, NULL otherwise).

coefficients

A list of G matrices (outcomes x polynomial terms): the class mean trajectories.

random_cov

A list of G random-effect covariance matrices (NULL for random = "none").

residual_var

A G x K matrix of residual variances.

proportions

Class proportions.

probabilities, assignments

Posterior class probabilities (persons x classes) and modal classes, in the order of ids.

ids

The person identifiers, in the order used by every person-level output (use it to align auxiliary variables for bch_robust).

weights

Per-person robustness weights, averaged over classes with the posterior probabilities. With robust_method = "huber" they lie in (0, 1] (1 = full weight). With robust_method = "t" they are the E-step weights (\nu + n_i)/(\nu + \delta_i), where n_i is the number of observed values and \delta_i the squared Mahalanobis distance of person i: they exceed 1 for persons closer to their class trajectory than expected and are small for outlying persons. All 1 when robust = FALSE.

random_effects

Posterior means of each person's random effects under his/her modal class.

fit

One-row data.frame: Classes, LogLik, Parameters, AIC, BIC (with \log N, N = persons), SABIC, Entropy, Min_Size, Max_Size, the penalty levels, and WAIC (MCMC).

penalty

Penalty settings and the selected zeros / fusions.

converged, iterations

EM convergence information.

mcmc_draws, mcmc_diagnostics, waic

MCMC output (see plot_mcmc_chains).

spec, internal, call_args, data

Model specification, standardized-scale estimates, arguments and data used for refitting.

Model

For person i in class g, the stacked vector y_i of all his/her observed values (all outcomes, all occasions) is

y_i = X_i \beta_g + Z_i b_i + e_i, \quad b_i \sim N(0, D_g / u_i), \quad e_i \sim N(0, R_{g} / u_i),

where X_i contains, for each observation of outcome k at time t, the polynomial basis (1, t, \dots, t^{degree}) in the block of outcome k; Z_i the first 1 (random = "intercept") or 2 (random = "slope") of those terms; R_g is diagonal with one residual variance per outcome; and u_i = 1 (Gaussian) or u_i \sim \mathrm{Gamma}(\nu/2, \nu/2) (multivariate t: Pinheiro, Liu & Wu, 2001). Marginally, y_i is Gaussian or multivariate t with mean X_i \beta_g and scale V_{ig} = Z_i D_g Z_i' + R_g.

Persons may be observed at different times, on different numbers of occasions, and not on every outcome at every occasion: each person contributes the exact likelihood of the values actually observed, so the estimates are maximum likelihood under missing-at-random missingness (including drop-out that depends on earlier observed values). Every outcome is standardized internally (all reported estimates are on the original scale); the time variable is used as supplied, so center it where the intercept should be interpreted (e.g. years since baseline).

Robust estimation

robust_method = "t"

A mixture of multivariate-t linear mixed models: the random effects and the errors of a person share the latent scale u_i, so a person whose trajectory is far from every class (or who has a few gross errors) is down-weighted as a whole through E[u_i] = (\nu + n_i) / (\nu + d_i), with d_i the Mahalanobis distance of y_i and n_i its length. It is a proper likelihood model, fitted by ECM with an ECME update of \nu, so information criteria, the bootstrapped likelihood ratio test (blrt_gmm_robust) and the BCH method are used as intended. Recommended.

robust_method = "huber"

Huber weights computed from d_i with a \chi^2_{n_i} cutoff (alpha) down-weight outlying persons in the M-step; an estimating-equation approach whose reported log-likelihood is the Gaussian one evaluated at the robust estimates (a heuristic for information criteria).

LASSO penalties

Two penalties on the class trajectories (fixed effects) are available, separately or together; both act on the standardized outcome scale and are expressed per person, i.e. the EM engine maximizes

\ell(\theta)/N - \lambda_{growth} \sum_{g,k,l \ge 1} |\beta_{gkl}| - \lambda_{diff} \sum_{g,k,l} |\beta_{gkl} - \bar\beta_{kl}|,

with \bar\beta_{kl} the (unweighted) mean over classes.

lambda_growth

Sparse trajectories: the growth terms (slope, quadratic, ...) of every class and outcome are shrunk toward 0, intercepts are not penalized. A coefficient set exactly to 0 means that the class does not change on that outcome (e.g. a "stable" class). This is the Lasso of Tibshirani (1996) applied to the fixed effects of a mixture of regressions (Khalili & Chen, 2007; Du et al., 2013).

lambda_diff

Sparse class differences: every coefficient of every class is shrunk toward the across-class mean, so that coefficients (or, with group_diff = TRUE, whole outcomes) on which the classes do not differ are fused. With group_diff = TRUE the penalty is \lambda_{diff} \sum_k w_k \| S_k (B_k - 1\bar\beta_k') \|, where S_k scales every coefficient by the square root of its information (so that intercepts and slopes are penalized on comparable scales; Simon & Tibshirani, 2012) and w_k is the adaptive weight; it removes outcome k from the class separation altogether when its group is set to zero – the grouped variable-selection penalty of Xie, Pan & Shen (2008), which extends the penalized model-based clustering of Pan & Shen (2007) used by robust_lpa's lambda.

By default (adaptive = TRUE) the penalties are adaptive (Zou, 2006; Wang & Leng, 2008 for the group version): every coefficient, deviation or outcome group is weighted by the reciprocal of its unpenalized estimate, which gives consistent selection and makes the penalty level interpretable – a term is set to zero approximately when its Wald statistic is below \sqrt{N\lambda} in absolute value, so lambda = z^2 / N corresponds to a threshold z (e.g. 4 / N for |z| < 2). The penalized fit is started from the unpenalized maximum-likelihood solution (which also supplies the weights). The penalized fixed-effects step of each iteration is solved exactly (to numerical tolerance) by ADMM (Boyd et al., 2011); the number of fixed-effect parameters used by AIC/BIC is the number of free coefficients left by the zeros and fusions (the generalized-Lasso degrees of freedom; Tibshirani & Taylor, 2011). Penalized estimates are shrunk toward zero / toward each other: relax = TRUE refits the model without penalty while keeping the selected zeros and fusions (the relaxed Lasso; Meinshausen, 2007), which is preferable for reporting and inference. Choose lambda_growth / lambda_diff with estimate_gmm_robust (BIC or cross-validation over persons). In the MCMC engine the same penalties become Bayesian-Lasso (Laplace) priors with rate N\lambda times the adaptive weight of each term (Park & Casella, 2008; Bayesian group Lasso of Kyung et al., 2010, when group_diff = TRUE), whose posterior mode is the EM penalized estimate; relax does not apply to the MCMC engine.

Estimation

The EM engine uses an alternating ECM algorithm (Meng & van Dyk, 1997): each iteration updates the class proportions and the class trajectories by (penalized) weighted generalized least squares with the random effects integrated out, then, after a fresh E-step, the residual variances, the random-effect covariances and (t model) \nu; every step increases the (penalized) observed-data log-likelihood. The iterations are accelerated by SQUAREM (Varadhan & Roland, 2008) with a monotonicity safeguard. Several starts are run (k-means on person-level least-squares trajectories by default) and the best is kept.

MCMC details

The Gibbs sampler allocates persons with the random effects and latent scales integrated out, then draws \nu (t model; random-walk Metropolis with the latent scales integrated out, Gamma(2, 0.1) prior), the latent scales, the class trajectories with the random effects integrated out, the random effects, the random-effect covariance matrices, the residual variances (inverse-gamma(1, 0.1) on the standardized scale) and the mixing proportions (Dirichlet(1, ..., 1)). The random-effect covariances have a scaled inverse-Wishart prior (O'Malley & Zaslavsky, 2008), D = \mathrm{diag}(a) \Psi \mathrm{diag}(a) with \Psi \sim IW(q + 1, I) per block and a_j \sim N(0, A_j^2) (A_j = 5 standardized units for intercepts, divided by the standard deviation of the time term for slopes), sampled by parameter expansion (Liu & Wu, 1999; Gelman et al., 2008), which avoids the slow mixing of small variance components. The fixed effects have vague normal priors, combined with the Bayesian-Lasso priors when penalties are requested. The chains start from the EM fit with perturbed trajectories, draws are relabeled to that EM solution by an optimal assignment of the class mean trajectories, and the first half of each chain is discarded. The Huber option is a heuristic in the MCMC engine, as in robust_lpa.

References

Boyd, S., Parikh, N., Chu, E., Peleato, B., & Eckstein, J. (2011). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning, 3(1), 1-122. doi:10.1561/2200000016

Du, Y., Khalili, A., Neslehova, J. G., & Steele, R. J. (2013). Simultaneous fixed and random effects selection in finite mixture of linear mixed-effects models. Canadian Journal of Statistics, 41(4), 596-616. doi:10.1002/cjs.11192

Gelman, A., van Dyk, D. A., Huang, Z., & Boscardin, W. J. (2008). Using redundant parameterizations to fit hierarchical models. Journal of Computational and Graphical Statistics, 17(1), 95-122. doi:10.1198/106186008X287337

Khalili, A., & Chen, J. (2007). Variable selection in finite mixture of regression models. Journal of the American Statistical Association, 102(479), 1025-1038. doi:10.1198/016214507000000590

Kyung, M., Gill, J., Ghosh, M., & Casella, G. (2010). Penalized regression, standard errors, and Bayesian lassos. Bayesian Analysis, 5(2), 369-411. doi:10.1214/10-BA607

Liu, J. S., & Wu, Y. N. (1999). Parameter expansion for data augmentation. Journal of the American Statistical Association, 94(448), 1264-1274. doi:10.1080/01621459.1999.10473879

Meinshausen, N. (2007). Relaxed Lasso. Computational Statistics & Data Analysis, 52(1), 374-393. doi:10.1016/j.csda.2006.12.019

Meng, X.-L., & van Dyk, D. (1997). The EM algorithm – an old folk-song sung to a fast new tune. Journal of the Royal Statistical Society: Series B, 59(3), 511-567. doi:10.1111/1467-9868.00082

Muthen, B., & Shedden, K. (1999). Finite mixture modeling with mixture outcomes using the EM algorithm. Biometrics, 55(2), 463-469. doi:10.1111/j.0006-341X.1999.00463.x

Nagin, D. S. (1999). Analyzing developmental trajectories: A semiparametric, group-based approach. Psychological Methods, 4(2), 139-157. doi:10.1037/1082-989X.4.2.139

O'Malley, A. J., & Zaslavsky, A. M. (2008). Domain-level covariance analysis for multilevel survey data with structured nonresponse. Journal of the American Statistical Association, 103(484), 1405-1418. doi:10.1198/016214508000000724

Pan, W., & Shen, X. (2007). Penalized model-based clustering with application to variable selection. Journal of Machine Learning Research, 8, 1145-1164.

Park, T., & Casella, G. (2008). The Bayesian Lasso. Journal of the American Statistical Association, 103(482), 681-686. doi:10.1198/016214508000000337

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. doi:10.1198/10618600152628059

Simon, N., & Tibshirani, R. (2012). Standardization and the group lasso penalty. Statistica Sinica, 22(3), 983-1001.

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B, 58(1), 267-288. doi:10.1111/j.2517-6161.1996.tb02080.x

Tibshirani, R. J., & Taylor, J. (2011). The solution path of the generalized lasso. The Annals of Statistics, 39(3), 1335-1371. doi:10.1214/11-AOS878

Varadhan, R., & Roland, C. (2008). Simple and globally convergent methods for accelerating the convergence of any EM algorithm. Scandinavian Journal of Statistics, 35(2), 335-353. doi:10.1111/j.1467-9469.2007.00585.x

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. doi:10.1080/01621459.1996.10476679

Wang, H., & Leng, C. (2008). A note on adaptive group lasso. Computational Statistics & Data Analysis, 52(12), 5277-5286. doi:10.1016/j.csda.2008.05.006

Xie, B., Pan, W., & Shen, X. (2008). Variable selection in penalized model-based clustering via regularization on grouped parameters. Biometrics, 64(3), 921-930. doi:10.1111/j.1541-0420.2007.00955.x

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

See Also

estimate_gmm_robust (number of classes and penalty selection), blrt_gmm_robust, plot_robust_gmm, bch_robust (relating the classes to baseline or distal variables), plot_mcmc_chains.

Examples

data(neuro_long)
set.seed(1)
fit <- robust_gmm(neuro_long, id = "ID", time = "Year",
                  outcomes = c("Memory", "Executive"), G = 2,
                  robust_method = "t", n_starts = 2)
fit
summary(fit)
plot_robust_gmm(fit)


# Sparse class differences (group LASSO by outcome), relaxed refit
# (adaptive weights: lambda = z^2 / N removes groups with |z| below ~z)
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", lambda_diff = 9 / N,
                    group_diff = TRUE, relax = TRUE, n_starts = 2)
fit_l$penalty$outcome_selected

# MCMC engine
fit_b <- robust_gmm(neuro_long, id = "ID", time = "Year",
                    outcomes = "Memory", G = 2, engine = "MCMC",
                    robust_method = "t", mcmc_iter = 400, n_chains = 2)
summary(fit_b)

Fit a Single Robust Latent Profile Analysis Model

Description

Estimates a Latent Profile Analysis (finite mixture) model that is robust to multivariate outliers and handles missing data by full-information maximum likelihood, using either an EM or an MCMC (Bayesian Lasso) engine. Two robust estimators are available (see "Robust estimation" below): Huber-type down-weighting (robust_method = "huber", the default) and a mixture of multivariate t distributions (robust_method = "t"), which is a proper likelihood-based model. The MCMC engine runs n_chains independent chains (4 by default) and reports Gelman-Rubin \hat{R}, effective sample size and WAIC. Set cores > 1 to run the EM engine's random restarts, or the MCMC engine's chains, in parallel.

Usage

robust_lpa(
  data,
  G,
  model = 6,
  engine = "EM",
  max_iter = 100,
  tol = 1e-06,
  n_starts = 5,
  lambda = 0,
  mcmc_iter = 2000,
  prior_laplace = 0.1,
  robust = TRUE,
  alpha = 0.05,
  n_chains = 4,
  cores = 1,
  robust_method = c("huber", "t"),
  nu = NULL,
  init = c("kmeans", "random")
)

Arguments

Value

A list with S3 class "robust_lpa" (see print.robust_lpa and summary.robust_lpa) containing:

engine

The estimation engine used.

robust_method

"none", "huber" or "t".

means

A list of length G with the estimated profile means.

covariances

A list of length G with the estimated profile covariance matrices (for robust_method = "t": the scale matrices; the covariance of a t component is nu / (nu - 2) times the scale matrix for nu > 2).

proportions

Numeric vector of length G with the estimated mixing proportions.

nu

The (estimated or fixed) t degrees of freedom; NULL unless robust_method = "t". Posterior median for the MCMC engine.

probabilities

An n x G matrix of posterior profile-membership probabilities.

weights

Numeric vector of length n: each observation's robustness weight, averaged over profiles with its posterior probabilities; small values flag outliers. Huber weights lie in (0, 1] (1 = full weight). With robust_method = "t" they are the E-step weights (\nu + p_i)/(\nu + \delta_i) (p_i observed variables, \delta_i squared Mahalanobis distance), which exceed 1 for observations close to the profile mean. All 1 when robust = FALSE.

fit

A one-row data.frame with Model, Profiles, LogLik, Parameters, AIC, BIC, SABIC, Entropy, Min_Size, Max_Size and, for the MCMC engine, WAIC. Parameters counts every free parameter; only with the EM engine and lambda > 0 are the mean components shrunk exactly to zero left out.

assignments

Integer vector of length n with the most likely profile for each observation.

converged, iterations

(EM engine) whether the selected start met tol before max_iter, and how many iterations it used.

mcmc_draws

(MCMC engine only) a list with chains (the raw per-chain draws, relabeled to a common profile ordering), n_chains, mcmc_iter, burnin, and nu_estimated.

mcmc_diagnostics

(MCMC engine only) a data.frame with one row per scalar parameter (Parameter, Rhat, ESS).

waic

(MCMC engine only) a list with WAIC, lppd and p_waic.

data

The numeric matrix actually fit.

call_args

A named list of every argument controlling this fit, used to refit the identical specification on new data (e.g. by bch_robust's bootstrap correction).

Robust estimation

robust_method = "huber"

At every M-step each observation's contribution to a profile's mean and covariance is down-weighted by a Huber weight computed from its squared Mahalanobis distance to that profile's current (already robust) estimates, so the down-weighting accumulates across iterations and converges to an iteratively reweighted M-estimator. This is an estimating-equation approach, not a likelihood: the reported LogLik (and hence AIC/BIC/SABIC and blrt_robust) is the Gaussian mixture log-likelihood evaluated at the robust estimates, which is a useful heuristic but not the maximized objective. Huber-type estimators also have a limited breakdown point: a large, compact cluster of outliers can still mask itself.

robust_method = "t"

Each profile is a multivariate t distribution (McLachlan & Peel, 1998; Peel & McLachlan, 2000), whose heavier tails automatically down-weight outlying observations through the latent-scale weights (nu + p) / (nu + d). This is a proper likelihood-based model fitted by ECME, so LogLik, AIC/BIC, the BLRT and the BCH method are all used as intended; it is also markedly more resistant than Huber weighting to gross outliers. nu is estimated by default and counts as one extra parameter.

In the MCMC engine, "t" is implemented as an exact sampler (latent Gamma scales; nu updated with the latent scales integrated out), whereas "huber" remains a heuristic that down-weights sufficient statistics and does not target a well-defined posterior; prefer "t" for Bayesian inference.

Missing data

Rows with missing values contribute their exact observed-data likelihood (the marginal density of their observed entries) to the E-step, and the M-step uses the exact EM treatment of incomplete multivariate data: missing entries are replaced by their conditional expectations given the observed ones and the conditional covariance is added to the scatter matrix (Ghahramani & Jordan, 1994; Liu & Rubin, 1995). The resulting estimates are maximum likelihood under missing-at-random (MAR) missingness. The MCMC engine uses the equivalent data-augmentation step (missing entries are drawn from their conditional distribution every sweep). Rows with no observed variables are allowed and are allocated according to the mixing proportions only.

MCMC details

The chains are initialized from a preliminary EM fit (same model and robust method), with the profile means of each chain independently perturbed by half a within-profile standard deviation so that the chains start from dispersed points (as the Gelman-Rubin diagnostic assumes). Priors: Laplace (Bayesian Lasso) on the means, Dirichlet(1, ..., 1) on the mixing proportions, inverse-gamma(1, 1) on variances (models 1, 2 and the scales of models 4-5), inverse-Wishart(p + 1, I) on full covariance matrices (models 3 and 6), a uniform distribution on correlation matrices (models 4-5), and Gamma(2, 0.1) on nu. All conditionals are sampled exactly except the correlation matrices of models 4-5, which use random-walk Metropolis steps (step sizes adapted during burn-in only). Label switching is resolved by relabeling every draw to the EM solution (pivotal reordering) with an exact optimal assignment on standardized mean distances. $mcmc_diagnostics reports the Gelman-Rubin \hat{R} and the effective sample size of every scalar parameter; $fit$WAIC is the widely applicable information criterion (Watanabe, 2010), computed from up to 1000 posterior draws and preferable to the plug-in AIC/BIC for comparing Bayesian fits.

Reproducibility

Every random quantity (EM starts, MCMC chains, bootstrap replicates) is generated from a seed drawn from R's random number stream before any parallel dispatch, so set.seed(1); robust_lpa(..., cores = 1) and set.seed(1); robust_lpa(..., cores = 4) give identical results on the same machine. Across operating systems, compilers or linear-algebra libraries, results may differ in the last digits (and, when two solutions are almost equally good, in the selected start or the order of the profiles).

References

Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457-472. doi:10.1214/ss/1177011136

Ghahramani, Z., & Jordan, M. I. (1994). Supervised learning from incomplete data via an EM approach. Advances in Neural Information Processing Systems, 6, 120-127.

Liu, C., & Rubin, D. B. (1995). ML estimation of the t distribution using EM and its extensions, ECM and ECME. Statistica Sinica, 5(1), 19-39.

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

Park, T., & Casella, G. (2008). The Bayesian Lasso. Journal of the American Statistical Association, 103(482), 681-686. doi:10.1198/016214508000000337

Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571-3594.

See Also

estimate_profiles_robust to fit and compare many G / model combinations at once, blrt_robust for a bootstrapped likelihood ratio test, bch_robust to relate profiles to distal outcomes, and plot_mcmc_chains to inspect MCMC chains.

Examples

data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "Attention", "Executive_Functions",
                                     "RT_Stroop", "RT_TMT")]))
set.seed(1)
fit <- robust_lpa(x, G = 2, model = 6, n_starts = 2, max_iter = 50)
fit           # print.robust_lpa(): concise overview of the fit
summary(fit)  # summary.robust_lpa(): profile means/sizes and fit indices

# Multivariate-t mixture: a likelihood-based robust alternative
fit_t <- robust_lpa(x, G = 2, model = 6, n_starts = 2, max_iter = 50,
                    robust_method = "t")
fit_t$nu
head(sort(fit_t$weights))  # smallest weights = most outlying observations

# Compare with classical (non-robust) estimation
fit_classical <- robust_lpa(x, G = 2, model = 6, n_starts = 2, max_iter = 50,
                            robust = FALSE)
sapply(fit_t$means, `[`, "RT_Stroop")
sapply(fit_classical$means, `[`, "RT_Stroop")


# MCMC engine: 4 chains (default), with a small mcmc_iter for speed
fit_mcmc <- robust_lpa(x, G = 2, model = 6, engine = "MCMC",
                       robust_method = "t", mcmc_iter = 300, n_chains = 4)
summary(fit_mcmc)  # includes Rhat / ESS ranges and WAIC

Auxiliary M-Step Function for Robust Estimation

Description

Computes a single profile's updated mean and (unconstrained) covariance matrix from its posterior membership probabilities, with optional robustness weighting and optional soft-thresholding (LASSO-type shrinkage) of the mean. Used internally by robust_lpa at every EM iteration, for every profile; not intended to be called directly by end users.

Usage

robust_m_step(
  data,
  z,
  alpha = 0.05,
  lambda = 0,
  robust = TRUE,
  mu = NULL,
  sigma = NULL,
  robust_method = c("huber", "t"),
  nu = 4,
  weights = NULL,
  patterns = NULL
)

Arguments

Value

A list with mean (numeric row vector), covariance (a symmetric matrix) and weights (the robustness weights used).

Robustness weights

The robustness weights are computed from each observation's squared Mahalanobis distance (on its observed entries) to the profile's current mean/covariance, i.e. the parameters that produced the posterior probabilities z. Because these current parameters are themselves the robust estimates from the previous iteration, the down-weighting compounds across EM iterations and converges to a fixed point (an iteratively reweighted M-estimator), instead of being recomputed each time from a non-robust starting point (which lets outliers mask themselves by inflating the covariance they are measured against).

robust_method = "huber"

Huber weights: 1 inside the 1 - alpha chi-squared quantile (with degrees of freedom equal to the number of observed variables of the row), sqrt(cutoff / d) beyond it. The covariance is the Huber-weighted scatter matrix.

robust_method = "t"

The E-step expectation of the latent scale of a multivariate-t distribution, (nu + p_obs) / (nu + d); together with the covariance update below this is the exact ECM M-step of a multivariate-t mixture (McLachlan & Peel, 2000).

Missing data

Missing entries are handled by the exact EM treatment of incomplete multivariate data (Ghahramani & Jordan, 1994; Liu & Rubin, 1995): each missing block is replaced by its conditional expectation given the observed entries under the profile's current parameters, and the corresponding conditional covariance is added to the scatter matrix. This gives maximum-likelihood estimates under missing-at-random (MAR) missingness. When no current parameters are supplied (initialization), pairwise available-case moments are used as the starting point for that single step.

References

Ghahramani, Z., & Jordan, M. I. (1994). Supervised learning from incomplete data via an EM approach. Advances in Neural Information Processing Systems, 6, 120-127.

Liu, C., & Rubin, D. B. (1995). ML estimation of the t distribution using EM and its extensions, ECM and ECME. Statistica Sinica, 5(1), 19-39.

McLachlan, G. J., & Peel, D. (2000). Finite Mixture Models. Wiley. doi:10.1002/0471721182

Calculate a Simple Robust (Trimmed) Mean

Description

Computes a one-step trimmed centroid: an observation is included in the average only if its Euclidean distance to the coordinate-wise median of data is below threshold. This is a quick, easy-to-reason-about robust location estimate, not an iterative M-estimator; for the full robust mixture-model estimation used elsewhere in this package, see robust_lpa.

Usage

robust_mean(data, threshold = 10)

Arguments

Details

threshold is a distance from the coordinate-wise median of data, not from the origin – so a sensible value depends on the scale and spread of your variables. A reasonable starting point is a small multiple of a typical per-variable standard deviation times sqrt(ncol(data)) (roughly the scale of a Euclidean distance across all variables); mahalanobis-based thresholds (as used internally by robust_lpa) account for correlation and scale automatically and are preferable when variables are on very different scales.

Value

A numeric vector representing the robust mean of the variables. If no observation falls within threshold of the median, returns a vector of zeros with a warning.

Examples

data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
r_mean <- robust_mean(x, threshold = 3)

# Print the calculated robust means
r_mean

Summarize a Fitted Robust Growth Mixture Model

Description

Summarize a Fitted Robust Growth Mixture Model

Usage

## S3 method for class 'robust_gmm'
summary(object, ...)

Arguments

Value

An object of class "summary.robust_gmm": a list with the class trajectories (a data.frame with one row per class and outcome), class sizes, variance components, fit indices, penalty information and, for the MCMC engine, the range of the convergence diagnostics.

See Also

robust_gmm

Examples

data(neuro_long)
set.seed(1)
fit <- robust_gmm(neuro_long, id = "ID", time = "Year", outcomes = "Memory",
                  G = 2, n_starts = 2)
summary(fit)

Summarize a Fitted Robust LPA Model

Description

Builds a compact summary of a model fitted by robust_lpa, limited to the information most people actually need to interpret a fit: per-profile means, profile sizes/mixing proportions, the headline fit indices, and, for the MCMC engine, the Gelman-Rubin \hat{R} / effective sample size convergence range. Returns an object of class "summary.robust_lpa" with its own print method (print.summary.robust_lpa), following the usual summary()/print(summary()) convention used throughout R (e.g. summary.lm). For the full per-parameter Rhat/ESS table, use object$mcmc_diagnostics directly.

Usage

## S3 method for class 'robust_lpa'
summary(object, ...)

Arguments

Value

An object of class "summary.robust_lpa", a list with engine, method (a description of the estimation method), converged/iterations (EM engine), model, G, n, means (a variables x profiles matrix), sizes (a data.frame of profile sizes/mixing proportions), fit (the one-row fit-indices data.frame, restricted to the headline columns, including WAIC for the MCMC engine), and, for the MCMC engine only, mcmc_info (chain configuration and convergence diagnostics).

See Also

robust_lpa, print.robust_lpa

Examples

data(neuro_data)
x <- scale(as.matrix(neuro_data[, c("Memory", "RT_Stroop")]))
fit <- robust_lpa(x, G = 2, model = 2, n_starts = 3, max_iter = 20)
summary(fit)