我正在使用0-1分布数据集进行glm次回归.它与ggplot2::geom_smooth进行得很好;这是我的代码:

library(ggplot2)
 
 set.seed(123)
 df <- transform(data.frame(Conc=runif(200, min=200, max=1000)),
 AE=rbinom(200, 1, prob=plogis((Conc - 600)/100)))
 
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")

enter image description here

现在,我想用Bootstrap绘制95%置信区间.我try 了boot包和ggplot2::mean_cl_boot包,但都失败了.

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

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() 

推荐答案

这是你的模型,

> 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')

你想要 bootstrap 的东西,

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

从 bootstrap 的pred个版本中,

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

您想要95%CI(显示百分位 bootstrap 方法).

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

要绘制它,您需要应用fit$family$linkinv()函数,正如不久前这answer中解释的@GavinSimpson,如下所示:

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

enter image description here

Gavin的答案提供了ggplot 2的方式.

PS:如果我在OP中使用样本数据,它看起来与那里的第一个图相似.


Data:

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

R相关问答推荐

有没有方法将琴弦完全捕捉到R中的多边形?

R形式的一维数字线/箱形图样式图表

矩阵%*%矩阵中的错误:需要数字/复杂矩阵/向量参数

根据选中三个复选框中的一个或两个来调整绘图

如何写一个R函数来旋转最后n分钟?

如何基于两个条件从一列中提取行

函数可以跨多个列搜索多个字符串并创建二进制输出变量

R -使用矩阵reshape 列表

通过初始的shiny 应用更新部署的shiny 应用的数据和参数,其中部署的应用程序显示为URL

将统计检验添加到GGPUBR中的盒图,在R

SHILINY中DT列的条件着色

自定义交互作用图的标签

按两个因素将观测值分组后计算单独的百分比

填充图例什么时候会有点?

如果极点中存在部分匹配,则替换整个字符串

对一个数据帧中另一个数据帧中的值进行计数

如何在一种 colored颜色 中设置数值变量的 colored颜色 和高于阈值的 colored颜色 点?

将字符变量出现次数不相等的字符框整形为pivot_wider,而不删除重复名称或嵌套字符变量

使用相对风险回归计算RR

Ggplot2水平线和垂直线的图例图标不匹配