SIMCA — Soft Independent Modelling of Class Analogies
| Task | Classification (class modelling) |
Method key (analysis_type) | simca |
| Prediction on new data | Yes — see Prediction |
| Library | Own class SIMCA in analysis.py, built on sklearn.decomposition.PCA, StandardScaler and scipy.stats.f |
When to use
- Each class is a well-defined, homogeneous group (for example, one product or one analyte), and you want to describe each class on its own.
- Classes have different internal structure: SIMCA fits a separate model to each of them, so one class may vary along other directions than another.
- Many correlated features, few samples per class.
How it works
SIMCA builds one PCA model per class and assigns a sample to the class whose model reconstructs it best.
- One PCA model per class. For class with samples and features, the class data are centred (and scaled, with Scale Features on) and a PCA with components is fitted.
- Orthogonal distance. A sample is projected onto the class model and reconstructed as . The squared orthogonal distance (the residual, in the original units) is
- F statistic. The distance is compared with the typical residual of the training samples of the class:
- Assignment. The sample goes to the class with the smallest .
Probabilities are a softmax over . For AUC, Chrometrica uses the distances directly (for two classes: ), because the softmax is almost flat on many data sets and would give an AUC of exactly 0.5.
:::note Differences from textbook SIMCA
- Classic SIMCA also accepts or rejects a sample for every class separately, with a critical value , so a sample can belong to several classes or to none. Here the critical value is computed and stored with the model, but not used for the decision: every sample gets exactly one class, the nearest one. Significance Level (α) therefore does not change the predictions in the current version.
- Only the orthogonal distance is used. The score distance (Hotelling inside the model plane) is not taken into account.
:::
Parameters
Defaults as set by the analysis dialog.
| Parameter | UI label | Default | Notes |
|---|---|---|---|
n_components | Number of Components | , = number of feature columns | Components of each class model; reduced per class to |
alpha | Significance Level (α) | 0.05 | See the note above |
scale_x | Scale Features | on | Autoscale the data of each class with the mean and standard deviation of that class |
feature_selection | Enable Feature Selection | off | See Feature selection. |
average_replicates | Average technical replicates | off | See Averaging replicates. |
Preprocessing
- Each class model has its own centring and scaling, fitted on the training rows of that class.
- Every class needs at least 2 rows in every training fold. If a fold cannot be fitted, the analysis does not stop: the fold gets accuracy 0 and all its test rows are predicted as the first training label. Such folds pull the CV metrics down; if CV accuracy looks surprisingly low, check the class sizes.
- With Group K-Fold, the number of folds is reduced to the number of samples if it is larger.
- Before the model: replicate averaging and feature selection, if they are on. See Order of steps.
Results and metrics
The common classification metrics are described in Metrics. AUC is computed from the distances as described above; for more than two classes it is the macro-averaged one-vs-rest AUC.
SIMCA-specific fields:
| Field | Meaning |
|---|---|
prediction_probabilities | Softmax of for every sample (final model) |
model_info | n_components, alpha, scale_x, n_classes |
Visualizations
- Confusion matrix
- Permutation plot
- Time series and heatmap of the input data
Source code
The model is the class SIMCA in chrometrica/analysis/analysis.py. Fitting
one PCA model per class:
analysis.py · SIMCA.fit() · lines 332–412
def fit(self, X, y):
"""Fit SIMCA models for each class with numerical stability"""
self.classes_ = np.unique(y)
if len(self.classes_) < 2:
raise ValueError("SIMCA requires at least 2 classes")
for class_label in self.classes_:
# Get data for this class
X_class = X[y == class_label]
if X_class.shape[0] < 2:
raise ValueError(
f"Class {class_label} has insufficient samples ({X_class.shape[0]}). Need at least 2 samples per class.")
# Create and fit PCA model for this class
if self.scale_x:
scaler = StandardScaler(with_mean=True, with_std=True)
X_scaled = scaler.fit_transform(X_class)
else:
scaler = None
X_scaled = X_class
# Ensure n_components is valid
n_comp = min(self.n_components, X_class.shape[0] - 1, X_class.shape[1])
if n_comp < 1:
n_comp = 1
pca = PCA(n_components=n_comp, random_state=0)
pca.fit(X_scaled)
# Compute critical F value with bounds checking
N = X_class.shape[0] # Number of samples
J = X_class.shape[1] # Number of features
K = n_comp
# Determine degrees of freedom parameter with bounds
a = min(J, N - 1) # More conservative approach
if a <= K:
a = K + 1 # Ensure a > K for valid F distribution
# Critical F value with bounds checking
try:
df1 = max(1, a - K) # Ensure degrees of freedom >= 1
df2 = max(1, (a - K) * (N - K - 1)) # Ensure degrees of freedom >= 1
f_crit = stats.f.ppf(1.0 - self.alpha, df1, df2)
# Ensure finite F critical value
if not np.isfinite(f_crit) or f_crit <= 0:
f_crit = 1.0 # Fallback value
except (ValueError, RuntimeError):
f_crit = 1.0 # Fallback value
# Precompute training residuals for efficiency and stability
X_train_scaled = scaler.transform(X_class) if scaler else X_class
X_train_pred = pca.inverse_transform(pca.transform(X_train_scaled))
X_train_pred = scaler.inverse_transform(X_train_pred) if scaler else X_train_pred
od_train_squared = np.sum((X_class - X_train_pred) ** 2, axis=1)
total_od_train = np.sum(od_train_squared)
# Avoid division by zero in denominator
denominator = total_od_train / max(1, (a - K) * (N - K - 1))
if denominator <= 0:
denominator = 1e-10 # Small positive value to avoid division by zero
# Store model components
self.models[class_label] = {
'scaler': scaler,
'pca': pca,
'N': N,
'J': J,
'K': K,
'a': a,
'f_crit': f_crit,
'X_train': X_class,
'denominator': denominator,
'total_od_train': total_od_train
}
return self
Distance to a class and assignment:
analysis.py · SIMCA._compute_distance() · lines 414–454
def _compute_distance(self, X, class_label):
"""Compute orthogonal distance for samples to a specific class"""
model = self.models[class_label]
distances = []
for i in range(X.shape[0]):
sample = X[i:i + 1]
# Transform sample
if model['scaler'] is not None:
X_scaled = model['scaler'].transform(sample)
else:
X_scaled = sample
# Reconstruct sample
try:
X_pred = model['pca'].inverse_transform(
model['pca'].transform(X_scaled)
)
if model['scaler'] is not None:
X_pred = model['scaler'].inverse_transform(X_pred)
# Calculate orthogonal distance with numerical stability
od_squared = np.sum((sample - X_pred) ** 2)
# F-test statistic with bounds checking
numerator = od_squared / max(1, model['a'] - model['K'])
f_stat = numerator / max(1e-10, model['denominator']) # Avoid division by zero
# Ensure finite F-statistic
if not np.isfinite(f_stat) or f_stat < 0:
f_stat = 1e10 # Large value for invalid cases
except (ValueError, np.linalg.LinAlgError):
# Handle numerical errors in PCA transformation
f_stat = 1e10 # Large value for invalid cases
distances.append(f_stat)
return np.array(distances)
analysis.py · SIMCA.predict() · lines 456–482
def predict(self, X):
"""Predict class membership for each sample with numerical stability"""
X = np.array(X)
if X.ndim == 1:
X = X.reshape(1, -1)
predictions = []
confidence_scores = []
for i in range(X.shape[0]):
sample = X[i:i + 1]
best_class = None
best_confidence = float('inf')
for class_label in self.classes_:
# Compute distance for this sample to current class
f_stat = self._compute_distance(sample, class_label)[0]
# Store the class with the lowest F-statistic (most likely to belong)
if f_stat < best_confidence:
best_confidence = f_stat
best_class = class_label
predictions.append(best_class)
confidence_scores.append(best_confidence)
return np.array(predictions)
Full source of the SIMCA class
analysis.py · SIMCA() · lines 322–539
class SIMCA:
"""Soft Independent Modeling of Class Analogies (SIMCA) implementation with numerical stability"""
def __init__(self, n_components=2, alpha=0.05, scale_x=True):
self.n_components = n_components
self.alpha = alpha
self.scale_x = scale_x
self.models = {}
self.classes_ = None
def fit(self, X, y):
"""Fit SIMCA models for each class with numerical stability"""
self.classes_ = np.unique(y)
if len(self.classes_) < 2:
raise ValueError("SIMCA requires at least 2 classes")
for class_label in self.classes_:
# Get data for this class
X_class = X[y == class_label]
if X_class.shape[0] < 2:
raise ValueError(
f"Class {class_label} has insufficient samples ({X_class.shape[0]}). Need at least 2 samples per class.")
# Create and fit PCA model for this class
if self.scale_x:
scaler = StandardScaler(with_mean=True, with_std=True)
X_scaled = scaler.fit_transform(X_class)
else:
scaler = None
X_scaled = X_class
# Ensure n_components is valid
n_comp = min(self.n_components, X_class.shape[0] - 1, X_class.shape[1])
if n_comp < 1:
n_comp = 1
pca = PCA(n_components=n_comp, random_state=0)
pca.fit(X_scaled)
# Compute critical F value with bounds checking
N = X_class.shape[0] # Number of samples
J = X_class.shape[1] # Number of features
K = n_comp
# Determine degrees of freedom parameter with bounds
a = min(J, N - 1) # More conservative approach
if a <= K:
a = K + 1 # Ensure a > K for valid F distribution
# Critical F value with bounds checking
try:
df1 = max(1, a - K) # Ensure degrees of freedom >= 1
df2 = max(1, (a - K) * (N - K - 1)) # Ensure degrees of freedom >= 1
f_crit = stats.f.ppf(1.0 - self.alpha, df1, df2)
# Ensure finite F critical value
if not np.isfinite(f_crit) or f_crit <= 0:
f_crit = 1.0 # Fallback value
except (ValueError, RuntimeError):
f_crit = 1.0 # Fallback value
# Precompute training residuals for efficiency and stability
X_train_scaled = scaler.transform(X_class) if scaler else X_class
X_train_pred = pca.inverse_transform(pca.transform(X_train_scaled))
X_train_pred = scaler.inverse_transform(X_train_pred) if scaler else X_train_pred
od_train_squared = np.sum((X_class - X_train_pred) ** 2, axis=1)
total_od_train = np.sum(od_train_squared)
# Avoid division by zero in denominator
denominator = total_od_train / max(1, (a - K) * (N - K - 1))
if denominator <= 0:
denominator = 1e-10 # Small positive value to avoid division by zero
# Store model components
self.models[class_label] = {
'scaler': scaler,
'pca': pca,
'N': N,
'J': J,
'K': K,
'a': a,
'f_crit': f_crit,
'X_train': X_class,
'denominator': denominator,
'total_od_train': total_od_train
}
return self
def _compute_distance(self, X, class_label):
"""Compute orthogonal distance for samples to a specific class"""
model = self.models[class_label]
distances = []
for i in range(X.shape[0]):
sample = X[i:i + 1]
# Transform sample
if model['scaler'] is not None:
X_scaled = model['scaler'].transform(sample)
else:
X_scaled = sample
# Reconstruct sample
try:
X_pred = model['pca'].inverse_transform(
model['pca'].transform(X_scaled)
)
if model['scaler'] is not None:
X_pred = model['scaler'].inverse_transform(X_pred)
# Calculate orthogonal distance with numerical stability
od_squared = np.sum((sample - X_pred) ** 2)
# F-test statistic with bounds checking
numerator = od_squared / max(1, model['a'] - model['K'])
f_stat = numerator / max(1e-10, model['denominator']) # Avoid division by zero
# Ensure finite F-statistic
if not np.isfinite(f_stat) or f_stat < 0:
f_stat = 1e10 # Large value for invalid cases
except (ValueError, np.linalg.LinAlgError):
# Handle numerical errors in PCA transformation
f_stat = 1e10 # Large value for invalid cases
distances.append(f_stat)
return np.array(distances)
def predict(self, X):
"""Predict class membership for each sample with numerical stability"""
X = np.array(X)
if X.ndim == 1:
X = X.reshape(1, -1)
predictions = []
confidence_scores = []
for i in range(X.shape[0]):
sample = X[i:i + 1]
best_class = None
best_confidence = float('inf')
for class_label in self.classes_:
# Compute distance for this sample to current class
f_stat = self._compute_distance(sample, class_label)[0]
# Store the class with the lowest F-statistic (most likely to belong)
if f_stat < best_confidence:
best_confidence = f_stat
best_class = class_label
predictions.append(best_class)
confidence_scores.append(best_confidence)
return np.array(predictions)
def predict_proba(self, X):
"""Predict class probabilities with numerical stability"""
X = np.array(X)
if X.ndim == 1:
X = X.reshape(1, -1)
n_samples = X.shape[0]
n_classes = len(self.classes_)
probabilities = np.zeros((n_samples, n_classes))
# Compute distances for all samples to all classes
distances = np.zeros((n_samples, n_classes))
for j, class_label in enumerate(self.classes_):
distances[:, j] = self._compute_distance(X, class_label)
# Convert distances to probabilities using softmax with numerical stability
for i in range(n_samples):
# Use negative distances (smaller distance = higher probability)
scores = -distances[i]
# Handle extreme values - clip to prevent overflow
scores = np.clip(scores, -500, 500) # exp(700) is about the limit for float64
# Shift for numerical stability
scores_shifted = scores - np.max(scores)
# Compute softmax
exp_scores = np.exp(scores_shifted)
sum_exp = np.sum(exp_scores)
if sum_exp > 0:
probabilities[i] = exp_scores / sum_exp
else:
# Fallback to uniform distribution
probabilities[i] = np.ones(n_classes) / n_classes
# Ensure no NaN/Inf values
if np.any(np.isnan(probabilities[i])) or np.any(np.isinf(probabilities[i])):
probabilities[i] = np.ones(n_classes) / n_classes
return probabilities
def decision_function(self, X):
"""Return distance to each class (smaller = more likely)"""
X = np.array(X)
if X.ndim == 1:
X = X.reshape(1, -1)
n_samples = X.shape[0]
n_classes = len(self.classes_)
distances = np.zeros((n_samples, n_classes))
for j, class_label in enumerate(self.classes_):
distances[:, j] = self._compute_distance(X, class_label)
return distances
The AUC score:
analysis.py · simca_scores_for_auc() · lines 542–567
def simca_scores_for_auc(model, X, is_binary):
"""Continuous SIMCA decision score suitable for ROC AUC.
``SIMCA.predict_proba`` is a softmax over clipped negative F-distances and
collapses to a near-constant vector (~0.5 per class) on many real datasets,
which makes ``roc_auc_score`` return exactly 0.5 regardless of the data.
This helper instead derives the score from ``decision_function`` (the raw
orthogonal-distance F statistics, one column per class in ``model.classes_``
order, i.e. ascending label).
Binary -> 1-D array, higher == more likely the positive (larger-label)
class: ``log d(neg) - log d(pos)``.
Multiclass -> 2-D array (n_samples, n_classes) == ``-decision_function``,
columns in ``model.classes_`` order, for one-vs-rest AUC.
"""
d = np.asarray(model.decision_function(X), dtype=float)
# F distances are >= 0; guard against inf / NaN from degenerate folds.
d = np.nan_to_num(d, nan=1e12, posinf=1e12, neginf=0.0)
d = np.clip(d, 0.0, 1e12)
if is_binary and d.ndim == 2 and d.shape[1] == 2:
eps = 1e-12
scores = np.log(d[:, 0] + eps) - np.log(d[:, 1] + eps)
return np.nan_to_num(scores, nan=0.0, posinf=1e12, neginf=-1e12)
return np.nan_to_num(-d, nan=0.0, posinf=1e12, neginf=-1e12)
How run_simca_analysis uses the class:
analysis.py · run_simca_analysis() · lines 2722–2917
n_components = parameters.get('n_components', 2)
alpha = parameters.get('alpha', 0.05)
scale_x = parameters.get('scale_x', True)
# ...
if cv_method == 'group':
groups = parameters.get('groups')
if groups is None:
raise ValueError("Groups parameter required for group cross-validation")
cv = GroupKFold(n_splits=min(cv_folds, len(np.unique(groups))))
elif cv_method == 'stratified':
cv = StratifiedKFold(n_splits=cv_folds, shuffle=True, random_state=42)
else:
cv = KFold(n_splits=cv_folds, shuffle=True, random_state=42)
for train_idx, test_idx in cv.split(X_processed, y, groups=groups if cv_method == 'group' else None):
X_train, X_test = X_processed[train_idx], X_processed[test_idx]
y_train, y_test = y[train_idx], y[test_idx]
try:
# Fit SIMCA on training data
simca_cv = SIMCA(n_components=n_components, alpha=alpha, scale_x=scale_x)
simca_cv.fit(X_train, y_train)
# Predict on test data
y_pred = simca_cv.predict(X_test)
# ---------- AUC computation (continuous SIMCA decision score) ----------
# NB: SIMCA.predict_proba is a softmax of clipped negative F-distances
# and collapses to a near-constant (~0.5) on many datasets, which makes
# roc_auc_score return exactly 0.5. Use decision_function distances
# instead (see simca_scores_for_auc).
y_score = simca_scores_for_auc(simca_cv, X_test, is_binary)
# ...
# Fit final model on all data
simca_final_model = None # stays None if the final fit fails (no model bundle)
try:
simca = SIMCA(n_components=n_components, alpha=alpha, scale_x=scale_x)
simca.fit(X_processed, y)
simca_final_model = simca
y_pred_full = simca.predict(X_processed)
y_pred_proba = simca.predict_proba(X_processed)
Full source of run_simca_analysis()
analysis.py · run_simca_analysis() · lines 2708–3068
def run_simca_analysis(X, y, parameters, cv_method, cv_folds):
"""Run SIMCA analysis with numerical stability, proper error handling, and optional feature selection"""
from sklearn.model_selection import StratifiedKFold, GroupKFold, KFold
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score, roc_auc_score
import numpy as np
if parameters.get('average_replicates', False):
groups = parameters.get('groups')
if groups is None:
raise ValueError("'groups' must be provided when averaging replicates is enabled.")
X, y, new_groups = average_replicates(X, y, groups, method='mean')
parameters['groups'] = new_groups
# Get parameters with defaults
n_components = parameters.get('n_components', 2)
alpha = parameters.get('alpha', 0.05)
scale_x = parameters.get('scale_x', True)
# Get overall classes for consistent confusion matrix dimensions
overall_classes = np.unique(y)
overall_classes_list = overall_classes.tolist()
n_overall_classes = len(overall_classes_list)
is_binary = (n_overall_classes == 2) # moved earlier
# Check for feature selection
feature_selection_params = get_feature_selection_params(parameters)
if feature_selection_params:
X_processed, selected_features, feature_scores = select_features(
X, y, feature_selection_params, 'classification'
)
else:
X_processed = X
selected_features = list(range(X.shape[1]))
feature_scores = None
# Enhanced cross-validation with metrics
cv_scores = []
cv_predictions = []
cv_true_labels = []
cv_class_reports = []
cv_confusion_matrices = []
cv_sensitivity_scores = []
cv_specificity_scores = []
# ---------- AUC initialisation ----------
cv_auc_scores = []
cv_scores_all = []
if cv_method == 'group':
groups = parameters.get('groups')
if groups is None:
raise ValueError("Groups parameter required for group cross-validation")
cv = GroupKFold(n_splits=min(cv_folds, len(np.unique(groups))))
elif cv_method == 'stratified':
cv = StratifiedKFold(n_splits=cv_folds, shuffle=True, random_state=42)
else:
cv = KFold(n_splits=cv_folds, shuffle=True, random_state=42)
for train_idx, test_idx in cv.split(X_processed, y, groups=groups if cv_method == 'group' else None):
X_train, X_test = X_processed[train_idx], X_processed[test_idx]
y_train, y_test = y[train_idx], y[test_idx]
try:
# Fit SIMCA on training data
simca_cv = SIMCA(n_components=n_components, alpha=alpha, scale_x=scale_x)
simca_cv.fit(X_train, y_train)
# Predict on test data
y_pred = simca_cv.predict(X_test)
# ---------- AUC computation (continuous SIMCA decision score) ----------
# NB: SIMCA.predict_proba is a softmax of clipped negative F-distances
# and collapses to a near-constant (~0.5) on many datasets, which makes
# roc_auc_score return exactly 0.5. Use decision_function distances
# instead (see simca_scores_for_auc).
y_score = simca_scores_for_auc(simca_cv, X_test, is_binary)
if y_score is not None:
try:
if is_binary:
auc_fold = roc_auc_score(y_test, y_score)
else:
auc_fold = roc_auc_score(y_test, y_score, multi_class='ovr', average='macro')
except ValueError:
# e.g. only one class present in this fold's y_test
auc_fold = None
cv_auc_scores.append(auc_fold)
# Keep the same score type in the aggregated-AUC accumulator
cv_scores_all.extend(np.asarray(y_score).tolist())
else:
cv_auc_scores.append(None)
# Store fold results
accuracy = accuracy_score(y_test, y_pred)
cv_scores.append(accuracy)
cv_predictions.extend(y_pred)
cv_true_labels.extend(y_test)
# Store per-fold classification report
fold_report = classification_report(y_test, y_pred, output_dict=True, zero_division=0)
cv_class_reports.append(fold_report)
# Create confusion matrix using overall classes to ensure consistent dimensions
cm = confusion_matrix(y_test, y_pred, labels=overall_classes_list)
cv_confusion_matrices.append(cm)
# Calculate sensitivity and specificity for binary classification using the full 2x2 matrix
if is_binary:
tn, fp, fn, tp = cm.ravel()
sensitivity = tp / (tp + fn) if (tp + fn) > 0 else 0
specificity = tn / (tn + fp) if (tn + fp) > 0 else 0
cv_sensitivity_scores.append(sensitivity)
cv_specificity_scores.append(specificity)
else:
cv_sensitivity_scores.append(None)
cv_specificity_scores.append(None)
except (ValueError, np.linalg.LinAlgError) as e:
print(f"CV fold failed: {e}")
cv_scores.append(0.0)
# Add fallback predictions
fallback_pred = [y_train[0]] * len(y_test) if len(y_train) > 0 else [overall_classes_list[0]] * len(y_test)
cv_predictions.extend(fallback_pred)
cv_true_labels.extend(y_test)
# Add fallback confusion matrix with overall classes
fallback_cm = np.zeros((n_overall_classes, n_overall_classes))
cv_confusion_matrices.append(fallback_cm)
# Fallback sensitivity/specificity
if is_binary:
cv_sensitivity_scores.append(0.0)
cv_specificity_scores.append(0.0)
else:
cv_sensitivity_scores.append(None)
cv_specificity_scores.append(None)
# AUC fallback
cv_auc_scores.append(None)
# ---------- After CV loop: aggregated AUC and summary ----------
auc_aggregated = None
if (cv_scores_all
and any(a is not None for a in cv_auc_scores)
and len(cv_scores_all) == len(cv_true_labels)
and len(np.unique(cv_true_labels)) > 1):
try:
if is_binary:
auc_aggregated = roc_auc_score(cv_true_labels, cv_scores_all)
else:
scores_matrix = np.array(cv_scores_all)
auc_aggregated = roc_auc_score(cv_true_labels, scores_matrix,
multi_class='ovr', average='macro')
except ValueError:
auc_aggregated = None
cv_auc_scores_filtered = [a for a in cv_auc_scores if a is not None]
cv_auc_mean = np.mean(cv_auc_scores_filtered) if cv_auc_scores_filtered else None
cv_auc_std = np.std(cv_auc_scores_filtered) if cv_auc_scores_filtered else None
# Calculate CV-based metrics
cv_precision_scores = []
cv_recall_scores = []
cv_f1_scores = []
for report in cv_class_reports:
if 'macro avg' in report:
cv_precision_scores.append(report['macro avg']['precision'])
cv_recall_scores.append(report['macro avg']['recall'])
cv_f1_scores.append(report['macro avg']['f1-score'])
# Calculate overall CV confusion matrix (aggregated)
cv_confusion_matrix_aggregated = confusion_matrix(cv_true_labels, cv_predictions,
labels=overall_classes_list).tolist()
if is_binary:
cm_agg = np.array(cv_confusion_matrix_aggregated)
tn, fp, fn, tp = cm_agg.ravel()
cv_sensitivity_aggregated = tp / (tp + fn) if (tp + fn) > 0 else 0
cv_specificity_aggregated = tn / (tn + fp) if (tn + fp) > 0 else 0
else:
cv_sensitivity_aggregated = None
cv_specificity_aggregated = None
if cv_confusion_matrices:
cv_confusion_matrix_avg, avg_matrix_classes = create_average_confusion_matrix(
cv_confusion_matrices, [overall_classes_list] * len(cv_confusion_matrices)
)
if is_binary:
cm_avg = np.array(cv_confusion_matrix_avg)
tn, fp, fn, tp = cm_avg.ravel()
cv_sensitivity_avg = tp / (tp + fn) if (tp + fn) > 0 else 0
cv_specificity_avg = tn / (tn + fp) if (tn + fp) > 0 else 0
else:
cv_sensitivity_avg = None
cv_specificity_avg = None
else:
cv_confusion_matrix_avg = []
avg_matrix_classes = []
cv_sensitivity_avg = None
cv_specificity_avg = None
# Fit final model on all data
simca_final_model = None # stays None if the final fit fails (no model bundle)
try:
simca = SIMCA(n_components=n_components, alpha=alpha, scale_x=scale_x)
simca.fit(X_processed, y)
simca_final_model = simca
y_pred_full = simca.predict(X_processed)
y_pred_proba = simca.predict_proba(X_processed)
final_classes = simca.classes_.tolist() if hasattr(simca, 'classes_') else overall_classes_list
except (ValueError, np.linalg.LinAlgError) as e:
print(f"Final model fitting failed: {e}")
n_classes = len(overall_classes_list)
y_pred_full = np.array([y[0]] * len(y)) if len(y) > 0 else np.array([overall_classes_list[0]] * len(y))
y_pred_proba = np.ones((len(y), n_classes)) / n_classes
final_classes = overall_classes_list
# Reorder the average confusion matrix to match final_classes if needed
if cv_confusion_matrix_avg and avg_matrix_classes and avg_matrix_classes != final_classes:
avg_to_final = {cls: idx for idx, cls in enumerate(avg_matrix_classes)}
n_classes = len(final_classes)
reordered_matrix = np.zeros((n_classes, n_classes))
for i, true_cls in enumerate(final_classes):
for j, pred_cls in enumerate(final_classes):
if true_cls in avg_to_final and pred_cls in avg_to_final:
orig_i = avg_to_final[true_cls]
orig_j = avg_to_final[pred_cls]
if (orig_i < len(cv_confusion_matrix_avg) and
orig_j < len(cv_confusion_matrix_avg[0])):
reordered_matrix[i, j] = cv_confusion_matrix_avg[orig_i][orig_j]
cv_confusion_matrix_avg = reordered_matrix.tolist()
avg_matrix_classes = final_classes
# Full model sensitivity and specificity
full_confusion_matrix = confusion_matrix(y, y_pred_full, labels=final_classes)
if is_binary and full_confusion_matrix.shape == (2, 2):
tn, fp, fn, tp = full_confusion_matrix.ravel()
full_sensitivity = tp / (tp + fn) if (tp + fn) > 0 else 0
full_specificity = tn / (tn + fp) if (tn + fp) > 0 else 0
else:
full_sensitivity = None
full_specificity = None
validation_message = None
if is_binary and full_sensitivity is not None and full_specificity is not None:
avg_sens_spec = (full_sensitivity + full_specificity) / 2
accuracy = accuracy_score(y, y_pred_full)
accuracy_diff = abs(accuracy - avg_sens_spec)
if accuracy_diff > 0.01:
validation_message = f"Note: Accuracy ({accuracy:.3f}) differs from (sensitivity + specificity)/2 ({avg_sens_spec:.3f}) by {accuracy_diff:.3f}"
# Ensure CV scores are valid
if not cv_scores or len(cv_scores) == 0:
cv_scores = [0.0] * cv_folds
# Sensitivity / specificity mean and std
cv_sensitivity_scores_filtered = [s for s in cv_sensitivity_scores if s is not None]
cv_sensitivity_mean = np.mean(cv_sensitivity_scores_filtered) if cv_sensitivity_scores_filtered else None
cv_sensitivity_std = np.std(cv_sensitivity_scores_filtered) if cv_sensitivity_scores_filtered else None
cv_specificity_scores_filtered = [s for s in cv_specificity_scores if s is not None]
cv_specificity_mean = np.mean(cv_specificity_scores_filtered) if cv_specificity_scores_filtered else None
cv_specificity_std = np.std(cv_specificity_scores_filtered) if cv_specificity_scores_filtered else None
# --- Full-data (resubstitution) AUC and accuracy, field parity with run_lda_analysis ---
# Use the continuous SIMCA decision score (decision_function), NOT predict_proba,
# which collapses to a constant and yields AUC == 0.5 (see simca_scores_for_auc).
full_auc = None
if is_binary:
try:
full_scores = simca_scores_for_auc(simca, X_processed, is_binary=True)
full_auc = roc_auc_score(y, full_scores)
except Exception:
full_auc = None
try:
full_accuracy = float(accuracy_score(y, y_pred_full))
except Exception:
full_accuracy = None
# Sanitize probabilities for JSON serialization
def sanitize_array(arr):
arr = np.array(arr)
arr = np.nan_to_num(arr, nan=0.0, posinf=1.0, neginf=0.0)
return arr.tolist()
results = {
'method': 'SIMCA',
'cv_scores': cv_scores,
'cv_mean': float(np.mean(cv_scores)) if cv_scores else 0.0,
'cv_std': float(np.std(cv_scores)) if cv_scores else 0.0,
'cv_precision': np.mean(cv_precision_scores) if cv_precision_scores else 0.0,
'cv_precision_std': np.std(cv_precision_scores) if cv_precision_scores else 0.0,
'cv_recall': np.mean(cv_recall_scores) if cv_recall_scores else 0.0,
'cv_recall_std': np.std(cv_recall_scores) if cv_recall_scores else 0.0,
'cv_f1': np.mean(cv_f1_scores) if cv_f1_scores else 0.0,
'cv_f1_std': np.std(cv_f1_scores) if cv_f1_scores else 0.0,
# Sensitivity
'cv_sensitivity': cv_sensitivity_mean,
'cv_sensitivity_std': cv_sensitivity_std,
'cv_sensitivity_aggregated': cv_sensitivity_aggregated,
'cv_sensitivity_avg': cv_sensitivity_avg,
'cv_sensitivity_scores': cv_sensitivity_scores,
# Specificity
'cv_specificity': cv_specificity_mean,
'cv_specificity_std': cv_specificity_std,
'cv_specificity_aggregated': cv_specificity_aggregated,
'cv_specificity_avg': cv_specificity_avg,
'cv_specificity_scores': cv_specificity_scores,
# AUC (new)
'cv_auc': sanitize_for_json(cv_auc_mean),
'cv_auc_std': sanitize_for_json(cv_auc_std),
'cv_auc_scores': sanitize_for_json(cv_auc_scores),
'cv_auc_aggregated': sanitize_for_json(auc_aggregated),
# Full-data AUC / accuracy (+ CV accuracy alias), field parity with LDA
'auc': sanitize_for_json(full_auc),
'accuracy': full_accuracy,
'cv_accuracy': float(np.mean(cv_scores)) if cv_scores else 0.0,
# Validation
'is_binary': is_binary,
'validation_message': validation_message,
# Confusion matrices
'cv_confusion_matrix': cv_confusion_matrix_avg,
'cv_confusion_matrices': [cm.tolist() for cm in cv_confusion_matrices],
'cv_confusion_matrix_aggregated': cv_confusion_matrix_aggregated,
'cv_class_reports': cv_class_reports,
'classification_report': classification_report(y, y_pred_full, output_dict=True, zero_division=0),
'confusion_matrix': full_confusion_matrix.tolist(),
'sensitivity': full_sensitivity,
'specificity': full_specificity,
'labels': y.tolist(),
'predictions': y_pred_full.tolist(),
'prediction_probabilities': sanitize_array(y_pred_proba),
'classes': final_classes,
'model_info': {
'n_components': n_components,
'alpha': alpha,
'scale_x': scale_x,
'n_classes': n_overall_classes
}
}
if feature_selection_params:
results['feature_selection'] = {
'selected_features': selected_features,
'feature_scores': feature_scores,
'method': feature_selection_params.get('method', 'anova'),
'k': feature_selection_params.get('k', 'all')
}
if is_binary and full_sensitivity is not None and full_specificity is not None:
results['metric_consistency_check'] = {
'accuracy': accuracy_score(y, y_pred_full),
'sensitivity': full_sensitivity,
'specificity': full_specificity,
'avg_sens_spec': (full_sensitivity + full_specificity) / 2,
'difference': abs(accuracy_score(y, y_pred_full) - (full_sensitivity + full_specificity) / 2)
}
# Final fitted model for inference on unknown samples (popped in run_analysis)
return attach_model_bundle(results, simca_final_model, selected_features, feature_selection_params)
Limitations and common pitfalls
- Nearest-class rule. A sample of an unknown kind is still assigned to the closest class. SIMCA in Chrometrica cannot say "none of the classes".
- Small classes. A class with few rows gives a PCA model with few components and a noisy residual scale.
- Number of components. Too few components leave class structure in the residuals; too many make every class model fit every sample.
References
- Wold S. Pattern recognition by means of disjoint principal components models. Pattern Recognition, 8, 127–139 (1976). doi:10.1016/0031-3203(76)90014-5
- Wold S., Sjöström M. SIMCA: a method for analyzing chemical data in terms of similarity and analogy. In: Chemometrics: Theory and Application, ACS Symposium Series 52, 243–282 (1977).
- Pomerantsev A. L. Acceptance areas for multivariate classification derived by projection methods. Journal of Chemometrics, 22, 601–609 (2008). doi:10.1002/cem.1147