Permutation test

A permutation test answers one question: could a model this good have been obtained by chance, from data in which the features say nothing about the class? It is the most direct check that a cross-validated result is real, and it is especially useful for small data sets, where a high CV accuracy can appear by luck.

In Chrometrica, the test runs when you build a permutation plot for a finished analysis (Permutation Plot (Validation) in the visualization dialog).

Available for: LDA, logistic regression, SVM, random forest, k-NN, PLS-DA, SIMCA. Binary classification only: the label column must have exactly two classes.

How the test works

1. Take the data of the analysis: the same rows, features, CV scheme and folds, replicate averaging and model parameters.
2. Compute the chosen metric for the original labels, twice: on the training data (the model fitted on all rows and tested on the same rows) and with cross-validation.
3. Shuffle the labels randomly, keeping the features in place. This destroys any relation between features and class while keeping everything else (class sizes, feature distributions) the same.
4. Repeat step 2 for the shuffled labels. Both metrics are computed against the shuffled labels: the labels the model was fitted on.
5. Repeat steps 3–4 Number of Permutations times.
6. Compare the original metric with the distribution of the shuffled ones and compute a p-value, separately for the training and the CV metric.

A shuffled model can still reach a high training metric by memorising random labels, but it should not reach a high CV metric. If the original CV metric is well above the shuffled ones, the model has learned something real.

Settings

UI label	Parameter	Default	Notes
Metric	metric	PLS-DA: R² / Q²; others: AUC	PLS-DA: R² / Q², Accuracy, AUC. Other methods: AUC, Accuracy.
Number of Permutations	n_permutations	200	50–500 in the dialog
Feature Selection in Permutations	feature_selection_mode	Repeated for each permutation	Only if the analysis used feature selection. See below.

The shuffles use a fixed seed (random_state = 42), so the same settings give the same result.

Metrics:

• Accuracy: share of correctly classified rows.
• AUC: area under the ROC curve, computed on the scores of all test rows pooled together. For SIMCA, the distance-based score is used (see SIMCA).
• R² / Q² (PLS-DA only): the classes are coded 0 and 1 and compared with the predicted class codes; R² on the training data, Q² with CV. See PLS-DA → Results.

Feature selection:

Option	feature_selection_mode	What happens
Repeated for each permutation (recommended)	per_permutation	Features are re-selected on the shuffled labels in every permutation, exactly as in the original analysis. The original and the shuffled models go through the same procedure: the p-value is fair.
Once, on original labels (as in analysis)	once	Features are selected once, on the original labels, and every shuffled model uses them. Faster, but the selected features "know" the original labels: the original model looks better than it is, and the p-value is too small.

Number of permutations

The smallest p-value the test can give is $1/(m+1)$, where $m$ is the number of permutations: 0.005 for 200 permutations, 0.01 for 100, 0.02 for 50. Use at least 100 permutations if you want to report $p < 0.01$.

The test fits $m + 1$ times as many models as the analysis itself (one CV run and one full fit per permutation). With leave-one-sample-out Group K-Fold on a large data set this can take a long time; start with 50–100 permutations.

Reading the p-value

For the original metric $s_0$ and the metrics $s_1, \dots, s_m$ of the permutations that could be computed,

$$p = \frac{\#\{\, i : s_i \ge s_0 \,\} + 1}{m + 1} .$$

This is the standard estimate (Phipson & Smyth, 2010), also used by sklearn.model_selection.permutation_test_score:

• the original labelling counts as one of the permutations, so $p$ is never 0;
• ties count against the model. Metrics such as accuracy take few distinct values and often reach 1.0; a shuffled model with the same value as the original one increases $p$.

Two p-values are reported:

• CV p-value: the one to use. A small value (for example, $p \le 0.05$) means that the cross-validated performance is unlikely to be due to chance.
• Train p-value: for comparison only. Flexible models (random forest, k-NN) often fit shuffled labels perfectly on the training data, so this value is often large even for a good model.

On the plot card, both p-values are shown in green when they are $\le 0.05$ and in red otherwise. With the default 200 permutations the smallest possible value is 0.005.

What the p-value does and does not tell you is explained in Interpreting permutation p-values.

Permutation plot

The plot shows each permutation as a point:

• x axis: the correlation between the shuffled and the original labels (1.0 for the original labels, around 0 for a thorough shuffle);
• y axis: the metric value, for the training data and for CV.

The original values are marked with stars at x = 1. A straight line is fitted through the training points and through the CV points; the legend shows their intercepts (the classic validation plot for PLS-DA). A good model has its original CV point clearly above the cloud of shuffled CV points. See Permutation plot for the plot options.

Source code

run_permutation_test in chrometrica/analysis/analysis.py. It is called by generate_permutation_plot in chrometrica/analysis/tasks.py, which rebuilds the data of the analysis and reads the model parameters from results['model_info'].

The estimator used for each analysis type:

analysis.py · run_permutation_test() · lines 3810–3867

def _build_estimator():
    """Build an estimator that mirrors what the corresponding
    run_*_analysis actually fits (see work/step1_plan.txt, 4.B / B0).

    Hyper-parameters are read from ``parameters``, populated from
    ``results['model_info']`` by generate_permutation_plot. Missing keys
    (analyses saved before the parameter was recorded) fall back to the
    defaults of the corresponding run_*_analysis.
    """
    at = analysis_type
    if at == 'logistic':
        from sklearn.linear_model import LogisticRegression
        est = LogisticRegression(
            random_state=42,
            max_iter=int(parameters.get('max_iter', 1000)),
            C=float(parameters.get('C', 1.0)),
        )
        if bool(parameters.get('scale_data', True)):
            from sklearn.preprocessing import StandardScaler
            from sklearn.pipeline import Pipeline
            return Pipeline([('scaler', StandardScaler()), ('logistic', est)])
        return est
    if at == 'svm':
        from sklearn.svm import SVC
        return SVC(
            random_state=42, probability=True,
            kernel=parameters.get('kernel', 'rbf'),
            C=float(parameters.get('C', 1.0)),
        )
    if at == 'random_forest':
        from sklearn.ensemble import RandomForestClassifier
        return RandomForestClassifier(
            random_state=42,
            n_estimators=int(parameters.get('n_estimators', 100)),
            max_depth=_optional_int(parameters.get('max_depth')),
        )
    if at == 'knn':
        from sklearn.neighbors import KNeighborsClassifier
        est = KNeighborsClassifier(
            n_neighbors=int(parameters.get('n_neighbors', 5)),
            weights=parameters.get('weights', 'uniform'),
            algorithm=parameters.get('algorithm', 'auto'),
            metric=parameters.get('knn_metric', 'minkowski'),
            p=int(parameters.get('p', 2)),
        )
        if bool(parameters.get('scale_data', True)):
            from sklearn.preprocessing import StandardScaler
            from sklearn.pipeline import Pipeline
            return Pipeline([('scaler', StandardScaler()), ('knn', est)])
        return est
    if at == 'simca':
        return SIMCA(
            n_components=int(parameters.get('n_components', 2)),
            alpha=float(parameters.get('alpha', 0.05)),
            scale_x=bool(parameters.get('scale_x', True)),
        )
    # default / 'lda'
    return LinearDiscriminantAnalysis(solver=parameters.get('solver', 'svd'))

The permutation loop and the p-value:

analysis.py · run_permutation_test() · lines 4077–4117

# ---- Original and permutations ----
orig_train, orig_cv = compute_metrics(y)
correlations = [1.0]
train_metrics = [orig_train]
cv_metrics = [orig_cv]

for _ in range(n_permutations):
    y_perm = rng.permutation(y)
    y_perm_num = np.array([class_to_num[val] for val in y_perm])
    corr, _ = pearsonr(y_num, y_perm_num)
    correlations.append(corr)

    train_perm, cv_perm = compute_metrics(y_perm)
    train_metrics.append(train_perm)
    cv_metrics.append(cv_perm)

# ---- Permutation p-values ----
# p = (b + 1) / (m + 1), b = number of permutations whose metric is greater
# than OR EQUAL to the original value, m = number of valid (non-None)
# permutation values (Phipson & Smyth, 2010; as in
# sklearn.permutation_test_score). The original labelling counts as one of
# the permutations, so p is never 0; ties count against the model (metrics
# are discrete and often saturate at 1.0). Computed independently for the
# train-data run and the cross-validation run; None if it cannot be computed.
def _perm_p_value(orig, perm_values):
    if orig is None:
        return None
    valid = [v for v in perm_values if v is not None]
    if not valid:
        return None
    return (sum(1 for v in valid if v >= orig) + 1) / (len(valid) + 1)

train_p_value = _perm_p_value(orig_train, train_metrics[1:])
cv_p_value = _perm_p_value(orig_cv, cv_metrics[1:])

return {
    'correlations': correlations,
    'train_metric_values': train_metrics,
    'cv_metric_values': cv_metrics,
    'train_p_value': train_p_value,
    'cv_p_value': cv_p_value,

Full source of run_permutation_test()

analysis.py · run_permutation_test() · lines 3769–4118

def run_permutation_test(X, y, parameters, cv_method, cv_folds,
                         n_permutations=200, random_state=None,
                         analysis_type='plsda', metric='r2',
                         feature_selection_params=None):
    """
    Perform permutation testing for binary classification models.

    Both the train (self-test) and the CV metric of every point are computed
    against the labels the model was fitted on (the permuted labels).

    Parameters
    ----------
    analysis_type : str
        'plsda', 'lda', 'logistic', 'svm', 'random_forest', 'knn', 'simca'.
        Every non-PLS-DA type is handled by the generic branch which builds
        the estimator via ``_build_estimator`` (mirrors run_*_analysis).
    metric : str
        For PLS-DA: 'r2', 'accuracy', 'auc'
        For every other analysis_type: 'accuracy', 'auc' ('r2' is rejected).
    feature_selection_params : dict or None
        None: X is used as is (features already selected by the caller, if any).
        dict: feature selection (``select_features``) is re-run on the full X
        for the original labels and for every permutation, before the CV loop.
    """
    from scipy.stats import pearsonr
    from sklearn.model_selection import StratifiedKFold, GroupKFold, KFold
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
    from sklearn.metrics import accuracy_score, roc_auc_score
    import numpy as np

    overall_classes = np.unique(y)
    if len(overall_classes) != 2:
        raise ValueError("Permutation plot is only defined for binary classification.")

    overall_classes_list = overall_classes.tolist()
    class_to_num = {overall_classes_list[0]: 0, overall_classes_list[1]: 1}
    y_num = np.array([class_to_num[val] for val in y])

    groups = parameters.get('groups', None) if cv_method == 'group' else None
    rng = np.random.RandomState(random_state)

    def _build_estimator():
        """Build an estimator that mirrors what the corresponding
        run_*_analysis actually fits (see work/step1_plan.txt, 4.B / B0).

        Hyper-parameters are read from ``parameters``, populated from
        ``results['model_info']`` by generate_permutation_plot. Missing keys
        (analyses saved before the parameter was recorded) fall back to the
        defaults of the corresponding run_*_analysis.
        """
        at = analysis_type
        if at == 'logistic':
            from sklearn.linear_model import LogisticRegression
            est = LogisticRegression(
                random_state=42,
                max_iter=int(parameters.get('max_iter', 1000)),
                C=float(parameters.get('C', 1.0)),
            )
            if bool(parameters.get('scale_data', True)):
                from sklearn.preprocessing import StandardScaler
                from sklearn.pipeline import Pipeline
                return Pipeline([('scaler', StandardScaler()), ('logistic', est)])
            return est
        if at == 'svm':
            from sklearn.svm import SVC
            return SVC(
                random_state=42, probability=True,
                kernel=parameters.get('kernel', 'rbf'),
                C=float(parameters.get('C', 1.0)),
            )
        if at == 'random_forest':
            from sklearn.ensemble import RandomForestClassifier
            return RandomForestClassifier(
                random_state=42,
                n_estimators=int(parameters.get('n_estimators', 100)),
                max_depth=_optional_int(parameters.get('max_depth')),
            )
        if at == 'knn':
            from sklearn.neighbors import KNeighborsClassifier
            est = KNeighborsClassifier(
                n_neighbors=int(parameters.get('n_neighbors', 5)),
                weights=parameters.get('weights', 'uniform'),
                algorithm=parameters.get('algorithm', 'auto'),
                metric=parameters.get('knn_metric', 'minkowski'),
                p=int(parameters.get('p', 2)),
            )
            if bool(parameters.get('scale_data', True)):
                from sklearn.preprocessing import StandardScaler
                from sklearn.pipeline import Pipeline
                return Pipeline([('scaler', StandardScaler()), ('knn', est)])
            return est
        if at == 'simca':
            return SIMCA(
                n_components=int(parameters.get('n_components', 2)),
                alpha=float(parameters.get('alpha', 0.05)),
                scale_x=bool(parameters.get('scale_x', True)),
            )
        # default / 'lda'
        return LinearDiscriminantAnalysis(solver=parameters.get('solver', 'svd'))

    def compute_metrics(y_cur):
        """Return (train_metric, cv_metric) for the labels y_cur."""
        if feature_selection_params:
            X_cur, _, _ = select_features(X, y_cur, feature_selection_params, 'classification')
        else:
            X_cur = X

        cv_splitter = None
        splitter_groups = None
        if cv_method == 'group':
            if groups is None:
                raise ValueError("Groups required for group CV.")
            cv_splitter = GroupKFold(n_splits=cv_folds)
            splitter_groups = groups
        elif cv_method == 'stratified':
            cv_splitter = StratifiedKFold(n_splits=cv_folds, shuffle=True, random_state=random_state)
        else:
            cv_splitter = KFold(n_splits=cv_folds, shuffle=True, random_state=random_state)

        cv_predictions = []
        cv_true_labels = []
        cv_scores_all = []  # for AUC

        if analysis_type == 'plsda':
            style = parameters.get('style', 'hard')
            alpha = parameters.get('alpha', 0.05)
            gamma = parameters.get('gamma', 0.01)
            scale_x = bool(parameters.get('scale_x', True))
            n_components = parameters.get('n_components', 7)

            for train_idx, test_idx in cv_splitter.split(X_cur, y_cur, groups=splitter_groups):
                X_train, X_test = X_cur[train_idx], X_cur[test_idx]
                y_train, y_test = y_cur[train_idx], y_cur[test_idx]

                model = PLSDA(
                    n_components=n_components,
                    style=style,
                    alpha=alpha,
                    gamma=gamma,
                    scale_x=scale_x,
                    random_state=random_state
                )
                model.fit(X_train, y_train)

                if style == 'hard':
                    y_pred = model.predict(X_test)
                else:
                    y_pred_list = model.predict(X_test)
                    y_pred = []
                    for pred in y_pred_list:
                        if pred and pred[0] != "NOT_ASSIGNED":
                            y_pred.append(pred[0])
                        else:
                            y_pred.append("UNKNOWN")
                cv_predictions.extend(y_pred)
                cv_true_labels.extend(y_test)

                # Collect probabilities for AUC if needed
                if metric == 'auc':
                    if hasattr(model, "predict_proba"):
                        y_prob = model.predict_proba(X_test)
                        pos_idx = overall_classes_list.index(overall_classes[1])
                        cv_scores_all.extend(y_prob[:, pos_idx].tolist())

            # --- Compute requested metric ---
            y_true_num = np.array([class_to_num.get(label, np.nan) for label in cv_true_labels])
            y_pred_num = np.array([class_to_num.get(pred, 0.5) for pred in cv_predictions])
            mask = ~np.isnan(y_true_num) & ~np.isnan(y_pred_num)
            y_true_num = y_true_num[mask]
            y_pred_num = y_pred_num[mask]

            if metric == 'r2':
                # Q² (cross-validated R²)
                cv_metric = None
                if len(y_true_num) > 0:
                    mean_y = np.mean(y_true_num)
                    ss_tot = np.sum((y_true_num - mean_y) ** 2)
                    if ss_tot > 0:
                        ss_res = np.sum((y_true_num - y_pred_num) ** 2)
                        cv_metric = 1 - ss_res / ss_tot
            elif metric == 'accuracy':
                cv_metric = accuracy_score(cv_true_labels, cv_predictions) if cv_predictions else None
            elif metric == 'auc':
                cv_metric = None
                if cv_scores_all and len(np.unique(cv_true_labels)) > 1:
                    try:
                        cv_metric = roc_auc_score(cv_true_labels, cv_scores_all)
                    except ValueError:
                        cv_metric = None
            else:
                raise ValueError(f"Unsupported metric: {metric}")

            # --- Train metric on full data ---
            model_full = PLSDA(
                n_components=n_components,
                style=style,
                alpha=alpha,
                gamma=gamma,
                scale_x=scale_x,
                random_state=random_state
            )
            model_full.fit(X_cur, y_cur)
            y_pred_full = model_full.predict(X_cur)
            if style == 'hard':
                y_pred_full_array = np.array(y_pred_full)
            else:
                y_pred_full_array = np.array([
                    pred[0] if (pred and pred[0] != "NOT_ASSIGNED") else "UNKNOWN"
                    for pred in y_pred_full
                ])

            if metric == 'r2':
                y_true_full = np.array([class_to_num.get(label, np.nan) for label in y_cur])
                y_pred_full_num = np.array([class_to_num.get(pred, 0.5) for pred in y_pred_full_array])
                mask_full = ~np.isnan(y_true_full) & ~np.isnan(y_pred_full_num)
                y_true_full = y_true_full[mask_full]
                y_pred_full_num = y_pred_full_num[mask_full]
                train_metric = None
                if len(y_true_full) > 0:
                    mean_y_full = np.mean(y_true_full)
                    ss_tot_full = np.sum((y_true_full - mean_y_full) ** 2)
                    if ss_tot_full > 0:
                        ss_res_full = np.sum((y_true_full - y_pred_full_num) ** 2)
                        train_metric = 1 - ss_res_full / ss_tot_full
            elif metric == 'accuracy':
                train_metric = accuracy_score(y_cur, y_pred_full_array)
            elif metric == 'auc':
                if hasattr(model_full, "predict_proba"):
                    y_prob_full = model_full.predict_proba(X_cur)
                    pos_idx = overall_classes_list.index(overall_classes[1])
                    pos_proba = y_prob_full[:, pos_idx]
                    train_metric = roc_auc_score(y_cur, pos_proba) if len(np.unique(y_cur)) > 1 else None
                else:
                    train_metric = None
            else:
                train_metric = None

        else:
            # Generic branch: lda / logistic / svm / random_forest / knn / simca.
            # The estimator is produced by _build_estimator() and mirrors the
            # model that the corresponding run_*_analysis actually fits.
            if metric not in ('accuracy', 'auc'):
                raise ValueError(
                    f"analysis_type '{analysis_type}' supports only 'accuracy' or 'auc' metric"
                )

            pos_idx = overall_classes_list.index(overall_classes[1])

            for train_idx, test_idx in cv_splitter.split(X_cur, y_cur, groups=splitter_groups):
                X_train, X_test = X_cur[train_idx], X_cur[test_idx]
                y_train, y_test = y_cur[train_idx], y_cur[test_idx]

                try:
                    model = _build_estimator()
                    model.fit(X_train, y_train)
                    y_pred = model.predict(X_test)
                except Exception:
                    # e.g. SIMCA singular covariance on a permuted fold -> skip fold
                    continue

                cv_predictions.extend(y_pred)
                cv_true_labels.extend(y_test)

                if metric == 'auc':
                    try:
                        if analysis_type == 'simca':
                            # predict_proba collapses to ~0.5; use decision_function score
                            s = np.asarray(simca_scores_for_auc(model, X_test, True))
                            cv_scores_all.extend(s.tolist())
                        elif hasattr(model, "predict_proba"):
                            y_prob = np.asarray(model.predict_proba(X_test))
                            cv_scores_all.extend(y_prob[:, pos_idx].tolist())
                    except Exception:
                        pass

            # --- CV metric ---
            if metric == 'accuracy':
                cv_metric = accuracy_score(cv_true_labels, cv_predictions) if cv_predictions else None
            else:  # auc
                cv_metric = None
                if (cv_scores_all
                        and len(cv_scores_all) == len(cv_true_labels)
                        and len(np.unique(cv_true_labels)) > 1):
                    try:
                        cv_metric = roc_auc_score(cv_true_labels, cv_scores_all)
                    except ValueError:
                        cv_metric = None

            # --- Train metric on full data ---
            train_metric = None
            try:
                model_full = _build_estimator()
                model_full.fit(X_cur, y_cur)
                if metric == 'accuracy':
                    y_pred_full = model_full.predict(X_cur)
                    train_metric = accuracy_score(y_cur, y_pred_full)
                else:  # auc
                    if analysis_type == 'simca' and len(np.unique(y_cur)) > 1:
                        s_full = simca_scores_for_auc(model_full, X_cur, True)
                        train_metric = roc_auc_score(y_cur, s_full)
                    elif hasattr(model_full, "predict_proba") and len(np.unique(y_cur)) > 1:
                        y_prob_full = np.asarray(model_full.predict_proba(X_cur))
                        train_metric = roc_auc_score(y_cur, y_prob_full[:, pos_idx])
            except Exception:
                train_metric = None

        return train_metric, cv_metric

    # ---- Original and permutations ----
    orig_train, orig_cv = compute_metrics(y)
    correlations = [1.0]
    train_metrics = [orig_train]
    cv_metrics = [orig_cv]

    for _ in range(n_permutations):
        y_perm = rng.permutation(y)
        y_perm_num = np.array([class_to_num[val] for val in y_perm])
        corr, _ = pearsonr(y_num, y_perm_num)
        correlations.append(corr)

        train_perm, cv_perm = compute_metrics(y_perm)
        train_metrics.append(train_perm)
        cv_metrics.append(cv_perm)

    # ---- Permutation p-values ----
    # p = (b + 1) / (m + 1), b = number of permutations whose metric is greater
    # than OR EQUAL to the original value, m = number of valid (non-None)
    # permutation values (Phipson & Smyth, 2010; as in
    # sklearn.permutation_test_score). The original labelling counts as one of
    # the permutations, so p is never 0; ties count against the model (metrics
    # are discrete and often saturate at 1.0). Computed independently for the
    # train-data run and the cross-validation run; None if it cannot be computed.
    def _perm_p_value(orig, perm_values):
        if orig is None:
            return None
        valid = [v for v in perm_values if v is not None]
        if not valid:
            return None
        return (sum(1 for v in valid if v >= orig) + 1) / (len(valid) + 1)

    train_p_value = _perm_p_value(orig_train, train_metrics[1:])
    cv_p_value = _perm_p_value(orig_cv, cv_metrics[1:])

    return {
        'correlations': correlations,
        'train_metric_values': train_metrics,
        'cv_metric_values': cv_metrics,
        'train_p_value': train_p_value,
        'cv_p_value': cv_p_value,
    }

References

• Phipson B., Smyth G. K. Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology, 9(1), Article 39 (2010). doi:10.2202/1544-6115.1585
• Ojala M., Garriga G. C. Permutation tests for studying classifier performance. Journal of Machine Learning Research, 11, 1833–1863 (2010). PDF
• Westerhuis J. A. et al. Assessment of PLSDA cross validation. Metabolomics, 4, 81–89 (2008). doi:10.1007/s11306-007-0099-6