plot predictions of a model with an interaction and a poly term

plot_poly_fits.R

library(tidyverse)
library(ggplot2)
library(visreg)

#some fake data
set.seed(35)
my_data <- data.frame(x1 = runif(100,-50,50), x2 = runif(100, -50, 50)) %>%
  mutate(y = rnorm(100, 0.001*x2*(x1 - x1^2), 100))

#let's see that data
ggplot(my_data,
       aes(x=x1, y=y, color=x2)) +
  geom_point()
  
#fit a model
mod <- lm(y ~ x2*poly(x1,2), data=my_data)
summary(mod)

#visualize
visreg(mod, xvar = "x1", by = "x2", gg=TRUE)
visreg(mod, xvar = "x2", by = "x1", gg=TRUE)

#if visreg didn't work
library(modelr)
#make a data with lots of x1
new_data <- my_data %>%
  data_grid(x1 = seq_range(x1, 100),
            x2 = seq_range(x2, 3))

#get predicted levels and CIs 
#note, this works for metafor objects, too
pred_data <- predict(mod, newdata = new_data,
                     interval="confidence") %>%
  as_tibble() %>%
  rename(y = fit) %>%
  bind_cols(new_data)

#now plot
ggplot(pred_data,
       aes(x=x1, y=y, ymin=lwr, ymax=upr)) +
  geom_ribbon(alpha=0.5) +
  geom_line() +
  facet_wrap(~x2)