Regression

Linear and generalized linear models. OLS via QR (CPU) or Cholesky (GPU). GLM via IRLS with Gaussian, Binomial, and Poisson families.

Linear and generalized linear models.

Usage:

from pystatistics import DataSource from pystatistics.regression import Design, fit

# OLS from DataSource ds = DataSource.from_file(“data.csv”) design = Design.from_datasource(ds, y=’target’) result = fit(design)

# OLS from arrays (convenience) result = fit(X, y)

# GLM (logistic regression) result = fit(X, y, family=’binomial’)

# GLM (Poisson regression) result = fit(X, y, family=’poisson’)

pystatistics.regression.fit(X_or_design, y=None, *, family=None, backend=None, solver=None, force=False, tol=1e-08, max_iter=25, names=None, l2=0.0, weights=None, offset=None, conf_level=0.95)[source]

Fit a linear or generalized linear model.

When family is None (default), fits ordinary least squares (LM) via QR decomposition or GPU Cholesky. When family is specified, fits a GLM via IRLS (Iteratively Reweighted Least Squares).

Accepts EITHER:

  1. A Design object (from DataSource or arrays)

  2. Raw X and y arrays (convenience)

Parameters:

  • X_or_design (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | Design) – Design object or X matrix

  • y (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – Response vector (required if X_or_design is array)

  • family (str | Family | None) – GLM family specification. None for OLS, or a string (‘gaussian’, ‘binomial’, ‘poisson’) or Family instance.

  • backend (Literal['auto', 'cpu', 'gpu', 'gpu_fp64'] | None) – Compute backend = (device, precision). Default None → ‘cpu’ (the R-reference path, validated for regulated-industry use). Values: ‘cpu’ (float64), ‘gpu’ (float32, CUDA or MPS), ‘gpu_fp64’ (float64, CUDA only), or ‘auto’ (GPU-fp32 if CUDA present, else CPU).

  • solver (Literal['qr', 'svd'] | None) – Numerical routine for the linear-model fit (family=None only): ‘qr’ (default) or ‘svd’. Not configurable on the GPU backend (which uses Cholesky on the normal equations) and not applicable to GLMs.

  • force (bool) – If True, proceed with the GPU float32 Cholesky path even when it is unreliable — for OLS, on ill-conditioned designs; for GLM, when IRLS does not converge in float32 (returns the possibly-inaccurate fit instead of raising). Has no effect on CPU backends.

  • tol (float) – Convergence tolerance for IRLS (GLM only). Default 1e-8 matches R’s glm.control().

  • max_iter (int) – Maximum IRLS iterations (GLM only). Default 25 matches R’s glm.control().

  • names (list[str] | None) – Optional list of predictor names for output labeling. If len(names) == p - 1 (one fewer than columns in X), “(Intercept)” is prepended automatically. If len(names) == p, used as-is.

  • l2 (float) – L2 (ridge) penalty strength. Default 0.0 (unpenalized). When > 0, fits a ridge-penalized model: predictors are standardized, the intercept is left unpenalized, and l2 is the penalty added on the standardized scale (matching MASS::lm.ridge’s lambda). A penalized fit does not report standard errors / t / p values (not valid for a biased estimator). See the ridge() convenience wrapper.

  • weights (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – Per-observation prior weights (n,), matching R’s lm(..., weights=) / glm(..., weights=). For OLS this is weighted least squares; for a GLM these are the IRLS prior weights. Must be non-negative and not all zero. None ⇒ unit weights. Not supported together with l2 > 0 (raises).

  • offset (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – Additive term in the linear predictor, η = Xβ + offset (n,), matching R’s glm(..., offset=). Used as-is, never estimated — e.g. log(exposure) for a Poisson rate model. None ⇒ no offset. Not supported together with l2 > 0 (raises).

  • conf_level (float)

Returns:

LinearSolution when family is None. GLMSolution when family is specified.

Return type:

LinearSolution | GLMSolution

Examples

# OLS with named output >>> result = fit(X, y, names=[‘albumin’, ‘copper’, ‘protime’]) >>> result.coef[‘copper’] 0.005255

# Logistic regression >>> result = fit(X, y, family=’binomial’)

# Poisson regression >>> result = fit(X, y, family=’poisson’)

pystatistics.regression.ridge(X_or_design, y=None, *, lam, family=None, backend=None, tol=1e-08, max_iter=25, names=None, weights=None, offset=None)[source]

Ridge (L2-penalized) regression — a thin wrapper over fit(..., l2=lam).

ridge(X, y, lam=λ) fits an L2-penalized linear model; pass family= for a penalized GLM (e.g. logistic ridge). Predictors are standardized and the intercept is unpenalized; lam is the penalty on the standardized scale, matching MASS::lm.ridge. Penalized fits do not report standard errors.

Equivalent to fit(X_or_design, y, family=family, l2=lam, ...); provided for discoverability. (Future: lasso / elastic_net wrappers will sit alongside this once the coordinate-descent solver exists.)

Parameters:
Return type:

LinearSolution | GLMSolution

class pystatistics.regression.Design(_X, _y, _n, _p, _source=None, _names=None)[source]

Bases: object

Regression design matrix specification.

Wraps a DataSource and provides X, y for regression. Immutable after construction.

Construction:

Design.from_datasource(ds, y=’target’) # X = all other columns Design.from_datasource(ds, x=[‘a’,’b’], y=’c’) # X = specified columns Design.from_datasource(ds) # Uses ds[‘X’] and ds[‘y’] Design.from_arrays(X, y) # Direct from arrays

Parameters:
classmethod from_datasource(source, *, x=None, y=None, terms=None)[source]

Build Design from DataSource.

Parameters:

  • source (DataSource) – The DataSource

  • x (str | list[str] | None) – Predictor column(s). If None and source has ‘X’, uses that. If None and y is specified, uses all columns except y.

  • y (str | None) – Response column. If None, uses ‘y’ from source.

  • terms (Sequence[Any] | None) – Structured term spec for categorical predictors and interactions (mutually exclusive with x). A list whose elements are bare column names (numeric main effects), C(name, ref=…) markers (categorical main effects), or tuples of those (interactions). When given, an intercept column is added automatically and the expanded column labels travel with the Design (see names). Requires y.

Returns:

Design ready for regression

Return type:

Design

Assumes good faith: garbage in, garbage out.

classmethod from_arrays(X, y)[source]

Build Design directly from arrays.

Parameters:
Return type:

Design

property X: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Design matrix (n x p).

property y: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Response vector (n,).

property n: int

Number of observations.

property p: int

Number of predictors.

property source: DataSource | None

Original DataSource, if available.

property names: tuple[str, ...] | None

Column labels aligned with X, when built from a term spec.

None for designs built from plain arrays/columns (the caller supplies names to fit() in that case).

supports(capability)[source]

Check if underlying data supports a capability.

Parameters:

capability (str)

Return type:

bool

XtX()[source]

Compute X’X (for standard errors).

Return type:

ndarray[tuple[Any, …], dtype[floating[Any]]]

Xty()[source]

Compute X’y.

Return type:

ndarray[tuple[Any, …], dtype[floating[Any]]]

class pystatistics.regression.C(name, ref=None)[source]

Bases: object

Mark a column as a categorical predictor (treatment / dummy coded).

Parameters:

  • name (str) – the column name in the DataSource.

  • ref (str | None) – the reference (baseline) level that is dropped. If None, the first level in sorted order is used as the baseline — matching R’s default factor ordering and the library’s ANOVA encoder.

A factor with k levels expands to k - 1 indicator columns labeled name[level] for each non-reference level.

name: str
ref: str | None = None
class pystatistics.regression.DataSource(_data, _capabilities, _metadata=<factory>)[source]

Bases: object

Universal data container. Domain-agnostic.

Construct via factory classmethods, not directly.

The lumber yard analogy: DataSource has data (logs). It doesn’t know or care what you’re building—furniture (regression), paper (MVN MLE), or two-by-fours (survival analysis).

Parameters:
keys()[source]

Return the names of all available arrays.

Returns:

frozenset of array names

Return type:

frozenset[str]

Example

>>> ds = DataSource.from_arrays(X=X, y=y)
>>> ds.keys()
frozenset({'X', 'y'})
property n_observations: int

Number of statistical units (rows).

property metadata: dict[str, Any]

Domain-agnostic metadata.

supports(capability)[source]

Check if this DataSource supports a capability.

Parameters:

capability (str) – Use constants from pystatistics.core.capabilities

Returns:

True if supported, False otherwise

Return type:

bool

Note

Unknown capabilities return False, never raise.

classmethod from_arrays(*, X=None, y=None, data=None, columns=None, **named_arrays)[source]

Construct from NumPy arrays.

Parameters:
Return type:

DataSource

classmethod from_file(path, *, columns=None)[source]

Construct from file (CSV, NPY).

Parameters:
Return type:

DataSource

classmethod from_dataframe(df, *, source_path=None)[source]

Construct from pandas DataFrame.

Numeric columns are stored as float64. Non-numeric columns (strings, objects, pandas categoricals) are preserved as-is rather than being force-cast — they are the raw material for categorical predictors, which the regression term builder encodes via C(…). Force-casting them to float would either crash or silently corrupt the data.

Parameters:

  • df (pd.DataFrame)

  • source_path (str | None)

Return type:

DataSource

classmethod from_tensors(*, X=None, y=None, **named_tensors)[source]

Construct from PyTorch tensors (already on GPU).

Parameters:

  • X (torch.Tensor | None)

  • y (torch.Tensor | None)

  • named_tensors (torch.Tensor)

Return type:

DataSource

to(device)[source]

Return a new DataSource with all arrays on the specified device.

Transfers the underlying materialized arrays to the given compute device ('cpu', 'cuda', 'cuda:0', 'mps', …) and returns a new DataSource. The original is unchanged — DataSources are immutable.

Intended workflow: pay the host↔device transfer once, reuse the resulting DataSource across many fits:

ds = DataSource.from_arrays(X=X, y=y)
gds = ds.to("cuda")                       # pay transfer once
pca(gds['X'], backend="gpu")              # no transfer
multinom(gds['y'], gds['X'], backend="gpu")  # no transfer

Without this, a stateless per-call API re-transfers X from host memory on every fit. Measured on a 1M × 100 FP32 matrix on an RTX 5070 Ti: per-call pageable H2D ≈ 66 ms (92% of total PCA wall time, which is only ~5 ms of actual compute). After .to("cuda"), each subsequent fit sees the 5 ms ceiling.

This method returns a NEW DataSource (Rule 5: no hidden state). The original CPU DataSource is untouched and still usable — for example, the 'cpu' backend can continue to operate on it while a sibling GPU DataSource drives the 'gpu' backend.

Parameters:

device (str) – PyTorch device string. Typical values: 'cpu', 'cuda', 'cuda:0', 'mps'. Any string torch.device() accepts is valid.

Returns:

A new DataSource whose arrays are torch.Tensor instances on the requested device (or numpy.ndarray instances if device='cpu'). Scalar metadata and array keys are preserved.

Raises:

  • ValidationError – If the DataSource is not materialized (streaming sources cannot be snapshotted to a device).

  • RuntimeError – If the requested device is unavailable (e.g. CUDA requested but no GPU present).

Return type:

DataSource

property device: str

Current compute device. 'cpu' for numpy-backed sources.

classmethod build(*args, **kwargs)[source]

Convenience factory that dispatches to appropriate from_* method.

Examples

DataSource.build(X=X, y=y) # from_arrays DataSource.build(“data.csv”) # from_file

Return type:

DataSource

class pystatistics.regression.LinearSolution(_result, _design, _names=None, _conf_level=0.95, _standard_errors=None, _t_statistics=None, _p_values=None)[source]

Bases: SolutionReprMixin

User-facing regression results.

Wraps the backend Result and provides convenient accessors for all regression outputs including standard errors and t-statistics.

Parameters:
property coefficients: ndarray[tuple[Any, ...], dtype[floating[Any]]]
property coef: dict[str, float]

Named coefficient mapping (like R’s coef() or statsmodels .params).

Returns dict mapping variable names to coefficient values. Falls back to ‘B[0]’, ‘B[1]’… when names are not available.

property residuals: ndarray[tuple[Any, ...], dtype[floating[Any]]]
property fitted_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]
property rss: float
property tss: float
property r_squared: float
property adjusted_r_squared: float
property residual_std_error: float
property standard_errors: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Standard errors of coefficients.

For CPU QR backend: uses R^{-1} from QR decomposition to compute (X’X)^{-1} = R^{-1} R^{-T}, matching R’s backsolve(R, diag(p)) exactly.

For GPU backend (no R available): uses np.linalg.inv(X’X). # NOT A FALLBACK: mathematically equivalent to QR path, # just a different computation route since GPU doesn’t store QR factors.

For rank-deficient matrices, aliased coefficients get NaN standard errors.

property t_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]

t-statistics for coefficients.

property p_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Two-sided p-values for coefficient t-tests.

property conf_level: float

Confidence level for conf_int (default 0.95).

property conf_int: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Wald confidence intervals for the coefficients, shape (p, 2).

coef ± t * se using the Student-t quantile at df_residual (the finite-sample reference for OLS, matching R’s confint.lm). Penalized (ridge) fits have NaN standard errors, so their intervals are NaN — a biased estimator has no valid Wald interval.

property rank: int
property df_residual: int
property info: dict[str, Any]
property timing: dict[str, float] | None
property backend_name: str
property warnings: tuple[str, ...]
summary()[source]

Generate R-style summary output.

Return type:

str

class pystatistics.regression.LinearParams(coefficients, residuals, fitted_values, rss, tss, rank, df_residual)[source]

Bases: object

Parameter payload for linear regression.

This is the immutable data computed by backends.

Parameters:
coefficients: ndarray[tuple[Any, ...], dtype[floating[Any]]]
residuals: ndarray[tuple[Any, ...], dtype[floating[Any]]]
fitted_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]
rss: float
tss: float
rank: int
df_residual: int
class pystatistics.regression.GLMSolution(_result, _design, _names=None, _conf_level=0.95, _standard_errors=None, _test_statistics=None, _p_values=None)[source]

Bases: SolutionReprMixin

User-facing GLM results.

Wraps the backend Result and provides convenient accessors for all GLM-specific outputs including deviance, AIC, and multiple residual types.

Parameters:
property coefficients: ndarray[tuple[Any, ...], dtype[floating[Any]]]
property coef: dict[str, float]

Named coefficient mapping (like R’s coef() or statsmodels .params).

property fitted_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Fitted values on the response scale (mu).

property linear_predictor: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Linear predictor eta = X @ beta.

property residuals_deviance: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Deviance residuals (signed).

property residuals_pearson: ndarray[tuple[Any, ...], dtype[floating[Any]]]

(y - mu) / sqrt(V(mu)).

Type:

Pearson residuals

property residuals_working: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Working residuals from the final IRLS iteration.

property residuals_response: ndarray[tuple[Any, ...], dtype[floating[Any]]]

y - mu.

Type:

Response residuals

property deviance: float

Residual deviance.

property null_deviance: float

Null deviance (intercept-only model).

property aic: float

Akaike Information Criterion.

property bic: float

Bayesian Information Criterion.

Derives −2·logL from the AIC using the same parameter count the AIC was built with (ic_param_count — the coefficients in rank plus any ML-estimated dispersion the AIC penalized: σ² for Gaussian, the shape for Gamma, θ for an auto-estimated negative binomial), then re-penalizes with log(n). Falls back to rank only for the legacy case where ic_param_count was not recorded.

property dispersion: float

Dispersion parameter.

Fixed at 1.0 for Binomial and Poisson. Estimated from data for Gaussian.

property rank: int
property df_residual: int
property df_null: int
property converged: bool
property n_iter: int
property family_name: str
property standard_errors: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Standard errors of coefficients.

Computed as sqrt(dispersion * diag((X’WX)^{-1})) where W are the final IRLS weights. Uses QR from the final iteration when available (CPU). GPU path uses direct (X’WX)^{-1} inversion, which is mathematically equivalent.

property z_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Wald test statistics (z for fixed dispersion, t for estimated).

property p_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Two-sided p-values.

Uses z-distribution for Binomial/Poisson (fixed dispersion=1). Uses t-distribution for Gaussian (estimated dispersion).

property conf_level: float

Confidence level for conf_int (default 0.95).

property conf_int: ndarray[tuple[Any, ...], dtype[floating[Any]]]

Wald confidence intervals for the coefficients, shape (p, 2).

coef ± q * se on the link scale, matching the reference distribution used for the p-values: the normal quantile for fixed-dispersion families (binomial/Poisson), the Student-t quantile at df_residual for estimated-dispersion families (Gaussian/Gamma). exp(conf_int) gives odds/rate-ratio intervals for log-link families.

property info: dict[str, Any]
property timing: dict[str, float] | None
property backend_name: str
property warnings: tuple[str, ...]
summary()[source]

Generate R-style GLM summary output.

Return type:

str

class pystatistics.regression.GLMParams(coefficients, fitted_values, linear_predictor, residuals_working, residuals_deviance, residuals_pearson, residuals_response, deviance, null_deviance, aic, dispersion, rank, df_residual, df_null, n_iter, converged, family_name, link_name, ic_param_count=None)[source]

Bases: object

Parameter payload for GLM (IRLS).

This is the immutable data computed by GLM backends.

Parameters:
coefficients: ndarray[tuple[Any, ...], dtype[floating[Any]]]
fitted_values: ndarray[tuple[Any, ...], dtype[floating[Any]]]
linear_predictor: ndarray[tuple[Any, ...], dtype[floating[Any]]]
residuals_working: ndarray[tuple[Any, ...], dtype[floating[Any]]]
residuals_deviance: ndarray[tuple[Any, ...], dtype[floating[Any]]]
residuals_pearson: ndarray[tuple[Any, ...], dtype[floating[Any]]]
residuals_response: ndarray[tuple[Any, ...], dtype[floating[Any]]]
deviance: float
null_deviance: float
aic: float
dispersion: float
rank: int
df_residual: int
df_null: int
n_iter: int
converged: bool
family_name: str
ic_param_count: int | None = None
class pystatistics.regression.Family(link=None)[source]

Bases: ABC

GLM family specification.

Defines the relationship between the mean and variance of the response distribution, along with a link function.

Parameters:

link (str | Link | None)

abstract property name: str
abstractmethod variance(mu)[source]

Variance function V(μ).

Parameters:

mu (ndarray[tuple[Any, ...], dtype[_ScalarT]])

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

abstractmethod deviance(y, mu, wt)[source]

Compute total deviance: 2 * Σ wt_i * d(y_i, μ_i).

The deviance is twice the difference between the saturated log-likelihood and the model log-likelihood.

Parameters:
Return type:

float

abstractmethod initialize(y, weights=None)[source]

Initialize μ from y for IRLS starting values.

Must return values in the valid range for the link function. weights are the per-observation prior weights (None ⇒ unit weights); only families whose R mustart depends on the prior weights (Binomial) consult them.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

property dispersion_is_fixed: bool

Whether the dispersion parameter is known a priori.

True for Binomial (φ=1) and Poisson (φ=1). False for Gaussian (φ=σ² estimated from data).

property n_ic_dispersion_params: int

Number of ML-estimated dispersion/shape parameters the information criteria penalize as free parameters, beyond the regression coefficients counted in rank.

R counts the dispersion of Gaussian (σ²) and Gamma (the shape) GLMs as a free parameter in both AIC and BIC — its logLik reports df = rank + 1. The fixed-dispersion families (Binomial, Poisson) and the fixed-θ negative binomial do not, so this returns 0 by default. aic() of the affected families adds the +2 for this parameter; recording the count here lets GLMSolution.bic re-penalize it with log(n) instead of leaving it at the AIC constant.

abstractmethod log_likelihood(y, mu, wt, dispersion)[source]

Compute the log-likelihood for AIC.

This must match R’s family$aic() / (-2) for consistency.

Parameters:
Return type:

float

aic(y, mu, wt, rank, dispersion)[source]

Compute AIC = -2 * loglik + 2 * rank.

For families with estimated dispersion (Gaussian), R computes AIC differently. Subclasses may override.

Parameters:
Return type:

float

class pystatistics.regression.Gaussian(link=None)[source]

Bases: Family

Gaussian (Normal) family. Default link: identity.

V(μ) = 1 Deviance = Σ wt_i * (y_i - μ_i)² (= RSS for identity link)

Parameters:

link (str | Link | None)

property name: str
variance(mu)[source]

Variance function V(μ).

Parameters:

mu (ndarray[tuple[Any, ...], dtype[_ScalarT]])

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

initialize(y, weights=None)[source]

Initialize μ from y for IRLS starting values.

Must return values in the valid range for the link function. weights are the per-observation prior weights (None ⇒ unit weights); only families whose R mustart depends on the prior weights (Binomial) consult them.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

deviance(y, mu, wt)[source]

Compute total deviance: 2 * Σ wt_i * d(y_i, μ_i).

The deviance is twice the difference between the saturated log-likelihood and the model log-likelihood.

Parameters:
Return type:

float

log_likelihood(y, mu, wt, dispersion)[source]

Compute the log-likelihood for AIC.

This must match R’s family$aic() / (-2) for consistency.

Parameters:
Return type:

float

aic(y, mu, wt, rank, dispersion)[source]

Compute AIC matching R’s gaussian family.

R’s gaussian()$aic uses MLE dispersion (dev/n, not dev/df) and adds +2 for the dispersion parameter. The formula is:

AIC = -2 * loglik(sigma_mle) + 2 + 2 * rank

where sigma_mle = sqrt(deviance / n), and the +2 comes from the gaussian family counting the dispersion as an extra parameter.

Parameters:
Return type:

float

property dispersion_is_fixed: bool

Whether the dispersion parameter is known a priori.

True for Binomial (φ=1) and Poisson (φ=1). False for Gaussian (φ=σ² estimated from data).

property n_ic_dispersion_params: int

Number of ML-estimated dispersion/shape parameters the information criteria penalize as free parameters, beyond the regression coefficients counted in rank.

R counts the dispersion of Gaussian (σ²) and Gamma (the shape) GLMs as a free parameter in both AIC and BIC — its logLik reports df = rank + 1. The fixed-dispersion families (Binomial, Poisson) and the fixed-θ negative binomial do not, so this returns 0 by default. aic() of the affected families adds the +2 for this parameter; recording the count here lets GLMSolution.bic re-penalize it with log(n) instead of leaving it at the AIC constant.

class pystatistics.regression.Binomial(link=None)[source]

Bases: Family

Binomial family. Default link: logit.

V(μ) = μ(1-μ) Deviance = 2 * Σ wt_i * [y_i log(y_i/μ_i) + (1-y_i) log((1-y_i)/(1-μ_i))]

Parameters:

link (str | Link | None)

property name: str
variance(mu)[source]

Variance function V(μ).

Parameters:

mu (ndarray[tuple[Any, ...], dtype[_ScalarT]])

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

initialize(y, weights=None)[source]

Initialize μ from y for IRLS starting values.

Must return values in the valid range for the link function. weights are the per-observation prior weights (None ⇒ unit weights); only families whose R mustart depends on the prior weights (Binomial) consult them.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

deviance(y, mu, wt)[source]

Compute total deviance: 2 * Σ wt_i * d(y_i, μ_i).

The deviance is twice the difference between the saturated log-likelihood and the model log-likelihood.

Parameters:
Return type:

float

log_likelihood(y, mu, wt, dispersion)[source]

Compute the log-likelihood for AIC.

This must match R’s family$aic() / (-2) for consistency.

Parameters:
Return type:

float

property dispersion_is_fixed: bool

Whether the dispersion parameter is known a priori.

True for Binomial (φ=1) and Poisson (φ=1). False for Gaussian (φ=σ² estimated from data).

class pystatistics.regression.Poisson(link=None)[source]

Bases: Family

Poisson family. Default link: log.

V(μ) = μ Deviance = 2 * Σ wt_i * [y_i log(y_i/μ_i) - (y_i - μ_i)]

Parameters:

link (str | Link | None)

property name: str
variance(mu)[source]

Variance function V(μ).

Parameters:

mu (ndarray[tuple[Any, ...], dtype[_ScalarT]])

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

initialize(y, weights=None)[source]

Initialize μ from y for IRLS starting values.

Must return values in the valid range for the link function. weights are the per-observation prior weights (None ⇒ unit weights); only families whose R mustart depends on the prior weights (Binomial) consult them.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

deviance(y, mu, wt)[source]

Compute total deviance: 2 * Σ wt_i * d(y_i, μ_i).

The deviance is twice the difference between the saturated log-likelihood and the model log-likelihood.

Parameters:
Return type:

float

log_likelihood(y, mu, wt, dispersion)[source]

Compute the log-likelihood for AIC.

This must match R’s family$aic() / (-2) for consistency.

Parameters:
Return type:

float

property dispersion_is_fixed: bool

Whether the dispersion parameter is known a priori.

True for Binomial (φ=1) and Poisson (φ=1). False for Gaussian (φ=σ² estimated from data).

class pystatistics.regression.GammaFamily(link=None)[source]

Bases: Family

Gamma family. Default link: inverse.

V(μ) = μ² Deviance = 2 * Σ wt_i * [(y_i - μ_i)/μ_i - log(y_i/μ_i)]

Used for positive continuous data with variance proportional to mean². Typical applications: cost data, survival times, insurance claims.

References

R: stats::Gamma() McCullagh & Nelder (1989), Ch. 8

Parameters:

link (str | Link | None)

property name: str
variance(mu)[source]

Variance function V(μ).

Parameters:

mu (ndarray[tuple[Any, ...], dtype[_ScalarT]])

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

initialize(y, weights=None)[source]

Initialize μ from y for IRLS starting values.

Must return values in the valid range for the link function. weights are the per-observation prior weights (None ⇒ unit weights); only families whose R mustart depends on the prior weights (Binomial) consult them.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

deviance(y, mu, wt)[source]

Compute total deviance: 2 * Σ wt_i * d(y_i, μ_i).

The deviance is twice the difference between the saturated log-likelihood and the model log-likelihood.

Parameters:
Return type:

float

log_likelihood(y, mu, wt, dispersion)[source]

Compute the log-likelihood for AIC.

This must match R’s family$aic() / (-2) for consistency.

Parameters:
Return type:

float

aic(y, mu, wt, rank, dispersion)[source]

Compute AIC matching R’s Gamma()$aic.

R’s Gamma family evaluates the AIC log-likelihood at a dispersion of dev / sum(wt) — the MLE of the dispersion under the Gamma distribution — NOT the moment estimate dev / df_residual that R reports in summary.glm and that PyStatistics stores in GLMParams.dispersion for standard errors. The two diverge whenever rank > 0, so the dispersion argument (which the solver derives from df_residual) must be ignored here and the AIC-specific dispersion recomputed internally.

R also counts the estimated dispersion/shape as a free parameter, adding +2 on top of 2 * rank. Concretely R computes Gamma()$aic as -2 * sum(wt * dgamma(y, 1/disp, scale=mu*disp, log=TRUE)) + 2 and glm.fit adds 2 * rank, giving the formula below. log_likelihood evaluated at disp is algebraically the weighted dgamma sum.

Parameters:
Return type:

float

property dispersion_is_fixed: bool

Whether the dispersion parameter is known a priori.

True for Binomial (φ=1) and Poisson (φ=1). False for Gaussian (φ=σ² estimated from data).

property n_ic_dispersion_params: int

Number of ML-estimated dispersion/shape parameters the information criteria penalize as free parameters, beyond the regression coefficients counted in rank.

R counts the dispersion of Gaussian (σ²) and Gamma (the shape) GLMs as a free parameter in both AIC and BIC — its logLik reports df = rank + 1. The fixed-dispersion families (Binomial, Poisson) and the fixed-θ negative binomial do not, so this returns 0 by default. aic() of the affected families adds the +2 for this parameter; recording the count here lets GLMSolution.bic re-penalize it with log(n) instead of leaving it at the AIC constant.

class pystatistics.regression.NegativeBinomial(theta=None, link=None)[source]

Bases: Family

Negative binomial family. Default link: log.

V(μ) = μ + μ²/θ (where θ is the dispersion parameter)

For fixed θ, this is a standard GLM with known variance function. When θ is unknown, it must be estimated via profile likelihood (see regression._nb_theta).

Parameters:

  • theta (float | None) – The dispersion parameter (> 0). Larger θ means less overdispersion; θ → ∞ recovers Poisson. If None, theta must be estimated externally (e.g., via fit(family=’negative.binomial’)).

  • link (str | Link | None) – Link function (default: log).

References

R: MASS::negative.binomial(), MASS::glm.nb() Venables & Ripley (2002), Modern Applied Statistics with S, Ch. 7.4

property name: str
variance(mu)[source]

Variance function V(μ).

Parameters:

mu (ndarray[tuple[Any, ...], dtype[_ScalarT]])

Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

initialize(y, weights=None)[source]

Initialize μ from y for IRLS starting values.

Must return values in the valid range for the link function. weights are the per-observation prior weights (None ⇒ unit weights); only families whose R mustart depends on the prior weights (Binomial) consult them.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[_ScalarT]]

deviance(y, mu, wt)[source]

Compute total deviance: 2 * Σ wt_i * d(y_i, μ_i).

The deviance is twice the difference between the saturated log-likelihood and the model log-likelihood.

Parameters:
Return type:

float

log_likelihood(y, mu, wt, dispersion)[source]

Compute the log-likelihood for AIC.

This must match R’s family$aic() / (-2) for consistency.

Parameters:
Return type:

float

property dispersion_is_fixed: bool

Whether the dispersion parameter is known a priori.

True for Binomial (φ=1) and Poisson (φ=1). False for Gaussian (φ=σ² estimated from data).