clophfit.fitting.residual_tests#

Formal residual diagnostics: the four assumptions a least-squares fit makes.

Each assumption has a picture that shows it failing and a test that measures it:

Assumption	What failure looks like	Test
Homoscedasticity	funnel / fan against prediction	Breusch-Pagan (Koenker), White
Independence	waves along the titration	Durbin-Watson, lag-1 r, runs
Normality	curved QQ, skewed histogram	Shapiro-Wilk, KS, Anderson-Darling
No outliers	isolated extreme points	externally studentized t, PRESS

Two things make these tests easy to misread on titration plates, and every function here reports accordingly:

They are pooled over hundreds of points, so any real effect, however small, reaches p ~ 0. Each test therefore also returns an effect size - R^2, the DW statistic, W, D, a rate or a ratio - and those, not the p-values, are what compare models.
The residuals belong to a fit. A plate spends several parameters per well, so residuals shrink and pick up induced correlation even when the noise is independent and calibrated. Durbin-Watson is not expected at exactly 2, and the Kolmogorov-Smirnov test against N(0, 1) mixes shape with that shrinkage (see residuals.qq_fit); the Lilliefors form, with location and scale estimated, tests shape alone.

The outlier and PRESS measures need each point’s leverage, computed per well from the weighted Jacobian of binding_1site() at the supplied parameters. For a Bayesian multi-well fit those are posterior means and the latent per-well pH offset is not counted, so leverage is an approximation there; for a classical per-well fit it is exact.

Functions#

`breusch_pagan`(z, yhat, x)	Test whether the squared residual follows the prediction or x (Breusch-Pagan).
`white_test`(z, yhat, x)	White's test: squared residual on ŷ, x, their squares and cross product.
`durbin_watson`(table, *[, z_col])	Durbin-Watson pooled within series: each well's label ordered along the titration.
`normality_tests`(z, *[, n_mc, seed])	Shapiro-Wilk, Kolmogorov-Smirnov against N(0, 1), Lilliefors and Anderson-Darling.
`durbin_watson_null`(table, params, *[, is_ph])	Durbin-Watson and lag-1 correlation expected of fitted residuals, per label.
`runs_test`(table, *[, z_col])	Wald-Wolfowitz runs test on residual signs, pooled within series.
`runs_test_null`(table, params, *[, is_ph, n_mc, seed])	Distribution of the pooled runs count for fitted residuals, per label.
`well_leverage`(table, params, *[, is_ph])	Leverage `h_ii` of every residual, from its well's weighted binding fit.
`studentized_outliers`(table, params, *[, is_ph, alpha])	Externally studentized residual of every point, and whether it is an outlier.
`press_statistics`(table, params, *[, is_ph])	PRESS: the prediction error of every point when it is left out of its well's fit.
`residual_tests`(table[, params, is_ph, n_mc, seed])	All four assumption checks, one row per label.

Module Contents#

clophfit.fitting.residual_tests.breusch_pagan(z, yhat, x)#

Test whether the squared residual follows the prediction or x (Breusch-Pagan).

The studentized form, n R^2 of z^2 regressed on the regressors, is used rather than the original, which assumes normal errors and rejects on kurtosis alone.

Parameters:

z (np.ndarray) – Standardised residuals.
yhat (np.ndarray) – Predictions.
x (np.ndarray) – Titration coordinate (pH or concentration).

Returns:

lm, df (2), p and r2 - the fraction of the squared residuals’ variance the fan explains, the effect size to compare.

Return type:

dict[str, float]

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> y = rng.uniform(100, 1000, 800)
>>> fan = breusch_pagan(rng.normal(0, y / 500), y, rng.uniform(5, 9, 800))
>>> fan["p"] < 1e-6
True

clophfit.fitting.residual_tests.white_test(z, yhat, x)#

White’s test: squared residual on ŷ, x, their squares and cross product.

Catches a variance that bends with the signal - high at both ends, say - which the linear Breusch-Pagan design cannot see.

Parameters:

z (np.ndarray) – Standardised residuals.
yhat (np.ndarray) – Predictions.
x (np.ndarray) – Titration coordinate.

Returns:

lm, df (5), p and r2.

Return type:

dict[str, float]

clophfit.fitting.residual_tests.durbin_watson(table, *, z_col='std_res')#

Durbin-Watson pooled within series: each well’s label ordered along the titration.

DW = sum((e_t - e_{t-1})^2) / sum(e_t^2) with the differences taken only inside a (well, label) series, so a jump between wells is never counted. Well below its null, neighbouring steps share an error (a wave, or a misplaced x); above it, they alternate. The null is not 2 for series this short: independent errors give 2 (n - 1) / n per series - 1.71 for a seven-step titration - which dw_expected reports. Fitted residuals are not independent even then: the fit makes neighbours alternate, and for a binding curve that lifts the null well above 2. Use durbin_watson_null() for the null of the fitted residuals.

Parameters:

table (pd.DataFrame) – Residual table with well, label, step (or x) and z_col.
z_col (str) – The residual column.

Returns:

dw, dw_expected (its value for independent errors of these series lengths), the pooled lag-1 correlation lag1_r, n_pairs and p, a two-sided normal-approximation p-value for lag1_r against 0.

Return type:

dict[str, float]

clophfit.fitting.residual_tests.normality_tests(z, *, n_mc=999, seed=0)#

Shapiro-Wilk, Kolmogorov-Smirnov against N(0, 1), Lilliefors and Anderson-Darling.

Shapiro-Wilk W measures shape only; it is the most powerful of the three, so on hundreds of points it rejects small departures - read W.
KS against N(0, 1) tests shape and scale: a calibrated fit still fails it when its residuals are shrunk by the parameters it spent.
Lilliefors is KS with location and scale estimated, p by Monte Carlo (scipy.stats.goodness_of_fit): shape only, like Shapiro-Wilk, but weighted towards the centre rather than the tails.
Anderson-Darling, also with location and scale estimated and p by Monte Carlo: shape only, weighted towards the tails, where a Student-t or mixture likelihood is meant to differ from the Normal.

Parameters:

z (np.ndarray) – Standardised residuals.
n_mc (int) – Monte Carlo samples for the Lilliefors and Anderson-Darling p-values.
seed (int) – Seed for those samples, so a report is reproducible.

Returns:

sw_w, sw_p, ks_n01_d, ks_n01_p, lillie_d, lillie_p, ad_a2, ad_p.

Return type:

dict[str, float]

clophfit.fitting.residual_tests.durbin_watson_null(table, params, *, is_ph=True)#

Durbin-Watson and lag-1 correlation expected of fitted residuals, per label.

With independent, calibrated errors the weighted residuals of a least-squares fit have covariance I - H, H the hat matrix of the well’s weighted Jacobian (both labels, one shared K). So, per series, E[sum (e_t - e_{t-1})^2] = tr(D'D (I - H)) and E[sum e_t^2] = tr(I - H) with D the first-difference operator along the steps; pooling those traces over wells gives the null DW of the pooled statistic, and the summed lag-1 covariances over the summed variances of the leading and trailing points the null lag-1 correlation. Plateau points carry the most leverage, so the series ends are the least variable, which is why the two differ. A multi-well Bayesian fit spends more than these parameters (the latent pH axis), so its true null lies further from the independent-error value.

Parameters:

table (pd.DataFrame) – Residual table with well, label, step (or x), x, raw_res and sigma.
params (Mapping[str, Mapping[str, float]]) – Well to fitted K/S0_<label>/S1_<label>.
is_ph (bool) – pH titration (else concentration).

Returns:

Indexed by label: dw_expected_fit and lag1_expected.

Return type:

pd.DataFrame

Examples

A straight line fitted to seven independent points: its residuals alternate, so the null is above the independent-error 12/7.

>>> import numpy as np, pandas as pd
>>> x = np.linspace(0, 1, 7)
>>> q, _ = np.linalg.qr(np.column_stack([np.ones(7), x]))
>>> s = np.eye(7) - q @ q.T
>>> d = np.diff(np.eye(7), axis=0)
>>> round(float(np.trace(d.T @ d @ s) / np.trace(s)), 2)
2.36

clophfit.fitting.residual_tests.runs_test(table, *, z_col='std_res')#

Wald-Wolfowitz runs test on residual signs, pooled within series.

A run is a stretch of residuals with one sign along a (well, label) series. Too few runs means neighbouring steps share their error, which is the wave Durbin-Watson also detects. The runs test uses signs only, so a few large residuals, which dominate DW’s sums of squares, cannot drive it. It suits Student-t and mixture fits, where such points are expected. Runs, their expectations and variances (conditional on each series’ sign counts) are summed over series, so short titrations still pool into a usable test. Like DW, the textbook null assumes independent errors. Fitted residuals alternate more, so a correct fit shows more runs than runs_expected. Use runs_test_null() for the fitted-residual null.

Parameters:

table (pd.DataFrame) – Residual table with well, label, step (or x) and z_col.
z_col (str) – The residual column.

Returns:

runs (observed), runs_expected and runs_z under independent errors, runs_p (two-sided, normal approximation) and n_series.

Return type:

dict[str, float]

Examples

Twenty series that each change sign once have 40 runs, far fewer than independent signs would give:

>>> import pandas as pd
>>> rows = [
...     {"well": w, "label": "1", "step": i, "std_res": 1.0 if i < 4 else -1.0}
...     for w in range(20)
...     for i in range(8)
... ]
>>> out = runs_test(pd.DataFrame(rows))
>>> out["runs"], out["runs_p"] < 1e-6
(40.0, True)

clophfit.fitting.residual_tests.runs_test_null(table, params, *, is_ph=True, n_mc=999, seed=0)#

Distribution of the pooled runs count for fitted residuals, per label.

With independent, calibrated errors, a least-squares fit leaves weighted residuals e = (I - H) u with u ~ N(0, I) and H the hat matrix of the well’s weighted Jacobian (both labels, one shared K), as in durbin_watson_null(). The runs count of that projection has no closed form, so it is simulated. Each Monte Carlo draw projects fresh noise through every well and sums the runs per label, which gives the null distribution of the statistic runs_test() reports.

Parameters:

table (pd.DataFrame) – Residual table with well, label, step (or x), x, raw_res and sigma.
params (Mapping[str, Mapping[str, float]]) – Well to fitted K/S0_<label>/S1_<label>.
is_ph (bool) – pH titration (else concentration).
n_mc (int) – Monte Carlo draws.
seed (int) – Seed for the draws.

Returns:

Indexed by label: runs_expected_fit and runs_sd_fit (mean and SD of the simulated pooled runs), and runs_null (the simulated counts, as an array, for the p-value in residual_tests()).

Return type:

pd.DataFrame

clophfit.fitting.residual_tests.well_leverage(table, params, *, is_ph=True)#

Leverage h_ii of every residual, from its well’s weighted binding fit.

Parameters:

table (pd.DataFrame) – Residual table with well, label, x, raw_res, sigma.
params (Mapping[str, Mapping[str, float]]) – Well to its fitted K, S0_<label>, S1_<label>.
is_ph (bool) – pH titration (else concentration).

Returns:

h_ii aligned with table’s index; NaN for wells without parameters.

Return type:

pd.Series

clophfit.fitting.residual_tests.studentized_outliers(table, params, *, is_ph=True, alpha=_OUTLIER_ALPHA)#

Externally studentized residual of every point, and whether it is an outlier.

Each well is its own regression: t_i follows Student-t with n - p - 1 degrees of freedom, and the cutoff is Bonferroni-corrected for the well’s n points, so a plate of clean wells flags ~alpha of its wells, not of its points.

Parameters:

table (pd.DataFrame) – Residual table with well, label, x, raw_res, sigma.
params (Mapping[str, Mapping[str, float]]) – Well to its fitted parameters.
is_ph (bool) – pH titration (else concentration).
alpha (float) – Family-wise level per well.

Returns:

table’s index with t_abs, dof, cutoff and outlier.

Return type:

pd.DataFrame

clophfit.fitting.residual_tests.press_statistics(table, params, *, is_ph=True)#

PRESS: the prediction error of every point when it is left out of its well’s fit.

e_(i) = e_i / (1 - h_ii) on each well’s linearised least-squares residuals (see well_leverage()), exact for a linear fit and first-order for a nonlinear one, so no refits are needed. Reported on the raw scale, in signal counts, which makes it comparable across noise models fitted to the same points: a lower PRESS is a mean model that predicts better.

Parameters:

table (pd.DataFrame) – Residual table with well, label, x, raw_res, sigma.
params (Mapping[str, Mapping[str, float]]) – Well to its fitted parameters.
is_ph (bool) – pH titration (else concentration).

Returns:

press (sum of squared deleted residuals), press_rms (its root mean, in counts), press_ratio (PRESS / SSE; 1 would mean no point pulls its own fit) and max_leverage.

Return type:

dict[str, float]

clophfit.fitting.residual_tests.residual_tests(table, params=None, *, is_ph=True, n_mc=999, seed=0)#

All four assumption checks, one row per label.

Parameters:

table (pd.DataFrame) – Canonical residual table: well, label, step or x, yhat, raw_res, sigma, std_res.
params (Mapping[str, Mapping[str, float]] | None) – Well to fitted K/S0_<label>/S1_<label>. Without it the leverage-based columns (studentized outliers, PRESS) are left NaN.
is_ph (bool) – pH titration (else concentration).
n_mc (int) – Monte Carlo samples for the Lilliefors and Anderson-Darling p-values and the fitted-residual runs null.
seed (int) – Seed for those samples.

Returns:

Indexed by label. Columns: n; bp_r2/bp_p and white_r2/white_p (homoscedasticity); dw, dw_expected, lag1_r, dw_p, and with params the fitted-residual nulls dw_expected_fit and lag1_expected; runs, runs_expected, runs_z, runs_p, and with params runs_expected_fit, runs_sd_fit and the Monte Carlo runs_p_fit (independence); sw_w, sw_p, ks_n01_d, ks_n01_p, lillie_d, lillie_p, ad_a2, ad_p (normality); t_max, outlier_rate, wells_with_outlier (studentized outliers); press, press_rms, press_ratio, max_leverage (PRESS).

Return type:

pd.DataFrame