How can I plot fitted models with observations whe

Question

I am learning how to apply a multilevel modelling to traditional ANOVA research designs. I would like to plot fitted regression lines for each treatment using ggplot2. I would like to plot regression lines based on models I fitted, not let ggplot2 to plot since I would like to see how estimations differ based on changing models. I am aware that I can calculate coefficients and slopes by myself, but because the model is relatively complicated, I am looking for ways to plot more easily.

This is the sample code for the type of research design I am dealing with. I found that the sjPlot package (http://www.strengejacke.de/sjPlot/sjp.lm/) provides very nice plots, which shows the regression lines for each treatment over each test timing and actual observations on a scatter plot. This is exactly what I would like to make with ggplot.

 require(tidyverse)
 require(sjPlot)
 set.seed (100)
 dat <- data_frame(
    participant_id = c(c(1:15), c(1:15)),
    treatment = c(sample (letters [1:3], 15, replace = T), sample (letters [1:3], 15, replace = T)),
    test_timing = c(sample(letters [1:3], 15, replace = T),sample(letters [1:3], 15, replace = T)),
    learning_gain = (runif(30, min = 1, max = 20))
)

fit <- lm (learning_gain ~ treatment * test_timing -1, data = dat)
sjp.lm(fit, type = "pred", 
vars = c("test_timing", "treatment"),facet.grid = F)

Or, something like this:

 sjp.lm(fit, type = "pred", 
    vars = c("test_timing", "treatment"),facet.grid = T)

I would really appreciate it if you could teach me how I can make a plot similar to this image using the ggplot2 package. Thank you!

Answer 1

Indeed, sjPlot just calls a function from ggeffects. ggeffects again just returns a data frame with the marginal effects / predicted values. So you can use this data frame to build your own ggplot-objects, or simply use the plot()-function to create a plot (a ggplot-object).

There is a package vignette that describes how to build own plots using ggplot: https://cran.r-project.org/web/packages/ggeffects/vignettes/marginaleffects.html

Furthermore, the help-files (see https://cran.r-project.org/web/packages/ggeffects/ggeffects.pdf) are full of examples on how to to create plots using ggplot. Maybe this helps you?

Here is an example for your specific case, but note that you will not see straight lines, because your variables are categorical:

require(tidyverse)
require(ggeffects)
set.seed (100)
dat <- data_frame(
  participant_id = c(c(1:15), c(1:15)),
  treatment = c(sample (letters [1:3], 15, replace = T), sample (letters [1:3], 15, replace = T)),
  test_timing = c(sample(letters [1:3], 15, replace = T),sample(letters [1:3], 15, replace = T)),
  learning_gain = (runif(30, min = 1, max = 20))
)

# fit model
fit <- lm (learning_gain ~ treatment * test_timing -1, data = dat)

# compute marginal effects
me <- ggpredict(fit, c("test_timing", "treatment"))

# see results
me

# A tibble: 9 x 5
      x predicted conf.low conf.high group
  <dbl>     <dbl>    <dbl>     <dbl> <fct>
1  1.00      9.81   -0.472      20.1 a    
2  1.00      9.29    4.15       14.4 b    
3  1.00     13.2     5.95       20.5 c    
4  2.00     14.1     8.93       19.2 a    
5  2.00     11.6     7.00       16.2 b    
6  2.00      7.23   -3.06       17.5 c    
7  3.00     13.7     7.77       19.6 a    
8  3.00     13.5     8.33       18.6 b    
9  3.00     10.9     6.66       15.1 c    

Line plot

# note that predictors are categorical, so no straight line
ggplot(me, aes(x = x, y = predicted, colour = group)) +
  geom_line()

Dot plot

# plot dots
ggplot(me, aes(x = x, y = predicted, colour = group)) +
  geom_point()

Dot plots with error bars

# add CI, need jitter
ggplot(me, aes(x = x, y = predicted, colour = group)) +
  geom_point(position = position_dodge(.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(.2))

Facets

# facets
ggplot(me, aes(x = x, y = predicted, colour = group)) +
  geom_point(position = position_dodge(.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(.2)) +
  facet_grid(~group)

plot()-function from the ggeffects package

# the simple way
plot(me)

Raw Observations

ggpredict() returns the raw data as attribute to the return value, so you can also plot the original observations as well:

raw <- attr(me, "rawdata")

ggplot(me, aes(x = x, y = predicted, colour = group)) +
  geom_point(position = position_dodge(.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(.2)) +
  # using jitter for observations here
  geom_jitter(data = raw, mapping = aes(x = x, y = response, colour = group))

Or, also the simple way:

plot(me, rawdata = TRUE)

See also ?plot.ggeffects for the different plot options, for instance, use jitter = FALSE to remove the jittering for the raw data points etc...

Answer 2

The emmeans package is pretty useful:

pr <- emmeans::emmeans(fit, ~treatment * test_timing)
ggplot(summary(pr), 
       aes(test_timing, emmean, color = treatment, ymin = lower.CL, ymax = upper.CL)) +
  geom_pointrange(position = position_dodge(0.2), size = 1)