statspai.postestimation

postestimation

Post-estimation tools for StatsPAI.

Provides: - margins(): Average Marginal Effects (AME), Marginal Effects at the Mean (MEM) - marginsplot(): Visualize marginal effects - test(): Wald / F test for linear restrictions (beta1 = beta2, joint significance) - lincom(): Linear combinations of coefficients with inference - nlcom(): Nonlinear combinations by the delta method - predictions() / comparisons() / slopes() (+ avg_ versions) / hypotheses(): marginaleffects-style quantities with delta-method inference

Equivalent to Stata's margins, test, lincom commands.

postestimation_contract

postestimation_contract(result: Any, *, data: Optional[DataFrame] = None, include_diagnostics: bool = True) -> Dict[str, Any]

Return the post-estimation actions supported by result.

Parameters:

Name	Type	Description	Default
result	object	Fitted StatsPAI result, estimator, or compatible object.	required
data	DataFrame	Analysis frame available for data-dependent actions. When supplied, the contract marks margins/predict paths as ready rather than merely method-available.	None
include_diagnostics	bool	Include scalar diagnostics from model_info / diagnostics.	True

Returns:

Type	Description
dict	Machine-readable contract with available, missing, recommended_next, and optional diagnostics keys.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> df = pd.DataFrame({
...     "x1": rng.normal(size=200),
...     "x2": rng.normal(size=200),
... })
>>> df["y"] = (1.0 + 0.5 * df["x1"] - 0.3 * df["x2"]
...            + rng.normal(size=200))
>>> res = sp.regress("y ~ x1 + x2", data=df)
>>> con = sp.postestimation_contract(res, data=df)
>>> con["result_type"]
'EconometricResults'
>>> "margins" in con["available"]   # data supplied -> ready
True
>>> con["has_data"]
True

lincom

lincom(result: Any, expression: str, alpha: float = 0.05) -> dict[str, Any]

Estimate a linear combination of coefficients with inference.

Parameters:

Name	Type	Description	Default
result	EconometricResults or CausalResult	Fitted model. After irr=True the combination is of the coefficients b (Stata's lincom without its irr option).	required
expression	str	Linear combination (Stata lincom syntax). Examples: - "x1 + x2" — beta_x1 + beta_x2 - "x1 - x2" — beta_x1 - beta_x2 - "2*x1 + 3*x2" — 2beta_x1 + 3beta_x2 - "_cons + x1" — intercept plus slope	required
alpha	float	Significance level.	0.05

Returns:

Type	Description
dict	estimate, se, statistic (t when the fit has residual degrees of freedom, else z; also under the historical key z), df (None for z), distribution ("t" or "z"), pvalue, ci and expression.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> x1, x2 = rng.normal(size=n), rng.normal(size=n)
>>> y = 1.0 + 0.5 * x1 + 0.5 * x2 + rng.normal(size=n)
>>> df = pd.DataFrame({"y": y, "x1": x1, "x2": x2})
>>> result = sp.regress("y ~ x1 + x2", data=df)
>>> out = sp.lincom(result, "x1 + x2")    # beta1 + beta2
>>> sorted(out)
['ci', 'df', 'distribution', 'estimate', 'expression',
 'pvalue', 'se', 'statistic', 'z']
>>> diff = sp.lincom(result, "x1 - x2")   # beta1 - beta2
>>> bool(diff["se"] > 0)
True

test

test(result: Any, hypothesis: Union[str, Sequence[str]]) -> dict[str, Any]

Wald test for linear restrictions on coefficients.

Parameters:

Name	Type	Description	Default
result	EconometricResults or CausalResult	Fitted model with .params and .std_errors; restrictions that involve several coefficients also need its covariance matrix. After a fit reported as incidence-rate ratios (irr=True) the restrictions are on the coefficients b, not on exp(b), as Stata's test after poisson, irr.	required
hypothesis	str or list of str	Hypothesis specification (Stata test syntax), or a list of them to be tested jointly (["x1 = 0", "x2 = 2*x3"]). A coefficient is named as the result names it, including names such as I(exper ** 2) or np.log(sales). Examples: - "x1 = 0" — beta_x1 = 0 - "x1 = x2" — beta_x1 = beta_x2 - "x1 = x2 = 0" — joint: beta_x1 = beta_x2 and beta_x2 = 0 - "x1 x2" — joint: beta_x1 = 0 and beta_x2 = 0 - "(x1 = 0) (x2 + x3 = 1)" — grouped restrictions - "x1 = 0, x2 + x3 = 1" — the same, in the comma-separated spelling of statsmodels' f_test - "x1 - 2*x2 = 0" — linear restriction; _cons names the intercept whatever the estimator calls it. - "2.race = 3.race", "1.union#c.grade" — Stata's factor-variable names for the coefficients the formula engine calls C(race)[T.2] and grade:C(union)[T.1]; "i.race" on its own tests every level of the factor jointly (Stata's testparm i.race).	required

Returns:

Type	Description
dict	statistic (F when the fit has residual degrees of freedom, else chi-squared), pvalue, df = (q, df_resid or None), chi2 (the Wald statistic), distribution ("F" or "chi2") and hypothesis. Linearly dependent restrictions are dropped, as Stata does; q counts the independent ones.

Raises:

Type	Description
MethodIncompatibility	Unknown coefficient names, restrictions without coefficients, inconsistent restrictions, or a multi-coefficient restriction on a result that carries no covariance matrix.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> x1, x2, x3 = rng.normal(size=n), rng.normal(size=n), rng.normal(size=n)
>>> y = 1.0 + 0.5 * x1 + 0.5 * x2 - 0.3 * x3 + rng.normal(size=n)
>>> df = pd.DataFrame({"y": y, "x1": x1, "x2": x2, "x3": x3})
>>> result = sp.regress("y ~ x1 + x2 + x3", data=df)
>>> out = sp.test(result, "x1 = x2")        # beta1 = beta2?
>>> sorted(out)
['chi2', 'df', 'distribution', 'hypothesis', 'pvalue', 'statistic']
>>> joint = sp.test(result, "x1 = x2 = 0")  # joint: beta1 = beta2 = 0?
>>> bool(joint["pvalue"] < 0.05)
True
>>> restr = sp.test(result, "x1 + x2 = 1")  # beta1 + beta2 = 1?

avg_comparisons

avg_comparisons(result: Any, variables: Any = None, newdata: Optional[DataFrame] = None, *, by: Optional[Union[str, Sequence[str]]] = None, comparison: str = 'difference', type: str = 'response', wts: Any = None, df: Any = None, conf_level: float = 0.95, vcov: bool = True) -> DataFrame

Average contrasts (marginaleffects::avg_comparisons).

The average over rows of the unit-level contrasts of :func:comparisons (comparison='ratioavg' gives the ratio of average predictions instead). For a binary treatment on the response scale this is the average treatment effect implied by the model (the g-formula).

Parameters are those of :func:comparisons.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(3)
>>> d = pd.DataFrame({"x": rng.normal(size=400), "d": rng.integers(0, 2, 400)})
>>> d["y"] = (rng.uniform(size=400) < 1 / (1 + np.exp(-(d.x + d.d)))).astype(int)
>>> fit = sp.logit("y ~ d + x", data=d)
>>> ate = sp.avg_comparisons(fit, "d")
>>> bool(0 < ate.estimate.iloc[0] < 1)
True

avg_predictions

avg_predictions(result: Any, newdata: Optional[DataFrame] = None, *, by: Optional[Union[str, Sequence[str]]] = None, type: str = 'response', wts: Any = None, df: Any = None, conf_level: float = 0.95, vcov: bool = True) -> DataFrame

Average prediction, overall or by= group (marginaleffects::avg_predictions).

Parameters are those of :func:predictions.

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=200), "g": rng.integers(0, 3, 200)})
>>> d["y"] = 1 + d.x + d.g + rng.normal(size=200)
>>> fit = sp.regress("y ~ x + C(g)", data=d)
>>> sp.avg_predictions(fit, by="g").shape[0]
3

avg_slopes

avg_slopes(result: Any, variables: Any = None, newdata: Optional[DataFrame] = None, *, by: Optional[Union[str, Sequence[str]]] = None, slope: str = 'dydx', type: str = 'response', wts: Any = None, df: Any = None, conf_level: float = 0.95, eps: float = 1e-05, vcov: bool = True) -> DataFrame

Average partial derivatives (marginaleffects::avg_slopes).

The average marginal effect; with wts=True and df='model' it is Stata's margins, dydx() after a weighted fit.

Parameters are those of :func:slopes.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(5)
>>> d = pd.DataFrame({"x": rng.normal(size=300)})
>>> d["y"] = (rng.uniform(size=300) < 1 / (1 + np.exp(-d.x))).astype(int)
>>> fit = sp.logit("y ~ x", data=d)
>>> ame = sp.avg_slopes(fit, "x")
>>> bool(0 < ame.estimate.iloc[0] < fit.params["x"])
True

comparisons

comparisons(result: Any, variables: Any = None, newdata: Optional[DataFrame] = None, *, by: Optional[Union[str, Sequence[str]]] = None, comparison: str = 'difference', type: str = 'response', wts: Any = None, df: Any = None, conf_level: float = 0.95, vcov: bool = True) -> DataFrame

Unit-level contrasts of predictions (marginaleffects::comparisons).

For each variable, predictions at a "high" and a "low" value of it, every other column held at the row's own value, contrasted per row:

• a factor (C(g), strings, categories): each level against the first;
• a numeric 0/1 variable: 1 against 0;
• another numeric variable: an increase of 1 from the row's value, or as given in variables=: a number (that increase), "sd" / "2sd" (mean +/- sd/2, mean +/- sd), "iqr" (third against first quartile), "minmax" (maximum against minimum) or a pair (low, high) of fixed values.

Parameters:

Name	Type	Description	Default
result	fitted result	As in :func:predictions.	required
variables	(str, list or dict)	Variables to contrast (default: every model variable); a dict maps a numeric variable to its change.	None
newdata	DataFrame	Rows to evaluate (default: the estimation sample).	None
by	str or list of str	Average within groups.	None
comparison	str	'difference', 'ratio', 'lnratio' (per row); 'ratioavg', 'lnratioavg' (ratio of the average predictions). When the contrasts are averaged (avg_comparisons or by=), 'ratio' and 'lnratio' mean 'ratioavg' and 'lnratioavg', as in marginaleffects.	'difference'
type	str	As in :func:predictions.	'response'
wts	str	As in :func:predictions.	'response'
df	str	As in :func:predictions.	'response'
conf_level	str	As in :func:predictions.	'response'
vcov	str	As in :func:predictions.	'response'

Returns:

Type	Description
DataFrame	term, contrast, the data's columns (unit level) or the by columns, and the estimate columns.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(2)
>>> d = pd.DataFrame({"x": rng.normal(size=300), "d": rng.integers(0, 2, 300)})
>>> d["y"] = 1 + 2 * d.d + d.x + rng.normal(size=300)
>>> fit = sp.regress("y ~ d + x", data=d)
>>> sp.comparisons(fit, "d").shape[0] == len(d)
True

hypotheses

hypotheses(result: Any, hypothesis: Union[str, Sequence[str], None] = None, *, df: Any = None, conf_level: float = 0.95, joint: Union[bool, Sequence[str]] = False, joint_test: str = 'f') -> DataFrame

Delta-method test of a function of estimates (marginaleffects::hypotheses).

result is a fitted model (the estimates are its coefficients) or a table returned by :func:predictions, :func:comparisons or :func:slopes (the estimates are its rows, with their joint covariance from the table's Jacobian). b1, b2 ... refer to the estimates by position; a model's coefficients can also be named. "lhs = rhs" tests lhs - rhs = 0; an expression without = is tested against 0.

Parameters:

Name	Type	Description	Default
result	fitted result or marginaleffects-style table		required
hypothesis	str or list of str	E.g. "b1 - b2 = 0", "b2 / b3 = 1", "x1 = 2 * x2". Only arithmetic and the functions sp.nlcom accepts. A list gives one row per hypothesis (or, with joint=True, one joint test of all of them). May be omitted when joint names the estimates.	None
df	(None, 'model' or float)	Reference distribution of the one-at-a-time tests (None: normal).	None
conf_level	float		0.95
joint	bool or list of str	Joint Wald test instead of one row per hypothesis. True with hypothesis tests all of them at once -- linear or not, as Stata testnl does (the restrictions' Jacobian enters W = g' (G V G')^-1 g); True without hypothesis tests that every estimate is 0; a list of coefficient names (or b1 ...) tests that those are 0 (marginaleffects::hypotheses(joint=)).	False
joint_test	('f', 'chisq')	"f": F = W / q on (q, df2) degrees of freedom, df2 the fit's residual degrees of freedom (marginaleffects' convention, also after a logit) or df when given; "chisq": W on q.	"f"

Returns:

Type	Description
DataFrame	One row per hypothesis: term (the hypothesis), estimate (lhs - rhs) and the inference columns. A joint test returns one row: hypothesis (the restrictions joined by "; "), statistic, p_value, df1 and (F) df2, test.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(6)
>>> d = pd.DataFrame({"x1": rng.normal(size=300), "x2": rng.normal(size=300)})
>>> d["y"] = 1 + d.x1 + d.x2 + rng.normal(size=300)
>>> fit = sp.regress("y ~ x1 + x2", data=d)
>>> h = sp.hypotheses(fit, "b2 = b3")
>>> bool(h.p_value.iloc[0] > 0.01)
True
>>> sp.hypotheses(fit, ["x1 = 1", "x2 = 1"]).shape[0]
2
>>> j = sp.hypotheses(fit, ["x1 = 1", "x2 = 1"], joint=True)
>>> list(j.columns[:3])
['hypothesis', 'statistic', 'p_value']

predictions

predictions(result: Any, newdata: Optional[DataFrame] = None, *, by: Optional[Union[str, Sequence[str]]] = None, type: str = 'response', wts: Any = None, df: Any = None, conf_level: float = 0.95, vcov: bool = True) -> DataFrame

Unit-level predictions with delta-method standard errors.

marginaleffects::predictions: one row per row of newdata (default: the estimation sample), with the data's columns and estimate, std_error, statistic, p_value, conf_low, conf_high. With by= the predictions are averaged within groups (what :func:avg_predictions does).

Parameters:

Name	Type	Description	Default
result	fitted result	From sp.regress, sp.logit, sp.probit, sp.cloglog, sp.poisson, sp.nbreg or sp.glm (identity, logit, probit, cloglog or log link).	required
newdata	DataFrame	Rows to predict for. Default: the estimation sample of the fitted data.	None
by	str or list of str	Average within groups of these columns.	None
type	('response', 'link')	Prediction scale.	'response'
wts	(None, True, str or array)	Weights for averages: None (unweighted, the marginaleffects default), True (the fit's weight column), a column name or an array.	None
df	(None, 'model' or float)	Reference distribution: None is normal (marginaleffects), 'model' the fit's own (t after regress), or a number.	None
conf_level	float		0.95
vcov	bool	False returns the estimates alone (no standard errors), as marginaleffects' vcov = FALSE. A switch, not a variance type: the variance is the fit's own. After a fit with absorbed fixed effects (sp.feols / fepois / feglm with | fe) the effects enter at their estimates. An average (or by= average) that depends on them after a Poisson fit (sp.fepois) has standard errors that count the effects as parameters, as Stata's margins after poisson y x i.g; unit-level quantities and other models then need vcov=False (marginaleffects refuses all of them on a fixest model). Contrasts and slopes of a linear model do not depend on the effects and keep their standard errors.	True

Returns:

Type	Description
DataFrame

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> d = pd.DataFrame({"x": rng.normal(size=300), "g": rng.integers(0, 2, 300)})
>>> d["y"] = (rng.uniform(size=300) < 1 / (1 + np.exp(-d.x))).astype(int)
>>> fit = sp.logit("y ~ x + g", data=d)
>>> p = sp.predictions(fit)
>>> bool(((p.estimate > 0) & (p.estimate < 1)).all())
True

slopes

slopes(result: Any, variables: Any = None, newdata: Optional[DataFrame] = None, *, by: Optional[Union[str, Sequence[str]]] = None, slope: str = 'dydx', type: str = 'response', wts: Any = None, df: Any = None, conf_level: float = 0.95, eps: float = 1e-05, vcov: bool = True) -> DataFrame

Unit-level partial derivatives (marginaleffects::slopes).

dY/dX of the prediction with respect to each numeric variable, differentiating through interactions and formula transformations; slope='eyex' / 'eydx' / 'dyex' give elasticities and semi-elasticities. Factor and 0/1 variables are reported as contrasts, as marginaleffects does.

Parameters:

Name	Type	Description	Default
result	Any	As in :func:comparisons.	required
variables	Any	As in :func:comparisons.	required
newdata	Any	As in :func:comparisons.	required
by	Any	As in :func:comparisons.	required
type	Any	As in :func:comparisons.	required
wts	Any	As in :func:comparisons.	required
df	Any	As in :func:comparisons.	required
conf_level	Any	As in :func:comparisons.	required
vcov	Any	As in :func:comparisons.	required
slope	('dydx', 'eyex', 'eydx', 'dyex')		'dydx'
eps	float	Relative step of the central difference used on transformed terms.	1e-5

Returns:

Type	Description
DataFrame

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(4)
>>> d = pd.DataFrame({"x": rng.normal(size=300)})
>>> d["y"] = 1 + 2 * d.x + 0.5 * d.x**2 + rng.normal(size=300)
>>> fit = sp.regress("y ~ x + I(x**2)", data=d)
>>> s = sp.slopes(fit, "x")
>>> bool(np.allclose(s.estimate, fit.params.iloc[1] + 2 * fit.params.iloc[2] * d.x))
True

contrast

contrast(result: Any, data: DataFrame, variable: str, method: str = 'r', reference: Any = None, alpha: float = 0.05, effect: str = 'difference', subset: Any = None) -> DataFrame

Compute contrasts of predictive margins across levels of a variable.

Equivalent to Stata's margins <var>, contrast(ar) / contrast(r).

Parameters:

Name	Type	Description	Default
result	EconometricResults	Fitted model result.	required
data	DataFrame	Estimation data.	required
variable	str	Categorical variable whose levels are contrasted.	required
method	str	Contrast type: 'r' (reference): each level vs reference level. 'ar' (adjacent): each level vs the previous level. 'gw' (grand-mean weighted): each level vs the weighted grand mean of all levels.	'r'
reference	scalar	Reference level when method='r'. Defaults to the smallest observed level.	None
alpha	float	Significance level.	0.05
effect	(difference, ratio, odds_ratio)	Scale of the contrast between two predictive margins m1 and m0. 'ratio' is m1 / m0 (a marginal risk ratio after a binary-outcome model) and 'odds_ratio' is [m1 / (1 - m1)] / [m0 / (1 - m0)], the marginal odds ratio, which unlike the coefficient of a logistic regression is a population-averaged quantity. Inference is by the delta method on the log scale: se is the standard error of the log ratio, the test is of a ratio of one, and the interval is the exponentiated log-scale interval (R marginaleffects comparison = "lnratioavg" / "lnoravg" with transform = exp). Not available with method='gw'.	'difference'
subset	str or boolean array-like	Average over a subpopulation of data: a boolean mask aligned with its rows, or an expression for data.eval such as "treat == 1". The model's factor levels still come from the whole of data, so with "treat == 1" the contrast between the levels of treat is the effect on the treated and with "treat == 0" the effect on the untreated (Stata margins, subpop(); R marginaleffects newdata = subset(...)). Passing an already subsetted frame as data is not the same: a level that is absent from it cannot be contrasted.	None

Returns:

Type	Description
DataFrame	Columns: contrast_label, contrast, se, z, pvalue, ci_lower, ci_upper.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n = 300
>>> group = rng.integers(0, 3, size=n).astype(float)
>>> df = pd.DataFrame({
...     "group": group,
...     "experience": rng.normal(10, 3, size=n),
... })
>>> df["wage"] = (10 + 2.0 * df["group"]
...               + 0.3 * df["experience"] + rng.normal(size=n))
>>> result = sp.regress("wage ~ group + experience", data=df)
>>> c = sp.contrast(result, data=df, variable="group",
...                 method="r", reference=0)
>>> c["contrast_label"].tolist()
['1.0 vs 0', '2.0 vs 0']

event_study_table

event_study_table(result: Any, *, regex: Optional[str] = None, label_fmt: str = 't={t}', include_reference: bool = False) -> _MarginsResult

Adapter that turns an event-study fit into a regtable input.

Two extraction paths:

1. CausalResult fast path — when result.model_info carries an "event_study" DataFrame (the canonical shape produced by :func:sp.event_study), read its relative_time / estimate / se / ci_lower / ci_upper / pvalue columns directly. Reference period (estimate==0 and se==0) is dropped by default (set include_reference=True to keep it visible).

2. Regex path — when regex is provided, scan result.params.index for coefficients matching the pattern and use the first capture group as the relative time.

Rows are sorted by relative time (so the rendered table reads t=-3, t=-2, …, t=+3 regardless of the underlying coefficient name encoding).

Parameters:

Name	Type	Description	Default
result	CausalResult or EconometricResults	Fitted model. The CausalResult fast path is used automatically when model_info['event_study'] is present.	required
regex	str	Pattern with one capture group that matches the relative time in coefficient names. Required when the CausalResult fast path is not applicable. Examples: r"^tau_(-?\d+)$", r"^lag(\d+)$", r"::(-?\d+)$".	None
label_fmt	str	Format string for the row label of each event-time bin. The {t} placeholder receives the integer relative time.	``"t={t}"``
include_reference	bool	Whether to render the reference period row (typically t=-1) where the estimate is identically zero.	False

Returns:

Type	Description
_MarginsResult	Duck-typed result accepted directly by :func:sp.regtable.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(7)
>>> rows = []
>>> for u in range(60):
...     tt = rng.choice([3, 4, 5, 100])  # 100 marks never-treated
...     for t in range(8):
...         post = 1.0 if (tt < 100 and t >= tt) else 0.0
...         y = 0.3 * u / 60 + 0.1 * t + 1.5 * post + rng.normal(scale=0.5)
...         rows.append({"unit": u, "time": t, "treat_time": tt, "y": y})
>>> panel = pd.DataFrame(rows)
>>> r = sp.event_study(panel, y="y", treat_time="treat_time",
...                    time="time", unit="unit", window=(-3, 3))
>>> tbl = sp.event_study_table(r)
>>> _ = sp.regtable(tbl, title="Event study", output="text")

margins_at

margins_at(result: Any, data: DataFrame, at: Dict[str, Any], alpha: float = 0.05, contrast: Optional[str] = None, subset: Any = None, predict: Optional[str] = None, category: Any = None, vcov: bool = True) -> DataFrame

Compute predictive margins at specific covariate values.

Equivalent to Stata's margins, at(experience=(1 5 10 15 20)).

For each combination of at values, every observation has the at variables set to those values while all other covariates stay at their observed levels. The prediction -- on the model's default scale: the linear prediction after regress, Pr(y=1) after logit / probit, the expected count after poisson -- is averaged over the estimation sample (weighted when the fit was) to give the predictive margin at that point, with delta-method SEs and t(df) / z inference as the fit's own.

Parameters:

Name	Type	Description	Default
result	EconometricResults	Fitted model result (must have .params and associated data).	required
data	DataFrame	Data to compute margins on.	required
at	dict	Mapping of variable names to lists/arrays of values. If multiple variables are given, the Cartesian product of all value lists is used. Example:: at={"experience": [1, 5, 10], "female": [0, 1]} produces 6 grid points. Every key must be a model variable.	required
alpha	float	Significance level for confidence intervals.	0.05
contrast	(None, first, adjacent)	Return differences between grid points instead of the margins: each point against the first one ('first') or against the point before it ('adjacent'). The standard error uses the covariance between the two margins, which the separate rows do not show. With at={"dose": [30, 60]} this is g-computation of the effect of moving everyone from 30 to 60, whatever shape the model gives the dose (R marginaleffects::avg_comparisons( variables = list(dose = c(30, 60)))).	None
subset	str or boolean array-like	Average over a subpopulation of data; see :func:margins.	None
predict	optional	The prediction after a limited-dependent-variable model and the outcome categories to report; see :func:margins. Multiple-outcome predictions add an outcome column.	None
category	optional	The prediction after a limited-dependent-variable model and the outcome categories to report; see :func:margins. Multiple-outcome predictions add an outcome column.	None
vcov	bool	False returns the margins without standard errors. After a fit with absorbed fixed effects the predictive margins depend on the estimated effects (and so do contrasts after a non-linear link); after a Poisson fit their standard errors count the effects as parameters (Stata), otherwise they are refused unless vcov=False; see :func:margins.	True

Returns:

Type	Description
DataFrame	One row per grid point with columns for each at variable, plus margin, se, ci_lower, ci_upper. With contrast, one row per comparison: the at columns of the point, versus (the point it is compared with), contrast, se, z, pvalue, ci_lower, ci_upper.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> df = pd.DataFrame({
...     "experience": rng.integers(0, 20, size=200).astype(float),
...     "female": rng.integers(0, 2, size=200).astype(float),
... })
>>> df["wage"] = (10.0 + 0.5 * df["experience"]
...               - 1.0 * df["female"] + rng.normal(size=200))
>>> result = sp.regress(
...     "wage ~ experience + female + experience:female", data=df
... )
>>> m = sp.margins_at(result, data=df, at={"experience": [1, 5, 10, 15, 20]})
>>> m.shape[0]
5
>>> m.columns.tolist()
['experience', 'margin', 'se', 'ci_lower', 'ci_upper']

margins_at_plot

margins_at_plot(margins_at_df: DataFrame, x: Optional[str] = None, by: Optional[str] = None, ax: Any = None, figsize: Tuple[float, float] = (8, 5), title: Optional[str] = None, xlabel: Optional[str] = None, ylabel: str = 'Predicted Value', palette: Optional[List[str]] = None) -> Tuple[Any, Any]

Plot predictive margins from margins_at() with confidence bands.

Parameters:

Name	Type	Description	Default
margins_at_df	DataFrame	Output from margins_at().	required
x	str	Variable to place on the x-axis. If None, inferred as the at-variable with the most unique values.	None
by	str	Variable to produce separate lines for (legend grouping).	None
ax	matplotlib Axes		None
figsize	tuple		(8, 5)
title	str		None
xlabel	str		None
ylabel	str		'Predicted Value'
palette	list of str	Colours for each by group.	None

Returns:

Type	Description
(fig, ax)

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(42)
>>> df = pd.DataFrame({
...     "experience": rng.integers(0, 20, size=200).astype(float),
...     "female": rng.integers(0, 2, size=200).astype(float),
... })
>>> df["wage"] = (10.0 + 0.5 * df["experience"]
...               - 1.0 * df["female"] + rng.normal(size=200))
>>> res = sp.regress(
...     "wage ~ experience + female + experience:female", data=df
... )
>>> m = sp.margins_at(
...     res, data=df, at={"experience": [1, 5, 10, 15, 20]}
... )
>>> fig, ax = sp.margins_at_plot(m, x="experience")
>>> fig.savefig("margins_at.png")

margins_table

margins_table(result: Any, data: Optional[DataFrame] = None, variables: Optional[List[str]] = None, at: Optional[Dict[str, Any]] = None, method: str = 'ame', eps: float = 1e-05, alpha: float = 0.05) -> _MarginsResult

Marginal-effects result that pipes straight into sp.regtable.

Thin adapter over :func:margins: computes the marginal effects (same kwargs), then wraps the resulting DataFrame in a :class:_MarginsResult exposing params / std_errors / tvalues / pvalues / conf_int_* so that sp.regtable(margins_table(model)) produces a publication-quality marginal-effects table — closing the "estimator → margins table" gap that previously required users to hand-build add_rows.

Mirrors the R workflow modelsummary(avg_slopes(model)).

Parameters:

Name	Type	Description	Default
result	EconometricResults	Fitted model — same input sp.margins accepts.	required
data	Optional[DataFrame]	Forwarded verbatim to :func:margins.	None
variables	Optional[DataFrame]	Forwarded verbatim to :func:margins.	None
at	Optional[DataFrame]	Forwarded verbatim to :func:margins.	None
method	Optional[DataFrame]	Forwarded verbatim to :func:margins.	None
eps	Optional[DataFrame]	Forwarded verbatim to :func:margins.	None
alpha	Optional[DataFrame]	Forwarded verbatim to :func:margins.	None

Returns:

Type	Description
_MarginsResult	A duck-typed result object usable directly as a positional argument to :func:sp.regtable.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> x = rng.normal(size=n)
>>> z = rng.normal(size=n)
>>> p = 1 / (1 + np.exp(-(0.5 * x - 0.4 * z)))
>>> df = pd.DataFrame({
...     "y": (rng.uniform(size=n) < p).astype(int),
...     "x": x,
...     "z": z,
... })
>>> m = sp.logit("y ~ x + z", data=df)
>>> mt = sp.margins_table(m)
>>> _ = sp.regtable(mt, output="text")
>>> sp.regtable(mt, output="latex", filename="margins.tex")

marginsplot

marginsplot(margins_df: DataFrame, ax: Any = None, figsize: Tuple[float, float] = (8, 5), color: str = '#2C3E50', title: Optional[str] = None) -> Tuple[Any, Any]

Plot marginal effects with confidence intervals.

Parameters:

Name	Type	Description	Default
margins_df	DataFrame	Output from margins().	required
ax	matplotlib Axes		None
figsize	tuple		(8, 5)
color	str		'#2C3E50'
title	str		None

Returns:

Type	Description
(fig, ax)

Examples:

>>> import statspai as sp
>>> df = sp.cps_wage()
>>> result = sp.regress("log_wage ~ education + experience + female",
...                     data=df)
>>> me = sp.margins(result, data=df)
>>> fig, ax = sp.marginsplot(me, title='AME of wage determinants')
>>> fig.savefig('margins.png')

pwcompare

pwcompare(result: Any, data: DataFrame, variable: str, adjust: str = 'none', alpha: float = 0.05) -> DataFrame

Pairwise comparisons of predictive margins across all levels.

Equivalent to Stata's pwcompare <var>.

Parameters:

Name	Type	Description	Default
result	EconometricResults	Fitted model result.	required
data	DataFrame	Estimation data.	required
variable	str	Categorical variable whose levels are compared pairwise.	required
adjust	str	P-value adjustment method: 'none': unadjusted. 'bonferroni': Bonferroni correction. 'sidak': Sidak correction. 'holm': Holm step-down procedure.	'none'
alpha	float	Significance level for (adjusted) confidence intervals.	0.05

Returns:

Type	Description
DataFrame	Columns: comparison, diff, se, z, pvalue, pvalue_adj, ci_lower, ci_upper.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> n = 300
>>> group = rng.integers(0, 3, size=n).astype(float)
>>> df = pd.DataFrame({
...     "group": group,
...     "experience": rng.normal(10, 3, size=n),
... })
>>> df["wage"] = (10 + 2.0 * df["group"]
...               + 0.3 * df["experience"] + rng.normal(size=n))
>>> result = sp.regress("wage ~ group + experience", data=df)
>>> pw = sp.pwcompare(result, data=df, variable="group",
...                   adjust="bonferroni")
>>> bool((pw["pvalue_adj"] >= pw["pvalue"]).all())
True