Chemistry: preparing a product-by-compound block#

A sensory-to-chemistry study has two tables with one row per product: a block of sensory attribute means, and a block of concentrations or integrated peak areas. Relating the second to the first with PLS is the easy part. Getting the second one ready is where the answers are quietly won or lost, and process_improve.chemistry is the four decisions that make up that work.

The order is fixed#

trim  ->  transform  ->  centre  ->  scale

Trim first, because a compound detected in two of forty products has no concentration worth transforming. Transform before centring, because the range ratio that decides between a log and a linear scale is a property of the raw values. Centre before scaling, because a scaling constant estimated around the wrong centre is the wrong constant.

from process_improve.chemistry import (
    apply_transform,
    center_and_scale,
    classify_zero_states,
    normalisation_check,
    trim_by_prevalence,
)

states = classify_zero_states(chem, lod={"linalool": 0.02})
totals, outside = normalisation_check(chem)

kept, dropped, presence = trim_by_prevalence(chem, min_nonzero=3)
transformed, applied = apply_transform(kept, lod={"linalool": 0.02})
scaled, constants = center_and_scale(transformed, presence[kept.columns])

scaled is now ready for PLS with scale=False, since it is already centred and scaled.

A zero is not self-describing#

A zero in a concentration table is one of two completely different things:

  • a rounded zero, also called left-censored: the compound is present, at a concentration the instrument could not resolve;

  • an essential zero: the compound is genuinely absent.

They want opposite handling. A rounded zero should be substituted with something below the detection limit and then modelled as a small quantity; an essential zero is a categorical fact and modelling it as a small quantity invents chemistry that is not there. The distinction cannot be recovered from an exported table, so classify_zero_states() does not guess:

>>> classify_zero_states(chem).set_index("compound")["zero_state"].unique()
array(['unknown'], dtype=object)

unknown is the default, and it is not a failure state: it records that nobody has yet said which of the other two applies. Never default to censored. Declaring a detection limit is the claim that non-detects lie below it, so passing one in lod classifies the compound as rounded; nothing else does that on your behalf.

A trimmed compound is not a discarded compound#

trim_by_prevalence() returns three frames, and the third is the interesting one:

kept, dropped, presence = trim_by_prevalence(chem, min_nonzero=3)

presence covers every compound, kept and dropped alike. For a rare compound the binary fingerprint often carries more than the concentration ever could: “this compound appears in exactly the three products the panel called green” is a finding, and a column of thirty-seven zeros and three numbers states it badly. A missing measurement stays NaN in that layer rather than becoming a zero, because never measured and measured-as-absent are different.

Choosing a transform#

choose_transform() reads the range ratio of the detected values: largest over smallest, ignoring zeros and missing cells. A compound spanning orders of magnitude is multiplicative and belongs on a log scale; one varying by a factor of two or three is additive and does not. In between, the honest answer is "ambiguous", and apply_transform() resolves it with a caller-chosen default rather than a coin toss:

>>> choose_transform(chem["limonene"])
'log'
>>> choose_transform(chem["ethanol"])
'linear'

A log-transformed compound has its non-detects substituted, by half the declared detection limit or half the smallest value seen, before the log is taken. Detected values survive exactly. A linear compound passes through untouched: a zero on a linear scale is a usable number and needs no substitution.

detected_only is off for a reason#

The rule that centring and scaling constants must not be estimated from imputed values is sound, and easy to over-apply. Where no imputation has happened, the zeros are real observations of “not detected”, and excluding them puts every one of those zeros many standard deviations below a centre computed from a handful of detected values. The column becomes effectively binary with a very large magnitude, and since PLS follows variance, the components then track how sparse a variable is rather than how it relates to the response. Every attribute comes back with the same handful of rare compounds at the top of its list, which reads as a finding and is an artefact.

So center_and_scale() defaults detected_only=False. Switch it on after apply_transform() has substituted non-detects for a log-scaled compound, and not before. The detected layer is a required argument either way, so the flag cannot be reached without the mask in hand.

Preprocessing a held-out fold honestly#

Both fitting steps have a replay partner:

train, test = chem.iloc[train_rows], chem.iloc[test_rows]

kept, _dropped, presence = trim_by_prevalence(train, min_nonzero=3)
train_t, applied = apply_transform(kept)
train_s, constants = center_and_scale(train_t, presence[kept.columns])

test_t = apply_fitted_transform(test[kept.columns], applied)
test_s = apply_fitted_center_scale(test_t, constants)

Without that pair, honest nested cross-validation is not possible: re-deriving the transform offsets and the scaling constants from the test rows would let those rows influence their own preprocessing, and the cross-validated score would then be measuring something other than out-of-sample performance.

The constants table names its column divisor and is divided by, not multiplied by. That is deliberate: the standalone scale() returns a multiplier while center() returns a subtrahend, and the two have been confused often enough that the new table leaves no room for it.

See also

Panel diagnostics: can this attribute be modelled? for the matching question on the sensory side: whether an attribute can be modelled as an intensity at all.

API#

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

Preprocessing a product-by-compound block of concentrations or peak areas.

One row per product, one column per compound, values being concentrations or integrated peak areas. Getting such a block ready for a PLS against a block of sensory attributes is four decisions, and the order they are taken in is fixed:

trim  ->  transform  ->  centre  ->  scale

Trimming first, because a compound seen in two of forty products has no concentration worth transforming. Transforming before centring, because the range ratio that decides between a log and a linear scale is a property of the raw values. Centring before scaling, because a scaling constant estimated around the wrong centre is the wrong constant.

Three things in here are easy to get wrong, and each is wrong quietly:

  • A zero is not self-describing. It is either a concentration below a detection limit (a rounded zero, the compound is there) or a compound that is genuinely absent (an essential zero). Those need opposite handling, and the distinction cannot be recovered from an exported table. classify_zero_states() therefore defaults to "unknown" and makes the caller declare. It never defaults to censored.

  • A trimmed compound is not a discarded compound. For a rare compound the binary fingerprint of where it appears at all often carries more than its concentration does, so trim_by_prevalence() returns a presence layer covering every compound, not only the kept ones.

  • Preprocessing constants must not see the test rows. Every fitting function here has an apply_fitted_* partner that replays the constants computed elsewhere, so held-out rows can be preprocessed with training constants alone. Without that pair, honest nested cross-validation is not possible: the transform offsets and the scaling constants would both have seen the row being predicted.

References

Martin-Fernandez, Barcelo-Vidal and Pawlowsky-Glahn, “Dealing with zeros and missing values in compositional data sets using nonparametric imputation”, Mathematical Geology, 35(3), 253-278, 2003.

van den Berg, Hoefsloot, Westerhuis and others, “Centering, scaling, and transformations: improving the biological information content of metabolomics data”, BMC Genomics, 7, 142, 2006, doi:10.1186/1471-2164-7-142.

process_improve.chemistry.preprocessing.ZERO_STATES: tuple[str, ...] = ('rounded', 'essential', 'unknown')#

it records that nobody has yet said which of the other two applies.

Type:

The three zero states. unknown is the default and is not a failure

process_improve.chemistry.preprocessing.TRANSFORM_RULES: tuple[str, ...] = ('log', 'linear', 'ambiguous')#

The transform rules choose_transform() can return.

process_improve.chemistry.preprocessing.SCALING_METHODS: tuple[str, ...] = ('autoscale', 'pareto')#

Scaling methods understood by center_and_scale().

process_improve.chemistry.preprocessing.classify_zero_states(chem, declared=None, lod=None)[source]#

Record, per compound, what a zero in that column is taken to mean.

A zero is either rounded (left-censored: the compound is present at a concentration below the detection limit) or essential (structurally absent: the compound is not there). The two want opposite handling, and no amount of looking at an exported table recovers which one applies, so this function does not guess. A compound nobody has spoken for is "unknown".

Never default to censored. Classifying a zero as left-censored asserts a latent value below a detection limit, which is a claim about the chemistry, not a convenience. Declaring a detection limit for a compound is that claim, so passing one in lod classifies the compound as "rounded".

Parameters:

  • chem (pandas.DataFrame) – Products (rows) by compounds (columns) of concentrations or peak areas.

  • declared (dict or None) – Explicit per-compound states, {compound: "rounded" | "essential" | "unknown"}. Takes precedence over lod. Compound names not in the block are an error, since a typo would otherwise silently leave the compound unknown.

  • lod (dict or None) – Per-compound limits of detection, {compound: float}. A compound with a finite, positive limit and no entry in declared is classified "rounded".

Returns:

One row per compound, in the block’s column order, with columns:

compound

The compound name.

zero_state

"rounded", "essential" or "unknown".

source

How the state was arrived at: "declared", "lod", or "default".

lod

The declared limit of detection, or NaN.

n_zero, n_nonzero, n_missing

Cell counts, so a compound with no zeros at all (whose state is therefore moot) is easy to spot.

Return type:

pandas.DataFrame

Raises:

ValueError – If the block is empty, a declared state is not one of ZERO_STATES, or declared / lod names a compound the block does not have.

Examples

>>> states = classify_zero_states(chem, lod={"linalool": 0.02})
>>> states.query("zero_state == 'unknown' and n_zero > 0")["compound"].tolist()
process_improve.chemistry.preprocessing.trim_by_prevalence(chem, min_nonzero=3)[source]#

Split the block by how often each compound is seen at all.

A compound detected in one or two products has no concentration worth modelling: any apparent relationship is a line through two points. It is not thereby uninformative, though, so nothing is thrown away. The third return value is a presence layer covering every compound, kept and dropped alike, which for a rare compound is often the more useful representation: “this compound appears in exactly the three products the panel called ‘green’” is a finding, and it is one a concentration column full of zeros states badly.

Parameters:

  • chem (pandas.DataFrame) – Products (rows) by compounds (columns).

  • min_nonzero (int, default 3) – Keep a compound when it is non-zero in at least this many products.

Returns:

(kept, dropped, presence) – kept and dropped partition the columns of chem, both keeping every row. presence has the same shape as chem and holds 1.0 where a compound was detected, 0.0 where it was not, and NaN where the measurement is missing: a missing measurement is not an absence.

Return type:

tuple of pandas.DataFrame

Raises:

ValueError – If the block is empty or min_nonzero is negative.

Examples

>>> kept, dropped, presence = trim_by_prevalence(chem, min_nonzero=3)
>>> presence[dropped.columns].sum().sort_values(ascending=False).head()
process_improve.chemistry.preprocessing.normalisation_check(chem, factor=1.8)[source]#

Report the row totals, and which of them sit outside a fold-band around the median.

A product-by-compound block that has been normalised to a constant sum has row totals that are all equal; one that has not may still be fine, but a row whose total is several fold away from its neighbours usually means something mechanical rather than chemical (a different injection volume, a dilution that was not recorded, an integration that dropped a peak). Either way the answer changes what the rest of the pipeline should do, so it is worth looking before transforming.

Parameters:

  • chem (pandas.DataFrame) – Products (rows) by compounds (columns).

  • factor (float, default 1.8) – Half-width of the accepted band, as a fold change: a row is reported when its total is above median * factor or below median / factor. Must be greater than 1.

Returns:

(totals, outside) – totals is the row sum for every product, NaN for a row with no measurements at all. outside is the subset of totals beyond the band, in the block’s row order; it is empty when every row is inside.

Return type:

tuple of pandas.Series

Raises:

ValueError – If the block is empty, factor is not greater than 1, or the median row total is not positive (which leaves no band to compare against).

Examples

>>> totals, outside = normalisation_check(chem)
>>> outside / totals.median()  # how far out, as a fold change
process_improve.chemistry.preprocessing.choose_transform(col, ratio_log=10.0, ratio_linear=3.0)[source]#

Decide whether a compound should be modelled on a log or a linear scale.

The decision is made on the range ratio of the detected values: the largest divided by the smallest, ignoring zeros and missing cells. A compound spanning orders of magnitude is multiplicative and belongs on a log scale; one varying by a factor of two or three is additive and does not. In between there is no evidence either way and the honest answer is "ambiguous", which apply_transform() resolves with a caller-chosen default rather than a coin toss.

Parameters:

  • col (pandas.Series) – One compound’s values across products.

  • ratio_log (float, default 10.0) – A range ratio at or above this chooses "log".

  • ratio_linear (float, default 3.0) – A range ratio at or below this chooses "linear".

Returns:

One of TRANSFORM_RULES. "linear" is returned when a log is not applicable at all (a negative value present, or fewer than two detected values to form a ratio from).

Return type:

str

Raises:

ValueError – If ratio_log is not strictly greater than ratio_linear, which would make the two rules overlap.

Examples

>>> choose_transform(chem["limonene"])
'log'
process_improve.chemistry.preprocessing.apply_transform(chem, lod=None, ambiguous='linear', ratio_log=10.0, ratio_linear=3.0)[source]#

Transform every compound by the rule choose_transform() picks for it.

A "log" compound has its non-detects replaced (by half the declared detection limit, or half the smallest value seen) and is then taken to base 10. A "linear" compound passes through untouched: a zero on a linear scale is a usable number and needs no substitution.

The substitution is an imputation, and it is the only imputation this module performs. It is what makes detected_only=True in center_and_scale() meaningful, and what makes it meaningless before this step has run.

Parameters:

  • chem (pandas.DataFrame) – Products (rows) by compounds (columns). Trim first: see trim_by_prevalence().

  • lod (dict or None) – Per-compound limits of detection, {compound: float}. Used to set the substitution value for a log-transformed compound.

  • ambiguous ({"linear", "log"}, default "linear") – What to do with a compound whose range ratio decides nothing. The default keeps the values as they are, which is the smaller claim.

  • ratio_log (float) – Passed to choose_transform().

  • ratio_linear (float) – Passed to choose_transform().

Returns:

(transformed, applied) – transformed has the shape and labels of chem. applied has one row per compound with columns compound, rule ("log" or "linear"), offset (the substitution value; 0.0 for a linear compound), range_ratio, and chosen_by ("range_ratio" or "ambiguous_default"). Feed it to apply_fitted_transform() to replay these decisions on held-out rows.

Return type:

tuple of pandas.DataFrame

Raises:

ValueError – If the block is empty, ambiguous is not "linear" or "log", lod names a compound the block does not have, or the two ratio thresholds overlap.

Examples

>>> transformed, applied = apply_transform(kept, lod={"linalool": 0.02})
>>> applied.query("rule == 'log'")["compound"].tolist()
process_improve.chemistry.preprocessing.apply_fitted_transform(chem, applied)[source]#

Replay a transform table computed elsewhere, without re-deriving it.

This is the half of apply_transform() that must be used on held-out rows. Re-deriving the rule and the offset from the test rows would let those rows influence their own preprocessing, and the cross-validated score would then be measuring something other than out-of-sample performance.

Parameters:
Returns:

The transformed block, with the labels of chem.

Return type:

pandas.DataFrame

Raises:

ValueError – If applied is missing a required column, names a rule other than "log" or "linear", or has no entry for some column of chem.

Examples

>>> train_t, applied = apply_transform(chem.iloc[train_rows])
>>> test_t = apply_fitted_transform(chem.iloc[test_rows], applied)
process_improve.chemistry.preprocessing.center_and_scale(transformed, detected, method='autoscale', *, detected_only=False)[source]#

Centre and scale a transformed block, returning the constants for replay.

"autoscale" divides by the standard deviation (ddof=1, matching MCUVScaler), giving every compound equal say. "pareto" divides by its square root, which keeps some of the original variance structure and is the usual choice when the large peaks are meant to stay large.

Parameters:

  • transformed (pandas.DataFrame) – Output of apply_transform().

  • detected (pandas.DataFrame) – A 1 / 0 / NaN layer of the same shape, marking which cells were actually observed; the presence return of trim_by_prevalence(), restricted to the same columns. Only read when detected_only=True, but required either way so that a caller cannot switch the flag on without having the mask to hand.

  • method ({"autoscale", "pareto"}, default "autoscale") – The divisor: the standard deviation, or its square root.

  • detected_only (bool, keyword-only, default False) –

    Compute the centring and scaling constants from the detected cells alone.

    Leave this off unless imputation has actually run. The rule it implements, that constants must not be estimated from imputed values, is sound and easy to over-apply. On a column whose zeros are real observations of “not detected”, excluding them puts every one of those zeros many standard deviations below a centre estimated from a handful of detected values. The column becomes effectively binary with a very large magnitude, and since PLS follows variance, the components then track how sparse a variable is rather than how it relates to the response: every attribute comes back with the same handful of rare compounds at the top of its list. Switch it on after apply_transform() has substituted non-detects for a log-scaled compound, and not before.

Returns:

(scaled, constants) – scaled has the labels of transformed. constants has one row per compound with columns compound, center, divisor, method and n_used. divisor is what the centred values were divided by, not a multiplier; feed the table to apply_fitted_center_scale() rather than reapplying it by hand.

Return type:

tuple of pandas.DataFrame

Raises:

ValueError – If the block is empty, method is not one of SCALING_METHODS, or detected does not line up with transformed.

Examples

>>> scaled, constants = center_and_scale(transformed, presence[transformed.columns])
>>> constants.sort_values("divisor").head()
process_improve.chemistry.preprocessing.apply_fitted_center_scale(transformed, constants)[source]#

Replay centring and scaling constants computed elsewhere.

The partner of center_and_scale(), for held-out rows. divisor is divided by, not multiplied by: it is the standard deviation (or its square root), the same quantity center_and_scale() divided the training rows by.

Parameters:
Returns:

The scaled block, with the labels of transformed.

Return type:

pandas.DataFrame

Raises:

ValueError – If constants is missing a required column, has no entry for some column of transformed, or carries a non-positive divisor.

Examples

>>> train_s, constants = center_and_scale(train_t, train_presence)
>>> test_s = apply_fitted_center_scale(test_t, constants)