Skip to main content

LDA — Linear Discriminant Analysis

TaskClassification
Method key (analysis_type)lda
Prediction on new dataYes — see Prediction
Librarysklearn.discriminant_analysis.LinearDiscriminantAnalysis

When to use​

LDA is a good first classifier:

  • the classes differ mainly in their mean colour;
  • the spread inside the classes is similar;
  • there are clearly more samples than features.

It is fast, has no parameters to tune, and gives a low-dimensional projection that you can look at on a score plot.

How it works​

LDA assumes that the samples of every class kk follow a normal distribution with their own mean μk\mu_k and a covariance matrix Σ\Sigma that is the same for all classes. A sample xx is assigned to the class with the largest discriminant score

δk(x)=x⊤Σ−1μk−12 μk⊤Σ−1μk+log⁡πk,\delta_k(x) = x^\top \Sigma^{-1} \mu_k - \tfrac{1}{2}\, \mu_k^\top \Sigma^{-1} \mu_k + \log \pi_k ,

where πk\pi_k is the share of class kk in the training data. The boundaries between classes are therefore straight lines (hyperplanes).

The same model also gives a projection: at most K−1K-1 discriminant axes (for KK classes) along which the ratio of the between-class variance to the within-class variance is the largest (Fisher's criterion). These axes are the "components" on the score plot.

Parameters​

Defaults as set by the analysis dialog. Everything else uses the scikit-learn defaults: no shrinkage, class priors taken from the class sizes.

ParameterUI labelDefaultNotes
solverSolversvdsvd: singular value decomposition, does not compute the covariance matrix and works when there are more features than samples. eigen: eigenvalue decomposition of the covariance matrices; can fail when the within-class covariance is singular (more features than samples). Any other value is replaced by svd.
feature_selectionEnable Feature SelectionoffSee Feature selection.
average_replicatesAverage technical replicatesoffSee Averaging replicates. Not available with the sample_name label column.

Preprocessing​

No scaling is applied. LDA does not need it: the result does not change when you multiply a feature by a constant.

Before the model: replicate averaging and feature selection, if they are on. The order of steps is described in Cross-validation → Order of steps.

Results and metrics​

The common classification metrics (CV accuracy, precision, recall, F1, sensitivity, specificity, AUC, confusion matrices) are described in Metrics. For binary problems, AUC uses the LDA decision function as the score.

LDA-specific fields:

FieldMeaning
transformed_dataCoordinates of every sample on the discriminant axes (final model)
loadingsscalings_: the weights of the features on each discriminant axis
explained_variance_ratioShare of the between-class variance explained by each axis
model_infoThe solver actually used

Visualizations​

Source code​

run_lda_analysis in chrometrica/analysis/analysis.py. The imports, the model in the CV loop and the final model:

analysis.py · run_lda_analysis() · lines 976–1178
def run_lda_analysis(X, y, parameters, cv_method, cv_folds):
    """Run LDA analysis with cross-validation and optional feature selection"""
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
    from sklearn.model_selection import StratifiedKFold, GroupKFold
    from sklearn.metrics import classification_report, confusion_matrix, accuracy_score, roc_auc_score
    import numpy as np
    # ...
    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]

        # Train and predict
        lda_cv = LinearDiscriminantAnalysis(solver=solver)
        lda_cv.fit(X_train, y_train)
        y_pred = lda_cv.predict(X_test)
    # ...
    # Fit final model on all data
    lda = LinearDiscriminantAnalysis(solver=solver)
    lda.fit(X_processed, y)
    X_lda = lda.transform(X_processed)
    y_pred_full = lda.predict(X_processed)
Full source of run_lda_analysis()
analysis.py · run_lda_analysis() · lines 976–1333
def run_lda_analysis(X, y, parameters, cv_method, cv_folds):
    """Run LDA analysis with cross-validation and optional feature selection"""
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
    from sklearn.model_selection import StratifiedKFold, GroupKFold
    from sklearn.metrics import classification_report, confusion_matrix, accuracy_score, roc_auc_score
    import numpy as np

    # ---------- Helper to sanitize JSON ----------
    def sanitize_for_json(obj):
        if isinstance(obj, dict):
            return {k: sanitize_for_json(v) for k, v in obj.items()}
        elif isinstance(obj, list):
            return [sanitize_for_json(v) for v in obj]
        elif isinstance(obj, float) and np.isnan(obj):
            return None
        elif isinstance(obj, (np.float32, np.float64)) and np.isnan(obj):
            return None
        elif isinstance(obj, np.ndarray):
            return sanitize_for_json(obj.tolist())
        else:
            return obj

    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

    # 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

    # 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)

    # Solver from the analysis dialog. 'lsqr' is not offered: it has no
    # transform(), which the score/loading plots need.
    solver = parameters.get('solver', 'svd')
    if solver not in ('svd', 'eigen'):
        logger.warning(f"LDA: unsupported solver '{solver}', using 'svd'")
        solver = 'svd'

    # --- AUC initialisation (only used if binary) ---
    cv_auc_scores = []       # AUC per fold
    cv_scores_all = []       # decision scores / probabilities for all test samples

    # Enhanced cross-validation with metrics
    cv_scores = []
    cv_predictions = []
    cv_true_labels = []
    cv_class_reports = []
    cv_confusion_matrices = []
    cv_sensitivity_scores = []
    cv_specificity_scores = []

    if cv_method == 'stratified':
        cv = StratifiedKFold(n_splits=cv_folds, shuffle=True, random_state=42)
    elif cv_method == 'group':
        groups = parameters.get('groups')
        if groups is None:
            raise ValueError("Groups parameter required for group cross-validation")
        cv = GroupKFold(n_splits=cv_folds)
    else:
        from sklearn.model_selection import KFold
        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]

        # Train and predict
        lda_cv = LinearDiscriminantAnalysis(solver=solver)
        lda_cv.fit(X_train, y_train)
        y_pred = lda_cv.predict(X_test)

        # --- AUC computation only for binary ---
        if is_binary:
            # Get scores (decision function or probability of positive class)
            if hasattr(lda_cv, "decision_function"):
                y_score = lda_cv.decision_function(X_test)
            elif hasattr(lda_cv, "predict_proba"):
                y_prob = lda_cv.predict_proba(X_test)
                y_score = y_prob[:, 1]   # positive class probability
            else:
                y_score = None

            if y_score is not None:
                # Attempt to compute AUC for this fold; if it fails (e.g. only one class), set to None
                try:
                    auc_fold = roc_auc_score(y_test, y_score)
                except ValueError:
                    auc_fold = None
                cv_auc_scores.append(auc_fold)
                # Store scores for aggregated AUC (always store, even if auc_fold is None)
                if isinstance(y_score, np.ndarray):
                    cv_scores_all.extend(y_score.tolist())
                else:
                    cv_scores_all.extend(y_score)
            else:
                cv_auc_scores.append(None)
        else:
            # Multiclass: skip AUC entirely
            cv_auc_scores.append(None)

        # Store fold results (unchanged)
        cv_scores.append(accuracy_score(y_test, y_pred))
        cv_predictions.extend(y_pred)
        cv_true_labels.extend(y_test)

        fold_report = classification_report(y_test, y_pred, output_dict=True, zero_division=0)
        cv_class_reports.append(fold_report)

        cm = confusion_matrix(y_test, y_pred, labels=overall_classes_list)
        cv_confusion_matrices.append(cm)

        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)

    # --- After CV loop: compute aggregated AUC only for binary ---
    auc_aggregated = None
    if is_binary and cv_scores_all:
        # Try to compute aggregated AUC from all scores; if it fails, set to None
        try:
            auc_aggregated = roc_auc_score(cv_true_labels, cv_scores_all)
        except ValueError:
            auc_aggregated = None
        # Ensure no NaN
        if auc_aggregated is not None and np.isnan(auc_aggregated):
            auc_aggregated = None

    # Compute mean/std of per-fold AUC (only binary)
    if is_binary:
        valid_auc = [a for a in cv_auc_scores if a is not None]
        cv_auc_mean = np.mean(valid_auc) if valid_auc else None
        cv_auc_std = np.std(valid_auc) if valid_auc else None
    else:
        cv_auc_mean = None
        cv_auc_std = None

    # ---- Continue with existing CV metrics (unchanged) ----
    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'])

    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
    lda = LinearDiscriminantAnalysis(solver=solver)
    lda.fit(X_processed, y)
    X_lda = lda.transform(X_processed)
    y_pred_full = lda.predict(X_processed)

    final_classes = lda.classes_.tolist()

    # Reorder average confusion matrix if needed
    if cv_confusion_matrix_avg 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

    # Explained variance ratio
    if hasattr(lda, 'explained_variance_ratio_'):
        explained_variance_ratio = lda.explained_variance_ratio_.tolist()
    else:
        try:
            if hasattr(lda, 'scalings_'):
                X_projected = lda.transform(X_processed)
                total_var = np.var(X_projected, axis=0, ddof=1)
                explained_variance_ratio = (total_var / np.sum(total_var)).tolist()
            else:
                n_components = X_lda.shape[1] if len(X_lda.shape) > 1 else 1
                explained_variance_ratio = [1.0 / n_components] * n_components
        except:
            n_components = X_lda.shape[1] if len(X_lda.shape) > 1 else 1
            explained_variance_ratio = [1.0 / n_components] * n_components

    cumulative_variance = np.cumsum(explained_variance_ratio).tolist()
    cv_confusion_matrices_lists = [cm.tolist() for cm in cv_confusion_matrices]
    loadings = lda.scalings_.tolist() if hasattr(lda, 'scalings_') else None

    # 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}"

    accuracy = accuracy_score(y, y_pred_full)

    auc = None
    if is_binary:
        try:
            if hasattr(lda, "predict_proba"):
                y_prob = lda.predict_proba(X_processed)
                # positive class is overall_classes[1]
                pos_idx = overall_classes_list.index(overall_classes[1])
                pos_proba = y_prob[:, pos_idx]
                auc = roc_auc_score(y, pos_proba)
        except Exception:
            auc = None

    # Mean and std for sensitivity/specificity
    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

    # Build results dictionary – AUC fields are sanitized
    results = {
        'method': 'LDA',
        'cv_scores': cv_scores,
        'cv_mean': np.mean(cv_scores) if cv_scores else 0.0,
        'cv_accuracy': np.mean(cv_scores) if cv_scores else 0.0,
        'cv_std': 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 metrics
        '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 metrics
        '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 metrics (sanitized; for multiclass they become None/empty)
        '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),
        # Validation flag
        'is_binary': is_binary,
        'validation_message': validation_message,
        # Existing metrics
        'cv_confusion_matrix': cv_confusion_matrix_avg,
        'cv_confusion_matrices': cv_confusion_matrices_lists,
        '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,
        'auc': auc,
        'accuracy': accuracy,
        'transformed_data': X_lda.tolist(),
        'loadings': loadings,
        'labels': y.tolist(),
        'predictions': y_pred_full.tolist(),
        'classes': final_classes,
        'explained_variance_ratio': explained_variance_ratio,
        'cumulative_variance': cumulative_variance,
        'model_info': {'solver': solver}
    }

    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)
        }

    results = sanitize_for_json(results)

    # Final fitted model for inference on unknown samples (popped in run_analysis)
    return attach_model_bundle(results, lda, selected_features, feature_selection_params)

Limitations and common pitfalls​

  • More features than samples. With few samples per class and many features, the covariance estimate is poor and LDA overfits. Use feature selection, or PLS-DA, which is built for this case.
  • At least two rows per class. A class with one row makes the within-class covariance undefined.
  • Linear boundaries only. If classes are separated by a curved boundary, try SVM or random forest.
  • Different spreads. If one class is much more variable than another, the shared-covariance assumption fails and accuracy drops.

References​