如何使用估计的边际均值生成可以自定义的交互图？

Question

我目前正在使用线性模型和方差分析来研究两个因素（温度、物种）对昆虫卵期持续时间的影响，如下：

## 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)

温度和物种之间的相互作用非常显着，因此我使用 emmeans 包中的 emmip() 函数绘制了一个简单的相互作用图，其中每个点都是估计的边际平均值，如下所示：

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

生成的图如下：

我想知道是否有人知道如何对此图进行以下编辑，或者知道如何生成相同图的替代方法，因为我知道 emmip() 函数不允许进行大量自定义。我想要做出的改变是：

1. 我想完全删除背景（所以它是纯白色的）
2. 我想将线条（以及图例中的）的点和颜色设为黑色，每条线都是不同的类型（例如实线和虚线），并且每条线的点都是不同的形状（例如 A 为圆形，A 为三角形） B).
3. 我希望置信区间为浅灰色
4. 在每个点周围，我希望有一个物种和温度的每种组合的所有数据点的抖动图。如果可能的话，由于重叠，每个物种的形状是否不同（或必要时不同的颜色）。

任何人可以提供的任何帮助将不胜感激。

Answer 1

你想要的是基础图形。

> warp.lm <- lm(breaks ~ wool*tension, data=warpbreaks)
> 
> library(emmeans)
> emm <- summary(emmeans(warp.lm, ~ wool*tension))
> 
> 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)),
+           col='grey')
+ )
> 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',
+              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), 
+        pch=as.integer(emm$wool), lty=as.integer(emm$wool) + 1)
> dev.off()

Answer 2

添加

plotit=FALSE

再次运行模型并保存结果。然后使用

emmip_ggplot

绘制模型。这允许您向绘图添加更多图层。我使用 emmip 包中的 auto.noise 数据集。

noise.lm <- lm(noise ~ size * type + side, data = auto.noise)

emmip.dat <- emmip(noise.lm, type ~ size, CIs = TRUE, plotit=FALSE)

emmip_ggplot(emmip.dat,  
       linearg = list(col="black", size=1), 
       CIarg = list(lwd = 1, alpha = 1, col="grey70"),
       dotarg = list(col="black", size=4)) +
  theme_classic() +
  theme(legend.key.width = unit(2,"cm")) +
  geom_jitter(data=mutate(auto.noise, tvar=type), aes(x=size, y=noise), col="black")