基于Bootstrap的R语言自定义函数实现GLM回归95%置信区间估算

通过参考建议，已实现R中基于Bootstrap方法为二项式GLM回归模型计算95%置信区间的代码，并封装为可复用的自定义函数，以下是相关实现：

原始实现代码 #

library(ggplot2)

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

fit <- glm(AE ~ Conc, family='binomial', data=df)
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, family='binomial', data=bdf)
  },
  simplify=FALSE
)

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

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

ggplot(data = ndf) +
  geom_line(aes(Conc, fit$family$linkinv(pred)))+
  geom_ribbon( aes(Conc, 
                 ymin = fit$family$linkinv(bpredci[, 1]),
                 ymax = fit$family$linkinv(bpredci[, 2])
                   ),
             fill = "grey", alpha = 0.5)

封装后的自定义函数代码 #

library(ggplot2)
library(dplyr) # 函数依赖dplyr包的select函数

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

mf <- function(mydata, x, y){
  # 拟合二项式GLM模型
  fit <- glm(as.formula(paste(y, "~", x)), family='binomial', data=mydata)
  # 构建预测数据集
  ndf <- mydata %>% select(all_of(x))
  # 基于连接函数获取预测值
  pred <- predict(fit, newdata=ndf, type='link') 

  # 执行Bootstrap重复抽样与模型拟合
  set.seed(42)
  bf <- replicate(
    999L, {
      bdf <- mydata[sample.int(nrow(mydata), replace=TRUE), ]
      glm(as.formula(paste(y, "~", x)), family='binomial', data=bdf)
    },
    simplify=FALSE
  )

  # 获取Bootstrap样本的预测值
  bpred <- sapply(bf, predict, newdata=ndf, type='link')

  # 计算95%置信区间
  ci <- \(mea, sd) mea + as.matrix(sd*(-qt(.025, Inf))) %*% cbind(-1, 1)
  bpredci <- ci(matrixStats::rowMeans2(bpred), matrixStats::rowSds(bpred))

  # 绘制回归曲线与置信区间
  ggplot(ndf) +
    geom_line(aes(.data[[x]], fit$family$linkinv(pred))) +
    geom_ribbon(aes(.data[[x]],
                  ymin = fit$family$linkinv(bpredci[, 1]),
                  ymax = fit$family$linkinv(bpredci[, 2])),
              fill = "grey",
              alpha = 0.5,
              inherit.aes = FALSE)
}

# 调用自定义函数示例
mf(mydata = df, "Conc", "AE")

内容的提问来源于stack exchange，提问作者Dylan Li