Interaction terms (provisional)#

Warning

The three functions in process_improve.interactions are unvalidated on real data. They exist to test conditions that a small observational study rarely meets, so this code may go unreached for a long time; its coverage is unit tests only, and nothing here has been run against a real product-by-compound block with a real sensory response. Treat the API as subject to change, and read any result from it as a hypothesis rather than a finding.

An interaction between two chemical predictors is a real thing to look for: a compound whose perceptual effect depends on the level of another is exactly what a linear additive model misses. Whether the data can support the claim is a separate matter, and usually the binding one.

Does the pair support a term at all?#

An interaction term says the effect of A depends on the level of B. Estimating that needs products where A is high and B is low, and products where the reverse holds, as well as the two agreeing corners. With one corner empty the term is fitted from three points of support and will report whatever noise lives there.

>>> from process_improve.interactions import pair_coverage
>>> covered, detail = pair_coverage(x["linalool"].to_numpy(), x["geraniol"].to_numpy())
>>> covered, detail["low_high"], detail["high_low"], detail["correlation"]
(False, 1, 0, 0.91)

Two variables that co-vary occupy only the agreeing corners and fail this check. That is the correct answer, not a defect to work around. Their interaction is not identifiable from these observations, and no amount of regularisation makes it so. The returned correlation is nearly always the explanation.

Each variable is split at its own median, so the marginal split is balanced by construction and only the joint distribution can fail.

Building the terms#

The order is not negotiable: transform, centre and scale, multiply, then re-centre and re-scale the products.

from process_improve.interactions import interaction_terms

terms, constants = interaction_terms(x_scaled, [("linalool", "geraniol")])

The second pass is not tidiness. The product of two standardised columns is not itself centred: for approximately bivariate normal parents its mean is the parents’ correlation \(r\) and its variance is \(1 + r^2\). Skipping the re-centring therefore leaks correlation into the intercept, and skipping the re-scaling hands the model a column with more variance than a genuine predictor would have, inflating exactly the pairs whose interactions deserve the least trust.

constants records center, divisor and parent_correlation per term, so a held-out block can be built with training constants: multiply the same parents, subtract center, divide by divisor. Parents that do not look centred and unit-variance draw a SpecificationWarning, since both the reasoning above and the reported correlation depend on them being so.

How stable is a selection?#

A selection made once on all the data is a selection made once. stability_selection() repeats it on complementary half-samples and reports how often each name comes back:

>>> stability_selection(select, x_scaled, sensory_means, n_iter=50).head(3)
    name  n_selected  n_subsamples  selection_frequency
0  ethyl_hexanoate         98           100                 0.98
1        linalool          61           100                 0.61
2         nonanal          14           100                 0.14

The halves are complementary: each split is used in both directions, so every product appears in exactly half the subsamples and the two runs of a split share no rows.

References#

Meinshausen and Buhlmann, “Stability selection”, Journal of the Royal Statistical Society: Series B, 72(4), 417-473, 2010, doi:10.1111/j.1467-9868.2010.00740.x.

Shah and Samworth, “Variable selection with error control: another look at stability selection”, Journal of the Royal Statistical Society: Series B, 75(1), 55-80, 2013, doi:10.1111/j.1467-9868.2011.01034.x.

API#

  1. Kevin Dunn, 2010-2026. MIT License.

Two-factor interaction terms, and whether the data can support them.

Warning

Provisional. These three functions are unvalidated on real data. They exist to test conditions that a small observational study rarely meets: a pair of predictors whose observations populate all four corners of their plane, and a selection procedure stable enough for its choices to mean something. Because those conditions are rarely met, this code may go unreached for a long time after shipping, and its coverage is unit tests only. Nothing here has been run against a real product-by-compound block with a real sensory response. Treat the API as subject to change, and read any result from it as a hypothesis rather than a finding.

An interaction between two chemical predictors is a real thing to look for: a compound whose perceptual effect depends on the level of another is exactly what a linear additive model misses. But an interaction term is only estimable when the data actually separates the four combinations of high and low, and in an observational set of thirty products it usually does not. pair_coverage() says so before a model is fitted rather than after.

The ordering inside interaction_terms() is not negotiable: transform, centre and scale, multiply, then re-centre and re-scale the products. The product of two standardised variables is not itself centred. For approximately bivariate normal columns its mean is the parents’ correlation \(r\) and its variance is \(1 + r^2\), so skipping the second centring leaks correlation into the intercept and systematically inflates the columns belonging to correlated pairs, exactly the pairs whose interactions are least trustworthy.

References

Meinshausen and Buhlmann, “Stability selection”, Journal of the Royal Statistical Society: Series B, 72(4), 417-473, 2010, doi:10.1111/j.1467-9868.2010.00740.x.

Shah and Samworth, “Variable selection with error control: another look at stability selection”, Journal of the Royal Statistical Society: Series B, 75(1), 55-80, 2013, doi:10.1111/j.1467-9868.2011.01034.x.

process_improve.interactions.interaction_terms(x_log, pairs)[source]#

Build centred, scaled product terms from centred, scaled predictors.

The order matters and is not negotiable: the parents must already be transformed, centred and scaled; this function multiplies them and then centres and scales the products again. A product of two standardised columns has mean \(r\) and variance \(1 + r^2\) (for approximately bivariate normal parents), so a term skipping the second pass carries its parents’ correlation into the intercept and arrives at the model with more variance than a genuine predictor would, inflating exactly the pairs whose interactions deserve the least trust.

Check pair_coverage() first. This function will happily build a term for a pair that supports nothing.

Parameters:

  • x_log (pandas.DataFrame) – Predictors, already transformed, centred and scaled (see process_improve.chemistry). A SpecificationWarning is raised when they do not look centred and unit-variance, because the reasoning above, and the parent_correlation column below, assume it.

  • pairs (sequence of (str, str)) – The pairs to build. Both names must be columns of x_log. A pair naming the same column twice gives a quadratic term, whose mean before re-centring is 1 rather than \(r\); that is allowed, but it is not an interaction.

Returns:

(terms, constants) – terms has one column per pair, named "a_x_b", on the rows of x_log, centred and unit-variance scaled. constants has one row per term with term, left, right, center, divisor and parent_correlation. divisor is divided by, not multiplied by, and the pair lets a held-out block be built with training constants: multiply the same parents, then subtract center and divide by divisor.

Return type:

tuple of pandas.DataFrame

Raises:

ValueError – If pairs is empty, names a column that is not in x_log, repeats a pair, or would produce a term name that collides with an existing column.

Examples

>>> terms, constants = interaction_terms(x_scaled, [("linalool", "geraniol")])
>>> terms.mean().abs().max() < 1e-12
True
process_improve.interactions.pair_coverage(x_a, x_b, min_per_corner=4)[source]#

Ask whether the observations populate all four corners of the (A, B) plane.

An interaction term claims that the effect of A depends on the level of B. Estimating that claim needs products where A is high and B is low, and products where the reverse holds, as well as the two agreeing corners. With one corner empty the term is fitted from three points of support and will report whatever the noise there suggests.

Two variables that co-vary populate only the agreeing corners and will fail this check. That is the correct answer, not a defect to work around: their interaction is not identifiable from these observations, and no amount of regularisation makes it so.

Parameters:

  • x_a (numpy.ndarray) – The two predictors, one value per product, the same length. Rows missing either value are dropped.

  • x_b (numpy.ndarray) – The two predictors, one value per product, the same length. Rows missing either value are dropped.

  • min_per_corner (int, default 4) – How many observations a corner needs before it counts as populated.

Returns:

(covered, detail) – covered is True when every corner holds at least min_per_corner observations. detail carries low_low, low_high, high_low, high_high (the corner counts, with the first word describing A), n (rows used), min_per_corner, threshold_a and threshold_b (the median splits), and correlation (Pearson’s r between the two, which is usually the explanation when the answer is False).

Return type:

tuple of (bool, dict)

Raises:

ValueError – If the two inputs have different lengths, or min_per_corner is below 1.

Notes

Each variable is split at its own median, so the marginal split is balanced by construction and only the joint distribution can fail. A value exactly at the median counts as low.

Examples

>>> covered, detail = pair_coverage(x["linalool"].to_numpy(), x["geraniol"].to_numpy())
>>> covered, detail["correlation"]
(False, 0.91)
process_improve.interactions.stability_selection(select, x, y, n_iter=100, seed=0)[source]#

Report how often each predictor is selected across complementary half-samples.

A selection made once on all the data is a selection made once. Repeating it on random halves, and reporting how often each name comes back, separates the choices the data supports from the ones that depended on which products happened to be in the set. The halves are complementary: each split is used in both directions, so every product appears in exactly half of the subsamples and the two runs of a split share no rows.

Parameters:

  • select (callable) – select(x, y) returning the names it selected, as any iterable of labels. Called 2 * n_iter times, so keep it cheap.

  • x (pandas.DataFrame) – Predictor block, one row per product. Its columns are the universe: a name the callable returns that is not a column here is an error.

  • y (pandas.DataFrame) – Response block, one row per product. Sub-sampled with x, row for row.

  • n_iter (int, default 100) – Number of complementary splits, so 2 * n_iter calls to select.

  • seed (int, default 0) – Seed for the splits, so a reported frequency can be reproduced.

Returns:

One row per column of x, sorted by frequency then name, with name, n_selected (subsamples that chose it), n_subsamples (always 2 * n_iter) and selection_frequency (the ratio).

Return type:

pandas.DataFrame

Raises:

ValueError – If the blocks disagree on rows, n_iter is below 1, there are fewer than four products to split, or select returns a name that is not a column of x.

Examples

>>> frequencies = stability_selection(select, x_scaled, sensory_means, n_iter=50)
>>> frequencies.query("selection_frequency > 0.6")["name"].tolist()