如何在R中绘制逻辑回归glm预测值和置信区间

Question

我有一个二项式 glm，其中包含一个存在/不存在的响应变量和一个具有9个级别的因子变量，如下所示：

data$y<-factor(data$y,levels=c(0,1),labels=c("absent","present"))
table(data$y,data$site_name)

          Andulay Antulang Basak Dauin Poblacion District 1 Guinsuan Kookoo's Nest Lutoban Pier Lutoban South Malatapay Pier
  absent        4        4     1                          0        3             1            5             5              2
  present       2        2     5                          6        3             5            1             1              4

model <- glm(y~site_name,data=data,binomial)

仅为简洁起见，跳过模型推断和验证，我该如何在一个箱线图中绘制每个站点出现概率及其置信区间？我想要的是类似于在R中绘制预测概率和置信区间所示的内容，但我希望用箱线图展示它，因为我的回归变量site_name是一个具有9个级别的因子，而不是一个连续变量。

我认为可以按以下方式计算必要的值（但不确定正确性）：

将模型系数转换为成功概率的函数：

calc_val <- function(x){return(round(1/(1+1/(exp(x))),3))}

基于模型的预测概率：

prob <- tapply(predict(model,type="response"),data$site_name,function(x){round(mean(x),3)})
means <- as.data.frame(prob)

预测概率的75%和95%置信区间：

ci <- cbind(confint(model,level=0.9),confint(model,level=0.5))
rownames(ci) <- gsub("site_name","",rownames(ci))
ci <- t(apply(ci,1,calc_val))

将所有内容整合到一张表格中

ci<-cbind(means,ci)
ci
                            prob   5 %  95 %  25 %  75 %   Pr(>|z|) stderr
Andulay                    0.333 0.091 0.663 0.214 0.469 0.42349216  0.192
Antulang                   0.333 0.112 0.888 0.304 0.696 1.00000000  0.192
Basak                      0.833 0.548 0.993 0.802 0.964 0.09916496  0.152
Dauin Poblacion District 1 1.000 0.000    NA 0.000 1.000 0.99097988  0.000
Guinsuan                   0.500 0.223 0.940 0.474 0.819 0.56032414  0.204
Kookoo's Nest              0.833 0.548 0.993 0.802 0.964 0.09916496  0.152
Lutoban Pier               0.167 0.028 0.788 0.130 0.501 0.51171512  0.152
Lutoban South              0.167 0.028 0.788 0.130 0.501 0.51171512  0.152
Malatapay Pier             0.667 0.364 0.972 0.640 0.903 0.25767454  0.192

所以我的问题有两个：

1. 概率和置信区间的计算是否正确？
2. 如何在箱线图中绘制这个数据？

编辑 这是一些通过 dput 提供的样本数据（也修改了上面的表格以匹配数据）：

# dput(data[c("y", "site_name")])
data <- structure(list(y = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L), .Label = c("absent", "present"), class = "factor"), site_name = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 9L, 9L, 9L, 9L, 9L, 9L, 4L, 4L, 4L, 4L, 4L, 4L, 8L, 8L, 8L, 8L, 8L, 8L, 7L, 7L, 7L, 7L, 7L, 7L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("Andulay", "Antulang", "Basak", "Dauin Poblacion District 1", "Guinsuan", "Kookoo's Nest", "Lutoban Pier", "Lutoban South", "Malatapay Pier"), class = "factor")), .Names = c("y", "site_name"), row.names = c(125L, 123L, 126L, 124L, 128L, 127L, 154L, 159L, 157L, 158L, 156L, 155L, 111L, 114L, 116L, 115L, 112L, 113L, 152L, 151L, 148L, 150L, 153L, 149L, 143L, 146L, 144L, 147L, 142L, 145L, 164L, 165L, 161L, 163L, 160L, 162L, 120L, 122L, 121L, 117L, 118L, 119L, 137L, 136L, 139L, 141L, 140L, 138L, 129L, 134L, 131L, 135L, 133L, 130L), class = "data.frame")
#

- Dolf Andringa

---

有没有可能提供一个可重现的例子呢？ - Ben Bolker

我添加了一些数据。我不知道我展示上面的表格是哪个物种，所以我更新了它们以匹配粘贴的数据。但问题更多地涉及技术性质。如何获得逻辑回归预测值和置信区间的箱线图？当然，除非我完全理解错了什么（这并非不可想象）。 - Dolf Andringa

你好，请提供你的数据，你能编辑你的问题，加入执行 dput(data[c("y", "site_name")]) 的结果吗？（希望其他人可以从你的问题中复制数据到他们的 R 会话中 - 在你发布的格式下无法做到这一点） - user20650

Answer 1

这是一个最基础、只包含基础功能的解决方案。

拟合模型：

mm <- glm(y~site_name,data=dd,family=binomial)

使用网站名称创建一个预测框架：

pframe <- data.frame(site_name=unique(dd$site_name))

使用标准误差预测（在logit/线性预测器比例尺上）

pp <- predict(mm,newdata=pframe,se.fit=TRUE)
linkinv <- family(mm)$linkinv ## inverse-link function

将预测结果、下限和上限结合起来，再进行反变换到概率尺度：

pframe$pred0 <- pp$fit
pframe$pred <- linkinv(pp$fit)
alpha <- 0.95
sc <- abs(qnorm((1-alpha)/2))  ## Normal approx. to likelihood
alpha2 <- 0.5
sc2 <- abs(qnorm((1-alpha2)/2))  ## Normal approx. to likelihood
pframe <- transform(pframe,
                    lwr=linkinv(pred0-sc*pp$se.fit),
                    upr=linkinv(pred0+sc*pp$se.fit),
                    lwr2=linkinv(pred0-sc2*pp$se.fit),
                    upr2=linkinv(pred0+sc2*pp$se.fit))

情节。

with(pframe,
{
    plot(site_name,pred,ylim=c(0,1))
    arrows(as.numeric(site_name),lwr,as.numeric(site_name),upr,
           angle=90,code=3,length=0.1)
})

作为一个箱线图：

with(pframe,
{
    bxp(list(stats=rbind(lwr,lwr2,pred,upr2,upr),
             n = rep(1,nrow(pframe)),
             conf = NA,
             out = NULL,
             group = NULL,
             names=as.character(site_name)))
})

有很多其他的方法可以做到这一点；我建议

library("ggplot2")
ggplot(pframe,aes(site_name,pred))+
     geom_pointrange(aes(ymin=lwr,ymax=upr))+
     geom_linerange(aes(ymin=lwr2,ymax=upr2),lwd=1.5)+
     coord_flip()

另一种解决方案是通过 y~site_name-1 拟合模型，在这种情况下，将为每个站点的概率分配单独的参数，并使用 profile()/confint() 找到置信区间；这比在上面的答案中依赖于参数/预测的抽样分布的正态性要稍微更准确一些。

- Ben Bolker