from __future__ import annotations
import warnings
from collections import defaultdict
from typing import TYPE_CHECKING, Any, Literal, NoReturn
import numpy as np
import pandas as pd
from scipy.stats import shapiro, t
from sklearn.utils import Bunch
if TYPE_CHECKING:
from collections.abc import Callable
__eps = np.finfo(np.float32).eps
[docs]
def t_value(p: float, v: float) -> float:
r"""
Return the value on the x-axis if you plot the cumulative t-distribution with a fractional
area of `p` (p is therefore a fractional value between 0 and 1 on the y-axis) and `v` is the
degrees of freedom.
Examples
--------
Since the cumulative distribution passes symmetrically through the x-axis at 0.0 for any
number of degrees of freedom
>>> t_value(0.5, v=10)
0.0
Zero fractional area under the curve is always at :math:`-\infty`:
>>> t_value(0.0, v=10)
-inf
100% fractional area is always at :math:`+\infty`:
>>> t_value(1.0, v=10)
inf
See also
--------
t_value_cdf: does the inverse of this function.
"""
return t.ppf(p, df=v)
[docs]
def t_value_cdf(z: float, v: float) -> float:
r"""
Return the fractional area under the cumulative t-distribution (y-axis value) at the t-value
`z` on the x-axis, with `v` degrees of freedom.
Examples
--------
The cumulative distribution is symmetric through the x-axis at 0.0 for any number of degrees
of freedom, so half of the area lies below zero:
>>> t_value_cdf(0.0, v=10)
0.5
Zero fractional area under the curve is at :math:`-\infty`:
>>> t_value_cdf(-np.inf, v=10)
0.0
100% fractional area is at :math:`+\infty`:
>>> t_value_cdf(np.inf, v=10)
1.0
See also
--------
t_value: does the inverse of this function.
"""
return t.cdf(z, df=v)
[docs]
def test_normality(x: np.ndarray | pd.Series) -> float:
"""
Check the p-value of the hypothesis that the data are from a normal distribution.
If the p-value is less than the chosen alpha level (e.g. 0.05 or 0.025), then there is
evidence that the data tested are NOT normally distributed.
On the other hand, if the p-value is greater than the chosen alpha level, then the null
hypothesis that the data came from a normally distributed population can not be rejected.
NOTE: it does not mean that the data are normally distributed, just that we have nothing better
to say about it. See the `Shapiro-Wilk test`_.
Implementation: Uses the Shapiro Wilk test directly taken from `scipy.stats.shapiro`_.
.. _Shapiro-Wilk test:
https://en.wikipedia.org/wiki/Shapiro-Wilk_test
.. _scipy.stats.shapiro:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html
"""
output = shapiro(x)
return output[1]
[docs]
def Sn(x: np.ndarray | pd.Series, constant: float = 1.1926) -> np.floating: # noqa: N802
"""
Compute a robust scale estimator. The Sn metric is an efficient alternative to MAD.
Parameters
----------
x : np.ndarray or pd.Series
A vector of values. NaN entries are ignored.
constant : float, optional
Multiplicative constant that makes the estimator consistent with iid values from a
Gaussian distribution with no outliers. Default is 1.1926.
Returns
-------
np.floating
A scalar value, the Sn estimate of spread. Returns NaN if all entries are missing,
and 0.0 when only a single non-missing value is supplied.
Notes
-----
Follows the Rousseeuw-Croux definition [2]_ exactly, including the high median
inside and the low median outside, and reproduces ``robustbase::Sn`` [1]_ (with
``finite.corr=TRUE``) for both odd and even sample sizes.
Disadvantages of MAD:
* It does not have particularly high efficiency for data that is in fact normal (37%).
In comparison, the median has 64% efficiency for normal data.
* The MAD statistic also has an implicit assumption of symmetry. That is, it measures the
distance from a measure of central location (the median).
References
----------
.. [1] https://cran.r-project.org/web/packages/robustbase/
.. [2] Rousseeuw, Peter J.; Croux, Christophe (December 1993),
"Alternatives to the Median Absolute Deviation", Journal of the American Statistical
Association, American Statistical Association, 88 (424): 1273-1283,
https://dx.doi.org/10.2307/2291267
"""
arr = np.asarray(x, dtype=np.float64)
n = np.sum(~np.isnan(arr))
if n == 0:
return np.float64(np.nan)
elif n == 1:
return np.float64(0.0)
# Remove NaNs and compute all pairwise absolute differences via broadcasting
clean = arr[~np.isnan(arr)]
diffs = np.abs(clean[:, None] - clean[None, :])
# Rousseeuw-Croux define Sn = c * lomed_i { himed_j |x_i - x_j| }, using the
# HIGH median inside and the LOW median outside. These are order statistics,
# not the averaging median numpy provides: himed is element floor(n/2) + 1
# and lomed is element floor((n + 1) / 2), one-based, of the sorted values.
# They coincide with np.median only for odd n, so using np.median for both
# (the previous behaviour) biased every even-n estimate downward: 0.443 in
# place of 0.886 at n = 2, and still 0.5 percent low at n = 1500.
inner_high_median = np.sort(diffs, axis=1)[:, n // 2]
outer_low_median = np.sort(inner_high_median)[(n + 1) // 2 - 1]
if n <= 9:
# Correction factors for n = 2 to 9:
correction = [0.743, 1.851, 0.954, 1.351, 0.993, 1.198, 1.005, 1.131][n - 2]
elif n % 2:
correction = n / (n - 0.9) # odd number of cases
else:
correction = 1.0
return constant * outer_low_median * correction
def _contains_nan(a: np.ndarray, nan_policy: str = "propagate") -> tuple[bool, str]:
"""From scipy.stats.stats."""
policies = ["propagate", "raise", "omit"]
if nan_policy not in policies:
raise ValueError(f"nan_policy must be one of {{{', '.join(f'{s!r}' for s in policies)}}}")
try:
# Calling np.sum to avoid creating a huge array into memory
# e.g. np.isnan(a).any()
with np.errstate(invalid="ignore"):
contains_nan = np.isnan(np.sum(a))
except TypeError:
# This can happen when attempting to sum things which are not
# numbers (e.g. as in the function `mode`). Try an alternative method:
try:
contains_nan = any(np.isnan(v) for v in a.ravel())
except TypeError: # pragma: no cover
# Don't know what to do. Fall back to omitting nan values and
# issue a warning.
contains_nan = False
nan_policy = "omit"
warnings.warn(
"The input array could not be properly checked for nan values. nan values will be ignored.",
RuntimeWarning,
stacklevel=2,
)
if contains_nan and nan_policy == "raise":
raise ValueError("The input contains nan values")
return (contains_nan, nan_policy)
[docs]
def ttest_independent(
sample_A: pd.Series,
sample_B: pd.Series,
conflevel: float = 0.995,
*,
equal_var: bool = True,
) -> dict:
"""Core calculation for a test of differences between the average of A and the average of B.
No checking of inputs.
Parameters
----------
sample_A : iterable
Vector of n_a measurements.
sample_B : iterable
Vector of n_b measurements.
conflevel : float
Value between 0 and 1 (closer to 1.0), that gives the level of confidence required for
the 2-sided test.
equal_var : bool
``True`` (the default) runs Student's pooled-variance test: the two samples are
assumed to share one variance, which is estimated from both and spent on
``n_A + n_B - 2`` degrees of freedom. ``False`` runs Welch's test: each sample
keeps its own variance and the degrees of freedom come from the
Welch-Satterthwaite approximation. Matches
:func:`scipy.stats.ttest_ind`'s parameter of the same name.
Returns
-------
dict
Outcomes from the statistical test. The keys of interest are
``"t value"`` (the test statistic; its ``"p value"`` is computed from
the t distribution on ``"Degrees of freedom"``), ``"Std error of
difference"`` (the denominator of that statistic, and the half-width
unit of the confidence interval), and ``"Pooled std dev"``.
``"Equal variance assumed"`` records which test was run;
``"Pooled std dev"`` is ``NaN`` under Welch, where no pooled estimate exists.
.. deprecated:: 1.70.0
``"z value"`` and ``"Pooled standard deviation"`` are misleading
names kept as aliases, and will be removed in 2.0. The statistic
is a t-statistic, not a z-statistic, and ``"Pooled standard
deviation"`` holds the standard *error* of the difference in
means, ``sqrt(svar * (1/nA + 1/nB))``, not the pooled standard
deviation ``sqrt(svar)``. Use ``"t value"`` and ``"Std error of
difference"`` instead.
Notes
-----
The default stays Student's so existing results do not move, but it is the weaker
choice. Pooling is only valid when the two variances really are equal; when they are
not, and especially when the larger variance sits with the smaller sample, Student's
test does not hold its nominal error rate, while Welch's does and costs almost
nothing when the variances *are* equal. Delacre, Lakens & Leys (2017) make the case
that Welch's should be the default, and it already is in R's ``t.test``:
https://rips-irsp.com/articles/10.5334/irsp.82
Testing the variances first and choosing on the outcome is worse than either, not
better: the pre-test spends its own error rate, and the conditional procedure has no
clean size. Pick ``equal_var`` from what is known about the measurement, not from the
data in hand.
"""
axis: Literal[0] = 0
v1, v2 = sample_A.var(axis=axis, ddof=1), sample_B.var(axis=axis, ddof=1)
n_A, n_B = sample_A.shape[axis], sample_B.shape[axis]
m1, m2 = sample_A.mean(), sample_B.mean()
d = m2 - m1
with np.errstate(divide="ignore", invalid="ignore"):
if equal_var:
df = n_A + n_B - 2.0
svar = ((n_A - 1) * v1 + (n_B - 1) * v2) / df
sd_z_variate = np.sqrt(svar * (np.divide(1.0, n_A) + np.divide(1.0, n_B)))
pooled_std_dev = np.sqrt(svar)
else:
# Welch: each sample contributes its own squared standard error, and the
# degrees of freedom are the Welch-Satterthwaite approximation, which
# interpolates between min(n_A, n_B) - 1 and n_A + n_B - 2 according to how
# unbalanced the two contributions are.
se2_A, se2_B = np.divide(v1, n_A), np.divide(v2, n_B)
sd_z_variate = np.sqrt(se2_A + se2_B)
df = np.divide(
(se2_A + se2_B) ** 2,
np.divide(se2_A**2, n_A - 1) + np.divide(se2_B**2, n_B - 1),
)
# No single variance is estimated, so there is no pooled standard deviation
# to report. NaN says that, where a number would imply a quantity the test
# never formed.
pooled_std_dev = np.nan
z_variate = np.divide(d, sd_z_variate)
ct = abs(t_value(p=(1 - conflevel) / 2.0, v=df))
confint_lo = d - ct * sd_z_variate
confint_hi = d + ct * sd_z_variate
return {
"Group A number": int(n_A),
"Group B number": int(n_B),
"Group A average": sample_A.mean(),
"Group B average": sample_B.mean(),
"z value": z_variate,
"ConfInt: Lo": confint_lo,
"ConfInt: Hi": confint_hi,
"p value": 2 * t_value_cdf(-np.abs(z_variate), df),
"Degrees of freedom": df,
"Equal variance assumed": bool(equal_var),
"Pooled standard deviation": sd_z_variate,
# Correctly named entries for the two mislabelled ones above, which are
# kept as deprecated aliases. "z value" is a t-statistic: its p value is
# computed from the t distribution on `df` degrees of freedom. "Pooled
# standard deviation" is the standard ERROR of the difference in means,
# sqrt(svar * (1/nA + 1/nB)); the pooled standard deviation is sqrt(svar).
"t value": z_variate,
"Std error of difference": sd_z_variate,
"Pooled std dev": pooled_std_dev,
}
def _apply_multiplicity_correction(output: pd.DataFrame, correction: str | None) -> pd.DataFrame:
"""Add corrected p-values to a table of pairwise comparisons.
A family of pairwise comparisons inflates the family-wise error rate: with
k groups there are k(k-1)/2 tests, so at k=5 ten raw p-values are compared
against 0.05 and the chance of at least one false positive is far above 5
percent. The raw ``p value`` column is always kept as-is; the corrected
values are added alongside so both are visible.
"""
if correction is None or output.empty:
return output
known = {"holm": holm_bonferroni, "bh": benjamini_hochberg, "benjamini-hochberg": benjamini_hochberg}
key = correction.strip().lower()
if key not in known:
raise ValueError(f"correction must be one of 'holm', 'bh', or None; got {correction!r}.")
adjusted = known[key](output["p value"].to_numpy())
output = output.copy()
output["p value (adjusted)"] = adjusted.p_adjusted
output["reject"] = adjusted.reject
output["correction"] = key
return output
[docs]
def ttest_independent_from_df( # noqa: PLR0913 - six is the honest width: two columns, a frame, and three switches
df: pd.DataFrame,
grouper_column: str,
values_column: str,
conflevel: float = 0.995,
correction: str | None = None,
*,
equal_var: bool = True,
) -> pd.DataFrame:
"""
Calculate the t-test for differences between two or more groups and returns a confidence
interval for the difference. The test is for UNPAIRED differences.
The dataframe `df` contains a `grouper_column` with 2 or more unique values (e.g. 'A' and 'B').
All unique values of the `grouper_column` are used, and t-tests are done between the values
in the `values_column`.
Args:
df (pd.DataFrame): Dataframe of the values and grouping variable.
grouper_column (str): Indicates which column will be grouped on.
values_column (str): Which column contains the numeric values to calculate the test on.
conflevel (float, optional): [description]. Defaults to 0.995.
correction (str | None, optional): "holm", "bh", or None. Defaults to None. See below.
equal_var (bool, optional): Forwarded to :func:`ttest_independent`. Defaults to True.
Equal variances: ``equal_var=True`` (the default) is Student's pooled-variance test;
``False`` is Welch's. Welch is the safer choice for a pairwise family, where the groups
have no particular reason to share one variance; see :func:`ttest_independent`'s Notes.
Multiplicity: by default the returned p-values are UNCORRECTED. With k groups this
runs k(k-1)/2 tests, so the chance of at least one false positive is well above the
per-test level. Pass ``correction="holm"`` to control the family-wise error rate (see
:func:`holm_bonferroni`) or ``correction="bh"`` to control the false discovery rate
(see :func:`benjamini_hochberg`). Either adds ``p value (adjusted)``, ``reject`` and
``correction`` columns; the raw ``p value`` column is always kept.
Output: Dataframe with columns containing the statistical outputs of the t-test, including:
1. Group "A" name
2. Group "B" name
3. Group "A" mean
4. Group "B" mean
5. z-value for the difference between group "B" minus group "A"
6. p-value for this z-value
7. Confidence interval low value for difference between group "B" minus group "A"
8. Confidence interval high value for difference between group "B" minus group "A"
Example:
df : has 3 levels in the grouper variable; ['Marco', 'Pete', 'Sam']
Output will have 3 rows:
Group A name Group B name
Marco Pete
Marco Sam
Pete Sam
"""
data_subset = df[[grouper_column, values_column]].copy()
data_subset = data_subset.dropna()
output = pd.DataFrame()
# Enumerate the groups from the CLEANED data. Taking them from `df` (the
# previous behaviour) kept groups whose values are all missing, and NaN
# group labels, each of which then produced a silent all-NaN comparison
# row against an empty sample.
groups = list(data_subset[grouper_column].unique())
while len(groups) > 0:
groupA_name = groups.pop(0)
for groupB_name in groups:
sample_A = data_subset[data_subset[grouper_column].eq(groupA_name)][values_column]
sample_B = data_subset[data_subset[grouper_column].eq(groupB_name)][values_column]
sample_A = sample_A.astype(np.float64)
sample_B = sample_B.astype(np.float64)
basic_stats = ttest_independent(sample_A, sample_B, conflevel, equal_var=equal_var)
basic_stats.update(
{
"Group A name": groupA_name,
"Group B name": groupB_name,
}
)
output = pd.concat([output, pd.DataFrame(basic_stats, index=[0])])
return _apply_multiplicity_correction(output, correction)
[docs]
def ttest_paired(differences: pd.Series, conflevel: float = 0.995) -> dict:
"""Core calculation for a test of paired differences.
Parameters
----------
differences : pd.Series
The paired differences (e.g. ``sample_A - sample_B``) for which the test is run.
conflevel : float, optional
Value between 0 and 1 (closer to 1.0), that gives the level of confidence required for
the 2-sided test. Default is 0.995.
Returns
-------
dict
Outcomes from the statistical test, with keys:
- ``"Differences mean"``: the mean of the input ``differences``.
- ``"z value"``: the test statistic (mean divided by its standard error).
- ``"ConfInt: Lo"``, ``"ConfInt: Hi"``: lower and upper bounds of the two-sided
confidence interval around the mean difference at the requested ``conflevel``.
- ``"p value"``: two-sided p-value from the t-distribution.
- ``"Degrees of freedom"``: ``n - 1`` for ``n`` paired observations.
- ``"Standard deviation"``: the **standard error of the mean difference**,
i.e. ``sample_std / sqrt(n)`` (the scale used to build the confidence interval),
NOT the sample standard deviation of ``differences``.
"""
# A paired t-test needs at least 2 differences (so dof >= 1 and the
# 1/n term inside sqrt is finite). SEC-24 (#273).
if differences.shape[0] < 2:
raise ValueError(f"ttest_paired requires at least 2 paired observations; got {differences.shape[0]}.")
diff_mean = differences.mean()
diff_svar = differences.std(ddof=1)
dof = differences.shape[0] - 1 # n-1 d
# By the central limit theorem, the `differences` values should be normally distributed
# with average (central value) given by mean of group A values, minus mean of group B values.
# Scale factor is the standard deviations of `differences`: estimated, therefore t-distribution
ct = abs(t_value(p=(1 - conflevel) / 2.0, v=dof))
with np.errstate(divide="ignore", invalid="ignore"):
sd_z_variate = diff_svar * np.sqrt(np.divide(1.0, dof + 1))
z_variate = np.divide(diff_mean, sd_z_variate)
confint_lo = diff_mean - ct * sd_z_variate
confint_hi = diff_mean + ct * sd_z_variate
return {
"Differences mean": diff_mean,
"z value": z_variate,
"ConfInt: Lo": confint_lo,
"ConfInt: Hi": confint_hi,
"p value": 2 * t_value_cdf(-np.abs(z_variate), dof),
"Degrees of freedom": dof,
"Standard deviation": sd_z_variate,
}
[docs]
def ttest_paired_from_df(
df: pd.DataFrame,
grouper_column: str,
values_column: str,
conflevel: float = 0.995,
correction: str | None = None,
) -> pd.DataFrame:
"""
Calculate the t-test for paired differences between two or more groups and returns a
confidence interval for the difference. The test is for PAIRED differences.
The differences is always defined as the A values minus the B values: after - before, or A - B.
The dataframe `df` contains a `grouper_column` with 2 or more unique values (e.g. 'A' and 'B').
All unique values of the `grouper_column` are used, and t-tests are done between the values
in the `values_column`.
When selecting the columns, the number of values per column must be the same.
Args:
df (pd.DataFrame): Dataframe of the values and grouping variable.
grouper_column (str): Indicates which column will be grouped on.
values_column (str): Which column contains the numeric values to calculate the test on.
conflevel (float, optional): [description]. Defaults to 0.995.
correction (str | None, optional): "holm", "bh", or None. Defaults to None. See below.
Multiplicity: by default the returned p-values are UNCORRECTED. With k groups this
runs k(k-1)/2 tests, so the chance of at least one false positive is well above the
per-test level. Pass ``correction="holm"`` to control the family-wise error rate (see
:func:`holm_bonferroni`) or ``correction="bh"`` to control the false discovery rate
(see :func:`benjamini_hochberg`). Either adds ``p value (adjusted)``, ``reject`` and
``correction`` columns; the raw ``p value`` column is always kept.
Output: Dataframe with columns containing the statistical outputs of the t-test, including:
1. Group A name
2. Group B name
3. Group A mean
4. Group B mean
5. Differences mean: average difference between the groups (not the same as the difference
of the averages from items 3 and 4 above)
6. z-value for the difference between group "A" minus group "B"
7. p-value for this z-value
8. Confidence interval low value for difference between group "A" minus group "B"
9. Confidence interval high value for difference between group "A" minus group "B"
"""
data_subset = df[[grouper_column, values_column]].copy()
data_subset = data_subset.dropna()
output = pd.DataFrame()
# Enumerate the groups from the CLEANED data. Taking them from `df` (the
# previous behaviour) kept groups whose values are all missing, and NaN
# group labels, each of which then produced a silent all-NaN comparison
# row against an empty sample.
groups = list(data_subset[grouper_column].unique())
while len(groups) > 0:
groupA_name = groups.pop(0)
for groupB_name in groups:
sample_A = data_subset[data_subset[grouper_column].eq(groupA_name)][values_column]
sample_B = data_subset[data_subset[grouper_column].eq(groupB_name)][values_column]
sample_A = sample_A.astype(np.float64)
sample_B = sample_B.astype(np.float64)
if sample_A.shape[0] != sample_B.shape[0]:
raise ValueError(
"Paired t-test requires both groups to have the same number of samples; "
f"got {sample_A.shape[0]} and {sample_B.shape[0]}."
)
differences = sample_A - sample_B.to_numpy() # only the values of one vector are needed!
basic_stats = ttest_paired(differences, conflevel)
basic_stats.update(
{
"Group A name": groupA_name,
"Group B name": groupB_name,
"Group A number": sample_A.shape[0],
"Group B number": sample_B.shape[0],
"Group A average": sample_A.mean(),
"Group B average": sample_B.mean(),
}
)
output = pd.concat([output, pd.DataFrame(basic_stats, index=[0])])
return _apply_multiplicity_correction(output, correction)
[docs]
def confidence_interval(df: pd.DataFrame, column_name: str, conflevel: float = 0.95, style: str = "robust") -> tuple:
"""
Calculate the confidence interval, returned as a tuple, for the `column_name` (str) in the
dataframe `df`, for a given confidence level `conflevel` (default: 0.95).
`style`: ['robust'; 'regular']: indicates which style of estimates to use for the center and
spread. Default: 'robust'
Missing values are ignored.
Raises
------
ValueError
If `style` is neither 'robust' nor 'regular', or if the column holds fewer than 2
non-missing values. An unrecognised `style` used to fall through to the classical
branch, so a typo returned a different interval with nothing to signal it (#561).
"""
known_styles = ("robust", "regular")
if style.lower() not in known_styles:
raise ValueError(f"style must be one of {known_styles}; got {style!r}.")
data = df[column_name]
n = data.count()
# A t-CI needs at least 2 non-missing observations to compute a
# spread and (n - 1) > 0 degrees of freedom. ``t.ppf(., -1)`` /
# ``spread / sqrt(0)`` silently yielded NaN / inf in earlier
# versions; reject up front. SEC-24 (#273).
if n < 2:
raise ValueError(
f"confidence_interval requires at least 2 non-missing values in column {column_name!r}; got {n}."
)
if style.lower() == "robust":
center = data.median()
# MAD (scale="normal") consistently estimates sigma, but the interval
# is for the MEDIAN, whose asymptotic standard error is
# sigma * sqrt(pi / 2) / sqrt(n), not sigma / sqrt(n). Without the
# sqrt(pi/2) factor the advertised 95% interval has roughly 87%
# coverage.
spread = median_absolute_deviation(data.to_numpy(), nan_policy="omit") * float(np.sqrt(np.pi / 2.0))
else:
center = data.mean()
spread = data.std()
c_t = t_value(1 - (1 - conflevel) / 2, n - 1)
return (center - c_t * spread / np.sqrt(n), center + c_t * spread / np.sqrt(n))
def _mad_1d(x: np.ndarray, center: Callable, nan_policy: str) -> float:
"""Taken from `scipy.stats.stats`.
Median absolute deviation for 1-d array x.
This is a helper function for `median_abs_deviation`; it assumes its
arguments have been validated already. In particular, x must be a
1-d numpy array, center must be callable, and if nan_policy is not
'propagate', it is assumed to be 'omit', because 'raise' is handled
in `median_abs_deviation`.
No warning is generated if x is empty or all nan.
"""
isnan = np.isnan(x)
if isnan.any():
if nan_policy == "propagate":
return np.nan
x = x[~isnan]
if x.size == 0:
# MAD of an empty array is nan.
return np.nan
# Edge cases have been handled, so do the basic MAD calculation.
med = center(x)
return np.median(np.abs(x - med))
[docs]
def biweight_midvariance(x: np.ndarray | pd.Series, nan_policy: str = "omit") -> float:
"""Return the Mosteller-Tukey robust scale (biweight midvariance) of ``x``.
The biweight midvariance is a robust, highly efficient estimator of the
variance: it down-weights observations far from the median and ignores
gross outliers entirely.
Parameters
----------
x : np.ndarray or pd.Series
One-dimensional sample of numeric values.
nan_policy : {"omit", "propagate"}, optional
``"omit"`` (default) drops missing values; ``"propagate"`` returns
``nan`` if any value is missing.
Returns
-------
float
The robust variance estimate. Returns ``0.0`` when the MAD is zero
(e.g. a constant sample), and ``nan`` for an empty sample.
References
----------
Mosteller and Tukey, *Data Analysis and Regression*, pp. 207-208, 1977.
"""
a = np.asarray(x, dtype=float).ravel()
isnan = np.isnan(a)
if isnan.any():
if nan_policy == "propagate":
return float("nan")
a = a[~isnan]
n = a.size
if n == 0:
return float("nan")
location = np.median(a)
spread_mad = np.median(np.abs(a - location))
if spread_mad == 0:
return 0.0
# c = 9 is the biweight MIDVARIANCE tuning constant (Mosteller & Tukey;
# also NIST DATAPLOT and astropy); c = 6, previously used here, is the
# biweight LOCATION constant and rejects far too much of the sample for
# a scale estimate to be consistent with sigma^2 on Gaussian data.
ui = (a - location) / (9.0 * spread_mad)
valid = ui**2 <= 1.0
a_valid, ui_valid = a[valid], ui[valid]
numerator = (a_valid - location) ** 2 * (1 - ui_valid**2) ** 4
denominator = (1 - ui_valid**2) * (1 - 5 * ui_valid**2)
return float(n * numerator.sum() / denominator.sum() ** 2)
[docs]
def holm_bonferroni(p_values: np.ndarray | pd.Series | list, alpha: float = 0.05) -> Bunch:
"""Holm-Bonferroni step-down correction for multiple comparisons.
Holm's method controls the family-wise error rate while being uniformly
more powerful than the plain Bonferroni correction. It is the recommended
post-hoc correction for a family of pairwise comparisons.
Parameters
----------
p_values : array-like
The raw (uncorrected) p-values of the individual comparisons.
alpha : float, optional
Family-wise significance level, by default 0.05.
Returns
-------
sklearn.utils.Bunch
A bunch with, in the same order as the input:
* ``p_adjusted``: the Holm-adjusted p-values.
* ``reject``: boolean array, ``True`` where the null hypothesis is
rejected at level ``alpha``.
* ``alpha``: the family-wise level used.
References
----------
Holm, "A simple sequentially rejective multiple test procedure",
Scandinavian Journal of Statistics, 6, 65-70, 1979.
"""
p = np.asarray(p_values, dtype=float).ravel()
m = p.size
if m == 0:
return Bunch(p_adjusted=np.array([]), reject=np.array([], dtype=bool), alpha=alpha)
order = np.argsort(p)
p_sorted = p[order]
# Step-down weights m, m-1, ..., 1; the running max enforces monotonicity.
weights = m - np.arange(m)
adjusted_sorted = np.minimum(np.maximum.accumulate(weights * p_sorted), 1.0)
p_adjusted = np.empty(m)
p_adjusted[order] = adjusted_sorted
return Bunch(p_adjusted=p_adjusted, reject=p_adjusted <= alpha, alpha=alpha)
[docs]
def benjamini_hochberg(p_values: np.ndarray | pd.Series | list, alpha: float = 0.05) -> Bunch:
"""Benjamini-Hochberg step-up correction controlling the false discovery rate.
The Benjamini-Hochberg (BH) procedure controls the expected proportion of
false positives among the rejected hypotheses (the false discovery rate),
rather than the family-wise error rate that :func:`holm_bonferroni`
controls. BH is the conventional choice when a family contains many
comparisons and a few false positives are tolerable, for example testing
one product effect across many sensory attributes.
Parameters
----------
p_values : array-like
The raw (uncorrected) p-values of the individual comparisons.
alpha : float, optional
Target false-discovery rate, by default 0.05.
Returns
-------
sklearn.utils.Bunch
A bunch with, in the same order as the input:
* ``p_adjusted``: the BH-adjusted p-values (q-values), monotone in the
rank order of the input p-values.
* ``reject``: boolean array, ``True`` where the hypothesis is rejected
at the target false-discovery rate ``alpha``.
* ``alpha``: the target false-discovery rate used.
References
----------
Benjamini and Hochberg, "Controlling the false discovery rate: a practical
and powerful approach to multiple testing", Journal of the Royal
Statistical Society B, 57, 289-300, 1995.
"""
p = np.asarray(p_values, dtype=float).ravel()
m = p.size
if m == 0:
return Bunch(p_adjusted=np.array([]), reject=np.array([], dtype=bool), alpha=alpha)
order = np.argsort(p)
p_sorted = p[order]
ranks = np.arange(1, m + 1)
# Step-up: q_(i) = min over k >= i of (m / k) * p_(k); the reverse running
# minimum enforces the required monotonicity.
scaled = (m / ranks) * p_sorted
adjusted_sorted = np.minimum(np.minimum.accumulate(scaled[::-1])[::-1], 1.0)
p_adjusted = np.empty(m)
p_adjusted[order] = adjusted_sorted
return Bunch(p_adjusted=p_adjusted, reject=p_adjusted <= alpha, alpha=alpha)
def _relative_spread(spread: float, center: float) -> float:
"""Return ``spread / center``, or NaN when the centre is (near) zero.
A relative standard deviation is undefined at a zero centre: the division
gives ``inf`` (or NaN for 0/0) and emits a RuntimeWarning through the
public API and the agent-facing tools. Report NaN instead, which is what
"undefined" means to every downstream consumer.
"""
if not np.isfinite(center) or abs(center) <= __eps:
return float("nan")
return float(spread / center)
[docs]
def summary_stats(x: np.ndarray | pd.Series, method: str = "robust") -> dict:
"""
Return summary statistics of the numeric values in vector ``x``.
Parameters
----------
x : numpy.ndarray or pandas.Series
A vector of univariate values to summarize.
method : str, optional
If ``"robust"`` (the default), the reported center is the median and
the spread is the Sn robust estimate; otherwise the mean and the
sample standard deviation are used.
Returns
-------
dict
A summary of the univariate vector. The most useful keys are
``"center"`` (a measure of the center, e.g. the median for the robust
method) and ``"spread"`` (a measure of the spread, e.g. the Sn robust
estimate for the robust method).
"""
if isinstance(x, pd.Series):
values = x.copy(deep=True).to_numpy()
elif isinstance(x, np.ndarray):
values = x.ravel()
else:
raise TypeError("Expecting a NumPy vector or Pandas series.")
x = values
out = {}
out["mean"] = np.nanmean(x)
out["std_ddof0"] = np.nanstd(x, ddof=0)
out["std_ddof1"] = np.nanstd(x, ddof=1)
out["rsd_classical"] = _relative_spread(out["std_ddof1"], out["mean"])
(
out["percentile_05"],
out["percentile_25"],
out["median"],
out["percentile_75"],
out["percentile_95"],
) = np.nanpercentile(x, [5, 25, 50, 75, 95])
out["iqr"] = out["percentile_75"] - out["percentile_25"]
out["min"] = np.nanmin(x)
out["max"] = np.nanmax(x)
out["N_non_missing"] = np.sum(~np.isnan(x))
if method.lower() == "robust":
out["center"], out["spread"] = out["median"], Sn(x)
if ((out["max"] - out["min"]) > 0) and (out["spread"] == 0) and (out["N_non_missing"] > 0):
# Don't fully trust the Sn() yet. It works strangely on quantized data when there is
# little variation. Replace the RSD with the classically calculated version in this
# very specific case. This example shows it: [99, 95, 95, 100, 100, 100, 100, 95, 100,
# 100, 100, 100, 105, 105, 100, 95, 105, 100, 95, 100]
out["center"], out["spread"] = out["mean"], out["std_ddof1"]
else:
out["center"], out["spread"] = out["mean"], out["std_ddof1"]
out["rsd"] = _relative_spread(out["spread"], out["center"])
return out
[docs]
def detect_outliers_esd(
x: np.ndarray | pd.Series, algorithm: str = "esd", max_outliers_detected: int = 1, **kwargs
) -> tuple[list[int], defaultdict[Any, Any]]:
"""
Return a list of indexes of points in the vector `x` which are likely outliers.
A second output (can be ignored) contains the details of the values used to make the decision.
Arguments:
x {list, sequence, NumPy vector/array} -- [A sequence, list or vector which can be
unravelled.]
Keyword Arguments:
algorithm -- Two algorithms are possible to detect outliers: (default: "esd")
'esd': Generalized ESD Test for Outliers.
If `max_outliers_detected=1` this is essentially Grubb's test.
For more details, please see:
https://www.itl.nist.gov/div898/handbook/eda/section3/eda35h3.htm
'cc-robust': Build a robust control-chart for the sequence `x` and
points should lie outside the +/- 3 sigma limits are considered
outliers.
Not Implemented Yet: left here as an idea for the future, but not
confirmed yet. Passing ``algorithm='cc-robust'`` (or any other value
besides ``'esd'``) silently returns an empty outlier list and an
empty details dict; it does not raise.
max_outliers_detected -- The maximum number of outliers that
should be detected, as required by the algorithms.
kwargs -- Algorithm dependent arguments. Defaults are shown here.
'esd':
'robust_variant' = False. When True, uses the median and MAD
for the center and the spread respectively. WARNING: the ESD
critical values (lambda_i) are derived for the classical
statistic max|x - mean| / std, which is bounded above by
(N-1)/sqrt(N); a MAD-scaled statistic has no such bound and is
routinely 2-3x larger on contaminated data, so the robust
variant declares far more outliers than the nominal alpha
implies (it can flag points in perfectly clean data). It was
previously the default; it is now opt-in and should be treated
as a screening heuristic, not a calibrated test.
'alpha' = 0.05. The significance level of the testing.
"""
algorithm = algorithm.strip().lower()
x = pd.Series(x).reset_index(drop=True)
max_outliers_detected = int(max_outliers_detected)
if algorithm == "esd":
"""
https://www.itl.nist.gov/div898/handbook/eda/section3/eda35h3.htm
Note: the two-sided test is implemented here. p = 1-alpha/(2*(n-i+1))
"""
# 0. Default settings. robust_variant defaults to False: the ESD
# critical values are calibrated for the mean/std statistic (see the
# docstring warning about the MAD-scaled variant).
robust_variant = bool(kwargs.get("robust_variant", False))
alpha = float(kwargs.get("alpha", 0.05))
if not 0.0 < alpha <= 1.0:
raise ValueError(f"alpha must lie in (0.0, 1.0], got {alpha}.")
n_observed = int(pd.Series(x).notna().sum())
# The ESD statistic for the i-th suspected outlier uses t on
# dof = N - i - 1 degrees of freedom, so testing r outliers needs
# r <= N - 2. Beyond that the critical value is undefined (t.ppf on
# dof <= 0 is NaN), every later iteration is silently inert, and an
# empty result would read as "no outliers" rather than "not testable".
if max_outliers_detected > 0 and max_outliers_detected > n_observed - 2:
raise ValueError(
f"max_outliers_detected ({max_outliers_detected}) cannot exceed the sample size minus 2 "
f"({n_observed - 2}); the generalized ESD test needs at least two observations beyond "
f"the outliers it tests for. The sample has {n_observed} non-missing observations."
)
# 1. Run K-S test first to check normality
# 2. https://www.itl.nist.gov/div898/handbook/eda/section3/eda35h3.htm
# `sample` shrinks by one point per iteration; keep it separate from
# the `x` parameter so its type stays a Series throughout.
sample: pd.Series = pd.Series(x).copy(deep=True)
extra_out: defaultdict[str, list[Any]] = defaultdict(list)
N = sample.shape[0] - pd.isna(sample).sum()
for k in range(max_outliers_detected):
i = k + 1
extra_out["i"].append(i)
if robust_variant:
variation = median_absolute_deviation(sample.to_numpy())
R = ((sample - sample.median()) / variation).abs()
else:
variation = sample.std()
R = ((sample - sample.mean()) / variation).abs()
# The ESD critical value is defined against the ORIGINAL sample
# size N (NIST / Rosner), so N is deliberately not refreshed here.
dof = N - i - 1
p = 1 - alpha / (2 * (N - i + 1))
t_s = t_value(p, dof)
lambda_i = (N - i) * t_s / np.sqrt((dof + np.power(t_s, 2)) * (dof + 2))
extra_out["lambda"].append(lambda_i)
g = R.max()
# The Grubbs p-value, by contrast, describes the sample actually in
# hand this iteration, which has shrunk by every point dropped so
# far. Using the original N here (the previous behaviour) reported
# a p-value for a larger sample than the statistic came from.
n_i = int(sample.notna().sum())
s = g**2 * n_i * (2 - n_i) / (g**2 * n_i - (n_i - 1) ** 2) # R-function Grubbs formula
p_value = 0 if s <= 0 else min(n_i * (1 - t_value_cdf(np.sqrt(s), n_i - 2)), 1)
# R is all-NaN when the spread is zero, so there is no candidate to
# drop. Record None rather than the -1 sentinel used before, which
# is not a label in the reset RangeIndex and raised KeyError from
# the drop below whenever the spread was non-zero (an infinity in
# the data reaches exactly that state).
R_i_idx: Any | None = None if pd.isna(g) else R.idxmax()
extra_out["R_i_idx"].append(R_i_idx)
extra_out["R_i"].append(g)
extra_out["p-value"].append(p_value)
if R_i_idx is not None and variation > __eps:
# The variation, if zero or small, will fail to drop this index
sample = sample.drop(R_i_idx)
try:
# NIST / Rosner: "the number of outliers is determined by finding
# the LARGEST i such that R_i > lambda_i". The crossings are not
# monotone in general (that is the masking effect this test
# exists to handle), so take the last crossing, not the first.
cutoff_i = np.where(
np.array(extra_out["R_i"]) - np.array(extra_out["lambda"]) >= 0,
)[0][-1]
except IndexError:
cutoff_i = -1
# The outlier indices are the points from the start of
# `extra_out['R_i_idx']`, up to, and including, the `cutoff_i` index
# in the list.
extra_out["cutoff"] = cutoff_i
outlier_index = [idx for idx in extra_out["R_i_idx"][0 : (cutoff_i + 1)] if idx is not None]
return outlier_index, extra_out
else:
return [], defaultdict(dict)
[docs]
def tietjen_moore_test( # noqa: PLR0913
x: np.ndarray | pd.Series,
n_outliers: int,
*,
two_sided: bool = True,
alpha: float = 0.05,
n_simulations: int = 10000,
random_state: int | None = None,
) -> Bunch:
"""Tietjen-Moore test for a *specified* number of outliers.
Tests the null hypothesis "there are no outliers" against the alternative
"the ``n_outliers`` most extreme observations are outliers". Unlike the
generalised ESD test, the number of suspected outliers must be fixed in
advance. The test statistic has no closed-form critical value, so it is
obtained by simulation under the normal null.
Parameters
----------
x : np.ndarray or pd.Series
One-dimensional sample. Missing values are dropped.
n_outliers : int
The number of suspected outliers to test for (1 <= n_outliers < N).
two_sided : bool, optional
If ``True`` (default) the test looks for outliers on either tail (the
observations with the largest absolute deviation from the mean). If
``False`` it tests only the ``n_outliers`` largest observations.
alpha : float, optional
Significance level, by default 0.05.
n_simulations : int, optional
Number of Monte-Carlo samples used to estimate the critical value.
random_state : int or None, optional
Seed for the simulation, for reproducibility.
Returns
-------
sklearn.utils.Bunch
With ``statistic``, ``critical_value``, ``reject`` (``True`` when the
outliers are significant), ``outlier_indices`` (positions in the
missing-value-removed sample), ``n_outliers`` and ``alpha``.
References
----------
Tietjen and Moore, "Some Grubbs-type statistics for the detection of
several outliers", Technometrics, 14, 583-597, 1972. See also the NIST
handbook, section 3.5.h.3.
"""
a = np.asarray(x, dtype=float).ravel()
a = a[~np.isnan(a)]
n = a.size
if not (1 <= n_outliers < n):
raise ValueError(f"n_outliers must be between 1 and {n - 1}, got {n_outliers}.")
def _statistic(sample: np.ndarray) -> float:
ranking = np.argsort(np.abs(sample - sample.mean())) if two_sided else np.argsort(sample)
kept = sample[ranking[: n - n_outliers]]
denominator = np.sum((sample - sample.mean()) ** 2)
return float(np.sum((kept - kept.mean()) ** 2) / denominator)
observed = _statistic(a)
# Simulate the null distribution: all observations i.i.d. standard normal.
rng = np.random.default_rng(random_state)
simulated = np.array([_statistic(rng.standard_normal(n)) for _ in range(n_simulations)])
critical_value = float(np.quantile(simulated, alpha))
outlier_indices = (
np.argsort(np.abs(a - a.mean()))[n - n_outliers :] if two_sided else np.argsort(a)[n - n_outliers :]
)
# A small statistic indicates that the removed points really are outliers.
return Bunch(
statistic=observed,
critical_value=critical_value,
reject=observed < critical_value,
outlier_indices=np.sort(outlier_indices),
n_outliers=n_outliers,
alpha=alpha,
)
[docs]
def distribution_fit(
x: np.ndarray | pd.Series,
distribution: str = "norm",
alpha: float = 0.05,
) -> Bunch:
"""Check how well a sample fits a named distribution.
Fits the parameters of the requested ``scipy.stats`` distribution by
maximum likelihood and runs a Kolmogorov-Smirnov goodness-of-fit test
(NIST handbook, section 3.5.7).
Parameters
----------
x : np.ndarray or pd.Series
One-dimensional sample. Missing values are dropped.
distribution : str, optional
Name of any continuous ``scipy.stats`` distribution, by default
``"norm"``.
alpha : float, optional
Significance level for the ``fits_well`` verdict, by default 0.05.
Returns
-------
sklearn.utils.Bunch
With ``distribution``, fitted ``parameters``, ``ks_statistic``,
``ks_pvalue``, ``fits_well`` (``True`` when the fit is not rejected at
level ``alpha``) and the sample size ``n``.
Notes
-----
Because the distribution parameters are estimated from the same data, the
KS p-value is conservative (the true Type-I error is smaller than
``alpha``); it remains a useful screening check.
"""
from scipy import stats as scipy_stats # noqa: PLC0415
a = np.asarray(x, dtype=float).ravel()
a = a[~np.isnan(a)]
dist = getattr(scipy_stats, distribution)
parameters = dist.fit(a)
# Compare against the frozen, fully-parameterised distribution. Passing the
# frozen ``.cdf`` (which applies loc/scale internally) rather than the
# distribution name plus ``args`` keeps this working across scipy versions,
# where the latter form can forward loc/scale positionally to the raw cdf.
ks_statistic, ks_pvalue = scipy_stats.kstest(a, dist(*parameters).cdf)
return Bunch(
distribution=distribution,
parameters=parameters,
ks_statistic=float(ks_statistic),
ks_pvalue=float(ks_pvalue),
fits_well=bool(ks_pvalue > alpha),
n=int(a.size),
)
[docs]
def variance_decomposition(df: pd.DataFrame, measured: str, repeat: str) -> dict:
"""
Given a DataFrame `df` of raw data, and an indication of which column is the `measured` value
column, and which is the `repeat` indicator, it will calculate the within and between replicate
standard deviation.
Example
Two measurements on day 1 ``[101, 102]`` and two measurements on day 2 ``[94, 95]``. The
between-day variation can already be expected to be much greater than the within-day variation.
>>> df = pd.DataFrame(data={'Result': [101, 102, 94, 95], 'Repeat': [1, 1, 2, 2]})
Result Repeat
0 101 1
1 102 1
2 94 2
3 95 2
>>> output = variance_decomposition(df, measured="Result", repeat="Repeat")
{'total_ms': 16.666667,
'total_dof': 3,
'within_ms': 0.5,
'within_stddev': 0.70711,
'within_dof': 2,
'between_ms': 49.0,
'between_stddev': 4.9244,
'between_dof': 1}
Notes
-----
* SSQ = sum of squares
* DOF = degrees of freedom
* MS = mean square = (sum of squares) / (degrees of freedom) = SSQ / DOF = variance
* ``between_ms`` is the raw ANOVA mean square, which estimates the sum of
the within-group variance and n0 times the between-group variance (n0 =
average group size), NOT the between-group variance itself.
``between_stddev`` reports the actual between-group variance COMPONENT,
the square root of ``max(0, (MS_between - MS_within) / n0)``. Earlier
versions reported ``sqrt(MS_between)``, which mixes the within-group
variance into the "between" number.
"""
# Overall statistics:
total_ms = max(0, df[measured].var())
total_dof = max(0, df[measured].count() - 1)
# Within a group: calculate the standard deviation (of variance) of each group. If you take
# the average of those variances, pooled, you get an estimate of the within-group variance.
within_ms = 0.0
within_dof = 0
for _, group in df.groupby(repeat):
dof_group_i = max(0, group[measured].count() - 1)
within_dof += dof_group_i
# handling missing data makes for messier code
within_ms = np.nansum([within_ms, group[measured].var() * dof_group_i])
within_ms = 0 if within_dof == 0 else within_ms / within_dof
# Between groups: the ANOVA relationship is such that SSQ(total) = SSQ(between) + SSQ(within).
# Therefore the between-group statistics are found by differencing
between_dof = max(0, total_dof - within_dof)
between_ms = 0 if between_dof == 0 else max(0.0, (total_ms * total_dof - within_ms * within_dof) / between_dof)
# MS_between estimates sigma_within^2 + n0 * sigma_between^2, so the
# between-group VARIANCE COMPONENT is (MS_between - MS_within) / n0,
# clipped at zero when the between mean square is smaller than the noise
# floor. n0 is the effective (average) group size; for unbalanced data the
# standard ANOVA n0 = (N - sum(n_i^2)/N) / (k - 1) is used.
group_sizes = np.array([group[measured].count() for _, group in df.groupby(repeat)], dtype=float)
n_total = float(group_sizes.sum())
k_groups = int((group_sizes > 0).sum())
n0 = (n_total - float(np.sum(group_sizes**2)) / n_total) / (k_groups - 1) if k_groups > 1 and n_total > 0 else 0.0
between_variance_component = max(0.0, (between_ms - within_ms) / n0) if n0 > 0 else 0.0
return {
"total_ms": total_ms,
"total_dof": total_dof,
"within_ms": within_ms,
"within_stddev": np.sqrt(within_ms),
"within_dof": within_dof,
"between_ms": between_ms,
"between_stddev": np.sqrt(between_variance_component),
"between_dof": between_dof,
}
# ---------------------------------------------------------------------------
# Migration helpers - old names raise helpful errors
# ---------------------------------------------------------------------------
_RENAMED = {
"median_abs_deviation": "median_absolute_deviation",
"normality_check": "test_normality",
"within_between_standard_deviation": "variance_decomposition",
"outlier_detection_multiple": "detect_outliers_esd",
"ttest_difference_calculate": "ttest_independent",
"ttest_paired_difference_calculate": "ttest_paired",
"ttest_difference": "ttest_independent_from_df",
"ttest_paired_difference": "ttest_paired_from_df",
}
def __getattr__(name: str) -> NoReturn:
"""Raise a helpful error when a renamed module attribute is accessed."""
if name in _RENAMED:
new = _RENAMED[name]
raise AttributeError(
f"{name!r} has been renamed to {new!r}. Use: from process_improve.univariate.metrics import {new}"
)
raise AttributeError(f"module {__name__!r} has no attribute {name!r}")