Visualizing adjusted means (LS means or Estimated Marginal means) procedure in ANCOVA (ggplot2)

Question

Trying to figure out how to visually demonstrate how adjusted means works in ANCOVA. There are some good published examples in the primary literature, but I haven't been able to replicate their visualizations with ggplot2. Examples I am trying to replicate:

Packard and Boardman 1999 (Figure 2) and Barrett 2011 (Figure 1)

library(ggplot2)
library(grid)
library(emmeans)
library(HH)
library(multcomp)

Example using 'litter' data:

data(litter)

gest.mean <- mean(litter$gesttime) #mean of the covariate

model1 <- lm(weight ~ gesttime * dose, data=litter)
pred1 <- predict(model1)

model2 <- lm(weight ~ gesttime + dose, data=litter)
pred2 <- predict(model2)

#plot different slopes
plot1 <- ggplot(data = cbind(litter, pred1),
   aes(gesttime, weight, color=dose)) + geom_point() +
geom_line(aes(y=pred1))+  #plots the predicted values (fitted line)
geom_vline(xintercept = gest.mean, linetype="dashed")+
labs(title = "Model1: Separate Slopes ANCOVA", subtitle = "model1 <- 
lm(weight ~ gesttime * dose, data=litter)")

#plot same slopes
plot2 <- ggplot(data = cbind(litter, pred2),
   aes(gesttime, weight, color=dose)) + geom_point() +
geom_line(aes(y=pred2))+
geom_vline(xintercept = gest.mean, linetype="dashed")+
labs(title = "Model2: Equal Slopes ANCOVA", subtitle = "model2 <- lm(weight ~ 
gesttime + dose, data=litter)")

#dashed vertical line shows the mean of covariate
#emmeans are calculated by adjusting points to mean of covariate along group specific slope

grid.newpage()
grid.draw(rbind(ggplotGrob(plot1), ggplotGrob(plot2), size = "last"))

summary(model1)
aov(model1)
summary(model2)
aov(model2)

#compare fits of model with interaction (sep. slopes) vs. model without (eq. slopes)
anova(model1,model2)

#EMmean post hocs to compare differences among four treatments at the grand mean of the covariate
#same as comparing intercepts when slopes are equal

#calculate model1 estimated marginal means (using interaction)
model1.emm <- emmeans(model1, "dose") #note that is gives warning message because sep slopes (interaction)
pairs(model1.emm)

#compare model1 marginal means (LS means)
plot(model1.emm, comparisons = TRUE)
CLD(model1.emm)

#calculate model2 estimated marginal means
model2.emm <- emmeans(model2, "dose")
pairs(model2.emm)

#compare model2 marginal means (LS means)
plot(model2.emm, comparisons = TRUE)
CLD(model2.emm)

#Just to show how EM means are used (intersect grand mean of covariate) 
plot3 <- ggplot(data = cbind(litter, pred2),
            aes(gesttime, weight, color=dose)) + geom_point() +
geom_line(aes(y=pred2))+
geom_vline(xintercept = gest.mean)+
geom_hline(yintercept = 28.87, linetype="dashed", color=c(1))+
geom_hline(yintercept = 29.33, linetype="dashed")+
geom_hline(yintercept = 30.56, linetype="dashed")+
geom_hline(yintercept = 32.35, linetype="dashed")+
labs(title = "Model2: Equal Slopes ANCOVA")

plot3

plot4 <-plot3 +
geom_segment(mapping=aes(x=gesttime, xend=gesttime+0.5, y=weight, 
yend=weight+0.5, colour = "dose"), arrow=arrow(), size=.25, color="blue")

plot4 
#obviously not what I wanted; individuals are not connected to mean of covariate (gesttime=22.08) along group-specific slope (sep. slopes) or common slope (eq. slopes)

Answer 1

You could do

library(emmeans)
plt = emmip(model2, dose ~ gesttime,
    cov.reduce = range)

So far, you have a ggplot object with the fitted lines. Now, get the adjusted means.

emmdat = as.data.frame(emmeans(model2, ~ dose*gesttime))

That data frame has the predictor values and EMMs that you need to plot. Add the appropriate ggplot() code to plot these points to plt, and display the result.

The same, using model1, will illustrate why you get a warning! The adjusted means compare differently at different gestation times.

Answer 2

Here is code to replicate the Barrett 2011 ANCOVA plot (Figure1). I follow the procedure of fitting an interaction first (separate slopes) and removing non-significant interaction to yield a minimum adequate model using equal slopes to fit adjusted values and adjusted means (LS means or EM means).

library(ggplot2)
library(dplyr)
library(grid)

#extract data from the Barrett 2011 paper
X <- c(11,21,30,41,52,65,71,77,8,17,29,42,51,64,72,79)
Y <- c(33,32,38,49,51,53,59,65,20,22,31,28,42,52,48,55)
Group <- as.factor(c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2))
data <-data.frame(X,Y,Group)

X.mean <- mean(data$X) #mean of the covariate

model1 <- lm(Y ~ X * Group, data=data)
pred1 <- predict(model1)

model2 <- lm(Y ~ X + Group, data=data)
pred2 <- predict(model2)

#plot different slopes
plot1 <- ggplot(data = cbind(data, pred1),
                aes(X, Y, color=Group)) + geom_point() +
  geom_line(aes(y=pred1))+  #plots the predicted values (fitted line)
  geom_vline(xintercept = X.mean, linetype="dashed", alpha = 0.15)+
  labs(title = "Model1: Separate Slopes ANCOVA", subtitle = "model1 <- lm(Y ~ X * Group, data=data)")+
  theme_classic()

#plot same slopes
plot2 <- ggplot(data = cbind(data, pred2),
                aes(X, Y, color=Group)) + geom_point() +
  geom_line(aes(y=pred2))+
  geom_vline(xintercept = X.mean, linetype="dashed", alpha = 0.15)+
  labs(title = "Model2: Equal Slopes ANCOVA", subtitle = "model2 <- lm(Y ~ X + Group, data=data)")+
  theme_classic()

grid.newpage()
grid.draw(rbind(ggplotGrob(plot1), ggplotGrob(plot2), size = "last"))

summary(model1)
anova(model1)
summary(model2)
anova(model2)
anova(model1, model2) #no sig. difference, drop interaction term and use simplest model (equal slopes)

plot3 <- ggplot(data = cbind(data, pred2),
                aes(X, Y, color=Group)) + 
  geom_point()+
  geom_line(aes(y=pred2))+
  #geom_vline(xintercept = X.mean, linetype="dashed", alpha = 0.45)+
  labs(title = "Model2: Equal Slopes ANCOVA", subtitle = "model2 <- lm(Y ~ X + Group, data=data)")+
  theme_classic()
plot3

#mutate to calc adjusted values of individuals
data <- data%>%mutate(adjY=Y-0.498*(X-X.mean)) 
#0.498 is the 'common slope' of model2; equal slopes ANCOVA 

plot4 <- ggplot(data = cbind(data, pred2),
                aes(X, Y, color=Group)) + 
  geom_point()+
  #geom_line(aes(y=pred2))+
  geom_vline(xintercept = X.mean, linetype="dashed", alpha = 0.45)+
  #labs(title = "Model2: Equal Slopes ANCOVA", subtitle = "model2 <- lm(Y ~ X + Group, data=data)")+
  geom_segment(aes(x=X, xend=X.mean, y=Y, yend=data$adjY), size=.25)+
  theme_classic()
plot4 

plot5 <-ggplot(data, aes(x=Group, y=adjY, color=Group))+
  geom_point()+
  stat_summary(geom="point", fun.y= "mean", shape = 8, color="black", size=5)+
  geom_hline(yintercept = mean(data$adjY), linetype="dashed", alpha = 0.45)+
  theme_classic()
plot5

grid.newpage()
grid.draw(rbind(ggplotGrob(plot3), ggplotGrob(plot4), ggplotGrob(plot5), size = "last"))

Barrett 2011 Figure 1