From c38aecf679860abd0c55f1bcd55c0b152ab155b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Simon=20M=C3=BCller?= Date: Thu, 28 May 2026 15:22:26 +0200 Subject: [PATCH] Influence diagnostics: dffits + influence masks (#27 batch 3) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add 4 influence diagnostics following the established pattern: - dffits_fit + pl_dffits: scaled per-row influence on fitted value. - influential_cooks_fit + pl_influential_cooks: boolean mask of observations above Cook's distance threshold (default 4/n). - influential_dffits_fit + pl_influential_dffits: boolean mask of observations above DFFITS threshold (default 2 * sqrt(p/n)). - high_leverage_points_fit + pl_high_leverage_points: boolean mask of observations with leverage > threshold (default 2*p/n). Takes only X (no y). Each fits OLS internally via fit_ols_for_residual_diag (except high_leverage_points which only builds the design matrix). Two new struct output schemas: dffits_output_dtype{dffits: List, n_observations} influence_mask_output_dtype{is_influential: List, n_influential, n_observations} Python (exprs/regression.py): dffits, influential_cooks, influential_dffits, high_leverage_points — all exported. Tests: tests/test_influence_diagnostics.py: 7 pytests covering shape and detection of planted high-influence points (large DFFITS, influential masks flag the planted indices, high_leverage_points catches extreme x values, custom thresholds tighten/loosen masks). tests/rust_api.rs: smoke tests in diagnostic_fits for all four *_fit entry points. Full Python suite: 448/448 pass (441 prior + 7 new). Rust integration: 12/12 pass. Co-Authored-By: Claude Opus 4.7 (1M context) --- python/polars_statistics/__init__.py | 8 + python/polars_statistics/exprs/__init__.py | 8 + python/polars_statistics/exprs/regression.py | 107 +++++++++++ src/expressions/regression.rs | 190 ++++++++++++++++++- tests/rust_api.rs | 37 ++++ tests/test_influence_diagnostics.py | 90 +++++++++ 6 files changed, 439 insertions(+), 1 deletion(-) create mode 100644 tests/test_influence_diagnostics.py diff --git a/python/polars_statistics/__init__.py b/python/polars_statistics/__init__.py index 2a2fb7f..6900d71 100644 --- a/python/polars_statistics/__init__.py +++ b/python/polars_statistics/__init__.py @@ -133,6 +133,10 @@ studentized_residuals, externally_studentized_residuals, residual_outliers, + dffits, + influential_cooks, + influential_dffits, + high_leverage_points, logistic_pearson_residuals, logistic_deviance_residuals, logistic_working_residuals, @@ -354,6 +358,10 @@ "studentized_residuals", "externally_studentized_residuals", "residual_outliers", + "dffits", + "influential_cooks", + "influential_dffits", + "high_leverage_points", "logistic_pearson_residuals", "logistic_deviance_residuals", "logistic_working_residuals", diff --git a/python/polars_statistics/exprs/__init__.py b/python/polars_statistics/exprs/__init__.py index 6fce7bd..56e8a02 100644 --- a/python/polars_statistics/exprs/__init__.py +++ b/python/polars_statistics/exprs/__init__.py @@ -89,6 +89,10 @@ studentized_residuals, externally_studentized_residuals, residual_outliers, + dffits, + influential_cooks, + influential_dffits, + high_leverage_points, logistic_pearson_residuals, logistic_deviance_residuals, logistic_working_residuals, @@ -260,6 +264,10 @@ "studentized_residuals", "externally_studentized_residuals", "residual_outliers", + "dffits", + "influential_cooks", + "influential_dffits", + "high_leverage_points", "logistic_pearson_residuals", "logistic_deviance_residuals", "logistic_working_residuals", diff --git a/python/polars_statistics/exprs/regression.py b/python/polars_statistics/exprs/regression.py index 95655df..d78b9cd 100644 --- a/python/polars_statistics/exprs/regression.py +++ b/python/polars_statistics/exprs/regression.py @@ -1124,6 +1124,113 @@ def _residual_diag_args(y, x, add_intercept): return [y.cast(pl.Float64), pl.lit(add_intercept, dtype=pl.Boolean), *x_exprs] +# ============================================================================ +# Influence diagnostics (issue #27 batch 3) +# ============================================================================ + + +def dffits( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + add_intercept: bool | None = None, + with_intercept: bool | None = None, +) -> pl.Expr: + """DFFITS per observation from an internal OLS fit. + + Measures the scaled change in fitted value if observation ``i`` were + dropped. |DFFITS_i| > 2 * sqrt(p / n) is the common influence cutoff. + + Returns a struct with ``dffits`` (List[float]) and ``n_observations``. + """ + add_intercept = _resolve_intercept(add_intercept, with_intercept) + return register_plugin_function( + plugin_path=LIB, + function_name="pl_dffits", + args=_residual_diag_args(y, x, add_intercept), + returns_scalar=True, + ) + + +def _influence_args(y, x, threshold, add_intercept): + if isinstance(y, str): + y = pl.col(y) + x_exprs = [(pl.col(xi) if isinstance(xi, str) else xi).cast(pl.Float64) for xi in x] + return [ + y.cast(pl.Float64), + pl.lit(add_intercept, dtype=pl.Boolean), + pl.lit(threshold, dtype=pl.Float64), + *x_exprs, + ] + + +def influential_cooks( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + threshold: float | None = None, + add_intercept: bool | None = None, + with_intercept: bool | None = None, +) -> pl.Expr: + """Boolean mask of observations with Cook's distance above threshold. + + Default threshold is ``4 / n`` (Belsley-Kuh-Welsch). Returns a struct + with ``is_influential`` (List[bool]), ``n_influential`` and + ``n_observations``. + """ + add_intercept = _resolve_intercept(add_intercept, with_intercept) + return register_plugin_function( + plugin_path=LIB, + function_name="pl_influential_cooks", + args=_influence_args(y, x, threshold, add_intercept), + returns_scalar=True, + ) + + +def influential_dffits( + y: Union[pl.Expr, str], + *x: Union[pl.Expr, str], + threshold: float | None = None, + add_intercept: bool | None = None, + with_intercept: bool | None = None, +) -> pl.Expr: + """Boolean mask of observations with |DFFITS| above threshold. + + Default threshold is ``2 * sqrt(p / n)``. Returns a struct with + ``is_influential`` (List[bool]), ``n_influential`` and ``n_observations``. + """ + add_intercept = _resolve_intercept(add_intercept, with_intercept) + return register_plugin_function( + plugin_path=LIB, + function_name="pl_influential_dffits", + args=_influence_args(y, x, threshold, add_intercept), + returns_scalar=True, + ) + + +def high_leverage_points( + *x: Union[pl.Expr, str], + threshold: float | None = None, + add_intercept: bool = True, +) -> pl.Expr: + """Boolean mask of high-leverage observations. + + Default threshold is ``2 * p / n``. Returns a struct with + ``is_influential`` (List[bool]), ``n_influential`` and ``n_observations``. + + Note: takes only feature columns (no y). + """ + x_exprs = [(pl.col(xi) if isinstance(xi, str) else xi).cast(pl.Float64) for xi in x] + return register_plugin_function( + plugin_path=LIB, + function_name="pl_high_leverage_points", + args=[ + pl.lit(add_intercept, dtype=pl.Boolean), + pl.lit(threshold, dtype=pl.Float64), + *x_exprs, + ], + returns_scalar=True, + ) + + def standardized_residuals( y: Union[pl.Expr, str], *x: Union[pl.Expr, str], diff --git a/src/expressions/regression.rs b/src/expressions/regression.rs index 293c147..b65f113 100644 --- a/src/expressions/regression.rs +++ b/src/expressions/regression.rs @@ -9,7 +9,8 @@ use serde::Deserialize; use anofox_regression::diagnostics::{ check_binary_separation, check_count_sparsity, compute_leverage, condition_diagnostic, - cooks_distance, externally_studentized_residuals, residual_outliers, standardized_residuals, + cooks_distance, dffits, externally_studentized_residuals, high_leverage_points, + influential_cooks, influential_dffits, residual_outliers, standardized_residuals, studentized_residuals, variance_inflation_factor, ConditionSeverity, SeparationCheck, SeparationType, }; @@ -221,6 +222,28 @@ fn outlier_mask_output_dtype(_input_fields: &[Field]) -> PolarsResult { Ok(Field::new("outliers".into(), DataType::Struct(fields))) } +/// Output type for DFFITS diagnostics (one value per row). +fn dffits_output_dtype(_input_fields: &[Field]) -> PolarsResult { + let fields = vec![ + Field::new("dffits".into(), DataType::List(Box::new(DataType::Float64))), + Field::new("n_observations".into(), DataType::UInt32), + ]; + Ok(Field::new("dffits_diag".into(), DataType::Struct(fields))) +} + +/// Output type for influence mask diagnostics: is_influential per row + counts. +fn influence_mask_output_dtype(_input_fields: &[Field]) -> PolarsResult { + let fields = vec![ + Field::new( + "is_influential".into(), + DataType::List(Box::new(DataType::Boolean)), + ), + Field::new("n_influential".into(), DataType::UInt32), + Field::new("n_observations".into(), DataType::UInt32), + ]; + Ok(Field::new("influence".into(), DataType::Struct(fields))) +} + // ============================================================================ // Helper Functions // ============================================================================ @@ -1630,6 +1653,171 @@ fn pl_residual_outliers(inputs: &[Series]) -> PolarsResult { residual_outliers_fit(inputs) } +// ============================================================================ +// Influence diagnostics (issue #27 batch 3) +// +// dffits, influential_cooks, influential_dffits, high_leverage_points. +// Each fits OLS internally (except high_leverage which only needs X). +// ============================================================================ + +fn dffits_output(values: Vec, n_obs: usize) -> PolarsResult { + let inner = Series::new("item".into(), values); + let d_s = Series::new("dffits".into(), &[inner]); + let n_s = Series::new("n_observations".into(), &[n_obs as u32]); + StructChunked::from_series("dffits_diag".into(), 1, [&d_s, &n_s].into_iter()) + .map(|ca| ca.into_series()) +} + +fn dffits_nan_output() -> PolarsResult { + dffits_output(vec![], 0) +} + +fn influence_mask_output( + mask: Vec, + n_influential: usize, + n_obs: usize, +) -> PolarsResult { + let inner = Series::new("item".into(), mask); + let m_s = Series::new("is_influential".into(), &[inner]); + let ni_s = Series::new("n_influential".into(), &[n_influential as u32]); + let nobs_s = Series::new("n_observations".into(), &[n_obs as u32]); + StructChunked::from_series("influence".into(), 1, [&m_s, &ni_s, &nobs_s].into_iter()) + .map(|ca| ca.into_series()) +} + +fn influence_mask_nan_output() -> PolarsResult { + influence_mask_output(vec![], 0, 0) +} + +/// DFFITS per row: scaled change in fitted value if observation i were dropped. +/// +/// Input contract: `[y, with_intercept (bool), x_0, ...]`. +pub fn dffits_fit(inputs: &[Series]) -> PolarsResult { + let (residuals, leverage, mse, n_params, n_rows) = match fit_ols_for_residual_diag(inputs) { + Some(v) => v, + None => return dffits_nan_output(), + }; + let d = dffits(&residuals, &leverage, mse, n_params); + let values: Vec = (0..d.nrows()).map(|i| d[i]).collect(); + dffits_output(values, n_rows) +} + +#[polars_expr(output_type_func=dffits_output_dtype)] +fn pl_dffits(inputs: &[Series]) -> PolarsResult { + dffits_fit(inputs) +} + +/// Influential observations by Cook's distance threshold (default 4/n). +/// +/// Input contract: `[y, with_intercept (bool), threshold (f64|null), x_0, ...]`. +pub fn influential_cooks_fit(inputs: &[Series]) -> PolarsResult { + if inputs.len() < 4 { + return influence_mask_nan_output(); + } + let threshold = inputs[2].f64()?.get(0); + // Build the cooks_distance contract from our own inputs. + let mut cd_inputs: Vec = Vec::with_capacity(2 + inputs.len() - 3); + cd_inputs.push(inputs[0].clone()); + cd_inputs.push(inputs[1].clone()); + cd_inputs.extend(inputs[3..].iter().cloned()); + + let (residuals, leverage, mse, n_params, n_rows) = match fit_ols_for_residual_diag(&cd_inputs) { + Some(v) => v, + None => return influence_mask_nan_output(), + }; + let cooks = cooks_distance(&residuals, &leverage, mse, n_params); + let idx = influential_cooks(&cooks, threshold); + let mut mask = vec![false; n_rows]; + for i in idx.iter() { + if *i < n_rows { + mask[*i] = true; + } + } + let n_inf = mask.iter().filter(|&&b| b).count(); + influence_mask_output(mask, n_inf, n_rows) +} + +#[polars_expr(output_type_func=influence_mask_output_dtype)] +fn pl_influential_cooks(inputs: &[Series]) -> PolarsResult { + influential_cooks_fit(inputs) +} + +/// Influential observations by DFFITS threshold (default 2 * sqrt(p/n)). +/// +/// Input contract: `[y, with_intercept (bool), threshold (f64|null), x_0, ...]`. +pub fn influential_dffits_fit(inputs: &[Series]) -> PolarsResult { + if inputs.len() < 4 { + return influence_mask_nan_output(); + } + let threshold = inputs[2].f64()?.get(0); + let mut cd_inputs: Vec = Vec::with_capacity(2 + inputs.len() - 3); + cd_inputs.push(inputs[0].clone()); + cd_inputs.push(inputs[1].clone()); + cd_inputs.extend(inputs[3..].iter().cloned()); + + let (residuals, leverage, mse, n_params, n_rows) = match fit_ols_for_residual_diag(&cd_inputs) { + Some(v) => v, + None => return influence_mask_nan_output(), + }; + let d = dffits(&residuals, &leverage, mse, n_params); + let idx = influential_dffits(&d, n_params, threshold); + let mut mask = vec![false; n_rows]; + for i in idx.iter() { + if *i < n_rows { + mask[*i] = true; + } + } + let n_inf = mask.iter().filter(|&&b| b).count(); + influence_mask_output(mask, n_inf, n_rows) +} + +#[polars_expr(output_type_func=influence_mask_output_dtype)] +fn pl_influential_dffits(inputs: &[Series]) -> PolarsResult { + influential_dffits_fit(inputs) +} + +/// High-leverage points (default threshold 2 * p / n). +/// +/// Input contract: `[with_intercept (bool), threshold (f64|null), x_0, ...]`. +/// No y is needed. +pub fn high_leverage_points_fit(inputs: &[Series]) -> PolarsResult { + if inputs.len() < 3 { + return influence_mask_nan_output(); + } + let with_intercept = inputs[0].bool()?.get(0).unwrap_or(true); + let threshold = inputs[1].f64()?.get(0); + let n_features = inputs.len() - 2; + let n_rows = inputs[2].len(); + if n_rows < 2 || n_features == 0 { + return influence_mask_nan_output(); + } + + // Build X matrix + let x = Mat::from_fn(n_rows, n_features, |row, col| { + inputs[2 + col] + .f64() + .ok() + .and_then(|ca| ca.get(row)) + .unwrap_or(0.0) + }); + let leverage = compute_leverage(&x, with_intercept); + let n_params = n_features + if with_intercept { 1 } else { 0 }; + let idx = high_leverage_points(&leverage, n_params, threshold); + let mut mask = vec![false; n_rows]; + for i in idx.iter() { + if *i < n_rows { + mask[*i] = true; + } + } + let n_inf = mask.iter().filter(|&&b| b).count(); + influence_mask_output(mask, n_inf, n_rows) +} + +#[polars_expr(output_type_func=influence_mask_output_dtype)] +fn pl_high_leverage_points(inputs: &[Series]) -> PolarsResult { + high_leverage_points_fit(inputs) +} + // ============================================================================ // GLM Expressions // ============================================================================ diff --git a/tests/rust_api.rs b/tests/rust_api.rs index ee43f74..82833d8 100644 --- a/tests/rust_api.rs +++ b/tests/rust_api.rs @@ -993,6 +993,43 @@ fn diagnostic_fits() { .unwrap_or(0); assert!(n_obs > 0); } + + // influence diagnostics (issue #27 batch 3). + { + let y_vals: Vec = x1 + .iter() + .enumerate() + .map(|(i, xi)| 0.5 + 1.5 * xi + 0.1 * (i as f64).sin()) + .collect(); + // dffits_fit: same contract as cooks_distance. + let dffits_inputs = vec![ + series_f64("y", &y_vals), + scalar_bool("with_intercept", true), + series_f64("x1", &x1), + series_f64("x2", &x2), + ]; + let _ = dffits_fit(&dffits_inputs).expect("dffits_fit failed"); + + // influential_cooks / influential_dffits: y, with_intercept, threshold, x... + let mask_inputs = vec![ + series_f64("y", &y_vals), + scalar_bool("with_intercept", true), + scalar_f64_null("threshold"), + series_f64("x1", &x1), + series_f64("x2", &x2), + ]; + let _ = influential_cooks_fit(&mask_inputs).expect("influential_cooks_fit failed"); + let _ = influential_dffits_fit(&mask_inputs).expect("influential_dffits_fit failed"); + + // high_leverage_points_fit: with_intercept, threshold, x... (no y) + let hlp_inputs = vec![ + scalar_bool("with_intercept", true), + scalar_f64_null("threshold"), + series_f64("x1", &x1), + series_f64("x2", &x2), + ]; + let _ = high_leverage_points_fit(&hlp_inputs).expect("high_leverage_points_fit failed"); + } } // ============================================================================= diff --git a/tests/test_influence_diagnostics.py b/tests/test_influence_diagnostics.py new file mode 100644 index 0000000..d640b41 --- /dev/null +++ b/tests/test_influence_diagnostics.py @@ -0,0 +1,90 @@ +"""Influence diagnostics tests (#27 batch 3).""" + +import numpy as np +import polars as pl +import pytest + +from polars_statistics import ( + dffits, + high_leverage_points, + influential_cooks, + influential_dffits, +) + + +@pytest.fixture +def outlier_data(): + """Linear data with planted high-influence points at indices [3, 17, 42].""" + rng = np.random.default_rng(1) + n = 100 + x = rng.normal(size=n) + y = 1.0 + 2.0 * x + rng.normal(scale=0.3, size=n) + # Make these both outliers AND high-leverage by yanking x as well. + for idx in (3, 17, 42): + x[idx] = 6.0 + y[idx] += 25.0 + return pl.DataFrame({"y": y, "x": x}) + + +class TestDffits: + def test_returns_one_per_row(self, outlier_data): + result = outlier_data.select(dffits("y", "x").alias("d")).item() + assert result["n_observations"] == 100 + assert len(result["dffits"]) == 100 + + def test_planted_indices_have_large_dffits(self, outlier_data): + result = outlier_data.select(dffits("y", "x").alias("d")).item() + for idx in (3, 17, 42): + assert abs(result["dffits"][idx]) > 1.0, ( + f"DFFITS at {idx} = {result['dffits'][idx]:.2f}" + ) + + +class TestInfluentialCooks: + def test_flags_planted_indices(self, outlier_data): + result = outlier_data.select( + influential_cooks("y", "x").alias("i") + ).item() + for idx in (3, 17, 42): + assert result["is_influential"][idx], f"missed Cook's at {idx}" + assert result["n_influential"] >= 3 + assert result["n_observations"] == 100 + + def test_custom_threshold(self, outlier_data): + """A higher threshold flags fewer points than the default.""" + default = outlier_data.select( + influential_cooks("y", "x").alias("i") + ).item() + strict = outlier_data.select( + influential_cooks("y", "x", threshold=1.0).alias("i") + ).item() + assert strict["n_influential"] <= default["n_influential"] + + +class TestInfluentialDffits: + def test_flags_planted_indices(self, outlier_data): + result = outlier_data.select( + influential_dffits("y", "x").alias("i") + ).item() + for idx in (3, 17, 42): + assert result["is_influential"][idx], f"missed DFFITS at {idx}" + assert result["n_observations"] == 100 + + +class TestHighLeveragePoints: + def test_flags_extreme_x_values(self, outlier_data): + """Planted x = 6.0 values should be flagged as high leverage.""" + result = outlier_data.select( + high_leverage_points("x").alias("h") + ).item() + for idx in (3, 17, 42): + assert result["is_influential"][idx], f"missed leverage at {idx}" + assert result["n_observations"] == 100 + + def test_no_y_argument_works(self): + """high_leverage_points takes only feature columns.""" + rng = np.random.default_rng(0) + df = pl.DataFrame({"x1": rng.normal(size=50), "x2": rng.normal(size=50)}) + result = df.select(high_leverage_points("x1", "x2").alias("h")).item() + assert result["n_observations"] == 50 + assert len(result["is_influential"]) == 50