0

The dataset I'm working with has an interaction term. I want to fit the model with the interaction term in the x-axis and the y.var in the y-axis. I tried following this example but didn't get too far on how to replicate it in ggplot2 (i.e. the plot style function doesn't work anymore..so I didn't know how to replicate the results).

temp <- rnorm(100, 2,1)set.seed(111)
temp <- rnorm(100, 3,1)
rainfall <- rnorm(100,5,1)
y.var <- rnorm(100, 2,1)
site <- rep(c("A","B","C","D"), each = 25)     

df <- data.frame(temp, rainfall, y.var, site)
df$site <- as.factor(as.character(df$site))

mod <- lmer(y.var ~ temp * rainfall + (1|site), data  = df)
summary(mod)
Rspacer
  • 2,369
  • 1
  • 14
  • 40
  • You want to replicate it in `ggplot`? The `plot_model` function does work. – Quinten May 17 '22 at 11:58
  • `ggplot2` is for plotting data. To plot model coefficients or model effects, I recommend `ggeffects` and `sjPlot`. Also, `lm` is not a mixed effects model, it's ordinary least squares. You'll want to use `lme4::lmer`. – Arthur May 17 '22 at 12:06
  • @Quinten `plot_model(fit, type = "pred", terms = c("temp:rainfall"))` Error in plot_model(fit, type = "pred", terms = c("temp:rainfall")) : could not find function "plot_model" – Rspacer May 17 '22 at 12:07
  • @Arthur See my updated df and code. I tried `ggPredict(model, c("temp * rainfall", "site"))` and I get an error that says... Error: $ operator not defined for this S4 class – Rspacer May 17 '22 at 12:18
  • Using the `ggeffects` package, `ggpredict(mod, terms=c("temp", "rainfall")) %>% plot()` should work. – DaveArmstrong May 17 '22 at 12:46

1 Answers1

4

If you want to plot the model in ggplot directly rather than using an extension package, you need to generate a data frame of predictions to plot. The benefit of doing it this way means you can also include your original data points on the plot.

Since you have y.var on the y axis, you are only left with one axis to plot two fixed effect variables. This means you will need to choose either rainfall or temperature for the x axis, and represent the other variable with another aesthetic such as color. In this example I'll use temperature for the x axis. Obviously to make the plot interpretable, we need to limit the number of "slices" of rainfall that we predict. Here I'll use 5.

The interaction effect is visible here as the change in slope of the line as rainfall increases.

library(geomtextpath)
library(lme4)

mod <- lmer(y.var ~ temp * rainfall + (1|site), data  = df)

newdf <- expand.grid(temp = seq(min(df$temp), max(df$temp), length = 100),
                     rainfall = seq(min(df$rainfall), max(df$rainfall), 
                                    length = 5))

newdf$y.var <- predict(mod, newdata = newdf)

ggplot(newdf, aes(x = temp, y = y.var, group = rainfall)) +
  geom_point(data = df, aes(shape = site, color = rainfall)) +
  geom_textline(aes(color = rainfall, label = round(rainfall, 2)), 
                hjust = 0.95) +
  scale_color_gradient(low = 'navy', high = 'red4') +
  theme_light(base_size = 16)

enter image description here

Allan Cameron
  • 147,086
  • 7
  • 49
  • 87
  • Great! Very helpful. Is there a way to specify the moisture levels. Say I want to only plot 1,3,5, and 7 or does the model preassign it? – Rspacer May 17 '22 at 15:20
  • 1
    @Biotechgeek the way to assign them is in the `expand.grid` call. Just change `rainfall = seq(..., length = 5)` to `rainfall = c(1, 3, 5, 7)` – Allan Cameron May 17 '22 at 15:27