Permutation nulls: is there anything here at all?#

At the sample sizes a sensory-to-chemistry study runs at, twenty or thirty products against fifty compounds, “is there anything here at all” is the hard question, and the obvious ways of answering it do not work.

Two ways that do not work#

A high in-sample R². With thirty products and three components, a good R² is close to guaranteed. It measures capacity, not evidence.

A count of variables exceeding VIP 1. This one is worth spelling out, because it looks like a statistic and is not. VIP is normalised so that

\[\sum_{j=1}^{K} \text{VIP}_j^2 = K\]

exactly, for any model, on any data:

>>> from process_improve.multivariate import PLS, vip
>>> model = PLS(n_components=2, scale=False).fit(x_scaled, y_scaled)
>>> float((vip(model) ** 2).sum()), x_scaled.shape[1]
(14.000000000000002, 14)

That identity is exactly what makes the familiar “VIP > 1” rule a sensible relative cut-off: the mean square is 1 by construction, so a score above 1 means above-average within this model. It also means the exceedance count describes the shape of the VIP distribution rather than the presence of a relationship. Permute the response and the count barely moves, so a null built on it has essentially no power and will report a false-discovery rate near 100% on data that does contain signal.

What does work#

Ask a permutation what it can achieve, and compare it with out-of-sample performance.

from process_improve.multivariate import check_predictive_signal

check_predictive_signal(chem, sensory_means, n_perm=999)

  attribute  q2_observed  q2_null_mean  q2_null_p95  p_value  n_permutations
0    fruity        0.612        -0.174        0.201    0.001             999
1     green        0.088        -0.166        0.213    0.284             999

fruity predicts held-out products far better than any reshuffling manages; green sits inside its own null. That is a distinction the VIP count above cannot draw.

Note the blocks go in unscaled. The default cross-validates PLS for you and re-derives the centring and scaling constants inside every fold, which is the second of the two traps below; scaling first would hand each fold constants computed from the row it is meant to be predicting.

This is expensive, and irreducibly so: every permutation refits once per fold, so leave-one-out on twenty products with n_perm=999 is twenty thousand fits. Develop at n_perm=50 and raise it for the number you intend to report.

Pass your own fit_predict to use a different model or a cheaper fold scheme:

def fit_predict(x, y):
    """Return out-of-sample predictions, one row per product."""
    ...  # your own CV loop; see the two responsibilities below

check_predictive_signal(chem, sensory_means, fit_predict, n_perm=999)

Two things check_predictive_signal() then leaves to you, inside fit_predict:

The cross-validation scheme is yours to choose. Leave-one-out is right when every product is precious, and needlessly expensive at hundreds of rows, where a permutation null multiplies the cost by n_perm.
Re-derive the response constants inside each fold too. It is easy to nest the predictor preprocessing and forget the response. Centre and scale the response on the training rows, predict, then back-transform with those same training constants before returning; otherwise the score is computed on a scale the fold never saw.

Whole rows are permuted, not columns. Permuting column by column would destroy the correlation structure among the attributes and inflate the null, making the test look more impressive than it is.

The p-value has a floor#

The p-value uses the (1 + count) / (n_perm + 1) form, which counts the observed statistic as one of its own null draws so it can never be exactly zero. The consequence is worth knowing before you choose n_perm: the smallest attainable p-value is 1 / (n_perm + 1). The default 500 permutations cannot report anything below 0.002, and a multiplicity correction over twenty attributes needs a floor well below the corrected threshold.

Testing the whole procedure, not one model#

count_discoveries_under_null() takes a callable that runs filtering, transformation, scaling and selection end to end, and counts discoveries under a permuted response:

>>> result = count_discoveries_under_null(select, chem, sensory_means, n_perm=200)
>>> result["observed"], result["null_mean"], result["null_to_observed_ratio"]
(9, 1.8, 0.2)

null_to_observed_ratio is not clipped. A value above 1 says shuffling found more than the real response did, which is the strongest evidence the procedure has nothing, and it is precisely the case worth seeing rather than rounding to a tidy-looking rate.

That ratio is a property of the procedure as run, which is the only version worth quoting: a procedure whose selection step is honest but whose filtering step peeked at the response has an FDR no formula recovers.

The response-independent steps are deliberately not hoisted out of the loop. Re-running them per permutation is a no-op when they really are response-independent, and hoisting them would be an assumption about your code. What is checked is determinism: a selector that disagrees with itself on two identical calls makes the FDR meaningless, and draws a SpecificationWarning rather than a silently contaminated number.

Recovering the expected class#

class_enrichment() asks a different kind of question, hypergeometrically: is a chemically expected class of compounds over-represented at the top of a ranking?

>>> class_enrichment(ranking_for_fruity, all_compounds, r"acetate|butanoate")
{'in_top': 5, 'class_size': 9, 'n_compounds': 61, 'n_drawn': 12,
 'p_value': 0.0004, 'matched': [...]}

At small sample sizes this is frequently the stronger evidence. Recovering the esters at the top of fruity is structure that noise does not produce, whereas a high R² on few products with several components very nearly is.

Note

Check where the ranking came from first. A one-component PLS has coefficient matrix outer(x_weights, y_loadings), so the absolute coefficients order identically for every attribute, and the attributes differ only in sign and magnitude. One-component solutions are common at small sample sizes, so treat this as the normal case rather than an edge case: an enrichment that looks attribute-specific may be one ranking reported many times.

API#

Kevin Dunn, 2010-2026. MIT License.

Permutation nulls that respond to signal, and a hypergeometric enrichment test.

At the sample sizes a sensory-to-chemistry study runs at, “is there anything here at all” is the hard question, and the obvious ways to answer it do not work. A high \(R^2\) on twelve products with three components is close to guaranteed. A count of how many variables exceed VIP 1 is not a test statistic at all: VIP is normalised so that \(\sum_j \text{VIP}_j^2 = K\) exactly, so that count describes the shape of the VIP distribution rather than the presence of a relationship, barely moves when the response is permuted, and yields a false-discovery rate near 100% on data that does contain signal.

What does work is asking a permutation what it can achieve:

check_predictive_signal() permutes the response across products, refits, and compares observed out-of-sample performance with what random reassignment reaches. Out-of-sample is what makes it bite; the same test on in-sample \(R^2\) would mostly measure model capacity. Ask this first: if it fails, nothing below can rescue the analysis.
count_discoveries_under_null() does the same for a whole selection procedure, counting discoveries under a permuted response, which audits the procedure as run rather than a model in isolation.
class_enrichment() asks whether a chemically expected class of compounds sits at the top of a ranking more often than chance allows. At small sample sizes this is frequently the stronger evidence: recovering the right class for an attribute is structure that noise does not produce. Its name is unchanged by the question-first naming used for the two functions above, because it already states its question, and because it tests a ranking you supply rather than a model this module fits: it sits alongside the pair, not on the same ladder.

References

Westerhuis, Hoefsloot, Smit and others, “Assessment of PLSDA cross validation”, Metabolomics, 4, 81-89, 2008, doi:10.1007/s11306-007-0099-6.

process_improve.multivariate._null.check_predictive_signal(x, y, fit_predict=None, n_components=2, n_perm=500, seed=0)[source]#

Test whether the predictors carry signal, against a shuffled response.

Ask this before asking which columns carry the signal. The question a permutation test can answer is “does this model predict held-out products better than it would if the responses were shuffled”, and unlike a VIP count (see vip()) the answer responds to whether signal is present.

Parameters:

x (pandas.DataFrame) – Predictor block, one row per product. Not permuted. Pass it unscaled when using the default fit_predict; see below.
y (pandas.DataFrame) – Response block, one row per product and one column per attribute.
fit_predict (callable, optional) –
fit_predict(x, y) returning out-of-sample predictions with one row per row of y and one column per attribute: a DataFrame, or an array of the same shape.

The default is leave-one-out cross-validated PLS with n_components, which is the right first thing to try and handles the two traps below for you. Supply your own to use a different model or a cheaper cross-validation scheme, in which case both traps become yours:
- The cross-validation scheme is yours to choose. Leave-one-out is right when every product is precious and needlessly expensive at hundreds of rows, where a permutation null multiplies the cost by n_perm. Pick deliberately.
- Re-derive the response constants inside each fold too. It is easy to nest the predictor preprocessing and forget the response. Centre and scale the response on the training rows, predict, then back-transform with those same training constants before the prediction is returned; otherwise the score is computed on a scale the fold never saw. Equivalently: do not scale x and y before calling, because the constants would then carry the held-out row.
n_components (int, default 2) – Latent components for the default fit_predict. Ignored when fit_predict is supplied. Capped at what the block supports.
n_perm (int, default 500) – Number of permutations. See the note on the p-value floor below.
seed (int, default 0) – Seed for the permutation order, so a reported p-value can be reproduced.

Returns:

One row per attribute (the columns of y), with columns:

attribute: The response name.
q2_observed: \(Q^2\) from fit_predict(x, y).
q2_null_mean, q2_null_p95: The mean and 95th percentile of \(Q^2\) under permutation. The 95th percentile is the more useful of the two: it is roughly the performance a shuffled response reaches one time in twenty.
p_value: The fraction of permutations reaching at least q2_observed, computed as (1 + count) / (n_perm + 1).
n_permutations: Permutations that produced a usable \(Q^2\).

Return type:

pandas.DataFrame

Raises:

ValueError – If the blocks disagree on rows, n_perm is below 1, or fit_predict returns the wrong shape.

Notes

This is expensive, and irreducibly so. Every permutation refits the model once per cross-validation fold, so the default leave-one-out scheme costs n_perm * n_products fits: on twenty products with the default 500 permutations that is ten thousand fits, a few minutes. The cost is the method, not the implementation - an out-of-sample null has to refit to be out-of-sample. Develop with n_perm around 50, then raise it for the number you intend to report, keeping the p-value floor below in mind.

The p-value uses the (1 + count) / (n_perm + 1) form, which counts the observed statistic as one of its own null draws so the result can never be exactly zero. That also puts a floor on it: the smallest attainable p-value is 1 / (n_perm + 1). The default 500 permutations cannot report anything below 0.002, so choose n_perm with the threshold you intend to apply in mind, and remember that a multiplicity correction over many attributes needs a floor well below the corrected threshold.

See also

count_discoveries_under_null: Audits a whole selection procedure rather than one model.

Examples

The default cross-validates PLS for you, on unscaled blocks:

>>> check_predictive_signal(chem, sensory_means, n_perm=999)

Supply your own when the model or the fold scheme has to change:

>>> def fit_predict(x, y):
...     return loo_predictions(x, y, n_components=2)  # your own CV loop
>>> check_predictive_signal(chem, sensory_means, fit_predict, n_perm=999)

process_improve.multivariate._null.class_enrichment(ranked, all_names, pattern, top_n=12)[source]#

Test whether a named class of compounds is over-represented at the top of a ranking.

At small sample sizes this is frequently stronger evidence than \(R^2\) or \(Q^2\). Recovering the chemically expected class for an attribute (the esters at the top of “fruity”, the pyrazines at the top of “roasted”) is structure that noise does not produce, whereas a high \(R^2\) on few products with several components very nearly is.

The test is hypergeometric: given n_compounds compounds of which class_size belong to the class, how surprising is it to draw in_top of them in the top top_n?

Note

Check where the ranking came from before reading much into a per-attribute result. A one-component PLS has a coefficient matrix outer(x_weights, y_loadings), so the absolute coefficients order identically for every attribute and the attributes differ only in sign and magnitude. One-component solutions are common at small sample sizes, so this is the normal case rather than an edge case: an enrichment that looks attribute-specific may be one ranking reported many times.

Parameters:

ranked (sequence of str) – Compound names, most important first. May be shorter than all_names.
all_names (sequence of str) – Every compound that could have been ranked: the population the draw comes from. Must not contain duplicates.
pattern (str) – Regular expression defining the class, matched against each name with re.search() (so a plain substring works).
top_n (int, default 12) – How much of the ranking counts as “the top”. Truncated to the length of ranked when it is shorter.

Returns:

With keys in_top (class members among the top n_drawn), class_size (class members in the population), n_compounds (population size), n_drawn (the effective number drawn), p_value (the hypergeometric upper-tail probability of seeing at least in_top), and matched (the class members found, in ranked order).

Return type:

dict

Raises:

ValueError – If all_names has duplicates, ranked contains a name not in all_names, top_n is below 1, or pattern is not a valid regular expression.

Examples

>>> class_enrichment(ranking_for_fruity, all_compounds, r"acetate|butanoate")
{'in_top': 5, 'class_size': 9, 'n_compounds': 61, 'n_drawn': 12, ...}

process_improve.multivariate._null.count_discoveries_under_null(select, x, y, n_perm=500, seed=0)[source]#

Count how many of these findings shuffling alone would have produced.

check_predictive_signal() tests one model; this tests the whole procedure that produced it. Pass a callable that runs filtering, transformation, scaling and selection end to end, and it is re-run in full on each permuted response. The average number of discoveries under permutation, divided by the number actually made, is an empirical false-discovery rate for the procedure as run, which is the only version of it worth quoting: a procedure whose selection step is honest but whose filtering step peeked at the response has an FDR that no formula recovers.

Parameters:

select (callable) – select(x, y) returning the names it selected: any iterable of hashable labels. It should be deterministic given its inputs; if it is not, the null counts pick up its own randomness and the reported FDR means nothing. A SpecificationWarning is raised when a repeat call on the observed data disagrees with the first.
x (pandas.DataFrame) – Predictor block, one row per product. Not permuted.
y (pandas.DataFrame) – Response block, one row per product. Whole rows are permuted, so each product keeps its full response vector.
n_perm (int, default 500) – Number of permutations.
seed (int, default 0) – Seed for the permutation order.

Returns:

With keys:

observed: Number of names selected on the real response.
null_mean, null_p95: Mean and 95th percentile of the selection count under permutation.
empirical_fdr: Deprecated since version 1.77.0: The old name for the ratio below, kept for one release and clipped to [0, 1] as it always was. Prefer null_to_observed_ratio, which is unclipped. Will be removed in 2.0.0.
null_to_observed_ratio: null_mean / observed; NaN when nothing was selected. Read it as “of the names this procedure returned, about this fraction is what shuffling alone would have produced”. It is not clipped: a value above 1 means shuffling found more than the real response did, which is the strongest evidence this procedure has nothing, and clipping it away would hide exactly the case worth seeing.
null_counts: The per-permutation counts, as an ndarray, for plotting the null.
selected: The names selected on the real response, in the order the callable returned them.

Return type:

dict

Raises:

ValueError – If the blocks disagree on rows or n_perm is below 1.