Interpretation of interaction plot from GLM

Question

I'm doing a moderation analysis in R. My dependent variable is a rate (Diffcount/totalcount), the independent variable is an index (continuous), and the moderator is a categorical variable (3 levels). I fitted a negative binomial regression model using glm.nb from MASS, and I used offset term like below:

glm.nb(Diffcount~Index1*factor3 + offset(log(totalcount)), data = dt)

I generated a interaction plot for this model, using interact_plot from interactions package. What I'm really confusing about is the interpretation of the y-axis. Does the value represent the rate (i.e., Diffcount/totalcount), a probability, or maybe something else? Thank you!

Answer 1

The y axis is the conditional expectation. It is interpreted depending on what you pass as covariates to the predict method for glm.nb. Let's see a small example.

I'm going to simulate data from

$$\log(E[y\mid x]) = -2 + 0.8x + 0.2z -0.3xz + \log(n) $$

for binary $z$, which is similar to your model. I will simulate the data to be from a poisson distribution despite us using the engative binomial.

library(tidyverse)

x <- seq(-2, 2, 0.1)
z <- letters[1:2]

dgrid <- crossing(x, z)
dgrid$n <- rpois(nrow(dgrid), 10)

X <- model.matrix(~x*z, data=dgrid)
eta <- X %*% c(-2, 0.8, 0.2, -0.3) + log(dgrid$n)
dgrid$y <- rpois(nrow(dgrid), exp(eta))

Let's now plot our predictions. To do so, we need to specify values of $x$, $z$, and $n$. Let's set $z=a$ and then compute the prediction over a grid of $x$ values for 2 values of $n$ (2 and 10). Plot is shown below

fit <- MASS::glm.nb(y~x*z + offset(log(n)), data=dgrid)
dpred <- crossing(x, z='a', n=c(2, 10))
dpred$mu <- predict(fit, newdata=dpred, type='response')

dpred %>% 
  ggplot(aes(x, mu, color=factor(n))) + 
  geom_line()

Two different curves. The curve depends on what the value of $n$ is -- generally, larger $n$ means larger values of the prediction.

Now, onto your questions.

What I'm really confusing about is the interpretation of the y-axis.

The y axis here and in your plot is the expected value for the outcome. It isn't a rate per se, it is an expected count variable. Why should prediction be a floating point number and not an integer when the outcome is an integer? It is the expected count, and it is more than reasonable to talk about the expected outcome being a float.

To get the expected rate, set $n=1$ when making the prediction call, as follows

dpred <- crossing(x, z='a', n=1)
dpred$mu <- predict(fit, newdata=dpred, type='response')

I imagine the library you've used to construct that plot might be setting $n=1$, but I can't be sure. If you can, verify this is actually the case before proceeding. Else, you can use a different library which makes this more transparent, like {marginaleffects}. Here is how I would make your plot in that library

library(marginaleffects)

plot_predictions(
  fit, 
  newdata = datagrid(n=1, x=seq(-2, 2, 0.1), by='z'),
  by=c('x', 'z')
)

With {marginaleffects}, you can specify the value of $n$ yourself. If you don't like the aesthetics of the plot, you can use the predictions function to get the data that would be used to create this plot.