pymmeans

Estimated marginal means (EMMs) for Python — a native implementation of R's emmeans package, with no R dependency.

Status: Beta (v0.2.7). API stable across the OLS / GLM / MixedLM / GEE / Cox / Beta surface; 90/100 strict parity with R emmeans (94/100 if partially-supported items count — see docs/r_parity_matrix.md) validated against tests/r_reference/ CSVs at atol=1e-4 (and tighter — see vs-r.md). 352 unit tests on the public surface (pytest), 54% line coverage. Minor API polish still possible; the numerical surface is frozen.

Install

pip install pymmeans              # from PyPI
pip install "pymmeans[plot]"      # add matplotlib for plot() / emmip()
pip install "pymmeans[tutorial]"  # add pysofra + jupyter for the showcase notebook

For local development:

git clone https://github.com/jturner-uofl/pymmeans.git
cd pymmeans
uv venv && uv pip install -e ".[dev,plot]"

Quickstart

Fit a model, then ask for marginal means, pairwise comparisons, and a plot:

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from pymmeans import emmeans, pairs, plot

rng = np.random.default_rng(0)
df = pd.DataFrame({
    "fertilizer": np.tile(np.repeat(["lo", "med", "hi"], 30), 2),
    "sunlight":   np.repeat(["shade", "sun"], 90),
})
df["growth"] = (
    df["fertilizer"].map({"lo": 0.5, "med": 1.2, "hi": 1.8})
    + (df["sunlight"] == "sun") * 0.4
    + rng.normal(0, 0.3, 180)
)
model = smf.ols("growth ~ fertilizer * sunlight", data=df).fit()

print(emmeans(model, "fertilizer"))
# Output (with the seed above):
#   fertilizer    emmean       SE     df  lower_cl  upper_cl
# 0         hi   2.0071   0.0376  174.0    1.9328    2.0813
# 1         lo   0.6811   0.0376  174.0    0.6069    0.7554
# 2        med   1.4399   0.0376  174.0    1.3657    1.5142

print(pairs(emmeans(model, "fertilizer")))
#   contrast  estimate      SE     df  t_ratio  p_value
# 0  hi - lo    1.3259  0.0532  174.0  24.9230      0.0
# 1 hi - med    0.5671  0.0532  174.0  10.6596      0.0
# 2 lo - med   -0.7588  0.0532  174.0 -14.2633      0.0
# (p_value=0.0 is pandas' default display of underflowed p < ~1e-15;
#  use `summary(pairs(...))` or `as_r_frame(...)` for "<.0001" format.)

ax = plot(emmeans(model, "fertilizer", by="sunlight"))

Interaction plot via emmip:

Showcase notebook

A full walkthrough of every analytical surface — EMMs, contrasts, multiplicity adjustments, mixed-model Kenward-Roger / Satterthwaite df, parametric-bootstrap LRT, and the ML adapter — is at examples/pymmeans_showcase.ipynb with a self-contained HTML render at examples/pymmeans_showcase.html.

What's in v0.1

Beyond R parity — ML adapter example

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from pymmeans import from_predict, ml_emmeans, ml_pairs, bootstrap_ci

# 0) Toy training data: 3 treatments × 2 sites + 2 numeric covariates
rng = np.random.default_rng(0)
df = pd.DataFrame({
    "treatment": np.repeat(["A", "B", "C"], 50),
    "site":      np.tile(np.repeat(["north", "south"], 25), 3),
    "age":       rng.uniform(20, 70, 150),
    "dose":      rng.uniform(0.0, 1.0, 150),
})
df["y"] = (
    df["treatment"].map({"A": 0.8, "B": 1.2, "C": 1.7})
    + 0.01 * df["age"]
    + 0.5 * df["dose"]
    + rng.normal(0, 0.1, 150)
)

# Pin a stable feature schema so sub-grid predictions always
# present the same columns sklearn saw at fit time.
_FEATURE_COLS = pd.get_dummies(
    df[["treatment", "site", "age", "dose"]],
).columns.tolist()

def featurize(d):
    return (
        pd.get_dummies(d[["treatment", "site", "age", "dose"]])
        .reindex(columns=_FEATURE_COLS, fill_value=0)
    )

X_train, y_train = featurize(df), df["y"]

# 1) Train any sklearn-style model
rf = RandomForestRegressor(random_state=0).fit(X_train, y_train)

# 2) Wrap as a pymmeans target via the predict callable
info = from_predict(
    predict_fn=lambda d: rf.predict(featurize(d)),
    data=df,
    factors=["treatment", "site"],
    numerics=["age", "dose"],
    refit_fn=lambda sample: (lambda fitted: lambda d: fitted.predict(featurize(d)))(
        RandomForestRegressor(random_state=0)
        .fit(featurize(sample), sample["y"])
    ),
)

# 3) Marginal means via prediction-surface averaging (g-computation)
em = ml_emmeans(info, "treatment")
#   treatment  emmean
# 0         A   1.531
# 1         B   1.939
# 2         C   2.464

# 4) Pairwise contrasts
ml_pairs(em)

# 5) Case-bootstrap CIs (refits on each resample for proper variance)
em_with_ci = bootstrap_ci(em, n_samples=500, kind="case", seed=0)
#   treatment  emmean     SE  lower_cl  upper_cl
# 0         A   1.531  0.025     1.483     1.582
# 1         B   1.939  0.022     1.899     1.983
# 2         C   2.464  0.028     2.402     2.509

R emmeans is restricted to models with a tractable β + V representation. pymmeans's prediction-surface averaging extends the marginal-effects workflow to any ML model with a .predict() method, sharing the same pairs / contrast / bootstrap_ci / effect_size / summary machinery as the linear-model path. The closest R-ecosystem analogue is marginaleffects, which also supports ML models; pymmeans brings the same population-average prediction workflow to Python without leaving the EMM toolbox.

Validation

352 public-surface unit tests pass (a small number of dependency-gated tests skip when optional packages are absent). 18 R-emmeans reference fits are cross-validated to atol=1e-4 (warpbreaks, pigs, ToothGrowth, InsectSprays, neuralgia binomial GLM response-scale, plus the broader reference suite). See tests/test_vs_r.py, tests/r_reference/, and docs/r_parity_matrix.md for the feature-level parity inventory.

Performance vs R

See docs/PERFORMANCE_REPORT.md for full numbers. Highlights:

pymmeans wins on the EMM side via analytic marginalization (computing L_marg as a Kronecker product of marginal factor codings, no grid materialization). Tukey uses a numerically correct hybrid of generalized Gauss-Laguerre (df < 400) and adaptive Gauss-Legendre (df ≥ 400, where SciPy's roots_genlaguerre overflows) — beating R at large k where Tukey is most painful, slightly slower at moderate k.

Documentation

A mkdocs-material site builds from docs/ (run mkdocs serve after pip install -e ".[docs]").

• docs/getting-started.md — worked examples
• docs/PERFORMANCE_REPORT.md — benchmark numbers
• docs/vs-r.md — R reference comparison and validation
• docs/r_parity_matrix.md — function-by-function map between pymmeans and R emmeans / lsmeans / pbkrtest
• docs/v0_2_roadmap.md — what's planned for the next release
• examples/pymmeans_showcase.ipynb — the full showcase notebook
• examples/basic_ols.py, examples/glm_logistic.py — runnable demos

License

GPL-3.0-or-later, matching the R emmeans package.