clophfit.fitting.residual_tests
===============================

.. py:module:: clophfit.fitting.residual_tests

.. autoapi-nested-parse::

   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 :func:`~clophfit.fitting.models.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
---------

.. autoapisummary::

   clophfit.fitting.residual_tests.breusch_pagan
   clophfit.fitting.residual_tests.white_test
   clophfit.fitting.residual_tests.durbin_watson
   clophfit.fitting.residual_tests.normality_tests
   clophfit.fitting.residual_tests.durbin_watson_null
   clophfit.fitting.residual_tests.runs_test
   clophfit.fitting.residual_tests.runs_test_null
   clophfit.fitting.residual_tests.well_leverage
   clophfit.fitting.residual_tests.studentized_outliers
   clophfit.fitting.residual_tests.press_statistics
   clophfit.fitting.residual_tests.residual_tests


Module Contents
---------------

.. py:function:: 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.

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

   :returns: ``lm``, ``df`` (2), ``p`` and ``r2`` - the fraction of the squared
             residuals' variance the fan explains, the effect size to compare.
   :rtype: dict[str, float]

   .. rubric:: 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


.. py:function:: 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.

   :param z: Standardised residuals.
   :type z: np.ndarray
   :param yhat: Predictions.
   :type yhat: np.ndarray
   :param x: Titration coordinate.
   :type x: np.ndarray

   :returns: ``lm``, ``df`` (5), ``p`` and ``r2``.
   :rtype: dict[str, float]


.. py:function:: 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
   :func:`durbin_watson_null` for the null of the fitted residuals.

   :param table: Residual table with ``well``, ``label``, ``step`` (or ``x``) and *z_col*.
   :type table: pd.DataFrame
   :param z_col: The residual column.
   :type z_col: str

   :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.
   :rtype: dict[str, float]


.. py:function:: 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.

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

   :returns: ``sw_w``, ``sw_p``, ``ks_n01_d``, ``ks_n01_p``, ``lillie_d``, ``lillie_p``,
             ``ad_a2``, ``ad_p``.
   :rtype: dict[str, float]


.. py:function:: 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.

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

   :returns: Indexed by label: ``dw_expected_fit`` and ``lag1_expected``.
   :rtype: pd.DataFrame

   .. rubric:: 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


.. py:function:: 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 :func:`runs_test_null` for the fitted-residual null.

   :param table: Residual table with ``well``, ``label``, ``step`` (or ``x``) and *z_col*.
   :type table: pd.DataFrame
   :param z_col: The residual column.
   :type z_col: str

   :returns: ``runs`` (observed), ``runs_expected`` and ``runs_z`` under independent
             errors, ``runs_p`` (two-sided, normal approximation) and ``n_series``.
   :rtype: dict[str, float]

   .. rubric:: 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)


.. py:function:: 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
   :func:`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 :func:`runs_test` reports.

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

   :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 :func:`residual_tests`).
   :rtype: pd.DataFrame


.. py:function:: well_leverage(table, params, *, is_ph = True)

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

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

   :returns: ``h_ii`` aligned with *table*'s index; NaN for wells without parameters.
   :rtype: pd.Series


.. py:function:: 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.

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

   :returns: *table*'s index with ``t_abs``, ``dof``, ``cutoff`` and ``outlier``.
   :rtype: pd.DataFrame


.. py:function:: 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 :func:`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.

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

   :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``.
   :rtype: dict[str, float]


.. py:function:: residual_tests(table, params = None, *, is_ph = True, n_mc = 999, seed = 0)

   All four assumption checks, one row per label.

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

   :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).
   :rtype: pd.DataFrame