statspai.diagnostics

diagnostics

Diagnostics and sensitivity analysis for StatsPAI.

Provides: - Oster (2019) coefficient stability bounds - McCrary (2008) density discontinuity test for RD manipulation

KitagawaResult dataclass

Bases: ResultProtocolMixin

Result of the Kitagawa (2015) LATE validity test.

Attributes:

Name	Type	Description
statistic	float	Test statistic (maximum violation of the inequality conditions).
p_value	float	Bootstrap p-value.
first_stage	float	First-stage effect: P(D=1|Z=1) - P(D=1|Z=0).
n_boot	int	Number of bootstrap replications used.
n_obs	int	Number of observations.
max_violation_treated	float	Maximum violation of the treated-complier CDF condition.
max_violation_untreated	float	Maximum violation of the untreated-complier CDF condition.

Methods:

Name	Description
summary	Print a formatted summary.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(4)
>>> z = rng.integers(0, 2, 400)
>>> d = ((0.3 + 0.5 * z + rng.normal(0, 0.2, 400)) > 0.5).astype(int)
>>> y = 1.0 + 0.8 * d + rng.normal(size=400)
>>> df = pd.DataFrame({"outcome": y, "treated": d, "assigned": z})
>>> res = sp.kitagawa_test(data=df, y="outcome", treatment="treated",
...                        instrument="assigned", n_boot=200, seed=42)
>>> type(res).__name__
'KitagawaResult'
>>> bool(0.0 <= res.p_value <= 1.0)
True

summary

summary() -> str

Return a formatted summary string.

RosenbaumResult dataclass

Bases: ResultProtocolMixin

Result container for :func:rosenbaum_bounds.

Examples:

>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> treated = rng.normal(1.0, 1.0, 80)
>>> control = rng.normal(0.0, 1.0, 80)
>>> res = sp.rosenbaum_bounds(treated, control)
>>> type(res).__name__
'RosenbaumResult'
>>> res.n_pairs
80

to_dict

to_dict() -> dict

JSON-safe dict of every field (agent-native serialization).

EvidenceFactorsResult dataclass

Bases: ResultProtocolMixin

Result of :func:evidence_factors.

Attributes:

Name	Type	Description
pvalue	float	Truncated product of the two factor p-value bounds.
pvalue_treated_vs_control1, pvalue_control2_vs_others	float	The two factors.
detail	DataFrame	Statistic, bounding moments and deviate of each factor.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = rng.normal(size=(300, 3))
>>> y[:, 0] += 0.8
>>> res = sp.evidence_factors(y, gamma=1.5, upsilon=1.5)
>>> type(res).__name__
'EvidenceFactorsResult'
>>> bool(res.pvalue < 0.05)
True

SensitivityTestResult dataclass

Bases: ResultProtocolMixin

Result of :func:rosenbaum_stratified and :func:noether_test.

Attributes:

Name	Type	Description
pvalue	float	Upper bound on the p-value at the first gamma.
gamma	float	The first (or only) Gamma analysed.
gamma_critical	float	The Gamma at which the bound equals alpha.
statistic	float	The test statistic.
detail	DataFrame	One row per Gamma.
diagnostics	dict	Counts that describe what entered the test.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> d = rng.normal(0.5, 1.0, 300)
>>> res = sp.noether_test(d, gamma=2.0)
>>> type(res).__name__
'SensitivityTestResult'
>>> res.diagnostics["n_pairs_used"]
101

WeakRobustResult

Bases: ResultProtocolMixin

Container holding the unified weak-IV-robust panel.

Returned by :func:sp.weakrobust. Supports dict-style lookup (panel["effective_F"]), .get(), .to_dict(), .to_frame() (Stata-style single table) and .summary().

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> z1 = rng.normal(size=n)
>>> z2 = rng.normal(size=n)
>>> u = rng.normal(size=n)
>>> educ = 0.7 * z1 + 0.5 * z2 + 0.5 * u + rng.normal(size=n)
>>> df = pd.DataFrame({
...     "wage": 1.0 + 0.4 * educ + u + rng.normal(size=n),
...     "educ": educ, "z1": z1, "z2": z2,
... })
>>> panel = sp.weakrobust(df, y="wage", endog="educ",
...                       instruments=["z1", "z2"],
...                       clr_simulations=2000, grid_size=101,
...                       random_state=0)
>>> type(panel).__name__
'WeakRobustResult'
>>> panel.endog
'educ'
>>> "effective_F" in panel.to_dict()
True

to_frame

to_frame() -> DataFrame

Return the Stata-style single-table summary.

WeightedRankPowerResult dataclass

Bases: ResultProtocolMixin

Result of :func:weighted_rank_power.

Attributes:

Name	Type	Description
power	DataFrame	Estimated power of the level-alpha sensitivity analysis at each Gamma.
mean, variance	float	Jackknife estimates of the mean of the per-block statistic and of the variance of that mean in a study of the observed size.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = rng.normal(size=(150, 3))
>>> y[:, 0] += 1.0
>>> res = sp.weighted_rank_power(y, gammas=[1.0, 2.0, 4.0])
>>> res.power["power"].is_monotonic_decreasing
True

WeightedRankResult dataclass

Bases: ResultProtocolMixin

Result of :func:weighted_rank.

Attributes:

Name	Type	Description
pvalue	float	Upper bound on the p-value at the first gamma.
gamma	float	The first (or only) value of Gamma analysed.
gamma_critical	float	The Gamma at which the bound equals alpha: the study is insensitive to biases smaller than this. 1.0 when the test does not reject even without hidden bias.
statistic, expectation, variance, deviate	float	The weighted rank statistic, the bounds on its null moments, and the standardised deviate, at the first gamma.
detail	DataFrame	One row per Gamma (and per phi when several are given).
estimate, conf_int	tuple of float, optional	With estimates=True, the range of Hodges-Lehmann point estimates and the confidence interval for an additive effect at the first gamma.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = rng.normal(size=(200, 3))
>>> y[:, 0] += 1.0
>>> res = sp.weighted_rank(y, gamma=2.0)
>>> type(res).__name__
'WeightedRankResult'
>>> res.n_blocks, res.block_size
(200, 3)
>>> bool(res.pvalue < 0.05)
True

diagnose_result

diagnose_result(result: Any, print_results: bool = True, alpha: float = 0.05) -> Dict[str, Any]

One-call diagnostic battery, auto-selected by method type.

Parameters:

Name	Type	Description	Default
result	EconometricResults or CausalResult	A fitted result from any StatsPAI estimator.	required
print_results	bool	If True, print formatted diagnostics to stdout.	True
alpha	float	Significance level for tests.	0.05

Returns:

Type	Description
dict	Keys depend on the method. Always includes 'method_type' and 'checks' (list of individual test result dicts).

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> rows = []
>>> for u in range(100):
...     treated = int(u < 50)
...     fe = rng.normal()
...     for post in (0, 1):
...         y = (1.0 + fe + 0.5 * post
...              + 0.8 * treated * post + rng.normal())
...         rows.append({"unit": u, "y": y,
...                      "treated": treated, "post": post})
>>> df = pd.DataFrame(rows)
>>> res = sp.did(df, y="y", treat="treated", time="post")
>>> out = sp.diagnose_result(res, print_results=False)
>>> out["method_type"]
'did'
>>> sorted(out.keys())
['checks', 'method_type', 'next_steps', 'violations']

binned_residuals

binned_residuals(x: Any, residuals: Optional[Any] = None, n_bins: Optional[int] = None, by: Optional[Any] = None, band: str = 'empirical') -> DataFrame

Average residuals within bins of the fitted value or a regressor.

The residuals of a binary regression take two values per fitted probability and a plot of them shows nothing. Averaged within bins holding equal numbers of observations they should scatter around zero inside +/- 2 sd / sqrt(n); a run of bins outside the band, or a trend across bins, shows where the mean model is wrong.

Parameters:

Name	Type	Description	Default
x	fitted model or array-like	A fitted model (sp.logit, sp.probit, sp.glm, sp.poisson, sp.bayes_regress): its fitted values are binned and its response residuals averaged. Or the values to bin on, with residuals.	required
residuals	array - like	Observed minus fitted. Required when x is an array.	None
n_bins	int	Number of bins. Default floor(sqrt(n)).	None
by	array - like	Bin on this variable instead of the fitted values of a model, e.g. a regressor, to see whether it needs a transformation.	None
band	('empirical', 'model')	'empirical': 2 sd / sqrt(n) of the residuals in the bin, the band of arm::binnedplot. Where the outcome hardly varies (fitted probabilities near 0 or 1) those residuals are nearly equal, the band collapses and the bin is flagged for no reason. 'model': 2 sqrt(sum p (1 - p)) / n, the standard error the fitted binary model implies for the bin average, which does not collapse. Needs a fitted model of a 0 / 1 outcome.	'empirical'

Returns:

Type	Description
DataFrame	One row per bin: xbar and ybar (means of the binning variable and of the residuals), n, x_lo, x_hi, two_se (2 sd / sqrt(n) of the residuals in the bin) and outside (|ybar| > two_se). attrs['share_outside'] is the share of bins outside the band; about 5 percent is expected from a correct model.

Warns:

Type	Description
UserWarning	Observations with a missing or infinite binning value or residual were dropped (with their number).

Notes

Bins are cut at the order statistics floor(n i / n_bins) of the binning variable and are closed on the right, so tied values stay together and bin sizes can differ by a few observations (the rule of R arm::binned.resids).

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=400)})
>>> df["y"] = (rng.uniform(size=400) < 1 / (1 + np.exp(-df["x"]))).astype(int)
>>> out = sp.binned_residuals(sp.logit("y ~ x", df))
>>> len(out)
20

References

gelman2006data, gelman2020regression

binned_residuals_plot

binned_residuals_plot(x: Any, residuals: Optional[Any] = None, n_bins: Optional[int] = None, by: Optional[Any] = None, ax: Any = None, xlabel: Optional[str] = None, band: str = 'empirical') -> Any

Plot of :func:binned_residuals with its +/- 2 se band.

Parameters are those of :func:binned_residuals (including band), plus an optional matplotlib ax and xlabel.

Returns:

Type	Description
Figure

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=400)})
>>> df["y"] = (rng.uniform(size=400) < 1 / (1 + np.exp(-df["x"]))).astype(int)
>>> fig = sp.binned_residuals_plot(sp.logit("y ~ x", df))

References

gelman2006data

confounder_adjust

confounder_adjust(effect: Any, *, confounder_outcome_effect: Any, exposure_confounder_effect: Any = None, exposed_prev: Any = None, unexposed_prev: Any = None, measure: str = 'coef', rare_outcome: bool = True) -> DataFrame

Adjust an observed effect for a specified unmeasured confounder.

Parameters:

Name	Type	Description	Default
effect	float or array - like	The observed estimate. Pass the point estimate and both confidence limits together to see all three move.	required
confounder_outcome_effect	float or array - like	Effect of the confounder on the outcome, on the scale of measure: a coefficient (per unit of the confounder) for 'coef', a ratio otherwise.	required
exposure_confounder_effect	float or array - like	Normal confounder with unit variance: the difference in its mean between the exposed and the unexposed.	None
exposed_prev	float or array - like	Binary confounder: its prevalence among the exposed and among the unexposed. Give both, instead of exposure_confounder_effect.	None
unexposed_prev	float or array - like	Binary confounder: its prevalence among the exposed and among the unexposed. Give both, instead of exposure_confounder_effect.	None
measure	('coef', 'rr', 'or', 'hr')	'coef': a difference (linear-model coefficient, risk difference, difference in means). 'rr', 'or', 'hr': risk, odds or hazard ratio.	'coef'
rare_outcome	bool	Odds and hazard ratios only. False converts the observed ratio and confounder_outcome_effect to approximate risk ratios first (sqrt for an odds ratio; VanderWeele's formula for a hazard ratio), which is advisable when the outcome is common (above roughly 15%). The output is then on the risk-ratio scale and effect_observed shows the converted value.	True

Returns:

Type	Description
DataFrame	One row per combination of inputs (arrays are broadcast), with effect_adjusted, effect_observed and the confounder specification used.

Notes

With d the mean difference of the confounder (exposed_prev - unexposed_prev when binary) and g its effect on the outcome:

• difference scale: adjusted = observed - g * d;
• ratio scale, normal confounder: adjusted = observed / g ** d;
• ratio scale, binary confounder: adjusted = observed * (1 + (g - 1) p0) / (1 + (g - 1) p1).

The formulas assume one confounder, independent of the measured covariates given the exposure, whose effect on the outcome is the same in both exposure groups.

Examples:

An observed coefficient of 6.58; a confounder one unit of which lowers the outcome by 2.3 and whose mean is 0.17 lower among the exposed:

>>> import statspai as sp
>>> out = sp.confounder_adjust(6.58, confounder_outcome_effect=-2.3,
...                            exposure_confounder_effect=-0.17)
>>> round(float(out['effect_adjusted'].iloc[0]), 3)
6.189

A risk ratio of 1.5 and a binary confounder (40% vs 10%) that raises risk by 80%:

>>> out = sp.confounder_adjust(1.5, confounder_outcome_effect=1.8,
...                            exposed_prev=0.4, unexposed_prev=0.1,
...                            measure='rr')
>>> round(float(out['effect_adjusted'].iloc[0]), 4)
1.2273

References

[@lin1998assessing], [@schlesselman1978assessing], [@mcgowan2022tipr]

confounder_tip

confounder_tip(effect: Any, *, confounder_outcome_effect: Any = None, exposure_confounder_effect: Any = None, exposed_prev: Any = None, unexposed_prev: Any = None, measure: str = 'coef', rare_outcome: bool = True) -> DataFrame

Unmeasured confounder that would move an effect to the null.

Leave exactly one property of the confounder unspecified and it is solved for; specify all of them and the result is how many independent confounders of that kind it would take.

Parameters:

Name	Type	Description	Default
effect	float or array - like	The observed estimate to tip. To ask what would make a result "not significant", pass the confidence limit closest to the null rather than the point estimate.	required
confounder_outcome_effect	float or array - like	Effect of the confounder on the outcome on the scale of measure.	None
exposure_confounder_effect	float or array - like	Normal confounder with unit variance: difference in its mean between the exposed and the unexposed.	None
exposed_prev	float or array - like	Binary confounder: prevalence among the exposed and unexposed.	None
unexposed_prev	float or array - like	Binary confounder: prevalence among the exposed and unexposed.	None
measure	('coef', 'rr', 'or', 'hr')	Scale of effect; the null is 0 for 'coef' and 1 for the ratios.	'coef'
rare_outcome	bool	See :func:confounder_adjust.	True

Returns:

Type

Description

DataFrame

effect_observed, the confounder specification with the solved value filled in, and n_unmeasured_confounders (1 when a property was solved for). A solved prevalence outside [0, 1] means no binary confounder with the stated properties can tip the effect; it is returned as NaN with a warning. So is the solved property of a confounder that has no effect on the outcome or no association with the exposure.

Examples:

How large a mean difference would a confounder with a -7 effect on the outcome need, to explain away a coefficient of 6.58?

>>> import statspai as sp
>>> tip = sp.confounder_tip(6.58, confounder_outcome_effect=-7)
>>> round(float(tip['exposure_confounder_effect'].iloc[0]), 2)
-0.94

How strongly would a binary confounder (50% vs 10%) have to raise risk to explain away a risk ratio of 1.2?

>>> tip = sp.confounder_tip(1.2, exposed_prev=0.5, unexposed_prev=0.1,
...                         measure='rr')
>>> round(float(tip['confounder_outcome_effect'].iloc[0]), 3)
1.526

References

[@lin1998assessing], [@mcgowan2022tipr], [@vanderweele2017sensitivity]

bias_factor

bias_factor(rr_eu: float, rr_ud: float) -> float

Confounding bias factor B (Ding & VanderWeele 2016).

Given the maximum risk ratios of an unmeasured confounder U with the exposure (rr_eu) and with the outcome (rr_ud), the factor by which they can bias a risk ratio is

B = (rr_eu * rr_ud) / (rr_eu + rr_ud - 1).

The E-value is the value e such that bias_factor(e, e) equals the observed risk ratio.

Examples:

>>> import statspai as sp
>>> round(sp.bias_factor(2.0, 2.0), 4)
1.3333
>>> round(sp.bias_factor(3.0, 3.0), 4)
1.8

evalue_from_result

evalue_from_result(result: Any, measure: str = 'SMD', rare: Optional[bool] = None, true: Optional[float] = None, *, sd: Optional[float] = None) -> Dict[str, Any]

Compute an E-value from a StatsPAI CausalResult object.

Parameters:

Name	Type	Description	Default
result	CausalResult	Result from any StatsPAI causal estimator exposing a scalar estimate (and ideally se / ci).	required
measure	str	How to interpret result.estimate. 'SMD' / 'MD' read it as a difference in means on the outcome's scale and divide it (with its SE and CI) by sd before converting to the risk-ratio scale; pass 'RR' / 'OR' / 'HR' for ratio estimates.	'SMD'
rare	Optional[bool]	Passed through to :func:evalue.	None
true	Optional[bool]	Passed through to :func:evalue.	None
sd	float	Standard deviation of the outcome, required for 'SMD' / 'MD': the E-value of a mean difference depends on it, and the raw ATT is not a standardised difference unless the outcome already has unit variance (then pass sd=1). Usually df[y].std(). Ignored for ratio measures.	None

Returns:

Type	Description
dict	Same structure as :func:evalue.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> rows = []
>>> for u in range(100):
...     treated = int(u < 50)
...     fe = rng.normal()
...     for post in (0, 1):
...         y = (1.0 + fe + 0.5 * post
...              + 0.8 * treated * post + rng.normal())
...         rows.append({"unit": u, "y": y,
...                      "treated": treated, "post": post})
>>> df = pd.DataFrame(rows)
>>> res = sp.did(df, y="y", treat="treated", time="post")
>>> ev = sp.evalue_from_result(res, measure="SMD", sd=df["y"].std())
>>> ev["measure"]
'SMD'
>>> ev["evalue_estimate"] > 1.0
True

evalue_rd

evalue_rd(n11: float, n10: float, n01: float, n00: float, true: float = 0.0, alpha: float = 0.05, grid: float = 0.0001) -> Dict[str, Any]

Exact E-value for a risk difference from a 2x2 table.

Uses the bound of Ding & VanderWeele (2016, Proposition A.8) on the population-standardised causal risk difference under confounding with bounding factor B,

RD_true >= (p1 - p0 * B) * (f + (1 - f) / B),

where p1 / p0 are the outcome risks in the exposed / unexposed and f is the exposure prevalence. The point E-value converts the bias factor that moves this bound to true (closed-form root of the quadratic in their Proposition A.11) to the E-value scale. The confidence-limit E-value uses the smallest bias factor on a grid 1, 1 + grid, 1 + 2*grid, ... at which the delta-method lower confidence limit of the bound reaches true (their equation A.5).

Parameters:

Name	Type	Description	Default
n11	float	Exposed cases and exposed non-cases.	required
n10	float	Exposed cases and exposed non-cases.	required
n01	float	Unexposed cases and unexposed non-cases.	required
n00	float	Unexposed cases and unexposed non-cases.	required
true	float	Reference risk difference (must be below the observed RD).	0.0
alpha	float	Significance level for the confidence-limit E-value.	0.05
grid	float	Step size of the bias-factor grid search for the CI E-value.	1e-4

Returns:

Type	Description
dict	evalue_estimate (point E-value), evalue_ci (E-value for the lower confidence limit; 1.0 when the interval already reaches true), rd (observed risk difference), bias_factor (the bias factor that moves the point estimate to true), measure ('RD') and true.

Notes

Clean-room implementation from the paper (spec: docs/dev/cleanroom/evalue.md); R EValue::evalues.RD (GPL-2) was used only as a black box and agrees to 1e-12. When no grid point up to the point-estimate bias factor qualifies (a coarse grid on a very precise table), that bias factor itself is used as the last grid point, so the answer still overshoots the exact root by less than one grid step; R returns NA there.

References

[@ding2016sensitivity]; [@vanderweele2017sensitivity].

Examples:

>>> import statspai as sp
>>> out = sp.evalue_rd(200, 150, 100, 250)
>>> round(out["evalue_estimate"], 6), round(out["evalue_ci"], 6)
(3.414214, 2.726306)

hausman_test

hausman_test(data: DataFrame, y: str, x: List[str], id: str, time: str, alpha: float = 0.05, sigmamore: bool = False) -> Dict[str, Any]

Hausman test for FE vs RE in panel data.

Equivalent to Stata's hausman fe re.

Parameters:

Name	Type	Description	Default
data	DataFrame	Panel data in long format.	required
y	str	Dependent variable.	required
x	list of str	Independent variables (time-varying).	required
id	str	Unit identifier.	required
time	str	Time period identifier.	required
alpha	float		0.05
sigmamore	bool	Stata's hausman fe re, sigmamore: both covariance matrices use the RE disturbance variance. Use it when the classical statistic is negative (recommendation == "inconclusive").	False

Returns:

Type	Description
dict	'statistic': Hausman chi² statistic 'df': degrees of freedom 'pvalue': p-value (nan when the statistic is negative) 'recommendation': 'FE', 'RE', or 'inconclusive' when V_FE - V_RE is not positive semi-definite and the statistic is negative (an AssumptionWarning is raised) 'psd_violation': whether V_FE - V_RE fails to be positive definite 'beta_fe', 'beta_re': coefficient vectors

Examples:

>>> import statspai as sp
>>> df = sp.dgp_panel(n_units=40, n_periods=8, seed=0)
>>> result = sp.hausman_test(df, y='y', x=['x'], id='unit', time='time')
>>> bool(result['recommendation'] in ('FE', 'RE'))
True
>>> bool(0.0 <= result['pvalue'] <= 1.0)
True

Notes

The test statistic is:

.. math:: H = (\hat{\beta}{FE} - \hat{\beta}{RE})' [V(\hat{\beta}{FE}) - V(\hat{\beta}{RE})]^{-1} (\hat{\beta}{FE} - \hat{\beta}{RE})

Under H0, H ~ χ²(k).

• Reject H0 (p < 0.05) → Use Fixed Effects
• Fail to reject → Use Random Effects (more efficient)

See Hausman (1978, Econometrica).

influence_measures

influence_measures(result: Any) -> DataFrame

Influence statistics of every observation of a linear regression.

With h the diagonal of the hat matrix, e the residual, s the root mean squared error on n - k degrees of freedom and s(i) the one of the fit without observation i:

============== ==================================================== leverage h rstandard e / (s sqrt(1 - h)) rstudent e / (s(i) sqrt(1 - h)) cooksd rstandard^2 h / (k (1 - h)) dfits rstudent sqrt(h / (1 - h)) welsch dfits sqrt((n - 1) / (1 - h)) covratio 1 / ((1 - h) ((n - k - 1 + rstudent^2) / (n - k))^k) stdp s sqrt(h), the standard error of the prediction stdf s sqrt(1 + h), of the forecast stdr s sqrt(1 - h), of the residual dfbeta_<x> the change of the coefficient of x when the observation is left out, in units of its standard error computed with s(i) ============== ====================================================

These are the definitions of Belsley, Kuh and Welsch (1980) that Stata's predict after regress and R's influence.measures follow.

Parameters:

Name	Type	Description	Default
result	EconometricResults	An unweighted fit of sp.regress.	required

Returns:

Type	Description
DataFrame	One row per observation of the estimation sample, in its order, with the columns above and fitted and residual.

Raises:

Type	Description
MethodIncompatibility	(a ValueError.) The fit is not a linear regression, or it was weighted: the stored design is unweighted and the residuals are not, so the statistics would mix the two.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=50)})
>>> df["y"] = 1 + 2 * df["x"] + rng.normal(size=50)
>>> out = sp.influence_measures(sp.regress("y ~ x", data=df))
>>> bool(np.isclose(out["leverage"].sum(), 2.0))
True

logit_gof

logit_gof(result: Any, groups: Optional[int] = None) -> Dict[str, Any]

Goodness of fit of a binary-outcome model.

Without groups, Pearson's chi-squared over the J covariate patterns, sum (y - m p)^2 / (m p (1 - p)), on J - k degrees of freedom. With groups = g, the Hosmer-Lemeshow test: the observations are sorted into g groups by the quantiles of the fitted probability (rows with equal probabilities stay together), and sum (O - E)^2 / (E (1 - E / n)) over the groups is referred to chi-squared with g - 2 degrees of freedom (Stata estat gof and estat gof, group(g)).

Parameters:

Name	Type	Description	Default
result	EconometricResults	A fit of sp.logit or sp.probit.	required
groups	int		None

Returns:

Type	Description
dict	statistic, df, pvalue, n, n_patterns or n_groups, and for the grouped test the table of observed and expected counts.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=300)})
>>> df["d"] = (rng.random(300) < 1 / (1 + np.exp(-df["x"]))).astype(int)
>>> out = sp.logit_gof(sp.logit("d ~ x", data=df), groups=10)
>>> int(out["df"])
8

logit_influence

logit_influence(result: Any) -> DataFrame

Pregibon's (1981) diagnostics after a logistic regression.

The statistics are computed for covariate patterns: the m rows that share the values of every regressor are one binomial observation with y successes and fitted probability p, and every row of a pattern gets the pattern's value. With V the model-based covariance of the coefficients:

============== ==================================================== residual Pearson, (y - m p) / sqrt(m p (1 - p)) hat m p (1 - p) x V x', the leverage rstandard residual / sqrt(1 - hat) deviance the signed square root of the pattern's deviance dx2 rstandard^2: the fall in Pearson's chi-squared when the pattern is left out (Hosmer-Lemeshow) ddeviance deviance^2 / (1 - hat): the fall in the deviance dbeta residual^2 hat / (1 - hat)^2: Pregibon's influence on the coefficients ============== ====================================================

These are the statistics of Stata's predict after logit / logistic.

Parameters:

Name	Type	Description	Default
result	EconometricResults	A fit of sp.logit.	required

Returns:

Type	Description
DataFrame	One row per observation of the estimation sample, with the columns above and pattern, m, successes, p.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.integers(0, 4, size=200)})
>>> df["d"] = (rng.random(200) < 0.2 + 0.15 * df["x"]).astype(int)
>>> out = sp.logit_influence(sp.logit("d ~ x", data=df))
>>> int(out["pattern"].nunique())
4

kitagawa_test

kitagawa_test(data: DataFrame, y: str, treatment: str, instrument: str, n_boot: int = 1000, n_grid: int = 100, seed: Optional[int] = None) -> KitagawaResult

Kitagawa (2015) specification test for the validity of LATE.

Tests whether the data are consistent with the LATE assumptions (independence, exclusion restriction, monotonicity) by checking that the implied complier potential-outcome CDFs are proper distribution functions.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
y	str	Outcome variable name.	required
treatment	str	Endogenous binary treatment variable (D).	required
instrument	str	Binary instrument variable (Z).	required
n_boot	int	Number of bootstrap replications for the p-value.	1000
n_grid	int	Number of grid points for evaluating the CDF conditions.	100
seed	int	Random seed for reproducibility.	None

Returns:

Type	Description
KitagawaResult	Test results including statistic, p-value, and first-stage effect.

Raises:

Type	Description
ValueError	If instrument or treatment are not binary, or first stage is zero.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> z = rng.integers(0, 2, size=800)  # binary instrument
>>> d = ((0.2 + 0.5 * z + rng.normal(scale=0.3, size=800)) > 0.5)
>>> d = d.astype(int)  # first stage: z shifts take-up
>>> y = 1.0 + 0.8 * d + rng.normal(size=800)
>>> df = pd.DataFrame(
...     {"outcome": y, "treated": d, "assigned": z})
>>> result = sp.kitagawa_test(
...     data=df,
...     y="outcome",
...     treatment="treated",
...     instrument="assigned",
...     n_boot=200,
...     seed=42,
... )
>>> type(result).__name__
'KitagawaResult'
>>> bool(0.0 <= result.p_value <= 1.0)
True

rosenbaum_bounds

rosenbaum_bounds(treated: Optional[Sequence[float]] = None, control: Optional[Sequence[float]] = None, *, data: Optional[DataFrame] = None, y: Optional[str] = None, treat: Optional[str] = None, pair_id: Optional[str] = None, method: str = 'wilcoxon', alternative: str = 'greater', gamma_grid: Optional[Sequence[float]] = None, alpha: float = 0.05, zero_method: str = 'pratt', estimates: bool = False) -> RosenbaumResult

Compute Rosenbaum bounds on a paired observational study.

You can pass either two parallel arrays (treated, control) of matched outcomes, or a long-format DataFrame with columns y, treat, pair_id.

Parameters:

Name	Type	Description	Default
treated	array - like	Outcome in the treated / control unit of each matched pair (same length). Ignored if data is provided.	None
control	array - like	Outcome in the treated / control unit of each matched pair (same length). Ignored if data is provided.	None
data	DataFrame	Long-format data. Must contain exactly two rows per pair_id, one with treat=1 and one with treat=0.	None
method	('wilcoxon', 'sign', 't')	Wilcoxon signed-rank bound (continuous), binomial sign test (robust / binary), or "t", the permutational t-test: the statistic is the sum of the positive pair differences, whose null distribution under bias Gamma has mean p * sum(|d|) and variance p (1 - p) * sum(d**2) with p = Gamma / (1 + Gamma). It is the test of the mean difference (R sensitivitymw::senmw(method="t") on pairs) and, unlike the rank tests, is driven by the largest differences.	"wilcoxon"
alternative	('greater', 'less', 'two-sided')	Direction of the alternative hypothesis for the treatment effect.	"greater"
gamma_grid	sequence of float	Gamma values (>= 1) over which to compute bounding p-values. Default: np.arange(1.0, 3.01, 0.1).	None
alpha	float	Significance level used to report gamma_critical.	0.05
zero_method	('pratt', 'wilcox')	Treatment of pairs with a zero difference in the Wilcoxon bound. "pratt" ranks them with the other pairs and gives them no weight (Rosenbaum's DOS2::senWilcox, Stata rbounds); "wilcox" drops them before ranking (rbounds::psens). Identical when no difference is zero. Ignored by the sign test, which always drops zeros.	"pratt"
estimates	bool	Also bound the effect itself (method="wilcoxon" only). For each Gamma the detail table gains the range of the Hodges-Lehmann point estimate of an additive effect, hl_lower to hl_upper, and the outer ends of its 1 - alpha confidence interval, ci_lower and ci_upper. At Gamma = 1 these are the usual Hodges-Lehmann estimate and interval. They are the columns t-hat+, t-hat-, CI+ and CI- of Stata rbounds; where the signed-rank statistic equals its target over a stretch of shifts, the middle of the stretch is reported, as there.	False

Returns:

Type	Description
RosenbaumResult	gamma_critical is the smallest Gamma in the grid at which the upper-bound p-value exceeds alpha. It is inf if the study is insensitive across the grid, and 1.0 if already sensitive.

See Also

weighted_rank : matched sets with several controls or several treated individuals, weights with higher design sensitivity, and the exact Gamma at which the bound crosses alpha. rosenbaum_stratified : strata of any size. noether_test : sign test on the pairs with the largest differences. amplify : reading Gamma as two odds ratios.

Examples:

Two parallel arrays of matched-pair outcomes (Wilcoxon bound):

>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> treated = rng.normal(1.0, 1.0, 80)
>>> control = rng.normal(0.0, 1.0, 80)
>>> res = sp.rosenbaum_bounds(treated, control)
>>> res.gamma_critical >= 1.0
True

A binary / robust sign-test bound on an insensitive example:

>>> import numpy as np
>>> t = np.array([2.0, 3.0, 4.0, 5.0])
>>> c = np.array([1.0, 1.0, 1.0, 1.0])
>>> bounds = sp.rosenbaum_bounds(t, c, method="sign", gamma_grid=[1.0, 2.0])
>>> bool(np.all(bounds.pvalue_upper >= bounds.pvalue_lower))
True

amplify

amplify(gamma: float, treatment_odds: Union[float, Sequence[float]]) -> Union[float, ndarray]

Rewrite a sensitivity parameter as two odds ratios.

A sensitivity analysis for matched pairs speaks of one number, Gamma: an unobserved covariate that could change the odds of treatment by that factor and that predicts the outcome perfectly. The same bound applies to every unobserved covariate that multiplies the odds of treatment by Lambda and the odds of a positive pair difference in outcomes by Delta, provided Gamma = (Lambda * Delta + 1) / (Lambda + Delta) [@rosenbaum2009amplification]. This function returns the Delta that goes with each Lambda.

Parameters:

Name	Type	Description	Default
gamma	float	Sensitivity parameter, > 1.	required
treatment_odds	float or sequence of float	Lambda, each > gamma.	required

Returns:

Type	Description
float or ndarray	Delta = (gamma * Lambda - 1) / (Lambda - gamma).

Examples:

>>> import statspai as sp
>>> sp.amplify(4, 7)
9.0
>>> sp.amplify(1.5, [2, 3, 4]).round(3).tolist()
[4.0, 2.333, 2.0]

References

[@rosenbaum2009amplification]

evidence_factors

evidence_factors(y: Any, *, gamma: float = 1.0, upsilon: float = 1.0, alternative: str = 'greater', trunc: float = 0.2, phi: Tuple[Any, Any] = ((8, 7, 8), (8, 8, 8)), scores: Sequence[float] = (1.0, 2.0, 5.0), block_scale: str = 'gap') -> EvidenceFactorsResult

Two evidence factors from blocks with a treated unit and two controls.

Each block holds one treated individual, one control from a first control group and one from a second. The first factor compares the treated individual with the first control, as a matched pair. The second compares the second control with the other two pooled. When the treatment has no effect the two p-values behave as if they came from unrelated studies, although they share data, and the biases that could explain away one comparison differ from those that could explain away the other [@rosenbaum2023second]. The factors are combined by the truncated product.

Parameters:

Name	Type	Description	Default
y	array-like of shape (I, 3)	Outcomes; columns are treated, first control, second control.	required
gamma	float	Bias allowed in the treated-versus-first-control comparison.	1.0
upsilon	float	Bias allowed in the comparison of the second control with the other two.	1.0
alternative	('greater', 'less')	Direction of the effect of treatment on the outcome.	"greater"
trunc	float	Truncation point of :func:truncated_product.	0.2
phi	pair	Block weights of the two factors, as in :func:weighted_rank.	((8, 7, 8), (8, 8, 8))
scores	sequence of 3 floats	Within-block scores of the second factor. Unequal steps guard against the dilution that pooling the treated individual with an unaffected control would otherwise cause.	(1, 2, 5)
block_scale	('gap', 'range')	Block dispersion used by the second factor; the first, on pairs, always uses the range.	"gap"

Returns:

Type	Description
EvidenceFactorsResult

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = rng.normal(size=(300, 3))
>>> y[:, 0] += 0.8
>>> res = sp.evidence_factors(y, gamma=1.5, upsilon=1.5)
>>> bool(res.pvalue <= max(res.pvalue_treated_vs_control1,
...                        res.pvalue_control2_vs_others))
True

References

[@rosenbaum2023second], [@rosenbaum2011some], [@zaykin2002truncated]

noether_test

noether_test(treated: Sequence[float], control: Optional[Sequence[float]] = None, *, f: float = 2.0 / 3.0, gamma: Union[float, Sequence[float]] = 1.0, alternative: str = 'greater', alpha: float = 0.05) -> SensitivityTestResult

Sign test on the matched pairs with the largest absolute differences.

Discards the fraction f of pairs whose absolute differences are smallest and applies the sign test to the rest [@noether1973some]. Small differences are nearly as likely to be positive as negative whether or not there is an effect, so they add noise that a small bias can exploit; dropping them makes the test less sensitive to unmeasured bias when the effect is not tiny [@rosenbaum2012exact]. The bound on the p-value is exact: under a bias of at most gamma the number of positive differences is at most binomial with success probability gamma / (1 + gamma).

Parameters:

Name	Type	Description	Default
treated	array - like	Treated outcome of each pair, or the treated-minus-control differences when control is omitted.	required
control	array - like	Control outcome of each pair.	None
f	float	Fraction of pairs set aside, 0 <= f < 1. Pairs are kept when the rank of their absolute difference, divided by the number of pairs, is at least f. 0 is the ordinary sign test.	2/3
gamma	float or sequence of float	Sensitivity parameter(s), >= 1.	1.0
alternative	('greater', 'less', 'two-sided')	"two-sided" doubles the smaller one-sided bound.	"greater"
alpha	float	Level for gamma_critical.	0.05

Returns:

Type	Description
SensitivityTestResult	statistic is the number of positive differences among the pairs used; diagnostics["n_pairs_used"] is how many those are.

Warns:

Type	Description
UserWarning	Pairs with a missing or infinite difference were dropped (with their number).

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> d = rng.normal(0.5, 1.0, 300)
>>> all_pairs = sp.noether_test(d, f=0, gamma=2.0)
>>> top_third = sp.noether_test(d, gamma=2.0)
>>> bool(top_third.gamma_critical > all_pairs.gamma_critical)
True

References

[@noether1973some], [@rosenbaum2012exact]

rosenbaum_stratified

rosenbaum_stratified(data: DataFrame, y: str, treat: str, strata: Optional[str] = None, *, gamma: Union[float, Sequence[float]] = 1.0, score: str = 'rank', alternative: str = 'greater', alpha: float = 0.05, bound: str = 'taylor') -> SensitivityTestResult

Rosenbaum bounds for a comparison within strata of any size.

Each stratum holds some treated individuals and some controls, in any numbers; the statistic is the sum of the treated individuals' scores. Under the null hypothesis of no effect, and allowing two individuals in the same stratum to differ in their odds of treatment by a factor of at most gamma, the function bounds the one-sided p-value [@rosenbaum2018sensitivity]. Matched pairs and matched sets are special cases; with strata=None it is the two-sample analysis of [@rosenbaum1990sensitivity].

Parameters:

Name	Type	Description	Default
data	DataFrame	One row per individual.	required
y	str	Outcome column (or, with score="raw", a column of scores).	required
treat	str	0/1 treatment indicator.	required
strata	str	Stratum identifier. None treats the sample as one stratum.	None
gamma	float or sequence of float	Sensitivity parameter(s), >= 1.	1.0
score	('rank', 'stratum_rank', 'aligned_rank', 'raw')	"rank" ranks the outcomes over the whole sample (the Wilcoxon rank sum in a single stratum). "stratum_rank" ranks within strata. "aligned_rank" subtracts the stratum mean and then ranks over the whole sample, the aligned ranks of Hodges and Lehmann. "raw" uses y itself, for the permutational t-test or for scores computed elsewhere.	"rank"
alternative	('greater', 'less', 'two-sided')		"greater"
alpha	float	Level for gamma_critical and for the Taylor bound.	0.05
bound	('taylor', 'separable')	"separable" takes, stratum by stratum, the hidden bias that maximises the null expectation; it is exact as the number of strata grows but may be slightly anti-conservative with a few large strata. "taylor" is a bound that holds at level alpha whatever the strata look like. Both are in detail. With a single stratum neither is needed: the candidate biases are searched for the largest p-value directly.	"taylor"

Returns:

Type	Description
SensitivityTestResult

Warns:

Type	Description
UserWarning	Rows with a missing value in a column used were dropped (with their number).

Notes

The moments of the statistic are computed exactly from Fisher's noncentral hypergeometric distribution. R's senstrat obtains them from the BiasedUrn package at its default precision of 1e-7 (method="BU"), or from a large-sample approximation (method="LS"); its method="RK" is exact and is the one this function agrees with to rounding.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> s = rng.integers(0, 6, 300)
>>> z = rng.binomial(1, 0.3 + 0.05 * s)
>>> yv = 0.6 * z + 0.3 * s + rng.normal(size=300)
>>> df = pd.DataFrame({"y": yv, "z": z, "s": s})
>>> res = sp.rosenbaum_stratified(df, "y", "z", "s", gamma=[1.0, 1.5],
...                               score="aligned_rank")
>>> res.diagnostics["n_strata"]
6
>>> bool(res.detail["pvalue"].is_monotonic_increasing)
True

References

[@rosenbaum2018sensitivity], [@rosenbaum1990sensitivity], [@gastwirth2000asymptotic]

truncated_product

truncated_product(pvalues: Sequence[float], trunc: float = 0.2) -> float

Combine independent p-values by the truncated product method.

The statistic is the product of the p-values that are at most trunc; its null distribution is known in closed form [@zaykin2002truncated]. trunc=1 is Fisher's method. Truncation keeps one large p-value from undoing several small ones, which matters when the p-values are upper bounds from a sensitivity analysis: such bounds tend to one as the allowed bias grows [@hsu2013effect].

Parameters:

Name	Type	Description	Default
pvalues	sequence of float	Independent p-values (or p-values whose joint distribution is stochastically larger than uniform on the unit cube, as with evidence factors).	required
trunc	float	Truncation point, 0 < trunc <= 1.	0.2

Returns:

Type	Description
float	The combined p-value; 1.0 when no p-value is at most trunc.

Examples:

>>> import statspai as sp
>>> round(sp.truncated_product([0.01, 0.3]), 4)
0.0399
>>> round(sp.truncated_product([0.01, 0.3], trunc=1.0), 4)  # Fisher
0.0204
>>> sp.truncated_product([0.5, 0.3])
1.0

References

[@zaykin2002truncated], [@hsu2013effect]

mccrary_test

mccrary_test(data: DataFrame, x: str, c: float = 0, bw: Optional[float] = None, n_bins: Optional[int] = None, alpha: float = 0.05, bin_width: Optional[float] = None) -> CausalResult

McCrary (2008) density discontinuity test for RD manipulation.

Tests the null hypothesis that the density of the running variable is continuous at the cutoff. A rejection suggests possible manipulation of the running variable.

A line-for-line port of R rdd::DCdensity (the R implementation of McCrary's DCdensity): a fine histogram with bin width 2 sd(x) n^(-1/2), a bandwidth from fourth-order polynomial fits to the bin heights on each side (3.348 (sigma^2 (c - l) / sum f''^2)^(1/5), averaged over the two sides), zero-padded bins, triangular-kernel local linear fits to the bin heights on each side, and theta = log f+ - log f- with se = sqrt(24 / (5 n h) (1/f+ + 1/f-)).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
x	str	Running variable name.	required
c	float	RD cutoff.	0
bw	float	Bandwidth of the local linear fits. Default: McCrary's rule above.	None
n_bins	int	Convenience alternative to bin_width: the bin width becomes (max(x) - min(x)) / n_bins. Ignored when bin_width is given.	None
alpha	float	Significance level.	0.05
bin_width	float	Histogram bin width (DCdensity's bin). Default: 2 sd(x) n^(-1/2).	None

Returns:

Type	Description
CausalResult	estimate: log density difference ln(f̂₊) - ln(f̂₋) pvalue: two-sided test for H₀: no discontinuity model_info['density_left']: estimated density from the left model_info['density_right']: estimated density from the right

Notes

Through 1.28.0 this was a different estimator under McCrary's name: ceil(sqrt(n/2)) bins of a width set by the narrower side, a Silverman-type bandwidth, a sandwich SE, and silent fallbacks to a density of 0.01 with SE 0.1 when a side had fewer than two bins in the window. :func:sp.rddensity (Cattaneo, Jansson & Ma) is the modern replacement; this function exists for comparability with McCrary.

References

mccrary2008manipulation

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> # Smooth running variable (no manipulation at the cutoff)
>>> df = pd.DataFrame({"score": rng.normal(0, 1, size=500)})
>>> result = sp.mccrary_test(df, x='score', c=0)
>>> bool(result.pvalue > 0.05)   # no manipulation -> fail to reject
True

oster_bounds

oster_bounds(data: Optional[DataFrame] = None, y: Optional[str] = None, treat: Optional[str] = None, controls: Optional[List[str]] = None, r_max: Optional[Union[float, str]] = None, delta: float = 1.0, beta_short: Optional[float] = None, r2_short: Optional[float] = None, beta_long: Optional[float] = None, r2_long: Optional[float] = None, alpha: float = 0.05, absorb: Optional[Union[str, List[str]]] = None, absorb_controls: Optional[Union[str, List[str]]] = None, moments: Optional[Dict[str, float]] = None, cluster: Optional[str] = None) -> OsterBoundsResult

Oster (2019) coefficient stability bounds.

Assesses robustness of a treatment effect estimate to omitted variable bias by computing how much unobservable selection (relative to observable selection) would be needed to explain away the result.

Can be called in two ways:

1. From data — runs short and long regressions internally::

   result = oster_bounds(data, y='wage', treat='training', controls=['age', 'education'])

2. From pre-computed statistics — supply β and R² directly::

   result = oster_bounds(beta_short=2.5, r2_short=0.15, beta_long=2.0, r2_long=0.45)

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data (required for method 1).	None
y	str	Outcome variable.	None
treat	str	Treatment variable.	None
controls	list of str	Control variables (included in the long but not short regression).	None
r_max	float or str	Maximum R² under full selection. Default: min(1.0, 1.3 × R²_long). Must not exceed 1; a value at or below R²_long is replaced by R²_long + 0.01 with a warning. A rule such as "1.3*r2_a" or "1.3*r2" multiplies the long regression's adjusted / plain R² (r2_a is Stata regress's, or xtreg, fe's with absorb=, the base empirical papers use with psacalc).	None
delta	float	Proportionality assumption: ratio of unobservable-to-observable selection. δ=1 means equal selection.	1.0
beta_short	float	Treatment coefficient from short regression (no controls).	None
r2_short	float	R² from short regression.	None
beta_long	float	Treatment coefficient from long regression (with controls).	None
r2_long	float	R² from long regression.	None
alpha	float	Significance level for the identified set.	0.05
absorb	str or list of str	Fixed effects in both regressions, absorbed, as after xtreg y d x, fe (the panel variable) or areg, absorb(); "id" or "id + firm". The coefficients and R² are then within, as psacalc computes them after xtreg, fe, and the exact solution is used.	None
absorb_controls	str or list of str	Fixed effects that are controls: absorbed in the long regression only, like i.ind#i.year dummies in xtreg y d x i.ind#i.year, fe. "ind^year + city^year" (# also accepted).	None
moments	dict	With the summary-statistics call, the three variances psacalc uses -- {"sigma_yy": ..., "sigma_xx": ..., "t_x": ...} (variance of y, of the treatment, and of the treatment residualised on the controls) -- so the exact solution replaces the first-order approximation.	None
cluster	str	With absorb=, the cluster variable of the xtreg, fe vce(cluster c) whose e(r2_a) the "...*r2_a" rule should reproduce: when the panel is nested in the cluster, Stata no longer charges the panel means in its degrees of freedom. Point estimates and the other inputs are unaffected.	None

Returns:

Type

Description

OsterBoundsResult

A dict (the keys below, unchanged) that also has summary() / glance() / to_latex() / to_word(). Keys: - beta_short, r2_short: short regression estimates - beta_long, r2_long: long regression estimates - r_max: maximum R² used - delta_star (also delta_for_zero): δ* such that adjusted β = 0 (robustness measure) -- the key sp.oster_delta uses in model_info['delta_star'] - beta_adjusted: bias-adjusted β under (δ, R_max) - beta_adjusted_alternatives: the other real roots of Oster's quadratic (δ = 1) / cubic (δ ≠ 1) -- exact method only - r2_a_long: adjusted R² of the long regression (NaN from summaries) - inputs: the psacalc inputs used (None for the approximation) - method: "exact" when called with data -- Oster's exact solution, which reproduces her Stata psacalc (delta and beta) to machine precision -- or "approximate" when only (β, R²) summaries are given, which uses Oster's first-order approximation β* ≈ β̃ − δ (β̊ − β̃)(R_max − R̃²)/(R̃² − R̊²) - identified_set: [min(β_adjusted, β_long), max(β_adjusted, β_long)] - robust: True if identified set excludes zero - interpretation: human-readable interpretation

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> # From data
>>> rng = np.random.default_rng(1)
>>> n = 300
>>> age = rng.normal(40, 10, n)
>>> education = rng.normal(13, 3, n)
>>> training = (rng.normal(size=n) + 0.3 * (education - 13) > 0).astype(float)
>>> wage = 20 + 2.0 * training + 0.1 * age + 0.5 * education + rng.normal(size=n)
>>> df = pd.DataFrame({"wage": wage, "training": training,
...                    "age": age, "education": education})
>>> result = sp.oster_bounds(df, y='wage', treat='training',
...                          controls=['age', 'education'])
>>> bool('delta_for_zero' in result)
True

>>> # From statistics (e.g., read from a paper)
>>> result = sp.oster_bounds(beta_short=2.5, r2_short=0.15,
...                          beta_long=2.0, r2_long=0.45)
>>> round(result['delta_for_zero'], 2)
8.89

diagnose

diagnose(data: DataFrame, y: str, x: List[str], print_results: bool = True) -> Dict[str, Any]

Comprehensive regression diagnostics in one call.

The tests are those of R's lmtest defaults; see :func:het_test and :func:reset_test for the Stata options that reproduce them.

Parameters:

Name	Type	Description	Default
data	DataFrame	Data used in the regression.	required
y	str	Dependent variable name.	required
x	list of str	Independent variable names (excluding constant).	required
print_results	bool	Print formatted output.	True

Returns:

Type	Description
dict	Keys: 'het_test', 'reset_test', 'vif', each containing test results.

Examples:

Run the full diagnostic suite on a Mincer-style wage regression:

>>> import statspai as sp
>>> df = sp.cps_wage()
>>> diag = sp.diagnose(df, y='log_wage',
...                    x=['education', 'experience'],
...                    print_results=False)
>>> sorted(diag.keys())
['het_test', 'reset_test', 'vif']

References

breusch1979simple, ramsey1969tests, white1980heteroskedasticity

het_test

het_test(data: DataFrame, y: str, x: List[str]) -> Dict[str, float]

Breusch-Pagan test for heteroskedasticity.

H0: Homoskedastic errors (constant variance). H1: Variance depends on regressors.

This is Koenker's studentized form, N R^2 from the regression of the squared residuals on every regressor: R's lmtest::bptest default and Stata's estat hettest, rhs iid. Stata's plain estat hettest uses the fitted values and assumes normal errors; sp.estat(fit, 'hettest', variables='fitted', version='normal') gives that one.

Parameters:

Name	Type	Description	Default
data	as in ``diagnose()``		required
y	as in ``diagnose()``		required
x	as in ``diagnose()``		required

Returns:

Type	Description
dict	'statistic', 'df', 'pvalue'.

References

Breusch, T.S. and Pagan, A.R. (1979). Econometrica, 47(5), 1287-1294. [@breusch1979simple]

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> df = pd.DataFrame({
...     "x1": rng.normal(size=200),
...     "x2": rng.normal(size=200),
... })
>>> df["y"] = (1.0 + 0.5 * df["x1"] - 0.3 * df["x2"]
...            + rng.normal(size=200))
>>> ht = sp.het_test(df, y="y", x=["x1", "x2"])
>>> int(ht["df"])
2
>>> ht["pvalue"] < 0.05   # homoskedastic data: do not reject H0
False

reset_test

reset_test(data: DataFrame, y: str, x: List[str], powers: int = 3) -> Dict[str, float]

Ramsey RESET test for functional form misspecification.

Tests whether nonlinear combinations of fitted values help explain the dependent variable. Rejection suggests the model is missing important nonlinearities.

With the default powers=3 this is R's lmtest::resettest; Stata's estat ovtest adds powers 2 to 4; powers=4 agrees with it to about 1e-6 (the two rescale the fitted values differently).

Parameters:

Name	Type	Description	Default
data	as in ``diagnose()``		required
y	as in ``diagnose()``		required
x	as in ``diagnose()``		required
powers	int	Include ŷ², ŷ³, ..., ŷ^powers in the auxiliary regression.	3

Returns:

Type	Description
dict	'statistic' (F), 'df1', 'df2', 'pvalue'.

References

Ramsey, J.B. (1969). JRSS-B, 31(2), 350-371. [@ramsey1969tests]

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> df = pd.DataFrame({
...     "x1": rng.normal(size=200),
...     "x2": rng.normal(size=200),
... })
>>> df["y"] = (1.0 + 0.5 * df["x1"] - 0.3 * df["x2"]
...            + rng.normal(size=200))
>>> rt = sp.reset_test(df, y="y", x=["x1", "x2"])
>>> rt["df1"], rt["df2"]
(2, 195)
>>> rt["pvalue"] > 0.05   # linear form is correct: do not reject
True

vif

vif(data: DataFrame, x: List[str]) -> DataFrame

Variance Inflation Factors for multicollinearity detection.

VIF > 10 is a common rule of thumb for concerning multicollinearity.

Equivalent to Stata's estat vif.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
x	list of str	Independent variables.	required

Returns:

Type	Description
DataFrame	Columns: 'variable', 'VIF', '1/VIF'.

Notes

VIF_j = 1 / (1 - R²_j) where R²_j is from regressing x_j on all other x variables.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> df = pd.DataFrame({
...     "x1": rng.normal(size=200),
...     "x2": rng.normal(size=200),
... })
>>> v = sp.vif(df, x=["x1", "x2"])
>>> list(v.columns)
['variable', 'VIF', '1/VIF']
>>> bool((v["VIF"] < 5).all())   # near-orthogonal: low VIFs
True

anderson_rubin_test

anderson_rubin_test(data: DataFrame, y: str, endog: Union[str, List[str]], instruments: List[str], exog: Optional[List[str]] = None, h0: Union[float, List[float]] = 0, alpha: float = 0.05, vcov: str = 'HC1', absorb: Optional[Union[str, List[str]]] = None, cluster: Optional[Union[str, List[str]]] = None, ar_vcov: str = 'classic') -> Dict[str, Any]

Anderson-Rubin (1949) test — size-correct under weak instruments.

Tests H0: β_endog = h0 and constructs a confidence set by inverting the test over a grid of candidate values. The AR test has correct size regardless of instrument strength and is the recommended inference procedure when F_eff < 10.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable.	required
endog	str or list of str	Endogenous regressor(s). With a list of two or more, the test is the joint AR F test of beta = h0 for all of them (at h0 = 0 Stata ivreg2's e(arf)), in the same ar_vcov / cluster / absorb forms; no confidence set, effective F or tF value is returned then (the AR set is a region in R^p -- use sp.weakrobust(..., project=[...]) for intervals on single coefficients).	required
instruments	list of str	Excluded instruments.	required
exog	list of str	Included exogenous controls.	None
h0	float or list of float	Null hypothesis value for the endogenous coefficient(s): one per endogenous regressor, or one number used for all.	0
alpha	float	Significance level.	0.05
vcov	(HC0, HC1, classic)	Variance estimator used for the Olea-Pflueger effective F reported alongside AR (one endogenous regressor).	'HC0'
absorb	str or list of str	High-dimensional fixed effects to partial out of y, the endogenous regressor, the instruments and the controls before the test — the same residualisation sp.iv(absorb=...) performs.	None
cluster	str or list of str	Cluster the AR statistic. Switches from the homoskedastic F-form to the (multiway) cluster-robust quadratic form (Z'e)' Omega^-1 (Z'e) with Omega the cluster-sum variance of the moment vector, referred to F(k_z, G-1). Without this, a panel AR test over-rejects exactly as a homoskedastic t-test would: the AR set is not automatically robust to serial correlation just because it is robust to weak instruments.	None
ar_vcov	(classic, HC0, HC1, HC2, HC3)	Variance behind the AR statistic itself when cluster is not given. 'classic' is the homoskedastic F of the regression of Y - b0 * D on the instruments (ivmodel::AR.test). The 'HC*' choices make the test robust to heteroskedasticity: the statistic is the Wald test that the instruments have zero coefficients in that regression, with the Eicker-Huber-White covariance of the chosen type, divided by k_z and referred to F(k_z, n - k_w - k_z). With one instrument it is the square of the robust t-ratio of the instrument, the Fieller-Anderson-Rubin interval of randomized encouragement designs. Under heteroskedasticity the classic statistic is not size-correct even with strong instruments. 'HC2' / 'HC3' need the leverages and are not available with absorb.	'classic'

Returns:

Type

Description

dict

'ar_stat' : AR F statistic at h0. 'ar_df' : Degrees of freedom (k_z, n - k_w - k_z). 'ar_pvalue' : P-value of AR test at h0. 'ar_ci' : (low, high) AR confidence set (grid-inverted). 'ar_vcov' : variance behind the AR statistic ('classic', an 'HC*' type, or 'cluster'). 'beta_2sls' : 2SLS point estimate. 'first_stage_F' : Classical first-stage F. 'effective_F' : Olea-Pflueger robust F_eff. 'tF_critical_value' : Lee et al. (2022) adjusted 5 % two-sided critical value at the effective F. None if alpha != 0.05 or there is more than one instrument. 'strength' : Instrument-strength interpretation. 'interpretation' : Human-readable summary.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> z1 = rng.normal(size=n)
>>> z2 = rng.normal(size=n)
>>> u = rng.normal(size=n)
>>> educ = 0.7 * z1 + 0.5 * z2 + 0.5 * u + rng.normal(size=n)
>>> df = pd.DataFrame({
...     "wage": 1.0 + 0.4 * educ + u + rng.normal(size=n),
...     "educ": educ, "z1": z1, "z2": z2,
... })
>>> result = sp.anderson_rubin_test(df, y='wage', endog='educ',
...                                 instruments=['z1', 'z2'])
>>> bool(0.0 <= result['ar_pvalue'] <= 1.0)
True

Notes

Under H0: β = β₀, regress Y - β₀·D on Z and W. The F-test on Z gives the AR statistic. The confidence set is constructed by collecting all β₀ values not rejected at level alpha. If the AR set is unbounded on one or both sides, the returned (low, high) will be ±inf accordingly.

effective_f_test

effective_f_test(data: DataFrame, endog: str, instruments: List[str], exog: Optional[List[str]] = None, vcov: str = 'HC1', absorb: Optional[Union[str, List[str]]] = None, cluster: Optional[Union[str, List[str]]] = None, *, y: Optional[str] = None, alpha: float = 0.05) -> Dict[str, Any]

Olea-Pflueger (2013) robust effective F statistic for weak instruments.

Computes the heteroskedasticity-robust effective F that is a pre-test for the concentration parameter of the first stage. Under homoskedasticity (vcov='classic') and a single instrument it reduces exactly to the standard first-stage F.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
endog	str	Endogenous regressor (single endogenous variable).	required
instruments	list of str	Excluded instruments.	required
exog	list of str	Included exogenous controls (a constant is added automatically).	None
vcov	(HC0, HC1, classic)	Variance estimator for the first-stage residuals: 'classic' — homoskedastic; F_eff equals first-stage F. 'HC0' — White heteroskedasticity-robust. 'HC1' — HC0 with small-sample correction n/(n-k). Ignored when cluster is given.	'HC0'
absorb	str or list of str	High-dimensional fixed effects to partial out of the endogenous regressor, the instruments and the controls before the first stage — the same residualisation sp.iv(absorb=...) performs, so the effective F describes the specification actually fitted. The absorbed degrees of freedom are charged to df_resid.	None
cluster	str or list of str	Cluster the first-stage moment variance. Omega becomes the (multiway, Cameron-Gelbach-Miller) cluster-sum meat with the ivreg2 finite-sample factor G_min/(G_min-1) * (n-1)/(n-K). This is the right diagnostic whenever the second stage is clustered: heteroskedasticity-only F_eff routinely overstates instrument strength under within-cluster correlation.	None
y	str	The outcome. Needed for the estimator-specific critical values: the worst-case bias of TSLS and of LIML depends on the reduced form as well as the first stage. Without it only the 'simplified' critical values are returned.	None
alpha	float	Level of the test of the null that the worst-case bias exceeds a fraction tau of the benchmark.	0.05

Returns:

Type

Description

dict

'F_eff' : Olea-Pflueger effective F. 'first_stage_F' : Classical first-stage F (for comparison). 'n_instruments' : Number of excluded instruments (k_z). 'n_obs' : Sample size. 'strength' : Interpretation string. 'critical_values' : {'simplified': {tau: c}, 'tsls': {...}, 'liml': {...}} for tau in 5%, 10%, 20% and 30%. Reject weak instruments (bias above tau of the benchmark) when F_eff exceeds c. 'simplified' uses the bound that holds for any estimator and data; 'tsls' and 'liml' use the worst-case Nagar bias computed from the data and are smaller. They are what Stata weakivtest prints, and agree with it to seven digits at 5%, 10% and 20%. At 30% weakivtest uses x = 3.33 where the definition gives 1 / 0.3; the exact value is used here. 'effective_df' : the matching K_eff, same layout. 'worst_case_bias' : {'tsls': B, 'liml': B} (with y). 'stock_yogo_10pct' : 23.1. Despite the key, kept for compatibility, this is the Montiel Olea-Pflueger rule of thumb for tau = 10% (the largest 'simplified' critical value over designs), not a Stock-Yogo value.

Notes

The formula (Andrews-Stock-Sun 2019, eq. 4.13) is

.. math::

F_{\text{eff}} = \frac{\hat\pi' (\tilde Z' \tilde Z)\hat\pi}
    {\mathrm{tr}\!\left(\hat\Omega\,(\tilde Z' \tilde Z)^{-1}\right)}

where :math:\tilde Z, \tilde D are residualized after partialling out the exogenous controls, :math:\hat\pi is the first-stage OLS coefficient vector, and :math:\hat\Omega = \sum_i \hat\eta_i^2 \tilde z_i \tilde z_i' is the HC meat.

Under homoskedasticity :math:\hat\Omega \approx \hat\sigma_\eta^2 \tilde Z'\tilde Z, so the trace collapses to :math:k_z \hat \sigma_\eta^2 and :math:F_{\text{eff}} reduces to the first-stage F.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> z1 = rng.normal(size=n)
>>> z2 = rng.normal(size=n)
>>> u = rng.normal(size=n)
>>> educ = 0.7 * z1 + 0.5 * z2 + 0.5 * u + rng.normal(size=n)
>>> df = pd.DataFrame({
...     "wage": 1.0 + 0.4 * educ + u + rng.normal(size=n),
...     "educ": educ, "z1": z1, "z2": z2,
... })
>>> res = sp.effective_f_test(df, endog='educ', instruments=['z1', 'z2'])
>>> bool(res['F_eff'] > res['stock_yogo_10pct'])  # strong instruments
True

tF_critical_value

tF_critical_value(first_stage_F: float, alpha: float = 0.05) -> float

Lee–McCrary–Moreira–Porter (2022, AER) tF adjusted critical value.

Returns the adjusted two-sided t-ratio critical value c(F) such that |t| > c(F) is a valid 1 - alpha test of the 2SLS coefficient in the just-identified model, robust to weak instruments.

Parameters:

Name	Type	Description	Default
first_stage_F	float	First-stage F statistic of the single excluded instrument, computed with the same variance estimator as the t-ratio it adjusts (the robust / cluster Wald F -- equal to the Olea-Pflueger effective F with one instrument). This is what ivDiag passes.	required
alpha	float	Significance level. Only 0.05 is implemented (the level LMMP tabulate).	0.05

Returns:

Type	Description
float	Adjusted critical value, linearly interpolated in sqrt(F) on the LMMP table exactly as ivDiag::tF does. 1.96 once F >= 10.3**2 = 106.09. inf below F = 4, the first tabulated point: c(F) diverges as F falls to 3.84, so the tF interval is unbounded there and the Anderson-Rubin set is the inference to report. (ivDiag instead clamps every F <= 4 to c(4) = 18.66; the two agree from F = 4 upward.)

Raises:

Type	Description
ValueError	If alpha is not 0.05.

Examples:

>>> import statspai as sp
>>> round(sp.tF_critical_value(10.0), 2)  # |t| must exceed 3.44 at F = 10
3.44
>>> sp.tF_critical_value(200.0)
1.96

References

lee2022valid, lal2024much

weakrobust

weakrobust(data: DataFrame, y: str, endog: Union[str, List[str]], instruments: List[str], exog: Optional[List[str]] = None, *, h0: Union[float, List[float]] = 0.0, alpha: float = 0.05, vcov: str = 'HC1', include_clr: bool = True, include_k: bool = True, clr_simulations: int = 20000, grid_size: int = 401, random_state: Optional[int] = None, clr_method: str = 'exact', cluster: Optional[Union[str, List[str]]] = None, test_vcov: Optional[str] = None, subset: Optional[Union[str, List[str]]] = None, project: Optional[Union[str, List[str]]] = None) -> WeakRobustResult

Stata-style unified weak-instrument-robust diagnostic panel.

Bundles in a single call:

1. Classical first-stage F and Olea–Pflueger (2013) effective F — instrument-strength pre-tests.
2. Kleibergen–Paap (2006) rk LM and Wald F — rank test of the reduced form, heteroskedasticity-robust Cragg-Donald analogue.
3. Anderson–Rubin (1949) F test of H0: β_endog = h0 and its weak-IV-robust confidence set.
4. Moreira (2003) Conditional LR (CLR) test and CLR confidence set by grid inversion (uniformly most powerful invariant for a single endogenous regressor).
5. Kleibergen (2002) K score test and K confidence set.
6. Lee–McCrary–Moreira–Porter (2022) tF adjusted critical value.

This entry point is the Python analogue of Stata 19's ivregress 2sls; estat weakrobust and weakiv / ivreg2 in Stata 17+, unifying tools that are scattered across ivmodel (R), linearmodels (Python), and the user-written packages weakiv, rivtest (Stata).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome.	required
endog	str or list of str	Endogenous regressor(s). With two or more the tests are joint tests of beta = h0 in Stata weakiv's form (see Notes); subset= / project= give intervals for single coefficients.	required
instruments	list of str	Excluded instruments.	required
exog	list of str	Included exogenous controls. An intercept is always added.	None
h0	float or list of float	Null value(s) for the endogenous coefficient(s): one per endogenous regressor, or one number used for all. All under-H0 tests (AR, CLR, K) are evaluated at beta = h0.	0.0
alpha	float	Significance level for the robust confidence sets.	0.05
vcov	(HC0, HC1, classic, cluster)	Used by the Olea–Pflueger effective F. 'cluster' (with cluster=) also makes the weak-IV-robust tests cluster-robust, i.e. it implies test_vcov='cluster'.	'HC0'
include_clr	bool	Also run the CLR test and invert it for a CLR confidence set.	True
include_k	bool	Also run the Kleibergen K score test and K confidence set.	True
clr_simulations	int	Monte-Carlo draws for the CLR null distribution; used only with clr_method='simulate'.	20 000
grid_size	int	Grid resolution used by AR/CLR/K confidence-set inversion.	401
random_state	int		None
clr_method	(exact, simulate)	How the CLR conditional null distribution is evaluated. 'exact' integrates it numerically (the closed form R ivmodel::CLR and Stata weakiv use), making the CLR p-value and set deterministic. 'simulate' is the pre-1.29 Monte-Carlo version. Used by the 'ivmodel' convention only.	'exact'
cluster	str	Cluster variable for the weak-IV-robust tests (one-way, as Stata weakiv ..., cluster()); implies test_vcov='cluster'. The effective F is clustered on it too.	None
test_vcov	(classic, robust, cluster)	Variance behind the AR / K / J / CLR tests, in Stata weakiv's form (see Notes). 'robust' is heteroskedasticity-robust (weakiv ..., robust). Default: 'cluster' when cluster is given, else 'classic'; leaving it unset with one endogenous regressor and no cluster keeps the ivmodel convention.	'classic'
subset	str or list of str	With several endogenous regressors: names whose coefficient gets a subset AR test (weakiv ..., subset()) at its h0 value and a confidence set, the other coefficients concentrated out by LIML. Referred to chi2(k - p + 1). i.i.d. errors only, as in Stata.	None
project	str or list of str	With several endogenous regressors: names whose coefficient gets the projection of the joint AR set (weakiv ..., project()): b1 is in the set if some value of the other coefficients is not rejected by the joint AR test, chi2(k). Conservative. Computed in closed form (LIML root) over the whole real line, where Stata searches a user-bounded grid. i.i.d. errors only.	None

Returns:

Type	Description
WeakRobustResult	Accepts .summary(), .to_frame() and dict-style lookup.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> z1 = rng.normal(size=n)
>>> z2 = rng.normal(size=n)
>>> u = rng.normal(size=n)
>>> educ = 0.7 * z1 + 0.5 * z2 + 0.5 * u + rng.normal(size=n)
>>> df = pd.DataFrame({
...     "wage": 1.0 + 0.4 * educ + u + rng.normal(size=n),
...     "educ": educ, "z1": z1, "z2": z2,
... })
>>> panel = sp.weakrobust(df, y='wage', endog='educ',
...                       instruments=['z1', 'z2'],
...                       clr_simulations=2000, grid_size=101,
...                       random_state=0)
>>> bool("effective_F" in panel.to_dict())
True

Notes

Two conventions. With one endogenous regressor and none of cluster / test_vcov / subset / project, the panel keeps the ivmodel convention (result['convention'] == 'ivmodel'): AR is an F statistic with n - k_w - k denominator degrees of freedom, CLR / K as R ivmodel. Otherwise ('weakiv') every statistic is in Stata weakiv's LM form -- AR chi2(k), K chi2(p), J chi2(k-p), the K-J combination (kwt = 0.8) and CLR with its conditioning rank statistic rk_stat -- and they reproduce weakiv to 1e-10 under i.i.d., robust and cluster variances. The CLR p-value is the exact conditional one: weakiv approximates it with a 20-node Simpson rule for one regressor (differences ~1e-7), and for several simulates it with AR and K drawn independently although AR = K + J, which is a different distribution (p 0.80 where the exact value is 0.75). The sets for one regressor invert the reported LM tests exactly; Stata's closed-form AR and K intervals invert the Wald form instead, so its own AR p-value at its interval endpoint is 0.053, not 0.05.

References

Anderson–Rubin (1949) Ann. Math. Stat. 20, 46-63. Kleibergen (2002) Econometrica 70, 1781-1803. Moreira (2003) Econometrica 71, 1027-1048. Kleibergen–Paap (2006) J. Econom. 133, 97-126. Olea–Pflueger (2013) JBES 31, 358-369. Lee–McCrary–Moreira–Porter (2022) AER 112, 3260-3290. [@anderson1949estimation]

weighted_rank_power

weighted_rank_power(y: Any, gammas: Sequence[float], *, data: Optional[DataFrame] = None, treat: Optional[str] = None, block: Optional[str] = None, treated: Any = 1, phi: PhiSpec = 'u868', sample_ratio: float = 1.0, alpha: float = 0.05) -> WeightedRankPowerResult

Estimated power of a weighted rank sensitivity analysis.

Uses the data in hand as a pilot. The sensitivity analysis at Gamma rejects when the statistic exceeds its bounding expectation by z_alpha bounding standard deviations; how often that happens depends on where the statistic is centred and how variable it is, both of which the jackknife over blocks estimates. Power is then read from the normal approximation, for a study with sample_ratio times as many blocks.

The estimate answers a design question (which weights would have served this study best, and how many blocks a similar study needs) and is not a substitute for the analysis itself: choosing phi by its estimated power on the same outcomes and then reporting that test ignores the selection. Use a list of phi in :func:weighted_rank for a choice that is accounted for.

Parameters:

Name	Type	Description	Default
y	Any	The block design, as in :func:weighted_rank.	required
data	Any	The block design, as in :func:weighted_rank.	required
treat	Any	The block design, as in :func:weighted_rank.	required
block	Any	The block design, as in :func:weighted_rank.	required
treated	Any	The block design, as in :func:weighted_rank.	required
gammas	sequence of float	Values of the sensitivity parameter.	required
phi	(str, tuple or callable)	Block weights, as in :func:weighted_rank.	"u868"
sample_ratio	float	Number of blocks in the planned study divided by the number observed.	1.0
alpha	float	Level of the one-sided test.	0.05

Returns:

Type	Description
WeightedRankPowerResult

Notes

With sample_ratio = 1 the result equals R's weightedRank::estPower. For other ratios that function divides the bounding standard deviation by the ratio where its square root is wanted (the variance of a mean over s times as many blocks is 1/s of the original), which overstates power when s > 1; the computation here follows the variance.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = rng.normal(size=(150, 3))
>>> y[:, 0] += 1.0
>>> small = sp.weighted_rank_power(y, gammas=[3.0])
>>> large = sp.weighted_rank_power(y, gammas=[3.0], sample_ratio=4.0)
>>> bool(large.power["power"][0] > small.power["power"][0])
True

References

[@rosenbaum2023bahadur]