2 Interaction term in a non-linear model

Published

December 7, 2017

2.1 The issue

In a linear model with interaction,

\[ E(y) = \beta_1 x_1 + \beta_2 x_2 + \beta_{12} x_1 x_2, \tag{2.1}\]

the cross partial derivative is simply the coefficient:

\[ \frac{\partial^2 E(y)}{\partial x_1 \partial x_2} = \beta_{12}. \tag{2.2}\]

In a nonlinear model — logit, Poisson, or any link function $F$ — the same coefficient $\beta_{12}$ sits on the link scale:

\[ F(E(y)) = \beta_1 x_1 + \beta_2 x_2 + \beta_{12} x_1 x_2. \tag{2.3}\]

The cross partial of $F(E(y))$ is still $\beta_{12}$. But the cross partial of $E(y)$ itself — the quantity most researchers want — is not:

\[ \frac{\partial^2 E(y)}{\partial x_1 \partial x_2} = \beta_{12} G'(\eta) + (\beta_{1} + \beta_{12} x_2)(\beta_2 + \beta_{12} x_1) G''(\eta), \tag{2.4}\]

where $G = F^{-1}$ and $\eta$ is the linear predictor. This is the point of Ai and Norton (2003): the probability-scale interaction effect depends on every coefficient, varies across the covariate space, and can differ in sign and significance from $\beta_{12}$.

The deeper issue is that nonlinearity alone generates a nonzero cross partial — even without an interaction term. Drop $\beta_{12}$ and the cross partial on $E(y)$ is still $\beta_1 \beta_2 G''(\eta) \ne 0$ whenever the link is curved. The marginal effect of $x_1$ depends on $x_2$ in any nonlinear model, interaction or not.

Greene (Econometric Analysis, 8th ed., sec. 17.3) and Buis (2010) draw the same conclusion: the model is specified on the link scale, and that is where coefficients and their tests are interpretable. In a logit, $\beta_{12}$ is the interaction on log-odds; exponentiating gives an odds-ratio interaction that is multiplicative. Partial effects on the probability scale are implications of the model, not parameters of it — they should be computed from the fitted model (via margins or marginaleffects), not read off the coefficient table.

2.2 A worked example

Let’s illustrate the Ai-Norton point directly: in a logit model with an interaction term, the cross partial derivative $\frac{\partial^2 E(y)}{\partial x_1 \partial x_2}$ (what we actually care about) is not the same object as the estimated coefficient $\beta_{12}$ on $x_1 x_2$ in the linear index. We simulate data with a true interaction effect on the log-odds scale, fit a logit model, and then use the marginaleffects package for two distinct things: first the marginal effect of $x_1$ on the probability scale at several values of $x_2$, and then the cross partial derivative itself.

We draw $n = 2000$ observations with $x_1, x_2 \sim N(0,1)$, independent of each other. The linear index is $\eta = -0.3 + 0.8x_1 + 0.5x_2 + 0.6\,x_1x_2$, the success probability is $p = \Lambda(\eta)$, and $y \sim \text{Bernoulli}(p)$. So the true interaction on the log-odds scale is $\beta_{12} = 0.6$, a single constant — which is exactly the object the coefficient table estimates, and exactly the object that does not transfer to the probability scale.

Code

library(dplyr)
library(marginaleffects)

set.seed(42)
n <- 2000
x1 <- rnorm(n)
x2 <- rnorm(n)
eta <- -0.3 + 0.8 * x1 + 0.5 * x2 + 0.6 * x1 * x2
p <- plogis(eta)
y <- rbinom(n, 1, p)
dat <- data.frame(y, x1, x2)

m <- glm(y ~ x1 * x2, data = dat, family = binomial)
summary(m)$coefficients

              Estimate Std. Error   z value     Pr(>|z|)
(Intercept) -0.2348658 0.05017155 -4.681255 2.851237e-06
x1           0.7792602 0.05930028 13.140920 1.918795e-39
x2           0.5296149 0.05537644  9.563903 1.133977e-21
x1:x2        0.6624955 0.06403580 10.345705 4.376826e-25

The coefficient on x1:x2 is 0.662 with a standard error of 0.064, against a true 0.6 — recovered to within one standard error, as it should be, since the logit is correctly specified. The other three are 0.779, 0.530 and $-0.235$ against true values of 0.8, 0.5 and $-0.3$.

That coefficient is $\hat\beta_{12}$, the interaction effect on the log-odds (linear index) scale — constant by construction. On the probability scale nothing is constant. Start with the marginal effect of $x_1$ itself, evaluated at three values of $x_2$:

Code

slopes(m, variables = "x1", by = "x2", newdata = datagrid(x2 = c(-1, 0, 1)))


 x2 Estimate Std. Error     z Pr(>|z|)     S    2.5 % 97.5 %
 -1   0.0253     0.0159  1.59    0.112   3.2 -0.00591 0.0565
  0   0.1919     0.0147 13.09   <0.001 127.6  0.16314 0.2206
  1   0.3538     0.0243 14.59   <0.001 157.7  0.30630 0.4014

Term: x1
Type: response
Comparison: dY/dX

The marginal effect of $x_1$ on $P(y=1)$ is 0.025 at $x_2=-1$, 0.192 at $x_2=0$, and 0.354 at $x_2=1$ — a fourteen-fold change across the three values, where the log-odds interaction is one fixed number. Notice too that at $x_2=-1$ the effect is not significant ($p = 0.112$) even though $\hat\beta_{12}$ is significant at better than $10^{-24}$. Already the probability-scale answer and the coefficient tell different stories.

The cross partial derivative is a further step: not the effect of $x_1$ at a given $x_2$, but how much that effect changes as $x_2$ moves. That is a difference of slopes, so it needs a hypothesis test on top of the by argument rather than a second call to slopes():

Code

slopes(m, variables = "x1", by = "x2",
       newdata = datagrid(x2 = c(-1, 1)),
       hypothesis = "b2 - b1 = 0")


 Hypothesis Estimate Std. Error  z Pr(>|z|)    S 2.5 % 97.5 %
    b2-b1=0    0.329     0.0299 11   <0.001 90.7  0.27  0.387

Type: response

The cross partial is 0.329 with a standard error of 0.030, $z = 11$ — that is the difference between the slope at $x_2=1$ and the slope at $x_2=-1$, $0.354 - 0.025$, with a proper standard error attached. It is large and highly significant, which is worth being honest about: in this example it agrees with $\hat\beta_{12}$ in significance, differing only in scale and magnitude. Note that it is nowhere near 0.662; the two numbers live on different scales and there is no reason to expect them to match. Ai and Norton’s warning is that the two quantities need not agree — not that they never do. The dependable lesson is the one visible in the table above: on the probability scale the answer depends on where you evaluate it, so a single coefficient cannot summarise it.

--- title: "Interaction term in a non-linear model" date: "2017-12-07" --- ## The issue In a linear model with interaction, $$ E(y) = \beta_1 x_1 + \beta_2 x_2 + \beta_{12} x_1 x_2, $$ {#eq-interaction-nonlinear-1} the cross partial derivative is simply the coefficient: $$ \frac{\partial^2 E(y)}{\partial x_1 \partial x_2} = \beta_{12}. $$ {#eq-interaction-nonlinear-2} In a nonlinear model — logit, Poisson, or any link function $F$ — the same coefficient $\beta_{12}$ sits on the *link* scale: $$ F(E(y)) = \beta_1 x_1 + \beta_2 x_2 + \beta_{12} x_1 x_2. $$ {#eq-interaction-nonlinear-3} The cross partial of $F(E(y))$ is still $\beta_{12}$. But the cross partial of $E(y)$ itself — the quantity most researchers want — is not: $$ \frac{\partial^2 E(y)}{\partial x_1 \partial x_2} = \beta_{12} G'(\eta) + (\beta_{1} + \beta_{12} x_2)(\beta_2 + \beta_{12} x_1) G''(\eta), $$ {#eq-interaction-nonlinear-5} where $G = F^{-1}$ and $\eta$ is the linear predictor. This is the point of [Ai and Norton (2003)](https://pdfs.semanticscholar.org/6285/8e64d9a337504d72cb862c4cc1e7fd27a7a0.pdf): the probability-scale interaction effect depends on every coefficient, varies across the covariate space, and can differ in sign and significance from $\beta_{12}$. The deeper issue is that nonlinearity alone generates a nonzero cross partial — even without an interaction term. Drop $\beta_{12}$ and the cross partial on $E(y)$ is still $\beta_1 \beta_2 G''(\eta) \ne 0$ whenever the link is curved. The marginal effect of $x_1$ depends on $x_2$ in any nonlinear model, interaction or not. Greene (*Econometric Analysis*, 8th ed., sec. 17.3) and [Buis (2010)](http://www.stata-journal.com/sjpdf.html?articlenum=st0194) draw the same conclusion: the model is specified on the link scale, and that is where coefficients and their tests are interpretable. In a logit, $\beta_{12}$ is the interaction on log-odds; exponentiating gives an odds-ratio interaction that is multiplicative. Partial effects on the probability scale are *implications* of the model, not parameters of it — they should be computed from the fitted model (via `margins` or `marginaleffects`), not read off the coefficient table. ## A worked example Let's illustrate the Ai-Norton point directly: in a logit model with an interaction term, the cross partial derivative $\frac{\partial^2 E(y)}{\partial x_1 \partial x_2}$ (what we actually care about) is not the same object as the estimated coefficient $\beta_{12}$ on $x_1 x_2$ in the linear index. We simulate data with a true interaction effect on the log-odds scale, fit a logit model, and then use the `marginaleffects` package for two distinct things: first the marginal effect of $x_1$ on the probability scale at several values of $x_2$, and then the cross partial derivative itself. We draw $n = 2000$ observations with $x_1, x_2 \sim N(0,1)$, independent of each other. The linear index is $\eta = -0.3 + 0.8x_1 + 0.5x_2 + 0.6\,x_1x_2$, the success probability is $p = \Lambda(\eta)$, and $y \sim \text{Bernoulli}(p)$. So the true interaction on the log-odds scale is $\beta_{12} = 0.6$, a single constant — which is exactly the object the coefficient table estimates, and exactly the object that does *not* transfer to the probability scale. ```{r} #| label: ai-norton-example #| cache: true #| warning: false #| message: false library(dplyr) library(marginaleffects) set.seed(42) n <- 2000 x1 <- rnorm(n) x2 <- rnorm(n) eta <- -0.3 + 0.8 * x1 + 0.5 * x2 + 0.6 * x1 * x2 p <- plogis(eta) y <- rbinom(n, 1, p) dat <- data.frame(y, x1, x2) m <- glm(y ~ x1 * x2, data = dat, family = binomial) summary(m)$coefficients ``` The coefficient on `x1:x2` is 0.662 with a standard error of 0.064, against a true 0.6 — recovered to within one standard error, as it should be, since the logit is correctly specified. The other three are 0.779, 0.530 and $-0.235$ against true values of 0.8, 0.5 and $-0.3$. That coefficient is $\hat\beta_{12}$, the interaction effect on the log-odds (linear index) scale — constant by construction. On the *probability* scale nothing is constant. Start with the marginal effect of $x_1$ itself, evaluated at three values of $x_2$: ```{r} #| label: ai-norton-slopes #| cache: true #| warning: false #| message: false slopes(m, variables = "x1", by = "x2", newdata = datagrid(x2 = c(-1, 0, 1))) ``` The marginal effect of $x_1$ on $P(y=1)$ is 0.025 at $x_2=-1$, 0.192 at $x_2=0$, and 0.354 at $x_2=1$ — a fourteen-fold change across the three values, where the log-odds interaction is one fixed number. Notice too that at $x_2=-1$ the effect is not significant ($p = 0.112$) even though $\hat\beta_{12}$ is significant at better than $10^{-24}$. Already the probability-scale answer and the coefficient tell different stories. The cross partial derivative is a *further* step: not the effect of $x_1$ at a given $x_2$, but how much that effect changes as $x_2$ moves. That is a difference of slopes, so it needs a hypothesis test on top of the `by` argument rather than a second call to `slopes()`: ```{r} #| label: ai-norton-cross-partial #| cache: true #| warning: false #| message: false slopes(m, variables = "x1", by = "x2", newdata = datagrid(x2 = c(-1, 1)), hypothesis = "b2 - b1 = 0") ``` The cross partial is 0.329 with a standard error of 0.030, $z = 11$ — that is the difference between the slope at $x_2=1$ and the slope at $x_2=-1$, $0.354 - 0.025$, with a proper standard error attached. It is large and highly significant, which is worth being honest about: in *this* example it agrees with $\hat\beta_{12}$ in significance, differing only in scale and magnitude. Note that it is nowhere near 0.662; the two numbers live on different scales and there is no reason to expect them to match. Ai and Norton's warning is that the two quantities *need not* agree — not that they never do. The dependable lesson is the one visible in the table above: on the probability scale the answer depends on where you evaluate it, so a single coefficient cannot summarise it.