`statspai.robustness`¶

robustness ¶

Robustness analysis tools.

spec_curve: Specification Curve Analysis (Simonsohn et al. 2020)
robustness_report: Automated battery of robustness checks
subgroup_analysis: Subgroup heterogeneity analysis with forest plot
refute: rerun an estimator on altered data whose answer is known

RefutationResult `dataclass` ¶

Bases: ResultProtocolMixin

Outcome of one refutation test.

Attributes:

Name	Type	Description
`method`	`str`	Which test was run.
`estimate`	`float`	The estimate on the data as they are.
`new_effect`	`float`	Mean of the estimates on the altered data.
`expected`	`float`	What the estimator should return on the altered data: zero for `'placebo_treatment'` and `'dummy_outcome'`, the original estimate for the other two.
`p_value`	`float`	Two-sided Monte Carlo p-value that the altered-data estimates are centred on `expected`: twice the smaller of the shares at or below and at or above it, with one added to numerator and denominator. Small means the estimates sit to one side of where they should be.
`refuted`	`bool`	`p_value < alpha`.
`interval`	`tuple of float`	Central `1 - alpha` range of the altered-data estimates.
`permutation_pvalue`	`float or None`	`'placebo_treatment'` only: the share of placebo estimates at least as large in absolute value as the original one. A randomisation p-value for the effect itself, not a refutation.
`simulations`	`ndarray`	The altered-data estimates.
`n_simulations, n_failed`	`int`	Reruns that returned an estimate, and reruns that raised.
`details`	`dict`	`first_error` when a rerun failed; `subset_fraction` for `'data_subset'`.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.normal(size=400)
>>> d = (x + rng.normal(size=400) > 0).astype(int)
>>> df = pd.DataFrame({'y': 2 * d + x + rng.normal(size=400), 'd': d, 'x': x})
>>> res = sp.refute(sp.aipw, df, y='y', treat='d', covariates=['x'],
...                 method='placebo_treatment', n_simulations=40, seed=0)
>>> isinstance(res, sp.RefutationResult), res.refuted
(True, False)

summary ¶

summary() -> str

A formatted account of the test.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(1)
>>> d = rng.integers(0, 2, 300)
>>> df = pd.DataFrame({'y': d + rng.normal(size=300), 'd': d,
...                    'x': rng.normal(size=300)})
>>> res = sp.refute(sp.aipw, df, y='y', treat='d', covariates=['x'],
...                 method='data_subset', n_simulations=40, seed=0)
>>> print(res.summary().splitlines()[1])
  Refutation: data_subset (random subsets of the sample)

RobustnessResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for robustness report results.

Returned by :func:robustness_report. Holds one row per robustness check in results_df (alternative controls, clustering, winsorising, trimming, subsamples), alongside the baseline_estimate / baseline_se for the key variable x. Call :meth:summary for a stability assessment.

Examples:

>>> import statspai as sp
>>> df = sp.cps_wage()
>>> report = sp.robustness_report(
...     df,
...     formula="log_wage ~ education + experience + tenure",
...     x="education",
...     winsor_levels=[0.01],
... )
>>> type(report).__name__
'RobustnessResult'
>>> report.x
'education'
>>> bool(report.n_checks > 0 and "check" in report.results_df.columns)
True

results_df `instance-attribute` ¶

results_df: DataFrame

One row per robustness check.

summary ¶

summary() -> str

Formatted text summary.

to_latex ¶

to_latex(*, caption: Optional[str] = 'Robustness Checks', label: Optional[str] = None) -> str

Export to LaTeX table.

plot ¶

plot(figsize: Optional[Tuple[float, float]] = None, title: Optional[str] = None, color: str = '#2C3E50', baseline_color: str = '#E74C3C') -> Any

Forest-plot style visualization of robustness checks.

Returns:

Type	Description
`fig, ax : matplotlib Figure and Axes`

FrontierSensitivityResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for frontier sensitivity analysis.

Returned by :func:copula_sensitivity, :func:survival_sensitivity, and :func:calibrate_confounding_strength.

Examples:

>>> import statspai as sp
>>> res = sp.copula_sensitivity(0.3, 0.1)
>>> isinstance(res, sp.FrontierSensitivityResult)
True
>>> res.method
'copula_gaussian'

SpecCurveResult `dataclass` ¶

Bases: ResultProtocolMixin

Holds all specification curve outputs.

Returned by :func:sp.spec_curve. Carries one row per specification in results_df plus summary statistics (median_estimate, share_significant, share_positive) and a .plot() method producing the two-panel Simonsohn-Simmons-Nelson figure.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> df = pd.DataFrame({
...     "education": rng.normal(12, 3, size=n),
...     "experience": rng.normal(10, 5, size=n),
...     "female": rng.integers(0, 2, size=n),
... })
>>> df["wage"] = (5.0 + 0.4 * df["education"] + 0.1 * df["experience"]
...               - 0.5 * df["female"] + rng.normal(size=n))
>>> res = sp.spec_curve(
...     data=df, y="wage", x="education",
...     controls=[[], ["experience"], ["experience", "female"]],
...     se_types=["nonrobust", "hc1"],
... )
>>> type(res).__name__
'SpecCurveResult'
>>> res.n_specs                     # 3 control sets x 2 SE types
6
>>> res.results_df.shape[0] == res.n_specs
True

References

simonsohn2020specification

results_df `instance-attribute` ¶

results_df: DataFrame

One row per specification with columns: spec_id, estimate, se, ci_lower, ci_upper, pvalue, significant, plus one column per choice dimension.

model_info `class-attribute` `instance-attribute` ¶

model_info: Dict[str, Any] = field(default_factory=dict)

multiplicity: the session count of specifications of y (this curve's included), see core/_spec_counter.py. dropped_specs: the specifications left out because their regressors are linearly dependent, one dict per specification with its choices and the redundant columns (collinear, reason).

summary ¶

summary(alpha: float = 0.05) -> str

Return a formatted text summary.

to_latex ¶

to_latex(*args: Any, caption: Optional[str] = 'Specification Curve Summary', label: Optional[str] = None) -> str

Export summary to a LaTeX table.

to_dataframe ¶

to_dataframe() -> DataFrame

Return the full results DataFrame.

plot ¶

plot(alpha: float = 0.05, color_sig: str = '#2C3E50', color_nonsig: str = '#BDC3C7', figsize: Optional[Tuple[float, float]] = None, title: Optional[str] = None, sort_by: str = 'estimate') -> Tuple[Any, Tuple[Any, Any]]

Draw the canonical two-panel specification curve plot.

Top panel: sorted point estimates with 95 % CIs, coloured by significance. Bottom panel: indicator matrix showing which analytical choices produced each specification.

Parameters:

Name	Type	Description	Default
`alpha`	`float`	Significance threshold for colouring.	`0.05`
`color_sig`	`str`	Colours for significant / non-significant estimates.	`'#2C3E50'`
`color_nonsig`	`str`	Colours for significant / non-significant estimates.	`'#2C3E50'`
`figsize`	`tuple`	Figure size (width, height). Auto-sized if None.	`None`
`title`	`str`	Title for the top panel.	`None`
`sort_by`	`str`	Column to sort specifications by. Default `'estimate'`.	`'estimate'`

Returns:

Type	Description
`fig, axes : matplotlib Figure and array of Axes`

SubgroupResult `dataclass` ¶

Bases: ResultProtocolMixin

Container for subgroup heterogeneity analysis.

Returned by :func:sp.subgroup_analysis. Holds one row per subgroup in results_df, the pooled overall_estimate / overall_se and per-grouping interaction Wald tests in het_tests.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> df = pd.DataFrame({
...     "education": rng.normal(12, 3, size=n),
...     "experience": rng.normal(10, 5, size=n),
...     "female": rng.integers(0, 2, size=n),
... })
>>> df["wage"] = (5.0 + 0.4 * df["education"] + 0.1 * df["experience"]
...               + 0.3 * df["female"] * df["education"]
...               + rng.normal(size=n))
>>> res = sp.subgroup_analysis(
...     data=df,
...     formula="wage ~ education + experience",
...     x="education",
...     by={"Gender": "female"},
... )
>>> type(res).__name__
'SubgroupResult'
>>> res.x
'education'
>>> "Gender" in res.het_tests
True

results_df `instance-attribute` ¶

results_df: DataFrame

Columns: group_var, group_val, estimate, se, ci_lower, ci_upper, pvalue, nobs, label, note. A subgroup whose regression cannot be run (too few rows, or x / a control constant or collinear within it) keeps its row with NaN estimates and the reason in note; note is empty for the others.

het_tests `instance-attribute` ¶

het_tests: Dict[str, Dict[str, float]]

Per group_var: chi2, pvalue, df.

to_latex ¶

to_latex(*, caption: Optional[str] = 'Subgroup Heterogeneity Analysis', label: Optional[str] = None) -> str

Export to LaTeX.

plot ¶

plot(figsize: Optional[Tuple[float, float]] = None, title: Optional[str] = None, color: str = '#2C3E50', overall_color: str = '#E74C3C') -> Any

Forest plot of subgroup estimates.

Returns:

Type	Description
`(fig, ax)`

SensitivityDashboard `dataclass` ¶

Result of a unified sensitivity analysis.

Always contains an e_value entry; other entries are optional depending on what the estimator provides.

to_dict ¶

to_dict(detail: str = 'agent') -> Dict[str, Any]

JSON-serialisable view, for MCP / agent callers.

Without this the generic serialiser fell back to hunting for estimate / se fields, found none on a dashboard, and emitted an empty payload — so the MCP sensitivity tool returned nothing but its own title and an agent calling it got no E-value, no Oster delta, no robustness value at all.

calibrate_confounding_strength ¶

calibrate_confounding_strength(estimate: float, se: float, *, observed_r2_outcome: float, observed_r2_treatment: float, dof: Optional[float] = None, alpha: float = 0.05, target_estimate: float = 0.0, multipliers: Optional[Sequence[float]] = None) -> FrontierSensitivityResult

Calibrate the strength of an unobserved confounder required to explain the observed effect to a target value.

Cinelli & Hazlett's (2020) benchmark bounds: a confounder Z that is k times as strong as an observed benchmark covariate X_j (k times its partial R² with the treatment and with the outcome) has at most

R²_{D~Z|X} = k R²_{D~Xj|X-j} / (1 - R²_{D~Xj|X-j}),
R²_{Y~Z|D,X} = ((√k + √r) / √(1 - r))² R²_{Y~Xj|D,X-j} / (1 - R²_{Y~Xj|D,X-j}) with r = k R²_{D~Xj}² / ((1 - k R²_{D~Xj})(1 - R²_{D~Xj})),

and biases the OLS coefficient by at most se · √dof · √(R²_{Y~Z|D,X} R²_{D~Z|X} / (1 - R²_{D~Z|X})). This is R sensemakr::ovb_bounds / ovb_partial_r2_bound / adjusted_estimate with kd = ky = k, generalised to an arbitrary target_estimate (the estimate is moved toward the target).

Parameters:

Name	Type	Description	Default
`estimate`	`float`	OLS coefficient of the treatment and its (classical) standard error.	required
`se`	`float`	OLS coefficient of the treatment and its (classical) standard error.	required
`observed_r2_outcome`	`float in (0, 1)`	Partial R² of the benchmark covariate with the outcome, given the treatment and the other covariates (sensemakr `r2yxj.dx`).	required
`observed_r2_treatment`	`float in (0, 1)`	Partial R² of the benchmark covariate with the treatment, given the other covariates (sensemakr `r2dxj.x`).	required
`dof`	`float`	Residual degrees of freedom of the regression. Required: the bias bound scales with `√dof`.	`None`
`alpha`	`float`	Level of the adjusted confidence interval, referred to t with `dof` degrees of freedom as R `sensemakr::ovb_bounds` reports it (sensemakr's stand-alone `adjusted_ci` uses `dof - 1`).	`0.05`
`target_estimate`	`float`	Effect value to explain away.	`0.0`
`multipliers`	`sequence of float`	Values of `k`; default `np.linspace(0.5, 5.0, 19)`. A `k` that implies `R²_{D~Z\|X} >= 1` is impossible and its row is NaN with `feasible = False` (sensemakr stops with an error instead).	`None`

Notes

Before 1.30 the bias bound omitted the √dof factor and scaled both partial R² linearly (k R²), so it was smaller than Cinelli and Hazlett's by roughly √dof -- a factor of ~20 at 400 observations -- and the function reported effects as robust that a confounder as strong as the benchmark explains away.

Examples:

>>> import statspai as sp
>>> res = sp.calibrate_confounding_strength(
...     0.3, 0.1, observed_r2_outcome=0.1, observed_r2_treatment=0.1,
...     dof=500)
>>> res.method
'calibrate_confounding_strength'
>>> list(res.curve.columns)[:3]
['multiplier', 'r2_outcome', 'r2_treatment']

References

Baitairian et al. (arXiv:2510.16560, 2025). [@baitairian2025calibrating] Cinelli & Hazlett (JRSS-B 2020). [@cinelli2020making]

copula_sensitivity ¶

copula_sensitivity(estimate: float, se: float, *, sigma_u: float = 1.0, sigma_y: float = 1.0, rho_grid: Optional[Sequence[float]] = None, alpha: float = 0.05) -> FrontierSensitivityResult

Gaussian-copula sensitivity to unobserved confounding.

Under a Gaussian copula with correlation rho between U (one latent unit-level confounder) and Y (outcome), the bias in an OLS / DML point estimate scales linearly with rho:

bias(rho) = rho * sigma_u * sigma_y / sigma_D²  ≈ rho * sigma_u * sigma_y

under the normalisation sigma_D = 1. The adjusted estimate is estimate - bias(rho); we sweep rho on a grid to find the breakpoint rho* that zeros the effect: the grid value of smallest |rho| whose adjusted interval covers zero or lies on the other side of it.

Parameters:

Name	Type	Description	Default
`estimate`	`float`		required
`se`	`float`		required
`sigma_u`	`float`	Standard deviations of the latent confounder and the outcome. With default values the bias coefficient is numerically equal to `rho`, matching Chernozhukov-Cinelli-Hazlett's "percentile scaling."	`1.0`
`sigma_y`	`float`	Standard deviations of the latent confounder and the outcome. With default values the bias coefficient is numerically equal to `rho`, matching Chernozhukov-Cinelli-Hazlett's "percentile scaling."	`1.0`
`rho_grid`	`sequence of float`	Correlation grid. Defaults to `np.linspace(-0.5, 0.5, 21)`.	`None`
`alpha`	`float`		`0.05`

Returns:

Type	Description
`FrontierSensitivityResult`

Examples:

>>> import statspai as sp
>>> res = sp.copula_sensitivity(0.3, 0.1)
>>> res.method
'copula_gaussian'
>>> int(len(res.curve))
21

References

Balgi, Braun, Peña & Daoud (arXiv:2508.08752, 2025). [@balgi2025sensitivity]

survival_sensitivity ¶

survival_sensitivity(log_hr: float, se_log_hr: float, *, gamma_grid: Optional[Sequence[float]] = None, baseline_survival_t: float = 0.5, alpha: float = 0.05) -> FrontierSensitivityResult

Nonparametric sensitivity for survival / hazard-ratio outcomes.

Extends Rosenbaum's Gamma bounds to hazard ratios and converts them into shifted survival differences at a chosen time t.

Given an observed log hazard ratio log_hr with SE se_log_hr, bound the log-HR at sensitivity parameter Γ; the worst case moves it toward the null (no effect), the best case away from it:

log_hr_worst(Γ) = log_hr − sign(log_hr) log(Γ)
log_hr_best(Γ)  = log_hr + sign(log_hr) log(Γ)

and translate the worst case into a survival shift at time t using the proportional-hazards identity S_1(t) = S_0(t) ^ exp(log_hr).

Parameters:

Name	Type	Description	Default
`log_hr`	`float`		required
`se_log_hr`	`float`		required
`gamma_grid`	`sequence of float`	Gamma (≥ 1) values. Defaults to `np.linspace(1.0, 3.0, 21)`.	`None`
`baseline_survival_t`	`float`	Baseline S_0(t) used to report Δ survival at time t.	`0.5`
`alpha`	`float`		`0.05`

Notes

breakpoint is the smallest grid Γ whose worst-case confidence interval contains 0 or lies on the other side of it (with a small SE the interval can cross 0 between two grid values). The continuous crossing is the confidence limit nearest the null on the hazard-ratio scale, Γ* = exp(|log_hr| - z se) (when it exceeds 1) -- the bias factor whose E-value Γ* + sqrt(Γ*(Γ* - 1)) R EValue::evalues.HR(rare = TRUE) reports for the confidence limit.

Before 1.30 the worst case was always log_hr - log(Γ): for a protective effect (log_hr < 0) it moved away from the null, so no Γ ever overturned the effect and every protective hazard ratio was reported as robust.

Examples:

>>> import statspai as sp
>>> res = sp.survival_sensitivity(0.4, 0.15)
>>> res.method
'survival_gamma'
>>> int(len(res.curve))
21

References

Hu & Westling (arXiv:2511.01412, 2025). [@hu2025nonparametric]

subgroup_analysis ¶

subgroup_analysis(data: DataFrame, formula: str, x: str, by: Dict[str, str], robust: str = 'hc1') -> SubgroupResult

Run subgroup heterogeneity analysis with forest plot.

Estimate the effect of x on y within each subgroup defined by the variables in by, and test for heterogeneity using interaction-based Wald tests.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame`	Analysis dataset.	required
`formula`	`str`	Regression formula, e.g. `"wage ~ education + experience"`.	required
`x`	`str`	Key explanatory variable.	required
`by`	`dict[str, str]`	Mapping of display name → column name for grouping. Example: `{'Gender': 'female', 'Region': 'region'}`.	required
`robust`	`str`	Standard error type for subgroup regressions.	`'hc1'`

Returns:

Type	Description
`SubgroupResult`	Container with `.plot()`, `.summary()`, `.to_latex()`, `.results_df`.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> df = pd.DataFrame({
...     "female": rng.integers(0, 2, size=n),
...     "region": rng.integers(0, 3, size=n),
...     "education": rng.normal(12, 3, size=n),
...     "experience": rng.normal(10, 5, size=n),
... })
>>> df["wage"] = (
...     5 + 0.4 * df["education"] + 0.1 * df["experience"]
...     - 0.5 * df["female"] + rng.normal(size=n)
... )
>>> result = sp.subgroup_analysis(
...     data=df,
...     formula="wage ~ education + experience",
...     x='education',
...     by={'Gender': 'female', 'Region': 'region'},
... )
>>> fig = result.plot()
>>> type(result).__name__
'SubgroupResult'

statspai.robustness¶

robustness ¶

RefutationResult dataclass ¶

summary ¶

RobustnessResult dataclass ¶

results_df instance-attribute ¶

summary ¶

to_latex ¶

plot ¶

FrontierSensitivityResult dataclass ¶

SpecCurveResult dataclass ¶

results_df instance-attribute ¶

model_info class-attribute instance-attribute ¶

summary ¶

to_latex ¶

to_dataframe ¶

plot ¶

SubgroupResult dataclass ¶

results_df instance-attribute ¶

het_tests instance-attribute ¶

to_latex ¶

plot ¶

SensitivityDashboard dataclass ¶

to_dict ¶

calibrate_confounding_strength ¶

copula_sensitivity ¶

survival_sensitivity ¶

subgroup_analysis ¶

`statspai.robustness`¶

RefutationResult `dataclass` ¶

RobustnessResult `dataclass` ¶

results_df `instance-attribute` ¶

FrontierSensitivityResult `dataclass` ¶

SpecCurveResult `dataclass` ¶

results_df `instance-attribute` ¶

model_info `class-attribute` `instance-attribute` ¶

SubgroupResult `dataclass` ¶

results_df `instance-attribute` ¶

het_tests `instance-attribute` ¶

SensitivityDashboard `dataclass` ¶