Regression#

Backwards-compatible re-exporter for process_improve.regression.

The implementation now lives in process_improve.regression._robust_regression (ENG-23 / #305): the renamed file makes filename-ranked tooling (Jump-to-File, fuzzy search, codecov reports) less ambiguous about which methods.py is being shown.

Every public name remains importable as before:

from process_improve.regression.methods import OLS, robust_regression, repeated_median_slope

class process_improve.regression.methods.OLS(fit_intercept=True, na_rm=True, conflevel=0.95, pi_resolution=50)[source]#

Bases: RegressorMixin, BaseEstimator

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.

coefficients_#

Fitted slope coefficients (excludes the intercept).

Type:: np.ndarray of shape (K,)

intercept_#

Fitted intercept (np.nan if fit_intercept is False).

Type:: float

standard_errors_#

Standard errors of coefficients_.

Type:: np.ndarray of shape (K,)

standard_error_intercept_#

Standard error of the intercept.

Type:: float

t_values_#

t-statistics for each coefficient.

Type:: np.ndarray of shape (K,)

t_value_intercept_#

t-statistic for the intercept.

Type:: float

p_values_#

Two-sided p-values for each coefficient.

Type:: np.ndarray of shape (K,)

p_value_intercept_#

p-value for the intercept.

Type:: float

conf_intervals_#

Lower and upper bounds of the coefficient confidence intervals.

Type:: np.ndarray of shape (K, 2)

conf_interval_intercept_#

Lower and upper bounds of the intercept confidence interval.

Type:: np.ndarray of shape (2,)

r2_#

Coefficient of determination.

Type:: float

adj_r2_#

Adjusted R-squared.

Type:: float

se_#

Residual standard error (sqrt of residual variance).

Type:: float

df_resid_#

Residual degrees of freedom.

Type:: int

df_model_#

Model degrees of freedom (number of slope coefficients).

Type:: int

f_statistic_#

F-statistic for the overall regression.

Type:: float

f_pvalue_#

p-value associated with the F-statistic.

Type:: float

fitted_values_#

In-sample predictions.

Type:: np.ndarray of shape (N,)

residuals_#

In-sample residuals (NaN at rows removed by na_rm).

Type:: np.ndarray of shape (N_original,)

leverage_#

Hat-matrix diagonal $h_i = \mathbf{x}_i^T (X^T X)^{-1} \mathbf{x}_i$, over the rows that were fitted. Matches statsmodels hat_matrix_diag.

Type:: np.ndarray of shape (N,)

influence_#

Cook’s distance per fitted row. Matches statsmodels cooks_distance.

Type:: np.ndarray of shape (N,)

pi_range_#

Columns are x-grid, lower bound, upper bound of the prediction interval. np.nan if not applicable.

Type:: np.ndarray of shape (pi_resolution, 3) or float

feature_names_in_#

Column names of the feature matrix.

Type:: list[str]

target_name_#

Name of the target variable.

Type:: str

n_samples_#

Number of samples used in the fit (after na_rm).

Type:: int

n_features_in_#

Number of input features.

Type:: int

is_fitted_#

Whether fit() has been called successfully.

Type:: bool

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

fit(X, y)[source]#

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 – Fitted estimator.

Return type:

OLS

predict(X)[source]#

Predict target values for X.

Parameters:: X (array-like of shape (N, K))
Returns:: y_pred
Return type:: np.ndarray of shape (N,)

prediction_interval(X, conflevel=None)[source]#

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:

A bunch with three length-M arrays: predicted (the point prediction), and lower / upper (the prediction-interval bounds).

Return type:

sklearn.utils.Bunch

summary()[source]#

Return an R-style summary(lm(...)) string for the fitted model.

Return type:: str

to_dict()[source]#

Return the legacy dictionary representation used by multiple_linear_regression.

Return type:: dict

set_score_request(*, sample_weight='$UNCHANGED$')#

Configure whether metadata should be requested to be passed to the score method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.
False: metadata is not requested and the meta-estimator will not pass it to score.
None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.
str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters:

sample_weight (str, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED) – Metadata routing for sample_weight parameter in score.
self (OLS)

Returns:

self – The updated object.

Return type:

object

process_improve.regression.methods.fit_robust_lm(x, y)[source]#

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

Parameters:

x (ndarray)
y (ndarray)

Return type:

ndarray

process_improve.regression.methods.multiple_linear_regression(X, y, fit_intercept=True, na_rm=True, conflevel=0.95, pi_resolution=50)[source]#

Linear regression of the N rows and K columns of matrix X onto the single column ‘y’.

Backwards-compatible wrapper around OLS. New code should use the 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)

Parameters:

X (ndarray | DataFrame)
y (ndarray | DataFrame | Series)
fit_intercept (bool)
na_rm (bool)
conflevel (float)
pi_resolution (int)

Return type:

dict

process_improve.regression.methods.repeated_median_slope(x, y, nowarn=False)[source]#

Robust slope calculation via Siegel’s repeated-median estimator.

https://en.wikipedia.org/wiki/Repeated_median_regression

An elegant (simple) method to compute the robust slope between a vector x and y. For each point i the median of the pairwise slopes (y[j] - y[i]) / (x[j] - x[i]) over all j != i is computed; the returned slope is the median of those per-point medians.

Parameters:

x (np.ndarray or sequence) – Independent variable. Coerced to a 1-D numpy array. Must have at least 3 elements (unless nowarn=True).
y (np.ndarray or sequence) – Dependent variable. Must have the same length as x.
nowarn (bool, optional) – If True, skip the minimum-length check and the SEC-19 cap on the number of points. Default False. The equal-length check and the finite-data check are always applied.

Returns:

The repeated-median estimate of the slope. Returns np.nan if all inner medians are undefined (e.g. all x values are equal).

Return type:

float

Raises:

ValueError – If x and y differ in length; if, unless nowarn=True, there are two points or fewer or more than the SEC-19 cap; or if fewer than three finite (x, y) pairs remain after non-finite pairs are dropped, so no repeated median can be formed.

Notes

Missing data: a pair is used only when both x[i] and y[i] are finite. Pairs holding a NaN or an infinity in either vector are dropped before any slope is formed, and the estimate is computed from the remainder. The dropping is pairwise and silent, so check the inputs first if the count of usable observations matters to you.

Cost: the pairwise slopes are formed at once as an N x N matrix, so time and memory grow as N squared. At the SEC-19 cap of 5,000 points the call needs about 425 MB of memory; at 1,000 points, 17 MB. Faster algorithms exist (Matousek, Mount, Netanyahu, 1998, https://link.springer.com/article/10.1007/PL00009190), but they are not used here.

process_improve.regression.methods.robust_regression(x, y, fit_intercept=True, na_rm=True, conflevel=0.95, nowarn=False, pi_resolution=50, *, intervals=None, n_boot=999, random_state=0)[source]#

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 (np.ndarray, pd.DataFrame or pd.Series) – Sequences of numerical values, paired by position.
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 process_improve._random.check_random_state(). The default seed makes a repeated call return the same intervals. Used only with intervals="bootstrap".

Returns:

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.

Return type:

dict

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 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".

process_improve.regression.methods.t_value(p, v)[source]#

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 $-\infty$:

>>> t_value(0.0, v=10)
-inf

100% fractional area is always at $+\infty$:

>>> t_value(1.0, v=10)
inf

Power of a regressor#

Box and Tidwell’s estimate of the power of a regressor.

A response that is curved in a regressor $x$ can be straight in a power of it, $x^\alpha$, with $\ln x$ in the place of $\alpha = 0$. Box and Tidwell (Technometrics, 1962) estimate the powers of the model

\[y = \beta_0 + \sum_j \beta_j x_j^{\alpha_j} + (\text{regressors that enter as they are}) + \varepsilon\]

with ordinary regressions only. Near a trial power $\alpha$,

\[x^{\alpha'} \approx x^{\alpha} + (\alpha' - \alpha)\, x^{\alpha} \ln x,\]

so adding the constructed variable $x^{\alpha}\ln x$ to a regression on $x^{\alpha}$ gives it the coefficient $\gamma = \beta(\alpha' - \alpha)$. Each step:

regresses $y$ on the regressors at the trial powers, which gives each $\hat\beta$;
adds the constructed variables, which gives each $\hat\gamma$;
moves each power to $\alpha + \hat\gamma / \hat\beta$.

From $\alpha = 1$, the regressor as it is, the first step is Box and Tidwell’s $\alpha_1 = 1 + \hat\gamma/\hat\beta$. Repeated, the steps are the Gauss-Newton method for the least-squares estimate of the powers; here a step is halved whenever it would raise the residual sum of squares. The standard error of each power comes from the same linearization at the estimate, as for any nonlinear least-squares parameter.

The constructed variable’s t value at $\alpha = 1$ tests whether the regressor needs a power other than 1 at all: it is the score test of a straight-line relation in that regressor.

Examples

A response that is straight in $\sqrt{x}$:

>>> import numpy as np
>>> x = np.arange(1.0, 11.0)
>>> found = box_tidwell(x, 2 + 3 * np.sqrt(x))
>>> round(found["power"]["x"], 6), found["df_resid"]
(0.5, 7)

process_improve.regression.box_tidwell.SERIES = 1e-06#: Below this size a power’s columns are summed from their series in it, exact to rounding there.

process_improve.regression.box_tidwell.HALVINGS = 40#: The most times a step is halved in search of a lower residual sum of squares.

process_improve.regression.box_tidwell.FLAT = 1.4901161193847656e-08#: A regressor whose straight-line part of the fit is below this share of the spread of y has no slope at its trial power: its power cannot be moved, whatever rounding leaves in its coefficient.

process_improve.regression.box_tidwell.box_tidwell(X, y, columns=None, max_iter=100, tol=1e-08)[source]#

Estimate the power of one or more regressors by Box and Tidwell’s method.

Parameters:

X (pandas.DataFrame, numpy.ndarray or Sequence[float]) – The regressors, a column each. A one-dimensional array is one regressor, named "x"; the columns of a two-dimensional array are named "x1", "x2", and so on.
y (numpy.ndarray or Sequence[float]) – The response.
columns (Sequence[str] or None) – The regressors whose power is estimated; each must be positive and take three or more different values. None estimates a power for every column. The other columns enter the model as they are, beside the intercept.
max_iter (int) – The most steps to take.
tol (float) – The steps stop when no power moves by more than this.

Returns:

power: the least-squares estimate of each power, by column name;
std_error: its standard error, $s$ times the square root of the power’s diagonal element of $(J^\top J)^{-1}$, with $J$ the derivatives of the fitted values in every parameter at the estimate;
linearity_t and linearity_p_value: the t value of each constructed variable in the regression with every power at 1, and its two-sided p value; a small p value says the relation is not a straight line in that regressor;
sse, the residual sum of squares at the estimate, and sse_linear, with every power at 1;
df_resid, the residual degrees of freedom: the rows less the intercept, the regressors’ coefficients and the powers; n, the rows;
iterations, the powers at the start (all 1) and after each step: the first step, unless it had to be halved, is Box and Tidwell’s $\alpha_1$.

Return type:

dict[str, Any]

Raises:

ValueError – If the data or the arguments cannot give a power: a name that is not a column, a column to transform with a value that is not positive or fewer than three different values, a value that is not a finite number, a response of another length or with no variation, too few rows to leave a residual degree of freedom, a regressor with no clear straight-line effect on the response, or powers that do not settle in max_iter steps. Powers that keep growing in size mean that no finite power fits best: the fit improves the more a regressor is reduced to its largest or smallest value.

Nonlinear least squares#

Least squares for a model that is nonlinear in its parameters, with the usual inference.

The model $y = f(x; \theta) + \varepsilon$ is fitted by minimizing the residual sum of squares $S(\theta) = \sum_u (y_u - f(x_u; \theta))^2$. Inference rests on the linearization of $f$ at the estimate $\hat\theta$, the matrix $J$ of the derivatives of the fitted values in the $p$ parameters (Box, Hunter and Hunter, section 14.6; Bates and Watts):

the covariance of the estimates is $s^2 (J^\top J)^{-1}$, with $s^2 = S(\hat\theta)/(n - p)$, which gives standard errors, t values and intervals;
the approximate joint confidence region for all the parameters is $\{\theta : S(\theta) \le S(\hat\theta)\,[1 + \tfrac{p}{n - p} F_{p,\,n-p}]\}$. The region follows the contours of $S$ itself, so it keeps the curved shape of the valley that the linearized intervals miss.

Examples

Box, Hunter and Hunter’s biochemical oxygen demand against incubation time, with $\eta = \beta_1 (1 - e^{-\beta_2 x})$:

>>> import numpy as np
>>> days = np.array([1, 2, 3, 5, 7, 10])
>>> demand = np.array([109, 149, 149, 191, 213, 224])
>>> fit = nonlinear_regression(lambda x, b1, b2: b1 * (1 - np.exp(-b2 * x)), days, demand, {"b1": 200, "b2": 0.5})
>>> round(fit["estimates"]["b1"], 1), round(fit["estimates"]["b2"], 4), round(fit["sse"], 2)
(213.8, 0.5472, 1168.01)

process_improve.regression.nonlinear.nonlinear_regression(model, x, y, start, conf_level=0.95)[source]#

Fit y = model(x, *parameters) by least squares, and give the inference of the linearized model.

Parameters:

model (Callable[..., Any]) – Called as model(x, theta_1, ..., theta_p), with the parameters in the order of start; returns the fitted value of each row of x.
x (numpy.ndarray, pandas.DataFrame or Sequence[float]) – The regressors, passed to model as they are.
y (numpy.ndarray or Sequence[float]) – The response, one value per row of x.
start (Mapping[str, float]) – The parameters’ names and their starting values, from which the search begins. A nonlinear model can have more than one minimum: start near the expected values.
conf_level (float) – The confidence level of the intervals and of the joint region.

Returns:

estimates, std_errors, t_values, p_values (two-sided, of each parameter against 0) and intervals (the estimate plus or minus $t_{n-p}$ standard errors), each by parameter name;
correlation: the correlations of the estimates, a nested dict;
sse ($S(\hat\theta)$), df_resid ($n - p$), s2, n;
fitted and residuals, in the order of the rows;
region_sse: $S(\hat\theta)\,[1 + \tfrac{p}{n - p} F_{p,\,n-p}]$, the contour of the residual sum of squares that bounds the joint confidence region.

Return type:

dict[str, Any]

Raises:

ValueError – If start names no parameter or gives one that is not a finite number, y has a value that is not a finite number, there are no more rows than parameters, model does not return a finite value for each row at the start, conf_level is not between 0 and 1, the search does not converge, or the estimates cannot be told apart (the derivatives of the fitted values in the parameters are linearly dependent at the estimate).