From 0a1f49c92fd0f3f5240118716abe062950360046 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Simon=20M=C3=BCller?= Date: Thu, 28 May 2026 16:43:34 +0200 Subject: [PATCH] docs: update for v0.5.0 additions (#34) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bring the mkdocs site and CHANGELOG in line with everything that landed across PRs #20–#33. docs/api/regression/diagnostics.md — full sweep adding ~20 new entries grouped by category: - Multicollinearity: vif, generalized_vif, high_vif_predictors - OLS residual battery: standardized / studentized / externally studentized residuals, residual_outliers - Influence / leverage: leverage, cooks_distance, dffits, influential_cooks, influential_dffits, high_leverage_points - GLM residuals: pearson / deviance / working for logistic and poisson - Goodness of fit: pearson_chi_squared_logistic / _poisson Other docs: - regression/linear.md: huber, pls - regression/glm.md: logistic_regression (sklearn-style); lambda_ kwarg note on all GLM expressions - regression/alm.md: 25 distributions; loss/link/role_trim/ extra_parameter kwargs - regression/summary-predict.md: quantile_summary, quantile_predict, isotonic_predict, lm_dynamic_predict - classes/linear.md: Huber, PLS classes - classes/glm.md: LogisticRegression; lambda_ on existing classes - outputs.md: 10 new struct schemas copied from src/expressions/regression.rs - index.md, getting-started.md: brief "Use from Rust" + diagnostics quick tour Drive-by fixes: - docs/API_REFERENCE.md and docs/api/README.md had broken ../README.md links and a dead TOC anchor that blocked `mkdocs build --strict`. CHANGELOG.md: new 0.5.0 section listing every Added / Changed item above with PR references. Also corrected the 0.4.0 entry's release date from 2025-01-17 to 2026-01-17 to match the published tag. `mkdocs build --strict` passes (only the intentional API_REFERENCE.md "not in nav" info remains). Co-Authored-By: Claude Opus 4.7 (1M context) --- CHANGELOG.md | 62 +++- docs/API_REFERENCE.md | 5 +- docs/api/README.md | 2 +- docs/api/classes/glm.md | 36 +++ docs/api/classes/linear.md | 63 ++++ docs/api/outputs.md | 135 +++++++++ docs/api/regression/alm.md | 56 +++- docs/api/regression/diagnostics.md | 400 ++++++++++++++++++++++++- docs/api/regression/glm.md | 38 ++- docs/api/regression/linear.md | 49 +++ docs/api/regression/summary-predict.md | 29 ++ docs/getting-started.md | 22 ++ docs/index.md | 12 + 13 files changed, 889 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f7aef6..b64ea20 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,66 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.4.0] - 2025-01-17 +## [0.5.0] - 2026-05-28 + +### Added + +- **Hybrid crate (Rust + Python)** — `polars-statistics` now builds as both a + `cdylib` (Python plugin) and an `rlib` (Rust dependency). Every Polars + expression has a public `_fit` Rust entry point in + `polars_statistics::expressions`. New `python` Cargo feature gates the + pyo3 / numpy / pymodels surface; downstream Rust crates use + `default-features = false`. (#13) +- **New regression model wrappers**: + - `Huber` — M-estimator robust to outliers; class + `huber()` expression. (#14) + - `LogisticRegression` — sklearn-style API with `predict_proba`, + `decision_function`, `score`, `penalty="l2"`, `C` kwarg. Distinct from + the existing `Logistic` wrapper. (#14) + - `PLS` — Partial Least Squares with `transform()` for the latent space; + class + `pls()` expression. (#19) +- **Penalized IRLS `lambda_` kwarg** on the GLM PyClass models — `PyLogistic`, + `PyPoisson`, `PyNegativeBinomial`, `PyTweedie`, `PyProbit`, `PyCloglog` + (the expression layer already supported it). (#15) +- **ALM expression parity with `PyALM`** — all 25 distributions reachable + from `ps.alm(...)`, plus `loss`, `link`, `role_trim`, and + `extra_parameter` kwargs. (#16) +- **Summary / predict completeness** — added the matching expressions for + families that previously only had a base fit (#18): + - `quantile_summary`, `quantile_predict` + - `isotonic_predict` + - `lm_dynamic_predict` +- **Diagnostics toolkit** (#17 + #27): + - Multicollinearity: `vif`, `generalized_vif`, `high_vif_predictors` + - OLS residual battery: `standardized_residuals`, + `studentized_residuals`, `externally_studentized_residuals`, + `residual_outliers` + - GLM residuals (logistic + Poisson): `*_pearson_residuals`, + `*_deviance_residuals`, `*_working_residuals` for each family + - Influence / leverage: `leverage`, `cooks_distance`, `dffits`, + `influential_cooks`, `influential_dffits`, `high_leverage_points` + - Goodness of fit: `pearson_chi_squared_logistic`, + `pearson_chi_squared_poisson` +- **Documentation**: + - Sweep of the README and mkdocs site to cover every v0.5.0 addition. + - New "Use from Rust" section + `examples/rust_wls.rs` walking through + the rlib path. + +### Changed + +- Updated `anofox-regression` dependency to v0.5.4 (introduces + `HuberRegressor` and the sklearn-style `LogisticRegression`) +- Updated `anofox-statistics` dependency to v0.4.1 +- Python `__version__` caught up from 0.3.0 → 0.5.0 (was lagging two + minor versions behind the wheel metadata) + +### Backwards compatibility + +Additive. `PyLogistic` is unchanged — the new sklearn-style +`LogisticRegression` is a separate class. The `alm()` expression's input +contract grew but all existing keyword-only callers continue to work +because the new kwargs default to `None` / `"likelihood"`. + +## [0.4.0] - 2026-01-17 ### Added @@ -197,6 +256,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - SIMD-optimized linear algebra via faer - Automatic parallelization for group operations +[0.5.0]: https://github.com/DataZooDE/polars-statistics/compare/v0.4.0...v0.5.0 [0.4.0]: https://github.com/DataZooDE/polars-statistics/compare/v0.3.0...v0.4.0 [0.3.0]: https://github.com/DataZooDE/polars-statistics/compare/v0.2.0...v0.3.0 [0.2.0]: https://github.com/DataZooDE/polars-statistics/compare/v0.1.0...v0.2.0 diff --git a/docs/API_REFERENCE.md b/docs/API_REFERENCE.md index 249acd8..8e5cea5 100644 --- a/docs/API_REFERENCE.md +++ b/docs/API_REFERENCE.md @@ -6,7 +6,7 @@ > **Version:** 0.4.0 > **Last Updated:** 2026-01-17 -Complete API reference for polars-statistics. For quick start examples, see the [README](../README.md). +Complete API reference for polars-statistics. For quick start examples, see the [README](https://github.com/DataZooDE/polars-statistics#readme). --- @@ -24,7 +24,6 @@ Complete API reference for polars-statistics. For quick start examples, see the - [TOST Equivalence Tests](#tost-equivalence-tests) - [Regression Models](#regression-models) - [Linear Models](#linear-models) - - [Robust Regression](#robust-regression) - [Diagnostics](#diagnostics) - [GLM Models](#glm-models) - [Augmented Linear Model (ALM)](#augmented-linear-model-alm) @@ -1978,7 +1977,7 @@ Struct { ## See Also -- [README](../README.md) - Quick start guide +- [README](https://github.com/DataZooDE/polars-statistics#readme) - Quick start guide - [Polars Documentation](https://docs.pola.rs/) - [faer](https://github.com/sarah-ek/faer-rs) - Linear algebra backend diff --git a/docs/api/README.md b/docs/api/README.md index 19b7c48..2a4945c 100644 --- a/docs/api/README.md +++ b/docs/api/README.md @@ -68,5 +68,5 @@ result.with_columns( ## See Also -- [README](../../README.md) - Quick start guide +- [README](https://github.com/DataZooDE/polars-statistics#readme) - Quick start guide - [Polars Documentation](https://docs.pola.rs/) diff --git a/docs/api/classes/glm.md b/docs/api/classes/glm.md index c3931b4..695f91a 100644 --- a/docs/api/classes/glm.md +++ b/docs/api/classes/glm.md @@ -2,6 +2,8 @@ Direct GLM model access outside of Polars expressions. +All GLM classes (`Logistic`, `Poisson`, `NegativeBinomial`, `Tweedie`, `Probit`, `Cloglog`) accept a `lambda_=0.0` kwarg for L2 (ridge) penalty applied inside the IRLS update. The sklearn-style [`LogisticRegression`](#logisticregression) class exposes the same penalty via `C = 1 / lambda_` and an explicit `penalty` choice. + ## Common Interface ```python @@ -44,6 +46,39 @@ probabilities = model.predict_proba(X_new) # Probability estimates --- +## LogisticRegression + +Sklearn-style logistic regression. Distinct from [`Logistic`](#logistic): uses inverse-strength regularization `C = 1 / lambda_` and an explicit `penalty` choice. + +```python +from polars_statistics import LogisticRegression + +model = LogisticRegression( + penalty: str = "l2", # "l2" or "none" + C: float = 1.0, # Inverse of regularization strength + threshold: float = 0.5, + with_intercept: bool = True, + max_iter: int = 100, + tol: float = 1e-8, + compute_inference: bool = True, + confidence_level: float = 0.95, +) +model.fit(X, y) # y: binary (0/1) + +# Methods +classes = model.predict(X_new) # 0/1 predictions +probs = model.predict_proba(X_new) # Probability estimates +scores = model.decision_function(X_new) # Linear scores (log-odds) +accuracy = model.score(X_new, y_new) # Mean accuracy + +# Properties +model.coefficients # np.ndarray +model.intercept # float or None +model.n_iter # int — IRLS iterations until convergence +``` + +--- + ## Poisson Poisson regression for count data. @@ -135,6 +170,7 @@ model.fit(X, y) # y: binary (0/1) | Class | Parameters | |-------|------------| | `Logistic` | `lambda_`, `with_intercept` | +| `LogisticRegression` | `penalty`, `C`, `threshold`, `with_intercept`, `max_iter`, `tol`, `compute_inference`, `confidence_level` | | `Poisson` | `lambda_`, `with_intercept` | | `NegativeBinomial` | `theta`, `estimate_theta`, `lambda_`, `with_intercept` | | `Tweedie` | `var_power`, `lambda_`, `with_intercept` | diff --git a/docs/api/classes/linear.md b/docs/api/classes/linear.md index bda9ce1..1628ee9 100644 --- a/docs/api/classes/linear.md +++ b/docs/api/classes/linear.md @@ -129,6 +129,67 @@ model.fit(X, y) --- +## Huber + +Huber M-estimator (robust regression). + +```python +from polars_statistics import Huber + +model = Huber( + epsilon: float = 1.35, + alpha: float = 0.0001, + with_intercept: bool = True, + max_iter: int = 100, + tol: float = 1e-5, +) +model.fit(X, y) + +# Properties +model.coefficients # np.ndarray +model.intercept # float or None +model.scale # float — robust scale estimate +model.epsilon # float +model.outliers # np.ndarray[bool] — down-weighted observations +model.n_outliers # int +model.r_squared # float +model.mse # float +model.rmse # float +model.n_observations # int +``` + +--- + +## PLS + +Partial Least Squares regression. + +```python +from polars_statistics import PLS + +model = PLS( + n_components: int = 2, + with_intercept: bool = True, + tol: float = 1e-6, + scale: bool = True, +) +model.fit(X, y) + +# Methods +predictions = model.predict(X_new) +latent = model.transform(X_new) # Project into latent space + +# Properties +model.coefficients # np.ndarray +model.intercept # float or None +model.n_components # int +model.explained_variance_ratio # np.ndarray (one entry per component) +model.r_squared # float +model.n_observations # int +``` + +--- + ## Quantile Quantile regression. @@ -179,6 +240,8 @@ model.fitted_values # np.ndarray | `WLS` | `with_intercept`, `compute_inference` | | `RLS` | `forgetting_factor`, `with_intercept` | | `BLS` | `lower_bound`, `upper_bound`, `with_intercept` | +| `Huber` | `epsilon`, `alpha`, `with_intercept`, `max_iter`, `tol` | +| `PLS` | `n_components`, `with_intercept`, `tol`, `scale` | | `Quantile` | `tau`, `with_intercept` | | `Isotonic` | `increasing` | diff --git a/docs/api/outputs.md b/docs/api/outputs.md index 9255abb..3812b32 100644 --- a/docs/api/outputs.md +++ b/docs/api/outputs.md @@ -88,6 +88,141 @@ Struct { --- +## VIF Output + +Returned by `vif`. + +``` +Struct { + terms: List[String], # "x1", "x2", ... + vif: List[Float64], + n_observations: UInt32, +} +``` + +--- + +## VIF Mask Output + +Returned by `high_vif_predictors`. + +``` +Struct { + is_high: List[Boolean], # one per input feature column + n_high: UInt32, + n_features: UInt32, +} +``` + +--- + +## GVIF Output + +Returned by `generalized_vif`. + +``` +Struct { + gvif: List[Float64], # one value per group + n_groups: UInt32, +} +``` + +--- + +## Leverage Output + +Returned by `leverage`. + +``` +Struct { + leverage: List[Float64], # one h_ii per input row + n_observations: UInt32, +} +``` + +--- + +## Cook's Distance Output + +Returned by `cooks_distance`. + +``` +Struct { + cooks_d: List[Float64], # one D_i per input row + n_observations: UInt32, +} +``` + +--- + +## DFFITS Output + +Returned by `dffits`. + +``` +Struct { + dffits: List[Float64], # one value per input row + n_observations: UInt32, +} +``` + +--- + +## Influence Mask Output + +Returned by `influential_cooks`, `influential_dffits`, `high_leverage_points`. + +``` +Struct { + is_influential: List[Boolean], + n_influential: UInt32, + n_observations: UInt32, +} +``` + +--- + +## Residual Diagnostics Output + +Returned by `standardized_residuals`, `studentized_residuals`, `externally_studentized_residuals`, and all GLM residual functions (`logistic_*_residuals`, `poisson_*_residuals`). + +``` +Struct { + residuals: List[Float64], # one residual per input row + n_observations: UInt32, +} +``` + +--- + +## Outlier Mask Output + +Returned by `residual_outliers`. + +``` +Struct { + is_outlier: List[Boolean], + n_outliers: UInt32, + n_observations: UInt32, +} +``` + +--- + +## Chi-Squared Output + +Returned by `pearson_chi_squared_logistic` and `pearson_chi_squared_poisson`. + +``` +Struct { + chi_squared: Float64, + df_resid: UInt32, + n_observations: UInt32, +} +``` + +--- + ## GLM Output Returned by `logistic`, `poisson`, `negative_binomial`, `tweedie`, `probit`, `cloglog`. diff --git a/docs/api/regression/alm.md b/docs/api/regression/alm.md index 3ffef66..9ab23b5 100644 --- a/docs/api/regression/alm.md +++ b/docs/api/regression/alm.md @@ -1,6 +1,6 @@ # Augmented Linear Model (ALM) -Flexible regression supporting 24+ distributions with automatic link function selection. +Flexible regression supporting 25 distributions with configurable link functions and loss criteria. ## `alm` @@ -9,6 +9,10 @@ ps.alm( y: Union[pl.Expr, str], *x: Union[pl.Expr, str], distribution: str = "normal", + link: str | None = None, # None = canonical link for the distribution + loss: str = "likelihood", # "likelihood" | "mse" | "mae" | "ham" | "role" + role_trim: float | None = None, # Trim fraction for the "role" loss + extra_parameter: float | None = None, with_intercept: bool = True, ) -> pl.Expr ``` @@ -17,22 +21,59 @@ ps.alm( **Example:** ```python -# Fit with Laplace distribution for robust estimation +# Robust regression via Laplace likelihood df.group_by("group").agg( ps.alm("y", "x1", "x2", distribution="laplace").alias("model") ) + +# Gamma regression with log link +df.group_by("group").agg( + ps.alm("y", "x1", distribution="gamma", link="log").alias("model") +) ``` --- -## Supported Distributions +## Supported Distributions (25) | Category | Distributions | |----------|---------------| -| Continuous | `normal`, `laplace`, `student_t`, `logistic` | -| Positive | `lognormal`, `loglaplace`, `gamma`, `inverse_gaussian`, `exponential` | -| Bounded (0,1) | `beta` | +| Continuous | `normal`, `laplace`, `student_t`, `logistic`, `asymmetric_laplace`, `generalised_normal`, `s` | +| Positive | `lognormal`, `loglaplace`, `logs`, `loggeneralisednormal`, `gamma`, `inverse_gaussian`, `exponential`, `folded_normal`, `rectified_normal` | +| Bounded (0,1) | `beta`, `logit_normal` | | Count | `poisson`, `negative_binomial`, `binomial`, `geometric` | +| Ordinal | `cumulative_logistic`, `cumulative_normal` | +| Transformed | `boxcox_normal` | + +--- + +## Link Functions + +When `link=None` ALM picks the canonical link for the distribution. Explicit options: + +| `link` | Inverse | Typical use | +|--------|---------|-------------| +| `"identity"` | `η` | Gaussian | +| `"log"` | `exp(η)` | Positive / count | +| `"logit"` | `1 / (1 + exp(-η))` | Binary, beta | +| `"probit"` | `Φ(η)` | Binary | +| `"inverse"` | `1 / η` | Gamma | +| `"sqrt"` | `η²` | Poisson | +| `"cloglog"` | `1 - exp(-exp(η))` | Binary (asymmetric) | + +--- + +## Loss Functions + +| `loss` | Description | +|--------|-------------| +| `"likelihood"` | Maximum likelihood (default) | +| `"mse"` | Mean squared error | +| `"mae"` | Mean absolute error | +| `"ham"` | Half-absolute-moment | +| `"role"` | Robust loss with trimming via `role_trim` (default 0.05) | + +`extra_parameter` is required by distributions that take an auxiliary parameter (e.g. degrees of freedom for `student_t`, shape for `generalised_normal`). --- @@ -45,9 +86,10 @@ df.group_by("group").agg( | Heavy tails | `student_t` | | Positive continuous | `lognormal`, `gamma` | | Right-skewed positive | `gamma`, `inverse_gaussian` | -| Proportions/rates | `beta` | +| Proportions/rates | `beta`, `logit_normal` | | Count data | `poisson`, `negative_binomial` | | Overdispersed counts | `negative_binomial` | +| Ordinal outcomes | `cumulative_logistic`, `cumulative_normal` | --- diff --git a/docs/api/regression/diagnostics.md b/docs/api/regression/diagnostics.md index acc8a52..0901204 100644 --- a/docs/api/regression/diagnostics.md +++ b/docs/api/regression/diagnostics.md @@ -1,6 +1,6 @@ # Regression Diagnostics -Tools for detecting multicollinearity, quasi-separation, and other model issues. +Tools for detecting multicollinearity, quasi-separation, influence and outliers, plus residual batteries for OLS and GLM fits. ## `condition_number` @@ -86,6 +86,393 @@ df.select(ps.check_count_sparsity("count", "x1", "x2")) --- +## Multicollinearity + +### `vif` + +Variance Inflation Factor per predictor: `VIF_j = 1 / (1 - R²_j)`. + +```python +ps.vif(*x: Union[pl.Expr, str]) -> pl.Expr +``` + +**Returns:** See [VIF Output](../outputs.md#vif-output) + +**Example:** +```python +df.select(ps.vif("x1", "x2", "x3").alias("vif")) +``` + +VIF > 5 (moderate) / > 10 (severe) typically flags problematic collinearity. The intercept is not included. + +--- + +### `generalized_vif` + +Generalized VIF (GVIF) for grouped predictors such as one-hot dummies. + +```python +ps.generalized_vif( + *x: Union[pl.Expr, str], + group_sizes: list[int], # e.g. [1, 1, 3] for two scalars + a 3-level dummy +) -> pl.Expr +``` + +**Returns:** See [GVIF Output](../outputs.md#gvif-output) + +**Example:** +```python +# x1, x2 are scalars; cat_a/cat_b/cat_c are dummies for one categorical +df.select( + ps.generalized_vif("x1", "x2", "cat_a", "cat_b", "cat_c", + group_sizes=[1, 1, 3]).alias("gvif") +) +``` + +For single-column groups GVIF coincides with `vif`. + +--- + +### `high_vif_predictors` + +Boolean mask of features whose VIF exceeds `threshold` (default 10). + +```python +ps.high_vif_predictors( + *x: Union[pl.Expr, str], + threshold: float = 10.0, +) -> pl.Expr +``` + +**Returns:** See [VIF Mask Output](../outputs.md#vif-mask-output) + +**Example:** +```python +df.select(ps.high_vif_predictors("x1", "x2", "x3", threshold=5.0)) +``` + +--- + +## OLS Residual Battery + +Each function fits OLS internally and returns one residual per row. + +### `standardized_residuals` + +`r_i / sqrt(MSE)`. + +```python +ps.standardized_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +**Example:** +```python +df.select(ps.standardized_residuals("y", "x1", "x2").alias("r_std")) +``` + +--- + +### `studentized_residuals` + +Internally studentized: `r_i / sqrt(MSE * (1 - h_ii))`. + +```python +ps.studentized_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +### `externally_studentized_residuals` + +Leave-one-out studentized residuals; t-distributed with `n - p - 1` df under the null. + +```python +ps.externally_studentized_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +### `residual_outliers` + +Boolean outlier mask based on `|studentized residual| > threshold` (default 2.0). + +```python +ps.residual_outliers( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + threshold: float = 2.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Outlier Mask Output](../outputs.md#outlier-mask-output) + +**Example:** +```python +df.select(ps.residual_outliers("y", "x1", "x2", threshold=2.5)) +``` + +--- + +## Influence / Leverage + +### `leverage` + +Hat-matrix diagonal `h_ii`. Sums to `p`; flags `h_ii > 2p/n` are conventional. + +```python +ps.leverage( + *x: Union[pl.Expr, str], + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Leverage Output](../outputs.md#leverage-output) + +**Example:** +```python +df.select(ps.leverage("x1", "x2").alias("h")) +``` + +--- + +### `cooks_distance` + +`D_i = (e_i² / (p · MSE)) · (h_ii / (1 - h_ii)²)`. Common cutoffs: `D_i > 4/n` or `D_i > 1`. + +```python +ps.cooks_distance( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Cook's Distance Output](../outputs.md#cooks-distance-output) + +**Example:** +```python +df.select(ps.cooks_distance("y", "x1", "x2").alias("d")) +``` + +--- + +### `dffits` + +Scaled change in fitted value when observation `i` is dropped. Cutoff: `|DFFITS_i| > 2 * sqrt(p / n)`. + +```python +ps.dffits( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [DFFITS Output](../outputs.md#dffits-output) + +--- + +### `influential_cooks` + +Boolean mask of observations with Cook's distance above `threshold` (default `4 / n`). + +```python +ps.influential_cooks( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + threshold: float | None = None, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Influence Mask Output](../outputs.md#influence-mask-output) + +--- + +### `influential_dffits` + +Boolean mask of observations with `|DFFITS|` above `threshold` (default `2 * sqrt(p / n)`). + +```python +ps.influential_dffits( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + threshold: float | None = None, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Influence Mask Output](../outputs.md#influence-mask-output) + +--- + +### `high_leverage_points` + +Boolean mask of high-leverage observations. Takes only feature columns; default threshold `2 * p / n`. + +```python +ps.high_leverage_points( + *x: Union[pl.Expr, str], + threshold: float | None = None, + add_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Influence Mask Output](../outputs.md#influence-mask-output) + +--- + +## GLM Residuals + +Each function fits the appropriate GLM internally and returns one residual per row. All accept `lambda_=0.0` for ridge-penalized IRLS. + +### `logistic_pearson_residuals` + +```python +ps.logistic_pearson_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +### `logistic_deviance_residuals` + +```python +ps.logistic_deviance_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +### `logistic_working_residuals` + +IRLS adjusted-dependent-variable residuals. + +```python +ps.logistic_working_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +### `poisson_pearson_residuals` + +```python +ps.poisson_pearson_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +### `poisson_deviance_residuals` + +```python +ps.poisson_deviance_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +### `poisson_working_residuals` + +```python +ps.poisson_working_residuals( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Residual Diagnostics Output](../outputs.md#residual-diagnostics-output) + +--- + +## Goodness of Fit + +### `pearson_chi_squared_logistic` + +`Σ pearson_residual²` from a logistic fit, plus residual degrees of freedom. For a well-specified model `X² / df_resid ≈ 1`. + +```python +ps.pearson_chi_squared_logistic( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Chi-Squared Output](../outputs.md#chi-squared-output) + +--- + +### `pearson_chi_squared_poisson` + +`Σ pearson_residual²` from a Poisson fit, plus residual degrees of freedom. + +```python +ps.pearson_chi_squared_poisson( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + lambda_: float = 0.0, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Chi-Squared Output](../outputs.md#chi-squared-output) + +--- + ## Diagnostic Workflow ```python @@ -94,19 +481,18 @@ import polars_statistics as ps # 1. Check multicollinearity before fitting cond = df.select(ps.condition_number("x1", "x2", "x3")) -print(cond) +vif = df.select(ps.vif("x1", "x2", "x3")) # 2. For binary outcomes, check separation sep = df.select(ps.check_binary_separation("y", "x1", "x2")) -if sep["has_separation"][0]: - print("Warning: separation detected") - # Use regularization +if sep["separation"][0]["has_separation"]: model = df.select(ps.logistic("y", "x1", "x2", lambda_=1.0)) else: model = df.select(ps.logistic("y", "x1", "x2")) -# 3. For count data, check sparsity -sparse = df.select(ps.check_count_sparsity("count", "x1", "x2")) +# 3. After fitting OLS, scan for outliers and influential points +df.select(ps.residual_outliers("y", "x1", "x2").alias("outliers")) +df.select(ps.influential_cooks("y", "x1", "x2").alias("influence")) ``` --- diff --git a/docs/api/regression/glm.md b/docs/api/regression/glm.md index 07a1327..9a7201e 100644 --- a/docs/api/regression/glm.md +++ b/docs/api/regression/glm.md @@ -1,6 +1,14 @@ # Generalized Linear Models (GLM) -GLM models for binary classification, count data, and other non-normal response distributions. All GLM models support optional Ridge regularization via the `lambda_` parameter. +GLM models for binary classification, count data, and other non-normal response distributions. + +## Penalized IRLS + +All GLM expressions (`logistic`, `poisson`, `negative_binomial`, `tweedie`, `probit`, `cloglog`) accept a `lambda_` kwarg (default `0.0`) that adds an L2 (ridge) penalty on the coefficients inside the IRLS update. Set `lambda_ > 0` to stabilize estimation under collinearity or quasi-separation. + +The sklearn-style `logistic_regression` exposes the same penalty via `C = 1 / lambda_` and an explicit `penalty` choice. + +--- ## `logistic` @@ -24,6 +32,34 @@ df.group_by("group").agg(ps.logistic("success", "x1", "x2").alias("model")) --- +## `logistic_regression` + +Sklearn-style logistic regression. Distinct from [`logistic`](#logistic): uses inverse-strength regularization `C = 1 / lambda_` and an explicit `penalty` choice. Returns the same [GLM Output](../outputs.md#glm-output) schema. + +```python +ps.logistic_regression( + y: Union[pl.Expr, str], # Binary (0/1) + *x: Union[pl.Expr, str], + penalty: str = "l2", # "l2" or "none" + C: float = 1.0, # Inverse of regularization strength + threshold: float = 0.5, + max_iter: int = 100, + tol: float = 1e-8, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [GLM Output](../outputs.md#glm-output) + +**Example:** +```python +df.group_by("group").agg( + ps.logistic_regression("success", "x1", "x2", penalty="l2", C=0.5).alias("model") +) +``` + +--- + ## `poisson` Poisson regression for count data. diff --git a/docs/api/regression/linear.md b/docs/api/regression/linear.md index 257810b..b93969b 100644 --- a/docs/api/regression/linear.md +++ b/docs/api/regression/linear.md @@ -126,6 +126,55 @@ ps.nnls( --- +## `huber` + +Huber M-estimator (robust regression). Down-weights observations with `|residual / scale| > epsilon`. + +```python +ps.huber( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + epsilon: float = 1.35, # Huber threshold (must be > 1.0) + alpha: float = 0.0001, # L2 ridge penalty in the WLS update + max_iter: int = 100, + tol: float = 1e-5, + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Linear Model Output](../outputs.md#linear-model-output) + +**Example:** +```python +df.group_by("group").agg(ps.huber("y", "x1", "x2").alias("model")) +``` + +--- + +## `pls` + +Partial Least Squares — projects `X` onto `n_components` latent variables maximizing covariance with `y`. Useful for highly collinear features. + +```python +ps.pls( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + n_components: int = 2, + tol: float = 1e-6, + scale: bool = True, # Scale X to unit variance before fitting + with_intercept: bool = True, +) -> pl.Expr +``` + +**Returns:** See [Linear Model Output](../outputs.md#linear-model-output) + +**Example:** +```python +df.group_by("group").agg(ps.pls("y", "x1", "x2", "x3", n_components=2).alias("model")) +``` + +--- + ## `quantile` Quantile regression for estimating conditional quantiles (e.g., median). diff --git a/docs/api/regression/summary-predict.md b/docs/api/regression/summary-predict.md index 2d04249..684ec12 100644 --- a/docs/api/regression/summary-predict.md +++ b/docs/api/regression/summary-predict.md @@ -22,10 +22,13 @@ ps.tweedie_summary(y, *x, var_power=1.5, with_intercept=True) -> pl.Expr ps.probit_summary(y, *x, with_intercept=True) -> pl.Expr ps.cloglog_summary(y, *x, with_intercept=True) -> pl.Expr ps.alm_summary(y, *x, distribution="normal", with_intercept=True) -> pl.Expr +ps.quantile_summary(y, *x, tau=0.5, with_intercept=True) -> pl.Expr ``` Formula variants also available: `ps.ols_formula_summary(formula, ...)`, etc. +**Note on `quantile_summary`:** Quantile regression has no analytic standard errors, so the `std_error`, `statistic`, and `p_value` columns are returned as NaN. The coefficient estimates themselves are still valid. + **Returns:** See [Summary Output](../outputs.md#summary-output) ### Example @@ -71,9 +74,35 @@ All models have prediction functions: - `logistic_predict`, `poisson_predict`, `negative_binomial_predict` - `tweedie_predict`, `probit_predict`, `cloglog_predict` - `alm_predict` +- `quantile_predict`, `isotonic_predict`, `lm_dynamic_predict` Formula variants also available: `ps.ols_formula_predict(formula, ...)`, etc. +**Models without analytic intervals.** The following predictors return `lower` / `upper` as NaN — point predictions only: + +| Function | Why no intervals | +|----------|------------------| +| `quantile_predict` | Quantile regression intervals require bootstrap (not yet wired through). | +| `isotonic_predict` | Step-function fit with no parametric uncertainty estimate. | +| `lm_dynamic_predict` | Predictions use time-averaged coefficients from the dynamic model. | + +```python +# Quantile regression predictions +df.with_columns( + ps.quantile_predict("y", "x1", "x2", tau=0.5).over("group").alias("p") +).unnest("p") + +# Isotonic step-function predictions (single feature) +df.with_columns( + ps.isotonic_predict("y", "x").alias("p") +).unnest("p") + +# Dynamic linear model predictions +df.with_columns( + ps.lm_dynamic_predict("y", "x1", "x2").over("group").alias("p") +).unnest("p") +``` + **Returns:** See [Prediction Output](../outputs.md#prediction-output) ### Parameters diff --git a/docs/getting-started.md b/docs/getting-started.md index e3bb1b5..31ef51a 100644 --- a/docs/getting-started.md +++ b/docs/getting-started.md @@ -94,6 +94,28 @@ df.group_by("group").agg( # Columns: term, estimate, std_error, statistic, p_value ``` +## Diagnostics Quick Tour + +Detect multicollinearity, outliers, and influential points before (or after) fitting: + +```python +import polars as pl +import polars_statistics as ps + +# 1. Multicollinearity — Variance Inflation Factor per predictor +df.select(ps.vif("x1", "x2", "x3").alias("vif")) + +# 2. Outliers — boolean mask from studentized residuals +df.select( + ps.residual_outliers("y", "x1", "x2", threshold=2.5).alias("outliers") +) + +# 3. Influence — Cook's distance per observation +df.select(ps.cooks_distance("y", "x1", "x2").alias("d")) +``` + +The full battery — leverage, DFFITS, GLM residuals (Pearson/deviance/working), Pearson χ² goodness of fit, GVIF for grouped predictors — is documented on the [Diagnostics page](api/regression/diagnostics.md). + ## Model Classes For direct model access outside Polars expressions: diff --git a/docs/index.md b/docs/index.md index 4f6790e..0fc3cc2 100644 --- a/docs/index.md +++ b/docs/index.md @@ -44,6 +44,18 @@ result.with_columns( ) ``` +## Use from Rust + +`polars-statistics` builds as both a Python extension (`cdylib`) and a Rust library (`rlib`). Other Rust crates can depend on it directly and call the same statistical and regression code that the Python plugin uses — no Python boundary, no FFI overhead. + +```toml +[dependencies] +polars = { version = "0.52", features = ["lazy", "partition_by"] } +polars-statistics = { version = "0.5", default-features = false } +``` + +`default-features = false` disables the `python` feature, so `pyo3` and `numpy` are not linked. Every Polars expression has a public Rust counterpart named `_fit` (e.g. `ols_fit`, `vif_fit`, `logistic_predict_fit`) under `polars_statistics::expressions`. See the [Use from Rust](https://github.com/DataZooDE/polars-statistics#use-from-rust) section of the README for a full runnable example and the complete list of `_fit` entry points. + ## Examples | Example | Description |