Rbase 中同一图中两个模型的阴影置信区间

Question

我有一个数据集，其中包含 142 个关于在四个站点（Site）处通过不同宽度间隙（LargDiscon）的个体交叉概率（nRF01）的观察结果，如下所示：

  LargDiscon nRF01            Site
1         90     0 Azay-sur-Thouet
2         79     0 Azay-sur-Thouet
3         55     0 Azay-sur-Thouet
4        125     0 Azay-sur-Thouet
5         40     0 Azay-sur-Thouet
6         28     0 Azay-sur-Thouet

我使用二项式误差分布进行了 GLM，以宽度间隙、位置（图中的彩色斜率）和两者之间的相互作用作为解释变量。我还需要整体模型（图中的黑色斜率），因为我们使用该模型进行其他分析。

gt <- glm(nRF01 ~ LargDiscon*Site, data=p, family="binomial") #model with four sites
gt2 <- glm(nRF01 ~ LargDiscon, data=p, family="binomial") #model overall

之后，我用站点绘制模型，但我还使用以下代码添加了整体坡度：

# Plot
par(cex=0.8, cex.lab = 1.5, cex.sub=1) # increase size font
par(cex.lab=1.7, cex.axis=1.9, cex.main=1.9,
    mar = c(5, 5, 5, 5))
plot(jitter(nRF01, 0.1) ~ LargDiscon, 
     data=p, t="n", 
     xlim=c(0, 170), 
     xaxs="i", 
     ylim=c(-0.05, 1.05), 
     xlab="Gap width (m)", 
     ylab="Probability of crossing", 
     main=" ")


lg <- seq(1, 170, 1)
nd <- data.frame(LargDiscon=rep(lg,times=4), Site=rep(levels(p$Site), each=length(lg))) #predictions
nd$pred <- as.vector(predict(gt, nd, type="response"))  # using model to predict values of points 
pred <- as.vector(predict(gt2,  list(LargDiscon=lg), type="response")) #pred of overall model

# add confidence interval
# how ?

abline(h=c(0,1), col="darkgrey")
abline(v=0, col="darkgrey")
col <- c("#FF6101", "#D90000", "#0093A5", "#902E98")
lty <- c("dotted", "dotdash", "twodash", "longdash")

for(i in 1:4) { 
    points(pred ~ LargDiscon, data=nd[nd$Site==levels(p$Site)[i], ], 
           lwd=2, 
           col=col[i],
           lty = lty[i],
           t="l")
    }
points(lg, pred, lwd=2, col="black", t="l")
points(jitter(nRF01, 0.08) ~ LargDiscon, data=p, pch=20, bg="black", col = c("#FF6101", "#D90000", "#0093A5", "#902E98"))
abline(h=0.5, col="darkgrey", lty=2)

但是我确实需要为每个模型添加置信区间，但我失败了。我尝试使用 ggplot 但也失败了。有人可以帮我在这张图中添加置信区间吗？

Answer 1

我用 ggplot 成功了。我已经包含了下面的脚本。这允许您在同一个图表上显示您的模型 + 整体模型。

# Create dataframe with fittedvalues
nd <- data.frame(LargDiscon=rep(lg,times=4), Site=rep(levels(p$Site), each=length(lg))) # predictions
## add the fitted values by predicting from the model for the new data
nd <- add_column(nd, fit = predict(gt, newdata = nd, type = 'response'))
nd <- add_column(nd, wrong_se = predict(gt, newdata = nd, type = 'response',
                                        se.fit = TRUE)$se.fit)
nd <- mutate(nd, wrong_upr = fit + (2 * wrong_se), wrong_lwr = fit - (2 * wrong_se))


# Same for the overall model
lg <- seq(1, 170, 1) #width of my explcative variable
nd2 <- data.frame(LargDiscon = lg, fit = as.vector(predict(gt2,  list(LargDiscon=lg), type="response")))
nd2 <- add_column(nd2, fit = predict(gt2, newdata = nd2, type = 'response'))
nd2 <- add_column(nd2, wrong_se = predict(gt2, newdata = nd2, type = 'response',
                                        se.fit = TRUE)$se.fit)
nd2 <- mutate(nd2, wrong_upr = fit + (2 * wrong_se), wrong_lwr = fit - (2 * wrong_se))
nd2$Site <- "Overall"


## plot it

ggplot() +
  # hline
  geom_hline(yintercept=1, linetype="dashed", color = "darkgrey")+
  geom_hline(yintercept=0, linetype="dashed", color = "darkgrey")+
  geom_hline(yintercept=0.5, linetype="dashed", color = "darkgrey")+
  #raw data
  geom_jitter(data = p, aes(x = LargDiscon, y = nRF01, fill = Site, color = Site), width = 0.005, height = 0.01, alpha = 0.5) +
  # model site
  geom_line(data = nd, aes(x = LargDiscon, y = fit, colour = Site), size = 1) +
  geom_ribbon(data = nd, aes(x = LargDiscon, y = fit, ymin = wrong_lwr, ymax = wrong_upr, fill = Site),
              alpha = 0.1,
              colour = NA) +
  # #overall model
  geom_line(data=nd2, aes(x=LargDiscon, y=fit, color = Site), size=1) +
  geom_ribbon(data = nd2, aes(x = LargDiscon, y = fit, ymin = wrong_lwr, ymax = wrong_upr),
              alpha = 0.1,
              colour = NA) +
  # geom_ribbon(data = nd2, aes(ymin = wrong_lwr, ymax = wrong_upr),
  #             alpha = 0.1,
  #             colour = NA) +
  # geom_line(nd2, aes(x = LargDiscon, y = fit), size = 1) +
  #scale_colour_discrete(name = 'Site') +
  scale_fill_manual(values = c("#FF6101", "#D90000", "#0093A5", "#902E98", "black"), aesthetics = c("color", "fill")) +
  labs(x = 'Gap width (m)', y = 'Probability of gap crossing') +
  theme_classic()