Code
import polars as pl
from polars import col
import seaborn as sns
from bossanova import model, load_datasetMarginal Effects & the Mixing Board
In the previous notebooks we built a toolkit for working with regression models — fitting, diagnosing, and running inference on parameter estimates. We asked: “How certain are we in \(\hat{\beta}\)?”
We learned that by default, the parameters reflect how the GLM can see and this is controlled by the design matrix. Remember for categorical variables (switches) we have a choice of a variety of different encoding schemes to choose from (dummy, sum, poly). This means that if you want the model to “see” that data differently you need to change the design matrix (i.e. use a different coding scheme, center variables, etc).
In this notebook we shift focus from the model’s internal wiring (parameters) to its predictions. We’ll learn to:
- Ask models questions directly using
.explore()— no more puzzling over what an estimate means - Compute marginal effects — predicted values at specific levels of a predictor
- Make comparisons — pairwise and interaction contrasts
Marginal What?
Consider a typical statistical setting where you ask a research question and selects a model to answer it. Perhaps the model is designed to capture a salient feature of the data-generating process, achieve a satisfactory level of goodness-of-fit, or control for confounders that could introduce bias in a causal estimate. You fits your model using an estimator, and obtain parameter estimates along with standard errors (uncertainty).
Even if a fitted model is relatively simple, the parameter estimates it generates may not map directly onto an estimand that directly informs you research question. The raw parameters may be very difficult to interpret substantively. This can become especially true when we meet Generalized Linear Models soon, e.g. logistic regression for a binary outcome. Typically we get estimates expressed as log-odds-ratios. Despite the simple nature of this model, its parameters are still expressed as complicated functions of probabilities and interactions don’t mean what we think they do!.
If interactions between simple probabilities are challenging to grasp, how can we expect to interpret the estimates produced by fancy regression models with splines, non-linearites, etc? How can we help our audience understand the results produced by mixed-effects Bayesian models, or by complex machine learning algorithms?
This is the main goal of marginal effects analysis. The big idea is that in most cases, analysts should not focus on the raw parameters of their models. Instead we should ask the model what it predicts for different values of our variables based on what its learned from the data. Instead of decoding what \(\hat{\beta}_2 = 0.188\) means, we ask: “What does the model predict when radio = 20 vs radio = 40?”
The Mixing Board
Remember the analogy we learned in class about how the GLM sees variables:
- Continuous: sliders
- Categorical: switches
Let’s make a complex mixer-board-esque regression model with multiple continuous (slider) and categorical (switch) explanatory variables:
A marginal effect is what happens to the output when you move one slider or flip one switch, holding everything else in place.
In bossanova each model has a .explore() method that provides an interface that lets you set the knobs and ask: “What does the model predict at these settings?”
It returns marginal estimated means for categorical variables and marginal estimated slopes for continuous variables — always with proper standard errors and confidence intervals.
The key advantage over the default .params:
marginal effects are invariant to coding schemes. Since we’re always working with model predictions, they don’t depend upon parameterization!
Warm-Up: Credit Dataset
Let’s start with the credit dataset and build up from simple to complex, so the .explore() syntax becomes second nature.
| Variable | Description |
|---|---|
| Balance | Average credit card balance ($) |
| Income | Annual income ($1000s) |
| Ethnicity | African American, Asian, Caucasian |
| Student | Yes / No |
Code
credit = load_dataset("credit")Continuous Predictor: The Slope
Let’s fit a simple model and use .explore() to ask it questions.
Code
m_income = model("Balance ~ Income", credit).fit()
m_income.params
shape: (2, 2)
| term | estimate |
|---|---|
| str | f64 |
| "Intercept" | 246.5 |
| "Income" | 6.048 |
.params tells us the slope is ~6.05 — for each additional $1000 in income, balance increases by about $6. Let’s verify this with .explore()
Code
m_income.explore("Income")
shape: (1, 2)
| Income | estimate |
|---|---|
| str | f64 |
| "marginal_effect" | 6.048 |
That 6.048 is the same as our \(\hat{\beta}_1\) from .params — for a simple model with one continuous predictor, the marginal slope is the parameter estimate. But as we’ll see, this equivalence breaks down with interactions.
The super-power of marginal effects estimation is that it allows you to ask the model all kinds of questions. For example, what if we wanted to zoom into predicted balance at specific income levels?
We can use the @[] syntax to zoom into specific values or levels of any of our predictors:
Code
m_income.explore("Income@[20, 50, 80]")
shape: (3, 2)
| Income | estimate |
|---|---|
| i64 | f64 |
| 20 | 367.5 |
| 50 | 548.9 |
| 80 | 730.4 |
We asked: “What does the model predict at Income = 20, 50, and 80?” The @ syntax pins the predictor at specific values.
For a continuous predictor, .explore() without @ returns the marginal slope — how much the prediction changes per unit increase.
Categorical Predictor: Estimated Marginal Means
Now let’s fit a model with a categorical predictor:
Code
m_eth = model("Balance ~ Ethnicity", credit).fit()
m_eth.params
shape: (3, 2)
| term | estimate |
|---|---|
| str | f64 |
| "Intercept" | 531.0 |
| "Ethnicity[Asian]" | -18.69 |
| "Ethnicity[Caucasian]" | -12.5 |
With treatment coding (the default), .params gives us: an Intercept (the mean for African American, the reference level), and differences from that reference. To get the actual predicted balance for Asian, we’d have to compute \(531.0 + (-18.69) = 512.3\).
Or… we just ask:
Code
m_eth.explore("Ethnicity")
shape: (3, 2)
| Ethnicity | estimate |
|---|---|
| str | f64 |
| "African American" | 531.0 |
| "Asian" | 512.3 |
| "Caucasian" | 518.5 |
These are the marginal estimated means — the model’s predicted Balance for each Ethnicity group. No reference level, no mental arithmetic. Just the predictions.
What about group differences? Treatment coding only gave us differences from the reference. To get all pairwise comparisons:
Code
# We can use .effects just like .params to get a nicer dataframe
m_eth.explore("pairwise(Ethnicity)").effects
shape: (3, 2)
| contrast | estimate |
|---|---|
| str | f64 |
| "EthnicityAsian - EthnicityAfri… | -18.69 |
| "EthnicityCaucasian - Ethnicity… | -12.5 |
| "EthnicityCaucasian - Ethnicity… | 6.184 |
Every pairwise comparison — with SEs, CIs, and p-values — regardless of how the model was coded internally. None of these differences are significant here (all \(p > .05\)), consistent with Ethnicity not being a strong predictor of Balance.
Multiple Predictors & Interactions
Now let’s add an interaction. With Balance ~ Income * Student, the relationship between Income and Balance can differ for students vs non-students:
Code
m_int = model("Balance ~ Income * Student", credit).fit()
m_int.params
shape: (4, 2)
| term | estimate |
|---|---|
| str | f64 |
| "Intercept" | 200.6 |
| "Income" | 6.218 |
| "Student[Yes]" | 476.7 |
| "Income:Student[Yes]" | -1.999 |
4 parameters — and interpreting them requires remembering what “the effect of Income when Student = No” means vs “the additional effect when Student = Yes.” Let’s use the mixing board instead.
EMMs — “What’s the predicted Balance for students vs non-students?”
We can use pairwise() like before or we can use Python style indexing to select levels within a variable to compare:
Code
m_int.explore("Student[Yes - No]")
shape: (1, 2)
| contrast | estimate |
|---|---|
| str | f64 |
| "StudentYes - StudentNo" | 386.3 |
Notice how this differs from the parameter estimate above!
That’s because we didn’t center Income, so the estimate assumes Income = 0 But if we center it, we’ll see that the estimate now matches the marginal effect
Code
model("Balance ~ center(Income) * Student", credit).fit().infer()
shape: (4, 8)
| term | estimate | se | ci_lower | ci_upper | statistic | df | p_value |
|---|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| "Intercept" | 481.8 | 20.64 | 441.2 | 522.4 | 23.34 | 396.0 | 2.2200e-16 |
| "center(Income)" | 6.218 | 0.5921 | 5.054 | 7.382 | 10.5 | 396.0 | 2.2200e-16 |
| "Student[Yes]" | 386.3 | 65.36 | 257.8 | 514.8 | 5.91 | 396.0 | 7.3780e-9 |
| "center(Income):Student[Yes]" | -1.999 | 1.731 | -5.403 | 1.404 | -1.155 | 396.0 | 0.2489 |
Students carry about $386 more in balance — averaged across all Income levels.
Another super-powerof marginal effects is helping you unpack complex interactions.
For example, we can see that the interaction between Income (centered) and Student is -1.99 above. But what does this mean and in which direction?
If you remember the default treatment/dummy-coding, then you’ll remember that “No” is the default reference group for Student. So this is the Income slope for Students - the Income slope for Non-Students.
But if we used a different encoding scheme this would change:
Code
# Sum/deviation coding
model("Balance ~ center(Income) * sum(Student)", credit).fit().infer()
shape: (4, 8)
| term | estimate | se | ci_lower | ci_upper | statistic | df | p_value |
|---|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| "Intercept" | 674.9 | 32.68 | 610.7 | 739.2 | 20.65 | 396.0 | 2.2200e-16 |
| "center(Income)" | 5.219 | 0.8656 | 3.517 | 6.92 | 6.029 | 396.0 | 3.7890e-9 |
| "Student[No]" | -193.1 | 32.68 | -257.4 | -128.9 | -5.91 | 396.0 | 7.3780e-9 |
| "center(Income):Student[No]" | 0.9996 | 0.8656 | -0.7022 | 2.701 | 1.155 | 396.0 | 0.2489 |
Notice that the estimate for the interaction change but the inference test result did not (same t-statistic and p-value).
You can imagine with more complex models with additional interactions (e.g. 3-way, 4-way) things will get really hard to interpret really fast!
But we can simplify things by relying on marginal effects to do the right thing. .explore() supports a formula syntax similar to the formula syntax you use when creating models so you can slice-and-dice your data to ask the model any number of questions:
Slopes by group — “How does the Income-Balance relationship differ?”
Code
m_int.explore("Income ~ Student")
shape: (2, 3)
| Student | Income | estimate |
|---|---|---|
| str | str | f64 |
| "No" | "marginal_effect" | 6.218 |
| "Yes" | "marginal_effect" | 4.219 |
The Income slope is ~$6.2/thousand for non-students vs ~$4.2/thousand for students. The interaction means income has a weaker effect on balance for students. This is why the original interaction estimate with the default dummy coding was -1.99! It reflects the difference in these two numbers!
You can combine this with the @[] syntax to focus on specific values.
The ~ separates the focal variable (what you want EMMs/slopes for) from the conditioning variables (what you want to hold fixed):
Code
# Compare students and non-students at specific income levels
m_int.explore("Student[Yes - No] ~ Income@[20, 50, 80]")
shape: (3, 3)
| Income | contrast | estimate |
|---|---|---|
| f64 | str | f64 |
| 20.0 | "StudentYes - StudentNo" | 436.7 |
| 50.0 | "StudentYes - StudentNo" | 376.7 |
| 80.0 | "StudentYes - StudentNo" | 316.7 |
Penguins - A Complex Model
The credit examples had at most one categorical and one continuous predictor. Now let’s build a model with multiple categorical and continuous predictors interacting, where .explore() really becomes essential.
We’ll model penguin body mass as a function of sex, species, and flipper length — with a full three-way interaction:
\[ \text{body\_mass\_g} \sim \text{sex} \times \text{species} \times \text{flipper\_length\_mm} \]
This model allows the relationship between flipper length and body mass to differ across every combination of sex and species.
Data Loading (& a quick lesson on missing values)
Before modeling, let’s check our data quality. A quick habit to build: always inspect your data for missing or unexpected values before fitting a model.
Code
penguins_raw = load_dataset("penguins")
# Check for nulls
penguins_raw.null_count()
shape: (1, 8)
| species | island | bill_length_mm | bill_depth_mm | flipper_length_mm | body_mass_g | sex | year |
|---|---|---|---|---|---|---|---|
| u32 | u32 | u32 | u32 | u32 | u32 | u32 | u32 |
| 0 | 0 | 2 | 2 | 2 | 2 | 0 | 0 |
Code
# Check unique values in the sex column
penguins_raw.select("sex").unique()
shape: (3, 1)
| sex |
|---|
| str |
| "female" |
| "NA" |
| "male" |
Code
# Check unique values in the sex column
penguins_raw.select("species").unique()
shape: (3, 1)
| species |
|---|
| str |
| "Gentoo" |
| "Adelie" |
| "Chinstrap" |
Notice something sneaky: sex has three unique values — "female", "male", and "NA". That "NA" is a string, not a true missing value! Polars’s .null_count() shows 0 nulls for sex because the string "NA" isn’t null — it’s just a text label.
This is a common gotcha in real datasets. If we tried to fit a model without catching this, sex would become a three-level factor, creating nonsense parameters.
The fix: explicitly filter out the string "NA", then drop any remaining null rows from other columns:
Code
penguins = penguins_raw.filter(col("sex") != "NA").drop_nulls()
penguins.shape(333, 8)333 clean rows with 2 sex levels. Now we can fit our model.
Fitting the 3-Way Interaction
Below we’ve setup the full model for you (with default coding scheme) and outputted the parameter estimates.
12 parameters. Good luck interpreting sex[male]:species[Gentoo]:flipper_length_mm = -23.27 from the parameter table alone.
Let’s use .explore() to try unpack the estimates. We’ll start by using .jointtest() and then use .explore() to zoom in.
Code
m_penguin = model("body_mass_g ~ sex * species * flipper_length_mm",penguins).fit()
m_penguin.params
shape: (12, 2)
| term | estimate |
|---|---|
| str | f64 |
| "Intercept" | 991.5 |
| "sex[male]" | -592.6 |
| "species[Chinstrap]" | 233.4 |
| "species[Gentoo]" | -3804.0 |
| "sex[male]:species[Chinstrap]" | -4744.0 |
| … | … |
| "sex[male]:flipper_length_mm" | 6.283 |
| "species[Chinstrap]:flipper_len… | -0.6514 |
| "species[Gentoo]:flipper_length… | 22.57 |
| "sex[male]:species[Chinstrap]:f… | 21.98 |
| "sex[male]:species[Gentoo]:flip… | -23.27 |
Code
# ANOVA
m_penguin.jointtest()
shape: (7, 5)
| term | df1 | df2 | f_ratio | p_value |
|---|---|---|---|---|
| str | i64 | f64 | f64 | f64 |
| "sex" | 1 | 321.0 | 0.200843 | 0.654344 |
| "species" | 2 | 321.0 | 1.325464 | 0.26713 |
| "flipper_length_mm" | 1 | 321.0 | 55.402216 | 9.0955e-13 |
| "sex:species" | 2 | 321.0 | 3.912632 | 0.020948 |
| "sex:flipper_length_mm" | 1 | 321.0 | 0.905627 | 0.341993 |
| "species:flipper_length_mm" | 2 | 321.0 | 1.726036 | 0.179636 |
| "sex:species:flipper_length_mm" | 2 | 321.0 | 3.679455 | 0.026307 |
Let’s try to understand this 3-way interaction (last-row above)
Estimated Marginal Means
“What’s the model’s predicted body mass for each species?”
Code
m_penguin.explore("species")
shape: (3, 2)
| species | estimate |
|---|---|
| str | f64 |
| "Adelie" | 3885.0 |
| "Chinstrap" | 3823.0 |
| "Gentoo" | 4703.0 |
Gentoo penguins are the heaviest (~4703g), with Adelie (~3885g) and Chinstrap (~3823g) close together. These EMMs are averaged across both sexes and all flipper lengths.
“What about by sex?”
Code
m_penguin.explore("pairwise(sex)")
shape: (1, 2)
| contrast | estimate |
|---|---|
| str | f64 |
| "sexmale - sexfemale" | 665.2 |
Males are about 666g heavier than females on average.
To get the full picture, we can cross species with sex:
Code
m_penguin.explore("species ~ sex[male - female]")
shape: (3, 3)
| species | rhs_contrast | estimate |
|---|---|---|
| str | str | f64 |
| "Adelie" | "sexmale - sexfemale" | 670.0 |
| "Chinstrap" | "sexmale - sexfemale" | 343.4 |
| "Gentoo" | "sexmale - sexfemale" | 843.4 |
Now we’re starting to see the interaction: the sex difference is much larger for Gentoo (~844g) than for Chinstrap (~343g) or Adelie (~670g).
Marginal Slopes
“How does body mass change with flipper length?”
Code
m_penguin.explore("flipper_length_mm")
shape: (1, 2)
| flipper_length_mm | estimate |
|---|---|
| str | f64 |
| "marginal_effect" | 22.89 |
Overall, each additional mm of flipper length is associated with ~22.9g more body mass. But this is the average slope — our 3-way interaction means it varies by group:
Code
m_penguin.explore("flipper_length_mm ~ species")
shape: (3, 3)
| species | flipper_length_mm | estimate |
|---|---|---|
| str | str | f64 |
| "Adelie" | "marginal_effect" | 15.8 |
| "Chinstrap" | "marginal_effect" | 26.14 |
| "Gentoo" | "marginal_effect" | 26.73 |
The flipper-mass slope is weakest for Adelie (~15.8 g/mm) and strongest for Gentoo (~26.7 g/mm) and Chinstrap (~26.1 g/mm). We can also break it out by sex:
Code
m_penguin.explore("flipper_length_mm ~ sex")
shape: (2, 3)
| sex | flipper_length_mm | estimate |
|---|---|---|
| str | str | f64 |
| "female" | "marginal_effect" | 19.96 |
| "male" | "marginal_effect" | 25.82 |
Code
m_penguin.explore("flipper_length_mm ~ species + sex")
shape: (6, 4)
| species | sex | flipper_length_mm | estimate |
|---|---|---|---|
| str | str | str | f64 |
| "Adelie" | "female" | "marginal_effect" | 12.66 |
| "Adelie" | "male" | "marginal_effect" | 18.94 |
| "Chinstrap" | "female" | "marginal_effect" | 12.01 |
| "Chinstrap" | "male" | "marginal_effect" | 40.27 |
| "Gentoo" | "female" | "marginal_effect" | 35.23 |
| "Gentoo" | "male" | "marginal_effect" | 18.24 |
Interaction Contrasts: “Difference of Differences”
The direction reversal raises a natural question: “Is the Chinstrap-Adelie difference between males and females?”
That’s an interaction contrast — a difference of differences. We use brackets on both sides:
Code
m_penguin.explore("species[Chinstrap - Adelie] ~ sex[male - female]").infer()
shape: (1, 9)
| contrast | rhs_contrast | estimate | se | ci_lower | ci_upper | statistic | df | p_value |
|---|---|---|---|---|---|---|---|---|
| str | str | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| "speciesChinstrap - speciesAdel… | "sexmale - sexfemale" | -326.6 | 146.7 | -615.2 | -37.9 | -2.226 | 321.0 | 0.02673 |
The Chinstrap-Adelie gap is ~327g more negative for males than females (\(p = .0267\)). This is the sex × species interaction, focused on exactly the comparison we care about.
Slope Contrasts: “Does the relationship differ?”
We can also compare how the flipper-mass slope differs between groups:
Code
# Does the flipper slope differ between Gentoo and Adelie?
m_penguin.explore("flipper_length_mm ~ species[Gentoo - Adelie]").infer()
m_penguin.effects
shape: (1, 9)
| flipper_length_mm | rhs_contrast | estimate | se | ci_lower | ci_upper | statistic | df | p_value |
|---|---|---|---|---|---|---|---|---|
| str | str | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| "marginal_effect" | "speciesGentoo - speciesAdelie" | 10.93 | 7.026 | -2.891 | 24.76 | 1.556 | 321.0 | 0.1207 |
shape: (1, 9)
| flipper_length_mm | rhs_contrast | estimate | se | ci_lower | ci_upper | statistic | df | p_value |
|---|---|---|---|---|---|---|---|---|
| str | str | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| "marginal_effect" | "speciesGentoo - speciesAdelie" | 10.93 | 7.026 | -2.891 | 24.76 | 1.556 | 321.0 | 0.1207 |
The flipper-mass slope is ~10.9 g/mm steeper for Gentoo than Adelie — each additional mm of flipper length “buys” more body mass in Gentoo penguins.
Summary
We now have three complementary tools for interrogating a fitted model:
| Question | Tool | What it tells you |
|---|---|---|
| “Is this predictor worth including?” | .jointtest() or compare() |
Yes/no — invariant to coding |
| “What are the model’s internal weights?” | .params |
Coefficients — depend on coding |
| *“What does the model predict for ___?“* | .explore() |
Predictions — invariant to coding |
| “How uncertain am I??” | .infer() |
Inference - after .fit() or .explore() |
Challenge: Exploring the Credit Dataset
Now it’s your turn. The credit dataset has many predictors — use the mixing board to investigate:
- Fit a model:
Balance ~ Income * Student + Ethnicity * Gender - Use
.jointtest()— which terms matter? - Compute EMMs for
StudentandEthnicity— what are the predicted Balances? - Compare Income slopes for students vs non-students — does Income affect Balance differently?
- Use
pairwise()or bracket contrasts to compare Ethnicity groups - Condition on specific Income values with
@— how do predictions change? - Write up anything you find as if you were writing a Results section
Code
# Your code hereCode
# Your code hereCode
# Your code hereCode
# Your code hereCode
# Your code hereYour summary here…