使用 dplyr 计算 95%-CI 的长度
Posted
技术标签:
【中文标题】使用 dplyr 计算 95%-CI 的长度【英文标题】：Calculating length of 95%-CI using dplyr 【发布时间】：2016-06-27 11:25:17 【问题描述】：上次我问如何计算一个变量 (procras) 的每个测量场合 (周) 的平均分数，该变量 (procras) 已为多个受访者重复测量。所以我的（简化的）长格式数据集如下所示（这里有两个学生，5 个时间点，没有分组变量）：
studentID week procras
1 0 1.4
1 6 1.2
1 16 1.6
1 28 NA
1 40 3.8
2 0 1.4
2 6 1.8
2 16 2.0
2 28 2.5
2 40 2.8
使用 dplyr 我会得到每个测量场合的平均分数
mean_data <- group_by(DataRlong, week)%>% summarise(procras = mean(procras, na.rm = TRUE))
看起来像这样，例如：
Source: local data frame [5 x 2]
occ procras
(dbl) (dbl)
1 0 1.993141
2 6 2.124020
3 16 2.251548
4 28 2.469658
5 40 2.617903
使用 ggplot2，我现在可以绘制随时间的平均变化，并且通过轻松调整 dplyr 的 group_data()，我还可以获得每个子组的平均值（例如，男性和女性每个场合的平均得分）。 现在，我想在 mean_data 表中添加一列，其中包括 95%-CI 的长度，围绕每个场合的平均得分。
http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/ 解释了如何获取和绘制 CI，但是一旦我想对任何子组执行此操作，这种方法似乎就会出现问题，对吧？那么有没有办法让 dplyr 也自动在 mean_data 中包含 CI（基于组大小等）？ 之后，将新值作为 CI 绘制到我希望的图表中应该相当容易。 谢谢。
【问题讨论】：
【参考方案1】：您可以使用mutate 中的一些额外功能手动完成summarise
library(dplyr)
mtcars %>%
group_by(vs) %>%
summarise(mean.mpg = mean(mpg, na.rm = TRUE),
sd.mpg = sd(mpg, na.rm = TRUE),
n.mpg = n()) %>%
mutate(se.mpg = sd.mpg / sqrt(n.mpg),
lower.ci.mpg = mean.mpg - qt(1 - (0.05 / 2), n.mpg - 1) * se.mpg,
upper.ci.mpg = mean.mpg + qt(1 - (0.05 / 2), n.mpg - 1) * se.mpg)
#> Source: local data frame [2 x 7]
#>
#> vs mean.mpg sd.mpg n.mpg se.mpg lower.ci.mpg upper.ci.mpg
#> (dbl) (dbl) (dbl) (int) (dbl) (dbl) (dbl)
#> 1 0 16.61667 3.860699 18 0.9099756 14.69679 18.53655
#> 2 1 24.55714 5.378978 14 1.4375924 21.45141 27.66287
【讨论】：
谢谢，这对我来说几乎是完美的，我也可以用 ggplot 绘制 CI。我唯一的问题是，无论他们是否失踪，n.mpg = n()) 总是给我相同的数字，即参与者总数（n = 566）。由于纵向设计，出现了 dropout，所以使用总 n 会使 CI 不准确，因为 SE 和 df 是错误的。我试图通过从 n() 参数中减去 'sum(as.numeric(is.na(DataRlong$procras)))' 来解决这个问题，但这会减去所有情况下丢失案例的总数。 我怎么能告诉 r 只从 n 中减去在相应测量场合丢失的案例？ 可能有更好的方法来做到这一点，但我已经定义了自己的函数来计算过去完整观察的数量。您可以定义一个函数nobs <- function(x) length(x[!is.na(x)]) 并将n() 替换为nobs(procras)。【参考方案2】：
我使用 gmodels 包中的 ci 命令：
library(gmodels)
your_db %>% group_by(gouping_variable1, grouping_variable2, ...)
%>% summarise(mean = ci(variable_of_interest)[1],
lowCI = ci(variable_of_interest)[2],
hiCI = ci(variable_of_interest)[3],
sd = ci (variable_of_interest)[4])
【讨论】：
这个ci()函数是从哪里来的？
抱歉，我没有提到 gmodels 包。我刚刚更新了回复。希望对你有帮助【参考方案3】：
如果你想使用boot 包的多功能性，我找到了this blog post useful（下面的代码灵感来自那里）
library(dplyr)
library(tidyr)
library(purrr)
library(boot)
set.seed(321)
mtcars %>%
group_by(vs) %>%
nest() %>%
mutate(boot_res = map(data,
~ boot(data = .$mpg,
statistic = function(x, i) mean(x[i]),
R = 1000)),
boot_res_ci = map(boot_res, boot.ci, type = "perc"),
mean = map(boot_res_ci, ~ .$t0),
lower_ci = map(boot_res_ci, ~ .$percent[[4]]),
upper_ci = map(boot_res_ci, ~ .$percent[[5]]),
n = map(data, nrow)) %>%
select(-data, -boot_res, -boot_res_ci) %>%
unnest(cols = c(n, mean, lower_ci, upper_ci)) %>%
ungroup()
#> # A tibble: 2 x 5
#> vs mean lower_ci upper_ci n
#> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 0 16.6 15.0 18.3 18
#> 2 1 24.6 22.1 27.3 14
由reprex package (v0.3.0) 于 2020 年 1 月 22 日创建
代码的一些解释：
与nest()嵌套时，会创建一个列表列（默认称为data），其中包含2个数据框，是整个mtcars的2个子集，按vs分组（包含2个唯一的值，0 和 1)。
然后，使用mutate() 和map()，我们将boot 包中的函数boot() 应用于列表列data，从而创建列表列boot_res。然后通过将boot.ci() 函数应用于boot_res 列表列等来创建boot_res_ci 列表列。
使用select()，我们删除不再需要的列表列，并通过取消嵌套和取消分组最终结果来休息。
不幸的是，该代码并不容易浏览，但它可以用于另一个示例。
使用broom::tidy()
刚刚意识到包broom 具有处理boot() 输出的方法的实现，正如here 指出的那样。这使得代码不那么冗长，输出更完整，包括统计的偏差和标准误差（这里的平均值）：
library(dplyr)
library(tidyr)
library(purrr)
library(broom)
library(boot)
set.seed(321)
mtcars %>%
group_by(vs) %>%
nest() %>%
mutate(boot_res = map(data,
~ boot(data = .$mpg,
statistic = function(x, i) mean(x[i]),
R = 1000)),
boot_tidy = map(boot_res, tidy, conf.int = TRUE, conf.method = "perc"),
n = map(data, nrow)) %>%
select(-data, -boot_res) %>%
unnest(cols = -vs) %>%
ungroup()
#> # A tibble: 2 x 7
#> vs statistic bias std.error conf.low conf.high n
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 0 16.6 -0.0115 0.843 15.0 18.3 18
#> 2 1 24.6 -0.0382 1.36 22.1 27.3 14
由reprex package (v0.3.0) 于 2020 年 1 月 22 日创建
data.table简洁的语法
但是请注意，通过使用 data.table 包而不是 dplyr，我得到了更简洁的语法：
library(data.table)
library(magrittr)
library(boot)
library(broom)
mtcars <- mtcars %>% copy %>% setDT
set.seed(321)
mtcars[, c(n = .N,
boot(data = mpg,
statistic = function(x, i) mean(x[i]),
R = 1000) %>%
tidy(conf.int = TRUE, conf.method = "perc")),
by = vs]
#> vs n statistic bias std.error conf.low conf.high
#> 1: 0 18 16.61667 -0.01149444 0.8425817 15.03917 18.26653
#> 2: 1 14 24.55714 -0.03822857 1.3633112 22.06429 27.32839
由reprex package (v0.3.0) 于 2020 年 1 月 23 日创建
使用 data.table 一次多个变量
library(data.table)
library(magrittr)
library(boot)
library(broom)
mtcars <- mtcars %>% copy %>% setDT
# Specify here the variables for which you want CIs
variables <- c("mpg", "disp")
# Function to get the CI stats, will be applied to each column of a subset of
# data (.SD)
get_ci <- function(varb, ...)
boot(data = varb,
statistic = function(x, i) mean(x[i]),
R = 1000) %>%
tidy(conf.int = TRUE, ...)
set.seed(321)
mtcars[, c(n = .N,
lapply(.SD, get_ci) %>%
rbindlist(idcol = "varb")),
by = vs, .SDcols = variables]
#> vs n varb statistic bias std.error conf.low conf.high
#> 1: 0 18 mpg 16.61667 -0.01149444 0.8425817 15.03917 18.26653
#> 2: 0 18 disp 307.15000 -1.49692222 23.1501247 261.18766 353.04416
#> 3: 1 14 mpg 24.55714 -0.03215714 1.3800432 21.86628 27.50551
#> 4: 1 14 disp 132.45714 0.32994286 14.9070552 104.45798 163.57344
由reprex package (v0.3.0) 于 2020-01-23 创建
【讨论】：
【参考方案4】：更新 tidyr 1.0.0
@Valentin 提供的所有解决方案都是可行的，但我想暗示一个新的替代方案，它对你们中的一些人来说更具可读性。它用一个名为unnest_wider 的相对较新的[tidyr 1.0.0][1] 函数替换了所有summarise 解决方案。
有了它，您可以将代码简化为以下内容：
mtcars %>%
nest(data = -"vs") %>%
mutate(ci = map(data, ~ MeanCI(.x$mpg, method = "boot", R = 1000))) %>%
unnest_wider(ci)
给出：
# A tibble: 2 x 5
vs data mean lwr.ci upr.ci
<dbl> <list> <dbl> <dbl> <dbl>
1 0 <tibble [18 × 10]> 16.6 14.7 18.5
2 1 <tibble [14 × 10]> 24.6 22.0 27.1
不用自举计算置信区间更简单：
mtcars %>%
nest(data = -"vs") %>%
mutate(ci = map(data, ~ MeanCI(.x$mpg))) %>%
unnest_wider(ci)
【讨论】：
【参考方案5】：对于正态分布：
library(dplyr)
mtcars %>%
group_by(vs) %>%
summarise(mean.mpg = mean(mpg, na.rm = TRUE),
sd.mpg = sd(mpg, na.rm = TRUE),
n.mpg = n()) %>%
mutate(se.mpg = sd.mpg / sqrt(n.mpg),
lower.ci.mpg = mean.mpg - qnorm(0.975) * se.mpg,
upper.ci.mpg = mean.mpg + qnorm(0.975) * se.mpg)
【讨论】：
以上是关于使用 dplyr 计算 95%-CI 的长度的主要内容，如果未能解决你的问题，请参考以下文章
如何计算两条 ROC 曲线之间的 AUC 差异（95% CI）？