The scenarios model — at vs atexog

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from pymargins import Margins

rng = np.random.default_rng(42)
n = 2000
df = pd.DataFrame({
    "age": rng.integers(20, 75, n),
    "female": rng.binomial(1, 0.52, n),
    "treated": rng.binomial(1, 0.40, n),
})
lp = -1.5 + 0.04 * df["age"] - 0.3 * df["female"] + 0.8 * df["treated"]
df["y"] = rng.binomial(1, 1 / (1 + np.exp(-lp)))

fit = smf.glm("y ~ age + female + treated", data=df,
              family=sm.families.Binomial()).fit()
m = Margins.log_scale(fit, at="overall")

Two knobs control where in covariate space a margin is evaluated. They live at different levels of the API.

at — session-level aggregation rule

at= is set at session construction. It controls the default evaluation rule for variables not otherwise pinned:

Value	Per-variable behavior
"overall"	take the observed value on every row, then average
"typical"	median for continuous, mode for discrete
"mean"	mean for all (errors on non-numeric)
"median"	median for all
"mode"	mode for all (errors on continuous)
dict	per-variable override, with _default as a fallback rule
callable	(data) -> 1-row DataFrame, fully bespoke

"overall" corresponds to Stata’s bare margins and gives AAP / AME. "typical" corresponds to margins, atmeans for mixed factor / continuous models and gives APM / MEM.

atexog — per-call counterfactual pins

atexog= is a per-call dict (or a list of dicts wrapped as scenarios=) that pins specific variables to specific values. A list-valued entry produces a Cartesian product (a grid). Variables not mentioned in atexog follow the session’s at= rule.

# AAP at age=25, 45, 65, averaging the rest over the sample
print(Margins.log_scale(fit, at="overall").predict(
    atexog={"age": [25, 45, 65]}
).summary())

# APR at age=25, 45, 65, others held at typical profile
print(Margins.log_scale(fit, at="typical").predict(
    atexog={"age": [25, 45, 65]}
).summary())

============================================================
             Margins Result (delta, level=0.95)             
============================================================
        estimate  std err         z  P>|z|  [95% Conf. Int.]
------------------------------------------------------------
age=25    0.4120   0.0453  -19.5680  0.000    0.3770, 0.4503
age=45    0.6124   0.0184  -26.7196  0.000    0.5908, 0.6349
age=65    0.7819   0.0173  -14.1802  0.000    0.7558, 0.8089
============================================================

n = 2000
Note: std err is on the inference scale; estimate and CI are on the reporting scale.
κ: max=0.062
Delta-vs-sim disagreement: 0.156%

============================================================
             Margins Result (delta, level=0.95)             
============================================================
        estimate  std err         z  P>|z|  [95% Conf. Int.]
------------------------------------------------------------
age=25    0.3009   0.0726  -16.5385  0.000    0.2610, 0.3469
age=45    0.5005   0.0381  -18.1675  0.000    0.4645, 0.5393
age=65    0.6999   0.0288  -12.3895  0.000    0.6615, 0.7405
============================================================

n = 2000
Note: std err is on the inference scale; estimate and CI are on the reporting scale.
κ: max=0.067
Delta-vs-sim disagreement: 0.962%

Why split it this way?

Because the aggregation choice is a methodological commitment (AME vs MEM is an argument; if you flip mid-analysis the audit trail should show it) but the counterfactual pins are not (you genuinely do want to evaluate the same AME at several age points).

This is the same logic behind keeping phi, vcov, level, and method session-level: the analytical posture belongs in the constructor; the analytical question belongs in the method call.

See The session pre-commitment.