statspai.inference

inference

Inference module for StatsPAI.

Provides robust inference methods for supported estimator results: - Wild Cluster Bootstrap (Cameron, Gelbach & Miller 2008) - Randomization Inference (Fisher 1935; Young 2019)

BiTestResult

Bases: ResultProtocolMixin

Outcome of :func:bitest.

Attributes:

Name	Type	Description
n_obs, successes	int	Number of trials and of successes.
p_null	float	Probability of success under the null.
estimate	float	Observed share of successes.
expected	float	n_obs * p_null.
pvalue	float	Two-sided p-value. pvalue_greater is Pr(k >= successes) and pvalue_less is Pr(k <= successes) under the null.
k_opposite	int or None	The outcome in the opposite tail at which the two-sided region starts (Stata's r(k_opp)); None when that tail is empty.
ci	tuple of float	Clopper-Pearson 1 - alpha interval for the share.

Examples:

>>> import statspai as sp
>>> res = sp.bitest(n=41, successes=25)
>>> round(res.pvalue, 4)
0.211

BootstrapResult dataclass

Bases: ResultProtocolMixin

Container for bootstrap inference results.

Returned by :func:bootstrap. Holds the point estimate, the bootstrap standard error, the confidence interval, the two-sided p-value, and the full bootstrap distribution.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"y": rng.normal(1.0, 1.0, 200),
...                    "x": rng.normal(0.0, 1.0, 200)})
>>> res = sp.bootstrap(df, lambda d: d["y"].mean() - d["x"].mean(),
...                    n_boot=200, seed=0)
>>> type(res).__name__
'BootstrapResult'
>>> res.n_boot
200
>>> res.ci_method
'percentile'
>>> bool(res.ci_lower < res.estimate < res.ci_upper)
True

MetaAnalysisResult dataclass

Bases: ResultProtocolMixin

Result of a summary-data meta-analysis.

Attributes:

Name	Type	Description
estimate	float	Pooled effect under the chosen model (random-effects by default).
se, ci	(float, tuple)	Standard error and confidence interval of the pooled effect.
p_value	float	Two-sided p-value for the pooled effect being zero.
method	str	"fixed" or "DL" (DerSimonian-Laird random effects).
fixed_estimate, fixed_se	float	Fixed-effect pooled estimate and SE (always reported).
random_estimate, random_se	float	Random-effects pooled estimate and SE (always reported).
tau2, q, q_df, q_pvalue, i2, h2	float	Between-study variance and heterogeneity statistics.
prediction_interval	tuple or None	Random-effects prediction interval (where a future study's true effect is expected to lie); None when fewer than 3 studies.
weights	ndarray	Per-study weights under the chosen model (normalised to sum to 1).

Examples:

>>> import statspai as sp
>>> effects = [0.30, 0.45, 0.20, 0.55, 0.35]
>>> se = [0.12, 0.18, 0.10, 0.22, 0.14]
>>> result = sp.meta_analysis(effects, se, method="DL")
>>> type(result).__name__
'MetaAnalysisResult'
>>> bool(0.0 <= result.i2 <= 100.0)  # I^2 heterogeneity is a percentage
True
>>> bool(result.estimate > 0)  # all studies show a positive effect
True

egger_test

egger_test() -> Dict[str, float]

Egger's regression test for funnel-plot asymmetry.

Regresses the standard normal deviate y_i / se_i on precision 1 / se_i; a non-zero intercept indicates small-study effects. Requires at least 3 studies.

forest_plot

forest_plot(ax: Any = None, **kwargs: Any) -> Any

Forest plot: per-study effects + CIs with the pooled diamond.

funnel_plot

funnel_plot(ax: Any = None, **kwargs: Any) -> Any

Funnel plot: effect size vs. standard error (precision).

PATEEstimator

Population Average Treatment Effect estimator.

Reweights an experimental (RCT) sample to a target population to recover the PATE, correcting the SATE for external validity. The convenience wrapper :func:pate builds this object and calls :meth:fit; construct it directly when you want to reuse the estimator.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n_exp = 300
>>> x_exp = rng.normal(0.0, 1.0, n_exp)
>>> d = rng.integers(0, 2, n_exp)
>>> y = 1.0 + 0.5 * x_exp + 2.0 * d + rng.normal(0, 1.0, n_exp)
>>> df_rct = pd.DataFrame({"x": x_exp, "treated": d, "outcome": y})
>>> df_pop = pd.DataFrame({"x": rng.normal(0.5, 1.0, 400)})  # target pop
>>> est = sp.PATEEstimator(
...     data_experiment=df_rct,
...     data_target=df_pop,
...     y="outcome",
...     treatment="treated",
...     covariates=["x"],
...     method="ipw",
...     n_boot=50,
...     seed=0,
... )
>>> result = est.fit()
>>> bool(np.isfinite(result.estimate))
True

fit

fit() -> CausalResult

Estimate PATE and return a CausalResult.

FisherResult

Bases: ResultProtocolMixin

Result container for Fisher's exact permutation test.

Attributes:

Name	Type	Description
statistic	float	Observed test statistic.
p_value	float	Two-sided permutation p-value.
p_one_sided	float	One-sided (greater) permutation p-value.
ci	tuple of float	For statistic='ate' the confidence interval for a constant effect, by inversion of the test. For the other statistics the central 1 - alpha range of the null distribution of the statistic, which is not an interval for the effect.
perm_dist	ndarray	Full permutation distribution of the test statistic.
statistic_type	str	Name of the test statistic used.
n_perm	int	Number of permutations performed.
n_obs	int	Number of observations.
n_treated	int	Number of treated units.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> d = rng.integers(0, 2, 80)
>>> y = 0.5 * d + rng.normal(size=80)
>>> df = pd.DataFrame({"outcome": y, "treated": d})
>>> res = sp.fisher_exact(data=df, y="outcome", treatment="treated",
...                       statistic="ate", n_perm=2000, seed=42)
>>> type(res).__name__
'FisherResult'
>>> bool(0.0 <= res.p_value <= 1.0)
True
>>> res.n_obs
80

summary

summary() -> str

Return a formatted summary string.

plot

plot(ax: Optional[Any] = None, figsize: Tuple[int, int] = (8, 5)) -> Tuple[Any, Any]

Plot the permutation distribution with observed value.

Parameters:

Name	Type	Description	Default
ax	matplotlib Axes	Axes to plot on. If None, creates a new figure.	None
figsize	tuple	Figure size if creating a new figure.	(8, 5)

Returns:

Type	Description
tuple of (fig, ax)

ClassicTestResult dataclass

Bases: ResultProtocolMixin

Result of one of the tests of this module.

Attributes:

Name	Type	Description
method	str	The test.
statistic	float	The test statistic named in statistic_name.
pvalue	float	Its p-value (two-sided where the test has sides).
df	(float or tuple, optional)	Degrees of freedom of the reference distribution.
n_obs	int	Rows used.
estimates	dict	The other numbers the test reports (rank sums, variances, the correlation, ...), by name.
table	(DataFrame, optional)	The table by group, where the test has one.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"y": [1.0, 2, 3, 4, 5, 6, 7, 8],
...                    "g": [0, 0, 0, 0, 1, 1, 1, 1]})
>>> res = sp.ranksum(df, "y", by="g")
>>> type(res).__name__, res.statistic_name, res.n_obs
('ClassicTestResult', 'z', 8)
>>> sorted(res.estimates)[:2]
['adjusted_variance', 'porder']

SDTestResult

Bases: ResultProtocolMixin

Outcome of :func:sdtest.

Attributes:

Name	Type	Description
estimate	float	The sample standard deviation (one sample) or the ratio of the two sample standard deviations.
statistic	float	(n - 1) s^2 / sd0^2 (chi-squared) for one sample, s1^2 / s2^2 (F) for two.
df	float or tuple of float	Degrees of freedom: n - 1, or (n1 - 1, n2 - 1).
pvalue, pvalue_less, pvalue_greater	float	Two-sided p-value (twice the smaller tail) and the one-sided ones (Ha: sd < sd0 / Ha: sd > sd0, or the same for the ratio against one).
ci	tuple of float	1 - alpha confidence interval for the standard deviation (one sample) or for the ratio of standard deviations (two samples).
groups	DataFrame	One row per sample with n, mean, se, sd, ci_lower and ci_upper (the interval is for the mean).

Examples:

>>> import statspai as sp
>>> res = sp.sdtest(n=10, sd=1.3 ** 0.5, sd0=2)
>>> round(res.statistic, 3), res.df
(2.925, 9.0)
>>> bool(res.pvalue_less < 0.05)
True

SuestResult dataclass

Bases: ResultProtocolMixin

Joint estimates of several OLS equations.

Attributes:

Name	Type	Description
params	Series	Coefficients, indexed "<equation>:<term>" (the equation name is its dependent variable).
vcov	DataFrame	Joint covariance of all coefficients.
n_obs	dict	Estimation sample size of each equation.
n_clusters	int or None

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> d = pd.DataFrame({"x": rng.normal(size=100)})
>>> d["y1"], d["y2"] = d.x + rng.normal(size=100), rng.normal(size=100)
>>> r = sp.suest(d, ["y1 ~ x", "y2 ~ x"])
>>> type(r).__name__, r.test_zero("x")["df"]
('SuestResult', 2)

wald

wald(R: Any, q: Optional[Sequence[float]] = None) -> Dict[str, float]

Wald test of R b = q (chi-squared, as Stata test).

R is a matrix over params (columns in params order) or a list of dicts mapping coefficient names to weights, one per restriction.

test_equal

test_equal(term: str, equations: Optional[Sequence[str]] = None) -> Dict[str, float]

Test that term has the same coefficient in every equation.

test_zero

test_zero(term: str, equations: Optional[Sequence[str]] = None) -> Dict[str, float]

Test that term is zero in every equation (joint significance).

TTestResult

Bases: ResultProtocolMixin

Outcome of :func:ttest.

Attributes:

Name	Type	Description
estimate	float	The sample mean (one sample), the mean of the differences (paired), or the difference between the two group means.
se, statistic, df	float	Standard error of estimate, the t statistic for H0: estimate == null and its degrees of freedom.
pvalue	float	Two-sided p-value. pvalue_less and pvalue_greater are the one-sided ones (Ha: diff < null and Ha: diff > null).
ci	tuple of float	1 - alpha confidence interval for estimate.
groups	DataFrame	One row per sample (and combined for two unpaired samples) with n, mean, se, sd, ci_lower and ci_upper.
method	str	Which test was run.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"g": np.repeat([0, 1], 50)})
>>> df["y"] = 1.0 + 0.5 * df.g + rng.normal(size=100)
>>> res = sp.ttest(df, "y", by="g", unequal=True)
>>> type(res).__name__
'TTestResult'
>>> list(res.groups.index)
['0', '1', 'combined']
>>> bool(res.ci[0] < res.estimate < res.ci[1])
True

ci

ci(data: DataFrame, variables: _Vars = None, *, stat: str = 'means', method: str = 'exact', alpha: float = 0.05) -> DataFrame

Confidence intervals for means, variances or proportions.

Equivalent to Stata's ci means, ci variances (sd for the standard deviation) and ci proportions.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
variables	str or list of str	Columns; every numeric column when omitted. Missing values are dropped variable by variable.	None
stat	('means', 'variances', 'sd', 'proportions')	'means': mean +/- t(n - 1) se. 'variances' / 'sd': the chi-squared interval (n - 1) s^2 / chi2, exact for a normal population and sensitive to departures from it. 'proportions': for 0 / 1 variables, by method.	'means'
method	('exact', 'wald', 'wilson', 'agresti', 'jeffreys')	stat='proportions' only. 'exact' (the default, as in Stata) is the Clopper-Pearson interval, which never covers less than 1 - alpha and is conservative; 'wilson' is the score interval, usually the better choice for reporting.	'exact'
alpha	float	1 - alpha is the confidence level.	0.05

Returns:

Type	Description
DataFrame	Indexed by variable, with n, the estimate (mean, variance, sd or proportion), se (means and proportions), ci_lower and ci_upper.

Raises:

Type	Description
MethodIncompatibility	An unknown stat or method, or a proportion asked of a variable that is not 0 / 1.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"y": rng.normal(2, 1, 100),
...                    "d": rng.integers(0, 2, 100)})
>>> means = sp.ci(df, "y")
>>> bool(means.loc["y", "ci_lower"] < 2 < means.loc["y", "ci_upper"])
True
>>> list(sp.ci(df, "d", stat="proportions", method="wilson").columns)
['n', 'proportion', 'se', 'ci_lower', 'ci_upper']

prtest

prtest(data: Optional[DataFrame] = None, y: Optional[str] = None, by: Optional[str] = None, *, other: Optional[str] = None, p: float = 0.5, n: Union[None, int, Sequence[int]] = None, proportion: Union[None, float, Sequence[float]] = None, alpha: float = 0.05) -> TTestResult

Large-sample z test on a proportion, or on the difference of two.

Equivalent to Stata's prtest / prtesti and (without continuity correction) R's prop.test. For a small sample use the exact :func:bitest.

============================================ =========================== Call Stata ============================================ =========================== sp.prtest(df, "d", p=0.4) prtest d == 0.4 sp.prtest(df, "d", by="g") prtest d, by(g) sp.prtest(df, "d", other="e") prtest d == e sp.prtest(n=50, proportion=0.52, p=0.4) prtesti 50 0.52 0.4 sp.prtest(n=(30, 45),prtesti 30 0.4 45 0.67proportion=(0.4, 0.67)) ============================================ ===========================

Parameters:

Name	Type	Description	Default
data	DataFrame		None
y	str	A 0 / 1 variable.	None
by	str	Grouping column with exactly two distinct non-missing values. The estimate is the proportion in the group with the smaller value minus that in the group with the larger one.	None
other	str	A second 0 / 1 variable, treated as an independent sample.	None
p	float	One sample: the proportion under the null. Not used with two samples, whose null is equal proportions.	0.5
n	number or pair of numbers	Summary statistics in place of data: observations and sample proportion of one sample, or pairs for two.	None
proportion	number or pair of numbers	Summary statistics in place of data: observations and sample proportion of one sample, or pairs for two.	None
alpha	float	1 - alpha is the confidence level of the intervals.	0.05

Returns:

Type	Description
TTestResult	statistic is z (df is infinite). With two samples the statistic uses the standard error pooled under the null (result.se_null); se and ci use the unpooled one.

Raises:

Type	Description
MethodIncompatibility	A variable that is not 0 / 1, by without exactly two groups, or a null or sample proportion outside (0, 1).

Notes

The normal approximation is poor when n p or n (1 - p) is small; the intervals here are Wald intervals and can leave [0, 1]. :func:ci with stat='proportions' gives exact and Wilson intervals.

Examples:

>>> import statspai as sp
>>> res = sp.prtest(n=50, proportion=0.52, p=0.4)
>>> round(res.statistic, 4)
1.7321
>>> round(res.pvalue, 4)
0.0833

sktest

sktest(data: DataFrame, variables: _Vars = None, *, adjust: bool = True) -> DataFrame

Skewness and kurtosis tests for normality, and their joint test.

One row per variable: the test that the skewness is that of a normal (zero), the test that the kurtosis is (three), and the joint chi-squared with two degrees of freedom. Equivalent to Stata's sktest.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
variables	str or list of str	Columns to test; every numeric column when omitted. Missing values are dropped variable by variable.	None
adjust	bool	Apply Royston's adjustment to the joint statistic (the one Stata's sktest documents), which keeps its size in small samples. False is the D'Agostino, Belanger and D'Agostino (1990) statistic as it stands (Stata's noadjust).	True

Returns:

Type	Description
DataFrame	Indexed by variable, with n, skewness, kurtosis (the moment coefficients, 3 for a normal), p_skew, p_kurt, chi2 and p_chi2.

Raises:

Type	Description
DataInsufficient	Fewer than 8 observations of a variable (the approximations are not defined below that), or no variation.

Notes

Rejecting says the data are unlikely under normality. With a large sample the test rejects for departures too small to matter, and with a small one it has little power; a plot of the distribution is the better guide to whether a departure matters.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"z": rng.normal(size=300),
...                    "e": rng.exponential(size=300)})
>>> out = sp.sktest(df)
>>> list(out.columns)
['n', 'skewness', 'kurtosis', 'p_skew', 'p_kurt', 'chi2', 'p_chi2']
>>> bool(out.loc["e", "p_chi2"] < 0.01 < out.loc["z", "p_chi2"])
True

References

[@dagostino1990suggestion]

swilk

swilk(data: DataFrame, variables: _Vars = None) -> DataFrame

Shapiro-Wilk W test for normality.

data may also be a Series or a one-dimensional array, which is tested as a single variable.

Equivalent to Stata's swilk and R's shapiro.test: W is the squared correlation of the ordered sample with the expected normal order statistics, and the p-value comes from Royston's (1992) normal approximation.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
variables	str or list of str	Columns to test; every numeric column when omitted. Missing values are dropped variable by variable.	None

Returns:

Type	Description
DataFrame	Indexed by variable, with n, W, V (Royston's index of departure from normality, 1 at the median of its null distribution; reported for n >= 12), z and pvalue.

Raises:

Type	Description
DataInsufficient	Fewer than 4 or more than 5,000 observations (the range the approximation covers), no variation, or an infinite value.

Notes

The same caution as for :func:sktest applies: in large samples a rejection can reflect a departure of no practical importance.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"e": rng.exponential(size=200)})
>>> out = sp.swilk(df, "e")
>>> bool(out.loc["e", "pvalue"] < 0.001)
True

References

[@royston1992approximating]

cor_test

cor_test(data: DataFrame, x: str, y: str, covariates: Optional[Sequence[str]] = None, *, alpha: float = 0.05) -> ClassicTestResult

Pearson correlation, or partial correlation, with test and interval.

Without covariates this is the correlation of x and y; with them, the correlation of what is left of each after a linear regression on the covariates. t = r sqrt(df / (1 - r^2)) with df = n - 2 - k (k linearly independent, non-constant covariates; a redundant one is named in a warning and not counted) is Student's t under no correlation and joint normality, and the interval is Fisher's: tanh(atanh(r) -/+ z / sqrt(n - 3 - k)).

A partial correlation of zero is what a causal graph implies for two nodes that a set of others d-separates, under linearity: conditioning on a mediator should remove a correlation, conditioning on a collider should create one. :meth:statspai.dag.DAG.test_implications runs this test for every implication of a graph.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
x	str	Numeric columns.	required
y	str	Numeric columns.	required
covariates	sequence of str	Numeric columns to partial out. Rows with a missing value in any column used are dropped.	None
alpha	float	One minus the confidence level.	0.05

Returns:

Type	Description
ClassicTestResult	statistic is the correlation; pvalue two-sided; estimates has t, ci_lower, ci_upper and n_covariates.

Notes

Matches R's cor.test (estimate, t, df, p-value and interval) and, with covariates, the t test of ggm::pcor.test. For rank correlations see :func:statspai.spearman and :func:statspai.ktau.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> x = rng.normal(size=500)
>>> m = x + rng.normal(size=500)           # x -> m -> y
>>> y = m + rng.normal(size=500)
>>> df = pd.DataFrame({"x": x, "m": m, "y": y})
>>> bool(sp.cor_test(df, "x", "y").pvalue < 0.001)
True
>>> bool(sp.cor_test(df, "x", "y", covariates=["m"]).pvalue > 0.05)
True

esize

esize(data: DataFrame, y: str, by: str, *, unequal: bool = False, alpha: float = 0.05) -> ClassicTestResult

Effect sizes for the difference between two group means.

Reports, with confidence intervals,

• Cohen's d: the difference in means over the pooled standard deviation;
• Hedges's g: d times the exact small-sample correction Gamma(m/2) / (sqrt(m/2) Gamma((m-1)/2)), m = n1 + n2 - 2, which removes the upward bias of d;
• Glass's Delta: the difference over the standard deviation of one group alone (delta1 the first group's, delta2 the second's), for when a treatment changes the spread as well as the mean and one group is the natural yardstick;
• the point-biserial correlation between the outcome and group membership.

The difference is the mean of the group with the smaller value of by minus the mean of the other, as in :func:statspai.ttest.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome.	required
by	str	Column with exactly two groups.	required
unequal	bool	Use Satterthwaite's degrees of freedom for the intervals of d, g and the point-biserial correlation (Stata's unequal). The estimates of d, g and Delta do not change, nor do Delta's intervals.	False
alpha	float	One minus the confidence level.	0.05

Returns:

Type	Description
ClassicTestResult	statistic is Cohen's d. table has one row per effect size with estimate, ci_lower and ci_upper; estimates carries the group sizes, means and standard deviations, the t statistic and its degrees of freedom.

Notes

The intervals invert the noncentral t distribution of the t statistic, so they are exact under normality and are not symmetric about the estimate. They agree with Stata's esize twosample to rounding. R's effsize::cohen.d uses a normal approximation unless noncentral = TRUE, and corrects g by 1 - 3 / (4m - 1), which differs from the exact factor in the seventh digit at m = 400.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"g": np.repeat(["a", "b"], 100)})
>>> df["y"] = rng.normal(size=200) + 0.5 * (df.g == "b")
>>> res = sp.esize(df, "y", by="g")
>>> list(res.table.index)
['cohens_d', 'hedges_g', 'glass_delta1', 'glass_delta2', 'point_biserial_r']
>>> bool(res.table.loc["cohens_d", "ci_upper"] < 0)   # a minus b
True

iv_wild_bootstrap

iv_wild_bootstrap(result: Any, data: DataFrame, cluster: str, variable: str, n_boot: int = 999, weight_type: str = 'rademacher', seed: Optional[int] = None, alpha: float = 0.05, beta0: float = 0.0, efficient: bool = True) -> Dict[str, Any]

WRE wild cluster bootstrap for a single (endogenous) 2SLS coefficient.

Parameters:

Name	Type	Description	Default
result	EconometricResults	A fitted sp.ivreg result (2SLS); must carry data_info['iv'].	required
data	DataFrame	The estimation data (row-aligned to the fitted sample).	required
cluster	str	Cluster variable name in data.	required
variable	str	Coefficient to test (H0: beta = beta0).	required
n_boot	int	Bootstrap controls (see wild_cluster_boot).	999
weight_type	int	Bootstrap controls (see wild_cluster_boot).	999
seed	int	Bootstrap controls (see wild_cluster_boot).	999
alpha	int	Bootstrap controls (see wild_cluster_boot).	999
beta0	int	Bootstrap controls (see wild_cluster_boot).	999
efficient	bool	If True, use the DM2010 restricted efficient reduced form (regress the endogenous regressor on the instruments and the restricted structural residual, exploiting the u-v correlation). If False, use the plain restricted reduced form (projection on instruments only).	True

Returns:

Type	Description
dict	beta_hat, se_cluster, t_stat, p_boot (symmetric, strict inequality as in boottest), ci_boot (the test-inversion confidence set), t_distribution, n_clusters, n_obs, n_boot (2**G when the Rademacher grid is enumerated), enumerated, weight_type and method ('WRE' / 'WUR').

Notes

After sp.feols with absorbed fixed effects and an IV part, the stored columns are residualised on the effects and each bootstrap draw is residualised again, so the bootstrap is that of the model with the effects as dummies. sp.wild_cluster_boot routes 2SLS results here.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n, g = 400, 12
>>> firm = rng.integers(0, g, n)
>>> z = rng.normal(size=n)
>>> u = rng.normal(size=g)[firm] + rng.normal(size=n)
>>> d = 0.6 * z + 0.5 * u + rng.normal(size=n)
>>> df = pd.DataFrame({"y": 1 + 0.3 * d + u, "d": d, "z": z, "firm": firm})
>>> fit = sp.ivreg("y ~ (d ~ z)", data=df, cluster="firm")
>>> out = sp.iv_wild_bootstrap(fit, df, "firm", "d", n_boot=199, seed=1)
>>> bool(0.0 <= out["p_boot"] <= 1.0)
True

References

roodman2019fast

cr2_se

cr2_se(result: EconometricResults, data: DataFrame, cluster: str, alpha: float = 0.05) -> EconometricResults

CR2 bias-corrected cluster-robust standard errors (Bell & McCaffrey 2002).

Applies a bias correction to the cluster-robust variance estimator that accounts for leverage and improves finite-sample performance. The correction factor uses the inverse square root of (I - H_gg), where H_gg is the within-cluster hat matrix block.

Parameters:

Name	Type	Description	Default
result	EconometricResults	A fitted regression result from sp.regress().	required
data	DataFrame	The original data used for estimation.	required
cluster	str	Name of the cluster variable in data.	required
alpha	float	Significance level for confidence intervals.	0.05

Returns:

Type	Description
EconometricResults	New results object with CR2-corrected standard errors.

Notes

The CR2 estimator replaces the raw residuals e_g with adjusted residuals:

e*_g = (I - H_gg)^{-1/2} e_g

where H_gg = X_g (X'X)^{-1} X_g' is the within-cluster block of the hat matrix. This removes the downward bias in the conventional cluster-robust variance estimator.

Each coefficient's p-value and confidence interval use a t distribution with that coefficient's own Bell-McCaffrey degrees of freedom (diagnostics['satterthwaite_dof']), the ones R estimatr and clubSandwich report. data_info['df_resid'] is their minimum.

See Bell & McCaffrey (2002) and Pustejovsky & Tipton (2018).

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> state = rng.integers(0, 10, size=n)
>>> x1 = rng.normal(size=n)
>>> y = 1.0 + 0.5 * x1 + rng.normal(size=n) + rng.normal(size=10)[state]
>>> df = pd.DataFrame({"y": y, "x1": x1, "state": state})
>>> res = sp.regress("y ~ x1", data=df, cluster="state")
>>> cr2 = sp.cr2_se(res, data=df, cluster="state")
>>> cr2.model_info["se_type"]
'CR2 (Bell-McCaffrey)'

jackknife_se

jackknife_se(result: EconometricResults, data: DataFrame, cluster: str, alpha: float = 0.05) -> EconometricResults

Leave-one-cluster-out jackknife standard errors.

Computes cluster jackknife variance by re-estimating the model G times, each time dropping one cluster. Valid when G is small and conventional cluster-robust SEs are unreliable.

Parameters:

Name	Type	Description	Default
result	EconometricResults	A fitted regression result from sp.regress().	required
data	DataFrame	The original data used for estimation.	required
cluster	str	Name of the cluster variable in data.	required
alpha	float	Significance level for confidence intervals.	0.05

Returns:

Type	Description
EconometricResults	New results object with jackknife standard errors, t-statistics based on effective DoF, and adjusted p-values.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> state = rng.integers(0, 8, size=n)
>>> x1 = rng.normal(size=n)
>>> x2 = rng.normal(size=n)
>>> y = 1.0 + 0.5 * x1 - 0.3 * x2 + rng.normal(size=n) + rng.normal(size=8)[state]
>>> df = pd.DataFrame({"y": y, "x1": x1, "x2": x2, "state": state})
>>> result = sp.regress("y ~ x1 + x2", data=df, cluster="state")
>>> jk = sp.jackknife_se(result, data=df, cluster="state")
>>> jk.model_info["se_type"]
'cluster jackknife'

Notes

The cluster jackknife variance is:

V_jk = ((G-1)/G) * sum_g (theta_g - theta_bar)^2

where theta_g is the estimate from dropping cluster g, and theta_bar is the mean of the leave-one-out estimates.

The effective degrees of freedom are G - 1, which gives wider confidence intervals than asymptotic cluster-robust SEs.

Each refit is least squares on the design the fit stored (its terms, its intercept or lack of one, and sqrt(w) scaling for a weighted fit). A result whose design cannot be recovered, one with absorbed fixed effects, and a cluster variable with a single cluster are refused; sp.jackknife with a statistic that refits the model covers those estimators.

wild_cluster_boot

wild_cluster_boot(result: EconometricResults, data: DataFrame, cluster: Union[str, List[str]], variable: Optional[Union[str, List[str]]] = None, n_boot: int = 999, weight_type: str = 'rademacher', seed: Optional[int] = None, alpha: float = 0.05, h0: Union[float, List[float]] = 0.0, confidence_set: bool = False, bootcluster: Optional[Union[str, List[str]]] = None) -> Dict[str, Any]

Wild cluster bootstrap test of one coefficient, or of several jointly.

The bootstrap is the one the fitted estimator defines -- the one Stata boottest runs after the same command:

• least squares (sp.regress, sp.feols, sp.hdfe_ols, a fixed-effects sp.panel): the WCR (wild cluster restricted) bootstrap of Cameron, Gelbach & Miller (2008) on the stored design, after checking that least squares on it reproduces the coefficients;
• sp.ivreg (2SLS): the WRE bootstrap of :func:statspai.iv_wild_bootstrap (endogenous coefficients);
• sp.logit / sp.probit / sp.poisson / sp.glm (those links): the score wild cluster bootstrap (Kline-Santos), bit-exact with boottest when the Rademacher grid is enumerated;
• sp.did (2 x 2): its own vce='wild', replayed from the call record (variable may be omitted).

Any other model raises: refitting it by least squares (as this function did before) bootstraps a different estimator.

Parameters:

Name	Type	Description	Default
result	EconometricResults	A fitted result (see above).	required
data	DataFrame	The original data used for estimation.	required
cluster	str	Name of the cluster variable.	required
variable	str	Name of the coefficient to test (H0: beta = h0). Required except after sp.did, where the ATT is tested.	None
n_boot	int	Number of bootstrap replications (use odd number).	999
weight_type	str	Bootstrap weight distribution: - 'rademacher': +/-1 with equal probability. - 'webb': Webb (2014) 6-point distribution for G < 12. - 'mammen': Mammen (1993) 2-point distribution.	'rademacher'
seed	int	Random seed for reproducibility.	None
alpha	float	Significance level for confidence interval.	0.05
h0	float	The value of the coefficient under the null, imposed when the bootstrap samples are generated (Stata boottest x = h0).	0.0
confidence_set	bool	Also return ci_inverted, the set of null values the test does not reject at level alpha: the confidence set Stata boottest prints. The test is repeated over null values with the same bootstrap weights (a few hundred times), so it is off by default. One coefficient, one-way clustering.	False
bootcluster	str or list of str	Least squares: the variable(s) whose clusters the bootstrap weights are drawn for (Stata boottest, bootcluster()). Default: the intersection of the cluster variables, boottest's default. As in boottest, bootstrap clusters that cut across one-way error clusters are replaced by their intersections with them (a subcluster bootstrap, with a warning); ones nested in no error clustering and coarser are refused.	None

Notes on multi-way and joint tests

cluster=[a, b] uses the multi-way variance sum_s (-1)^(|s|+1) c_s V_s with each term's own G_s/(G_s - 1) (n - 1)/(n - k), as boottest does; draws whose variance is not positive have no statistic and are left out of the p-value (n_boot_valid). variable=[x1, x2] with h0=[r1, r2] tests the restrictions jointly by the Wald statistic F = W / q (statistic) with all of them imposed on the bootstrap data, p = #{F* > F} / B (boottest (x1 = r1) (x2 = r2)). Both are least-squares only and return no interval.

Returns:

Type

Description

dict

Keys (least squares; the other bootstraps return the subset they define and a method key naming the bootstrap): - method: which bootstrap ran - beta_hat: point estimate - se_cluster: conventional cluster-robust SE - t_stat: observed t-statistic - p_boot: symmetric two-sided bootstrap p-value, #{|t*| > |t|} / B (strict inequality, as in boottest) - ci_boot: bootstrap percentile-t confidence interval - t_distribution: array of bootstrap t-statistics - n_clusters: number of clusters - n_boot: number of replications used (2**G when the Rademacher grid is enumerated, i.e. when 2**G <= n_boot) - n_boot_requested / enumerated - ci_inverted: with confidence_set=True, the interval of null values with a bootstrap p-value of at least alpha. The p-value is a step function of the null value and each endpoint is the jump where it falls below alpha. boottest interpolates between the points of a grid instead, and its endpoints differ from these by up to about 2% of a standard error.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> state = rng.integers(0, 8, size=n)
>>> x1 = rng.normal(size=n)
>>> x2 = rng.normal(size=n)
>>> y = 1.0 + 0.5 * x1 - 0.3 * x2 + rng.normal(size=n) + rng.normal(size=8)[state]
>>> df = pd.DataFrame({"y": y, "x1": x1, "x2": x2, "state": state})
>>> result = sp.regress("y ~ x1 + x2", data=df, cluster="state")
>>> boot = sp.wild_cluster_boot(result, data=df, cluster="state",
...                             variable="x1", n_boot=199, seed=0)
>>> bool(0.0 <= boot["p_boot"] <= 1.0)
True

Notes

The WCR bootstrap imposes H0: beta_variable = 0 when constructing bootstrap samples. Under few clusters, this yields much better size control than unrestricted wild bootstrap or cluster-robust SEs.

For G < 12, use weight_type='webb' per Webb (2014).

cluster_robust_se

cluster_robust_se(X: ndarray, resid: ndarray, clusters: Union[ndarray, List[ndarray]], **kwargs: Any) -> ndarray

Return cluster-robust standard errors (diagonal sqrt of vcov).

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> firm = rng.integers(0, 10, size=n)
>>> year = rng.integers(0, 8, size=n)
>>> X = np.column_stack([np.ones(n), rng.normal(size=n)])
>>> resid = rng.normal(size=n)
>>> se = sp.cluster_robust_se(X, resid, [firm, year])
>>> se.round(3).tolist()
[0.061, 0.011]

cr3_jackknife_vcov

cr3_jackknife_vcov(X: ndarray, y: ndarray, cluster: ndarray) -> ndarray

CR3 cluster-jackknife variance (delete-one-cluster).

Drops each cluster in turn, refits OLS on the reduced sample, and forms the empirical variance of the leave-one-out coefficient vector. Scales by (G-1)/G. Often more conservative than CR1/CR2 and useful when clusters are unbalanced (Bell-McCaffrey 2002).

Complexity is O(G · k²) refits; suitable when G <= a few hundred.

Parameters:

Name	Type	Description	Default
X	(ndarray, shape(n, k))		required
y	(ndarray, shape(n))		required
cluster	(ndarray, shape(n))		required

Returns:

Type	Description
(ndarray, shape(k, k))

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> firm = rng.integers(0, 10, size=n)
>>> X = np.column_stack([np.ones(n), rng.normal(size=n)])
>>> y = X @ np.array([1.0, 0.5]) + rng.normal(size=n)
>>> V = sp.cr3_jackknife_vcov(X, y, firm)
>>> V.shape
(2, 2)

multiway_cluster_vcov

multiway_cluster_vcov(X: ndarray, resid: ndarray, clusters: Union[ndarray, List[ndarray]], df_adjust: bool = True, n_params: Optional[int] = None, psd_correct: bool = True) -> ndarray

Compute N-way cluster-robust variance of an OLS coefficient vector.

Parameters:

Name	Type	Description	Default
X	(ndarray, shape(n, k))	Design matrix used in the regression.	required
resid	(ndarray, shape(n))	OLS residuals.	required
clusters	ndarray or list of ndarrays	One or more cluster variables, one per dimension. Non-numeric labels are supported.	required
df_adjust	bool	If True, apply the G/(G-1) * (n-1)/(n-k) CR1 finite-sample correction per component variance. If False, uses raw sandwich (useful when the caller has already degreed-freedom adjusted).	True
n_params	int	Override for the k used in DOF adjustment; useful when FEs have been absorbed (pass total absorbed DOF here).	None
psd_correct	bool	Project V onto PSD cone by zeroing negative eigenvalues.	True

Returns:

Type	Description
(ndarray, shape(k, k))	Multiway-cluster-robust variance-covariance matrix.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> firm = rng.integers(0, 10, size=n)
>>> year = rng.integers(0, 8, size=n)
>>> X = np.column_stack([np.ones(n), rng.normal(size=n)])
>>> resid = rng.normal(size=n)
>>> V = sp.multiway_cluster_vcov(X, resid, [firm, year])
>>> V.shape
(2, 2)

ppi_mean

ppi_mean(*, y: Any, yhat: Any, yhat_unlabeled: Any, tune: bool = True, alpha: float = 0.05) -> CausalResult

Prediction-powered estimate of a population mean.

Combines n human-labeled outcomes with N model-predicted outcomes on unlabeled rows. Valid regardless of prediction quality; with tune=True (PPI++) the interval is asymptotically never wider than the classical labeled-only interval.

Parameters:

Name	Type	Description	Default
y	(array - like, shape(n))	Gold-standard (human) outcomes on the labeled sample.	required
yhat	(array - like, shape(n))	Model predictions on the same labeled rows.	required
yhat_unlabeled	(array - like, shape(N))	Model predictions on the unlabeled rows.	required
tune	bool	Use the PPI++ power-tuning weight λ ∈ [0, 1]. False fixes λ = 1 (the original PPI estimator).	True
alpha	float	CI level (1 - alpha confidence).	0.05

Returns:

Type	Description
CausalResult	estimate / se / ci / pvalue for the mean; model_info carries lambda, the classical labeled-only estimate and SE, the imputed-only mean, the labeled-sample correlation between y and yhat, and the effective variance ratio versus classical inference.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(7)
>>> theta = 2.0
>>> n, N = 100, 5000
>>> y_all = theta + rng.normal(0, 1, n + N)
>>> f_all = y_all + rng.normal(0, 0.5, n + N)   # decent predictions
>>> r = sp.ppi_mean(y=y_all[:n], yhat=f_all[:n],
...                 yhat_unlabeled=f_all[n:])
>>> bool(r.ci[0] < theta < r.ci[1])
True
>>> bool(r.se < np.std(y_all[:n], ddof=1) / np.sqrt(n))  # beats classical
True

References

Angelopoulos et al. (2023), Science 382(6671), 669–674 — doi:10.1126/science.adi6000. PPI++: arXiv:2311.01453.

ppi_ols

ppi_ols(*, X: Any, y: Any, yhat: Any, X_unlabeled: Any, yhat_unlabeled: Any, target: Optional[Union[str, int]] = None, tune: bool = True, add_intercept: bool = True, alpha: float = 0.05) -> CausalResult

Prediction-powered OLS with a labeled audit sample.

Runs OLS of the model-predicted outcome on the unlabeled covariates, then rectifies every coefficient with the labeled sample's OLS(y) − OLS(yhat) contrast. Standard errors combine the two independent sandwich variances; with tune=True each coefficient gets its own PPI++ weight λ_j ∈ [0, 1].

Parameters:

Name	Type	Description	Default
X	(DataFrame or array, shape(n, p))	Covariates for the labeled rows.	required
y	(array - like, shape(n))	Gold-standard outcomes for the labeled rows.	required
yhat	(array - like, shape(n))	Model-predicted outcomes for the same labeled rows.	required
X_unlabeled	(DataFrame or array, shape(N, p))	Covariates for the unlabeled rows (same columns as X).	required
yhat_unlabeled	(array - like, shape(N))	Model-predicted outcomes for the unlabeled rows.	required
target	str or int	Which coefficient is the headline estimate. Defaults to the first non-intercept term.	None
tune	bool	Per-coefficient PPI++ power tuning. False fixes λ = 1.	True
add_intercept	bool	Prepend a constant column to both design matrices.	True
alpha	float	CI level (1 - alpha confidence).	0.05

Returns:

Type	Description
CausalResult	Headline coefficient in estimate / se / ci / pvalue; the full coefficient table (with per-term λ and the classical labeled-only comparison) in detail.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(3)
>>> n, N = 150, 4000
>>> x_all = rng.normal(size=n + N)
>>> y_all = 1.0 + 2.0 * x_all + rng.normal(0, 1, n + N)
>>> f_all = y_all + rng.normal(0, 0.5, n + N)
>>> r = sp.ppi_ols(
...     X=pd.DataFrame({"x": x_all[:n]}), y=y_all[:n], yhat=f_all[:n],
...     X_unlabeled=pd.DataFrame({"x": x_all[n:]}),
...     yhat_unlabeled=f_all[n:])
>>> bool(r.ci[0] < 2.0 < r.ci[1])
True
>>> list(r.detail["term"])
['const', 'x']

References

Angelopoulos et al. (2023), Science 382(6671), 669–674 — doi:10.1126/science.adi6000. PPI++: arXiv:2311.01453.

fisher_exact

fisher_exact(data: DataFrame, y: str, treatment: str, statistic: str = 'ate', controls: Optional[List[str]] = None, n_perm: int = 10000, stratify: Optional[str] = None, cluster: Optional[str] = None, seed: Optional[int] = None, alpha: float = 0.05) -> FisherResult

Fisher's exact randomization test with enhanced features.

Computes a permutation-based p-value under the sharp null hypothesis and, for the difference in means, the confidence interval for a constant effect obtained by inverting the test.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
y	str	Outcome variable name.	required
treatment	str	Binary treatment variable name (0/1).	required
statistic	str	Test statistic to use: - 'ate': Average treatment effect (difference in means). - 'ks': Kolmogorov-Smirnov statistic. - 'rank_sum': Wilcoxon rank-sum statistic. - 't': studentized difference in means (unequal variances).	'ate'
controls	list of str	Control variables for covariate-adjusted inference. When provided, the test statistic is computed on residuals from regressing Y on controls.	None
n_perm	int	Number of random permutations. When the design admits no more than n_perm distinct assignments (complete, cluster or stratified randomization) they are all enumerated instead, so the p-value is exact (the R ri2 / randomizr rule), and n_perm on the result reports the number of assignments.	10000
stratify	str	Variable for stratified permutation (permute within strata).	None
cluster	str	Variable for cluster-level randomization (permute by cluster).	None
seed	int	Random seed for reproducibility.	None
alpha	float	Significance level for the test-inversion confidence interval.	0.05

Returns:

Type	Description
FisherResult	Object with statistic, p_value, ci, perm_dist, and methods for summary() and plot().

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> d = rng.integers(0, 2, 80)
>>> y = 0.5 * d + rng.normal(size=80)
>>> df = pd.DataFrame({"outcome": y, "treated": d})
>>> result = sp.fisher_exact(
...     data=df, y="outcome", treatment="treated",
...     statistic="ate", n_perm=2000, seed=42)
>>> type(result).__name__
'FisherResult'
>>> bool(0.0 <= result.p_value <= 1.0)
True
>>> result.n_obs
80

Notes

Under the sharp null H0: Y_i(1) = Y_i(0) for all i, treatment assignment is the only source of randomness. The p-value is the proportion of permuted statistics at least as extreme as observed.

The confidence interval inverts the permutation test: it is the set of tau_0 for which the test of Y - tau_0 * D does not reject at level alpha. It is computed on the same assignments as the p-value (the difference in means is linear in the outcome, so no new draws are needed) and each end is found by bisection; with controls the shifted outcome is residualized, as the null requires. An end is infinite, with a warning, when no hypothesised effect is rejected, which happens when the design has too few assignments for alpha. Through 1.38.0 the ends were read off a 101-point grid spanning six standard deviations of the outcome, each point tested on its own 500 draws, and with controls the unresidualized tau_0 * D was subtracted from the residualized outcome. For statistic='t', 'ks' and 'rank_sum' the interval is the same inversion, in outcome units: the statistic is recomputed on Y - tau_0 * D for a stored set of assignments. Through 1.38.0 those statistics returned the 2.5% and 97.5% percentiles of the null distribution of the statistic, which is not an interval for the effect.

For stratified experiments, treatment is permuted within strata. For cluster-randomized experiments, treatment is permuted at the cluster level; with both cluster and stratify whole clusters are permuted within their stratum (the treatment and the stratum must then be constant within clusters).

See Young (2019, QJE) for why randomization p-values should accompany asymptotic p-values in experimental papers.

ri_test

ri_test(data: DataFrame, y: str, treat: str, stat: str = 'diff_means', n_perms: int = 1000, cluster: Optional[str] = None, seed: Optional[int] = None, alpha: float = 0.05, strata: Optional[str] = None, covariates: Optional[List[str]] = None, absorb: Optional[Union[str, List[str]]] = None, interact: Optional[str] = None) -> Dict[str, Any]

Randomization inference p-value.

Computes the Fisher exact p-value by permuting the treatment vector and recalculating the test statistic.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable.	required
treat	str	Treatment. 'diff_means' / 't' / 'ks' need a binary 0/1 indicator; 'ols' / 'ols_t' and callables also take a continuous one (an exposure or dose), whose observed values are permuted. Under cluster it must be constant within clusters (the permuted unit is the cluster); a treatment that switches on within a cluster, such as exposure x post, is passed as treat='exposure', interact='post'.	required
stat	str or callable	Test statistic: - 'diff_means': difference in means (Y_bar_1 - Y_bar_0) - 'ks': Kolmogorov-Smirnov statistic - 't': the difference in means divided by its unequal-variance (Neyman) standard error. Unlike 'diff_means', the test based on it also has the right size in large samples for the weak null of a zero average effect when effects vary across units - 'rank_sum': Wilcoxon rank sum of the treated outcomes (mid-ranks), centred at its null mean n_1 (n + 1) / 2 - 'lin': coefficient on the treatment in the regression of y on the treatment, the centred covariates and their interactions with the treatment (:func:sp.lm_lin) - 'lin_t': that coefficient divided by its HC2 standard error, the covariate-adjusted statistic with the same large-sample guarantee as 't' - 'ols': coefficient on the treatment in the OLS regression of y on the treatment, covariates and (with strata) stratum fixed effects -- the regression-adjusted statistic of stratified experiments (R ri2::conduct_ri(y ~ Z + x + block)) - 'ols_t': its t-statistic (HC1, or CR1 by cluster) - A callable f(Y, D) -> float for custom statistics.	'diff_means'
n_perms	int	Number of random permutations. Use 10000+ for publications. When the design admits no more than n_perms distinct assignments (C(n, n_1), or C(G, G_1) under cluster) every assignment is enumerated instead and the p-value is exact -- the rule R ri2 / randomizr::obtain_permutation_matrix apply.	1000
cluster	str	Cluster-level permutation (permute treatment at cluster level).	None
seed	int	Random seed.	None
alpha	float		0.05
strata	str	Randomization strata (blocks): treatment is re-randomized within each stratum, keeping its number of treated units (or clusters, with cluster). Permuting across strata when the experiment assigned within them tests the wrong null distribution and loses power (UCT, QJE 2016: education p 0.70 unstratified against 0.17).	None
covariates	list of str	Controls for stat='ols' / 'ols_t' / 'lin' / 'lin_t'.	None
absorb	str or list of str	Fixed effects for stat='ols' / 'ols_t', as reghdfe's absorb() (column names, "a^b" for combinations; singletons dropped). They are swept out by the sp.hdfe_ols absorber (Frisch-Waugh-Lovell). Under cluster the regressor of any permutation is R v with R the swept cluster indicators (times interact) and v the permuted cluster values, so only G columns are swept once and each permutation is a small matrix product -- the statistic equals the sp.hdfe_ols coefficient (or its cluster t) refitted on the permuted treatment. Without cluster every permutation re-sweeps the regressor.	None
interact	str	For stat='ols' / 'ols_t': the regressor is treat * data[interact] while treat is what is permuted -- the intensity DiD exposure x post with exposure permuted across clusters.	None

Returns:

Type	Description
dict	'observed': observed test statistic 'p_value': two-sided randomization p-value 'p_one_sided': one-sided (greater) p-value 'n_perms': number of permutations used (the number of distinct assignments when they were enumerated) 'exact': True when every assignment of the design was enumerated 'perm_distribution': array of permuted statistics

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> d = rng.integers(0, 2, 80)
>>> y = 0.5 * d + rng.normal(size=80)
>>> df = pd.DataFrame({"outcome": y, "treatment": d})
>>> result = sp.ri_test(df, y='outcome', treat='treatment',
...                     n_perms=2000, seed=42)
>>> sorted(result.keys())
['exact', 'n_perms', 'observed', 'p_one_sided', 'p_value', 'perm_distribution']
>>> bool(0.0 <= result['p_value'] <= 1.0)
True
>>> result['n_perms']
2000

Notes

Under the sharp null H0: Y_i(1) = Y_i(0) for all i, the treatment assignment is the only source of randomness. The RI p-value is the fraction of permutation statistics at least as extreme as the observed statistic.

For cluster-randomized experiments, treatment is permuted at the cluster level (all units in a cluster get the same permuted status).

See Young (2019, QJE) for why RI p-values should be reported alongside asymptotic p-values in experimental papers.

ksmirnov

ksmirnov(data: DataFrame, y: str, by: str) -> ClassicTestResult

Two-sample Kolmogorov-Smirnov test of equal distributions.

D is the largest difference between the two empirical distribution functions. The one-sided p-values are exp(-2 m n / (m + n) D^2) and the combined one is the limiting Kolmogorov distribution evaluated at D sqrt(m n / (m + n)); both are large-sample approximations, as printed by Stata's ksmirnov (which also offers an exact p-value for small samples, not computed here).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
by	str	Grouping column with exactly two values.	required

Returns:

Type	Description
ClassicTestResult	statistic is the combined D. table has one row per direction: the first group's distribution function above the second's (its values tend to be smaller) and below it.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"y": [1.0, 2, 3, 4, 5, 6, 7, 8],
...                    "g": [0, 0, 0, 0, 1, 1, 1, 1]})
>>> sp.ksmirnov(df, "y", by="g").statistic
1.0

ktau

ktau(data: DataFrame, x: str, y: str) -> ClassicTestResult

Kendall's rank correlations tau-a and tau-b and the score test.

The score S is the number of concordant pairs less the number of discordant ones; tau-a = S / (n (n - 1) / 2) and tau-b divides S by the geometric mean of the numbers of pairs not tied on x and not tied on y. The variance of S under independence is corrected for ties (Kendall 1975), and z = (|S| - 1) / sd carries a continuity correction, as in Stata's ktau.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
x	str		required
y	str		required

Returns:

Type	Description
ClassicTestResult	statistic is tau-b; estimates holds tau_a, score, se_score and z.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"x": [1.0, 2, 3, 4, 5], "y": [2.0, 1, 4, 3, 5]})
>>> round(sp.ktau(df, "x", "y").statistic, 2)
0.6

kwallis

kwallis(data: DataFrame, y: str, by: str) -> ClassicTestResult

Kruskal-Wallis test that several groups have the same distribution.

H = 12 / (N (N + 1)) sum(R_j^2 / n_j) - 3 (N + 1) with R_j the rank sum of group j, referred to chi-squared with k - 1 degrees of freedom. The statistic corrected for ties divides H by 1 - sum(t^3 - t) / (N^3 - N); it is the one R's kruskal.test reports and the second one Stata's kwallis prints.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
by	str	Grouping column.	required

Returns:

Type	Description
ClassicTestResult	statistic / pvalue are the tie-corrected ones; estimates['chi2_unadjusted'] and ['p_unadjusted'] the others.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"y": [1.0, 2, 3, 4, 5, 6, 7, 8, 9],
...                    "g": [1, 1, 1, 2, 2, 2, 3, 3, 3]})
>>> round(sp.kwallis(df, "y", by="g").statistic, 3)
7.2

median_test

median_test(data: DataFrame, y: str, by: str) -> ClassicTestResult

Median test: do the groups have the same share of values above the overall median?

Each observation is classed as greater than the median of the pooled sample or not, and Pearson's chi-squared test of independence is applied to the resulting table. With two groups the continuity corrected statistic is reported as well (Stata median).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
by	str		required

Returns:

Type	Description
ClassicTestResult	table is the 2 x k table of counts.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"y": [1.0, 2, 3, 4, 5, 6, 7, 8],
...                    "g": [0, 0, 0, 0, 1, 1, 1, 1]})
>>> sp.median_test(df, "y", by="g").statistic
8.0

oneway

oneway(data: DataFrame, y: str, by: str, compare: Optional[str] = None) -> ClassicTestResult

One-way analysis of variance.

The F test that the group means are equal, Bartlett's test that the group variances are equal and, on request, every pairwise difference of means with a p-value adjusted for the number of comparisons (Stata oneway, R aov / bartlett.test).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
by	str		required
compare	(None, 'bonferroni', 'sidak', 'scheffe')	Adjustment of the pairwise comparisons. Each difference is the mean of the row group less the mean of the column group, tested with the pooled within-group variance.	None

Returns:

Type	Description
ClassicTestResult	statistic is F with df = (k - 1, N - k). estimates holds the sums of squares, mean squares, bartlett_chi2 / bartlett_p and, with compare, the tables differences and pvalues.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"y": [1.0, 2, 3, 4, 5, 6, 7, 8, 9],
...                    "g": [1, 1, 1, 2, 2, 2, 3, 3, 3]})
>>> sp.oneway(df, "y", by="g").statistic
27.0

ranksum

ranksum(data: DataFrame, y: str, by: str) -> ClassicTestResult

Wilcoxon rank-sum (Mann-Whitney) test that two groups have the same distribution.

The statistic is the sum of the ranks of the first group (the lower value of by) in the pooled sample, ties ranked by their mean rank. Its variance under the null is n1 n2 (N + 1) / 12 less the correction for ties n1 n2 sum(t^3 - t) / (12 N (N - 1)), and z = (T - n1 (N + 1) / 2) / sd is referred to the standard normal without a continuity correction (Stata ranksum; R's wilcox.test(correct = FALSE) gives the same p-value).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	The variable to compare.	required
by	str	Grouping column with exactly two values.	required

Returns:

Type	Description
ClassicTestResult	statistic is z. estimates holds the rank sums, the variances and porder, the estimate of P(y1 > y2) + P(y1 = y2)/2.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"y": [1.0, 3, 5, 7, 2, 4, 6, 9, 11],
...                    "g": [0, 0, 0, 0, 1, 1, 1, 1, 1]})
>>> res = sp.ranksum(df, "y", by="g")
>>> round(res.statistic, 3)
-0.98

robvar

robvar(data: DataFrame, y: str, by: str) -> ClassicTestResult

Robust tests of equal variances across groups.

W0 is Levene's F statistic on the absolute deviations from the group means; W50 replaces the mean by the median and W10 by the 10 percent trimmed mean, here the mean of the values between the group's 10th and 90th percentiles (Brown and Forsythe). Each is referred to F with k - 1 and N - k degrees of freedom (Stata robvar).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
by	str		required

Returns:

Type	Description
ClassicTestResult	statistic / pvalue are those of W0; estimates holds W50, W10 and their p-values.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"y": [1.0, 2, 3, 4, 10, 20, 30, 40],
...                    "g": [0, 0, 0, 0, 1, 1, 1, 1]})
>>> round(sp.robvar(df, "y", by="g").statistic, 3)
9.624

signrank

signrank(data: DataFrame, y: str, other: Optional[str] = None, value: float = 0.0) -> ClassicTestResult

Wilcoxon matched-pairs signed-rank test.

The differences y - other (or y - value) are ranked by absolute size, ties and zeros included, and the ranks of the positive ones are summed. Under the null of a distribution symmetric about zero the sum has mean n (n + 1) / 4 less the share of the zeros and variance n (n + 1) (2n + 1) / 24, less sum(t^3 - t) / 48 for ties among the nonzero differences and n0 (n0 + 1) (2 n0 + 1) / 24 for the n0 zeros (Stata signrank).

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
other	str	The paired variable. Without it y is compared with value.	None
value	float		0

Returns:

Type	Description
ClassicTestResult	statistic is z.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"before": [10.0, 12, 9, 14, 11, 13],
...                    "after": [12.0, 15, 9, 13, 15, 18]})
>>> res = sp.signrank(df, "after", other="before")
>>> res.n_obs
6

spearman

spearman(data: DataFrame, x: str, y: str) -> ClassicTestResult

Spearman's rank correlation and the test that it is zero.

The correlation of the ranks (ties get their mean rank). The p-value refers t = rho sqrt((n - 2) / (1 - rho^2)) to Student's t with n - 2 degrees of freedom, as Stata's spearman and R's cor.test(method = "spearman", exact = FALSE) do.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
x	str		required
y	str		required

Returns:

Type	Description
ClassicTestResult	statistic is rho itself; estimates['t'] the t statistic.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({"x": [1.0, 2, 3, 4, 5], "y": [2.0, 1, 4, 3, 5]})
>>> round(sp.spearman(df, "x", "y").statistic, 2)
0.8

ztest

ztest(data: Optional[DataFrame] = None, y: Optional[str] = None, by: Optional[str] = None, *, other: Optional[str] = None, mu: float = 0.0, sd: Union[float, Sequence[float]] = 1.0, n: Union[None, int, Sequence[int]] = None, mean: Union[None, float, Sequence[float]] = None, alpha: float = 0.05) -> TTestResult

z test on a mean, or on a difference of two means, with known sd.

The population standard deviation is taken as given, so the statistic is standard normal whatever the sample size. Equivalent to Stata's ztest / ztesti. When the standard deviation is estimated from the same data, use :func:ttest.

============================================ ========================= Call Stata ============================================ ========================= sp.ztest(df, "y", mu=20, sd=6) ztest y == 20, sd(6) sp.ztest(df, "y", by="g", sd=6) ztest y, by(g) sd(6) sp.ztest(df, "y", by="g", sd=(5, 7)) ztest y, by(g) sd1(5) sd2(7) sp.ztest(n=10, mean=88, sd=0.71, mu=85) ztesti 10 88 0.71 85 ============================================ =========================

Parameters:

Name	Type	Description	Default
data	DataFrame		None
y	str	The variable whose mean is tested.	None
by	str	Grouping column with exactly two distinct non-missing values; the estimate is the mean in the group with the smaller value minus the mean in the group with the larger one.	None
other	str	A second variable, treated as an independent sample (each column uses its own non-missing rows). A paired z test needs the known correlation of the two and is not offered.	None
mu	float	The null value of the mean, or of the difference.	0.0
sd	float or pair of floats	The known population standard deviation; a pair gives one per sample.	1.0
n	number or pair of numbers	Summary statistics in place of data: observations and sample mean of one sample, or pairs for two.	None
mean	number or pair of numbers	Summary statistics in place of data: observations and sample mean of one sample, or pairs for two.	None
alpha	float	1 - alpha is the confidence level of the intervals.	0.05

Returns:

Type	Description
TTestResult	With df infinite and statistic the z statistic.

Raises:

Type	Description
MethodIncompatibility	When neither data nor summary statistics are given, when by does not have exactly two groups, or when sd is not positive.

Examples:

Scores were N(85, 0.5). After tutoring, ten tests average 88. Did the mean rise?

>>> import statspai as sp
>>> res = sp.ztest(n=10, mean=88, sd=0.5 ** 0.5, mu=85)
>>> round(res.statistic, 2)
13.42
>>> bool(res.pvalue_greater < 0.05)   # Ha: mean > 85
True

vcov

vcov(result: Any, vce: Any, *, data: Optional[DataFrame] = None) -> Any

Re-report a fitted result under another standard-error choice.

Re-runs the call that produced result with vce= replaced and returns the new result. Every other argument -- weights, fixed effects, covariates, seeds -- is the one you passed. Because the estimator itself computes the new standard errors, they are exactly what sp.<estimator>(..., vce=vce) returns, including its small-sample factors and its refusals (an estimator that does not implement a kind raises rather than substituting another).

Parameters:

Name	Type	Description	Default
result	fitted result	Any result from an estimator that accepts vce=.	required
vce	str	Stata-style request, e.g. "robust", "hc3", "cluster state", "cluster state year" (where multiway clustering is offered), "bootstrap". Which kinds a given estimator accepts is listed in sp.function_schema(name).	required
data	DataFrame	The frame the model was fitted on. Needed only when that frame no longer exists (results hold it weakly).	None

Returns:

Type	Description
result	A new result of the same type. model_info['vcov_refit'] records the request and whether the point estimates moved.

Raises:

Type	Description
MethodIncompatibility	When the result has no call record, the estimator has no vce=, the requested kind is not available for it, or replaying the original call does not reproduce the original estimates.

Notes

The original call is replayed once, unchanged, before vce= is swapped; its estimates must match the result's to a relative 1e-10. The point estimates may legitimately move under the new request when the estimator uses the cluster variable for more than the variance (rows with a missing cluster value are dropped; cross-fitting folds and cluster-Lasso loadings follow the clusters); a warning says so.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=200),
...                    "g": rng.integers(0, 40, 200)})
>>> df["y"] = 1 + 2 * df["x"] + rng.normal(size=200)
>>> fit = sp.regress("y ~ x", data=df)
>>> clustered = sp.vcov(fit, "cluster g")
>>> direct = sp.regress("y ~ x", data=df, vce="cluster g")
>>> bool(np.allclose(clustered.std_errors, direct.std_errors, rtol=0, atol=0))
True

vcov_compare

vcov_compare(result: Any, vce: Union[str, Sequence[Any]], *, data: Optional[DataFrame] = None, terms: Optional[Iterable[str]] = None) -> DataFrame

Standard errors of one fit under several vce= choices, side by side.

The first block is the result as fitted; each further block re-runs the call with one vce (see :func:statspai.vcov). A request the estimator does not implement is reported in the note column instead of stopping the table.

Parameters:

Name	Type	Description	Default
result	fitted result		required
vce	str or list	One request or several, e.g. ["robust", "cluster state", "cluster state year"].	required
data	DataFrame	As in :func:statspai.vcov.	None
terms	iterable of str	Keep only these coefficients.	None

Returns:

Type	Description
DataFrame	One row per (request, term): vce, term, estimate, se, ci_lower, ci_upper, pvalue, se_ratio (to the fitted SE of the same term) and note.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(1)
>>> df = pd.DataFrame({"x": rng.normal(size=300),
...                    "g": rng.integers(0, 30, 300)})
>>> df["y"] = 1 + df["x"] + rng.normal(size=300)
>>> fit = sp.regress("y ~ x", data=df)
>>> tab = sp.vcov_compare(fit, ["hc1", "hc3", "cluster g"], terms=["x"])
>>> list(tab["vce"])
['(as fitted)', 'hc1', 'hc3', 'cluster g']

wild_cluster_bootstrap

wild_cluster_bootstrap(data: DataFrame, y: str, x: List[str], cluster: str, test_var: Optional[str] = None, h0: float = 0, n_boot: int = 999, weight_type: str = 'rademacher', seed: Optional[int] = None, alpha: float = 0.05) -> Dict[str, Any]

Wild Cluster Bootstrap p-value and confidence interval.

Provides valid inference when the number of clusters is small, where conventional cluster-robust standard errors are unreliable.

Implements the WCR (Wild Cluster Restricted) bootstrap from Cameron, Gelbach & Miller (2008), with Rademacher weights (default) or Webb (2014) 6-point weights for very few clusters.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
y	str	Outcome variable.	required
x	list of str	All regressors (including the variable being tested).	required
cluster	str	Cluster variable name.	required
test_var	str	Variable to test. Default: last variable in x.	None
h0	float	Null hypothesis value for the test variable coefficient.	0
n_boot	int	Number of bootstrap replications. With Rademacher weights and 2**G <= n_boot the full grid of 2**G sign vectors is enumerated instead (the rule Stata boottest and R fwildclusterboot apply), so the p-value is exact and n_boot in the result reports 2**G.	999
weight_type	str	Bootstrap weight distribution: - 'rademacher': ±1 with equal probability. Standard choice. - 'webb': 6-point distribution from Webb (2014). Recommended when G < 12 clusters. - 'mammen': Mammen (1993) 2-point distribution.	'rademacher'
seed	int	Random seed for reproducibility (unused when the grid is enumerated).	None
alpha	float	Significance level for the confidence interval.	0.05

Returns:

Type

Description

dict

Keys: - beta_hat: OLS point estimate of test_var - se_cluster: conventional cluster-robust SE - t_stat: t-statistic under H0 - p_boot: symmetric two-sided bootstrap p-value, #{|t*| > |t|} / B (strict inequality, as in boottest) - ci_boot: bootstrap percentile-t confidence interval, beta_hat minus the quantiles of t* times se_cluster. It is a valid interval (95% coverage in simulations with 8 and with 20 clusters), but it is not the set of null values the test does not reject: the t* distribution is drawn once, under h0, so with few clusters it can contain values at which p_boot would be below alpha. For the interval that is dual to the test, the one boottest reports, use :func:sp.wild_cluster_ci_inv. - n_clusters: number of clusters - n_boot: number of replications used (2**G if enumerated) - n_boot_requested: the n_boot argument - enumerated: whether the full Rademacher grid was used - weight_type: weight distribution used - recommendation: human-readable guidance

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> education = rng.integers(8, 18, n)
>>> experience = rng.integers(0, 30, n)
>>> wage = (2.0 + 0.4 * education + 0.1 * experience
...         + rng.normal(0, 2, n))
>>> df = pd.DataFrame({'wage': wage, 'education': education,
...                    'experience': experience,
...                    'state': rng.integers(0, 12, n)})
>>> result = sp.wild_cluster_bootstrap(
...     df, y='wage', x=['education', 'experience'],
...     cluster='state', test_var='education', n_boot=199, seed=0)
>>> print(f"Bootstrap p = {result['p_boot']:.4f}")
>>> bool(result['ci_boot'][0] < result['ci_boot'][1])
True

Notes

The WCR bootstrap imposes the null H0: β_test = h0 when generating bootstrap samples. This yields better size control than the unrestricted wild bootstrap, especially with few clusters.

For G < 12 clusters, Webb (2014) weights are recommended over Rademacher. The function will auto-warn if Rademacher is used with very few clusters.

See Cameron, Gelbach & Miller (2008, REStat) Section III for the full algorithm description.

subcluster_wild_bootstrap

subcluster_wild_bootstrap(data: DataFrame, y: str, x: List[str], cluster: str, subcluster: Optional[str] = None, test_var: Optional[str] = None, h0: float = 0.0, n_boot: int = 999, weight_type: str = 'webb', seed: Optional[int] = None, alpha: float = 0.05) -> Dict[str, Any]

Subcluster wild cluster bootstrap for few-treated-clusters.

When treatment varies within cluster (or when you want to re-expand the randomization at a finer grain), sign-flips happen at the sub-cluster level. SEs are still computed clustered at the coarse cluster level.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str		required
x	list of str		required
cluster	str	Primary cluster column (for SE computation).	required
subcluster	str or None	Finer grouping at which sign-flips occur. If None, every observation is its own subcluster (pure Rademacher at obs level).	None
test_var	str or None	Parameter to test; default last element of x.	None
h0	float	Null value.	0.0
n_boot	int	Bootstrap replications.	999
weight_type	('rademacher', 'webb', 'mammen')	Distribution of sign flips. 'webb' (6-point) recommended when treatment has <= 5 treated clusters.	'rademacher'
seed	int		None
alpha	float		0.05

Returns:

Type	Description
dict with ``p_boot`` (symmetric, ``#{|t*| > |t|} / B`` as in
``boottest``), ``ci_boot``, ``beta_hat``, ``t_stat``, ``se_cluster``,
``n_subclusters``, ``n_boot`` (``2**S`` when Rademacher weights and
``2**S <= n_boot`` trigger full enumeration), ``enumerated`` and
``recommendation``.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> cl = rng.integers(0, 10, size=n)
>>> treat = (cl < 3).astype(float)
>>> x1 = rng.normal(size=n)
>>> y = 1.0 + 0.4 * treat + 0.5 * x1 + rng.normal(size=n)
>>> df = pd.DataFrame({"y": y, "treat": treat, "x1": x1, "cl": cl})
>>> out = sp.subcluster_wild_bootstrap(
...     df, y="y", x=["treat", "x1"], cluster="cl",
...     test_var="treat", n_boot=199, weight_type="webb", seed=0,
... )
>>> 0.0 <= out["p_boot"] <= 1.0
True

wild_cluster_ci_inv

wild_cluster_ci_inv(data: DataFrame, y: str, x: List[str], cluster: str, test_var: Optional[str] = None, n_boot: int = 999, weight_type: str = 'webb', alpha: float = 0.05, grid_size: int = 41, grid_span: float = 6.0, seed: Optional[int] = None) -> Dict[str, Any]

Confidence interval via bootstrap p-value inversion.

Runs the wild cluster bootstrap repeatedly across a grid of null values and finds the boundary where the two-sided bootstrap p-value equals alpha. This yields a WCR-inverted CI with better small-cluster coverage than the percentile-t CI.

The grid is centered on the OLS point estimate with half-width grid_span * se_cluster and grid_size evenly-spaced points. It only brackets each endpoint: the bootstrap p-value is a step function of the null value, and the bracket is then bisected to machine precision to locate the jump where p falls below alpha. With Rademacher weights and 2**G <= n_boot the bootstrap grid is enumerated, so the interval does not depend on a seed. Stata boottest interpolates its endpoints between grid points instead of locating the jump; on a 12-cluster example with the same enumerated draws the two sets differ by 0.2% and 1.9% of a standard error.

Shares the data / model arguments of :func:subcluster_wild_bootstrap; the grid-specific parameters are documented below.

Parameters:

Name	Type	Description	Default
grid_size	int	Number of null-value grid points to evaluate (odd preferred).	41
grid_span	float	Half-width of the search grid in units of cluster-robust SE.	6.0

Returns:

Type	Description
dict with ``ci``, ``p_grid`` (grid_size,), ``h0_grid``, ``beta_hat``,
``se_cluster``.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> cl = rng.integers(0, 10, size=n)
>>> treat = (cl < 3).astype(float)
>>> x1 = rng.normal(size=n)
>>> y = 1.0 + 0.4 * treat + 0.5 * x1 + rng.normal(size=n)
>>> df = pd.DataFrame({"y": y, "treat": treat, "x1": x1, "cl": cl})
>>> out = sp.wild_cluster_ci_inv(
...     df, y="y", x=["treat", "x1"], cluster="cl",
...     test_var="treat", n_boot=199, grid_size=21, seed=0,
... )
>>> out["p_grid"].shape
(21,)
>>> len(out["ci"])
2