Demos

Two end-to-end walkthroughs ship in the repository root. Each one is a single self-contained script you can run with python demo_<name>.py after installing smmargins.

Williams-style logit walkthrough

demo_margins.py reproduces, on a simulated voting dataset, every core statistic in Richard Williams’ Margins01 notes:

1. Adjusted predictions at specific values (APR / margins, at(...))

2. Adjusted predictions at means (APM / margins, atmeans)

3. Average adjusted predictions (AAP / margins)

4. Marginal effects at representative values (MER / margins, dydx(..) at(..))

5. Marginal effects at means (MEM / margins, dydx(..) atmeans)

6. Average marginal effects (AME / margins, dydx(..))

7. Discrete contrasts for a categorical variable

8. Williams’ classic interaction example: AME of age separately by sex

Highlights from the script

Fit a logit with an interaction:

fit = smf.logit(
    "voted ~ age + income + C(educ) + female + age:female",
    data=df,
).fit()
M = Margins(fit)

APR — predictions at policy-relevant ages, averaging everything else over the sample:

M.predict(atexog={"age": [25, 45, 65]})

MER, MEM, and AME for age — these can differ meaningfully in nonlinear models with interactions:

M.dydx("age", atexog={"age": [25, 45, 65]})   # MER
M.dydx("age", at="mean")                      # MEM
M.dydx("age")                                 # AME

Discrete AME for a multi-level factor with an explicit reference level:

M.dydx("educ", reference="college")

Williams’ interaction lesson — same model, AME of age for each sex:

M.dydx("age", atexog={"female": [0, 1]})

Full source

"""
demo_margins.py
===============

Walkthrough of the core analyses in Richard Williams' *Margins01* notes
(https://academicweb.nd.edu/~rwilliam/stats/Margins01.pdf), implemented
on top of StatsModels + patsy + the ``smmargins`` package.

We'll reproduce, in turn:

  1. Adjusted predictions at specific values (APR / "margins, at(...)")
  2. Adjusted predictions at means (APM / "margins, atmeans")
  3. Average adjusted predictions (AAP / "margins")
  4. Marginal effects at representative values (MER / "margins, dydx(..) at(..)")
  5. Marginal effects at means (MEM / "margins, dydx(..) atmeans")
  6. Average marginal effects (AME / "margins, dydx(..)")
  7. Discrete changes for categorical variables
  8. The interaction example that Williams uses to motivate AME
"""

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

from smmargins import Margins

pd.options.display.width = 120
pd.options.display.float_format = "{: .4f}".format

# ---------------------------------------------------------------------------
# Simulate a binary-outcome dataset with structure similar to Williams' notes
# ---------------------------------------------------------------------------
rng = np.random.default_rng(7)
N = 5_000
df = pd.DataFrame(
    {
        "age":    rng.normal(45, 12, N).clip(18, 90),
        "income": rng.lognormal(10.5, 0.4, N),          # ~36k median
        "educ":   rng.choice(["hs", "college", "grad"], N, p=[0.4, 0.4, 0.2]),
        "female": rng.integers(0, 2, N),
    }
)
eta = (
    -4.0
    + 0.05 * df["age"]
    + 0.00001 * df["income"]
    + 0.8 * (df["educ"] == "college")
    + 1.4 * (df["educ"] == "grad")
    + 0.3 * df["female"]
    - 0.0004 * df["age"] * (df["female"])        # interaction
)
df["voted"] = (rng.uniform(0, 1, N) < 1 / (1 + np.exp(-eta))).astype(int)

print("Sample:")
print(df.head(3), "\n")

# ---------------------------------------------------------------------------
# Fit a logit with an interaction, like the Williams example
# ---------------------------------------------------------------------------
fit = smf.logit(
    "voted ~ age + income + C(educ) + female + age:female",
    data=df,
).fit(disp=False)
print("=" * 80)
print("Fitted logit")
print("=" * 80)
print(fit.summary().tables[1])
print()

# `analytic=True` is the default: the outer ∂g/∂β goes through
# `family.link.inverse_deriv` for any GLM (Logit/Probit/Poisson/...) and
# the identity link for OLS/WLS/GLS, falling back to central finite
# differences only when the link derivative isn't available. Set
# `analytic=False` to force FD; you'll get the same answers (see the
# parity check at the bottom of this file) but pay p extra forward
# predict() calls per statistic.
M = Margins(fit)

# ---------------------------------------------------------------------------
# 1. Adjusted predictions at representative values (APR)
#    Stata: margins, at(age=(25 45 65))
# ---------------------------------------------------------------------------
print("=" * 80)
print("1. APR  (predict at age=25,45,65; everything else at sample values)")
print("=" * 80)
print(M.predict(atexog={"age": [25, 45, 65]}))
print()

# ---------------------------------------------------------------------------
# 2. Adjusted prediction at means (APM)  vs  average adjusted prediction (AAP)
# ---------------------------------------------------------------------------
print("=" * 80)
print("2. APM  (margins, atmeans)   vs   AAP  (margins)")
print("=" * 80)
print("APM:"); print(M.predict(at="mean"))
print("\nAAP:"); print(M.predict())
print()

# ---------------------------------------------------------------------------
# 3. Marginal effect: MER vs MEM vs AME for `age`
#    (Williams points out these three can differ meaningfully in nonlinear
#    models with interactions)
# ---------------------------------------------------------------------------
print("=" * 80)
print("3. d Pr(voted)/d age : MER (at age=25,45,65),  MEM, and AME")
print("=" * 80)
print("MER (at age=25,45,65):")
print(M.dydx("age", atexog={"age": [25, 45, 65]}))
print("\nMEM (at means of everything):")
print(M.dydx("age", at="mean"))
print("\nAME (averaged over the sample):")
print(M.dydx("age"))
print()

# ---------------------------------------------------------------------------
# 4. Discrete contrast for the categorical variable `educ`
# ---------------------------------------------------------------------------
print("=" * 80)
print("4. Discrete AME for educ  (each level vs 'college' as reference)")
print("=" * 80)
print(M.dydx("educ", reference="college"))
print()

# ---------------------------------------------------------------------------
# 5. Discrete change for the dummy `female`  (auto-detected as discrete)
# ---------------------------------------------------------------------------
print("=" * 80)
print("5. AME for female (0/1 dummy):  Pr(voted|female=1) - Pr(voted|female=0)")
print("=" * 80)
print(M.dydx("female"))
print()

# ---------------------------------------------------------------------------
# 6. Interaction-sensitivity: marginal effect of age, separately for men/women
#    This is Williams' classic motivating example: the interaction coefficient
#    alone tells you little about what the marginal effect actually is for any
#    given subpopulation.
# ---------------------------------------------------------------------------
print("=" * 80)
print("6. AME of age, separately by sex  (Williams' interaction illustration)")
print("=" * 80)
print(M.dydx("age", atexog={"female": [0, 1]}))
print()

# ---------------------------------------------------------------------------
# 7. Adjusted predictions, age by sex — table suitable for plotting
# ---------------------------------------------------------------------------
print("=" * 80)
print("7. Predicted Pr(voted) over age, for each sex")
print("=" * 80)
tbl = M.predict(atexog={"age": list(range(20, 91, 10)), "female": [0, 1]})
print(tbl)

# ---------------------------------------------------------------------------
# 8. Analytic vs FD: same answers, faster path
#    Logit exposes `family.link.inverse_deriv`, so the analytic outer
#    Jacobian is used by default. Toggling `analytic=False` reroutes
#    every statistic through central finite differences — useful as a
#    sanity check or when working with a custom Link subclass that
#    doesn't implement inverse_deriv.
# ---------------------------------------------------------------------------
print()
print("=" * 80)
print("8. Analytic vs FD — same numbers, taken via different paths")
print("=" * 80)
M_fd = Margins(fit, analytic=False)
ame_an = M.dydx("age")
ame_fd = M_fd.dydx("age")
print(f"AME(age) analytic : est={ame_an.estimate[0]: .8f}  se={ame_an.se[0]: .8f}")
print(f"AME(age) FD       : est={ame_fd.estimate[0]: .8f}  se={ame_fd.se[0]: .8f}")
print(f"max abs diff      : "
      f"est {abs(ame_an.estimate[0] - ame_fd.estimate[0]): .2e}, "
      f"se {abs(ame_an.se[0] - ame_fd.se[0]): .2e}")

Healthcare-style 2x2 difference-in-differences

demo_did.py answers a clinical question:

Is there a rate difference of condition \(X\) between groups A and B, with or without a preexisting condition \(Y\)?

The script fits a logit on simulated patient data and reports, on the probability scale:

• 4 cell predictions \(P(X \mid \text{group}, Y)\)

• 2 simple effects \(P(X \mid B, Y) - P(X \mid A, Y)\) at each \(Y\)

• 1 difference-in-differences (whether the A-vs-B gap depends on \(Y\))

All with delta-method standard errors and confidence intervals. The DiD here is not the coefficient on the group:Y interaction — that coefficient is on the log-odds scale, while the clinical question is about probabilities. This is Ai & Norton (2003) in practice; see Mathematical motivation for the derivation.

Highlights from the script

Fit and call did():

fit = smf.logit(
    "condition_X ~ C(group) + preexist_Y + C(group):preexist_Y "
    "+ age + female",
    data=df,
).fit()
M = Margins(fit)

did = M.did("group", "preexist_Y",
            group_levels=["A", "B"],
            condition_levels=[0, 1])
print(did)              # cells + simple effects + DiD

Same DiD at one specific patient profile (60-year-old male):

M.did("group", "preexist_Y",
      group_levels=["A", "B"], condition_levels=[0, 1],
      atexog={"age": 60, "female": 0})

Plot-ready cell table:

tbl = did.cells.summary()        # estimate / SE / CI per cell

Full source

"""
demo_did.py
===========

Difference-in-differences example directly matching the question:

    "Is there a rate difference of condition X between group A and B,
     with or without preexisting condition Y?"

We fit a logit for P(X=1) on group (A/B), preexisting Y (0/1), their
interaction, and control covariates. Then we use Margins.did() to get,
on the *probability* scale:

  * 4 cell predictions      P(X | group, Y)
  * 2 simple effects        P(X|B,Y) - P(X|A,Y)   at each Y
  * 1 DiD                   (simple effect at Y=1) - (simple effect at Y=0)

All with delta-method standard errors and CIs.

The DiD here is the "does the A-vs-B gap depend on Y?" question.  It is
NOT the coefficient on group×Y (that's on the log-odds scale); on the
probability scale you have to go through the inverse link — which is
exactly what Margins.did() does.
"""
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

from smmargins import Margins

pd.options.display.width = 140
pd.options.display.float_format = "{: .4f}".format

# ---------------------------------------------------------------------------
# Simulate patient-level data
# ---------------------------------------------------------------------------
rng = np.random.default_rng(42)
N = 6_000
df = pd.DataFrame({
    "group":        rng.choice(["A", "B"], N, p=[0.55, 0.45]),
    "preexist_Y":   rng.integers(0, 2, N),             # 0 = no Y, 1 = has Y
    "age":          rng.normal(55, 15, N).clip(18, 95),
    "female":       rng.integers(0, 2, N),
})

# True data-generating process:
#   * baseline rate of X depends on age, sex, and Y
#   * group B has a modest additive bump in the log-odds
#   * the group effect is AMPLIFIED among patients with preexisting Y
#     (this is the thing we want to detect)
eta = (
    -3.5
    + 0.04 * df["age"]
    - 0.3 * df["female"]
    + 0.5 * (df["group"] == "B")
    + 1.1 * df["preexist_Y"]
    + 0.8 * (df["group"] == "B") * df["preexist_Y"]   # interaction
)
df["condition_X"] = (rng.uniform(0, 1, N) < 1 / (1 + np.exp(-eta))).astype(int)

print("Raw sample rates of condition X by cell:")
print(df.groupby(["group", "preexist_Y"])["condition_X"].mean().round(4))
print()

# ---------------------------------------------------------------------------
# Fit the logit with the group × preexist_Y interaction + controls
# ---------------------------------------------------------------------------
fit = smf.logit(
    "condition_X ~ C(group) + preexist_Y + C(group):preexist_Y + age + female",
    data=df,
).fit(disp=False)

print("=" * 84)
print("Logit model (coefficients are on the log-odds scale)")
print("=" * 84)
print(fit.summary().tables[1])
print()

# ---------------------------------------------------------------------------
# DiD on the *probability* (response) scale — what the clinical question asks
# ---------------------------------------------------------------------------
# Margins(fit) uses analytic outer Jacobians via family.link.inverse_deriv
# when available (Logit qualifies), falling back to central finite
# differences otherwise. did() reuses predict()'s machinery, so it
# inherits the analytic path automatically. Set Margins(fit,
# analytic=False) to force FD if you ever want to cross-check.
M = Margins(fit)
did = M.did("group", "preexist_Y",
            group_levels=["A", "B"],
            condition_levels=[0, 1])

print("=" * 84)
print("DiD on the probability scale, averaged over age and sex distribution")
print("=" * 84)
print(did)

# ---------------------------------------------------------------------------
# Interpretation
# ---------------------------------------------------------------------------
pA0 = did.cells.estimate[0]   # group=A, Y=0
pA1 = did.cells.estimate[1]   # group=A, Y=1
pB0 = did.cells.estimate[2]   # group=B, Y=0
pB1 = did.cells.estimate[3]   # group=B, Y=1
se_simple_Y0 = did.simple_effects.se[0]
se_simple_Y1 = did.simple_effects.se[1]
did_est, did_se = did.did.estimate[0], did.did.se[0]

print()
print("=" * 84)
print("Plain-language summary")
print("=" * 84)
print(f"Condition X rate, group A, no preexisting Y : {pA0:.3%}")
print(f"Condition X rate, group A, with Y           : {pA1:.3%}")
print(f"Condition X rate, group B, no preexisting Y : {pB0:.3%}")
print(f"Condition X rate, group B, with Y           : {pB1:.3%}")
print()
print(f"Rate difference (B - A) among NO-Y patients  : "
      f"{(pB0 - pA0):+.3%}  (SE {se_simple_Y0:.3%})")
print(f"Rate difference (B - A) among WITH-Y patients: "
      f"{(pB1 - pA1):+.3%}  (SE {se_simple_Y1:.3%})")
print()
print(f"Difference-in-differences                   : "
      f"{did_est:+.3%}  (SE {did_se:.3%})")
print(f"  -> the B-vs-A gap is {abs(did_est):.3%} larger among patients "
      f"with preexisting Y.")
print(f"  -> 95% CI: ({did.did.ci_lower[0]:+.3%}, {did.did.ci_upper[0]:+.3%})")
print(f"  -> p-value: {did.did.pvalues[0]:.4g}")

# ---------------------------------------------------------------------------
# Sensitivity: DiD at a specific patient profile (e.g. 60-year-old male)
# ---------------------------------------------------------------------------
print()
print("=" * 84)
print("DiD at a specific profile: 60-year-old male")
print("=" * 84)
did_profile = M.did(
    "group", "preexist_Y",
    group_levels=["A", "B"], condition_levels=[0, 1],
    atexog={"age": 60, "female": 0},
)
print(did_profile.did)

# ---------------------------------------------------------------------------
# Bonus: plottable table of cell predictions with CIs
# ---------------------------------------------------------------------------
print()
print("=" * 84)
print("Cells with 95% CIs (suitable for a plot)")
print("=" * 84)
tbl = did.cells.summary().copy()
print(tbl)

# If you wanted to plot:
#   import matplotlib.pyplot as plt
#   fig, ax = plt.subplots()
#   for g in ["A", "B"]:
#       sub = tbl[tbl.index.str.contains(f"group={g}")]
#       ax.errorbar([0, 1],
#                   sub["prediction"].values,
#                   yerr=(sub["prediction"] - sub["[95% CI lo]"]).values,
#                   marker="o", label=f"group {g}", capsize=4)
#   ax.set_xticks([0, 1]); ax.set_xticklabels(["no Y", "with Y"])
#   ax.set_ylabel("P(condition X)"); ax.legend(); plt.show()