statspai.bounds

bounds

Bounds and Partial Identification for causal effects.

When point identification is not possible (e.g., due to sample selection, non-compliance, or missing data), bounds provide informative intervals for the treatment effect.

Functions:

Name	Description
- **Lee Bounds** : Bounds on ATE under sample selection
- **Manski Bounds** : Worst-case bounds with minimal assumptions
- **Horowitz-Manski Bounds** : Tighter bounds conditioning on covariates
- **IV Bounds** : Bounds under imperfect instruments
- **Oster Delta** : Coefficient stability / identified set
- **Selection Bounds** : Lee bounds with covariates
- **Breakdown Frontier** : Assumption robustness frontier

References

Lee, D. S. (2009). Training, Wages, and Sample Selection: Estimating Sharp Bounds on Treatment Effects. RES, 76(3), 1071-1102. [@lee2009training]

Manski, C. F. (1990). Nonparametric Bounds on Treatment Effects. AER P&P, 80(2), 319-323. [@manski1990nonparametric]

Horowitz, J. L. & Manski, C. F. (2000). Nonparametric Analysis of Randomized Experiments with Missing Covariate and Outcome Data. JASA, 95(449), 77-84. [@horowitz2000nonparametric]

Nevo, A. & Rosen, A. M. (2012). Identification with Imperfect Instruments. Review of Economics and Statistics, 94(3), 659-671. doi:10.1162/REST_a_00171. [@nevo2012identification]

Oster, E. (2019). Unobservable Selection and Coefficient Stability. JBES, 37(2), 187-204. [@oster2019unobservable]

Masten, M. A. & Poirier, A. (2021). Salvaging Falsified Instrumental Variable Models. Econometrica, 89(3), 1449-1469. [@masten2021salvaging]

BalkePearlResult dataclass

Bases: ResultProtocolMixin

Result of :func:statspai.balke_pearl.

Fields: lower, upper, width, lower_monotone, upper_monotone, joint_probs, n_obs, detail.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> z = rng.integers(0, 2, 400)
>>> d = (z & (rng.random(400) < 0.8)).astype(int)
>>> y = (rng.random(400) < 0.3 + 0.4 * d).astype(int)
>>> res = sp.balke_pearl(pd.DataFrame({"y": y, "d": d, "z": z}),
...                      y="y", treat="d", instrument="z")
>>> isinstance(res, sp.BalkePearlResult)
True

MLBoundsResult dataclass

Bases: ResultProtocolMixin

ATE bounds produced by :func:ml_bounds.

Examples:

>>> import numpy as np, pandas as pd, statspai as sp
>>> rng = np.random.default_rng(0)
>>> df = pd.DataFrame({"x": rng.normal(size=300), "d": rng.integers(0, 2, 300)})
>>> df["y"] = df.x + 0.5 * df.d + rng.normal(size=300)
>>> res = sp.ml_bounds(df, y="y", treat="d", covariates=["x"],
...                    n_splits=2, n_bootstrap=20, random_state=0)
>>> isinstance(res, sp.MLBoundsResult)
True

center_shift

center_shift() -> float

Midpoint difference from the plain-means Manski estimate.

Both estimate the same interval (by iterated expectations the covariate-averaged bounds are the Manski bounds), so this is finite-sample estimation noise of the cross-fitted nuisances and tends to zero as n grows. It is not identifying information.

BoundsResult dataclass

Bases: ResultProtocolMixin

Result container for partial identification / bounds estimation.

Attributes:

Name	Type	Description
lower	float	Lower bound of the identified set.
upper	float	Upper bound of the identified set.
se_lower	float	Standard error of the lower bound.
se_upper	float	Standard error of the upper bound.
ci_lower	tuple	Confidence interval for the lower bound (lower_lo, lower_hi).
ci_upper	tuple	Confidence interval for the upper bound (upper_lo, upper_hi).
method	str	Name of the bounding method.
alpha	float	Significance level used for confidence intervals.
n_obs	int	Number of observations used.
model_info	dict	Additional method-specific information.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> age = rng.integers(20, 60, n).astype(float)
>>> education = rng.integers(8, 18, n).astype(float)
>>> trained = rng.integers(0, 2, n)
>>> wage = 10 + 0.3 * age + 0.5 * education + 5 * trained + rng.normal(0, 3, n)
>>> df = pd.DataFrame(
...     {"wage": wage, "trained": trained,
...      "age": age, "education": education})
>>> res = sp.horowitz_manski(
...     data=df, y="wage", treatment="trained",
...     covariates=["age", "education"], n_boot=50, random_state=0)
>>> isinstance(res, sp.BoundsResult)
True
>>> bool(res.lower <= res.upper)
True
>>> bool(res.width >= 0)
True

width property

width: float

Width of the identified set.

midpoint property

midpoint: float

Midpoint of the identified set.

parameter_ci

parameter_ci() -> Tuple[float, float]

Imbens-Manski (2004) confidence interval for the parameter.

[lower - c se_lower, upper + c se_upper] with c solving Phi(c + width / max(se)) - Phi(-c) = 1 - alpha: between the one-sided z_{1-alpha} (wide set) and the two-sided z_{1-alpha/2} (point identified). The same interval sp.lee_bounds reports. Until 2026-10-27 the printed interval was labelled Imbens-Manski but used z_{1-alpha/2}, which covers the whole identified set and over-covers the parameter.

includes_zero

includes_zero() -> bool

Check whether the identified set includes zero.

summary

summary() -> str

Return a formatted summary string.

plot

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

Interval plot showing the identified set with confidence intervals.

Parameters:

Name	Type	Description	Default
ax	matplotlib Axes	Axes to draw on; created if None.	None
**kwargs	Any	Passed to ax.errorbar.	{}

Returns:

Type	Description
Axes

lee_bounds

lee_bounds(data: DataFrame, y: str, treat: str, selection: Optional[str] = None, covariates: Optional[List[str]] = None, n_bootstrap: int = 500, alpha: float = 0.05, random_state: int = 42, se_method: str = 'bootstrap', trimming: str = 'quantile', weights: Optional[str] = None) -> CausalResult

Compute Lee (2009) bounds for ATE under sample selection.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data (including units with missing outcomes).	required
y	str	Outcome variable (may have NaN for selected-out units).	required
treat	str	Binary treatment variable (0/1).	required
selection	str	Binary selection/retention indicator (1 = observed, 0 = missing). None (Stata leebounds without select()) counts a unit as selected when its outcome is observed. Rows with a missing treatment or selection indicator, or selected with a missing outcome, are dropped with a warning, as leebounds drops them.	None
covariates	list of str	Discrete covariates for Lee's (2009) tightened bounds (Stata leebounds, tight()): each covariate cell is trimmed by its own share and the trimmed means are averaged over the covariate distribution of the always-observed (the retained units of the less-retained arm). Continuous covariates must be binned first. Cell-level results are in model_info['cells']. Analytic SEs of tightened bounds follow leebounds' approximation and need trimming='leebounds'; the bootstrap works with every trimming.	None
n_bootstrap	int	Bootstrap iterations for inference.	500
alpha	float	Significance level.	0.05
random_state	int		42
se_method	(bootstrap, analytic)	Standard errors of the two bounds. 'analytic' is the asymptotic variance of Lee (2009, Proposition 3) with sample analogues; with trimming='leebounds' it equals Stata leebounds' default vce(analytic). 'bootstrap' resamples units n_bootstrap times.	'bootstrap'
trimming	(quantile, exact, leebounds)	How the p = |p1 - p0| / max(p1, p0) share is trimmed from the arm with higher retention. 'quantile' is Lee's (2009) sample estimator: keep every retained outcome at or beyond the sample p-quantile (Stata _pctile definition). 'exact' keeps exactly (1 - p) n units of mass: the observations tied at the quantile share the mass still needed, the rule Tauchmann describes for leebounds with tied outcomes. On a binary or coarsely discrete outcome use 'exact': the quantile rule keeps every observation tied at the quantile, so it trims too little and the two bounds collapse (a warning says so when that happens). 'leebounds' reproduces Stata leebounds (Tauchmann) number for number, floating-point artefacts included -- use it to match a published leebounds table (see Notes).	'quantile'
weights	str	Sampling (probability) weights, Stata leebounds [pw=]: the retention rates, the trimming quantile, the kept mass, the means and the analytic variance are all weighted (weights rescaled to mean one). Needs trimming='leebounds' (which also holds the weights in single precision, as leebounds does) or 'exact'; not available with covariates. The bootstrap resamples units together with their weights. One difference from leebounds: when no observation sits at the threshold of the tail that keeps the upper outcomes, leebounds' standard error of that bound uses the unweighted variance of the untrimmed arm; here that variance is weighted for both bounds (tests/reference_parity/test_lee_bounds_weights_Stata_parity.py rebuilds Stata's number).	None

Notes

Before 1.29 StatsPAI kept floor((1 - p) n) whole observations, one fewer than Lee's quantile rule whenever p n is not an integer, so each bound was off by one observation's worth of mass.

Run as a black box, leebounds behaves as if the trimming percentage and the quantile threshold were rounded to the 16 or fewer significant digits Stata keeps when it stores a number in a macro. When the rounded threshold still equals an observed outcome, the tied observations share the remaining mass (trimming='exact'); when it does not -- common for continuous outcomes stored as float -- every outcome at or beyond the rounded value is kept, which can drop the quantile observation itself. On the UCT replication (QJE 2016, Table III assets) that moves the lower bound from -2.51 (Lee's rule) to Stata's -3.38. trimming='leebounds' reproduces this; docs/dev/cleanroom/leebounds.md records the experiments.

Returns:

Type	Description
CausalResult	estimate = midpoint of bounds. ci = Imbens-Manski confidence interval for the identified set. model_info contains 'lower_bound' and 'upper_bound'.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n = 600
>>> training = rng.integers(0, 2, size=n)
>>> # employment (selection) is higher among the treated -> differential attrition
>>> employed = (rng.uniform(size=n)
...             < np.where(training == 1, 0.8, 0.6)).astype(int)
>>> wage = 10.0 + 2.0 * training + rng.normal(size=n)
>>> wage = np.where(employed == 1, wage, np.nan)   # wage missing if not employed
>>> df = pd.DataFrame({'wage': wage, 'training': training, 'employed': employed})
>>> result = sp.lee_bounds(df, y='wage', treat='training',
...                        selection='employed', n_bootstrap=100)
>>> bool(result.model_info['lower_bound'] <= result.model_info['upper_bound'])
True

manski_bounds

manski_bounds(data: DataFrame, y: str, treat: str, y_lower: Optional[float] = None, y_upper: Optional[float] = None, assumption: str = 'none', alpha: float = 0.05, n_bootstrap: int = 500, random_state: int = 42) -> CausalResult

Compute Manski (1990) worst-case bounds on ATE.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
y	str	Outcome variable.	required
treat	str	Binary treatment variable (0/1).	required
y_lower	float	Known lower bound of the outcome. If None, uses observed min.	None
y_upper	float	Known upper bound of the outcome. If None, uses observed max.	None
assumption	str	Additional assumption: - 'none': no assumptions (widest bounds) - 'mtr': Monotone Treatment Response (Y(1) >= Y(0) for all): [0, worst-case upper]. - 'mts': Monotone Treatment Selection with positive selection, E[Y(d) | D=1] >= E[Y(d) | D=0] (Manski & Pepper 2000): [worst-case lower, E[Y|D=1] - E[Y|D=0]]. - 'mts_mtr': MTS and MTR jointly: [0, E[Y|D=1] - E[Y|D=0]]. Raises ValueError when E[Y|D=1] < E[Y|D=0]: the two assumptions are then jointly refuted and the identified set is empty. Before 1.29, 'mts' returned the 'mts_mtr' interval (i.e. silently imposed MTR as well) and, when the naive difference was negative, swapped its endpoints instead of reporting the empty set. These match the "worst case", "MTS positive selection" and "MTS and MTR positive selection" rows of Stata's tebounds at zero misclassification.	'none'
alpha	float		0.05
n_bootstrap	int		500
random_state	int		42

Returns:

Type	Description
CausalResult

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(2)
>>> n = 600
>>> training = rng.integers(0, 2, size=n)
>>> employed = (rng.uniform(size=n)
...            < np.where(training == 1, 0.7, 0.5)).astype(int)
>>> df = pd.DataFrame({'employed': employed, 'training': training})
>>> result = sp.manski_bounds(df, y='employed', treat='training',
...                           y_lower=0, y_upper=1, n_bootstrap=100)
>>> bool(result.model_info['lower_bound'] <= result.model_info['upper_bound'])
True

breakdown_frontier

breakdown_frontier(estimate: float, se: float, assumption: str = 'parallel_trends', max_violation: float = 0.1, n_grid: int = 100, alpha: float = 0.05) -> BoundsResult

Masten-Poirier (2021) breakdown frontier for qualitative conclusions.

For a given point estimate and standard error, computes how much an identifying assumption can be violated before a qualitative conclusion (e.g., "positive treatment effect") breaks down.

Parameters:

Name	Type	Description	Default
estimate	float	Point estimate of the treatment effect.	required
se	float	Standard error of the estimate.	required
assumption	str	Label for the identifying assumption being relaxed. Currently supports a generic linear violation model applicable to 'parallel_trends', 'exclusion_restriction', or 'selection_on_observables'.	'parallel_trends'
max_violation	float	Maximum magnitude of the assumption violation to explore.	0.1
n_grid	int	Grid resolution for the frontier.	100
alpha	float		0.05

Returns:

Type	Description
BoundsResult	lower/upper give the identified set at the maximum violation. model_info contains 'breakdown_point' (the violation magnitude at which the conclusion reverses) and 'frontier_grid' / 'frontier_bounds' for plotting.

Examples:

>>> import statspai as sp
>>> result = sp.breakdown_frontier(
...     estimate=0.05, se=0.02,
...     assumption="parallel_trends",
...     max_violation=0.1,
... )
>>> round(result.model_info["breakdown_point"], 4)
0.05
>>> (round(result.lower, 4), round(result.upper, 4))
(-0.05, 0.15)
>>> ax = result.plot()

horowitz_manski

horowitz_manski(data: DataFrame, y: str, treatment: str, covariates: List[str], y_lower: Optional[float] = None, y_upper: Optional[float] = None, n_boot: int = 500, alpha: float = 0.05, random_state: int = 42) -> BoundsResult

Horowitz-Manski (2000) bounds conditioning on covariates.

Averages the conditional worst-case (Manski) bounds over the covariate strata. For the ATE with a common outcome range [y_lower, y_upper] these bounds are linear in the cell means, so the average equals the unconditional worst-case bounds exactly (sp.manski_bounds(assumption='none')); conditioning tightens nothing unless the strata carry further restrictions.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
y	str	Outcome variable.	required
treatment	str	Binary treatment variable (0/1).	required
covariates	list of str	Covariates to condition on (discretised via quartiles for continuous variables).	required
y_lower	float	Known lower bound of Y. Defaults to observed min.	None
y_upper	float	Known upper bound of Y. Defaults to observed max.	None
n_boot	int	Bootstrap replications.	500
alpha	float		0.05
random_state	int		42

Returns:

Type	Description
BoundsResult

Notes

Bootstrap replicates that fail are dropped from the SE/CI computation; a RuntimeWarning reports the failed fraction and model_info['n_boot_failed'] records the count.

Examples:

>>> import statspai as sp, numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> trained = rng.binomial(1, 0.5, n)
>>> wage = 50 + 10 * trained + rng.normal(0, 15, n)
>>> df = pd.DataFrame({"wage": wage, "trained": trained,
...                    "age": rng.normal(40, 10, n)})
>>> result = sp.horowitz_manski(
...     data=df, y="wage", treatment="trained", covariates=["age"],
...     y_lower=float(df["wage"].min()), y_upper=float(df["wage"].max()),
... )
>>> result.method
'Horowitz-Manski Bounds (2000)'
>>> bool(result.lower < 0 < result.upper)   # worst-case bounds straddle zero
True

iv_bounds

iv_bounds(data: DataFrame, y: str, treatment: str, instrument: str, controls: Optional[List[str]] = None, assumption: str = 'monotone_iv', alpha: float = 0.05, n_boot: int = 500, random_state: int = 42, xu_sign: str = 'positive') -> BoundsResult

Nevo-Rosen (2012) bounds on a linear IV slope with an imperfect instrument.

The model is y = beta x + W delta + u with one endogenous regressor x (treatment), exogenous controls W and one instrument z that may be correlated with u. Nevo and Rosen replace exogeneity of z with two assumptions:

• A3 (same direction): corr(x, u) and corr(z, u) have the same sign;
• A4 (less endogenous): |corr(z, u)| <= |corr(x, u)|.

Under A3 and A4, V(1) = sd(z) x - sd(x) z is a second instrument whose bias has the opposite sign, and the sharp identified set for beta has the IV estimates with z and with V(1) as endpoints (working paper, Proposition 2; Proposition 3 with controls). It is a closed interval when the two instruments bound beta from opposite sides, and a half-line otherwise (typically when z and x are positively correlated). Which side is open follows from the sign of corr(x, u), which the data cannot tell and xu_sign states.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable.	required
treatment	str	Endogenous regressor (continuous or binary).	required
instrument	str	Imperfect instrument (one variable).	required
controls	list of str	Exogenous regressors W, partialled out of y, x and z (they serve as their own instruments).	None
assumption	(monotone_iv, less_than_late)	'monotone_iv': Nevo-Rosen A3 and A4, the sharp set above. 'less_than_late': deprecated. It returned Wald +/- |OLS - Wald|, an interval with no source in the literature that can exclude the true slope under the assumption it was described with; it is kept, unchanged up to rounding, until the next minor release.	'monotone_iv'
alpha	float		0.05
n_boot	int		500
random_state	int		42
xu_sign	(positive, negative)	The assumed sign of corr(x, u): 'positive' when OLS is believed to overstate beta. It decides which IV estimate is the upper and which the lower endpoint.	'positive'

Returns:

Type	Description
BoundsResult	lower / upper (-inf / inf for an open side). model_info holds beta_iv_z, beta_iv_v1 (the two endpoints), ols_estimate, wald_estimate (= beta_iv_z), closed and first_stage_f (the F statistic of z in the regression of x on z and the controls, homoskedastic).

Notes

Before 2026-10-27 'monotone_iv' returned [min(OLS, IV), max(OLS, IV)]. When z and x are positively correlated both estimates are biased in the same direction and the true slope lies outside that interval: with corr(x, u) > 0 it lies below both. first_stage_f was n corr(z, x)^2.

Bootstrap replicates that fail are dropped from the SE/CI computation; a RuntimeWarning reports the failed fraction and model_info['n_boot_failed'] records the count. An open side has a NaN standard error.

The formulas are those of the working paper (NBER w14434 and cemmap CWP16/08, Section 3, both read 2026-10-27); the published version is Review of Economics and Statistics 94(3), 659-671.

References

nevo2012identification

Examples:

>>> import statspai as sp, numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 2000
>>> z, v, e = rng.normal(size=(3, n))
>>> u = 0.2 * z + 0.8 * v + 0.5 * e      # z is mildly endogenous
>>> x = -z + v                            # negatively correlated with z
>>> df = pd.DataFrame({"y": x + u, "x": x, "z": z})
>>> result = sp.iv_bounds(data=df, y="y", treatment="x", instrument="z",
...                       n_boot=50)
>>> result.method
'Nevo-Rosen IV Bounds (monotone_iv)'
>>> bool(result.lower <= 1.0 <= result.upper)     # true slope 1
True

oster_delta

oster_delta(data: DataFrame, y: str, x_base: List[str], x_controls: List[str], r_max: float = 1.3, delta_range: Tuple[float, float] = (-2.0, 2.0), n_grid: int = 200, alpha: float = 0.05, n_boot: int = 500, random_state: int = 42) -> BoundsResult

Oster (2019) coefficient stability bounds and delta* computation.

Computes the identified set for the treatment coefficient beta under the assumption that selection on unobservables is proportional to selection on observables (delta), and R-squared would be at most r_max if all unobservables were included.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable.	required
x_base	list of str	Key treatment/variable(s) of interest.	required
x_controls	list of str	Additional controls whose inclusion tightens identification.	required
r_max	float	Maximum R-squared assumption. A value in (R_full, 1] is used as R_max itself. A value above 1 cannot be an R-squared and is read as a multiplier of the controlled R-squared, so the default 1.3 gives Oster's recommended R_max = min(1, 1.3 * R_full); <= 0 means the same. A value in (0, R_full] falls back to that recommendation with a warning.	1.3
delta_range	tuple	Range of proportional selection parameter delta.	(-2, 2)
n_grid	int	Grid points for delta in the identified set computation.	200
alpha	float		0.05
n_boot	int		500
random_state	int		42

Returns:

Type

Description

BoundsResult

lower/upper give the identified set for beta at delta=1 (equal selection) and the given r_max. model_info['beta_star_delta1'] and model_info['delta_star'] are Oster's exact solutions (the quadratic / cubic of her Proposition 2 as solved by her Stata psacalc, which they reproduce to machine precision). The first entry of x_base is the treatment; further x_base entries enter both the short and the long regression (psacalc's mcontrol()).

Notes

Bootstrap replicates that fail are dropped from the SE/CI computation; a RuntimeWarning reports the failed fraction and model_info['n_boot_failed'] records the count.

Examples:

>>> import statspai as sp, numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 500
>>> education = rng.normal(12, 3, n)
>>> experience = rng.normal(10, 5, n)
>>> tenure = rng.normal(5, 3, n)
>>> wage = (2 * education + 1.5 * experience + 0.8 * tenure
...         + rng.normal(0, 5, n))
>>> df = pd.DataFrame({"wage": wage, "education": education,
...                    "experience": experience, "tenure": tenure})
>>> result = sp.oster_delta(
...     data=df, y="wage",
...     x_base=["education"],
...     x_controls=["experience", "tenure"],
...     r_max=1.3,
... )
>>> result.method
'Oster (2019) Coefficient Stability Bounds'
>>> bool(result.lower <= result.upper)
True

selection_bounds

selection_bounds(data: DataFrame, y: str, treatment: str, selection: str, covariates: Optional[List[str]] = None, method: str = 'conditional', n_boot: int = 500, alpha: float = 0.05, random_state: int = 42) -> BoundsResult

Lee (2009) bounds for ATE under sample selection, optionally conditioning on covariates for tighter bounds.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Outcome variable (may have NaN when selection=0).	required
treatment	str	Binary treatment (0/1).	required
selection	str	Binary indicator: 1 = outcome observed, 0 = missing.	required
covariates	list of str	Covariates to condition on for tighter (conditional) bounds.	None
method	(conditional, unconditional)	'conditional': Lee's (2009) tightened bounds. Each covariate stratum is trimmed by its own share and the trimmed means are averaged over the covariate distribution of the always-observed (the selected units of the less-retained arm). Continuous covariates (float, or more than 10 values) are binned into quartiles first. Equals sp.lee_bounds(covariates=<strata>, trimming='quantile'). 'unconditional': Lee bounds ignoring covariates, equal to sp.lee_bounds(trimming='quantile').	'conditional'
n_boot	int		500
alpha	float		0.05
random_state	int		42

Returns:

Type	Description
BoundsResult

Notes

Since 2026-10-27 both methods run the sp.lee_bounds code. Before, the unconditional bounds kept floor((1 - p) n) observations (one fewer than Lee's quantile rule when p n is not whole, the rule sp.lee_bounds dropped in 1.29), and the conditional bounds weighted strata by their share of the whole sample and dropped strata missing an arm, instead of weighting by the always-observed. Rows with a missing treatment / selection indicator, or selected with a missing outcome, are dropped with a warning, as in sp.lee_bounds.

Bootstrap replicates that fail (or yield NaN bounds, e.g. a stratum losing all treated or selected units in a resample) are dropped from the SE/CI computation; a RuntimeWarning reports the failed fraction and model_info['n_boot_failed'] records the count.

Examples:

>>> import statspai as sp, numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 500
>>> trained = rng.binomial(1, 0.5, n)
>>> employed = rng.binomial(1, 0.8, n)
>>> wage = 50 + 10 * trained + rng.normal(0, 10, n)
>>> df = pd.DataFrame({"wage": np.where(employed == 1, wage, np.nan),
...                    "trained": trained, "employed": employed,
...                    "age": rng.normal(40, 10, n)})
>>> result = sp.selection_bounds(
...     data=df, y="wage", treatment="trained",
...     selection="employed",
...     covariates=["age"],
...     method="conditional",
... )
>>> result.method
'Lee (2009) Bounds (conditional)'
>>> bool(result.lower <= result.upper)
True