from __future__ import annotations
import warnings
from typing import Any, Literal
import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.exceptions import NotFittedError
from sklearn.utils import Bunch
from .._random import check_random_state
from .._warnings import user_stacklevel
from ..univariate.metrics import _check_probability, t_value
__eps = np.finfo(np.float32).eps
[docs]
def fit_robust_lm(x: np.ndarray, y: np.ndarray) -> np.ndarray:
"""
Fits a robust linear model between Numpy vectors `x` and `y`, with an
intercept. Returns a length-2 array ``[intercept, slope]`` (the
``params`` attribute returned by ``statsmodels.RLM``); no extra checking
on data consistency is done.
See also: regression.repeated_median_slope
"""
rlm_model = sm.RLM(y, np.vstack([np.ones(x.size), x.ravel()]).T, M=sm.robust.norms.HuberT())
rlm_results = rlm_model.fit()
return rlm_results.params
def _drop_non_finite_pairs(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""
Drop ``(x, y)`` pairs holding a NaN or an infinity in either vector.
A NaN in either vector makes every slope through that point undefined, so the inner
``np.nanmedian`` of :func:`repeated_median_slope` used to receive an all-NaN list: it
emitted "All-NaN slice encountered" and returned NaN, which the outer median then
quietly dropped. The slope was therefore computed from whichever points happened to be
clean, with nothing in the return value recording the omission. (#558)
Parameters
----------
x, y : np.ndarray
Equal-length 1-D arrays; the caller refuses vectors of different lengths.
Returns
-------
tuple[np.ndarray, np.ndarray]
The finite subset of each input.
Raises
------
ValueError
If fewer than three finite pairs remain.
"""
finite = np.isfinite(x) & np.isfinite(y)
if finite.all():
return x, y
n_finite = int(finite.sum())
if n_finite < 3:
raise ValueError(
f"repeated_median_slope requires at least 3 finite (x, y) pairs; got {n_finite} "
f"of {x.size}. Non-finite pairs (NaN or infinite in either vector) are omitted."
)
return x[finite], y[finite]
def _repeated_median(x: np.ndarray, y: np.ndarray) -> float:
"""
Siegel's repeated-median slope of finite, equal-length 1-D arrays, without input checks.
Row ``i`` of an N x N matrix holds the slopes ``(y[j] - y[i]) / (x[j] - x[i])``, with
NaN on the diagonal and wherever ``x[j] == x[i]``. Sorting each row puts its NaN last,
so the row's median sits at fixed positions among its defined slopes. A row with no
defined slope has no median and is left out of the outer one; with none at all the
slope is NaN. At the SEC-19 cap of 5,000 points the matrix and the x differences it is
divided by take 400 MB together.
Parameters
----------
x, y : np.ndarray
Finite 1-D arrays of the same length.
Returns
-------
float
The median over the rows of each row's median slope.
"""
x = np.asarray(x, dtype=float)
y = np.asarray(y, dtype=float)
dx = x[None, :] - x[:, None]
slopes = y[None, :] - y[:, None]
with np.errstate(divide="ignore", invalid="ignore"):
np.divide(slopes, dx, out=slopes)
slopes[dx == 0] = np.nan
del dx
slopes.sort(axis=1)
n_defined = np.count_nonzero(~np.isnan(slopes), axis=1)
rows = np.flatnonzero(n_defined)
if rows.size == 0:
return float("nan")
count = n_defined[rows]
# (a + b) / 2 is how np.median averages the middle pair, and gives a when a == b.
row_medians = (slopes[rows, (count - 1) // 2] + slopes[rows, count // 2]) / 2
return float(np.median(row_medians))
def _bootstrap_lines(
x: np.ndarray, y: np.ndarray, fit_intercept: bool, n_boot: int, rng: np.random.Generator
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Refit the robust line to ``n_boot`` resamples of the ``(x, y)`` pairs.
Each resample draws N pairs with replacement and is fitted as :func:`robust_regression`
fits the data: the repeated-median slope, then the median intercept. A resample whose
x values are all equal has no slope and is left out.
Parameters
----------
x, y : np.ndarray
Finite 1-D arrays of the same length.
fit_intercept : bool
If ``False``, every resampled line passes through the origin.
n_boot : int
The number of resamples.
rng : np.random.Generator
The source of every random draw.
Returns
-------
slopes, intercepts, errors : np.ndarray
One value per resample that has a slope. ``errors`` is a normal draw scaled by
1.4826 times the median absolute residual of that resample: the error of one new
observation, on a scale that outliers move as little as they move the line.
"""
if not isinstance(n_boot, (int, np.integer)) or n_boot < 1:
raise ValueError(f"n_boot must be a positive integer; got {n_boot!r}.")
picks = rng.integers(0, x.size, size=(n_boot, x.size))
slopes = np.array([_repeated_median(x[pick], y[pick]) for pick in picks])
has_slope = ~np.isnan(slopes)
if not has_slope.any():
# Only a few points and a few resamples can get here. One NaN each makes every interval NaN.
return np.full(1, np.nan), np.full(1, np.nan), np.full(1, np.nan)
x_b, y_b, slopes = x[picks[has_slope]], y[picks[has_slope]], slopes[has_slope]
lines = slopes[:, None] * x_b
intercepts = np.median(y_b - lines, axis=1) if fit_intercept else np.zeros(slopes.size)
scale = 1.4826 * np.median(np.abs(y_b - lines - intercepts[:, None]), axis=1)
return slopes, intercepts, scale * rng.standard_normal(slopes.size)
[docs]
def robust_regression( # noqa: PLR0913, PLR0915
x: np.ndarray | pd.DataFrame | pd.Series,
y: np.ndarray | pd.DataFrame | pd.Series,
fit_intercept: bool = True,
na_rm: bool = True,
conflevel: float = 0.95,
nowarn: bool = False,
pi_resolution: int = 50,
*,
intervals: Literal["bootstrap", "least_squares"] | None = None,
n_boot: int = 999,
random_state: int | np.random.Generator | None = 0,
) -> dict:
"""
Perform the Simple robust regression analysis between `x` and `y` variables.
Simple robust regression between an `x` and a `y` using the `repeated_median_slope` method
to calculate the slope. The intercept is the median intercept, when using that slope and the
provided `x` and `y` values, or forced to zero if fit_intercept=False.
Parameters
----------
x, y : np.ndarray, pd.DataFrame or pd.Series
Sequences of numerical values, paired by position.
fit_intercept : bool, default True
If True, fits an intercept term. If False, forces regression through origin.
na_rm : bool, default True
If True, removes all observations with one or more missing values.
conflevel : float, default 0.95
Confidence level for the confidence and prediction intervals, strictly between 0 and 1.
nowarn : bool, default False
If True, skips the minimum-size check and the SEC-19 cap on the number of points.
Users should ensure data validity beforehand.
pi_resolution : int, default 50
The number of x values, evenly spaced over the range of the data, at which
``pi_range`` gives the prediction interval.
intervals : {"bootstrap", "least_squares"} or None, default None
How ``conf_intervals``, ``conf_interval_intercept`` and ``pi_range`` are found.
``"bootstrap"`` takes percentiles over ``n_boot`` resamples of the ``(x, y)`` pairs,
each refitted with the repeated median, so the intervals resist outliers as the
slope does. ``"least_squares"`` applies the least-squares formulas to the residuals
around the robust line: one gross outlier widens those intervals even when it leaves
the slope where it was. ``None`` means ``"least_squares"`` until 3.0 and
``"bootstrap"`` from 3.0. Until then, a call whose intervals the change would move
issues a ``FutureWarning`` and returns its text as ``intervals_note`` (see Notes).
Passing either value silences it.
n_boot : int, default 999
The number of bootstrap resamples. Used only with ``intervals="bootstrap"``.
random_state : int, np.random.Generator or None, default 0
Seeds the bootstrap resamples; see
:func:`process_improve._random.check_random_state`. The default seed makes a
repeated call return the same intervals. Used only with ``intervals="bootstrap"``.
Returns
-------
dict
The fit, with these keys::
N: the number of observations used to fit the model
coefficients: a length-1 list containing the regression slope
intercept: the fitted intercept on the full-fit path
(0 when ``fit_intercept=False``); ``np.nan``
on the degenerate early-return path (see
below), regardless of ``fit_intercept``
standard_errors: a length-1 list containing the standard error of the slope
standard_error_intercept: standard error for the intercept (np.nan if fit_intercept=False)
R2: RegSS / (RegSS + RSS), around the robust line. Around a line
that is not the least-squares one, RegSS + RSS is not the
total sum of squares, so this is not 1 - RSS / TSS (the R2
of ``multiple_linear_regression``)
SE: the model's standard error
x_ssq: the sum of squares of (x - mean(x))
k: the number of model parameters (2 if fit_intercept else 1)
fitted_values: the N predicted values, one per row in y
residuals: the N residuals
t_value: the t-values for the standard errors, one per
coefficient (same meaning as in
``multiple_linear_regression``)
conf_intervals: 1 row x 2 columns (lower, upper) confidence interval for
the single slope coefficient, found as ``intervals`` says
conf_interval_intercept: (lower, upper) confidence interval for the intercept, found
as ``intervals`` says (np.nan if fit_intercept=False)
pi_range: pi_resolution rows x 3 columns: x, then the lower and upper
prediction limits for one new observation at that x, found
as ``intervals`` says
leverage: the hat-matrix diagonal (leverage) for each observation
influence: Cook-style influence values for each observation
intervals_note: only when ``intervals`` is not given and the default's change
in 3.0 would change the intervals: the text of the
``FutureWarning`` that says so
On the degenerate early-return path (fewer than three usable observations
in ``x`` or ``y``), the returned dictionary is the initialization stub:
``N`` is ``None``, ``x_ssq`` is absent, and ``leverage``, ``influence``,
``pi_range``, ``t_value``, ``fitted_values``, ``residuals``, and the
coefficient/interval arrays hold ``np.nan``.
Notes
-----
By default the three interval keys are least-squares intervals, as they were up to 2.2.2.
2.3.0 made the bootstrap the default without the announcement that the deprecation policy
requires for a changed default, so code that kept running got other numbers; 2.3.1
restored least squares. The default becomes ``"bootstrap"`` in 3.0, because its intervals
resist outliers as the slope does. Until then a call that the change would affect (one
with spread in x, a finite ``SE`` and residuals that are more than rounding) warns and
carries ``intervals_note``. A perfect fit stays quiet: every resample lies on its line,
so both kinds of interval have zero width.
A ``conflevel`` outside (0, 1), NaN included, raises ``ValueError`` before any work is
done. Up to 2.4.0 a level between -1 and 0 gave a reversed interval, and one with no t
value gave NaN intervals.
With ``intervals="bootstrap"`` the three interval keys are percentile
intervals over the resamples, each refitted as above. A resample whose x values are all
equal has no slope and is left out. The prediction limits add to each resampled line a
normal error scaled by 1.4826 times that resample's median absolute residual, so an
outlier widens them no more than it moves the line. Without outliers the slope interval
is a little wider than the least-squares one, because it measures the spread of the
repeated-median slope, which uses the data less efficiently than least squares does.
``SE``, the standard errors, ``t_value``, ``R2``, ``leverage`` and ``influence`` use the
least-squares formulas, applied to the residuals around the robust line, whichever
``intervals`` is: one gross outlier, which leaves the slope almost where it was, inflates
the residual sum of squares and so raises ``SE`` and lowers ``R2``.
The bootstrap refits the line ``n_boot`` times, and each fit grows as N squared (see
:func:`repeated_median_slope`). The default 999 resamples take about 0.1 s for 100
points, 3 s for 500 and over 10 s for 1,000; for more points, lower ``n_boot`` or use
``intervals="least_squares"``.
"""
if intervals not in (None, "bootstrap", "least_squares"):
raise ValueError(f"intervals must be 'bootstrap' or 'least_squares'; got {intervals!r}.")
_check_probability("conflevel", conflevel)
out: dict[str, Any] = {
"N": None,
"coefficients": [
np.nan,
],
"intercept": np.nan,
"standard_errors": [
np.nan,
],
"standard_error_intercept": np.nan,
"R2": np.nan,
"SE": np.nan,
"fitted_values": np.nan,
"residuals": np.nan,
"t_value": np.nan,
"conf_intervals": np.array([[np.nan, np.nan]]),
"conf_interval_intercept": np.array([np.nan, np.nan]),
"pi_range": np.nan,
"leverage": np.nan,
"k": 1,
"influence": np.nan,
}
# Data pre-processing: handle both Pandas and NumPy -> use Pandas internally for X and y
x_ = pd.DataFrame(x, copy=True) if isinstance(x, np.ndarray) else pd.DataFrame(x.values, copy=True)
y_ = pd.DataFrame(y.ravel(), copy=True) if isinstance(y, np.ndarray) else pd.DataFrame(y.values, copy=True)
# Removing missing values:
if na_rm:
missing_idx = y_.isna().any(axis=1) | x_.isna().any(axis=1)
x_ = x_.loc[~missing_idx, :]
y_ = y_.loc[~missing_idx]
# CASE when there is no data, or only 2 data point (repeate median slope needs more than 2)
if (y_.size <= 2) or (x_.size <= 2):
return out
x = x_.values.ravel()
y = y_.values.ravel()
# initialize statistical variables
k_params = 2 if fit_intercept else 1
dof_resid = len(x) - k_params
# Calculate robust regression
slope = repeated_median_slope(x, y, nowarn=nowarn)
intercept = np.nanmedian(y - slope * x) if fit_intercept else 0.0
mean_x, mean_y = np.mean(x), np.mean(y)
out["N"] = len(x)
out["intercept"] = intercept
out["coefficients"] = [
slope,
]
out["fitted_values"] = intercept + slope * x
out["residuals"] = y - out["fitted_values"]
# For robust method, calculate this way, since no guarantee both RegSS or RSS are < TSS.
# So ensure this way that the TSS = RegSS + RSS, and R2 is the ratio of RegSS/TSS
# https://learnche.org/pid/least-squares-modelling/least-squares-model-analysis
regression_ssq = np.sum(np.power(out["fitted_values"] - mean_y, 2))
residual_ssq = np.sum(out["residuals"] * out["residuals"])
total_ssq = regression_ssq + residual_ssq
out["R2"] = regression_ssq / total_ssq
out["SE"] = np.sqrt(residual_ssq / dof_resid)
out["x_ssq"] = np.sum(np.power(x - mean_x, 2))
out["k"] = k_params
# t-critical value for confidence intervals
c_t = t_value(1 - (1 - conflevel) / 2, dof_resid)
# Prediction intervals
pi_range = np.linspace(np.min(x), np.max(x), pi_resolution)
pi_y_pred = out["intercept"] + out["coefficients"][0] * pi_range
# Handle degenerate case where x has no variation
if out["x_ssq"] < __eps:
out["standard_error_intercept"] = np.nan
out["conf_interval_intercept"] = np.array([np.nan, np.nan])
out["standard_errors"] = [
np.nan,
]
out["pi_range"] = np.vstack([pi_range, pi_y_pred, pi_y_pred]).T
# With no variation in x every observation carries the same weight:
# the leverage is exactly 1/N. The previous code divided the ~zero
# deviations by the ~zero x_ssq (0/0), poisoning the leverage and
# the influence values computed from it with NaN/inf - inside the
# very branch that exists to guard against the degenerate x.
out["leverage"] = np.full(int(out["N"]), 1.0 / out["N"])
else:
out["standard_errors"] = [
out["SE"] * 1 / np.sqrt(out["x_ssq"]),
]
if fit_intercept:
out["standard_error_intercept"] = SE_b0 = out["SE"] * np.sqrt(1 / out["N"] + (mean_x) ** 2 / out["x_ssq"])
out["conf_interval_intercept"] = np.array(
[out["intercept"] - c_t * SE_b0, out["intercept"] + c_t * SE_b0],
)
var_y = (out["SE"] ** 2) * (1 + 1 / out["N"] + (pi_range - np.mean(x)) ** 2 / out["x_ssq"])
out["leverage"] = 1 / out["N"] + np.power(x - mean_x, 2) / out["x_ssq"]
else: # 1 / out["N"] is the term for the uncertainty for the intercept
out["standard_error_intercept"] = np.nan
out["conf_interval_intercept"] = np.array([np.nan, np.nan])
# The model without intercept centers around the origin
x_ssq_origin = np.sum(x**2)
var_y = (out["SE"] ** 2) * (1 + (pi_range**2) / x_ssq_origin)
out["leverage"] = np.power(x, 2) / x_ssq_origin
std_y = np.sqrt(var_y)
lower = pi_y_pred - c_t * std_y
upper = pi_y_pred + c_t * std_y
out["pi_range"] = np.vstack([pi_range, lower, upper]).T
# The t-statistic per coefficient: coefficient / standard error. This
# matches `multiple_linear_regression`, so the two functions agree on what
# the "t_value" key means. It is NaN on the degenerate-x path, where the
# standard error is itself undefined. Note `c_t` is a different quantity:
# the critical value used to set the interval widths.
out["t_value"] = np.array([out["coefficients"][0] / out["standard_errors"][0]])
out["conf_intervals"] = np.array(
[
[
out["coefficients"][0] - c_t * out["standard_errors"][0],
out["coefficients"][0] + c_t * out["standard_errors"][0],
],
]
)
# With no spread in x, or a non-finite value kept in (SE is then not finite), the NaN intervals above stand
# whichever intervals is, so the default's change in 3.0 changes nothing.
has_intervals = out["x_ssq"] >= __eps and np.isfinite(out["SE"])
# Nor does it for a perfect fit, whose residuals are rounding errors of a few units in the last place of the
# line's largest term (at most 6 over 4,000 random exact lines): every resample lies on that line, so both
# kinds of interval have zero width.
largest_term = max(np.abs(y).max(), abs(intercept), np.abs(slope * x).max())
perfect_fit = out["SE"] <= 16 * np.finfo(float).eps * largest_term
if intervals is None and has_intervals and not perfect_fit:
out["intervals_note"] = (
"process_improve.regression.robust_regression's default intervals='least_squares' is deprecated since "
"2.3.1 and will change to intervals='bootstrap' in 3.0, whose percentile-bootstrap intervals resist "
"outliers as the slope does. Pass intervals='least_squares' to keep these numbers, or "
"intervals='bootstrap' for the new ones; either silences this warning."
)
warnings.warn(out["intervals_note"], FutureWarning, stacklevel=user_stacklevel())
if intervals == "bootstrap" and has_intervals:
slopes, intercepts, errors = _bootstrap_lines(x, y, fit_intercept, n_boot, check_random_state(random_state))
tails = [50 * (1 - conflevel), 50 * (1 + conflevel)]
out["conf_intervals"] = np.percentile(slopes, tails)[None, :]
if fit_intercept:
out["conf_interval_intercept"] = np.percentile(intercepts, tails)
new_observations = intercepts[:, None] + slopes[:, None] * pi_range + errors[:, None]
out["pi_range"] = np.column_stack([pi_range, *np.percentile(new_observations, tails, axis=0)])
if out["SE"] < __eps:
out["influence"] = out["residuals"] * 0.0
else:
out["influence"] = (
np.power(out["residuals"] / ((1 - out["leverage"]) * out["SE"]), 2) * out["leverage"] / out["k"]
)
return out
_RENAMED = {"simple_robust_regression": "robust_regression"}
def __getattr__(name: str) -> None:
"""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.regression.methods import {new}"
)
raise AttributeError(f"module {__name__!r} has no attribute {name!r}")
[docs]
def multiple_linear_regression(
X: np.ndarray | pd.DataFrame,
y: np.ndarray | pd.DataFrame | pd.Series,
fit_intercept: bool = True,
na_rm: bool = True,
conflevel: float = 0.95,
pi_resolution: int = 50,
) -> dict:
"""
Linear regression of the N rows and K columns of matrix `X` onto the single column 'y'.
Backwards-compatible wrapper around :class:`OLS`. New code should use the
:class:`OLS` estimator directly, which exposes the same statistics as
sklearn-style attributes and prints an R-like ``summary(lm(...))``.
Notes and limitations:
* does not handle weighting
* N >= K at least as many rows as columns in X
Returns a dictionary of outputs. Keys always present::
N: number of observations actually used to fit
coefficients: a vector of K coefficients, one for each column in X;
shape ``[np.nan]`` on the degenerate/unfitted path
intercept: the fitted intercept; ``np.nan`` when
``fit_intercept=False`` or on the
degenerate/unfitted path
standard_errors: a vector of K standard errors, one per column in X
standard_error_intercept: standard error for the intercept
R2: the R^2 value
R2_regression_based: R^2 computed as ``RegSS / TSS`` (added post-fit)
R2_residual_based: R^2 computed as ``1 - RSS / TSS`` (added post-fit)
k: the number of model parameters (added post-fit)
SE: the model's standard error
fitted_values: the N predicted values, one per row in y
residuals: the N residuals
t_value: the t-values for the standard errors
conf_intervals: K rows x 2 columns (lower, upper) confidence intervals;
shape ``[np.nan, np.nan]`` on the degenerate/unfitted path
conf_interval_intercept: (lower, upper) confidence interval for the intercept
(added post-fit)
Keys present whenever the fit succeeds (populated for any K)::
leverage: hat-matrix diagonal
influence: Cook-style influence values
Keys present only for single-feature ``X`` (and only when ``fit_intercept``
is True and there is enough non-degenerate data)::
x_ssq: sum of squares of the centred predictor
pi_range: prediction interval above and below, over the
range of the predictor (``pi_resolution`` points)
"""
model = OLS(
fit_intercept=fit_intercept,
na_rm=na_rm,
conflevel=conflevel,
pi_resolution=pi_resolution,
)
return model.fit(X, y).to_dict()
[docs]
class OLS(RegressorMixin, BaseEstimator):
r"""Ordinary Least Squares regression with statistical diagnostics.
A scikit-learn-compatible estimator that fits an OLS model and exposes
inferential statistics (standard errors, t-values, p-values, confidence
intervals, F-statistic) and influence diagnostics (leverage, Cook's
distance). Calling ``print(model)`` after fitting renders a summary
similar to R's ``summary(lm(...))``.
Parameters
----------
fit_intercept : bool, default=True
If True, fits an intercept term. If False, the regression is forced
through the origin.
na_rm : bool, default=True
If True, drops rows with one or more missing values before fitting.
conflevel : float, default=0.95
Confidence level for confidence and prediction intervals.
pi_resolution : int, default=50
Number of grid points at which to compute prediction intervals over
the range of x. Only used when X has a single column and an intercept
is fitted.
Attributes
----------
coefficients_ : np.ndarray of shape (K,)
Fitted slope coefficients (excludes the intercept).
intercept_ : float
Fitted intercept (np.nan if ``fit_intercept`` is False).
standard_errors_ : np.ndarray of shape (K,)
Standard errors of ``coefficients_``.
standard_error_intercept_ : float
Standard error of the intercept.
t_values_ : np.ndarray of shape (K,)
t-statistics for each coefficient.
t_value_intercept_ : float
t-statistic for the intercept.
p_values_ : np.ndarray of shape (K,)
Two-sided p-values for each coefficient.
p_value_intercept_ : float
p-value for the intercept.
conf_intervals_ : np.ndarray of shape (K, 2)
Lower and upper bounds of the coefficient confidence intervals.
conf_interval_intercept_ : np.ndarray of shape (2,)
Lower and upper bounds of the intercept confidence interval.
r2_ : float
Coefficient of determination.
adj_r2_ : float
Adjusted R-squared.
se_ : float
Residual standard error (sqrt of residual variance).
df_resid_ : int
Residual degrees of freedom.
df_model_ : int
Model degrees of freedom (number of slope coefficients).
f_statistic_ : float
F-statistic for the overall regression.
f_pvalue_ : float
p-value associated with the F-statistic.
fitted_values_ : np.ndarray of shape (N,)
In-sample predictions.
residuals_ : np.ndarray of shape (N_original,)
In-sample residuals (NaN at rows removed by ``na_rm``).
leverage_ : np.ndarray of shape (N,)
Hat-matrix diagonal :math:`h_i = \mathbf{x}_i^T (X^T X)^{-1} \mathbf{x}_i`, over
the rows that were fitted. Matches ``statsmodels`` ``hat_matrix_diag``.
influence_ : np.ndarray of shape (N,)
Cook's distance per fitted row. Matches ``statsmodels`` ``cooks_distance``.
pi_range_ : np.ndarray of shape (pi_resolution, 3) or float
Columns are x-grid, lower bound, upper bound of the prediction
interval. ``np.nan`` if not applicable.
feature_names_in_ : list[str]
Column names of the feature matrix.
target_name_ : str
Name of the target variable.
n_samples_ : int
Number of samples used in the fit (after ``na_rm``).
n_features_in_ : int
Number of input features.
is_fitted_ : bool
Whether ``fit()`` has been called successfully.
Examples
--------
>>> import numpy as np
>>> from process_improve.regression.methods import OLS
>>> rng = np.random.default_rng(0)
>>> X = rng.standard_normal((50, 2))
>>> y = X @ [1.5, -2.0] + 0.5 + 0.1 * rng.standard_normal(50)
>>> model = OLS().fit(X, y)
>>> print(model) # doctest: +SKIP
Call:
OLS(fit_intercept=True, na_rm=True, conflevel=0.95)
...
See Also
--------
multiple_linear_regression : Backwards-compatible function returning a dict.
robust_regression : Robust regression via repeated-median slope.
"""
def __init__(
self,
fit_intercept: bool = True,
na_rm: bool = True,
conflevel: float = 0.95,
pi_resolution: int = 50,
) -> None:
self.fit_intercept = fit_intercept
self.na_rm = na_rm
self.conflevel = conflevel
self.pi_resolution = pi_resolution
[docs]
def fit( # noqa: C901, PLR0912, PLR0915
self,
X: np.ndarray | pd.DataFrame | pd.Series,
y: np.ndarray | pd.DataFrame | pd.Series,
) -> OLS:
"""Fit the OLS model.
Parameters
----------
X : array-like of shape (N, K)
Feature matrix. Pandas and NumPy inputs are both accepted.
y : array-like of shape (N,) or (N, 1)
Target vector.
Returns
-------
self : OLS
Fitted estimator.
"""
_check_probability("conflevel", self.conflevel)
# Preserve original input metadata (length, names) before we drop NA rows.
if isinstance(X, pd.DataFrame):
X_full = X.copy()
# Default integer column names (the typical case for numpy → DataFrame)
# are remapped to x1, x2, ... so the printed formula and table match R's lm().
if list(X_full.columns) == list(range(X_full.shape[1])):
X_full.columns = [f"x{i + 1}" for i in range(X_full.shape[1])]
elif isinstance(X, pd.Series):
X_full = X.to_frame(name=str(X.name) if X.name is not None else "x1")
else:
X_arr = np.asarray(X)
if X_arr.ndim == 1:
X_arr = X_arr.reshape(-1, 1)
X_full = pd.DataFrame(X_arr, columns=[f"x{i + 1}" for i in range(X_arr.shape[1])])
if isinstance(y, pd.Series):
y_full = pd.DataFrame(y.values.ravel())
target_name = y.name if y.name is not None else "y"
elif isinstance(y, pd.DataFrame):
y_full = pd.DataFrame(y.values.ravel())
target_name = y.columns[0] if len(y.columns) else "y"
else:
y_full = pd.DataFrame(np.asarray(y).ravel())
target_name = "y"
# X and y are matched positionally (row i of X with row i of y). Reset X
# to a clean RangeIndex so that a non-default index on the input X (for
# example a date index, or a sliced-out subset of a larger DataFrame)
# cannot misalign with the freshly built, RangeIndex-ed ``y_full``.
# Without this, statsmodels raises "indices for endog and exog are not
# aligned" once a design matrix is handed to ``sm.OLS`` below.
X_full = X_full.reset_index(drop=True)
n_original = y_full.shape[0]
self.feature_names_in_ = [str(c) for c in X_full.columns]
self.target_name_ = str(target_name)
self.n_features_in_ = X_full.shape[1]
# Initialize all fitted attributes to a "not enough data" state.
self._init_empty_attributes(n_original)
# Missing-value handling.
missing_idx = y_full.isna().any(axis=1)
if self.na_rm:
missing_idx = y_full.isna().any(axis=1) | X_full.isna().any(axis=1)
X_ = X_full.loc[~missing_idx, :]
y_ = y_full.loc[~missing_idx]
# Degenerate input: too few rows.
if y_.size <= 1 or X_.size <= 1:
self.is_fitted_ = False
return self
n_samples, n_features = X_.shape
k_params = n_features + (1 if self.fit_intercept else 0)
if n_samples < k_params:
raise ValueError("N >= K: You need at least as many rows as there are columns to fit a linear regression.")
if n_samples != y_.size:
raise ValueError(f"X and y must have the same number of rows: got {n_samples} and {y_.size}.")
# Build the design matrix.
design = X_.copy()
if self.fit_intercept:
design.insert(0, "__constant__", 1.0)
results = sm.OLS(y_, design).fit()
alpha = 1.0 - self.conflevel
conf = np.asarray(results._results.conf_int(alpha), dtype=float)
# se^2 * (X'X)^-1 for the design (intercept first, if fitted): needed to
# form prediction intervals at arbitrary new x via prediction_interval().
self._cov_params_ = np.asarray(results.cov_params(), dtype=float)
# statsmodels returns Series (indexed by column name) when fed a DataFrame,
# so always go through ``.values`` for positional access.
params = np.asarray(results.params.values if hasattr(results.params, "values") else results.params, dtype=float)
bse = np.asarray(results.bse.values if hasattr(results.bse, "values") else results.bse, dtype=float)
tvalues = np.asarray(
results.tvalues.values if hasattr(results.tvalues, "values") else results.tvalues, dtype=float
)
pvalues = np.asarray(
results.pvalues.values if hasattr(results.pvalues, "values") else results.pvalues, dtype=float
)
# Populate attributes.
self.n_samples_ = n_samples
self.df_resid_ = int(results.df_resid)
self.df_model_ = int(results.df_model)
if self.fit_intercept:
self.intercept_ = float(params[0])
self.standard_error_intercept_ = float(bse[0])
self.t_value_intercept_ = float(tvalues[0])
self.p_value_intercept_ = float(pvalues[0])
self.conf_interval_intercept_ = conf[0, :]
self.coefficients_ = params[1:]
self.standard_errors_ = bse[1:]
self.t_values_ = tvalues[1:]
self.p_values_ = pvalues[1:]
self.conf_intervals_ = conf[1:, :]
else:
self.intercept_ = float("nan")
self.standard_error_intercept_ = float("nan")
self.t_value_intercept_ = float("nan")
self.p_value_intercept_ = float("nan")
self.conf_interval_intercept_ = np.array([np.nan, np.nan])
self.coefficients_ = params
self.standard_errors_ = bse
self.t_values_ = tvalues
self.p_values_ = pvalues
self.conf_intervals_ = conf
self.fitted_values_ = np.asarray(results.fittedvalues, dtype=float)
self.r2_ = float(results.rsquared)
self.adj_r2_ = float(results.rsquared_adj)
self.se_ = float(np.sqrt(results.scale))
self.f_statistic_ = float(results.fvalue) if results.fvalue is not None else float("nan")
self.f_pvalue_ = float(results.f_pvalue) if results.f_pvalue is not None else float("nan")
# Residuals are written into a vector of the original shape (NaN where rows were dropped).
self.residuals_ = np.full(n_original, np.nan)
self.residuals_[~missing_idx.to_numpy()] = results.resid
# R^2 reported two ways for backwards compatibility with the dict API.
mean_y = float(np.mean(y_.values))
total_ssq = float(np.sum(np.power(y_.values - mean_y, 2)))
regression_ssq = float(np.sum(np.power(self.fitted_values_ - mean_y, 2)))
residual_ssq = float(np.nansum(self.residuals_ * self.residuals_))
self.r2_regression_based_ = regression_ssq / total_ssq if total_ssq > 0 else float("nan")
self.r2_residual_based_ = 1.0 - residual_ssq / total_ssq if total_ssq > 0 else float("nan")
# Leverage and Cook's distance, for any number of predictors. Both come from
# the hat-matrix diagonal of the model matrix statsmodels fitted, so they hold
# for multiple regression and with or without an intercept. The pseudo-inverse
# keeps a rank-deficient model matrix from raising.
exog = np.asarray(results.model.exog, dtype=float)
self.leverage_ = np.einsum("ij,ji->i", exog, np.linalg.pinv(exog))
eps = np.finfo(np.float32).eps
residual_room = 1.0 - self.leverage_
self.influence_ = np.zeros_like(self.leverage_)
if self.se_ >= eps:
# A row with leverage 1 sets its own fitted value, so it has no residual to
# redistribute and Cook's distance is not defined by the usual ratio. Zero is
# the reading for that row alone; every other row keeps its own value.
# np.asarray: results.resid is a Series when y arrives as one, and the
# documented type of both diagnostics is a plain array.
resid = np.asarray(results.resid, dtype=float)
movable = np.abs(residual_room) >= eps
self.influence_[movable] = (
np.power(resid[movable] / (residual_room[movable] * self.se_), 2) * self.leverage_[movable] / k_params
)
# A prediction-interval grid and the x sum of squares only mean something for a
# single predictor with an intercept.
self.x_ssq_ = float("nan")
self.pi_range_: np.ndarray | float = float("nan")
if n_features == 1:
x_col = X_.iloc[:, 0].to_numpy()
mean_x = float(np.mean(x_col))
x_ssq = float(np.sum(np.power(x_col - mean_x, 2)))
self.x_ssq_ = x_ssq
# No condition on se_: a near-perfect fit has a band of near-zero width, which is
# its prediction interval, and `prediction_interval()` gives the same at any x.
if self.fit_intercept and x_ssq > 0:
pi_range = np.linspace(np.min(x_col), np.max(x_col), self.pi_resolution)
pi_y_pred = self.intercept_ + self.coefficients_[0] * pi_range
var_y = (self.se_**2) * (1 + 1.0 / n_samples + (pi_range - mean_x) ** 2 / x_ssq)
std_y = np.sqrt(var_y)
c_t = t_value(1 - (1 - self.conflevel) / 2, n_samples - 2)
self.pi_range_ = np.vstack([pi_range, pi_y_pred - c_t * std_y, pi_y_pred + c_t * std_y]).T
self._k_ = k_params
self.is_fitted_ = True
return self
[docs]
def predict(
self,
X: np.ndarray | pd.DataFrame | pd.Series,
) -> np.ndarray:
"""Predict target values for ``X``.
Parameters
----------
X : array-like of shape (N, K)
Returns
-------
y_pred : np.ndarray of shape (N,)
"""
if not getattr(self, "is_fitted_", False):
raise NotFittedError("OLS must be fitted before calling predict().")
if isinstance(X, pd.DataFrame):
X_arr = X.to_numpy(dtype=float)
elif isinstance(X, pd.Series):
X_arr = X.to_numpy(dtype=float).reshape(-1, 1)
else:
X_arr = np.asarray(X, dtype=float)
if X_arr.ndim == 1:
X_arr = X_arr.reshape(-1, 1)
n_features = X_arr.shape[1]
if n_features != self.n_features_in_:
raise ValueError(
f"X has {n_features} feature(s), but OLS was fitted with {self.n_features_in_} feature(s)."
)
intercept = 0.0 if not self.fit_intercept else self.intercept_
return intercept + X_arr @ self.coefficients_
[docs]
def prediction_interval(
self,
X: np.ndarray | pd.DataFrame | pd.Series | float,
conflevel: float | None = None,
) -> Bunch:
"""Prediction interval for new observations at arbitrary ``X``.
Unlike the ``pi_range_`` attribute - which is evaluated on a fixed grid
spanning the training data - this method evaluates the prediction
interval at any predictor value(s) supplied by the caller, including
points outside the training range.
Parameters
----------
X : array-like of shape (M, K), (K,) or scalar
New predictor value(s). A scalar or 1-D array is interpreted as a
list of points when the model has a single feature, or as a single
multi-feature point otherwise.
conflevel : float or None, default=None
Confidence level for the interval. Defaults to the model's own
``conflevel``.
Returns
-------
sklearn.utils.Bunch
A bunch with three length-M arrays: ``predicted`` (the point
prediction), and ``lower`` / ``upper`` (the prediction-interval
bounds).
"""
assert getattr(self, "is_fitted_", False), "OLS must be fitted before calling prediction_interval()."
cl = self.conflevel if conflevel is None else conflevel
_check_probability("conflevel", cl)
X_arr = X.to_numpy(dtype=float) if isinstance(X, pd.Series | pd.DataFrame) else np.asarray(X, dtype=float)
if X_arr.ndim == 0:
X_arr = X_arr.reshape(1, 1)
elif X_arr.ndim == 1:
# Single-feature model: a 1-D input is a list of scalar points.
# Otherwise it is interpreted as one multi-feature point.
X_arr = X_arr.reshape(-1, 1) if self.n_features_in_ == 1 else X_arr.reshape(1, -1)
if X_arr.shape[1] != self.n_features_in_:
msg = f"Expected {self.n_features_in_} feature(s) per row, got {X_arr.shape[1]}."
raise ValueError(msg)
# Design rows must match the fitted design: intercept column first.
if self.fit_intercept:
design = np.column_stack([np.ones(X_arr.shape[0]), X_arr])
full_params = np.concatenate([[self.intercept_], self.coefficients_])
else:
design = X_arr
full_params = self.coefficients_
predicted = design @ full_params
# Variance of the estimated mean response: row @ (se^2 (X'X)^-1) @ row.
# The prediction interval adds se^2 for the new observation's own noise.
var_mean = np.einsum("ij,jk,ik->i", design, self._cov_params_, design)
std_pred = np.sqrt(self.se_**2 + var_mean)
c_t = t_value(1 - (1 - cl) / 2, self.df_resid_)
half_width = c_t * std_pred
return Bunch(predicted=predicted, lower=predicted - half_width, upper=predicted + half_width)
[docs]
def summary(self) -> str:
"""Return an R-style ``summary(lm(...))`` string for the fitted model."""
if not getattr(self, "is_fitted_", False):
return "OLS model has not been fitted (or fit failed due to insufficient data)."
return self._format_summary()
[docs]
def to_dict(self) -> dict:
"""Return the legacy dictionary representation used by ``multiple_linear_regression``."""
out: dict[str, Any] = {
"N": None,
"coefficients": [np.nan],
"intercept": np.nan,
"standard_errors": [np.nan],
"standard_error_intercept": np.nan,
"R2": np.nan,
"SE": np.nan,
"fitted_values": np.nan,
"residuals": np.nan,
"t_value": np.nan,
"conf_intervals": [np.nan, np.nan],
}
if not getattr(self, "is_fitted_", False):
return out
out["N"] = self.n_samples_
out["intercept"] = self.intercept_
out["coefficients"] = self.coefficients_
out["standard_errors"] = self.standard_errors_
out["standard_error_intercept"] = self.standard_error_intercept_
out["t_value"] = self.t_values_
out["conf_intervals"] = self.conf_intervals_
out["conf_interval_intercept"] = self.conf_interval_intercept_
out["R2"] = self.r2_
out["SE"] = self.se_
out["fitted_values"] = self.fitted_values_
out["residuals"] = self.residuals_
out["R2_regression_based"] = self.r2_regression_based_
out["R2_residual_based"] = self.r2_residual_based_
out["k"] = self._k_
if not np.isnan(self.x_ssq_):
out["x_ssq"] = self.x_ssq_
# A fitted model has at least two rows, so these are always real per-row arrays,
# never the one-element NaN placeholders an unfitted model holds.
out["leverage"] = self.leverage_
out["influence"] = self.influence_
if isinstance(self.pi_range_, np.ndarray):
out["pi_range"] = self.pi_range_
return out
def __repr__(self) -> str:
"""Return the R-style summary if fitted, otherwise the sklearn parameter repr."""
if not getattr(self, "is_fitted_", False):
return super().__repr__()
return self._format_summary()
# ------------------------------------------------------------------
# Internal helpers
# ------------------------------------------------------------------
def _init_empty_attributes(self, n_original: int) -> None:
"""Set all fitted attributes to NaN/empty so attribute access is safe even when fit fails."""
self.is_fitted_ = False
self.n_samples_ = 0
self.df_resid_ = 0
self.df_model_ = 0
self.intercept_ = float("nan")
self.standard_error_intercept_ = float("nan")
self.t_value_intercept_ = float("nan")
self.p_value_intercept_ = float("nan")
self.conf_interval_intercept_ = np.array([np.nan, np.nan])
self.coefficients_ = np.array([np.nan])
self.standard_errors_ = np.array([np.nan])
self.t_values_ = np.array([np.nan])
self.p_values_ = np.array([np.nan])
self.conf_intervals_ = np.array([[np.nan, np.nan]])
self.fitted_values_ = np.full(n_original, np.nan)
self.residuals_ = np.full(n_original, np.nan)
self.r2_ = float("nan")
self.adj_r2_ = float("nan")
self.se_ = float("nan")
self.f_statistic_ = float("nan")
self.f_pvalue_ = float("nan")
self.r2_regression_based_ = float("nan")
self.r2_residual_based_ = float("nan")
self.x_ssq_ = float("nan")
self.leverage_ = np.array([np.nan])
self.influence_ = np.array([np.nan])
self.pi_range_ = float("nan")
self._cov_params_ = np.array([[np.nan]])
self._k_ = 0
@staticmethod
def _signif_code(p: float) -> str:
if np.isnan(p):
return " "
if p < 0.001:
return "***"
if p < 0.01:
return "** "
if p < 0.05:
return "* "
if p < 0.1:
return ". "
return " "
@staticmethod
def _fmt_number(value: float, width: int = 10) -> str:
if np.isnan(value):
return f"{'NaN':>{width}}"
if value == 0:
return f"{0.0:>{width}.5f}"
magnitude = abs(value)
if magnitude < 1e-3 or magnitude >= 1e6:
return f"{value:>{width}.4e}"
return f"{value:>{width}.5f}"
@staticmethod
def _fmt_pvalue(p: float, width: int = 10) -> str:
if np.isnan(p):
return f"{'NaN':>{width}}"
if p < 1e-4:
return f"{p:>{width}.3e}"
return f"{p:>{width}.5f}"
def _format_summary(self) -> str:
# Call line: show the parameters that affected the fit.
call = f"OLS(fit_intercept={self.fit_intercept}, na_rm={self.na_rm}, conflevel={self.conflevel})"
target = self.target_name_ or "y"
feature_str = " + ".join(self.feature_names_in_) if self.feature_names_in_ else "1"
formula_line = (
f" formula: {target} ~ {feature_str}" if self.fit_intercept else f" formula: {target} ~ {feature_str} + 0"
)
# Residual quantiles.
resid = self.residuals_[~np.isnan(self.residuals_)]
q = np.quantile(resid, [0.0, 0.25, 0.5, 0.75, 1.0]) if resid.size else np.array([np.nan] * 5)
resid_header = f"{'Min':>10}{'1Q':>11}{'Median':>11}{'3Q':>11}{'Max':>11}"
resid_values = (
f"{self._fmt_number(q[0]):>10}"
f"{self._fmt_number(q[1]):>11}"
f"{self._fmt_number(q[2]):>11}"
f"{self._fmt_number(q[3]):>11}"
f"{self._fmt_number(q[4]):>11}"
)
# Coefficient table.
rows: list[tuple[str, float, float, float, float]] = []
if self.fit_intercept:
rows.append(
(
"(Intercept)",
self.intercept_,
self.standard_error_intercept_,
self.t_value_intercept_,
self.p_value_intercept_,
)
)
for name, est, se, tv, pv in zip(
self.feature_names_in_,
self.coefficients_,
self.standard_errors_,
self.t_values_,
self.p_values_,
strict=False,
):
rows.append((str(name), float(est), float(se), float(tv), float(pv)))
name_width = max((len(r[0]) for r in rows), default=11)
name_width = max(name_width, 11)
coef_header = f"{'':<{name_width}} {'Estimate':>10} {'Std. Error':>10} {'t value':>8} {'Pr(>|t|)':>10}"
coef_lines = []
for name, est, se, tv, pv in rows:
line = (
f"{name:<{name_width}} "
f"{self._fmt_number(est):>10} "
f"{self._fmt_number(se):>10} "
f"{self._fmt_number(tv, width=8):>8} "
f"{self._fmt_pvalue(pv):>10} {self._signif_code(pv)}"
)
coef_lines.append(line)
# F-statistic line.
if np.isnan(self.f_statistic_) or np.isnan(self.f_pvalue_):
f_line = "F-statistic: not available"
else:
f_line = (
f"F-statistic: {self.f_statistic_:.4g} on "
f"{self.df_model_} and {self.df_resid_} DF, "
f"p-value: {self.f_pvalue_:.4g}"
)
return (
"Call:\n"
f" {call}\n"
f"{formula_line}\n"
"\n"
"Residuals:\n"
f"{resid_header}\n"
f"{resid_values}\n"
"\n"
"Coefficients:\n"
f"{coef_header}\n" + "\n".join(coef_lines) + "\n"
"---\n"
"Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n"
"\n"
f"Residual standard error: {self.se_:.4g} on {self.df_resid_} degrees of freedom\n"
f"Multiple R-squared: {self.r2_:.4g},\tAdjusted R-squared: {self.adj_r2_:.4g}\n"
f"{f_line}"
)
# ENG-23 (#305): explicit ``__all__`` so the thin re-exporter ``methods.py``
# can do ``from ._robust_regression import *`` without triggering CodeQL's
# py/polluting-import warning.
__all__ = [
"OLS",
"fit_robust_lm",
"multiple_linear_regression",
"repeated_median_slope",
"robust_regression",
"t_value",
]
