如何引导 glm 回归，估计 95% 置信区间并绘制它？

Question

我正在使用 0-1 分布数据集进行

glm

回归。与

ggplot2::geom_smooth

相处得很顺利；这是我的代码：

library(ggplot2)
 
 set.seed(123)
 df = data.frame(Conc = runif(200, min = 200, max = 1000),    
                 AE = rbinom(200, 1, 0.5))
 
ggplot(df, aes(x = Conc, y = AE)) +   
  geom_jitter(height = 0.05, alpha = 0.5) +   
  geom_smooth(method = "glm", formula = y ~ log(x),
              method.args = list(family = "binomial"),
              color = "grey10")

现在，我想用 bootstrap 绘制 95% 置信区间。我尝试使用

boot

包和

ggplot2::mean_cl_boot

，但都失败了。

我没有保留所有不起作用的代码，但这是我尝试过的最新代码。说实话，我从其他答案中复制并尝试了这些代码，但我并不完全理解这些代码。

lm_coeffs = function(x, y) {
  coeffs = coefficients(lm(y~log(x)))
  tibble(C = coeffs[1], m=coeffs[2])
}

nboot = 1000

mtboot = lapply (seq_len(nboot), function(i) 
  df %>%
    slice_sample(prop=1, replace=TRUE) %>% 
    summarise(tibble(lm_coeffs(Conc, AE))))
mtboot = do.call(rbind, mtboot)

ggplot(df, aes(Conc, AE)) +
  geom_abline(aes(intercept=C, slope=m), data = mtboot,
              size=0.3, alpha=0.3, color='forestgreen') +
  geom_point() 

Answer 1

这是你的模型，

> fit <- glm(AE ~ Conc, df, family='binomial')
> ndf <- data.frame(Conc=seq(min(df$Conc), max(df$Conc), length.out=1e3))
> pred <- predict(fit, newdata=ndf, type='link')

您想要引导，

> set.seed(42)
> bf <- replicate(
+   999L, {
+     bdf <- df[sample.int(nrow(df), replace=TRUE), ]
+     glm(AE ~ Conc, bdf, family='binomial')
+   },
+   simplify=FALSE
+ )

并且从引导的

pred

ictions，

> bpred <- sapply(bf, predict, newdata=ndf, type='link')

您想要 95% CI。

> ci <- \(x, sd) x + as.matrix(sd*(-qt(.025, Inf))) %*% cbind(-1, 1)
> bpredci <- ci(matrixStats::rowMeans2(bpred), matrixStats::rowSds(bpred))

要绘制它，您需要应用

fit$family$linkinv()

函数，正如 @GavinSimpson 不久前在此 answer 中所解释的那样，如下所示：

> plot(df)
> lines(fit$family$linkinv(pred), col=4)
> for (j in 1:2) lines(fit$family$linkinv(bpredci[, j]), col=4, lty=2)

Gavin的答案提供了ggplot2方式。

PS：如果我在 OP 中使用样本数据，它看起来与那里的第一个图类似。

---

数据：

> set.seed(42)
> df <- transform(
+   data.frame(Conc=runif(200, min=200, max=1000)),
+   AE=rbinom(200, 1, prob=plogis((Conc - 600)/100))
+ )