【问题标题】:Calculating length of 95%-CI using dplyr使用 dplyr 计算 95%-CI 的长度
【发布时间】: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 绘制到我希望的图表中应该相当容易。 谢谢。

【问题讨论】:

    标签: r ggplot2 linechart confidence-interval trend


    【解决方案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 &lt;- 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)
        

        【讨论】:

          猜你喜欢
          • 2019-06-26
          • 1970-01-01
          • 1970-01-01
          • 1970-01-01
          • 1970-01-01
          • 2021-11-19
          • 2021-05-30
          • 2020-11-23
          • 1970-01-01
          相关资源
          最近更新 更多