使用 dplyr 计算 95%-CI 的长度
Calculating length of 95%-CI using dplyr
上次我问如何计算一个变量 (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 table 中添加一列，其中包括每次事件平均得分周围 95%-CI 的长度。
http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/ 解释了如何获取和绘制 CIs，但是一旦我想对任何子组执行此操作，这种方法似乎就会出现问题，对吗？那么有没有办法让 dplyr 也自动在 mean_data 中包含 CI（基于组大小等）？
之后，我希望将新值作为 CIs 绘制到图表中应该相当容易。
谢谢。
您可以使用 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
我使用 ci 命令来自 gmodels 包：
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])
如果你想使用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() 输出，正如 指出的那样。这使得代码不那么冗长，输出更完整，包括统计数据的偏差和标准误差（此处为平均值）：
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 年 1 月 23 日创建
更新 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)
正态分布：
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)
上次我问如何计算一个变量 (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 table 中添加一列，其中包括每次事件平均得分周围 95%-CI 的长度。
http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/ 解释了如何获取和绘制 CIs，但是一旦我想对任何子组执行此操作，这种方法似乎就会出现问题，对吗？那么有没有办法让 dplyr 也自动在 mean_data 中包含 CI（基于组大小等）？ 之后，我希望将新值作为 CIs 绘制到图表中应该相当容易。 谢谢。
您可以使用 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
我使用 ci 命令来自 gmodels 包：
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])
如果你想使用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() 输出，正如
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 年 1 月 23 日创建更新 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)
正态分布：
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)