Audience: graduate students and researchers in ecology who fit GL(M)Ms and want to interpret interaction terms correctly.
Format: short lecture blocks, and practice exercises.
Learning objectives
By the end of this workshop you should be able to:
Explain what an interaction term means for a categorical × numerical model and for a numerical × numerical model
Distinguish conditional effects from marginal effects, and choose the right one for a research question
Use marginaleffects (predictions(), slopes(), avg_predictions(), avg_slopes(), hypotheses()) to extract and test the quantities you actually care about
Learning objectives (continued)
Visualize interaction effects (slope-by-moderator plots, five-number-summary predictions, response surfaces) and know their trade-offs
Recognize that a non-linear link function (logit, log) creates an implicit interaction even when the model formula has none, and convert a raw GLM coefficient to a probability-scale or count-scale effect
Interpret a hurdle/zero-inflated model’s two components separately and combined, and recognize that two-way interactions are symmetric
Actual learning objective
Model summary outputs rarely answer a research question and definitely do not communicate the results in a comprehensible way.
Setup
library(tidyverse) # Data manipulation and plottinglibrary(glmmTMB) # Generalized linear (mixed) modelslibrary(marginaleffects) # Compute marginal/conditional effectslibrary(patchwork) # Combine multiple plotslibrary(GGally) # Pairwise plotslibrary(distributional) # Simulate mixture/gamma/normal distributions (Part 3)library(ggdist) # stat_slab() for predicted-distribution plots (Part 3)library(tinytable) # tt(): simple, nice-looking results tableslibrary(performance) # r2(), r2_tjur(): model R-squared measures# Set ggplot2 themeblack_theme <-theme(axis.line =element_line(linewidth =2, lineend ="round", color ="white"),panel.grid =element_blank(),panel.background =element_rect(fill ="black", color =NA),axis.ticks =element_blank(),axis.text =element_text(size =21, face ="bold", color ="white"),axis.title =element_text(size =21, face ="bold", color ="white"),plot.title =element_text(size =25, face ="bold", color ="white"),plot.background =element_rect(fill ="black", color =NA),legend.background =element_rect(fill ="black", color =NA),legend.text =element_text(size =21, color ="white"),legend.title =element_text(size =21, face ="bold", color ="white"),legend.position ="bottom")theme_set(black_theme)# Set global default for geom_point shapeupdate_geom_defaults("point", list(shape =1, size =3, alpha =0.8))# silence marginaleffects warnings about glmmTMBoptions(marginaleffects_safe =FALSE)set.seed(333)
Coding style & functions used throughout
A few conventions and workhorse functions you’ll see on repeat:
tibble() + rnorm()/rbinom()/rpois()/runif(): how we simulate data with a known true effect, so we can check whether our model recovers it
glmmTMB(y ~ x1 * x2, data = ..., family = ...): our model-fitting workhorse; * expands to x1 + x2 + x1:x2
marginaleffects verbs:
predictions() / avg_predictions() — model-implied values of \(Y\)
slopes() / avg_slopes() — model-implied effects (derivatives) of a predictor on \(Y\)
datagrid() — build a covariate grid to evaluate predictions/slopes at (e.g., fivenum, specific values, mean)
by = — keep results split by a grouping variable instead of averaging over it
hypotheses() — test linear combinations of estimates (e.g., “is slope A different from slope B?”)
set.seed(): called once per script so simulated data (and thus every number you see) is reproducible
Part 1 — Categorical × Numerical Interactions
What is an interaction?
An interaction means the effect of one predictor depends on the value of another.
No interaction: snow depth affects weight the same way for every sheep
Interaction: the effect of snow depth on weight differs by sex
For a categorical moderator (e.g., sex) and a continuous predictor (e.g., snow depth): the interaction means the slope differs between categories.
Our focus today isn’t the algebra — it’s using marginaleffects to find, test, and plot “it depends.”
Conditional vs. marginal effects
Conditional effects: the effect of a predictor at specific values of other predictors (e.g., the slope of snow depth for males)
Marginal effects: the effect averaged across the distribution of other predictors (e.g., the average slope of snow depth across the whole population)
Both are legitimate; they answer different questions. Interactions are exactly where the distinction starts to matter.
Case study: Bighorn sheep
Ram Mountain, Alberta. Question: how does winter snow depth relate to spring mass in Bighorn sheep, and does it differ by sex? Males and females segregate in winter — females use lower, less snowy habitat.
Built-in features: different intercepts, different slopes, and different snow-depth distributions by sex.
Simpson’s Paradox
Pooling males and females together and fitting a single trend line can hide — or even reverse — the true relationship.
Simpson’s Paradox
Code
ggplot(bighorn_population, aes(x = snow_depth, y = weight)) +geom_point(aes(color = sex)) +geom_smooth(method ="lm", se =FALSE, color ="red", linewidth =2)
Ignoring sex hides the true within-sex trends.
Pooled vs. within-sex
Pooled: snow depth looks positively associated with weight — even though both sexes lose weight with more snow.
Why? Sex is a confounder
Within each sex: more snow → lower weight (true negative relationship)
Between sexes: males see more snow and are heavier (confounded)
Pooled data compares heavy, high-snow males to light, low-snow females → spurious positive trend
Within-sex trends
Code
ggplot(bighorn_population) +geom_point(aes(x = snow_depth, y = weight, color = sex)) +geom_smooth(aes(x = snow_depth, y = weight), method ="lm", se =FALSE, color ="red", linewidth =2) +geom_smooth(aes(x = snow_depth, y = weight, color = sex), method ="lm", se =FALSE)
Within-sex trends are negative, as expected.
Exercise 1
Using bighorn_population:
Fit weight ~ sex * snow_depth with glmmTMB(). Call it mod_mass_true.
From the summary, identify which coefficient represents the male slope of snow depth. (Hint: it’s not printed directly — you have to add two rows together.)
In one sentence, explain why a naive weight ~ snow_depth model (ignoring sex) would mislead a manager trying to understand the effect of a harsh winter. Can you think of an example of such situation in your own field or research?
(~5 minutes — discuss with a neighbor)
Sampling from the population
In practice we can’t measure everyone. Two sampling scenarios:
Using the model that was fit to the entire population and the one with a biased sample, use the avg_predictions() function along with a balanced and unbalanced prediction grid to create a plot comparing the marginal predictions between the true population, biased population predictions and corrected population predictions. There should be three trend lines.
Interactions → conditional effects that depend on other variables
Simpson’s Paradox is a real risk when pooling across confounded groups
Marginal ≠ conditional; both are valid but answer different questions
Raw model coefficients rarely equal your research quantity of interest — use marginaleffects
Marginal effects inherit your sample’s covariate distribution unless you override it with newdata
Part 2 — Numerical × Numerical Interactions
Two numerical predictors
With two continuous predictors, an interaction means: the slope of one predictor changes as the other predictor changes — not a separate slope per group, but a sliding one.
A raw coefficient is only “the slope when the other predictor equals zero” — meaningful only if zero is a meaningful value. This motivates centering predictors on their mean, so the intercept and main-effect coefficients describe average conditions.
Case study: petal length
Simulated flowers; petal length as a function of temperature and rain, with an interaction.
Predictors are mean-centered: the intercept becomes the predicted response at average conditions rather than at (biologically meaningless) zero temperature/rain.
Every term has a small p-value — but what do these numbers mean for our research question? Time to unpack them.
Slope of one predictor at a single value of the other
Code
# Effect of temperature when rain_centered = 0 (i.e., average rain)slopes(mod_flower_1, newdata =datagrid(rain_centered =0),variables ="temp_centered")
predictions(mod_flower_1,newdata =datagrid(temp_centered =seq(min(flower_measurements_1$temp_centered),max(flower_measurements_1$temp_centered), length.out =30),rain_centered = fivenum )) |>as_tibble() |>mutate(rain =round(rain_centered +mean(rain), 2)) |>ggplot() +geom_line(aes(y = estimate, x = temp_centered +mean(temperature), color =as.factor(rain)),lineend ="round", linewidth =1.5) +geom_ribbon(aes(ymin = conf.low, ymax = conf.high,x = temp_centered +mean(temperature), fill =as.factor(rain)), alpha =0.2) +labs(y ="Predicted petal length", x ="Temperature", color ="Rain", fill ="Rain") +coord_cartesian(ylim =c(0, 50))
Response surfaces — use with caution
Code
predictions(mod_flower_1,newdata =datagrid(temp_centered =seq(min(flower_measurements_1$temp_centered),max(flower_measurements_1$temp_centered), length.out =30),rain_centered =seq(min(flower_measurements_1$rain_centered),max(flower_measurements_1$rain_centered), length.out =30) )) |>as_tibble() |>mutate(rain = rain_centered +mean(rain), temp = temp_centered +mean(temperature)) |>ggplot(aes(x = rain, y = temp, z = estimate)) +stat_contour_filled() +labs(x ="Rain", y ="Temperature", fill ="Predicted petal length")
Response surfaces: a caution
Contour bins are arbitrary and hide uncertainty — good for exploration, risky for publication-quality inference. Prefer line-and-ribbon plots when you can.
Exercise 3
Using mod_flower_1:
Compute the slope of temperature across the observed range of rain (mirror the rain-across-temperature plot above).
At what value of rain does the temperature effect become non-significant (CI crosses zero), if at all?
Rewrite the research question this model could answer as a single sentence, identifying which variable is the exposure and which is the moderator.
(~10 minutes)
Interactions beyond the mean: distributional regression
Interactions don’t have to act only on the average response. glmmTMB’s dispformula lets predictors affect the standard deviation too — useful when variance itself is ecologically meaningful (e.g., reduced trait variance under stress).
We can interpret the interactive effects of rain and temperature on the standard deviation of petal length in the same way that we would for the average. Examples are shown in the Numerical x Numerical Interactions tutorial online. I want to display another way to communicate the results. Knowing that the two parameters of the Normal distribution act together to generate the observations, we can predict new distributions of observations for different predictor values.
Predicting new observations
Code
predictions(mod_flower_measurements_3,newdata =datagrid(temp_centered = fivenum,rain_centered = fivenum ),type ="conditional") |>as_tibble() |>select(temp_centered, rain_centered, estimate) |>rename(mu = estimate) |>left_join(predictions(mod_flower_measurements_3,newdata =datagrid(temp_centered = fivenum,rain_centered = fivenum ),type ="disp") |>as_tibble() |>select(temp_centered, rain_centered, estimate) |>rename(sd = estimate)) |>rowwise() |>mutate(temp_centered =round(temp_centered),rain_centered =round(rain_centered),pred_obs =dist_normal(mu = mu, sd = sd)) -> predicted_observations# Build facet labels dynamically from the actual rounded temp_centered values# (hardcoding them would silently show NA facet labels if the simulated data changes)temp_labels <-setNames(paste0("temp = ", sort(unique(predicted_observations$temp_centered))),sort(unique(predicted_observations$temp_centered)))ggplot(predicted_observations, aes(xdist = pred_obs, color =as.factor(rain_centered))) +stat_slab(fill =NA, normalize ="xy") +scale_color_brewer(palette ="Accent") +facet_wrap(~as.factor(temp_centered), scales ="free", labeller =as_labeller(temp_labels)) +labs(x ="Predicted petal length", color ="Rain (centered)") +theme(axis.line.y =element_blank(),axis.text.y =element_blank(),axis.title.y =element_blank(),strip.text =element_text(color ="black", face ="bold", size =12))
Part 2 takeaways
With two numerical predictors, the interaction coefficient is the rate of change of one slope per unit of the other — not a separate slope for each group
Centering predictors makes main-effect coefficients interpretable as “effect at average conditions”
Plot slopes across the moderator’s range, not just the coefficient table
Five-number-summary prediction plots are usually more honest than response surfaces (they retain uncertainty)
Interactions can act on variance, not just the mean, via distributional regression (dispformula)
Part 3 — Link Functions Create Implicit Interactions
A different kind of interaction
So far, every interaction we’ve seen was written explicitly into the model formula (X1 * X2).
Non-linear link functions (logit, log, …) introduce a second, implicit source of interaction: even a model with no interaction term at all can show predictors whose response-scale effect depends on the values of other predictors.
This happens because the link function bends a straight line (on the link scale) into a curve (on the response scale) — and how much a fixed push on the link scale moves the response depends on where you are on that curve.
Case study: species occurrence
Simulated presence/absence of a species across 200 plots, as a function of tree cover, rain, and temperature — no interaction terms, fit with a binomial GLM (logit link):
We simulated the effect of tree cover to be 0.45 on the probability scale, but the model reports ~3.76. Did the model get it wrong?
No — it’s reporting a different scale
The logit-link GLM draws a straight line through the log-odds of occurrence, then bends that line into an S-shape to get back to a probability.
The logit link, visualized
Code
tibble(log_odds =seq(-6, 6, length.out =400)) |>mutate(probability =plogis(log_odds)) -> logistic_curvetibble(log_odds =c(0, -4.2, 4.2),probability =plogis(log_odds),label =c("steep middle:\na small push here\nmoves probability a lot","flat end:\na push here barely\nmoves probability","flat end:\na push here barely\nmoves probability")) -> annotggplot(logistic_curve, aes(x = log_odds, y = probability)) +geom_line(color ="white", linewidth =1.5, lineend ="round") +geom_point(data = annot, color ="white", size =3) +geom_text(data = annot,aes(label = label),color ="white",size =4.5,fontface ="bold",lineheight =0.9,vjust =c(-0.4, -0.6, 1.6),hjust =c(0.5, 0.5, 0.5) ) +scale_y_continuous(limits =c(-0.1, 1.05)) +labs(x ="log-odds scale (what the model's coefficient describes)",y ="probability scale",title ="The logit link bends a straight line into an S-shape" )
The logit link bends a straight line into an S-shape.
Reading the S-curve
Near the middle (p ≈ 0.5) the curve is steep: a small log-odds push moves probability a lot. Near the ends (p ≈ 0 or 1) the curve flattens: the same push barely moves probability.
No manual derivative needed — marginaleffects converts to the scale we want.
slopes() in this case asks: “at these covariate values, how much does predicted probability change for a small increase in tree cover?” — a probability-scale answer, directly comparable to how we simulated the data.
The tree-cover slope on the probability scale changes across rain × temp, though no interaction term exists in the model.
Why the slope moves around
Because the amount of “bending” depends on where you are on the S-curve, tree cover’s probability-scale effect depends on rain and temp too — even though we never told the model to include that interaction.
A note on measurement scale (Spake et al. 2023)
Whether an interaction “exists” can depend on which scale you look at: additive (probability, or log-odds as a straight difference) vs. multiplicative (ratios, or the log scale itself).
mod_1 has no interaction term, yet:
on the log-odds scale: no interaction, by construction
on the probability scale: an interaction appears, purely from the S-curve’s shape
Getting the scale question wrong is one way real interactions get mis-detected or mis-sized (a “Type-D” error in Spake et al.’s terminology).
Response vs. link scale, across each moderator
Code
slopes_by_temp_response <-slopes(mod_1,newdata =datagrid(temp =seq(min(occ_data$temp), max(occ_data$temp), length.out =30),rain =mean(occ_data$rain)),variables ="tree_cov", by ="temp", type ="response") |>as_tibble()slopes_by_temp_link <-slopes(mod_1,newdata =datagrid(temp =seq(min(occ_data$temp), max(occ_data$temp), length.out =30),rain =mean(occ_data$rain)),variables ="tree_cov", by ="temp", type ="link") |>as_tibble()slopes_by_rain_response <-slopes(mod_1,newdata =datagrid(rain =seq(min(occ_data$rain), max(occ_data$rain), length.out =30),temp =mean(occ_data$temp)),variables ="tree_cov", by ="rain", type ="response") |>as_tibble()slopes_by_rain_link <-slopes(mod_1,newdata =datagrid(rain =seq(min(occ_data$rain), max(occ_data$rain), length.out =30),temp =mean(occ_data$temp)),variables ="tree_cov", by ="rain", type ="link") |>as_tibble()plot_slope <-function(data, x, y_lab, scale_title) {ggplot(data) +geom_hline(yintercept =0, color ="blue", linewidth =1, linetype =2) +geom_ribbon(aes(x = .data[[x]], ymin = conf.low, ymax = conf.high), alpha =0.4, fill ="lightgrey") +geom_line(aes(x = .data[[x]], y = estimate), linewidth =1.5, lineend ="round", color ="white") +labs(y = y_lab, x = x, title = scale_title)}(plot_slope(slopes_by_temp_response, "temp", "Slope of tree_cov", "Response scale") +plot_slope(slopes_by_temp_link, "temp", NULL, "Link scale")) /(plot_slope(slopes_by_rain_response, "rain", "Slope of tree_cov", "Response scale") +plot_slope(slopes_by_rain_link, "rain", NULL, "Link scale"))
Slope of tree_cov across temp (top) and rain (bottom), on the probability (left) vs. log-odds (right) scale. On the link scale the slope is flat — no explicit interaction was fit.
Reading these panels
Link (log-odds) scale: the slope of tree_cov is flat (~1.98) across both moderators — exactly as it should be, since mod_1’s formula has no interaction term at all.
Response (probability) scale: the slope visibly curves with temp (and more weakly with rain) — the implicit interaction made concrete, purely from the logit link’s curvature.
With three predictors, a focal variable’s slope can no longer be shown across “the” moderator in one 2D plot — each panel above holds the other remaining predictor fixed at its mean. Alternatives: facet by the third variable, use small multiples at representative (e.g. five-number-summary) values, or report a single avg_slopes() that integrates over the observed distribution instead of fixing it.
Predicted probability of occurrence across tree_cov, averaged over the observed rain and temp values (counterfactual averaging).
avg_predictions()/avg_slopes() give one number (or curve) per focal value, integrating over rain and temp as observed in the sample — a marginal, “average effect” alternative to picking a single moderator value to condition on.
Exercise 4
Using mod_1:
Plot the predicted trend of occurence as a function of tree cover for the unique combinations of the five number summary of rain and temperature. Hint: use the predictions() function.
Plot the prediction for the average effect of temperature.
Average predicted probability of occurrence across temp, integrating over the observed distribution of rain and tree_cov.
Reviewer 2 wants a table
Code
# Average slopes on the probability scale for all three predictorsslopes_tbl_mod1 <-avg_slopes(mod_1, variables =c("tree_cov", "rain", "temp"), type ="response") |>as_tibble() |>select(term, estimate, conf.low, conf.high, p.value) |>mutate(across(c(estimate, conf.low, conf.high), ~round(.x, 3)),p.value =if_else(p.value <0.001, "<0.001", as.character(round(p.value, 3))) ) |>mutate(CI_95 = glue::glue("[{conf.low}, {conf.high}]")) |>select(Predictor = term, Slope = estimate, `95% CI`= CI_95, `p-value`= p.value)# Tjur's R2 as an extra row in the same 4-column layoutr2_row <-tibble(Predictor ="Tjur's R\u00b2",Slope =as.character(round(r2_tjur(mod_1), 3)),`95% CI`="",`p-value`="")mod1_tbl <-bind_rows( slopes_tbl_mod1 |>mutate(across(everything(), as.character)), r2_row)mod1_tbl |>tt(caption ="Average slopes of occurrence on the probability scale, and Tjur's R2 for mod_1") |>style_tt(i =nrow(slopes_tbl_mod1) +1, bold =TRUE, line ="t", line_color ="black", line_width =0.1)
Average slopes of occurrence on the probability scale, and Tjur's R2 for mod_1
Predictor
Slope
95% CI
p-value
rain
-0.008
[-0.04, 0.025]
0.639
temp
0.023
[0.015, 0.031]
<0.001
tree_cov
0.583
[0.371, 0.795]
<0.001
Tjur's R²
0.272
Case study: species abundance
The same phenomenon shows up with count data. Abundance (when present) as a function of the same three predictors, fit with a Poisson GLM (log link), again no interaction terms:
tibble(log_count =seq(-2, 6, length.out =400)) |>mutate(count =exp(log_count)) -> exp_curvetibble(log_count =c(-1.5, 2, 4.5),count =exp(log_count),label =c("flat end:\na push here barely\nmoves the expected count","moderate slope:\na push here moves\nthe expected count some","steep end:\na push here moves\nthe expected count a lot")) -> annot_expggplot(exp_curve, aes(x = log_count, y = count)) +geom_line(color ="white", linewidth =1.5, lineend ="round") +geom_point(data = annot_exp, color ="white", size =3) +geom_text(data = annot_exp,aes(label = label),color ="white",size =4.5,fontface ="bold",lineheight =0.9,vjust =c(-0.4, -0.6, 1.6),hjust =c(0.5, 0.5, 0.5) ) +labs(x ="log scale (what the model's coefficient describes)",y ="count scale",title ="The log link bends a straight line into an exponential curve" )
From log scale to counts
The log-scale tree_cov coefficient is ~0.73, not 4.5 — the same phenomenon as the logit link, just a different curve.
Two ways to report a log-link coefficient
Rate ratio: a constant multiplicative effect, valid everywhere — “each unit increase in tree cover multiplies expected abundance by this factor, regardless of rain or temp”
Local slope: an additive, count-scale effect (“this many more individuals”), but only valid near the specific covariate values used
Computing the rate ratio
Code
exp(fixef(mod_2)$cond["tree_cov"]) # rate ratio: ~2.08, i.e. a ~108% increase per unit tree cover
tree_cov
2.010566
Both describe the same underlying model — pick the one that matches your audience (percent change vs. a concrete headcount).
Exercise 5
Using mod_2:
Compute the rate ratio and the local count-scale slope for temp, holding rain at its minimum (a “drought” scenario) and tree_cov at its mean.
Repeat holding tree_cov at a low vs. a high value (e.g., the 10th and 90th percentiles) — does the rate ratio change? Does the local slope?
Explain, in your own words, why one of these two summaries stays constant across tree_cov and the other doesn’t.
(~7 minutes)
Exercise 5: answer (task 1)
Code
# Count-scale (additive) slope of temp under drought (rain = min), tree_cov at its meantemp_count_drought <-avg_slopes( mod_2,variables ="temp",newdata =datagrid(rain =min(abund_data$rain), tree_cov =mean(abund_data$tree_cov)),type ="response") |>as_tibble()# Rate ratio (exponentiated log-scale) slope of temp at the same pointtemp_link_drought <-avg_slopes( mod_2,variables ="temp",newdata =datagrid(rain =min(abund_data$rain), tree_cov =mean(abund_data$tree_cov)),type ="link") |>as_tibble()temp_drought_tbl_ex6 <- temp_count_drought |>select(count_est = estimate, count_low = conf.low, count_high = conf.high) |>bind_cols( temp_link_drought |>mutate(rr =exp(estimate), rr_low =exp(conf.low), rr_high =exp(conf.high)) |>select(rr, rr_low, rr_high, p.value) ) |>mutate(across(c(count_est, count_low, count_high, rr, rr_low, rr_high), ~round(.x, 3)),p.value =if_else(p.value <0.001, "<0.001", as.character(round(p.value, 3))),count_CI = glue::glue("[{count_low}, {count_high}]"),rr_CI = glue::glue("[{rr_low}, {rr_high}]") ) |>select(`Count-scale slope`= count_est, `(95% CI)`= count_CI,`Rate ratio`= rr, `(95% CI) `= rr_CI, `p-value`= p.value)temp_drought_tbl_ex6 |>tt(caption ="Effect of temperature on abundance under drought (rain at its minimum), at mean tree cover")
Effect of temperature on abundance under drought (rain at its minimum), at mean tree cover
Count-scale slope
(95% CI)
Rate ratio
(95% CI)
p-value
0.121
[0.082, 0.159]
1.018
[1.012, 1.024]
<0.001
Exercise 5: answer (task 2)
Code
# Low vs. high tree cover: 10th and 90th percentileslow_tree_q <-quantile(abund_data$tree_cov, 0.1)high_tree_q <-quantile(abund_data$tree_cov, 0.9)# Count-scale slope of temp at each tree_cov percentile, under drought (rain = min)temp_count_scale_ex6 <-avg_slopes( mod_2,variables ="temp",by ="tree_cov",newdata =datagrid(tree_cov =c(low_tree_q, high_tree_q), rain =min(abund_data$rain)),type ="response") |>as_tibble()# Rate ratio (exponentiated log-scale) slope of temp at the same two pointstemp_link_scale_ex6 <-avg_slopes( mod_2,variables ="temp",by ="tree_cov",newdata =datagrid(tree_cov =c(low_tree_q, high_tree_q), rain =min(abund_data$rain)),type ="link") |>as_tibble()temp_tree_tbl_ex6 <- temp_count_scale_ex6 |>select(tree_cov, count_est = estimate, count_low = conf.low, count_high = conf.high) |>left_join( temp_link_scale_ex6 |>mutate(rr =exp(estimate), rr_low =exp(conf.low), rr_high =exp(conf.high)) |>select(tree_cov, rr, rr_low, rr_high, p.value),by ="tree_cov" ) |>mutate(tree_cov_label =if_else(tree_cov == low_tree_q, "Low tree cover (p10)", "High tree cover (p90)"),across(c(count_est, count_low, count_high, rr, rr_low, rr_high), ~round(.x, 3)),p.value =if_else(p.value <0.001, "<0.001", as.character(round(p.value, 3))),count_CI = glue::glue("[{count_low}, {count_high}]"),rr_CI = glue::glue("[{rr_low}, {rr_high}]") ) |>select(`Tree cover`= tree_cov_label, `Count-scale slope`= count_est, `(95% CI)`= count_CI,`Rate ratio`= rr, `(95% CI) `= rr_CI, `p-value`= p.value)temp_tree_tbl_ex6 |>tt(caption ="Effect of temperature on abundance under drought (rain at its minimum), at low vs. high tree cover")
Effect of temperature on abundance under drought (rain at its minimum), at low vs. high tree cover
Tree cover
Count-scale slope
(95% CI)
Rate ratio
(95% CI)
p-value
Low tree cover (p10)
0.087
[0.058, 0.115]
1.018
[1.012, 1.024]
<0.001
High tree cover (p90)
0.166
[0.113, 0.22]
1.018
[1.012, 1.024]
<0.001
The rate ratio stays essentially constant across tree cover (it’s a fixed multiplicative effect, \[e^\beta\], independent of covariate values), while the count-scale slope grows substantially from low to high tree cover, since it’s rescaled by the much larger expected count \[\lambda\] there.
Case study: occurrence and abundance together
In reality, a species first has to be present before we can measure its abundance. A hurdle model combines two sub-models fit jointly:
A binomial (logit-link) component: does a plot clear the “hurdle” of zero? (occurrence)
A truncated count component: given it’s non-zero, how many individuals? (abundance when present)
Two ways to interpret: each part separately (what drives occurrence vs. abundance), or combined on the population scale (expected abundance across all plots = P(occurrence) × E[abundance | present]).
A simple hurdle example
# Simple hurdle simulation: abundance_pres = counts if present; occurence = presence;# abundance_pop = observed count in every plot (zero wherever the species is absent)tibble(tree_cov =generate(tree_cover, n_plots)) |>unnest(everything()) |>mutate(abundance_pres =rpois(n_plots, lambda = abundance_intercept + abundance_tree_effect * tree_cov ),occurence =rbinom(n_plots, size =1,prob = occurence_intercept + occurence_tree_effect * tree_cov),abundance_pop = abundance_pres * occurence) -> simple_pop_data# Hurdle model: main formula is the (zero-truncated Poisson) abundance-when-present part;# zi = ~ tree_cov models the probability of a structural zero (absence) with a logit linkmod_hu_simple <-glmmTMB(abundance_pop ~ tree_cov,zi =~ tree_cov,family = truncated_poisson,data = simple_pop_data)
Population vs. present-only predictions
Code
# type = "response": population-scale expected abundance across ALL plots (P(present) x count)predictions(mod_hu_simple, type ="response") |>mutate(type ="response") -> pred_resp_simple# type = "zprob": predicted probability of a structural zero (i.e. NOT occurring)predictions(mod_hu_simple, type ="zprob") -> pred_prob_simple# type = "conditional": expected abundance given the species is presentpredictions(mod_hu_simple, type ="conditional") |>mutate(type ="conditional") -> pred_cond_simplebind_rows(pred_resp_simple, pred_cond_simple) -> pred_abun_simpleplot_pred_abun_simple <-ggplot(pred_abun_simple) +geom_ribbon(aes(x = tree_cov, ymin = conf.low, ymax = conf.high, fill = type), alpha =0.2) +geom_line(aes(x = tree_cov, y = estimate, color = type), linewidth =1.5, lineend ="round") +labs(x ="Tree cover", y ="Abundance", color ="Prediction type", fill ="Prediction type")# zprob is P(structural zero), so 1 - zprob is the occurrence probabilityplot_pred_prob_simple <-ggplot(pred_prob_simple) +geom_ribbon(aes(x = tree_cov, ymin =1- conf.high, ymax =1- conf.low), fill ="white", alpha =0.2) +geom_line(aes(x = tree_cov, y =1- estimate), color ="white", linewidth =1.5, lineend ="round") +labs(x ="Tree cover", y ="Occurence") plot_pred_abun_simple / plot_pred_prob_simple
Top: predicted abundance, population-scale (response) vs. when present (conditional). Bottom: predicted occurrence probability (1 - zprob).
Predicted abundance when present (conditional) is always higher than population-scale abundance (response), since occurrence probability is below 100%. This is an implicit interaction between the model’s two parts: tree cover’s effect on population-scale abundance depends on its effect on both occurrence and abundance-when-present.
Simulating with an explicit interaction
Prey and competitors interact in their effect on abundance (when present); tree cover drives occurrence only. Predictors are simulated through the model’s own inverse link (exp(), plogis()) so they stay in a valid range:
zi = ~ scl_tree_cov: the zero-inflation (occurrence) part. The conditional formula’s scl_prey * scl_competitors is an explicit interaction, layered on top of the log link’s implicit one.
The interaction on the count scale
Code
slopes(mod_hu, variables ="scl_prey",newdata =datagrid(scl_competitors =seq(min(pop_data$scl_competitors),max(pop_data$scl_competitors), length.out =30)),type ="response") |># unscale back to raw competitor counts for an interpretable x-axismutate(competitors = scl_competitors *sd(pop_data$competitors) +mean(pop_data$competitors)) |>ggplot(aes(x = competitors, y = estimate)) +geom_hline(yintercept =0, color ="blue", linewidth =1, linetype =2) +geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha =0.4, fill ="lightgrey") +geom_line(linewidth =1.5, lineend ="round", color ="white") +labs(y ="Slope of scl_prey", x ="Competitors",title ="Count-scale effect of prey availability")
Count-scale effect of prey weakens (and reverses) as competitors increase.
The interaction on the link scale
Code
slopes(mod_hu, variables ="scl_prey",newdata =datagrid(scl_competitors =seq(min(pop_data$scl_competitors),max(pop_data$scl_competitors), length.out =30)),type ="link") |># unscale back to raw competitor counts for an interpretable x-axismutate(competitors = scl_competitors *sd(pop_data$competitors) +mean(pop_data$competitors)) |>ggplot(aes(x = competitors, y = estimate)) +geom_hline(yintercept =0, color ="blue", linewidth =1, linetype =2) +geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha =0.4, fill ="lightgrey") +geom_line(linewidth =1.5, lineend ="round", color ="white") +labs(y ="Slope of scl_prey", x ="Competitors",title ="Link-scale effect of prey")
Link-scale effect of prey is a straight, non-zero-crossing line in competitors — the explicit interaction, undistorted by the log link.
Unlike the count-scale slide, this one is linear in competitors by construction — the model’s formula (scl_prey * scl_competitors) is fit additively on the log scale, so the explicit interaction shows up here as a straight line, with no additional curvature from the log link itself.
Isolating the interaction: rate ratios
To isolate something closer to “the interaction itself,” we move to the rate-ratio scale: for each one-SD increase in competitors, by what factor does the effect of one extra SD of prey get multiplied? That’s exp() of the interaction coefficient, and because it lives on the model’s own additive (log) scale, it doesn’t depend on where competitors are evaluated — unlike the count-scale slope.
Code
rate_ratio_interaction_tbl <-avg_slopes(mod_hu, variables ="scl_prey",by ="scl_competitors",newdata =datagrid(scl_competitors =c(-1, 0, 1)),type ="link") |>as_tibble() |>mutate(rate_ratio =exp(estimate), rr_low =exp(conf.low), rr_high =exp(conf.high)) |>select(scl_competitors, rate_ratio, rr_low, rr_high, p.value) |>mutate(across(c(rate_ratio, rr_low, rr_high), ~round(.x, 3)),p.value =if_else(p.value <0.001, "<0.001", as.character(round(p.value, 3))) ) |>mutate(CI_95 = glue::glue("[{rr_low}, {rr_high}]")) |>select(`Competitors (scaled)`= scl_competitors, `Rate ratio for scl_prey`= rate_ratio, `95% CI`= CI_95, `p-value`= p.value)rate_ratio_interaction_tbl |>tt(caption ="Rate ratio of scl_prey's effect on abundance, evaluated at three levels of scl_competitors")
Rate ratio of scl_prey's effect on abundance, evaluated at three levels of scl_competitors
Competitors (scaled)
Rate ratio for scl_prey
95% CI
p-value
-1
1.227
[1.143, 1.318]
<0.001
0
0.953
[0.861, 1.055]
0.349
1
0.740
[0.615, 0.89]
0.001
Prey’s rate ratio drops from >1 (increasing abundance) at low competitor density to <1 (decreasing abundance) at high competitor density.
Predicted abundance, by competitor level
Code
predictions(mod_hu,newdata =datagrid(scl_prey =seq(min(pop_data$scl_prey), max(pop_data$scl_prey), length.out =30),scl_competitors = fivenum ),type ="conditional") |>mutate(prey = scl_prey *sd(pop_data$prey) +mean(pop_data$prey),competitors =round(scl_competitors *sd(pop_data$competitors) +mean(pop_data$competitors), 1)) |>ggplot() +geom_ribbon(aes(x = prey, ymin = conf.low, ymax = conf.high, fill =as.factor(competitors)), alpha =0.2) +geom_line(aes(x = prey, y = estimate, color =as.factor(competitors)), linewidth =1.5, lineend ="round") +labs(y ="Predicted abundance\n(when present)", x ="Prey", color ="Competitors", fill ="Competitors")
Abundance-when-present predictions across prey, one line per competitor five-number-summary level.
At low competitor density, abundance climbs steeply with prey; by the time competitors are abundant, the prey curve is nearly flat — the same shrinking slope we saw before, now shown as predicted trends rather than slopes.
Population-scale picture: low vs. high tree cover
Code
# Locate the two peaks of tree_cov's bimodal distribution, as scaled valuesd <-density(pop_data$tree_cov)peaks <- d$x[which(diff(sign(diff(d$y))) ==-2) +1]low_tree_pop <-round(peaks[1], 2)high_tree_pop <-round(peaks[length(peaks)], 2)low_tree_pop_scl <- (low_tree_pop -mean(pop_data$tree_cov)) /sd(pop_data$tree_cov)high_tree_pop_scl <- (high_tree_pop -mean(pop_data$tree_cov)) /sd(pop_data$tree_cov)predictions(mod_hu,newdata =datagrid(scl_prey =seq(min(pop_data$scl_prey), max(pop_data$scl_prey), length.out =30),scl_competitors = fivenum,scl_tree_cov =c(low_tree_pop_scl, high_tree_pop_scl) ),type ="response") |>mutate(prey = scl_prey *sd(pop_data$prey) +mean(pop_data$prey),competitors =round(scl_competitors *sd(pop_data$competitors) +mean(pop_data$competitors), 1),`Tree cover`=if_else(scl_tree_cov == low_tree_pop_scl, "Low tree cover", "High tree cover")) |>ggplot() +geom_ribbon(aes(x = prey, ymin = conf.low, ymax = conf.high, fill =as.factor(competitors)), alpha =0.2) +geom_line(aes(x = prey, y = estimate, color =as.factor(competitors)), linewidth =1.5, lineend ="round") +facet_wrap(~`Tree cover`) +labs(y ="Predicted population abundance", x ="Prey", color ="Competitors", fill ="Competitors")
Population-scale (response) predictions: prey x competitors, faceted by low vs. high tree cover.
Within each panel, the prey × competitors shape is the same. But the high-tree-cover panel sits above the low-tree-cover panel overall — more plots are occupied there. Tree cover mostly shifts how many plots have the species at all; prey and competitors jointly shape abundance within occupied plots.
A symmetric interaction, viewed from the other predictor
Code
predictions(mod_hu,newdata =datagrid(scl_competitors =seq(min(pop_data$scl_competitors), max(pop_data$scl_competitors), length.out =30),scl_prey = fivenum ),type ="conditional") |>mutate(competitors = scl_competitors *sd(pop_data$competitors) +mean(pop_data$competitors),prey =round(scl_prey *sd(pop_data$prey) +mean(pop_data$prey), 1)) |>ggplot() +geom_ribbon(aes(x = competitors, ymin = conf.low, ymax = conf.high, fill =as.factor(prey)), alpha =0.2) +geom_line(aes(x = competitors, y = estimate, color =as.factor(prey)), linewidth =1.5, lineend ="round") +labs(y ="Predicted abundance\n(when present)", x ="Competitors", color ="Prey", fill ="Prey")
Same interaction, axes swapped: competitors on x, one line per prey five-number-summary level.
A two-way interaction is mathematically symmetric: scl_prey:scl_competitors is the same coefficient as scl_competitors:scl_prey. This view makes a different point than the prey-on-x-axis plot: regardless of prey availability, enough competitors drive abundance down to roughly the same low floor. Neither plot is “more correct” — reporting only one risks showing half the story (Spake et al. 2023).
Abundance-when-present, averaged over the full observed range of the other predictor.
datagrid(..., grid_type = "counterfactual") pairs each focal value with every observed row, then averages — a true average over the sample rather than a slope evaluated at the other predictors’ means. Averaged over all observed competitor density, prey still shows a positive effect; averaged over all observed prey, competitors show a diminishing shape.
Population-scale average, by tree cover
Code
avg_pred_prey_resp <-avg_predictions(mod_hu,newdata =datagrid(model = mod_hu,scl_prey =seq(min(pop_data$scl_prey), max(pop_data$scl_prey), length.out =30),scl_tree_cov =c(low_tree_pop_scl, high_tree_pop_scl),grid_type ="counterfactual"),by =c("scl_prey", "scl_tree_cov"), type ="response") |>mutate(prey = scl_prey *sd(pop_data$prey) +mean(pop_data$prey),`Tree cover`=if_else(scl_tree_cov == low_tree_pop_scl, "Low tree cover", "High tree cover"))avg_pred_comp_resp <-avg_predictions(mod_hu,newdata =datagrid(model = mod_hu,scl_competitors =seq(min(pop_data$scl_competitors), max(pop_data$scl_competitors), length.out =30),scl_tree_cov =c(low_tree_pop_scl, high_tree_pop_scl),grid_type ="counterfactual"),by =c("scl_competitors", "scl_tree_cov"), type ="response") |>mutate(competitors = scl_competitors *sd(pop_data$competitors) +mean(pop_data$competitors),`Tree cover`=if_else(scl_tree_cov == low_tree_pop_scl, "Low tree cover", "High tree cover"))p_avg_prey_resp <- avg_pred_prey_resp |>ggplot(aes(x = prey, y = estimate, color =`Tree cover`, fill =`Tree cover`)) +geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha =0.2, color =NA) +geom_line(linewidth =1.5, lineend ="round") +labs(y ="Average predicted\npopulation abundance", x ="Prey")p_avg_comp_resp <- avg_pred_comp_resp |>ggplot(aes(x = competitors, y = estimate, color =`Tree cover`, fill =`Tree cover`)) +geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha =0.2, color =NA) +geom_line(linewidth =1.5, lineend ="round") +labs(y ="Average predicted\npopulation abundance", x ="Competitors")p_avg_prey_resp + p_avg_comp_resp
Population-scale average curves, one pair per tree-cover peak.
The shapes barely change between tree-cover panels — the prey × competitors interaction lives entirely in the abundance-when-present part. What changes is the overall height: high-tree-cover populations sit consistently above low-tree-cover ones, a pure occupancy effect, cleanly separated from the abundance-when-present shape.
A full rate-ratio table, by model component
Code
# Conditional-model rate ratio for scl_prey, at three levels of scl_competitorsprey_rr_tbl <-avg_slopes(mod_hu, variables ="scl_prey", by ="scl_competitors",newdata =datagrid(scl_competitors =c(-1, 0, 1)),type ="link") |>mutate(rate_ratio =exp(estimate), rr_low =exp(conf.low), rr_high =exp(conf.high),Term = glue::glue("scl_prey (at scl_competitors = {scl_competitors} SD)")) |>select(Term, rate_ratio, rr_low, rr_high, p.value)# Conditional-model rate ratio for scl_competitors, at three levels of scl_preycomp_rr_tbl <-avg_slopes(mod_hu, variables ="scl_competitors", by ="scl_prey",newdata =datagrid(scl_prey =c(-1, 0, 1)),type ="link") |>mutate(rate_ratio =exp(estimate), rr_low =exp(conf.low), rr_high =exp(conf.high),Term = glue::glue("scl_competitors (at scl_prey = {scl_prey} SD)")) |>select(Term, rate_ratio, rr_low, rr_high, p.value)# The interaction itself: ratio of rate ratios per 1-SD increase in the moderatorhu_confint <-confint(mod_hu)hu_coef_cond <-summary(mod_hu)$coefficients$condhu_coef_zi <-summary(mod_hu)$coefficients$ziinteraction_rr_tbl <-tibble(Term ="scl_prey:scl_competitors (interaction)",rate_ratio =exp(hu_coef_cond["scl_prey:scl_competitors", "Estimate"]),rr_low =exp(hu_confint["cond.scl_prey:scl_competitors", "2.5 %"]),rr_high =exp(hu_confint["cond.scl_prey:scl_competitors", "97.5 %"]),p.value = hu_coef_cond["scl_prey:scl_competitors", "Pr(>|z|)"])# Zero-inflation-model odds ratio for scl_tree_cov (no moderator in this model)zi_rr_tbl <-tibble(Term ="scl_tree_cov (zero-inflation)",rate_ratio =exp(hu_coef_zi["scl_tree_cov", "Estimate"]),rr_low =exp(hu_confint["zi.scl_tree_cov", "2.5 %"]),rr_high =exp(hu_confint["zi.scl_tree_cov", "97.5 %"]),p.value = hu_coef_zi["scl_tree_cov", "Pr(>|z|)"])hu_rr_tbl <-bind_rows( prey_rr_tbl |>mutate(Component ="Conditional (abundance, when present)"), comp_rr_tbl |>mutate(Component ="Conditional (abundance, when present)"), interaction_rr_tbl |>mutate(Component ="Conditional (abundance, when present)"), zi_rr_tbl |>mutate(Component ="Zero-inflation (occupancy, as an odds ratio)")) |>relocate(Component) |>mutate(across(c(rate_ratio, rr_low, rr_high), ~round(.x, 3)),p.value =if_else(p.value <0.001, "<0.001", as.character(round(p.value, 3))),CI_95 = glue::glue("[{rr_low}, {rr_high}]") ) |>select(Component, Term, `Rate/odds ratio`= rate_ratio, `95% CI`= CI_95, `p-value`= p.value)hu_rr_tbl |>tt(caption ="Standardized rate ratios (conditional model) and odds ratio (zero-inflation model) for mod_hu") |>style_tt(i =which(hu_rr_tbl$Component != dplyr::lag(hu_rr_tbl$Component, default ="")) -1,j =1:5, line ="b", line_color ="black", line_width =0.05) |>style_tt(fontsize =0.6)
Standardized rate ratios (conditional model) and odds ratio (zero-inflation model) for mod_hu
Component
Term
Rate/odds ratio
95% CI
p-value
Conditional (abundance, when present)
scl_prey (at scl_competitors = -1 SD)
1.227
[1.143, 1.318]
<0.001
Conditional (abundance, when present)
scl_prey (at scl_competitors = 0 SD)
0.953
[0.861, 1.055]
0.349
Conditional (abundance, when present)
scl_prey (at scl_competitors = 1 SD)
0.740
[0.615, 0.89]
0.001
Conditional (abundance, when present)
scl_competitors (at scl_prey = -1 SD)
0.650
[0.565, 0.748]
<0.001
Conditional (abundance, when present)
scl_competitors (at scl_prey = 0 SD)
0.505
[0.456, 0.559]
<0.001
Conditional (abundance, when present)
scl_competitors (at scl_prey = 1 SD)
0.392
[0.341, 0.45]
<0.001
Conditional (abundance, when present)
scl_prey:scl_competitors (interaction)
0.776
[0.705, 0.855]
<0.001
Zero-inflation (occupancy, as an odds ratio)
scl_tree_cov (zero-inflation)
0.898
[0.81, 0.996]
0.041
Prey’s rate ratio crosses 1 somewhere around average competitor density (why its raw summary() estimate looked like noise); competitors suppress abundance more as prey increases; the interaction row (~0.81) is the single number behind both patterns, and it’s the same number whichever predictor you treat as the moderator — the symmetry point from a few slides back, visible directly in the table. Tree cover’s odds ratio in the zero-inflation part tells a simpler, interaction-free story about occupancy alone.
Part 3 takeaways
Non-linear link functions (logit, log) create an implicit interaction: every predictor’s response-scale effect depends on the values of all other predictors, even with no interaction term in the formula
Raw GLM coefficients live on the link scale; marginaleffects::slopes() converts them to the natural (probability, count) scale for you
Multiplicative summaries (rate ratios, odds ratios) stay constant across the range of other predictors; additive, response-scale summaries do not
Whether an interaction “exists” can depend on the scale you’re measuring on (Spake et al. 2023) — this applies just as much to implicit link-function interactions as to explicit interaction terms
Hurdle/zero-inflated models separate “does it occur” from “how much, given it occurs” — interpret each part separately, or combine them (type = "response") for a population-scale answer
Interactions are mathematically symmetric: report (or at least check) both “A moderated by B” and “B moderated by A”
Exercise 6 / Homework
Now apply what we learned today to your data!
Wrap-up and further reading
Full write-ups with additional detail and 3D visualizations: Part 1, Part 2, Part 3
