Bayesian marginal effects — posterior draws via from_posterior

Every inference path in pymargins ultimately needs a distribution over the parameter vector β. The delta method assumes that distribution is normal with covariance vcov; Krinsky–Robb simulates draws from that same normal; the bootstrap resamples it. But if you fit the model in a Bayesian framework — PyMC, NumPyro, Stan, bambi — you already have the posterior as a bank of MCMC draws, and there is no reason to throw it away and re-approximate it as normal.

Margins.from_posterior takes those draws directly. It runs the simulation inference path, but instead of drawing β from N(β̂, V̂) it uses your draws as-is. The marginal effect is computed once per draw, and the spread of those values is the posterior of the marginal effect — a genuine credible interval on the outcome scale, with no delta-method linearization anywhere.

This demo uses the Mroz female-labor-force data and a logit. To keep the documentation build dependency-free, we generate the posterior draws with a Laplace approximation — a multivariate normal centred at the MLE with the model’s estimated covariance. Real MCMC draws from PyMC or NumPyro are an array of exactly the same shape (n_draws, n_params) and plug into the same call; the only thing that changes is where the draws come from. See the closing note for the drop-in PyMC/NumPyro snippet.

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 linearmodels.datasets import mroz

from pymargins import Margins

cols = ["inlf", "nwifeinc", "educ", "exper", "age", "kidslt6", "kidsge6"]
df = mroz.load()[cols].dropna().copy()
df["expersq"] = df["exper"] ** 2

fit = smf.glm(
    "inlf ~ nwifeinc + educ + exper + expersq + age + kidslt6 + kidsge6",
    data=df,
    family=sm.families.Binomial(),
).fit()
print(fit.summary().tables[1])

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.4255      0.860      0.495      0.621      -1.261       2.112
nwifeinc      -0.0213      0.008     -2.535      0.011      -0.038      -0.005
educ           0.2212      0.043      5.091      0.000       0.136       0.306
exper          0.2059      0.032      6.422      0.000       0.143       0.269
expersq       -0.0032      0.001     -3.104      0.002      -0.005      -0.001
age           -0.0880      0.015     -6.040      0.000      -0.117      -0.059
kidslt6       -1.4434      0.204     -7.090      0.000      -1.842      -1.044
kidsge6        0.0601      0.075      0.804      0.422      -0.086       0.207
==============================================================================

1. Build the posterior draw bank

The Laplace posterior is N(β̂, V̂). We draw 4 000 samples; the columns are in the same order as fit.params, which is exactly the order the adapter expects.

beta_hat = np.asarray(fit.params)
V_hat = np.asarray(fit.cov_params())

rng = np.random.default_rng(0)
posterior_draws = rng.multivariate_normal(beta_hat, V_hat, size=4000)
print("draw bank shape:", posterior_draws.shape)  # (n_draws, n_params)

draw bank shape: (4000, 8)

2. Open a posterior session

from_posterior defaults method="simulation" and n_sim to the number of draws. The point estimate defaults to the posterior mean; pass point_estimate= to report a different summary (e.g. the MLE or the posterior median).

m_bayes = Margins.from_posterior(fit, posterior_draws, at="overall")
print("method:", m_bayes.method, " draws:", m_bayes.n_sim)

method: simulation  draws: 4000

3. AME credible interval vs delta-method interval

Ask the posterior session for the average marginal effect of education on participation. The interval it returns is a 95 % credible interval — the 2.5th and 97.5th percentiles of the AME across the posterior draws. We put it side by side with the ordinary delta-method confidence interval from a standard session.

m_delta = Margins.linear_scale(fit, at="overall")

ame_bayes = m_bayes.dydx("educ")
ame_delta = m_delta.dydx("educ")

print("Posterior (credible interval):")
print(ame_bayes.summary())
print("\nDelta method (confidence interval):")
print(ame_delta.summary())

Posterior (credible interval):
===========================================================
          Margins Result (simulation, level=0.95)          
===========================================================
      estimate  std err  statistic  P>|z|  [95% Conf. Int.]
-----------------------------------------------------------
educ    0.0395   0.0072     0.0395  0.000    0.0247, 0.0528
===========================================================

n = 753
κ: 0.055

Delta method (confidence interval):
========================================================
           Margins Result (delta, level=0.95)           
========================================================
      estimate  std err       z  P>|z|  [95% Conf. Int.]
--------------------------------------------------------
educ    0.0395   0.0073  5.4144  0.000    0.0252, 0.0538
========================================================

n = 753
κ: 0.054
Delta-vs-sim disagreement: 2.267%

Under a flat prior the Laplace posterior is, by construction, the same normal the delta method assumes — so the two intervals here agree to within Monte-Carlo error. That agreement is the point: it shows the posterior path is calibrated against the analytic one in the case where they should coincide. The value of from_posterior shows up when they should not coincide — an informative or skewed prior, a small sample where the posterior is non-normal, or a nonlinear estimand whose posterior is asymmetric. There the credible interval follows the true posterior shape and the delta interval cannot.

4. The full posterior of an effect, not just an interval

Because every draw produces a complete marginal effect, you are not limited to a symmetric interval — you can report any posterior summary. Here is the AME of an additional young child (a discrete contrast) as a posterior, with its probability of being negative.

from pymargins import pairwise

scen, w = pairwise("kidslt6", [1, 0])
child_effect = m_bayes.contrasts(scenarios=scen, contrasts=w)
print(child_effect.summary())

draws = np.asarray(child_effect.draws)
print(f"\nPosterior mean : {draws.mean():.4f}")
print(f"P(effect < 0)  : {(draws < 0).mean():.3f}")
print(f"90% credible   : [{np.percentile(draws, 5):.4f}, "
      f"{np.percentile(draws, 95):.4f}]")

================================================================
            Margins Result (simulation, level=0.95)             
================================================================
           estimate  std err  statistic  P>|z|  [95% Conf. Int.]
----------------------------------------------------------------
kidslt6=1   -0.2697   0.0347    -0.2697  0.000  -0.3321, -0.1982
================================================================

n = 753
κ: 0.057

Posterior mean : -0.2668
P(effect < 0)  : 1.000
90% credible   : [-0.3228, -0.2100]

P(effect < 0) is a direct posterior probability — the kind of statement a Bayesian analysis is supposed to deliver and that a p-value cannot. It comes for free once the estimand is evaluated on each draw.

Plugging in real MCMC draws

Nothing above is specific to the Laplace approximation. If you fit the same logit in PyMC or NumPyro, you obtain a posterior array of shape (n_draws, n_params) and pass it to the same constructor:

# PyMC
import pymc as pm
with pm.Model() as model:
    # ... priors and likelihood for the same design matrix ...
    idata = pm.sample()
# stack chains × draws into (n_draws, n_params), columns in fit.params order
posterior_draws = idata.posterior[coef_names].to_array() \
    .transpose("chain", "draw", "variable") \
    .values.reshape(-1, len(coef_names))

m_bayes = Margins.from_posterior(fit, posterior_draws, at="overall")
m_bayes.dydx("educ")   # AME with a real posterior credible interval

The only requirements are that the draw columns line up with the adapter’s coefficient order (the order of fit.params) and that fit is the same specification you sampled. Everything downstream — predict, dydx, contrasts, evaluate, subgroups, plotting — behaves identically to any other session.

Where to next