Package 'synergyMixR'

Description

Collects and combines results from individual PBS array job tasks.

Usage

aggregate_partial_results(
  output_dir,
  pattern = "result_.*\\.rds$",
  remove_partial = FALSE,
  verbose = TRUE
)

Arguments

Value

Data frame with combined results, or NULL if no results found.

Examples

## Not run: 
# After PBS jobs complete
results <- aggregate_partial_results("results_sweep_modelsel")

# View results
head(results)

## End(Not run)

Description

Finds the optimal mapping from estimated cluster labels to true cluster labels using the Hungarian algorithm. This resolves the label switching problem inherent in clustering.

Usage

align_clusters_to_truth(z_true, z_hat)

Arguments

Details

The Hungarian algorithm (via clue::solve_LSAP) finds the assignment that maximizes the overlap between estimated and true clusters. This is essential for fair evaluation of clustering methods where label identities are arbitrary.

When K_hat != K_true, the algorithm handles the mismatch by:

• If K_hat < K_true: Some true clusters may not have a match

• If K_hat > K_true: Some estimated clusters map to the same true cluster

Value

A list with:

mapping

Named integer vector where names are estimated labels and values are corresponding true labels.

z_aligned

Integer vector of aligned estimated labels (same length as z_hat, with labels remapped to match z_true).

accuracy

Numeric scalar: proportion of correctly matched assignments after alignment.

Returns a list with all NA values if inputs are invalid.

Examples

## Not run: 
z_true <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
z_hat <- c(3, 3, 3, 1, 1, 1, 2, 2, 2)  # Labels are permuted
result <- align_clusters_to_truth(z_true, z_hat)
result$accuracy  # Should be 1.0 (perfect after alignment)

## End(Not run)

Description

After running rotate_mfa_model or any rotation, each factor's sign can be arbitrary. This function enforces a rule such as: "The variable with the largest absolute loading in each factor must be positive." Then, if that factor's column is flipped, we also flip the corresponding column in the factor scores, to keep the model consistent.

Usage

align_factor_signs_mfa(mfa_fit, method = c("max_abs"), cluster_ids = NULL)

Arguments

Details

- For each cluster k, we loop over each factor (column). We find the row index i_max that has the largest absolute loading. If Lambda[k][i_max, factor] < 0, we multiply that entire column by -1. - If factor_scores is present, we also multiply the corresponding column in each subject's factor score by -1 to maintain consistency. - This doesn't change the model fit or BIC, just the sign convention.

Value

The same mfa_fit object, but with possibly flipped signs in Lambda columns and factor_scores columns.

Description

Similar to align_factor_signs_mfa, but for a Mixture PCA model whose loadings are in W[[k]].

Usage

align_factor_signs_mpca(mpca_fit, method = c("max_abs"), cluster_ids = NULL)

Arguments

Details

- For each cluster k, we loop over each factor (column) in W[[k]]. If the entry with largest absolute loading is negative, we multiply that entire column by -1. - If factor_scores is present, we also multiply the corresponding column in each subject's factor scores by -1.

Value

The same mpca_fit object, but with possibly flipped signs in W columns and factor_scores columns (if present).

Description

Performs k-means clustering on subject-level features to obtain initial cluster assignments for mixture model fitting.

Usage

assign_by_kmeans(features, K)

Arguments

Details

This function uses the standard k-means algorithm with 10 random starts (nstart=10) to find a stable clustering solution. The resulting cluster assignments can be used as initial values for EM algorithms in mixture models.

Value

An integer vector of length N containing cluster labels (1..K). Each element indicates which cluster the corresponding subject is assigned to.

Examples

## Not run: 
# Create example features
features <- matrix(rnorm(100), nrow = 10, ncol = 10)

# Assign to 3 clusters
z_init <- assign_by_kmeans(features, K = 3)
table(z_init)  # Check cluster sizes

## End(Not run)

Description

Applies a bandpass filter to EMG signals to extract the frequency band of interest (typically 20-450 Hz for surface EMG).

Usage

bandpass_filter(
  data,
  sampling_rate,
  low_cutoff = 20,
  high_cutoff = 450,
  filter_order = 4,
  filter_type = "butterworth",
  zero_phase = TRUE
)

Arguments

Details

The bandpass filter removes:

• Low-frequency noise and motion artifacts (below low_cutoff)

• High-frequency noise (above high_cutoff)

Standard EMG bandwidth is 20-450 Hz for surface EMG. The lower cutoff helps remove motion artifacts, while the upper cutoff prevents aliasing and removes high-frequency noise.

Value

Filtered EMG signal, same structure as input.

Examples

## Not run: 
# Generate sample EMG-like signal
fs <- 1000  # 1000 Hz sampling rate
t <- seq(0, 1, by = 1/fs)
emg <- rnorm(length(t)) + 0.1 * sin(2 * pi * 10 * t)  # with low-freq noise

# Apply bandpass filter
filtered <- bandpass_filter(emg, sampling_rate = fs)

## End(Not run)

Description

Computes bootstrap confidence intervals for the loading matrices in mixture PCA or mixture FA models. Uses Procrustes alignment to ensure bootstrap samples are comparable across resamples.

Usage

bootstrap_loading_ci(
  list_of_data,
  K,
  r,
  B = 200,
  alpha = 0.05,
  mc_cores = 1,
  model_type = c("MPCA", "MFA"),
  max_iter = 50,
  n_init = 1,
  use_kmeans_init = TRUE,
  seed = NULL,
  verbose = TRUE
)

Arguments

Details

The bootstrap procedure: 1. Fit the model to the full dataset to obtain reference loadings 2. For each bootstrap iteration b = 1, ..., B: a. Resample subjects with replacement b. Fit the model to the bootstrap sample c. Align bootstrap loadings to reference using Procrustes rotation 3. Compute percentile confidence intervals from aligned samples

Procrustes alignment is essential because factor loadings are only identifiable up to rotation and sign. Without alignment, bootstrap samples would not be comparable, leading to inflated variance estimates.

Value

A list with the following elements:

loadings_median

List of K matrices (M x r), median loadings for each cluster

loadings_lower

List of K matrices (M x r), lower CI bound for each cluster

loadings_upper

List of K matrices (M x r), upper CI bound for each cluster

loadings_se

List of K matrices (M x r), standard errors for each cluster

bootstrap_samples

Array of dimension (B x K x M x r) containing all aligned bootstrap samples (optional, if return_samples=TRUE)

reference_fit

The reference model fit used for alignment

n_successful

Number of successful bootstrap iterations

alpha

The significance level used

Procrustes Alignment

For each cluster k, the bootstrap loading matrix W_k(b) is aligned to the reference W_k using orthogonal Procrustes: W_aligned(b) = W_k(b) R* where R* minimizes ||W_k - W_k(b) R||_F subject to R'R = I.

The optimal R is found via SVD: if W_k' W_k(b) = UDV', then R* = VU'.

See Also

mixture_pca_em_fit, mfa_em_fit

Examples

## Not run: 
# Simulate some data
sim <- simulate_dynamic_synergy_data(N = 50, K = 2, r = 3, M = 8, T_each = 100)

# Compute bootstrap CIs
boot_ci <- bootstrap_loading_ci(
  list_of_data = sim$X_list,
  K = 2,
  r = 3,
  B = 100,
  alpha = 0.05,
  mc_cores = 4
)

# View median loadings for cluster 1
print(boot_ci$loadings_median[[1]])

# View 95% CI width for first synergy, first muscle, cluster 1
width <- boot_ci$loadings_upper[[1]][1,1] - boot_ci$loadings_lower[[1]][1,1]
print(width)

## End(Not run)

Description

Uses the final hard assignments mfa_fit$z, reconstruct each subject's data.

Usage

calc_reconstruction_error_mixtureFA(list_of_data, mfa_fit)

Arguments

Value

SSE (sum of squared errors).

Description

Compute SSE for a MixturePCA model

Usage

calc_reconstruction_error_mixturePCA(list_of_data, mpca_fit)

Arguments

Value

SSE (sum of squared errors).

Description

Compute Total Reconstruction SSE for Single FA

Usage

calc_reconstruction_error_singleFA(list_of_data, fa_res)

Arguments

Value

A numeric value of the sum of squared errors (SSE).

Description

Compute Total Reconstruction SSE for Single PCA

Usage

calc_reconstruction_error_singlePCA(list_of_data, pca_res)

Arguments

Value

SSE (sum of squared errors).

Description

This function concatenates all matrices in list_of_data into one big matrix, calculates the global mean, and then sums the squared differences from that mean. The result can be used to compute VAF (Variance Accounted For).

Usage

calc_total_SST(list_of_data)

Arguments

Value

A numeric scalar, the total sum of squares (SST).

Description

Convenience function that extracts responsibilities from a fitted model and computes calibration against true labels.

Usage

calibrate_model(fit, z_true, n_bins = 10)

Arguments

Value

A calibration_result object

Examples

## Not run: 
# Fit model
fit <- mixture_pca_em_fit(list_of_data, K = 3, r = 2)

# Compute calibration
calib <- calibrate_model(fit, z_true)
print(calib)

## End(Not run)

Description

Checks whether the future ecosystem packages are installed.

Usage

check_future_available(batchtools = FALSE)

Arguments

Value

Logical. TRUE if packages are available.

Description

Checks the status of a PBS job.

Usage

check_pbs_job_status(job_id)

Arguments

Value

List with job status information.

Examples

## Not run: 
status <- check_pbs_job_status("12345")

## End(Not run)

Description

This function performs diagnostic checks to investigate SSE (Sum of Squared Errors) behavior in mixture model fits, particularly when SSE increases monotonically with latent dimension r. It examines three key aspects: EM convergence quality, noise variance (Psi) distribution, and VAF-SSE correlation.

Usage

check_SSE_diagnostics(
  fit_results,
  list_of_data,
  mixtureFA_fit = NULL,
  mixturePCA_fit = NULL,
  output_dir = NULL,
  psi_small_thresh = 0.01,
  psi_large_thresh = 100
)

Arguments

Details

The function performs three main diagnostic checks:

1. EM Convergence Quality Check

Since the C++ EM implementation doesn't return a convergence trace, this check examines the final model quality by comparing log-likelihoods and checking for degenerate solutions (e.g., empty clusters, extreme parameter values).

2. Psi (Noise Variance) Distribution Check

Examines the distribution of noise variance parameters in both MFA and MPCA models. Extreme values (very small < 0.01 or very large > 100) may indicate:

• Overfitting (Psi approaches 0)

• Poor model fit (Psi approaches infinity)

• Numerical instability

Creates histograms showing the distribution of Psi values across all clusters and muscles.

3. VAF-SSE Correlation Verification

Since VAF = 1 - SSE/SST, there should be a strong negative correlation (approximately -1) between VAF and SSE across methods. Deviations from this expected relationship may indicate:

• Scaling inconsistencies in SSE calculation

• Different SST values being used across methods

• Numerical precision issues

Creates a scatter plot of VAF vs SSE with the expected theoretical relationship.

Value

A list with diagnostic results:

convergence_ok

Logical; TRUE if convergence appears stable

psi_ok

Logical; TRUE if Psi values are in reasonable range

correlation_ok

Logical; TRUE if VAF-SSE correlation is strongly negative

warnings

Character vector of warning messages

psi_summary_mfa

Summary statistics for MFA Psi values

psi_summary_mpca

Summary statistics for MPCA Psi values

vaf_sse_correlation

Correlation coefficient between VAF and SSE

Examples

## Not run: 
# Example 1: Basic usage with fit_all_methods_on_data()
# The fitted models are stored as attributes and automatically extracted
sim_data <- simulate_mixture_data(N = 30, K = 3, r = 3, M = 8, T_each = 200, seed = 123)
fit_results <- fit_all_methods_on_data(
  list_of_data = sim_data$list_of_data,
  z_true = sim_data$z_true,
  N = 30, K = 3, r = 3, M = 8,
  seed = 123
)

# Run diagnostics (fitted models extracted from attributes automatically)
diagnostics <- check_SSE_diagnostics(
  fit_results = fit_results,
  list_of_data = sim_data$list_of_data,
  output_dir = "diagnostics"
)

# Check results
print(diagnostics$warnings)
print(diagnostics$psi_summary_mfa)
print(diagnostics$vaf_sse_correlation)

# Example 2: Pass fitted models explicitly
# Useful when you have fitted models from other sources
mfa_fit <- mfa_em_fit(sim_data$list_of_data, K = 3, r = 3)
mpca_fit <- mixture_pca_em_fit(sim_data$list_of_data, K = 3, r = 3)

# Create a minimal results data frame
results_df <- data.frame(
  Method = c("MixtureFA", "MixturePCA"),
  SSE = c(100, 120),
  VAF = c(0.85, 0.80)
)

# Run diagnostics with explicit fitted models
diagnostics2 <- check_SSE_diagnostics(
  fit_results = results_df,
  list_of_data = sim_data$list_of_data,
  mixtureFA_fit = mfa_fit,
  mixturePCA_fit = mpca_fit,
  output_dir = "diagnostics2",
  psi_small_thresh = 0.005,  # Custom threshold
  psi_large_thresh = 200     # Custom threshold
)

# Example 3: VAF-SSE correlation only (no fitted models)
# Useful for quick checks without EM convergence/Psi diagnostics
diagnostics3 <- check_SSE_diagnostics(
  fit_results = fit_results,
  list_of_data = sim_data$list_of_data
)
# Note: EM convergence and Psi diagnostics will be skipped with a warning

## End(Not run)

Description

This function enables pure method-level benchmarking by generating data using only the v2 dynamic synergy simulator and applying all six analysis methods (SingleFA, SinglePCA, TwoStep_FA, TwoStep_PCA, MixtureFA, MixturePCA). This complements compare_static_vs_dynamic by focusing on method performance rather than model comparison.

Usage

compare_baselines(
  N = 50,
  K = 2,
  r = 2,
  M = 6,
  T_each = 100,
  cluster_sep_spatial = 1,
  cluster_sep_temporal = 1,
  cluster_sep_stability = 1,
  seed = 123,
  n_init = 1,
  use_kmeans_init = FALSE,
  mc_cores = 1,
  method_pca = "EM",
  align_basis = TRUE,
  refine_assignments = TRUE
)

Arguments

Details

This function generates data using simulate_dynamic_synergy_data with three orthogonal control axes (spatial, temporal, stability), then applies all six analysis methods to evaluate their performance. This is useful for:

• Method benchmarking on realistic dynamic synergy data

• Parameter sensitivity analysis for method performance

• Integration with run_param_sweep for systematic evaluation

The function complements compare_static_vs_dynamic by focusing on method performance rather than model comparison. Together, these functions provide comprehensive evaluation of both models and methods.

Value

A data frame with the following columns:

Method

Analysis method name (e.g., "MixtureFA", "SinglePCA")

BIC

Bayesian Information Criterion

ARI

Adjusted Rand Index (clustering accuracy)

SSE

Sum of Squared Errors (reconstruction error)

VAF

Variance Accounted For (1 - SSE/totalSST)

model_type

Always "Dynamic_v2" for this function

seed

Random seed used for data generation

cluster_sep_spatial

Spatial separation parameter

cluster_sep_temporal

Temporal separation parameter

cluster_sep_stability

Stability separation parameter

N

Number of subjects

K

Number of clusters

r

Number of synergies

M

Number of muscles

T_each

Time series length

See Also

compare_static_vs_dynamic for model comparison, run_param_sweep for systematic parameter sweep, simulate_dynamic_synergy_data for data generation

Examples

## Not run: 
# Basic method comparison on v2 data
results <- compare_baselines(
  N = 50, K = 2, r = 2, M = 6, T_each = 100,
  cluster_sep_spatial = 1.0,
  cluster_sep_temporal = 1.0,
  cluster_sep_stability = 1.0,
  seed = 123
)
print(results)

# Compare methods with high spatial separation
results_spatial <- compare_baselines(
  N = 50, K = 2, r = 2, M = 6, T_each = 100,
  cluster_sep_spatial = 2.0,
  cluster_sep_temporal = 0.5,
  cluster_sep_stability = 1.0,
  seed = 456
)

# Visualize method performance
library(ggplot2)
ggplot(results, aes(x = Method, y = ARI)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Method Performance on v2 Dynamic Data")

## End(Not run)

Description

This function uses the v2 dynamic synergy simulator (simulate_dynamic_synergy_data), which provides three orthogonal control axes: spatial structure, temporal coordination, and stability-adaptability. This allows for systematic evaluation of how different motor control strategies affect muscle synergy structures.

Usage

compare_static_vs_dynamic(
  model_types = c("Dynamic_v2"),
  N = 30,
  K = 3,
  r = 3,
  M = 8,
  T_each = 200,
  cluster_sep_spatial = 1,
  cluster_sep_temporal = 1,
  cluster_sep_stability = 1,
  seed = 123,
  output_dir = "results/compare_dynamic_v2",
  n_init = 1,
  use_kmeans_init = FALSE,
  mc_cores = 1
)

Arguments

Details

The function generates data using the Dynamic_v2 model and applies six analysis methods (SingleFA, SinglePCA, TwoStep_FA, TwoStep_PCA, MixtureFA, MixturePCA) to evaluate their performance.

The function saves results to an RDS file in the specified output directory with filename format results_YYYY-MM-DD.rds.

Data is generated using simulate_dynamic_synergy_data, which creates time-varying synergy activation patterns with three orthogonal control axes: (1) Spatial structure, (2) Temporal coordination, and (3) Stability-Adaptability.

The generated data is analyzed using six methods, allowing comparison of method performance on realistic dynamic synergy data.

Value

A data frame with the following columns:

model_type

Character string: "Dynamic_v2"

Method

Analysis method name (e.g., "MixtureFA", "SinglePCA")

BIC

Bayesian Information Criterion

ARI

Adjusted Rand Index (clustering accuracy)

SSE

Sum of Squared Errors (reconstruction error)

VAF

Variance Accounted For (1 - SSE/totalSST)

seed

Random seed used for data generation

cluster_sep_spatial

Spatial separation parameter

cluster_sep_temporal

Temporal separation parameter

cluster_sep_stability

Stability separation parameter

N

Number of subjects

K

Number of clusters

r

Number of synergies

M

Number of muscles

T_each

Time series length

Examples

## Not run: 
# Basic comparison with v2 simulator
results <- compare_static_vs_dynamic(
  N = 30, K = 3, r = 3, M = 8, T_each = 200,
  cluster_sep_spatial = 1.0,
  cluster_sep_temporal = 1.0,
  cluster_sep_stability = 1.0,
  seed = 123
)
print(results)

# Compare with high spatial separation
results_spatial <- compare_static_vs_dynamic(
  N = 30, K = 3, r = 3, M = 8, T_each = 200,
  cluster_sep_spatial = 2.0,
  cluster_sep_temporal = 0.5,
  cluster_sep_stability = 1.0,
  seed = 456
)

# Visualize results
library(ggplot2)
ggplot(results, aes(x = Method, y = ARI)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Method Performance on Dynamic v2 Data")

## End(Not run)

Description

Computes the Adjusted Rand Index between true and predicted cluster assignments. This is a NA-safe wrapper around mclust::adjustedRandIndex.

Usage

compute_ari(true_cluster, pred_cluster)

Arguments

Details

The Adjusted Rand Index (ARI) measures the similarity between two clusterings, adjusted for chance. It is commonly used to evaluate clustering accuracy when ground truth is available.

Value

Numeric scalar: ARI value in [-1, 1], where 1 indicates perfect agreement. Returns NA if inputs are invalid or mclust is not available.

Examples

## Not run: 
true_cluster <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
pred_cluster <- c(1, 1, 2, 2, 2, 2, 3, 3, 3)
ari <- compute_ari(true_cluster, pred_cluster)

## End(Not run)

Description

This function takes a data frame with a file column (typically from read_all_results), reads each result file to extract fitted models, computes VAF for MFA and MPCA, and joins these back to the input data frame.

Usage

compute_ari_vaf_table(df)

Arguments

Value

The input data frame enriched with VAF_MFA and VAF_PCA columns. If models are not available in a file, VAF values will be NA.

Examples

## Not run: 
df <- read_all_results("results")
df <- compute_ari_vaf_table(df)
head(df)

## End(Not run)

Description

Unified interface for computing BIC for both MFA and MPCA models.

Usage

compute_BIC(loglik, K, r, M, N_obs, model_type)

Arguments

Details

BIC = -2 * loglik + p * log(N_obs), where p is the number of parameters. The parameter count differs between MFA and MPCA due to different model structures.

Value

A numeric scalar representing the BIC value (lower is better)

See Also

compute_BIC_mfa, compute_BIC_mpca

Examples

## Not run: 
# Compute BIC for an MFA model
bic <- compute_BIC(loglik, K = 2, r = 3, M = 8, N_obs = 1000, model_type = "MFA")

## End(Not run)

Description

Given the log-likelihood from compute_logLik_mfa, and known K, r, M, N, this function computes BIC = -2 * log(L) + nu * log(N), where nu is the parameter count for the MFA model. The corrected formula is nu = K * (M * r + r * (r-1)/2 + M + M) + (K - 1).

Usage

compute_BIC_mfa(loglik, K, r, M, N_obs)

Arguments

Details

The parameter count includes:

• M*r: factor loadings (Lambda)

• r*(r-1)/2: factor covariances (off-diagonal elements)

• M: unique variances (diagonal Psi)

• M: cluster means (mu)

• K-1: mixing proportions

This corrected formula properly penalizes larger r values, improving r estimation stability.

Note: N_obs should typically be the total number of time points (sum of T_i) rather than the number of subjects, as BIC penalizes based on sample size.

Value

A numeric BIC value.

Description

Given the log-likelihood from compute_logLik_mpca, and the number of clusters K, dimension r, etc., this function computes the corrected BIC: BIC = -2 * log(L) + nu * log(N), where nu is the parameter count for MPCA. The corrected formula is K*(M*r + r*(r-1)/2 + M + 1) + (K-1), i.e. for each cluster we have (M*r) for the directions, r*(r-1)/2 for factor covariances, M for the mean, and 1 for sigma^2. Then (K-1) for mixing proportions.

Usage

compute_BIC_mpca(loglik, K, r, M, N_obs)

Arguments

Details

The parameter count includes:

• M*r: loading matrix (W or P)

• r*(r-1)/2: factor covariances (off-diagonal elements)

• M: cluster means (mu)

• 1: isotropic noise variance (sigma2)

• K-1: mixing proportions

This corrected formula properly penalizes larger r values, improving r estimation stability.

Note: N_obs should typically be the total number of time points (sum of T_i) rather than the number of subjects, as BIC penalizes based on sample size.

Value

A numeric BIC value.

Description

Computes the Bayesian Information Criterion (BIC) for a fitted MFA or MPCA model. This is a unified wrapper that dispatches to the appropriate method-specific function.

Usage

compute_bic_wrapper(method, fit, K, r, M, N_obs)

Arguments

Details

BIC is computed as $-2 \log L + p \log(N)$, where L is the likelihood, p is the number of parameters, and N is the sample size.

Value

Numeric scalar: BIC value. Lower is better. Returns NA if computation fails.

Examples

## Not run: 
fit <- mfa_em_fit(list_of_data, K = 3, r = 2)
N_obs <- sum(sapply(list_of_data, nrow))
M <- ncol(list_of_data[[1]])
bic <- compute_bic_wrapper("MFA", fit, K = 3, r = 2, M = M, N_obs = N_obs)

## End(Not run)

Description

Computes calibration statistics by binning predicted probabilities and comparing to observed frequencies. This helps diagnose whether the model's probabilistic predictions are reliable.

Usage

compute_calibration(resp, z_true, n_bins = 10)

Arguments

Details

Calibration is computed for each cluster separately: 1. Bin the predicted probabilities into n_bins equal-width bins 2. For each bin, compute mean predicted probability and observed frequency 3. A well-calibrated model has mean_predicted close to mean_observed

The Expected Calibration Error (ECE) summarizes overall calibration:

$ECE = \sum_{b=1}^{B} \frac{|S_b|}{N} |acc(S_b) - conf(S_b)|$

where S_b is the set of samples in bin b, acc is accuracy (observed frequency), and conf is confidence (mean predicted probability).

Value

A list with class "calibration_result" containing:

calibration_df

Data frame with columns: bin_midpoint, mean_predicted, mean_observed, n_samples, cluster

ece

Expected Calibration Error (weighted average bin error)

mce

Maximum Calibration Error (worst bin)

n_bins

Number of bins used

K

Number of clusters

N

Number of subjects

See Also

plot_calibration, compute_ece

Examples

## Not run: 
# Fit a mixture model
fit <- mixture_pca_em_fit(list_of_data, K = 3, r = 2)

# Compute calibration using true labels
calib <- compute_calibration(fit$resp, z_true)

# View ECE
print(calib$ece)

# Plot calibration curve
plot_calibration(calib)

## End(Not run)

Description

Given a result with a vector z of length N (the cluster assignments), this function returns a frequency table, i.e., how many subjects are assigned to each cluster.

Usage

compute_cluster_sizes(res)

Arguments

Value

A table of cluster frequencies.

Examples

## Not run: 
tab <- compute_cluster_sizes(fit_mfa)
print(tab)

## End(Not run)

Description

Like compute_cluster_sizes, but for a Mixture PCA model result.

Usage

compute_cluster_sizes_mpca(res)

Arguments

Value

A table of frequencies of cluster assignments.

Description

Computes subspace similarity scores for each cluster, using the aligned mapping between estimated and true clusters. This provides a detailed assessment of how well each cluster's subspace is recovered.

Usage

compute_clusterwise_subspace_score(true_Lambda_list, est_Lambda_list, mapping)

Arguments

Details

For each true cluster k:

1. Find the estimated cluster that maps to k (using the mapping)

2. Compute subspace similarity between true_Lambda_list[[k]] and the corresponding estimated basis

3. If no estimated cluster maps to k, the score is NA

This function is NA-safe: if estimation failed for some clusters (NULL or invalid matrices), those clusters receive NA scores and are excluded from the mean calculation.

Value

A list with:

scores

Named numeric vector of subspace similarity scores for each true cluster (names are cluster indices as characters).

mean_score

Numeric scalar: mean of non-NA scores.

n_matched

Integer: number of true clusters with a valid estimated match.

Returns a list with NA values if inputs are invalid.

Examples

## Not run: 
# Generate data and fit model
data <- generate_synergy_data(N = 30, K = 3, r = 2, seed = 123)
fit <- mfa_em_fit(data$list_of_data, K = 3, r = 2)

# Align clusters
alignment <- align_clusters_to_truth(data$true_cluster, fit$z)

# Compute cluster-wise scores
scores <- compute_clusterwise_subspace_score(
  data$true_Lambda_list,
  fit$Lambda,
  alignment$mapping
)

## End(Not run)

Description

A convenience function to compute only the ECE from responsibilities and true labels.

Usage

compute_ece(resp, z_true, n_bins = 10)

Arguments

Details

ECE measures the average absolute difference between predicted probabilities and observed frequencies. A value of 0 indicates perfect calibration; values above 0.1 suggest significant miscalibration.

Guidelines for interpretation:

• ECE < 0.02: Excellent calibration

• ECE 0.02-0.05: Good calibration

• ECE 0.05-0.10: Moderate calibration

• ECE > 0.10: Poor calibration

Value

Numeric scalar, the ECE value. Lower is better (0 = perfect calibration).

See Also

compute_calibration

Examples

## Not run: 
ece <- compute_ece(fit$resp, z_true)
cat(sprintf("ECE: %.3f\n", ece))

## End(Not run)

Description

Computes the entropy term used in ICL (Integrated Complete-data Likelihood). Higher entropy indicates more uncertain cluster assignments.

Usage

compute_entropy(resp)

Arguments

Details

The entropy is computed as:

$H = -\sum_{i=1}^{N} \sum_{k=1}^{K} r_{ik} \log(r_{ik})$

where $r_{ik}$ is the responsibility of subject i for cluster k. A small constant (1e-16) is added to avoid log(0).

Value

A numeric scalar representing the entropy: -sum(resp * log(resp)).

Examples

## Not run: 
# Responsibilities from MFA fit
entropy <- compute_entropy(fit_mfa$resp)

## End(Not run)

Description

Given a single cluster's parameters (mu, Lambda, Psi) and data X, this function returns the estimated factor scores.

Usage

compute_factor_scores(X, mu, Lambda, Psi)

Arguments

Details

This follows the standard posterior mean formula in factor analysis: $\mathrm{E}[Z|X] = (I + \Lambda^T \Psi^{-1} \Lambda)^{-1} \Lambda^T \Psi^{-1} (X - \mu)$.

Value

A (T x r) matrix of factor scores.

See Also

compute_global_vaf_mfa

Examples

## Not run: 
# Suppose we have X (T x M), mu (M), Lambda (M x r), Psi (M x M).
scores <- compute_factor_scores(X, mu, Lambda, Psi)
head(scores)

## End(Not run)

Description

Given a fitted MFA model (with z, mu, Lambda, Psi for each cluster), this function computes the posterior factor scores Z_i(t) for each subject i and stores them in model$factor_scores[[i]].

Usage

compute_factor_scores_mfa(mfa_fit, list_of_data)

Arguments

Details

Internally, it calls compute_factor_scores_optimized for each subject, using that subject's assigned cluster k = z[i].

Value

The same mfa_fit object, but with a new element factor_scores (a list of length N), where each factor_scores[[i]] is a (T_i x r) matrix of factor scores for subject i.

Examples

## Not run: 
mfa_fit <- compute_factor_scores_mfa(mfa_fit, list_of_data)
# Now mfa_fit$factor_scores[[i]] is available for each subject i

## End(Not run)

Description

Given a fitted MPCA model (with z, mu, W, sigma2 for each cluster), this function computes the posterior factor (PC) scores Z_i(t) for each subject i and stores them in mpca_fit$factor_scores[[i]].

Usage

compute_factor_scores_mpca(mpca_fit, list_of_data)

Arguments

Details

We use compute_factor_scores_mpca_single for each subject, with the appropriate cluster parameters. This assumes isotropic noise sigma2[k] for cluster k.

Value

The same mpca_fit object, but with a new element factor_scores, where factor_scores[[i]] is a (T_i x r) matrix of PC scores for subject i.

Examples

## Not run: 
mpca_fit <- compute_factor_scores_mpca(mpca_fit, list_of_data)
mpca_fit$factor_scores[[1]]  # (T_1 x r) PC scores

## End(Not run)

Description

Given a single subject's data X, and the cluster parameters (mu, W, sigma2), this function computes the posterior factor scores (T x r) under a PPCA assumption with isotropic noise sigma^2 * I.

Usage

compute_factor_scores_mpca_single(X, mu, W, sigma2)

Arguments

Details

We use the posterior mean formula for PPCA with isotropic noise: z(t) = inv(W^T W + sigma2 * I_r) * W^T * (x(t) - mu). This is applied row by row in X, but we can do it in a vectorized manner.

If your noise model is not purely sigma2 * I, you will need a more general approach. Also note that if T=1, the result might be a (1 x r) matrix or (r)-vector, so be sure to keep it as a matrix if desired.

Value

A (T x r) matrix of factor (PC) scores for this subject.

Description

Given a single subject's data X and the cluster's parameters (mu, Lambda, Psi), this function computes the posterior factor scores (of size T x r) in a standard Factor Analysis approach.

Usage

compute_factor_scores_optimized(X, mu, Lambda, Psi)

Arguments

Details

We use the standard formula for Factor Analysis posterior means:

$\mathbf{Z}_i(t) = \left( \mathbf{I}_r + \Lambda^\top \Psi^{-1} \Lambda \right)^{-1} \Lambda^\top \Psi^{-1} \left[\mathbf{x}_i(t) - \mu\right].$

The code is implemented in a vectorized manner. If $r = \mathrm{ncol}(\Lambda)$, then the output dimension is (T x r).

Value

A (T x r) matrix of factor scores Z_i(t), where t = 1..T.

Examples

## Not run: 
X_i <- list_of_data[[i]]  # (T x M)
# from model$mu[[k]], model$Lambda[[k]], model$Psi[[k]]
Z_i <- compute_factor_scores_optimized(X_i, mu_k, Lambda_k, Psi_k)
dim(Z_i)  # (T x r)

## End(Not run)

Description

This function takes a fitted MFA model (with elements z, mu, Lambda, Psi) and reconstructs all subjects' data from the factor scores, computing the ratio of explained variance ($1 - \mathrm{SSE}/\mathrm{SST}$).

Usage

compute_global_vaf_mfa(list_of_data, res)

Arguments

Details

For each subject i, we identify its cluster assignment k_i, compute factor scores via compute_factor_scores, reconstruct Xhat from the factor model, accumulate the sum of squares of residuals (SSE), and the sum of squares total (SST).

Value

A numeric scalar for the global variance accounted for, $1 - \mathrm{SSE}/\mathrm{SST}$.

See Also

compute_factor_scores

Examples

## Not run: 
gvaf <- compute_global_vaf_mfa(my_data_list, fit_mfa)
cat("Global VAF:", gvaf, "\n")

## End(Not run)

Description

Similar to compute_global_vaf_mfa, but uses a Mixture PCA model structure. Given a model with z, mu, P, D, Psi, we reconstruct each subject's data from the principal component scores, accumulate SSE, and compare to total sum of squares (SST).

Usage

compute_global_vaf_mpca(list_of_data, res)

Arguments

Details

The user must define how to reconstruct data from PCA parameters. Typically, for each cluster k, we have W = P * D, so we either compute posterior scores or do a direct projection.

Value

A numeric scalar, $1 - \mathrm{SSE}/\mathrm{SST}$ across all subjects and time steps.

Description

Unified interface for computing ICL (Integrated Complete-data Likelihood) for either MFA or MPCA models.

Usage

compute_ICL(loglik, K, r, M, N_obs, resp, model_type)

Arguments

Value

A numeric ICL value.

See Also

compute_ICL_mfa, compute_ICL_mpca

Description

ICL (Integrated Complete-data Likelihood) is a model selection criterion that penalizes uncertain cluster assignments more heavily than BIC. This makes it more conservative for selecting the number of clusters K.

Usage

compute_ICL_mfa(loglik, K, r, M, N_obs, resp)

Arguments

Details

ICL is computed as:

$ICL = BIC + 2 \times Entropy$

where Entropy = -sum(resp * log(resp)).

The entropy penalty discourages models where cluster assignments are uncertain, which typically happens when K is overestimated (extra clusters with low membership). This makes ICL more conservative than BIC for K selection.

Value

A numeric ICL value.

References

Biernacki, C., Celeux, G., & Govaert, G. (2000). Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(7), 719-725.

See Also

compute_BIC_mfa, compute_entropy

Examples

## Not run: 
icl <- compute_ICL_mfa(fit$logLik, K = 2, r = 3, M = 8,
                       N_obs = 1000, resp = fit$resp)

## End(Not run)

Description

ICL (Integrated Complete-data Likelihood) is a model selection criterion that penalizes uncertain cluster assignments more heavily than BIC. This makes it more conservative for selecting the number of clusters K.

Usage

compute_ICL_mpca(loglik, K, r, M, N_obs, resp)

Arguments

Details

ICL is computed as:

$ICL = BIC + 2 \times Entropy$

where Entropy = -sum(resp * log(resp)).

Value

A numeric ICL value.

See Also

compute_BIC_mpca, compute_entropy

Description

Unified interface for computing log-likelihood for both MFA and MPCA models. Automatically detects model type from the fit object structure.

Usage

compute_logLik(list_of_data, fit, model_type = NULL)

Arguments

Value

A numeric scalar representing the total log-likelihood

See Also

compute_logLik_mfa, compute_logLik_mpca

Examples

## Not run: 
# Works with either MFA or MPCA fit objects
ll <- compute_logLik(list_of_data, fit)

# Or explicitly specify model type
ll <- compute_logLik(list_of_data, fit, model_type = "MFA")

## End(Not run)

Description

Given a fitted MFA model (with K clusters, each containing Lambda[[k]], mu[[k]], Psi[[k]], and mixing proportion pi[k]), this function computes the standard mixture log-likelihood:

$\sum_{i=1}^{N} \log\left( \sum_{k=1}^{K} \pi_k \prod_{t=1}^{T_i} \mathcal{N}(x_{i,t} \mid \mu_k, \Sigma_k) \right),$

where Sigma_k = Lambda_k * Lambda_k^T + Psi_k.

Usage

compute_logLik_mfa(list_of_data, mfa_fit)

Arguments

Details

Internally uses rowwise multivariate normal densities. You will need mvtnorm::dmvnorm or an equivalent. Ensure each Sigma_k is well-defined and invertible if needed.

Value

A numeric scalar, the total log-likelihood.

Description

Given a fitted Mixture PCA model with $P[[k]], $mu[[k]], $sigma2[k], and $pi[k], we compute

$\sum_{i=1}^{N} \log\Bigl(\sum_{k=1}^{K} \pi_k \prod_{t=1}^{T_i}\mathcal{N}\bigl(x_{i,t}\mid \mu_k,\Sigma_k\bigr)\Bigr)$

where Sigma_k = W_k * W_k^T + sigma2_k * I. A typical approach is to reconstruct Sigma_k from P[[k]] if needed (or W[[k]]).

Usage

compute_logLik_mpca(list_of_data, mpca_fit)

Arguments

Details

We assume each cluster k has Sigma_k = W_k * W_k^T + sigma2_k * I. You must assemble that to call mvtnorm::dmvnorm row by row.

Value

A numeric scalar, the total mixture log-likelihood.

Description

Computes Maximum Voluntary Contraction (MVC) reference values from MVC trial data for use in EMG normalization.

Usage

compute_mvc_reference(
  mvc_data,
  method = "peak",
  window_ms = 500,
  top_percent = 5,
  sampling_rate = NULL
)

Arguments

Value

Numeric vector of MVC values, one per channel.

Examples

# From single MVC trial
mvc_trial <- matrix(abs(rnorm(1000 * 4)), ncol = 4)
mvc_ref <- compute_mvc_reference(mvc_trial)

# From multiple MVC trials (takes maximum across trials)
mvc_trials <- list(
  trial1 = matrix(abs(rnorm(1000 * 4)), ncol = 4),
  trial2 = matrix(abs(rnorm(1000 * 4)), ncol = 4)
)
mvc_ref <- compute_mvc_reference(mvc_trials)

Description

Unified interface for computing reconstruction error for both MFA and MPCA models. Automatically detects model type from the fit object structure.

Usage

compute_reconstruction_error(
  list_of_data,
  fit,
  model_type = NULL,
  metric = c("sse", "rmse")
)

Arguments

Value

A numeric scalar representing the reconstruction error

Examples

## Not run: 
# Compute SSE
sse <- compute_reconstruction_error(list_of_data, fit)

# Compute RMSE
rmse <- compute_reconstruction_error(list_of_data, fit, metric = "rmse")

## End(Not run)

Description

Computes the temporal smoothness penalty for factor scores using second-order differences. This penalty measures how "wiggly" the factor activation profiles are, with higher values indicating more high-frequency noise.

Usage

compute_smoothness_penalty(factor_scores)

Arguments

Details

The smoothness penalty is computed as:

$\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$

where $z_{i,j,t}$ is the j-th factor score for subject i at time t.

This is equivalent to summing the squared second-order differences across all time series. When r is too large, models tend to fit high-frequency noise, producing wiggly factors with high smoothness penalties.

Value

A numeric scalar representing the total smoothness penalty across all subjects and factors.

Examples

## Not run: 
# Compute factor scores first
mfa_fit <- compute_factor_scores_mfa(mfa_fit, list_of_data)
# Then compute smoothness
smooth_val <- compute_smoothness_penalty(mfa_fit$factor_scores)

## End(Not run)

Description

Computes the similarity between two subspaces spanned by basis matrices using QR decomposition and SVD. This metric measures how well the estimated subspace recovers the true subspace, independent of rotation and scaling.

Usage

compute_subspace_similarity(L_true, L_hat)

Arguments

Details

The subspace similarity is computed as follows:

1. Orthonormalize both bases using QR decomposition

2. Compute the SVD of Q_true' * Q_hat

3. Return the mean of squared singular values (up to min_dim)

This metric is invariant to rotation and scaling of the basis vectors, making it suitable for comparing factor analysis solutions where rotation indeterminacy exists.

When r_true != r_hat, the metric uses min(r_true, r_hat) dimensions, which allows comparison even when the estimated rank differs from truth.

Value

Numeric scalar in [0, 1], where 1 indicates identical subspaces. Returns NA if inputs are invalid (NULL, empty, or contain NA/NaN/Inf).

Examples

## Not run: 
# Same subspace (rotated)
L_true <- matrix(c(1, 0, 0, 1), 2, 2)
R <- matrix(c(cos(pi/4), -sin(pi/4), sin(pi/4), cos(pi/4)), 2, 2)
L_hat <- L_true %*% R
compute_subspace_similarity(L_true, L_hat)  # Should be ~1

# Orthogonal subspaces
L_true <- matrix(c(1, 0, 0, 0), 2, 2)
L_hat <- matrix(c(0, 0, 0, 1), 2, 2)
compute_subspace_similarity(L_true, L_hat)  # Should be ~0

## End(Not run)

Description

Unified interface for computing global VAF for both MFA and MPCA models. Automatically detects model type from the fit object structure.

Usage

compute_VAF(list_of_data, fit, model_type = NULL)

Arguments

Details

VAF is computed as 1 - SSE/SST, where SSE is the sum of squared errors between observed and reconstructed data, and SST is the total sum of squares.

Value

A numeric scalar representing the global VAF (0 to 1)

See Also

compute_global_vaf_mfa, compute_global_vaf_mpca

Examples

## Not run: 
# Works with either MFA or MPCA fit objects
vaf <- compute_VAF(list_of_data, fit)

# Or explicitly specify model type
vaf <- compute_VAF(list_of_data, fit, model_type = "MFA")

## End(Not run)

Description

This function computes the Variance Accounted For (VAF) for a fitted MFA model across all subjects, reconstructing each subject's data via factor scores and summing squared errors.

Usage

compute_VAF_mfa(list_of_data, mfa_fit, z)

Arguments

Value

A numeric scalar representing the total VAF over all subjects.

Examples

## Not run: 
sim <- simulate_mixture_data(N = 50, K = 3, r = 2, M = 6, T_each = 100)
fit <- mfa_em_fit(sim$list_of_data, K = 3, r = 2)
vaf <- compute_VAF_mfa(sim$list_of_data, fit, fit$z)
print(vaf)

## End(Not run)

Description

This function computes the Variance Accounted For (VAF) for a fitted Mixture PCA model.

Usage

compute_VAF_mpca(list_of_data, mpca_fit, z)

Arguments

Value

A numeric scalar for the overall VAF.

Examples

## Not run: 
sim <- simulate_mixture_data(N = 50, K = 3, r = 2, M = 6, T_each = 100)
fit <- mixture_pca_em_fit(sim$list_of_data, K = 3, r = 2)
vaf <- compute_VAF_mpca(sim$list_of_data, fit, fit$z)
print(vaf)

## End(Not run)

Description

Computes the global Variance Accounted For (VAF) for a fitted MFA or MPCA model. This is a unified wrapper that dispatches to the appropriate method-specific function.

Usage

compute_vaf_wrapper(method, list_of_data, fit)

Arguments

Details

VAF is computed as $1 - SSE/SST$, where SSE is the sum of squared reconstruction errors and SST is the total sum of squares.

Value

Numeric scalar: VAF value in [0, 1], where 1 indicates perfect reconstruction. Returns NA if computation fails.

Examples

## Not run: 
fit <- mfa_em_fit(list_of_data, K = 3, r = 2)
vaf <- compute_vaf_wrapper("MFA", list_of_data, fit)

## End(Not run)

Description

Computes the Very Simple Structure (VSS) criterion, which evaluates how well factor loadings exhibit simple structure (each variable loads strongly on one factor, weakly on others).

Usage

compute_vss_criterion(
  Lambda,
  threshold = 0.5,
  method = c("threshold", "ratio"),
  rotate = c("none", "varimax")
)

Arguments

Details

The VSS criterion evaluates simple structure in factor loadings. Two methods are available:

**Threshold method (default):** For each variable (muscle), count how many factors have |loading| > threshold. Simple structure is achieved when this count is exactly 1. The VSS score is the proportion of variables across all clusters that exhibit simple structure.

**Ratio method:** For each variable (muscle), compute the concentration ratio max(lambda^2)/sum(lambda^2), which measures how concentrated the loadings are on a single factor. The VSS score is the average concentration ratio across all variables. This method is more robust when loadings are small in absolute magnitude, as it doesn't depend on an arbitrary threshold.

This implementation is inspired by the psych::VSS() function, adapted for mixture models where each cluster has its own loading matrix.

Value

A numeric scalar representing the VSS criterion. Higher values indicate better simple structure. Returns NA if rotation fails or loadings are degenerate.

References

Revelle, W., & Rocklin, T. (1979). Very Simple Structure: An Alternative Procedure For Estimating The Optimal Number Of Interpretable Factors. Multivariate Behavioral Research, 14(4), 403-414.

Examples

## Not run: 
# For an MFA model
vss_val <- compute_vss_criterion(mfa_fit$Lambda, threshold = 0.5)

## End(Not run)

Description

Counts how many tasks have completed by checking partial result files.

Usage

count_completed_tasks(output_dir, n_total = NULL)

Arguments

Value

List with completion statistics.

Description

Returns preset configurations for common EMG preprocessing scenarios.

Usage

emg_preset_config(preset = "default")

Arguments

Value

A list containing preprocessing configuration for each step.

Examples

# Get default configuration
config <- emg_preset_config("default")

# Modify for custom needs
config$notch$freq <- 60  # Change to 60 Hz
config$envelope$cutoff <- 10  # Higher cutoff for faster movements

Description

Comprehensive evaluation of model selection results, computing all relevant metrics including ARI, VAF, BIC, logLik, and selection accuracy.

Usage

evaluate_model_selection(
  list_of_data,
  true_cluster,
  K_true,
  r_true,
  select_result,
  runtime_sec = NA_real_
)

Arguments

Value

A list with evaluation metrics:

K_hat

Estimated number of clusters.

r_hat

Estimated factor dimension.

K_correct

Logical: K_hat == K_true.

r_correct

Logical: r_hat == r_true.

Kr_correct

Logical: both K and r correct.

ARI

Adjusted Rand Index.

VAF

Variance Accounted For.

logLik

Log-likelihood.

BIC

Bayesian Information Criterion.

runtime_sec

Runtime in seconds.

Examples

## Not run: 
# Generate data
data <- generate_synergy_data(N = 30, K = 3, r = 3, seed = 123)

# Run model selection
start_time <- Sys.time()
select_result <- select_K_r(data$list_of_data, method = "MFA", strategy = "grid")
runtime <- as.numeric(difftime(Sys.time(), start_time, units = "secs"))

# Evaluate
eval_result <- evaluate_model_selection(
  data$list_of_data, data$true_cluster,
  K_true = 3, r_true = 3,
  select_result, runtime_sec = runtime
)

## End(Not run)

Description

This function extracts cluster assignment labels from a fitted MFA or MPCA model object. It searches through common field names used by different model implementations.

Usage

extract_cluster_labels(fit)

Arguments

Details

The function searches for cluster labels in the following fields (in order):

• z: Standard field for MFA/MPCA models

• cluster: Alternative naming

• clusters: Plural form

• class: Classification result

• z_hat: Estimated assignments

• cl: Short form

• hard_cluster: Hard assignment from soft clustering

• labels: Generic label field

If the model has a resp (responsibility) matrix but no hard assignments, the function computes hard assignments using MAP (maximum a posteriori).

Value

An integer vector of cluster assignments (1 to K), or NA if cluster labels cannot be extracted.

Examples

## Not run: 
fit <- mfa_em_fit(list_of_data, K = 3, r = 2)
pred_cluster <- extract_cluster_labels(fit)

## End(Not run)

Description

Extracts the amplitude envelope from rectified EMG signals using various methods.

Usage

extract_envelope(
  data,
  sampling_rate,
  method = "lowpass",
  cutoff = 6,
  window_ms = 50,
  filter_order = 4,
  zero_phase = TRUE
)

Arguments

Details

The envelope represents the amplitude modulation of the EMG signal, showing muscle activation patterns over time. Common approaches:

• Low-pass filter: Applies a Butterworth low-pass filter to the rectified signal. Cutoff frequency typically 4-10 Hz.

• RMS: Computes root mean square in a sliding window. Good for preserving amplitude information.

• Moving average: Simple smoothing with a sliding window. Fast but may introduce artifacts at window edges.

Value

EMG envelope signal, same structure as input.

Examples

## Not run: 
# Generate and process EMG
fs <- 1000
raw_emg <- rnorm(1000)
rect_emg <- rectify(raw_emg)

# Low-pass envelope (most common)
env_lp <- extract_envelope(rect_emg, sampling_rate = fs, method = "lowpass")

# RMS envelope
env_rms <- extract_envelope(rect_emg, sampling_rate = fs,
                            method = "rms", window_ms = 100)

## End(Not run)

Description

Converts the summary from select_K_r to a flat data frame, useful for saving and plotting.

Usage

extract_selection_summary(select_result)

Arguments

Value

A data frame with columns K, r, logLik, BIC, and optionally stage (for twostage strategy).

Description

This function applies PCA to the concatenated data across all subjects, then for each subject i, computes the mean of the top r_dim principal component scores as a single feature vector. This is useful for initializing cluster assignments in mixture models.

Usage

extract_subject_features_by_singlePCA(list_of_data, r_dim = 2)

Arguments

Details

The function first concatenates all subject data into a single matrix, performs PCA on this combined dataset, and then computes the mean of the PC scores for each subject's time points. This provides a low-dimensional representation of each subject that can be used for clustering initialization.

Value

A numeric matrix of shape (N x r_dim), where N=length(list_of_data). Each row contains the mean PC scores for one subject.

Examples

## Not run: 
# Create example data
N <- 10
M <- 8
list_of_data <- lapply(1:N, function(i) matrix(rnorm(100 * M), 100, M))

# Extract 2D features for each subject
features <- extract_subject_features_by_singlePCA(list_of_data, r_dim = 2)
dim(features)  # Should be 10 x 2

## End(Not run)

Description

Helper function to apply consistent facet_grid() styling for cluster-factor plots.

Usage

facet_mixsynergy(rows = NULL, cols = NULL, scales = "fixed", ...)

Arguments

Value

A ggplot2 facet_grid object

Examples

## Not run: 
library(ggplot2)
ggplot(data, aes(x, y)) + 
  geom_line() + 
  facet_mixsynergy(rows = "Cluster", cols = "Factor", scales = "free_y")

## End(Not run)

Description

This function concatenates all subjects' data into one large matrix, then performs a naive factor analysis. It returns a simple BIC measure and the fitted loadings, mean, and diagonal Psi.

Usage

fit_single_factor_analysis(list_of_data, r = 2)

Arguments

Value

A list with elements (loadings, mu, diagPsi, logLik, BIC).

Examples

## Not run: 
sim <- simulate_mixture_data(N = 50, K = 3, r = 2, M = 6, T_each = 100)
fa_res <- fit_single_factor_analysis(sim$list_of_data, r = 2)
print(fa_res$BIC)

## End(Not run)

Description

Concatenates all data into one matrix, performs SVD, keeps r principal components, then computes a naive PPCA logLik and BIC.

Usage

fit_single_pca(list_of_data, r = 2)

Arguments

Value

A list with (P, mu, sigma2, logLik, BIC).

Examples

## Not run: 
sim <- simulate_mixture_data(N = 50, K = 3, r = 2, M = 6, T_each = 100)
pca_res <- fit_single_pca(sim$list_of_data, r = 2)
print(pca_res$BIC)

## End(Not run)

Description

This function implements a two-step baseline approach for comparison with mixture-based methods:

1. Step 1: Individual Decomposition - Apply PCA or FA to each subject's data separately to extract r basis vectors per subject.

2. Step 1.5: Basis Alignment - Align basis vectors across subjects using the Hungarian algorithm to ensure corresponding components are in the same order, improving clustering accuracy.

3. Step 2: Clustering - Normalize the basis vectors, concatenate them across all subjects, and apply k-means clustering. Assign each subject to a cluster based on majority vote of their basis vectors.

Usage

fit_two_step_baseline(
  data_list,
  r,
  K,
  method = c("PCA", "FA"),
  seed = 123,
  z_true = NULL,
  refine_assignments = FALSE,
  align_basis = TRUE
)

Arguments

Details

The algorithm addresses sign indeterminacy by normalizing each basis vector so that its maximum absolute value element is positive. Tie-breaking in majority vote is handled by computing the mean distance from a subject's basis vectors to each tied cluster centroid and selecting the closer one.

Step 1.5 (basis alignment) uses the Hungarian algorithm to align basis vectors across subjects. This ensures that corresponding components (e.g., the first component of each subject) represent similar patterns, which is crucial for accurate clustering. Without alignment, the arbitrary ordering of components from individual PCA/FA can lead to poor clustering performance.

Step 4 (re-optimization) fits cluster-specific models in standardized space to ensure consistency with Step 1. This improves reconstruction metrics (SSE/VAF) but does not change ARI unless refine_assignments = TRUE.

If refine_assignments = TRUE, a single reassignment pass is performed after fitting cluster-specific models. Each subject is reassigned to the cluster that minimizes reconstruction error. This can improve ARI while maintaining the two-step nature of the baseline (unlike iterative EM which would drift toward the mixture model).

For FA decomposition, the function uses psych::fa() if available. If the psych package is not installed or if FA fails for a subject, the function falls back to PCA with a warning.

Value

A list with:

z_est

Integer vector of estimated cluster assignments (1..K).

W_list

List of normalized basis matrices, one per subject.

W_centroids

Matrix of cluster centroids in basis space.

mu_list

List of centering vectors used for standardization, one per subject.

sd_list

List of scaling vectors used for standardization, one per subject.

cluster_models

List of cluster-specific models (each with W and n_members).

method

Character string indicating the method used (e.g., "TwoStep_PCA").

BIC

Set to NA (k-means does not have a canonical BIC).

ARI

Adjusted Rand Index if z_true is provided, otherwise NA.

Examples

## Not run: 
sim <- simulate_mixture_data(N = 50, K = 3, r = 2, M = 6, T_each = 100)
result <- fit_two_step_baseline(
  sim$list_of_data, r = 2, K = 3, method = "PCA",
  seed = 123, z_true = sim$z_true
)
print(result$ARI)

## End(Not run)

Description

Generates a PBS array job script for submitting parameter sweep tasks.

Usage

generate_pbs_script(
  n_tasks,
  output_dir,
  pbs_config = list(),
  sweep_type = c("modelsel", "baselines"),
  sweep_config = list()
)

Arguments

Value

Character. Path to the generated PBS script.

Examples

## Not run: 
script_path <- generate_pbs_script(
  n_tasks = 100,
  output_dir = "results",
  pbs_config = list(
    queue = "batch",
    walltime = "02:00:00",
    r_module = "R/4.3.0"
  )
)

## End(Not run)

Description

This function provides a standardized interface for generating v2 dynamic synergy data suitable for model selection experiments. It wraps simulate_dynamic_synergy_data and returns data in a format optimized for parameter sweep experiments with model selection.

Usage

generate_synergy_data(
  N = 30,
  K = 3,
  r = 3,
  M = 8,
  T_each = 200,
  cluster_sep_spatial = 1,
  cluster_sep_temporal = 1,
  cluster_sep_stability = 1,
  freq_range = c(0.5, 3),
  amp_range = c(0.5, 1.5),
  base_level = 0.05,
  rectify_output = TRUE,
  seed
)

Arguments

Details

This function is designed for use in parameter sweep experiments where model selection (K, r) is performed on generated data. The key differences from simulate_dynamic_synergy_data are:

• Returns list_of_data (renamed from X_list) for consistency with fitting functions

• Returns true_cluster (renamed from z_true) for consistency with evaluation functions

• Includes a params list for tracking generation parameters

• Seed is required (not optional) to ensure reproducibility

Value

A list with the following components:

list_of_data

List of length N; each element is a (T_each x M) matrix of generated EMG signals. This is the standard input format for MFA/MPCA fitting.

true_cluster

Integer vector of length N indicating true cluster assignments (1 to K). Named as true_cluster for consistency.

true_Lambda_list

List of length K; each element is a (M x r) matrix of true cluster-specific synergy loadings. Essential for subspace recovery evaluation. This is the canonical name for ground truth basis matrices.

params

List containing all generation parameters for reproducibility. Includes both canonical short names (T, sep_spatial, sep_temporal, sep_stability) and full names (T_each, cluster_sep_spatial, etc.) for compatibility. Also includes r_mode = "global" for future extension.

Lambda_list

List of length K; each element is a (M x r) matrix of cluster-specific synergy loadings. Kept for backward compatibility; identical to true_Lambda_list.

C_list

List of length N; each element is a (T_each x r) matrix of synergy activation signals (for visualization/validation).

spatial_params

List containing spatial structure parameters for each cluster.

temporal_params

List containing temporal coordination parameters for each cluster.

stability_params

List containing stability-adaptability parameters for each cluster.

See Also

simulate_dynamic_synergy_data for the underlying simulation function, run_param_sweep_modelsel for parameter sweep experiments.

Examples

## Not run: 
# Generate data for model selection experiment
data <- generate_synergy_data(
  N = 30, K = 3, r = 3, M = 8, T_each = 200,
  cluster_sep_spatial = 1.0,
  cluster_sep_temporal = 1.0,
  cluster_sep_stability = 1.0,
  seed = 123
)

# Access the data for fitting
list_of_data <- data$list_of_data
true_cluster <- data$true_cluster

# Fit MFA model
fit <- mfa_em_fit(list_of_data, K = 3, r = 3)

# Evaluate clustering accuracy
ari <- mclust::adjustedRandIndex(true_cluster, fit$z)

## End(Not run)

Description

After calling select_optimal_K_r_mfa (or a similar function), you have $all_models containing fits for various (K, r). This function extracts the model for a user-specified pair (K_target, r_target).

Usage

get_model_by_K_r(selection_obj, K_target, r_target)

Arguments

Value

The fitted model (the same structure returned by mfa_em_fit), or NULL if no match found.

Examples

## Not run: 
sel_obj <- select_optimal_K_r_mfa(my_data, Kvec=2:4, rvec=1:3)
best_model <- get_model_by_K_r(sel_obj, K_target=3, r_target=2)

## End(Not run)

Description

Given a selection_obj from select_optimal_K_r_mpca, extracts the model for the specified (K_target, r_target) if it exists.

Usage

get_model_by_K_r_mpca(selection_obj, K_target, r_target)

Arguments

Value

The fitted model (same structure as from mixture_pca_em_fit), or NULL.

Description

This function returns a small data.frame that describes the parameter ranges used in the simulation (e.g., K, r, N, M, T_each, etc.).

Usage

make_table1()

Value

A data.frame with columns Parameter and Values.

Examples

## Not run: 
table1 <- make_table1()
print(table1)

## End(Not run)

Description

This function filters the data for certain separation/noise conditions (e.g., sep=1.0 and noise=1.0) and computes group-wise means and standard deviations for ARI and BIC. It returns a small summary table.

Usage

make_table2(df)

Arguments

Value

A data.frame summarizing ARI/BIC means and sds, grouped by (K_true, r_true, N, M).

Examples

## Not run: 
df <- read_all_results("results")
table2 <- make_table2(df)
print(table2)

## End(Not run)

Description

This function calls the C++ function mfaTimeseriesCpp() to perform a Mixture Factor Analysis EM algorithm for fixed K and r. By default, it runs a single pass with the internal initialization from C++. However, if you specify multiple initial attempts (via n_init>1) and/or use_kmeans_init=TRUE, this function will try several different initial cluster assignments (in parallel using mclapply), then return the best solution (maximizing the final log-likelihood).

Usage

mfa_em_fit(
  list_of_data,
  K,
  r,
  max_iter = 50,
  iter_fa = 20,
  tol = 0.001,
  n_init = 1,
  use_kmeans_init = TRUE,
  kmeans_rdim = r,
  mc_cores = 1,
  seed = NULL,
  verbose = TRUE,
  n_threads = 0,
  cpp_seed = 0,
  progress_callback = NULL,
  ...
)

Arguments

Details

The single-run approach in C++ already does a built-in initialization if z_init is not provided. Therefore, if n_init=1 and use_kmeans_init=FALSE, we just call mfaTimeseriesCpp once (like the original design). Otherwise:

1. If use_kmeans_init=TRUE, we do one run where we assign clusters by k-means on some subject-level features (extracted by PCA).

2. We also sample n_init random assignments (each subject assigned to a random cluster).

3. For each assignment, we call mfa_em_fit_cpp_singleInit.

4. We compare all results' logLik and keep the best one.

This provides a more robust initialization scheme to avoid poor local maxima.

Value

A list with the same elements as a single run: z, pi, mu, Lambda, Psi, logLik, resp. It corresponds to the best solution (i.e. highest log-likelihood) among all tried inits.

See Also

mixture_pca_em_fit for MPCA fitting, select_r for selecting optimal r, plot_cluster_synergy_loadings_mfa for visualization

Examples

## Not run: 
# Suppose we have a list_of_data for N=100 subjects, each (T_i x M=8).
# We want to fit K=3, r=2, and try 5 random inits plus one k-means init:
fit <- mfa_em_fit(
  list_of_data, K=3, r=2,
  max_iter=50, iter_fa=20, tol=1e-3,
  n_init=5,
  use_kmeans_init=TRUE,
  kmeans_rdim=2,
  mc_cores=2,
  seed=123
)
print(fit$logLik)
head(fit$z)

## End(Not run)

Description

This helper function calls the C++ routine mfaTimeseriesCpp() with a user-specified cluster assignment (z_init) and then computes the final log-likelihood and responsibilities in R.

Usage

mfa_em_fit_cpp_singleInit(
  list_of_data,
  K,
  r,
  z_init,
  max_iter = 50,
  nIterFA = 20,
  tol = 0.001,
  verbose = TRUE,
  n_threads = 0,
  seed = 0,
  progress_callback = NULL,
  em_type = c("hard", "soft")
)

Arguments

Value

A list with elements:

z

Hard cluster assignments, length N.

pi

Cluster mixing proportions, length K.

mu

List of length K, each a mean vector (M x 1).

Lambda

List of length K, each an M x r factor loading matrix.

Psi

List of length K, each an M x M diagonal noise matrix.

logLik

Final log-likelihood computed in R.

resp

An N x K matrix of responsibilities.

Examples

## Not run: 
# Suppose we have list_of_data, K=3, r=2:
z_init <- sample.int(3, size=length(list_of_data), replace=TRUE)
fit_one <- mfa_em_fit_cpp_singleInit(list_of_data, K=3, r=2,
                                     z_init=z_init)
fit_one$logLik

## End(Not run)

Description

This function calls the C++ function mpcaTimeseriesCpp() to perform a Mixture PCA (PPCA) EM algorithm for fixed K and r, then converts the W matrices into (P, D), creates Psi, and computes the final log-likelihood/responsibilities in R. By default, it runs a single pass with the built-in initialization in C++. However, if n_init > 1 or use_kmeans_init=TRUE, it performs multiple initializations in parallel and returns the best-fitting result (highest log-likelihood).

Usage

mixture_pca_em_fit(
  list_of_data,
  K,
  r,
  max_iter = 50,
  iter_pca = 20,
  tol = 0.001,
  method = "EM",
  n_init = 1,
  use_kmeans_init = TRUE,
  kmeans_rdim = r,
  mc_cores = 1,
  seed = NULL,
  verbose = TRUE,
  n_threads = 0,
  cpp_seed = 0,
  progress_callback = NULL,
  ...
)

Arguments

Details

If n_init=1 and use_kmeans_init=FALSE, this function runs exactly one pass with the default initialization in C++ (i.e. subject i is assigned to cluster (i % K) + 1). Otherwise:

1. If use_kmeans_init=TRUE, we do a run where we apply k-means to some subject-level PCA features to get z_init, then run EM once.

2. We generate n_init random initial assignments, run EM for each in parallel, and collect the results.

3. We compare all solutions by final log-likelihood and pick the best.

Value

A list with elements:

See Also

mfa_em_fit for MFA fitting, select_r for selecting optimal r, plot_cluster_synergy_loadings_mpca for visualization

Examples

## Not run: 
# Suppose we have a list_of_data, K=3, r=2, and we want to try 5 random inits + k-means:
fit <- mixture_pca_em_fit(
  list_of_data, K=3, r=2, max_iter=50, iter_pca=20, tol=1e-3, method="EM",
  n_init=5, use_kmeans_init=TRUE, kmeans_rdim=2, mc_cores=2
)
print(fit$logLik)
head(fit$z)

## End(Not run)

Description

This helper function calls the updated C++ routine mpcaTimeseriesCpp() with a user-specified cluster assignment (z_init), then converts the returned W into (P, D) and constructs Psi, and finally computes the final log-likelihood and responsibilities in R.

Usage

mixture_pca_em_fit_cpp_singleInit(
  list_of_data,
  K,
  r,
  max_iter = 50,
  nIterPCA = 20,
  tol = 0.001,
  method = "EM",
  z_init,
  verbose = TRUE,
  n_threads = 0,
  seed = 0,
  progress_callback = NULL,
  em_type = c("hard", "soft")
)

Arguments

Value

A list with elements:

Examples

## Not run: 
# Suppose we have list_of_data, K=3, r=2, method="EM", plus some z_init:
z_init <- sample.int(3, size=length(list_of_data), replace=TRUE)
fit_one <- mixture_pca_em_fit_cpp_singleInit(
  list_of_data, K=3, r=2, max_iter=50, nIterPCA=20, tol=1e-3,
  method="EM", z_init=z_init
)
print(fit_one$logLik)
# Now we can call compute_factor_scores_mpca(fit_one, list_of_data) because fit_one$W exists

## End(Not run)

Description

Normalizes EMG signals using various methods to allow comparison across channels, subjects, or sessions.

Usage

normalize_emg(data, method = "max", reference = NULL, per_channel = TRUE)

Arguments

Details

Normalization is essential for:

• Comparing EMG across different muscles

• Comparing EMG across different subjects

• Comparing EMG across different recording sessions

MVC normalization (method = "mvc") is the gold standard in EMG research, expressing activation as a percentage of maximum voluntary contraction. Requires separate MVC trials for reference.

Max normalization (method = "max") is useful when MVC data is not available, normalizing to the peak value in the data.

Value

Normalized EMG signal, same structure as input.

Examples

# Max normalization (0 to 1 range)
emg <- matrix(rnorm(1000), ncol = 4)
norm_emg <- normalize_emg(emg, method = "max")

# MVC normalization (% MVC)
mvc_values <- c(100, 150, 80, 120)  # MVC for each channel
norm_emg <- normalize_emg(emg, method = "mvc", reference = mvc_values)

# Z-score normalization
norm_emg <- normalize_emg(emg, method = "zscore")

Description

Applies a notch (band-stop) filter to remove power line interference (50 Hz or 60 Hz) and optionally its harmonics.

Usage

notch_filter(
  data,
  sampling_rate,
  notch_freq = 50,
  Q = 30,
  harmonics = TRUE,
  zero_phase = TRUE
)

Arguments

Details

Power line interference (50/60 Hz) is a common artifact in EMG recordings. The notch filter creates a narrow band-stop at the specified frequency. When harmonics = TRUE, harmonics (100/120 Hz, 150/180 Hz, etc.) are also removed.

Value

Filtered EMG signal with power line noise removed.

Examples

## Not run: 
# Remove 50 Hz power line noise
clean_emg <- notch_filter(emg, sampling_rate = 1000, notch_freq = 50)

# Remove 60 Hz (Americas) with harmonics
clean_emg <- notch_filter(emg, sampling_rate = 1000, notch_freq = 60, harmonics = TRUE)

## End(Not run)

Description

Visualizes temporal activation curves per cluster by showing representative activation patterns from subjects in each cluster.

Usage

plot_activation_patterns(sim, n_subjects_per_cluster = 3, theme_base_size = 14)

Arguments

Value

A ggplot object showing temporal activation curves per cluster.

Examples

## Not run: 
sim <- simulate_dynamic_synergy_data(N=30, K=3, r=3, M=8, T_each=200)
plot_activation_patterns(sim)

## End(Not run)

Description

Creates heatmaps of basis vectors (W matrices) for each subject and component, showing how well the Two-step alignment procedure worked. Compares aligned basis vectors with true synergies (if available) to diagnose alignment failures.

Usage

plot_aligned_basis_heatmaps(
  W_list,
  W_true = NULL,
  z_true = NULL,
  z_est = NULL,
  max_subjects = 30,
  theme_base_size = 10
)

Arguments