Comprehensive Guide to K and r Model Selection

Introduction

This vignette provides a comprehensive guide to selecting the optimal number of clusters (K) and factors (r) in mixture-based muscle synergy analysis. While individual selection methods are documented elsewhere, this guide focuses on:

Two-stage K,r selection: A stable approach for joint optimization
Post-hoc evaluation: Tools for comparing multiple model candidates
Integrated workflow: End-to-end decision-making process for analysts

Load Package and Example Data

library(synergyMixR)

# Generate synthetic data
sim_data <- generate_synergy_data(
  N = 60, K = 3, r = 4, M = 8,
  T_each = 100,
  cluster_sep_spatial = 1.0,
  seed = 42
)

list_of_data <- sim_data$list_of_data
z_true <- sim_data$true_cluster
cat("True: K =", 3, ", r =", 4, "\n")

Part 1: Two-Stage K,r Selection

Why Two-Stage?

In mixture models, K (number of clusters) and r (number of factors) have strong interactions. Simultaneously searching over all (K, r) combinations can be:

Computationally expensive: O(|Kvec| × |rvec|) model fits
Unstable: Different initializations may prefer different (K, r) pairs
Confounded: Increasing K can compensate for inadequate r, and vice versa

The two-stage approach addresses these issues by:

First fixing r and optimizing K
Then fixing K at the optimal value and optimizing r

Algorithm

Stage 1: K Selection (r fixed)
├── Fix r = r_fixed (default: 2)
├── For each K in Kvec:
│   └── Fit model, compute BIC
└── Select best_K = argmin(BIC)

Stage 2: r Selection (K fixed)
├── Fix K = best_K
├── For each r in rvec:
│   └── Fit model, compute BIC
└── Select best_r = argmin(BIC)

Usage Example

# Two-stage selection for MFA
result_2stage <- select_optimal_K_r_mfa_twostage(
  list_of_data,
  Kvec = 2:5,      # Candidates for K
  rvec = 1:6,      # Candidates for r
  r_fixed = 2,     # Fixed r for Stage 1
  max_iter = 50,
  n_init = 3,
  use_kmeans_init = TRUE
)

# Results
cat("=== Two-Stage Selection Results ===\n")
cat("Stage 1 (K search with r =", 2, "):\n")
print(result_2stage$stage1_summary)

cat("\nStage 2 (r search with K =", result_2stage$best_K, "):\n")
print(result_2stage$stage2_summary)

cat("\nFinal selection: K =", result_2stage$best_K, ", r =", result_2stage$best_r, "\n")

For MPCA Models

# Two-stage selection for MPCA
result_2stage_mpca <- select_optimal_K_r_mpca_twostage(
  list_of_data,
  Kvec = 2:5,
  rvec = 1:6,
  r_fixed = 2,
  method = "EM"
)

cat("MPCA: K =", result_2stage_mpca$best_K, ", r =", result_2stage_mpca$best_r, "\n")

Choosing r_fixed

The r_fixed parameter in Stage 1 affects the K selection. Guidelines:

r_fixed	When to Use
2 (default)	Safe default for most applications
3-4	When you expect moderate synergy complexity
1	When data is very noisy or sample size is small

If you have prior knowledge about r (e.g., from literature or global r estimation), use that value:

# Use global r estimate as r_fixed
global_result <- select_r_global_mfa(list_of_data, rvec = 2:6)
r_estimated <- global_result$r_global

result_informed <- select_optimal_K_r_mfa_twostage(
  list_of_data,
  Kvec = 2:5,
  rvec = 1:6,
  r_fixed = r_estimated  # Use estimated r
)

Comparison: Two-Stage vs Grid Search

Method	Complexity	Stability	Optimality
Grid Search	O(	Kvec	×
Two-Stage	O(	Kvec	+

Recommendation: - Use two-stage for initial exploration or when computation is limited - Use grid search for final model selection when resources allow

Part 2: Post-Hoc Evaluation

After running a grid search, you may want to compare multiple (K, r) candidates beyond just BIC. The post-hoc evaluation functions provide additional metrics.

posthoc_mfa_evaluation()

This function computes additional metrics for all models from a grid search:

# First, run grid search
grid_result <- select_optimal_K_r_mfa(
  list_of_data,
  Kvec = 2:4,
  rvec = 2:5,
  n_init = 3
)

# Post-hoc evaluation
df_posthoc <- posthoc_mfa_evaluation(list_of_data, grid_result)

# View results (sorted by BIC)
print(df_posthoc)

Output Columns

Column	Description	Good Value
`K`	Number of clusters	-
`r`	Number of factors	-
`logLik`	Log-likelihood	Higher is better
`BIC`	Bayesian Information Criterion	Lower is better
`GlobalVAF`	Variance Accounted For	> 0.85
`minClusterSize`	Smallest cluster size	> N/(2K)
`maxClusterSize`	Largest cluster size	< N - N/(2K)

Interpreting Cluster Sizes

Extreme cluster size imbalance can indicate:

Overfitting: One cluster captures outliers
Wrong K: True K may be different
Poor initialization: Try more n_init

# Check for balanced clusters
df_balanced <- df_posthoc[
  df_posthoc$minClusterSize >= 10 &  # At least 10 subjects per cluster
  df_posthoc$GlobalVAF >= 0.85,      # Adequate reconstruction
]

cat("Candidates with balanced clusters and good VAF:\n")
print(df_balanced)

get_model_by_K_r()

Retrieve a specific model for detailed analysis:

# Get a model that's not BIC-optimal but has other desirable properties
# e.g., higher VAF, more balanced clusters
alt_model <- get_model_by_K_r(grid_result, K_target = 3, r_target = 4)

if (!is.null(alt_model)) {
  # Analyze this model
  alt_model <- compute_factor_scores_mfa(alt_model, list_of_data)

  # Visualize synergy patterns
  plot_cluster_synergy_loadings_mfa(alt_model, plot_type = "heatmap")

  # Check factor score smoothness
  smoothness <- compute_smoothness_penalty(alt_model$factor_scores)
  cat("Smoothness penalty:", smoothness, "\n")
}

Comparing Multiple Candidates

# Compare top 3 candidates by different criteria
top_by_bic <- df_posthoc[order(df_posthoc$BIC), ][1:3, ]
top_by_vaf <- df_posthoc[order(-df_posthoc$GlobalVAF), ][1:3, ]

cat("Top 3 by BIC:\n")
print(top_by_bic[, c("K", "r", "BIC", "GlobalVAF")])

cat("\nTop 3 by VAF:\n")
print(top_by_vaf[, c("K", "r", "BIC", "GlobalVAF")])

# If different, investigate both

Part 3: Integrated Workflow

This section presents the recommended end-to-end workflow for analysts.

Workflow Diagram

┌─────────────────────────────────────────────────────────────────┐
│ Phase 1: Preliminary Investigation                              │
├─────────────────────────────────────────────────────────────────┤
│ 1. Global r estimation (K-independent)                         │
│    select_r_global_mfa(list_of_data, rvec = 2:6)               │
│    → Obtain r_global as baseline                               │
│                                                                 │
│ 2. Visualize CV curve                                          │
│    plot(result, type = "global_cv")                            │
└────────────────────┬────────────────────────────────────────────┘
                     │
                     ▼
┌─────────────────────────────────────────────────────────────────┐
│ Phase 2: Joint K,r Selection                                    │
├─────────────────────────────────────────────────────────────────┤
│ Option A: Grid Search (comprehensive)                           │
│   select_optimal_K_r_mfa(Kvec = 1:5, rvec = r_global ± 1)      │
│                                                                 │
│ Option B: Two-Stage (efficient)                                 │
│   select_optimal_K_r_mfa_twostage(Kvec = 2:5, rvec = 1:6,      │
│                                    r_fixed = r_global)         │
│                                                                 │
│ Option C: Hybrid r-Selection (for fixed K)                      │
│   select_r_mixture(K = 3, rvec = 2:6, mode = "hybrid")         │
└────────────────────┬────────────────────────────────────────────┘
                     │
                     ▼
┌─────────────────────────────────────────────────────────────────┐
│ Phase 3: Diagnostics and Visualization                          │
├─────────────────────────────────────────────────────────────────┤
│ 1. Summary statistics                                           │
│    print(result); summary(result)                               │
│                                                                 │
│ 2. Visualizations                                               │
│    plot(result, type = "overview")                              │
│    plot_cluster_synergy_loadings_mfa(best_model)               │
│    plot_all_factor_scores_mfa(best_model, list_of_data)        │
│                                                                 │
│ 3. Check for issues                                             │
│    - Boundary r values                                          │
│    - Unbalanced clusters                                        │
│    - Non-smooth factor scores                                   │
└────────────────────┬────────────────────────────────────────────┘
                     │
                     ▼
┌─────────────────────────────────────────────────────────────────┐
│ Phase 4: Post-Hoc Evaluation and Final Decision                 │
├─────────────────────────────────────────────────────────────────┤
│ 1. Compare multiple candidates                                   │
│    df <- posthoc_mfa_evaluation(list_of_data, grid_result)     │
│                                                                 │
│ 2. Retrieve alternative models                                   │
│    alt <- get_model_by_K_r(grid_result, K = 3, r = 4)          │
│                                                                 │
│ 3. Final decision based on:                                     │
│    - BIC vs VAF trade-off                                       │
│    - Cluster balance                                            │
│    - Physiological interpretability                             │
│    - Factor score smoothness                                    │
└─────────────────────────────────────────────────────────────────┘

Complete Example

# ========================================
# Phase 1: Preliminary Investigation
# ========================================

cat("=== Phase 1: Global r Estimation ===\n")
global_result <- select_r_global_mfa(
  list_of_data,
  rvec = 2:6,
  criterion = "CV"
)
r_global <- global_result$r_global
cat("Estimated r_global:", r_global, "\n")

# ========================================
# Phase 2: Joint K,r Selection
# ========================================

cat("\n=== Phase 2: Grid Search ===\n")
# Focus search around r_global
rvec_focused <- max(1, r_global - 1):min(8, r_global + 2)

grid_result <- select_optimal_K_r_mfa(
  list_of_data,
  Kvec = 1:5,
  rvec = rvec_focused,
  n_init = 5,
  use_kmeans_init = TRUE
)

best_K <- grid_result$best_model_info$K
best_r <- grid_result$best_model_info$r
cat("BIC-optimal: K =", best_K, ", r =", best_r, "\n")

# ========================================
# Phase 3: Diagnostics
# ========================================

cat("\n=== Phase 3: Diagnostics ===\n")
best_model <- grid_result$best_model_info$model

# Compute factor scores
best_model <- compute_factor_scores_mfa(best_model, list_of_data)

# Check cluster sizes
cluster_sizes <- table(best_model$z)
cat("Cluster sizes:", paste(cluster_sizes, collapse = ", "), "\n")

# Check VAF
vaf <- compute_global_vaf_mfa(list_of_data, best_model)
cat("Global VAF:", round(vaf, 4), "\n")

# Check smoothness
smoothness <- compute_smoothness_penalty(best_model$factor_scores)
cat("Smoothness penalty:", round(smoothness, 2), "\n")

# Visualize
plot_cluster_synergy_loadings_mfa(best_model, plot_type = "heatmap")

# ========================================
# Phase 4: Post-Hoc Evaluation
# ========================================

cat("\n=== Phase 4: Post-Hoc Comparison ===\n")
df_posthoc <- posthoc_mfa_evaluation(list_of_data, grid_result)
print(head(df_posthoc, 10))

# Final decision checklist
cat("\n=== Final Decision Checklist ===\n")
cat("[ ] VAF >= 0.85:", vaf >= 0.85, "\n")
cat("[ ] BIC is optimal or near-optimal\n")
cat("[ ] Clusters are reasonably balanced:",
    min(cluster_sizes) >= length(list_of_data) / (2 * best_K), "\n")
cat("[ ] r is not at boundary:",
    best_r > min(rvec_focused) & best_r < max(rvec_focused), "\n")
cat("[ ] Synergy loadings are interpretable: (visual inspection)\n")
cat("[ ] Factor scores are smooth: (visual inspection)\n")

Decision Checklist

Use this checklist for final model selection:

Check	How to Verify	Good Status
VAF sufficient	`compute_global_vaf_mfa()`	>= 0.85-0.90
BIC near optimal	`grid_result$summary`	Within 10 of minimum
Clusters balanced	`table(model$z)`	min >= N/(2K)
r not at boundary	Check `best_r` vs `rvec`	Not at edges
Loadings interpretable	`plot_cluster_synergy_loadings_mfa()`	Clear structure
Scores smooth	`plot_all_factor_scores_mfa()`	No high-frequency noise

When to Reconsider

Scenario 1: r at boundary

if (best_r == max(rvec_focused)) {
  cat("Warning: r at upper boundary. Expand rvec.\n")
  # Re-run with expanded range
}

Scenario 2: Unbalanced clusters

if (min(cluster_sizes) < 5) {
  cat("Warning: Very small cluster. Consider:\n")
  cat("  - Reducing K\n")
  cat("  - Checking for outliers\n")
  cat("  - Using more initializations\n")
}

Scenario 3: Low VAF

if (vaf < 0.80) {
  cat("Warning: Low VAF. Consider:\n")
  cat("  - Increasing r\n")
  cat("  - Checking data quality\n")
  cat("  - Reducing noise in preprocessing\n")
}

Summary

This guide covered three key aspects of model selection:

Two-Stage Selection: Efficient, stable approach for joint K,r optimization
- Use select_optimal_K_r_mfa_twostage()
- Set r_fixed based on global r estimate or domain knowledge
Post-Hoc Evaluation: Compare candidates beyond BIC
- Use posthoc_mfa_evaluation() for comprehensive metrics
- Use get_model_by_K_r() to retrieve specific models
Integrated Workflow: Four-phase decision process
- Phase 1: Preliminary investigation (global r)
- Phase 2: Joint K,r selection (grid or two-stage)
- Phase 3: Diagnostics and visualization
- Phase 4: Post-hoc comparison and final decision

Key Recommendations:

Start with global r estimation to inform search range
Use focused grid search around estimated r
Always check cluster balance and VAF
Visualize synergy patterns for interpretability
Document your decision rationale

Session Info

sessionInfo()