The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
This release adds robust growth mixture models for longitudinal data
and is also a correctness release: every estimator in the package was
re-checked against a known reference (closed-form maximum likelihood,
mclust, lme4/nlme, or simulation
with known truth), and the problems found are fixed below. Several of
them change numerical results relative to 1.0.0; re-running analyses
fitted with 1.0.0 is recommended.
robust_gmm() fits growth mixture
models (finite mixtures of linear mixed-effects models; Verbeke &
Lesaffre, 1996; Muthen & Shedden, 1999) and, with
random = "none", latent class growth analysis (Nagin,
1999), for one or several outcomes measured repeatedly on unbalanced
occasions (long-format data; persons may have different visit times and
missing outcomes). Each class has a polynomial mean trajectory per
outcome (degree); within classes, persons have correlated
random intercepts or intercepts and slopes (random,
re_structure), with covariance and residual variances
shared or class-specific (re_cov, resid_var).
Every person contributes the exact likelihood of the values actually
observed, so estimation is valid under missing-at-random drop-out.
robust_method = "t"; the t linear mixed model of Pinheiro,
Liu & Wu, 2001, with nu estimated), which down-weight
whole outlying persons and keep a proper likelihood, or Huber
weights.lme4/nlme exactly
(checked in the tests).plot_mcmc_chains()).lambda_growth shrinks the growth terms of every class
toward zero (identifying classes that are stable on an outcome);
lambda_diff shrinks the class trajectories toward each
other, element-wise or by outcome (group_diff = TRUE, which
selects the outcomes that differentiate the classes; Xie, Pan &
Shen, 2008; the groups are standardized by the information of their
coefficients, Simon & Tibshirani, 2012). Penalties are adaptive by
default (Zou, 2006; Wang & Leng, 2008), so that
lambda = z^2 / N sets to zero terms whose Wald statistic is
below about z; the penalized step is solved exactly by
ADMM, the BIC counts only the free coefficients,
relax = TRUE refits the selected model without penalty
(relaxed Lasso), and the MCMC engine uses the corresponding Bayesian
(group, adaptive) Lasso priors.estimate_gmm_robust() fits a grid of
numbers of classes (and random-effect specifications), optionally tuning
the LASSO penalties by BIC or by K-fold cross-validation over
persons.blrt_gmm_robust(): bootstrapped
likelihood ratio test for the number of classes, simulating from the
fitted null model on the observed visit schedule and missingness
pattern.plot_robust_gmm() plots the class
trajectories over the observed individual trajectories;
print() and summary() methods.bch_robust() accepts
robust_gmm() fits, to relate the trajectory classes to
baseline characteristics or distal outcomes (one value per person, in
the order of $ids); its bootstrap correction resamples
whole persons.neuro_long (400
simulated persons, up to six annual visits, three outcomes, drop-out and
gross errors; its generating script is installed in
inst/scripts/) and a second vignette,
vignette("robust-growth-mixture", package = "RobustLPA").Parameter count, AIC, BIC and SABIC of
robust_lpa(). 1.0.0 left out of the count every
mean with absolute value below 1e-5, even without a LASSO penalty. On
standardized data with one profile all means are 0, so p
parameters were missing (model 1 on the five neuro_data
variables: 5 instead of 10; model 6: 15 instead of 20) and the BIC of
G = 1 was about 27.6 too low, favouring one profile. Means
are now left out only when the EM LASSO (lambda > 0) has
set them exactly to zero.
Parallel workers on Windows. With
cores > 1 the worker processes loaded RobustLPA from the
default library, which could hold another version than the one used in
the session (e.g. a new version installed in a separate library),
silently mixing code from the two. They now load the session’s copy, and
a warning is given if they cannot.
BCH (bch_robust()) now estimates the
classification error matrix correctly. 1.0.0 computed
D[g, c] as the average posterior of profile g
among observations assigned to c, i.e. P(C = g | assigned
c), the reverse of the P(assigned s | C = t) required by Bolck, Croon
& Hagenaars
(2004) / Vermunt (2010), and indexed the weight matrix
W = D^-1 with transposed subscripts. The two errors cancel
only for symmetric classification error. With unequal, overlapping
profiles (80/20, entropy about 0.55; true difference in the distal mean
= 1) the corrected profile-2 mean had a median bias of +0.44 and RMSE of
1.9 in 1.0.0 (worse than the uncorrected comparison), versus -0.05 and
0.19 now. Rows of $Classification_Matrix now sum to
1.
Constrained models (1-5) could return an unconstrained
solution. The EM initialization did not apply the
variance-covariance constraints, and the “keep the best iterate” logic
could then return those unconstrained starting values. For example,
model = 1 on neuro_data returned non-diagonal,
profile-specific covariances with a log-likelihood (-1395.4)
higher than the model-1 maximum (-1469.5, as found by
mclust EEI), which also invalidated BIC comparisons across
models. Classical fits now reproduce the mclust
maximum-likelihood solutions for models 1, 2, 3 and 6.
Missing data are now handled by exact maximum likelihood. The E-step already used the observed-data (FIML) density, but the M-step used pairwise available-case moments, which are biased under MAR. The M-step now uses the exact EM for incomplete data (conditional expectations of the missing entries plus their conditional covariance; Ghahramani & Jordan, 1994). Example: with 30% MAR missingness the estimated mean of the incomplete variable was -0.37 (truth 0) in 1.0.0 and now matches the closed-form ML solution. The MCMC engine now uses the equivalent data augmentation step.
Numerical underflow for outlying observations. Densities are now computed on the log scale throughout (log-sum-exp), so extreme observations receive valid posterior probabilities. In 1.0.0 an extreme outlier could get a posterior row of all zeros and a random profile.
MCMC models 1 and 5 now match the EM definitions. 1.0.0’s sampler used a single variance shared by all variables (a spherical model) for models 1 and 5; they now use variable-specific variances, as documented and as in the EM engine.
Models 4 and 5 are now fitted properly by both
engines. The EM engine used an approximate M-step that made the
likelihood decrease and stopped at poor solutions; it now uses
conditional-maximization steps that never decrease the likelihood. On
neuro_data (classical fit, G = 2), model 4 improves from
log-likelihood -1358.7 (83% classification accuracy) to -1315.4 (94%),
and model 5 from -1387.6 (52%) to -1318.8 (89%). The MCMC engine now
samples the scales exactly (slice sampling) and the correlation matrices
by Metropolis steps, instead of an approximate conditional that drifted
away from the posterior.
Huber weights now accumulate across EM iterations. In 1.0.0 each M-step recomputed the weights against a fresh non-robust estimate, so outliers could mask themselves; the weights are now computed against the current robust estimates (an iteratively reweighted M-estimator). Because the Huber fixed point is not a maximizer of the Gaussian likelihood, the final iterate is returned and no log-likelihood-decrease warnings are issued for Huber fits.
bch_robust(correction = "bootstrap") refits now run
with cores = 1 inside each replicate (previously they
inherited the original fit’s cores, causing nested
parallelism).
Profile labels (MCMC draws, bootstrap refits) are now matched by an exact optimal assignment (Hungarian algorithm) on standardized means, instead of a greedy nearest-mean rule. MCMC draws are relabeled to a fixed pivot (the preliminary EM solution).
estimate_profiles_robust() warns instead of silently
ignoring tune_lasso = TRUE with the MCMC engine, validates
its inputs, reports why a model failed, and uses the same parallel
backend as the rest of the package.
robust_lpa(robust_method = "t"), with nu
estimated by default or fixed): a likelihood-based robust model (Peel
& McLachlan, 2000) fitted by ECME in the EM engine and by an exact
Gibbs sampler (latent scales; nu updated with the scales
integrated out) in the MCMC engine. AIC/BIC, the BLRT and BCH are all
used as intended with this model, and it is markedly more resistant than
Huber weighting to gross outliers (in a test with 5% gross outliers on
data with unit variances, the Huber variance estimates were 2.3-3.6 –
3.5-6.1 in 1.0.0 – while the t fit estimated nu = 2.4 and a core scale
of 0.7-0.8).$weights: every fit now reports each observation’s
robustness weight (small values flag outliers).$fit$WAIC, $waic), starts its chains from a
preliminary EM fit with dispersed perturbations, and adapts its
Metropolis step sizes during burn-in.blrt_robust() simulates the null data from the fitted
model family (multivariate t for robust_method = "t").set.seed() now gives
identical results for any value of cores (one seed is drawn
per work unit before dispatch).init = "kmeans"; init = "random" restores
random soft partitions).$converged and
$iterations.testthat test suite checks every estimator against
closed-form or reference solutions.Major release: a new Bayesian MCMC estimation engine alongside the existing EM engine, an optional classical (non-robust) estimation mode, LASSO regularization, a new BCH auxiliary-variable method, parallel computing throughout, S3 print/summary methods, a new bundled dataset, and a new package vignette.
robust_lpa(engine = "MCMC") estimates the same six
variance-covariance parameterizations via Gibbs sampling, with a
Bayesian Lasso (Laplace-prior) option for the profile means
(prior_laplace) and the same Huber-weighted
robust/classical toggle as the EM engine
(robust/alpha). Runs multiple independent
chains by default (n_chains = 4), aligns each chain’s
arbitrarily-ordered profile labels to a common ordering before pooling,
and reports classic Gelman-Rubin and effective sample size convergence
diagnostics ($mcmc_diagnostics, via the new
coda dependency) for every scalar parameter, warning if any
exceeds the standard 1.1 R-hat threshold.
plot_mcmc_chains() draws multi-chain trace plots via
bayesplot.robust_lpa()
gains a robust argument (default TRUE) and
alpha, so either engine can be run with Huber-weighted
robust estimation (as in 0.1.0) or with classical (non-robust) maximum
likelihood, using the same engine, missing-data handling, and variance-
covariance parameterizations either way.lambda), and
estimate_profiles_robust(tune_lasso = TRUE) selects
lambda by k-fold cross-validation (k_folds,
lambda_grid).cores argument
runs independent work in parallel throughout the package: EM random
restarts or MCMC chains within a single robust_lpa() call;
the model/profile grid (and cross-validation folds) in
estimate_profiles_robust(); bootstrap replicates in
blrt_robust(); and bootstrap correction replicates in
bch_robust(). Uses parallel::mclapply() on
macOS/Linux and a parallel::makeCluster() PSOCK backend on
Windows.bch_robust() implements the Bolck-Croon-Hagenaars (2004)
three-step method for relating fitted profiles to a continuous
auxiliary/distal outcome variable, correcting for classification error.
Optionally (correction = "bootstrap") adds a nonparametric
bootstrap correction for classification uncertainty in the step-1 model,
reporting bootstrap standard errors, confidence intervals, and a Wald
chi-square test for the profile means
($Bootstrap_Correction) – a practical approximation to the
Bakk, Oberski & Vermunt (2014) sandwich correction, preferable to
the base $ANOVA_Table F-test (which treats the
classification matrix as fixed/known) for publication-grade
inference.blrt_robust() (present since 0.1.0) now supports both
estimation engines (engine = "EM" or "MCMC"),
replicates the observed data’s FIML missingness pattern in every
simulated bootstrap sample, and runs its bootstrap replicates in
parallel via cores.robust_lpa() now returns an object of class
"robust_lpa" with print() and
summary() methods: print() is a compact
overview (engine, model, profile count, N, headline fit indices, mixing
proportions), and summary() adds per-profile means, profile
sizes, and, for the MCMC engine, the Gelman-Rubin/effective sample size
range – both focused on what’s needed for a first read, rather than
dumping the full fitted object.robust_lpa() fits now store $data (the numeric
matrix used for estimation) and $call_args (every argument
controlling the fit), so other functions can refit the identical
specification on new or resampled data; used internally by
bch_robust()’s bootstrap correction.model, more n_starts, or fewer
profiles).neuro_data (n =
250, two groups) is a new synthetic dataset with a genuine group-level
difference in covariance structure (not only in means) and calibrated,
variable-magnitude outlier contamination, so that model/profile
selection, robust-vs-classical comparisons, and classification recovery
all behave as a realistic worked example throughout the package’s
documentation. See ?neuro_data for details and
data-raw/generate_neuro_data.R for the full generative
code.vignette("RobustLPA")) walks through the full workflow:
choosing a variance-covariance model, EM vs. MCMC estimation, robust
vs. classical estimation, LASSO regularization, model/profile selection,
the bootstrapped likelihood ratio test, and BCH auxiliary-variable
analysis, on neuro_data.parallel dependency under
Suggests; reference DOIs throughout the documentation use
the \doi{} Rd macro, as recommended by CRAN.RcppArmadillo) to natively
handle missing data.blrt_robust() for Bootstrapped Likelihood
Ratio Tests.plot_robust_lpa() for publication-ready
visualizations using ggplot2.These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.