r - 如何为一个样本比例绘制 95% 置信区间图

Question

qplot(x="SignUp", y=0.07, ymin=Lower_Level, ymax=Upper_Level, ylim=c(0,1), geom = "pointrange")+coord_flip() +
  
  ylab("SignUp Proportion")+geom_hline(yintercept=Upper_Level)+geom_hline(yintercept=Lower_Level)

这就是我设法绘制的。但我想要类似下图的东西。置信区间为 0.084 和 0.0551。样本比例为0.07

标签： rggplot2sampleconfidence-interval

Answer 1

我想你可以像这样显示估计概率的 95% 置信区间：

首先，从 1 和 0 的数据帧开始，代表您在样本中的“成功”和“失败”率。在这里，您的数字表明 1500 次成功中大约有 105 次，所以我们这样做：

df <- data.frame(x = c(rep(1, 105), rep(0, 1395)))

现在我们拟合一个逻辑回归，截距是我们估计的唯一参数：

mod <- coef(summary(glm(x ~ 1, family = binomial, data = df)))

mod
#>              Estimate Std. Error  z value      Pr(>|z|)
#> (Intercept) -2.586689  0.1011959 -25.5612 4.122466e-144

这里的估计应该与给定的估计和标准误差呈正态分布（在对数赔率尺度上），因此我们可以通过执行以下操作来获取适当范围内的密度值：

xvals <- seq(mod[1] - 3 * mod[2], mod[1] + 3 * mod[2], 0.01)
yvals <- dnorm(xvals, mod[1], mod[2])

现在我们将 x 值从对数赔率转换为概率：

pxvals <- exp(xvals)/(1 + exp(xvals))

我们还将创建一个向量来标记这些值是否在估计值的 1.96 标准差内：

level <- ifelse(xvals < mod[1] - 1.96 * mod[2], "lower",
          ifelse(xvals > mod[1] + 1.96 * mod[2], "upper", "estimate"))

现在我们将所有这些都放在一个数据框中并绘制：

plot_df <- data.frame(xvals, yvals, pxvals, level)

library(ggplot2)

ggplot(plot_df, aes(pxvals, yvals, fill = level)) +
  geom_area(alpha = 0.5) +
  geom_vline(xintercept = exp(mod[1])/(1 + exp(mod[1])), linetype = 2) + 
  scale_fill_manual(values = c("gray70", "deepskyblue4", "deepskyblue4"),
                    guide = guide_none()) +
  scale_x_continuous(limits = c(0.03, 0.13), breaks = 3:12/100,
                     name = "probability") +
  theme_bw()

请注意，因为我们已经转换了 x 轴，所以这不再是真正的密度图。结果，y 轴变得有些随意，但该图仍然准确地显示了概率估计的 95% 置信区间。

---

编辑

glm如果方法看起来太复杂，这是一种替代方法。它使用二项分布来获得 95% 的置信区间。您只需提供人口规模和“成功”数量

library(ggplot2)

population <- 1500
actual_successes <- 105
test_successes <- 1:300

density <- dbinom(test_successes, population, actual_successes/population)
probs   <- pbinom(test_successes, population, actual_successes/population)
label   <- ifelse(probs < 0.025, "low", ifelse(probs > 0.975, "high", "CI"))

ggplot(data.frame(probability = test_successes/population, density, label),
       aes(probability, density, fill = label)) +
  geom_area(alpha = 0.5) +
  geom_vline(xintercept = actual_successes/population, linetype = 2) + 
  scale_fill_manual(values = c("gray70", "deepskyblue4", "deepskyblue4"),
                    guide = guide_none()) +
  scale_x_continuous(limits = c(0.03, 0.13), breaks = 3:12/100,
                     name = "probability") +
  theme_bw()