Take regression data from Motulsky and make a 2 panel plot in base R (plot) and ggplot (qplot)

Data set

From Motulsky’s Intuitive Biostatistics, 2nd ed, Chapter 32.

Make dataframe

#Add raw data to vectors
insulin.sensitivity <- c(250, 220, 145, 115, 230, 200, 330, 400,370,260,270,530,375)

C20.22.fattyacids <- c(17.9,18.3,18.3, 18.4, 18.4, 20.2, 20.3, 21.8,21.99,22.1,23.1,24.2,24.4)

#make dataframe
insulin.dat <- data.frame(insulin.sensitivity,
                          C20.22.fattyacids)

Regression

#mak model
insulin.mod <- lm(insulin.sensitivity ~ C20.22.fattyacids,
      data = insulin.dat)

#predictions from model
mod.preds <- predict(insulin.mod)

#add predictions to DF
insulin.dat$mod.preds <- mod.preds

#CIs from model
CIs <- predict(insulin.mod, interval = "confidence" )

#Add CIs to dataframe
insulin.dat <- cbind(insulin.dat,CIs)

Make Motulsky Chapter 35 Figure 35.1

2 panel plot in plot()

2-panel plot in ggplot

Load packages

library(cowplot)
## Loading required package: ggplot2
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
library(ggplot2)

See here for information on the “complot” package

Make plot

#Make the left plot
qplot1 <- qplot(y = insulin.sensitivity,
      x = C20.22.fattyacids,
      data = insulin.dat) +
  geom_hline(yintercept = mean(insulin.dat$insulin.sensitivity), 
             color = "blue", 
             size = 1.2) +
  ggtitle("Null hypothesis (Ho)") +
  xlab("%C20-22 Fatty Acids") +
  ylab( expression( "Insulin Sensitivity (mg/" ~~ m^{-2}~"/min)"))


#Make the right plot
qplot2 <- qplot(y = insulin.sensitivity,
      x = C20.22.fattyacids,
      ylab = "",
      data = insulin.dat) +
  geom_smooth(method = "lm") +
  ggtitle("Regression Model (Ha)")+
  xlab("%C20-22 Fatty Acids")


plot_grid(qplot1, 
          qplot2,
          #labels=c("A", "B"), 
          ncol = 2,
          nrow = 1)