Plot percentage change figure with 95% CI and stats

Question

I am planning to reproduce the attached figure, but I have no clue how to do so:

Let´s say I would be using the CO2 example dataset, and I would like to plot the relative change of the Uptake according to the Treatment. Instead of having the three variables in the example figure, I would like to show the different Plants grouped for each day/Type.

So far, I managed only to get this bit of code, but this is far away from what it should look like.

aov1 <- aov(CO2$uptake~CO2$Type+CO2$Treatment+CO2$Plant)
plot(TukeyHSD(aov1, conf.level=.95))

Axes should be switched, and I would like to add statistical significant changes indicated with letters or stars.

Answer 1

You can do this by building it in base R - this should get you started. See comments in code for each step, and I suggest running it line by line to see what's being done to customize for your specifications:

Set up data

# Run model
aov1 <- aov(CO2$uptake ~ CO2$Type + CO2$Treatment + CO2$Plant)

# Organize plot data
aov_plotdata <- data.frame(coef(aov1), confint(aov1))[-1,] # remove intercept
aov_plotdata$coef_label <- LETTERS[1:nrow(aov_plotdata)] # Example labels

Build plot

#set up plot elements
xvals <- 1:nrow(aov_plotdata)
yvals <- range(aov_plotdata[,2:3])

# Build plot
plot(x = range(xvals), y = yvals, type = 'n', axes = FALSE, xlab = '', ylab = '') # set up blank plot
  points(x = xvals, y = aov_plotdata[,1], pch = 19, col = xvals) # add in point estimate
  segments(x0 = xvals, y0 = aov_plotdata[,2], y1 =  aov_plotdata[,3], lty = 1, col = xvals) # add in 95% CI lines
  axis(1, at = xvals, label = aov_plotdata$coef_label) # add in x axis
  axis(2, at = seq(floor(min(yvals)), ceiling(max(yvals)), 10)) # add in y axis
  segments(x0=min(xvals), x1 = max(xvals), y0=0, lty = 2) #add in midline
  legend(x = max(xvals)-2, y = max(yvals), aov_plotdata$coef_label, bty = "n", # add in legend
         pch = 19,col = xvals, ncol = 2)