Metrics Reference: Selection Criteria for K and r

Yusuke Matsui

2026-08-14

Introduction

This reference document provides detailed explanations of all metrics used for selecting the optimal number of factors (r) and clusters (K) in mixture-based muscle synergy analysis. Each metric has specific strengths and limitations, making it important to understand when and how to use each one.

Quick Reference Table

Metric	Purpose	Direction	Best For	Limitations
CV	Generalization	Lower	Global r estimation	Computationally expensive
VAF	Reconstruction	Higher	Minimum r threshold	Cannot detect overfitting
BIC	Model complexity	Lower	Joint K,r selection	Sensitive to sample size
Smoothness	Temporal quality	Lower	EMG-specific validation	Domain-specific
VSS	Interpretability	Higher	Post-rotation evaluation	Requires rotation
PA	Dimension estimation	Count	Exploratory analysis	Conservative

Detailed Metric Descriptions

---

Cross-Validation (CV)

Theoretical Background

Cross-validation estimates the model’s ability to generalize to unseen data. In the context of factor analysis, we use leave-one-subject-out CV:

For each subject $i$: 1. Fit model on $N-1$ subjects (excluding $i$) 2. Compute error on held-out subject $i$ 3. Average errors across all subjects

Mathematical Formulation

Negative Log-Likelihood CV: $$\text{CV}_{\text{NegLogLik}} = \frac{1}{N} \sum_{i=1}^{N} -\log p(X_i | \hat{\theta}_{-i})$$

Frobenius Norm CV: $$\text{CV}_{\text{Frob}} = \frac{1}{N} \sum_{i=1}^{N} \|X_i - \hat{X}_{i,-i}\|_F^2$$

where $\hat{\theta}_{-i}$ are parameters estimated without subject $i$.

Implementation

library(synergyMixR)

# CV-based r selection
result_cv <- select_r(
  list_of_data,
  method = "cv",
  rvec = 2:6,
  model_type = "MFA"
)

# View results
print(result_cv)
cat("Selected r:", result_cv$r_selected, "\n")

# Visualize
plot(result_cv$rvec, result_cv$cv_errors, type = "b",
     xlab = "r", ylab = "CV Error",
     main = "Cross-Validation Curve")
abline(v = result_cv$r_selected, col = "red", lty = 2)

When to Use

• Primary use case: Global r estimation before considering cluster structure
• Data requirements: Sufficient subjects (N >= 15-20 recommended)
• Computational note: Most expensive metric; consider reducing n_init for speed

Interpretation Guidelines

CV Curve Pattern	Interpretation	Action
Clear minimum	Strong signal for optimal r	Use minimum
Flat curve	r has little effect	Consider smallest r
Monotonic decrease	May need larger r range	Expand rvec
Noisy curve	Unstable estimation	Increase n_init or subjects

---

Variance Accounted For (VAF)

Theoretical Background

VAF measures the proportion of total variance in the data explained by the model. It is analogous to $R^2$ in regression.

Mathematical Formulation

$$\text{VAF} = 1 - \frac{\text{SSE}}{\text{SST}} = 1 - \frac{\sum_{i=1}^{N} \|X_i - \hat{X}_i\|_F^2}{\sum_{i=1}^{N} \|X_i - \bar{X}\|_F^2}$$

where: - $X_i$ is the observed data for subject $i$ - $\hat{X}_i$ is the reconstructed data - $\bar{X}$ is the grand mean across all subjects

Implementation

# Compute VAF for a fitted model
vaf <- compute_global_vaf_mfa(list_of_data, mfa_fit)
cat("Global VAF:", round(vaf, 4), "\n")

# VAF-based r selection
result_vaf <- select_r_by_vaf_mfa(
  list_of_data,
  K = 3,
  rvec = 2:6,
  vaf_threshold = 0.90  # Select smallest r achieving 90% VAF
)

cat("Selected r:", result_vaf$best_r, "\n")
print(result_vaf$summary)

VAF Thresholds in Literature

Domain	Typical Threshold	Reference
Muscle synergy (locomotion)	0.90	Ivanenko et al., 2004
Muscle synergy (reaching)	0.85-0.90	d’Avella et al., 2006
General EMG analysis	0.80-0.85	Tresch et al., 2006

Caution: VAF Alone is Insufficient

VAF increases monotonically with r. Choosing r based only on VAF will favor overfitting:

# Demonstrate VAF's monotonic increase
for (r in 1:6) {
  fit <- mfa_em_fit(list_of_data, K = 3, r = r, n_init = 1)
  vaf_r <- compute_global_vaf_mfa(list_of_data, fit)
  cat("r =", r, ": VAF =", round(vaf_r, 4), "\n")
}
# Output: VAF always increases with r

Solution: Use VAF as a threshold constraint, not the primary selection criterion.

---

Bayesian Information Criterion (BIC)

Theoretical Background

BIC balances model fit (likelihood) against model complexity (number of parameters). It penalizes complex models to prevent overfitting.

Mathematical Formulation

$$\text{BIC} = -2 \log L(\hat{\theta}) + p \log(n)$$

where: - $L(\hat{\theta})$ is the maximized likelihood - $p$ is the number of free parameters - $n$ is the sample size (total time points)

For MFA: $$p = K \cdot (M \cdot r + 2M) + (K - 1)$$

• $M \cdot r$: loadings per cluster
• $2M$: means and diagonal noise per cluster
• $K - 1$: mixing proportions

Implementation

# Compute BIC for a single model
loglik <- mfa_fit$logLik
bic <- compute_BIC_mfa(loglik, K = 3, r = 4, M = 8, N_total_rows = 6000)

# Grid search with BIC selection
grid_result <- select_optimal_K_r_mfa(
  list_of_data,
  Kvec = 1:5,
  rvec = 1:5
)

# View BIC surface
print(grid_result$summary)

# Visualize BIC heatmap
library(ggplot2)
ggplot(grid_result$summary, aes(x = r, y = K, fill = BIC)) +
  geom_tile() +
  scale_fill_viridis_c(direction = -1) +
  geom_point(data = grid_result$best_model_info,
             aes(x = r, y = K), color = "red", size = 4)

BIC vs AIC

Criterion	Penalty	Tendency
BIC	$p \log(n)$	More parsimonious (prefers simpler models)
AIC	$2p$	Less conservative

For muscle synergy analysis with moderate sample sizes, BIC is generally preferred.

Sample Size Sensitivity

BIC’s penalty term $p \log(n)$ grows with sample size. With very large datasets: - BIC may select very simple models - Consider AIC or cross-validation as alternatives

---

Smoothness Penalty

Theoretical Background

The smoothness penalty is designed for time-series data like EMG. It measures how “wiggly” the factor scores are over time. When r is too large, models tend to fit high-frequency noise, producing non-smooth factors.

Mathematical Formulation

Using second-order differences: $$\text{Smoothness} = \sum_{i=1}^{N} \sum_{j=1}^{r} \sum_{t=3}^{T_i} (z_{i,j,t} - 2z_{i,j,t-1} + z_{i,j,t-2})^2$$

This is equivalent to summing the squared discrete second derivatives of all factor score time series.

Implementation

# Compute factor scores first
mfa_fit <- compute_factor_scores_mfa(mfa_fit, list_of_data)

# Compute smoothness
smoothness <- compute_smoothness_penalty(mfa_fit$factor_scores)
cat("Smoothness penalty:", smoothness, "\n")

# Compare smoothness across r values
for (r in 2:6) {
  fit <- mfa_em_fit(list_of_data, K = 3, r = r, n_init = 1)
  fit <- compute_factor_scores_mfa(fit, list_of_data)
  smooth_r <- compute_smoothness_penalty(fit$factor_scores)
  cat("r =", r, ": Smoothness =", round(smooth_r, 2), "\n")
}

Interpretation

Smoothness Level	Interpretation
Low	Smooth, physiologically plausible activations
High	Noisy, potentially overfitting high-frequency components

Typical Pattern

As r increases: 1. Initially: smoothness may stay low (capturing true structure) 2. Eventually: smoothness increases sharply (fitting noise)

Look for the “elbow” where smoothness begins to increase rapidly.

Muscle Synergy Context

In muscle synergy analysis, factor scores represent muscle activation patterns that should: - Vary smoothly over the movement cycle - Not contain high-frequency oscillations (unless physiologically meaningful) - Reflect coordinated motor control

---

Very Simple Structure (VSS)

Theoretical Background

VSS evaluates the interpretability of factor loadings. The concept of “simple structure” (Thurstone, 1947) suggests that each variable (muscle) should load primarily on one factor (synergy).

Mathematical Formulation

Threshold Method: $$\text{VSS}_{\text{threshold}} = \frac{1}{K \cdot M} \sum_{k=1}^{K} \sum_{m=1}^{M} \mathbb{1}\left[\sum_{j=1}^{r} \mathbb{1}(|\lambda_{k,m,j}| > \tau) = 1\right]$$

Counts the proportion of muscles with exactly one loading above threshold $\tau$.

Ratio Method: $$\text{VSS}_{\text{ratio}} = \frac{1}{K \cdot M} \sum_{k=1}^{K} \sum_{m=1}^{M} \frac{\max_j(\lambda_{k,m,j}^2)}{\sum_j \lambda_{k,m,j}^2}$$

Measures how concentrated each muscle’s loading is on its dominant factor.

Implementation

# Compute VSS
vss_thresh <- compute_vss_criterion(mfa_fit$Lambda, threshold = 0.5, method = "threshold")
vss_ratio <- compute_vss_criterion(mfa_fit$Lambda, method = "ratio")

cat("VSS (threshold):", round(vss_thresh, 4), "\n")
cat("VSS (ratio):", round(vss_ratio, 4), "\n")

# VSS after rotation (recommended)
rotated_fit <- rotate_mfa_model(mfa_fit, rotation = "varimax")
vss_rotated <- compute_vss_criterion(rotated_fit$Lambda, method = "ratio")
cat("VSS after varimax:", round(vss_rotated, 4), "\n")

Rotation Importance

Unrotated factor loadings often do not exhibit simple structure. Apply varimax (or another rotation) before computing VSS:

# Without rotation - often lower VSS
vss_unrot <- compute_vss_criterion(mfa_fit$Lambda, rotate = "none")

# With rotation - typically higher VSS
vss_rot <- compute_vss_criterion(mfa_fit$Lambda, rotate = "varimax")

cat("Unrotated VSS:", round(vss_unrot, 4), "\n")
cat("Rotated VSS:", round(vss_rot, 4), "\n")

Interpretation Guidelines

VSS Score	Interpretation
> 0.8	Excellent simple structure
0.6 - 0.8	Good simple structure
0.4 - 0.6	Moderate; some complex loadings
< 0.4	Poor simple structure; difficult to interpret

---

Parallel Analysis (PA)

Theoretical Background

Parallel analysis compares observed eigenvalues to those expected from random data. Factors are retained only if their eigenvalues exceed the random baseline.

Algorithm

1. Compute eigenvalues of pooled covariance matrix: $\lambda_1 \geq \lambda_2 \geq ... \geq \lambda_M$
2. Generate B random permutations of the data
3. For each permutation, compute eigenvalues
4. Calculate 95th percentile of random eigenvalues: $\lambda^{95}_1, ..., \lambda^{95}_M$
5. Select r = max{j : $\lambda_j > \lambda^{95}_j$}

Implementation

# Global MFA-based r estimation
result_pa <- select_r_global_mfa(
  list_of_data,
  rvec = 1:8
)

cat("Selected r:", result_pa$r_global, "\n")

# View results
print(result_pa)

Characteristics

Strengths: - Objective, data-driven threshold - No arbitrary VAF cutoff needed - Good for exploratory analysis

Limitations: - Conservative (may underestimate r) - Ignores cluster structure - Assumes single population

---

Metric Comparison and Selection Strategy

Single Metric Usage

Metric	When to Use Alone
CV	Global r estimation with sufficient data
BIC	Joint K,r selection via grid search
PA	Quick exploratory analysis

Recommended Combinations

Conservative Selection (Avoid Overfitting)

# Prioritize BIC and smoothness
result <- select_r_multicriteria_mfa(
  list_of_data, K = 3, rvec = 2:6,
  weights = c(BIC = 2, VAF = 1, smoothness = 2, VSS = 1)
)

Interpretability Focus

# Prioritize VSS
result <- select_r_multicriteria_mfa(
  list_of_data, K = 3, rvec = 2:6,
  weights = c(BIC = 1, VAF = 1, smoothness = 1, VSS = 2)
)

Balanced Approach (Recommended Default)

result <- select_r_multicriteria_mfa(
  list_of_data, K = 3, rvec = 2:6,
  vaf_threshold = 0.85,
  weights = c(BIC = 1, VAF = 1, smoothness = 1, VSS = 1)
)

Decision Framework

┌─────────────────────────────────────────────────────────────┐
│ Step 1: Establish baseline with global CV                   │
│ r_cv = select_r(method = "cv", ...)$r_selected               │
└────────────────────┬────────────────────────────────────────┘
                     │
                     ▼
┌─────────────────────────────────────────────────────────────┐
│ Step 2: Check VAF threshold                                  │
│ If VAF(r_cv) < 0.85 → consider larger r                     │
│ If VAF(r_cv) > 0.95 → may be overfitting, check smoothness  │
└────────────────────┬────────────────────────────────────────┘
                     │
                     ▼
┌─────────────────────────────────────────────────────────────┐
│ Step 3: Validate with secondary metrics                      │
│ - Check smoothness (should not spike)                       │
│ - Check VSS after rotation (should be reasonable)           │
└────────────────────┬────────────────────────────────────────┘
                     │
                     ▼
┌─────────────────────────────────────────────────────────────┐
│ Step 4: Joint K,r selection with BIC                        │
│ Use focused rvec around validated r                         │
└─────────────────────────────────────────────────────────────┘

Summary

Key Takeaways

1. No single metric is perfect - each has strengths and limitations
2. CV is the gold standard for generalization but is computationally expensive
3. VAF should be a constraint, not the primary criterion
4. BIC works well for joint K,r selection with moderate sample sizes
5. Smoothness is EMG-specific - leverage temporal structure of your data
6. VSS requires rotation - always apply varimax before evaluation
7. PA is conservative - useful for exploratory analysis

Recommended Workflow

1. Start with global CV or PA for r estimation
2. Set VAF >= 0.85 as minimum threshold
3. Use multi-criteria selection with balanced weights
4. Validate final selection with smoothness and VSS
5. Visualize loadings and scores for interpretability

References

• Ivanenko, Y. P., et al. (2004). Five basic muscle activation patterns account for muscle activity during human locomotion. Journal of Physiology.
• d’Avella, A., et al. (2006). Control of fast-reaching movements by muscle synergy combinations. Journal of Neuroscience.
• Tresch, M. C., et al. (2006). Matrix factorization algorithms for the identification of muscle synergies. Journal of Neurophysiology.
• Revelle, W., & Rocklin, T. (1979). Very simple structure: An alternative procedure for estimating the optimal number of interpretable factors. Multivariate Behavioral Research.
• Horn, J. L. (1965). A rationale and test for the number of factors in factor analysis. Psychometrika.

Session Info

sessionInfo()