使用dplyr计算95%置信区间的长度

Question

使用dplyr计算95%置信区间的长度

rggplot2linechartconfidence-intervaltrend

16

上次我问如何计算在多个被试中重复测量某个变量（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％置信区间的长度。 http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/解释了如何获得和绘制置信区间，但是对于任何子组进行此操作似乎会出现问题，对吧？因此，有没有办法让dplyr自动将CI（基于组大小等）包含在mean_data中呢？然后，希望将新值作为置信区间轻松地绘制到图表中。谢谢。

- Rasul89

6个回答

14

我使用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])

- carfisma

4

如果您想使用boot包的多功能性，我发现这篇博客文章很有用（下面的代码受到其启发）

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-01-22创建}

代码解释:

当使用nest()时，默认会创建一个名为data的列表列，其中包含两个数据框，分别是mtcars按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

^{2020年01月22日由reprex package (v0.3.0)创建}

`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包 (v0.3.0)于2020-01-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日创建的。}

- Valentin_Ștefan

2

更新tidyr 1.0.0

@Valentin提供的所有解决方案都是可行的，但我想提示一种新的替代方案，它对于一些人来说更易读。它使用一个相对较新的[tidyr 1.0.0][1]函数 called unnest_wider来替换所有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

如果不使用bootstrap方法计算置信区间，可以更简单地使用以下方法：

mtcars %>% 
  nest(data = -"vs") %>%
  mutate(ci = map(data, ~ MeanCI(.x$mpg))) %>% 
  unnest_wider(ci)

- Agile Bean

0

为了帮助其他像我一样发现这篇文章有用但仍在寻找调整方法的人，我在此添加一个答案。

这里提供了一种基于@sboysel回答并使用他的“nobs”函数的替代手动解决方案。如果您想要在数据内部的组之间以及跨越多个变量进行总结（更改 across()以适应您的数据-在这里，它被编码为以特定字符串开头的变量），则此解决方案非常有用：

output1 <- your_data_frame %>% 
  dplyr::group_by(your_grouping_variable) %>% 
  dplyr::summarise(across(starts_with("your_string"),
                          .fns = list(
                            mean = ~mean(.x, na.rm = TRUE), 
                            sd = ~sd(.x, na.rm = TRUE), 
                            se = ~sd(.x, na.rm = TRUE)/sqrt(length(.x)),
                            n = ~nobs(.x),
                            ci_l = ~mean(.x, na.rm = TRUE) - (1.96 * sd(.x, na.rm = TRUE)/sqrt(nobs(.x))),
                            ci_u = ~mean(.x, na.rm = TRUE) + (1.96 * sd(.x, na.rm = TRUE)/sqrt(nobs(.x))))))

或者，如@carfisma所说，使用gmodels包中的ci来编写更简洁的代码：

output2 <- your_data_frame%>% 
  dplyr::group_by(your_grouping_variable) %>% 
  dplyr::summarise(across(starts_with("your_string"),
                          .fns = list(
                            mean = ~ci(.x, na.rm=TRUE)[1],
                            se = ~ci(.x, na.rm=TRUE)[4],
                            n = ~nobs(.x),
                            ci_l = ~ci(.x, na.rm=TRUE)[2], # confidence level default is 0.95
                            ci_u = ~ci(.x, na.rm=TRUE)[3])))

请注意，ci（）输出的第4个元素为标准误差，而不是像carfisma解决方案中所暗示的那样为sd。

使用dplyr版本1.0.10和gmodels 2.18.1.1

- user39683

0

对于正态分布：

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)

- SysEng

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- sboysel · Accepted Answer

您可以使用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

使用dplyr计算95%置信区间的长度

使用 broom::tidy()

data.table简洁语法

使用data.table一次处理多个变量

更新tidyr 1.0.0

使用 `broom::tidy()`

`data.table`简洁语法