Replication scripts (archive)¶

The notebook-style dataset demos under Demos — end-to-end analyses are the recommended entry point. The script-style demos archived here ship as standalone .py files in demo/ and are kept for reference and exact numerical reproduction against the Williams (2012) paper.

Williams (2012) — Adjusted predictions and marginal effects ¶

demo/williams_2012_demo.py replicates, on a simulated NHANES-like dataset, every core analysis from Richard Williams’ Stata Journal paper:

Logit model with factor variables
Adjusted Predictions at the Means (APM)
Average Adjusted Predictions (AAP)
Predictions at Representative Values (APR)
Marginal Effects at the Means (MEM) — continuous
Average Marginal Effects (AME) — continuous
Discrete changes for dummy variables
Marginal Effects at Representative Values (MER)
OLS model for comparison

"""
Demo: Replicating analyses from Richard Williams (2012),
"Using the Margins Command to Estimate and Interpret Adjusted Predictions
and Marginal Effects", Stata Journal 12(2): 308-331.

DOI: 10.1177/1536867X1201200209

This script generates synthetic NHANES-like data and demonstrates the
core analyses from the paper using statsmodels + pymargins.

Analyses covered:
  1. Logit model with factor variables
  2. Adjusted Predictions at the Means (APM)
  3. Average Adjusted Predictions (AAP)
  4. Predictions at Representative Values
  5. Marginal Effects at the Means (MEM) — continuous
  6. Average Marginal Effects (AME) — continuous
  7. Discrete changes for dummy variables
  8. Marginal Effects at Representative Values (MER)
  9. OLS model for comparison
"""

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

from pymargins import Margins

# ---------------------------------------------------------------------------
# 0. Generate synthetic NHANES-like data
# ---------------------------------------------------------------------------


def make_data(n=5000, seed=42):
    """Synthetic health-survey data mimicking NHANES II structure."""
    rng = np.random.default_rng(seed)

    # Demographics
    female = rng.binomial(1, 0.52, size=n)
    black = rng.binomial(1, 0.11, size=n)
    age = rng.integers(20, 75, size=n)

    # Age groups (1=20-29, 2=30-39, ..., 6=70+)
    agegrp = pd.cut(
        age, bins=[19, 29, 39, 49, 59, 69, 100], labels=[1, 2, 3, 4, 5, 6]
    ).astype(int)

    # BMI (correlated with age and sex)
    bmi = 22 + 0.15 * age + 1.5 * female + rng.normal(0, 4, size=n)
    bmi = np.clip(bmi, 15, 50)

    # Diabetes risk (logit link)
    lp = (
        -4.0
        + 0.55 * black
        + 0.10 * female
        + 0.06 * age
        + 0.03 * bmi
        + 0.5 * (agegrp == 2).astype(float)
        + 0.9 * (agegrp == 3).astype(float)
        + 1.4 * (agegrp == 4).astype(float)
        + 2.0 * (agegrp == 5).astype(float)
        + 2.6 * (agegrp == 6).astype(float)
    )
    diabetes = rng.binomial(1, 1 / (1 + np.exp(-lp)))

    # Continuous outcome (e.g., blood pressure) for OLS demo
    bp = (
        110
        + 0.4 * age
        + 2.5 * black
        + 1.2 * female
        + 0.5 * bmi
        + rng.normal(0, 8, size=n)
    )

    df = pd.DataFrame(
        {
            "diabetes": diabetes,
            "bp": bp,
            "black": black,
            "female": female,
            "age": age,
            "agegrp": agegrp,
            "bmi": bmi,
        }
    )
    return df


if __name__ == "__main__":
    df = make_data(n=5000, seed=42)
    print("=" * 70)
    print(f"Synthetic NHANES-like data (n = {len(df):,})")
    print("=" * 70)
    print(df.describe())
    print()

    # =====================================================================
    # 1. Logit model with factor variables (Williams §2)
    # =====================================================================
    print("=" * 70)
    print("1. LOGIT MODEL: diabetes ~ C(black) + C(female) + C(agegrp) + bmi + age")
    print("=" * 70)

    fit_logit = smf.glm(
        "diabetes ~ C(black) + C(female) + C(agegrp) + bmi + age",
        data=df,
        family=sm.families.Binomial(),
    ).fit()
    print(fit_logit.summary())
    print()

    # =====================================================================
    # 2. Adjusted Predictions at the Means (APM)  (Williams §3.1)
    # =====================================================================
    # In Stata: margins C(agegrp), atmeans
    # All covariates held at their typical values (continuous = median,
    # categorical = mode).  For models with factor variables this is the
    # safest pymargins equivalent because patsy requires integer category
    # codes; setting a dummy to its mean proportion (0.107) is not a valid
    # level for patsy even though it is mathematically correct.
    print("=" * 70)
    print("2. ADJUSTED PREDICTIONS AT THE MEANS (APM)")
    print("   Stata equivalent: margins agegrp, atmeans")
    print("=" * 70)

    m_apm = Margins.log_scale(fit_logit, at="typical")
    apm = m_apm.predict(
        atexog={"agegrp": list(range(1, 7))},
    )
    print(apm.summary(stars=True))
    print()

    # =====================================================================
    # 3. Average Adjusted Predictions (AAP)  (Williams §3.2)
    # =====================================================================
    # In Stata: margins C(agegrp)
    # Predictions averaged over the observed distribution of other covariates.
    print("=" * 70)
    print("3. AVERAGE ADJUSTED PREDICTIONS (AAP)")
    print("   Stata equivalent: margins agegrp")
    print("=" * 70)

    m_aap = Margins.log_scale(fit_logit, at="overall")
    aap = m_aap.predict(
        atexog={"agegrp": list(range(1, 7))},
    )
    print(aap.summary(stars=True))
    print()

    # =====================================================================
    # 4. Predictions at Representative Values  (Williams §3.3)
    # =====================================================================
    # In Stata: margins, at(age=(20 50 70)) atmeans
    # Evaluate predictions for specific ages, holding others at means.
    print("=" * 70)
    print("4. PREDICTIONS AT REPRESENTATIVE VALUES")
    print("   Stata equivalent: margins, at(age=(20 50 70)) atmeans")
    print("=" * 70)

    m_repr = Margins.log_scale(fit_logit, at="typical")
    repr_pred = m_repr.predict(
        atexog={"age": [20, 50, 70]},
    )
    print(repr_pred.summary(stars=True))
    print()

    # =====================================================================
    # 5. Marginal Effects at the Means (MEM) — continuous  (Williams §4.1)
    # =====================================================================
    # In Stata: margins, dydx(age) atmeans
    print("=" * 70)
    print("5. MARGINAL EFFECTS AT THE MEANS (MEM) — age")
    print("   Stata equivalent: margins, dydx(age) atmeans")
    print("=" * 70)

    m_mem = Margins.log_scale(fit_logit, at="typical")
    mem_age = m_mem.dydx("age")
    print(mem_age.summary(stars=True))
    print()

    # =====================================================================
    # 6. Average Marginal Effects (AME) — continuous  (Williams §4.2)
    # =====================================================================
    # In Stata: margins, dydx(age)
    print("=" * 70)
    print("6. AVERAGE MARGINAL EFFECTS (AME) — age")
    print("   Stata equivalent: margins, dydx(age)")
    print("=" * 70)

    m_ame = Margins.log_scale(fit_logit, at="overall")
    ame_age = m_ame.dydx("age")
    print(ame_age.summary(stars=True))
    print()

    # =====================================================================
    # 7. Discrete changes for dummy variables  (Williams §4.3)
    # =====================================================================
    # For binary/discrete variables, Stata margins computes the discrete
    # change (difference in predicted probability when flipping the dummy
    # from 0 to 1), not a derivative.
    #
    # In Stata: margins, dydx(black)
    # In pymargins: use contrasts() with two scenarios.
    print("=" * 70)
    print("7. DISCRETE CHANGE (MARGINAL EFFECT) — black")
    print("   Stata equivalent: margins, dydx(black)")
    print("=" * 70)

    m_disc = Margins.log_scale(fit_logit, at="overall")
    disc_black = m_disc.contrasts(
        scenarios=[
            {"atexog": {"black": 1}, "label": "black=1"},
            {"atexog": {"black": 0}, "label": "black=0"},
        ],
        contrasts=[+1, -1],
    )
    print(disc_black.summary(stars=True))
    print()

    # Discrete change at the means
    print("-" * 70)
    print("7b. DISCRETE CHANGE AT THE MEANS — female")
    print("    Stata equivalent: margins, dydx(female) atmeans")
    print("-" * 70)

    m_disc_mem = Margins.log_scale(fit_logit, at="typical")
    disc_female = m_disc_mem.contrasts(
        scenarios=[
            {"atexog": {"female": 1}, "label": "female=1"},
            {"atexog": {"female": 0}, "label": "female=0"},
        ],
        contrasts=[+1, -1],
    )
    print(disc_female.summary(stars=True))
    print()

    # =====================================================================
    # 8. Marginal Effects at Representative Values (MER)  (Williams §4.4)
    # =====================================================================
    # In Stata: margins, dydx(black) at(female=(0 1))
    # Compute the discrete change for black, separately for females and males.
    print("=" * 70)
    print("8. MARGINAL EFFECTS AT REPRESENTATIVE VALUES (MER)")
    print("   Stata equivalent: margins, dydx(black) at(female=(0 1))")
    print("=" * 70)

    m_mer = Margins.log_scale(fit_logit, at="overall")
    mer_female0 = m_mer.contrasts(
        scenarios=[
            {"atexog": {"black": 1, "female": 0}, "label": "black=1, female=0"},
            {"atexog": {"black": 0, "female": 0}, "label": "black=0, female=0"},
        ],
        contrasts=[+1, -1],
    )
    mer_female1 = m_mer.contrasts(
        scenarios=[
            {"atexog": {"black": 1, "female": 1}, "label": "black=1, female=1"},
            {"atexog": {"black": 0, "female": 1}, "label": "black=0, female=1"},
        ],
        contrasts=[+1, -1],
    )
    print("Effect of black for MALES (female=0):")
    print(mer_female0.summary(stars=True))
    print()
    print("Effect of black for FEMALES (female=1):")
    print(mer_female1.summary(stars=True))
    print()

    # =====================================================================
    # 9. OLS model — compare with logit  (Williams §5)
    # =====================================================================
    print("=" * 70)
    print("9. OLS MODEL: bp ~ C(black) + C(female) + age + bmi")
    print("=" * 70)

    fit_ols = smf.ols(
        "bp ~ C(black) + C(female) + age + bmi",
        data=df,
    ).fit()
    print(fit_ols.summary())
    print()

    # MEM for OLS (linear model: MEM = AME = coefficient)
    m_ols = Margins.linear_scale(fit_ols, at="typical")
    mem_ols_age = m_ols.dydx("age")
    print("MEM of age in OLS model (should match coefficient):")
    print(mem_ols_age.summary(stars=True))
    print()

    # AME for OLS (same as MEM in linear models)
    m_ols_ame = Margins.linear_scale(fit_ols, at="overall")
    ame_ols_age = m_ols_ame.dydx("age")
    print("AME of age in OLS model:")
    print(ame_ols_age.summary(stars=True))
    print()

    # =====================================================================
    # 10. LaTeX / HTML export demos
    # =====================================================================
    print("=" * 70)
    print("10. EXPORT EXAMPLES")
    print("=" * 70)

    print("--- LaTeX output for AME of age ---")
    print(ame_age.to_latex(stars=True, caption="AME of Age on Diabetes Risk"))
    print()

    print("--- HTML output for AME of age ---")
    print(ame_age.to_html(stars=True, caption="AME of Age on Diabetes Risk"))
    print()

    print("Demo complete.")

Williams (2012) — Inference scales ¶

demo/williams_2012_demo_scales.py reruns the same analyses on the log, logit, and lift inference scales, demonstrating how a single session locks in \((\phi, \phi^{-1})\) for the whole audit trail and how the κ diagnostic flags when the chosen scale is poorly linearized.

"""
Demo: Scales and contrasts — RR, direct ratio, and lift.

Uses the same synthetic NHANES-like data as williams_2012_demo.py.
Demonstrates three ways to quantify the effect of a binary covariate:

  1. Risk Ratio (RR)      — log_scale:  exp(log(p1) - log(p0))
  2. Direct Ratio         — linear_scale + evaluate:  p1 / p0
  3. Lift (RR - 1)        — log_scale result minus 1

Reference implementations:
  - StatsModels: manual point estimates
  - R marginaleffects: delta-method SEs and CIs
"""

import jax

jax.config.update("jax_enable_x64", True)

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

from pymargins import Margins


def make_data(n=5000, seed=42):
    """Synthetic health-survey data mimicking NHANES II structure."""
    rng = np.random.default_rng(seed)
    female = rng.binomial(1, 0.52, size=n)
    black = rng.binomial(1, 0.11, size=n)
    age = rng.integers(20, 75, size=n)
    agegrp = pd.cut(
        age, bins=[19, 29, 39, 49, 59, 69, 100], labels=[1, 2, 3, 4, 5, 6]
    ).astype(int)
    bmi = 22 + 0.15 * age + 1.5 * female + rng.normal(0, 4, size=n)
    bmi = np.clip(bmi, 15, 50)
    lp = (
        -4.0
        + 0.55 * black
        + 0.10 * female
        + 0.06 * age
        + 0.03 * bmi
        + 0.5 * (agegrp == 2).astype(float)
        + 0.9 * (agegrp == 3).astype(float)
        + 1.4 * (agegrp == 4).astype(float)
        + 2.0 * (agegrp == 5).astype(float)
        + 2.6 * (agegrp == 6).astype(float)
    )
    diabetes = rng.binomial(1, 1 / (1 + np.exp(-lp)))
    bp = (
        110
        + 0.4 * age
        + 2.5 * black
        + 1.2 * female
        + 0.5 * bmi
        + rng.normal(0, 8, size=n)
    )
    return pd.DataFrame(
        {
            "diabetes": diabetes,
            "bp": bp,
            "black": black,
            "female": female,
            "age": age,
            "agegrp": agegrp,
            "bmi": bmi,
        }
    )


def _print_section(title):
    print("=" * 70)
    print(title)
    print("=" * 70)


if __name__ == "__main__":
    df = make_data(n=5000, seed=42)

    fit_logit = smf.glm(
        "diabetes ~ C(black) + C(female) + C(agegrp) + bmi + age",
        data=df,
        family=sm.families.Binomial(),
    ).fit()

    # =====================================================================
    # 1. Risk Ratio via log_scale
    # =====================================================================
    _print_section("1. RISK RATIO — log_scale")
    print("   Inference: log(p1) - log(p0)   |   Reported: exp(...)")
    print()

    m_rr = Margins.log_scale(fit_logit, at="overall", kappa_threshold=float("inf"))
    rr_black = m_rr.contrasts(
        scenarios=[
            {"atexog": {"black": 1}, "label": "black=1"},
            {"atexog": {"black": 0}, "label": "black=0"},
        ],
        contrasts=[+1, -1],
    )
    print(rr_black.summary(stars=True))
    print()

    rr_female = m_rr.contrasts(
        scenarios=[
            {"atexog": {"female": 1}, "label": "female=1"},
            {"atexog": {"female": 0}, "label": "female=0"},
        ],
        contrasts=[+1, -1],
    )
    print(rr_female.summary(stars=True))
    print()

    # =====================================================================
    # 2. Direct Ratio via linear_scale + evaluate()
    # =====================================================================
    _print_section("2. DIRECT RATIO — linear_scale + evaluate(p[0]/p[1])")
    print("   Inference: p1 / p0 directly   |   Delta-method SE on ratio scale")
    print()

    m_ratio = Margins.linear_scale(
        fit_logit, at="overall", kappa_threshold=float("inf")
    )
    ratio_black = m_ratio.evaluate(
        scenarios=[
            {"atexog": {"black": 1}, "label": "black=1"},
            {"atexog": {"black": 0}, "label": "black=0"},
        ],
        compose=lambda p: p[0] / p[1],
    )
    print(ratio_black.summary(stars=True))
    print()

    ratio_female = m_ratio.evaluate(
        scenarios=[
            {"atexog": {"female": 1}, "label": "female=1"},
            {"atexog": {"female": 0}, "label": "female=0"},
        ],
        compose=lambda p: p[0] / p[1],
    )
    print(ratio_female.summary(stars=True))
    print()

    # =====================================================================
    # 3. True Lift = RR - 1
    # =====================================================================
    _print_section("3. TRUE LIFT (RR - 1)")
    print("   Computed from log_scale RR:  lift = RR - 1")
    print()

    lift_black_est = float(rr_black.estimate) - 1.0
    lift_black_ci = (
        float(rr_black.conf_int_lower) - 1.0,
        float(rr_black.conf_int_upper) - 1.0,
    )
    print(
        f"black:  lift = {lift_black_est:.4f}   CI = ({lift_black_ci[0]:.4f}, {lift_black_ci[1]:.4f})"
    )

    lift_female_est = float(rr_female.estimate) - 1.0
    lift_female_ci = (
        float(rr_female.conf_int_lower) - 1.0,
        float(rr_female.conf_int_upper) - 1.0,
    )
    print(
        f"female: lift = {lift_female_est:.4f}   CI = ({lift_female_ci[0]:.4f}, {lift_female_ci[1]:.4f})"
    )
    print()

    # =====================================================================
    # Compact reference table
    # =====================================================================
    _print_section("REFERENCE TABLE — All Scales")
    print(f"{'Scale':<20} {'black':>12} {'female':>12}")
    print("-" * 46)
    print(
        f"{'RR (log_scale)':<20} {float(rr_black.estimate):>12.4f} {float(rr_female.estimate):>12.4f}"
    )
    print(
        f"{'Direct ratio':<20} {float(ratio_black.estimate):>12.4f} {float(ratio_female.estimate):>12.4f}"
    )
    print(f"{'True lift (RR-1)':<20} {lift_black_est:>12.4f} {lift_female_est:>12.4f}")
    print()
    print("Demo complete.")

Cross-checks ¶

The companion files demo/williams_2012_demo_statsmodels.py, demo/williams_2012_demo_marginaleffects_r.R, and demo/williams_2012_demo_scales_marginaleffects_r.R reproduce the same numbers using statsmodels.get_margeff and the R marginaleffects package. These exist to pin pymargins to externally agreed answers at the precision both tools agree on.

Replication scripts (archive)¶

Williams (2012) — Adjusted predictions and marginal effects ¶

Williams (2012) — Inference scales ¶

Cross-checks ¶

pymargins

Navigation

Related Topics

Replication scripts (archive)¶

Williams (2012) — Adjusted predictions and marginal effects¶

Williams (2012) — Inference scales¶

Cross-checks¶

Williams (2012) — Adjusted predictions and marginal effects ¶

Williams (2012) — Inference scales ¶

Cross-checks ¶