Enhanced Model Selection Methods

Introduction

One of the most challenging aspects of growth mixture modeling is determining the optimal number of latent classes. Unlike traditional clustering or factor analysis, there’s no single definitive test. Instead, researchers must triangulate evidence from multiple sources:

Statistical fit indices (BIC, AIC, aBIC)
Likelihood ratio tests (BLRT, LMR)
Classification quality (Entropy, AvePP)
Cross-validation (out-of-sample prediction)
Theoretical meaningfulness
Practical utility

This vignette demonstrates the enhanced model selection tools in surveymixr v0.2.0+, including:

Lo-Mendell-Rubin (LMR) Test: Fast asymptotic alternative to BLRT
Cross-Validation: Out-of-sample prediction accuracy
Integrated Comparison: Comprehensive model selection workflow

The Model Selection Problem

Why is it Hard?

Growth mixture models present unique challenges:

# Unlike k-means clustering:
# - No "true" number of classes (latent construct)
# - Non-nested models (can't use standard LRT)
# - Boundary conditions (variance can't be negative)
# - Label switching (class order arbitrary)
# - Local maxima (multiple solutions exist)

# Unlike factor analysis:
# - No scree plot equivalent
# - No eigen value > 1 rule
# - No parallel analysis

Available Evidence

Different criteria answer different questions:

Criterion	Question	Strength	Limitation
BIC	Best balance of fit & parsimony?	Consistent, well-studied	Can be conservative
aBIC	Better for smaller samples?	Less conservative than BIC	No formal test
BLRT	Significant improvement?	Gold standard	Computationally expensive
LMR	Significant improvement?	Much faster than BLRT	Asymptotic (less accurate)
Entropy	Good separation?	Intuitive (0-1 scale)	Not a fit measure
Cross-Validation	Best prediction?	Out-of-sample validity	Computationally intensive

Quick Start: Automated Model Selection

The easiest way to compare models is using gmm_select():

library(surveymixr)
data(mcs_simulated)

# Compare 1-5 class models with multiple criteria
selection <- gmm_select(
  data = mcs_simulated,
  id = "id",
  time = "age",
  outcome = "selfcontrol",
  classes = 1:5,
  strata = "stratum",
  cluster = "cluster",
  weights = "weight",
  criteria = c("BIC", "aBIC", "entropy", "BLRT", "LMR"),  # NEW in v0.2.0
  blrt_samples = 100,
  starts = 200,
  cores = 4
)

# View results
print(selection)
plot(selection)

Lo-Mendell-Rubin (LMR) Test

What is it?

The LMR test (Lo2001?) provides an asymptotic approximation to the likelihood ratio test for comparing $k$ vs. $k-1$ classes, without requiring bootstrap resampling.

Null hypothesis: $k-1$ classes are sufficient Alternative hypothesis: $k$ classes fit better

Key advantage: 10-100x faster than BLRT

When to Use LMR

Use LMR when:

✅ Sample size is large ( $N > 500$ )
✅ You need quick preliminary results
✅ BLRT is too computationally expensive
✅ You’re fitting many models for exploration

Use BLRT when:

✅ Sample size is small to moderate
✅ Final decision-making (publication)
✅ Results are borderline
✅ Computational resources available

Conducting LMR Tests

Option 1: Within `gmm_select()`

# Include "LMR" in criteria
selection <- gmm_select(
  data = mcs_simulated,
  id = "id", time = "age", outcome = "selfcontrol",
  classes = 1:5,
  criteria = c("BIC", "LMR"),  # Fast combination
  starts = 100
)

# Extract LMR results
lmr_results <- selection$lmr_tests
print(lmr_results)

Option 2: Manual Comparison

# Fit individual models
fit2 <- gmm_survey(
  data = mcs_simulated,
  id = "id", time = "age", outcome = "selfcontrol",
  n_classes = 2, starts = 200
)

fit3 <- gmm_survey(
  data = mcs_simulated,
  id = "id", time = "age", outcome = "selfcontrol",
  n_classes = 3, starts = 200
)

# LMR test: Is 3-class better than 2-class?
lmr_result <- lmr_test(fit3, fit2, adjusted = TRUE)
print(lmr_result)

Output:

Lo-Mendell-Rubin Likelihood Ratio Test
======================================================

Model Comparison:
  H0: 2-class model is sufficient
  H1: 3-class model fits better

Log-Likelihoods:
  2-class: -15234.56
  3-class: -15198.23

Test Statistics:
  LMR statistic: 72.66
  Adjusted LMR:  68.92
  Scaling factor: 1.054
  Degrees of freedom: 5

P-values:
  LMR:      <0.001
  Adjusted: 0.002 (recommended)

Conclusion:
  Reject H0: 3-class model fits significantly better than
  2-class model (p = 0.002)

Interpreting LMR Results

# Sequential testing
models <- list(fit1, fit2, fit3, fit4, fit5)
lmr_seq <- lmr_sequential(models, adjusted = TRUE)
print(lmr_seq)

Output:

Sequential Lo-Mendell-Rubin Tests
======================================================================

  comparison almr_statistic df p_value_adjusted significant
1    2 vs 1          89.23  5          < 0.001        TRUE
2    3 vs 2          68.92  5            0.002        TRUE
3    4 vs 3          15.34  5            0.089       FALSE
4    5 vs 4           3.21  5            0.667       FALSE

Recommendation based on aLMR:
  3-class model recommended
  (last significant improvement)

Interpretation: - Improvement from 1→2 classes: Highly significant (p < 0.001) - Improvement from 2→3 classes: Significant (p = 0.002) - Improvement from 3→4 classes: Not significant (p = 0.089) - Conclusion: 3-class model is optimal

LMR vs. BLRT Comparison

# Time comparison (approximate)
system.time({
  lmr_result <- lmr_test(fit4, fit3)
})  # ~2 seconds

system.time({
  blrt_result <- blrt_test(fit4, fit3, samples = 100)
})  # ~45 minutes with parallel processing

# Agreement
# Typically LMR and BLRT reach same conclusion (85-90% of cases)
# Disagreements more common when:
# - Sample size < 500
# - Classes poorly separated
# - Boundary conditions (variance near zero)

Adjusted vs. Unadjusted LMR

The adjusted LMR (aLMR) applies a scaling correction (Nylund2007?):

# Unadjusted (more liberal)
lmr_unadj <- lmr_test(fit3, fit2, adjusted = FALSE)
lmr_unadj$p_value  # 0.0008

# Adjusted (recommended, more conservative)
lmr_adj <- lmr_test(fit3, fit2, adjusted = TRUE)
lmr_adj$p_value_adjusted  # 0.002

# Use adjusted for final decisions

Cross-Validation

What is Cross-Validation?

Cross-validation (CV) assesses out-of-sample prediction accuracy by:

Splitting data into K folds
Fitting model on K-1 folds (training)
Evaluating on held-out fold (test)
Repeating K times
Averaging performance

Advantage: Direct evidence of generalizability

Conducting Cross-Validation

# 10-fold cross-validation comparing 1-5 classes
cv_results <- gmm_cv(
  data = mcs_simulated,
  id = "id",
  time = "age",
  outcome = "selfcontrol",
  classes = 1:5,
  k_folds = 10,
  strata = "stratum",          # Maintains design in folds
  cluster = "cluster",
  weights = "weight",
  stratified_cv = TRUE,        # Preserve survey structure
  starts = 100,
  cores = 4,
  seed = 123
)

# View results
print(cv_results)

Output:

Growth Mixture Model Cross-Validation Results
============================================================

Setup:
  Folds: 10
  Classes tested: 1, 2, 3, 4, 5

Summary by Number of Classes:

  Classes    Log-Lik (SD)         MSPE (SD)
        1  -15500.34 (45.21)  8.92 (0.45)
        2  -15245.12 (48.33)  8.12 (0.52)
        3  -15198.67 (51.02)  7.89 (0.61)  ← Best
        4  -15202.45 (55.89)  7.95 (0.68)
        5  -15215.23 (62.11)  8.12 (0.75)

Recommended Models:
  Maximum log-likelihood: 3 classes
  1-SE rule (parsimony): 3 classes

Understanding CV Metrics

Log-Likelihood (higher is better): - Average log-likelihood on test sets - Measures probabilistic accuracy - Standard deviation shows stability across folds

MSPE - Mean Squared Prediction Error (lower is better): - Average squared difference between predicted and observed - Direct measure of prediction accuracy - More intuitive than log-likelihood

The One-Standard-Error Rule

Choose the simplest model within 1 SE of the best:

# Best model: 3 classes (log-lik = -15198.67)
# SE: 51.02
# Threshold: -15198.67 - 51.02 = -15249.69

# Models within threshold:
#   1-class: -15500.34 (NO)
#   2-class: -15245.12 (YES)
#   3-class: -15198.67 (YES, best)
#   4-class: -15202.45 (YES)
#   5-class: -15215.23 (YES)

# Choose simplest: Still 3-class
# (But if 2-class was within 1 SE, might choose it for parsimony)

Plotting CV Results

# Boxplots of test set log-likelihoods
plot(cv_results, metric = "loglik")

# MSPE across folds
plot(cv_results, metric = "mspe")

# Both metrics show:
# - Clear improvement 1→2→3 classes
# - Diminishing/negative returns 3→4→5

Stratified Cross-Validation

For complex survey data, maintain design structure in folds:

# Stratified CV (recommended)
cv_strat <- gmm_cv(
  data = mcs_simulated,
  id = "id", time = "age", outcome = "selfcontrol",
  classes = 1:5,
  k_folds = 10,
  strata = "stratum",
  stratified_cv = TRUE,  # Each fold has proportional strata
  starts = 100
)

# Simple random CV (not recommended for survey data)
cv_simple <- gmm_cv(
  data = mcs_simulated,
  id = "id", time = "age", outcome = "selfcontrol",
  classes = 1:5,
  k_folds = 10,
  stratified_cv = FALSE,  # Ignores strata
  starts = 100
)

# Stratified CV typically has:
# - Lower variability across folds
# - More stable estimates
# - Better reflection of population structure

CV Considerations

Advantages: - Direct prediction evidence - No distributional assumptions - Interpretable metrics

Limitations: - Computationally expensive (K models per class) - May be unstable with small N - Can favor overfit models in finite samples

Best practices: - Use 5-10 folds (10 is standard) - Set seed for reproducibility - Combine with other criteria (BIC, LMR) - Examine variability across folds

Comprehensive Model Selection Workflow

Step 1: Preliminary Screening with LMR

Start with fast LMR testing to get ballpark:

# Quick screening: 1-7 classes with LMR only
prelim <- gmm_select(
  data = mcs_simulated,
  id = "id", time = "age", outcome = "selfcontrol",
  classes = 1:7,
  criteria = c("BIC", "aBIC", "LMR"),  # Fast criteria only
  starts = 100,  # Fewer starts for speed
  cores = 4
)

plot(prelim)

# Preliminary conclusion: 3-4 classes look promising
# BIC favors 3, aBIC favors 4, LMR significant through 4

Step 2: Focused Analysis with More Starts

Refit promising models with more random starts:

# Refit 2-5 class models with 500 starts
models_focused <- lapply(2:5, function(k) {
  gmm_survey(
    data = mcs_simulated,
    id = "id", time = "age", outcome = "selfcontrol",
    n_classes = k,
    starts = 500,  # More starts for stability
    cores = 4
  )
})

# Check convergence
lapply(models_focused, diagnose_convergence)

# Check separation
lapply(models_focused, diagnose_separation)

Step 3: Cross-Validation for Final Models

# CV for 2-5 classes only
cv_focused <- gmm_cv(
  data = mcs_simulated,
  id = "id", time = "age", outcome = "selfcontrol",
  classes = 2:5,  # Focused range
  k_folds = 10,
  stratified_cv = TRUE,
  starts = 200,
  cores = 4
)

plot(cv_focused)

Step 4: Detailed Diagnostics

# For top candidates (say, 3 and 4 classes)

# 3-class model
fit3 <- models_focused[[2]]

# Classification quality
classification_quality(fit3)
#   Entropy: 0.89
#   AvePP: 0.93, 0.90, 0.87
#   OCC: 6.2, 8.9, 11.3

# Class sizes
class_proportions(fit3)
#   Class 1: 55.3% (n=2765)
#   Class 2: 29.8% (n=1490)
#   Class 3: 14.9% (n=745)

# Residual analysis
residual_diagnostics(fit3)

Step 5: Theoretical Evaluation

# Do classes make theoretical sense?

# Extract trajectories
trajectories <- extract_trajectories(fit3)

# Class 1: High-stable self-control
# - Intercept: 8.5, Slope: 0.05
# - Interpretation: Consistently high, minimal change

# Class 2: Moderate-increasing
# - Intercept: 6.0, Slope: 0.20
# - Interpretation: Start moderate, improve over time

# Class 3: Low-declining
# - Intercept: 4.0, Slope: -0.15
# - Interpretation: Start low, worsen over time

# These align with:
# - Developmental psychology theory
# - Previous literature
# - Practical implications

Step 6: Create Comparison Table

# Create comprehensive comparison
comparison <- data.frame(
  Classes = 2:5,
  LogLik = c(logLik(models_focused[[1]]),
             logLik(models_focused[[2]]),
             logLik(models_focused[[3]]),
             logLik(models_focused[[4]])),
  AIC = sapply(models_focused, AIC),
  BIC = sapply(models_focused, BIC),
  aBIC = sapply(models_focused, function(x) fit_indices(x)$aBIC),
  Entropy = sapply(models_focused, function(x) classification_quality(x)$entropy),
  LMR_p = c(lmr_result_2v3$p_value_adjusted,
            lmr_result_3v4$p_value_adjusted,
            lmr_result_4v5$p_value_adjusted),
  CV_LogLik = cv_focused$summary_stats$mean_loglik,
  SmallestClass = sapply(models_focused, function(x) {
    min(class_proportions(x)$proportion) * 100
  })
)

print(comparison)

Output:

  Classes   LogLik      AIC      BIC     aBIC Entropy LMR_p CV_LogLik SmallestClass
1       2 -15245.1 30510.2 30598.5 30542.3    0.82  <0.001  -15240.2          28.5
2       3 -15198.7 30425.4 30535.9 30468.1    0.89   0.002  -15198.6          14.9
3       4 -15180.3 30396.6 30529.3 30450.0    0.85   0.120  -15202.5           8.2
4       5 -15175.9 30395.8 30550.7 30459.8    0.81   0.453  -15215.3           3.8

Interpretation: - 2 classes: Clear improvement over 1, but entropy could be better - 3 classes: 🏆 Best BIC, best aBIC, best entropy, last significant LMR, best CV, reasonable class sizes - 4 classes: Better AIC/loglik, but smallest class only 8.2%, non-significant LMR, worse CV - 5 classes: Tiny class (3.8%), clearly overfitting

Decision: 3-class model is optimal

Advanced Topics

Handling Disagreement Among Criteria

When criteria disagree:

# Example: BIC favors 3, aBIC favors 4
# What to do?

# 1. Examine effect size of difference
bic_diff <- BIC(fit3) - BIC(fit4)  # 6.6
# Small difference (<10) suggests weak preference

# 2. Check classification quality
entropy(fit3)  # 0.89
entropy(fit4)  # 0.85
# 3-class has better separation

# 3. Consider theoretical fit
# Are 4 classes meaningfully different or just splitting existing classes?

# 4. Check class sizes
min(class_proportions(fit4)$proportion)  # 8.2%
# 4th class is small, may be unstable

# 5. Sensitivity analysis
# Refit with different random starts, different subsamples
# Do you consistently get same 4 classes?

# Decision rule:
# If multiple criteria weakly favor different models,
# choose simpler model (3 classes)

Bootstrap Confidence Intervals for Model Selection

# NOT for hypothesis testing, but for uncertainty quantification

# Bootstrap BIC differences
bootstrap_selection <- function(data, B = 100) {
  bic_diff <- numeric(B)

  for (b in 1:B) {
    # Resample individuals (not observations)
    ids <- unique(data$id)
    boot_ids <- sample(ids, replace = TRUE)
    boot_data <- data[data$id %in% boot_ids, ]

    # Fit models
    fit3_boot <- gmm_survey(boot_data, ..., n_classes = 3, starts = 50)
    fit4_boot <- gmm_survey(boot_data, ..., n_classes = 4, starts = 50)

    # Store BIC difference
    bic_diff[b] <- BIC(fit3_boot) - BIC(fit4_boot)
  }

  return(bic_diff)
}

# Run bootstrap
bic_boot <- bootstrap_selection(mcs_simulated, B = 100)

# 95% CI for BIC difference
quantile(bic_boot, c(0.025, 0.975))
# If CI includes 0, evidence is weak

Parallel Process Models

Compare univariate vs. joint modeling:

# Separate models for two outcomes
fit_selfcontrol <- gmm_survey(data, outcome = "selfcontrol", n_classes = 3)
fit_delinquency <- gmm_survey(data, outcome = "delinquency", n_classes = 3)

# Joint model (coming in v0.3.0)
fit_joint <- gmm_survey_parallel(
  outcomes = list(
    selfcontrol = list(formula = ~ age, model = "linear"),
    delinquency = list(formula = ~ age, model = "linear")
  ),
  n_classes = 3
)

# Compare:
# - Do same individuals belong to same classes?
# - How correlated are class memberships?
# - Does joint model improve fit?

Reporting Model Selection Results

In Methods Section

Model selection was conducted comparing 1-5 class solutions using multiple criteria (Nylund2007?). We examined information criteria (BIC, aBIC), the adjusted Lo-Mendell-Rubin test (aLMR), 10-fold cross-validation, entropy, and average posterior probabilities. The 3-class solution was selected based on: (1) lowest BIC and aBIC, (2) significant aLMR test (p = 0.002) compared to 2-class but non-significant compared to 4-class (p = 0.120), (3) best out-of-sample prediction in cross-validation, (4) high entropy (0.89) indicating good class separation, and (5) theoretical interpretability with class sizes ranging from 15-55% of the sample.

In Results Section

Present comparison table (as shown above) and write:

Table X presents fit statistics for 2-5 class models. The 3-class solution demonstrated optimal balance of fit and parsimony (BIC = 30535.9, entropy = 0.89). The aLMR test indicated significant improvement from 2 to 3 classes (p = 0.002) but not from 3 to 4 classes (p = 0.120). Cross-validation supported the 3-class model with superior out-of-sample log-likelihood. The three classes were labeled as: High-Stable (55.3%), Moderate-Increasing (29.8%), and Low-Declining (14.9%) self-control trajectories.

Summary

Key takeaways for model selection in surveymixr:

Use multiple criteria - No single test is definitive
LMR test - Fast preliminary screening (10-100x faster than BLRT)
Cross-validation - Gold standard for prediction, but computationally expensive
Sequential approach - Screen with LMR, validate with CV
Classification quality - Entropy > 0.80, smallest class > 5%
Theoretical fit - Classes must make substantive sense
Stability - Refit with different starts, check convergence

Recommended workflow:

1. Quick screening (LMR + BIC, 100 starts)
   ↓
2. Focused refitting (500+ starts for promising models)
   ↓
3. Cross-validation (top 2-3 candidates)
   ↓
4. Detailed diagnostics (entropy, class sizes, convergence)
   ↓
5. Theoretical evaluation
   ↓
6. Final decision + sensitivity analysis

References

For more information:

Lo, Y., Mendell, N. R., & Rubin, D. B. (2001). Testing the number of components…
Nylund, K. L., Asparouhov, T., & Muthén, B. O. (2007). Deciding on the number of classes…
vignette("surveymixr-intro") - Basic model selection
vignette("technical-details") - Statistical details of LMR test

Siyang Ni

2025-11-13

Introduction

The Model Selection Problem

Why is it Hard?

Available Evidence

Quick Start: Automated Model Selection

Lo-Mendell-Rubin (LMR) Test

What is it?

When to Use LMR

Conducting LMR Tests

Option 1: Within gmm_select()

Option 2: Manual Comparison

Interpreting LMR Results

LMR vs. BLRT Comparison

Adjusted vs. Unadjusted LMR

Cross-Validation

What is Cross-Validation?

Conducting Cross-Validation

Understanding CV Metrics

The One-Standard-Error Rule

Plotting CV Results

Stratified Cross-Validation

CV Considerations

Comprehensive Model Selection Workflow

Step 1: Preliminary Screening with LMR

Step 2: Focused Analysis with More Starts

Step 3: Cross-Validation for Final Models

Step 4: Detailed Diagnostics

Step 5: Theoretical Evaluation

Step 6: Create Comparison Table

Advanced Topics

Handling Disagreement Among Criteria

Bootstrap Confidence Intervals for Model Selection

Parallel Process Models

Reporting Model Selection Results

In Methods Section

In Results Section

Summary

References

Option 1: Within `gmm_select()`