How can I use lm() with poly() and ggeffect() to plot model effects when variable in quadratic term is not being plotted

Question

I am trying to plot model effects using the ggeffects package. Usually this is no problem. However, I have a model in which there is an interaction between a quadratic term and a first order term (Ex: var_a^2:var_b). In my actual model there are two quadratic terms and one is involved in an interaction with a third term. The example below shows the issue with two variables to keep the example as simple as possible.

If I use the poly() function as follows everything works fine; the model works fine and I can plot model effects. However, I would like to be able to remove he var_a^2:var_b interaction from the model.

# Create a data frame with explanatory variables
dat <- expand.grid(var_a = c(-1, 0, 1),
                   var_b = c(-1, 0, 1))

# Replicate so there plenty of data points for modelling
dat <- bind_rows(dat, dat)

# Create response value
dat$resp <- dat$var_a + 4 * dat$var_a^2 + dat$var_a * dat$var_b + rnorm(nrow(dat), 0, 0.2)


mod1 <- lm(resp ~ poly(var_a, 2) * var_b, data = dat)
summary(mod1)

plot(ggeffect(mod1, c('var_a', 'var_b')))
plot(ggeffect(mod1, c('var_b')))

Here's the plot of the var_a:var_b interaction

Here's a plot of effect var_b alone. I wouldn't normally plot var_b without var_a since they are involved in an interaction, but this is just to show my issue. It works here, but not below.

As stated above, I would like to remove the var_a^2:var_b interaction. As far as I know, to do that I need to separate the first and second order poly() terms as done below.

# This is the same as mod 1 without the var_a^2:var_b interaction
# I've split up the poly() terms so I can remove any unwanted terms
mod2 <- lm(resp ~ poly(var_a, 1) * var_b + poly(var_a, 2)[, 2], data = dat)
summary(mod2)

I can plot the var_a:var_b effect and get the essentially the same result as above.

plot(ggeffect(mod2, c('var_a', 'var_b')))

But I can no longer plot the var_b effect alone.

plot(ggeffect(mod2, 'var_b'))

For what it's worth, I believe the ggeffect() function is not working because the effects package on which it relies sets all variables other than the one being plotted (var_b) to a constant value which does not allow the quadratic effect of var_a to be calculated.

Is there a way to use the poly() function in a way that allows me to remove unwanted interactions from the model while still being able to use ggeffect() to plot model effects?

I tried to add ggeffects and effects tags to this question but don't have enough points to create a new tag. If anyone would like to do that, please do.

Answer 1

In your particular example, you could create the 2nd poly-variable before fitting the model:

library(ggeffects)

dat <- expand.grid(var_a = c(-1, 0, 1),
                   var_b = c(-1, 0, 1))

# Replicate so there plenty of data points for modelling
dat <- rbind(dat, dat)

# Create response value
dat$resp <- dat$var_a + 4 * dat$var_a^2 + dat$var_a * dat$var_b + rnorm(nrow(dat), 0, 0.2)

dat$var_a_p1 <- poly(dat$var_a, 2)[,2]
mod2 <- lm(resp ~ poly(var_a, 1) * var_b + var_a_p1, data = dat)
plot(ggemmeans(mod2, c('var_b')))
#> Loading required namespace: emmeans
#> NOTE: Results may be misleading due to involvement in interactions
#> Loading required namespace: ggplot2

^{Created on 2020-12-08 by the reprex package (v0.3.0)}