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

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, col="black",  
                    linearg = list(size=1), 
                    CIarg = list(lwd = 1, alpha = 1, col="grey70"),
                    dotarg = list(size=4)) +
  theme_classic() +
  theme(legend.key.width = unit(2,"cm"))