rplotinteractionemmeansjitter

How do I produce a customizable interaction plot using the estimated marginal means?

I am currently looking at how the effect of two factors (temperature, species) on the duration of the egg stage in insects using a linear model and ANOVA analysis as follows:

## Model
type3.egg <- list(Temperature=contr.sum, Species=contr.sum)
model <- lm(Duration.egg ~ Temperature*Species, data=egg.2, contrasts=type3.egg)
summary(model)
library(car)
Anova(model, type=3)

The interaction between temperature and species was significant so I plotted a simple interaction plot using the emmip() function in the package emmeans where each point is the estimated marginal mean as follows:

library(emmeans)
emmip(model, Species ~ Temperature, type="response", xlab="Temperature (°C)", 
      ylab="Duration of the egg stage (Days)", CI=TRUE)

The plot that was produced is as follows: enter image description here

I was wondering if anyone knew how I could make the following edits to this plot or knew of an alternative method to producing this same plot as I am aware that the emmip() function does not allow a great deal of customization. The changes I would like to make are:

  1. I would like to remove the background completely (so it is plain white)
  2. I would like to make the points and colors of the lines (and in the legend) black with each line being a different type (e.g. solid and dashed) and the points for each line being different shapes (e.g. circles for A, triangles for B).
  3. I would like the confidence intervals to be a light grey color
  4. Around each point I would like there to be a jitter plot of all of the data points for each combination of species and temperature. If possible, due to overlap, could these be different shapes for each species (or different colors if necessary).

Any help anyone can provide would be greatly appreciated.

Solution

  • What you want is base graphics. Use emmeans::emmeans to compute EMM beforehand.

    Example (using warpbreaks data):

    > warp.lm <- lm(breaks ~ wool*tension, data=warpbreaks)
> 
> (emm <- summary(emmeans::emmeans(warp.lm, ~ wool*tension), infer=c(TRUE, TRUE)))
 wool tension emmean   SE df lower.CL upper.CL t.ratio p.value
 A    L         44.6 3.65 48     37.2     51.9  12.218  <.0001
 B    L         28.2 3.65 48     20.9     35.6   7.739  <.0001
 A    M         24.0 3.65 48     16.7     31.3   6.581  <.0001
 B    M         28.8 3.65 48     21.4     36.1   7.891  <.0001
 A    H         24.6 3.65 48     17.2     31.9   6.734  <.0001
 B    H         18.8 3.65 48     11.4     26.1   5.149  <.0001
> 
> png('foo.png', 600, 400)
> par(mar=c(4.5, 4.5, 2, 2))
> with(warpbreaks, 
+      plot(x=jitter(as.integer(as.factor(tension))), y=breaks, 
+           xlab='Temperature (°C)', ylab='Duration of the egg stage (Days)',
+           pch=as.integer(as.factor(wool)), xaxt='n',
+           col='grey')
+      )
> with(warpbreaks, axis(1, at=seq_along(levels(tension)), labels=levels(tension)))
> for (i in seq_along(levels(emm$wool))) {
+   with(subset(emm, wool == levels(emm$wool)[i]),
+        lines(x=as.integer(tension) + switch(i, -.05, .05), 
+              y=emmean, type='o', cex=1.5,
+              lty=as.integer(wool), pch=as.integer(wool))
+   )
+ }
> arrows(as.integer(as.factor(emm$tension)) + c(-.05, .05), emm$lower.CL,
+        as.integer(as.factor(emm$tension)) + c(-.05, .05), emm$upper.CL, 
+        code=3, angle=90, length=.05, col='black', lwd=2, lty=1:2)
> legend('topright', title='Species', levels(emm$wool), cex=1.5,
+                pch=as.integer(emm$wool), lty=as.integer(emm$wool))
> dev.off()

    enter image description here