ddply 用置信区间总结？答案

【问题标题】：ddply summarise with Confidence interval?ddply 用置信区间总结？
【发布时间】：2020-06-19 05:33:56
【问题描述】：

我正在尝试使用 ddply 总结我的数据，并且我正在尝试找到一种方法来总结数据同时反映可靠性。

这是我的数据集的描述。

  BSTN ASTN  BSEC ASTN1 BSTN2 ASTN2 BSTN3 ASTN3 BSTN4 ASTN4 BSTN5 TFtime Ttime    ID
1 1001 1003 69551  1703  1703     0     0     0     0     0     0    399  2933 35404
2 1001 1006 69664  1703  1703     0     0     0     0     0     0    399  2284 35405
3 1001 1701 66606  1703  1703     0     0     0     0     0     0    118  1750 35406
4 1001 1701 66600  1703  1703     0     0     0     0     0     0    118  1750 35406
5 1001 1701 66601  1703  1703     0     0     0     0     0     0    118  1750 35406
6 1001 1703 69434     0     0     0     0     0     0     0     0      0  1005 35407

我想要的输出是总结按“ASTN”和“BSTN”分组的 Ttime 和 TFtime 的值。

对于 "Ttime" 和 "TFtime" 的平均值，我想反映 95% 的置信区间。 所以计算 95 内的 "Ttime" 和 "TFtime" 的平均值% 边界。 如果有多个 BSTN~ASTN 组合，我将如何使用 ddply 执行此过程。

以下是我使用并希望修改的代码。

Routetable<-ddply(A,c(.(BSTN),.(ASTN1),.(BSTN2),.(ASTN2),.(BSTN3),.(ASTN3),.(BSTN4),.(ASTN4),.(BSTN5),.(ASTN)), 
                    summarise, count=length(BSTN),mean=mean(Ttime),TFtimemean=mean(TFtime))

【问题讨论】：

标签： r dplyr plyr confidence-interval

【解决方案1】：

更新答案

我不确定，但我猜你真正想要做的是过滤所有大于/小于mean(x) -/+ 2*sd(x) 的值，并按每个组过滤。以下方法可以做到这一点。在ggplot2s Diamond 数据集的情况下，它保留了所有值的大约 97%，并且只是去除了极端值。

library(tidyverse)

diamonds %>% 
  group_by(cut, color) %>% 
  mutate(across(c(x,y,z),
                list(low = ~ mean(.x, na.rm = TRUE) - 2 * sd(.x, na.rm = TRUE),
                     high = ~ mean(.x, na.rm = TRUE) + 2 * sd(.x, na.rm = TRUE))
                )
         ) %>% 
  filter(x >= x_low & x <= x_high,
         y >= x_low & y <= y_high,
         z >= z_low & z <= z_high)
#> # A tibble: 52,299 x 16
#> # Groups:   cut, color [35]
#>    carat cut   color clarity depth table price     x     y     z x_low x_high
#>    <dbl> <ord> <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl>  <dbl>
#>  1 0.23  Ideal E     SI2      61.5    55   326  3.95  3.98  2.43  3.51   6.92
#>  2 0.21  Prem~ E     SI1      59.8    61   326  3.89  3.84  2.31  3.52   7.65
#>  3 0.290 Prem~ I     VS2      62.4    58   334  4.2   4.23  2.63  3.86   9.12
#>  4 0.31  Good  J     SI2      63.3    58   335  4.34  4.35  2.75  4.14   8.62
#>  5 0.24  Very~ I     VVS1     62.3    57   336  3.95  3.98  2.47  3.92   8.62
#>  6 0.26  Very~ H     SI1      61.9    55   337  4.07  4.11  2.53  3.66   8.30
#>  7 0.23  Very~ H     VS1      59.4    61   338  4     4.05  2.39  3.66   8.30
#>  8 0.3   Good  J     SI1      64      55   339  4.25  4.28  2.73  4.14   8.62
#>  9 0.23  Ideal J     VS1      62.8    56   340  3.93  3.9   2.46  3.88   8.76
#> 10 0.31  Ideal J     SI2      62.2    54   344  4.35  4.37  2.71  3.88   8.76
#> # ... with 52,289 more rows, and 4 more variables: y_low <dbl>, y_high <dbl>,
#> #   z_low <dbl>, z_high <dbl>

^{由reprex package (v0.3.0) 于 2020 年 6 月 23 日创建}

旧答案

借助更好的示例数据，我们可以实现更加程序化的方法。例如，我使用 ggplot2s diamonds 数据集。在下面的代码中查看我的 cmets。

library(tidyverse)

diamonds %>% 
  # set up your groups
  nest_by(cut, color) %>%  
  # define in `across` for which variables you want means and conf int to be calculated 
  mutate(ttest = list(summarise(data, across(c(x,y,z), 
                                          ~ broom::tidy(t.test(.x))))),
         ttest = list(unpack(ttest, c(x, y, z), names_sep = "_") %>% 
   # select only the estimates and conf intervalls                          
                        select(contains("estimate"), contains("conf")))) %>% 
  unnest(ttest)
#> # A tibble: 35 x 12
#> # Groups:   cut, color [35]
#>    cut   color      data x_estimate y_estimate z_estimate x_conf.low x_conf.high
#>    <ord> <ord> <list<tb>      <dbl>      <dbl>      <dbl>      <dbl>       <dbl>
#>  1 Fair  D     [163 × 8]       6.02       5.96       3.84       5.89        6.15
#>  2 Fair  E     [224 × 8]       5.91       5.86       3.72       5.80        6.02
#>  3 Fair  F     [312 × 8]       5.99       5.93       3.79       5.89        6.09
#>  4 Fair  G     [314 × 8]       6.17       6.11       3.96       6.06        6.28
#>  5 Fair  H     [303 × 8]       6.58       6.50       4.22       6.47        6.69
#>  6 Fair  I     [175 × 8]       6.56       6.49       4.19       6.43        6.70
#>  7 Fair  J     [119 × 8]       6.75       6.68       4.32       6.55        6.95
#>  8 Good  D     [662 × 8]       5.62       5.63       3.50       5.55        5.69
#>  9 Good  E     [933 × 8]       5.62       5.63       3.50       5.56        5.68
#> 10 Good  F     [909 × 8]       5.69       5.71       3.54       5.63        5.76
#> # … with 25 more rows, and 4 more variables: y_conf.low <dbl>,
#> #   y_conf.high <dbl>, z_conf.low <dbl>, z_conf.high <dbl>

^{由reprex package (v0.3.0) 于 2020 年 6 月 19 日创建}

如果您想根据均值的置信度 iIntervalls 过滤观察结果，您可以调整我上面的方法，如下所示。请注意，这与过滤每个变量的顶部和底部 2.5 % 不同，您会丢失大量数据。

library(tidyverse)

diamonds %>% 
  nest_by(cut, color) %>% 
  mutate(ttest = summarise(data, across(c(x,y,z), 
                                             ~ broom::tidy(t.test(.x)))) %>% 
         unpack(c(x,y,z), names_sep = "_")) %>% 
  unpack(ttest) %>% 
  select(cut, color, data, contains("estimate"), contains("conf")) %>% 
  rowwise(cut, color) %>% 
  mutate(data = list(filter(data,
                       x >= x_conf.low & x <= x_conf.high,
                       y >= x_conf.low & y <= y_conf.high,
                       z >= z_conf.low & z <= z_conf.high))) %>% 
  unnest(data)
#> # A tibble: 322 x 19
#> # Groups:   cut, color [30]
#>    cut   color carat clarity depth table price     x     y     z x_estimate
#>    <ord> <ord> <dbl> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>      <dbl>
#>  1 Fair  D      0.91 SI2      62.5    66  3079  6.08  6.01  3.78       6.02
#>  2 Fair  D      0.9  SI2      65.7    60  3205  5.98  5.93  3.91       6.02
#>  3 Fair  D      0.9  SI2      64.7    59  3205  6.09  5.99  3.91       6.02
#>  4 Fair  D      0.95 SI2      64.4    60  3384  6.06  6.02  3.89       6.02
#>  5 Fair  D      0.9  SI2      64.9    57  3473  6.03  5.98  3.9        6.02
#>  6 Fair  D      0.9  SI2      64.5    61  3473  6.1   6     3.9        6.02
#>  7 Fair  D      0.9  SI1      64.5    61  3689  6.05  6.01  3.89       6.02
#>  8 Fair  D      0.91 SI1      64.7    61  3730  6.06  5.99  3.9        6.02
#>  9 Fair  D      0.9  SI2      64.6    59  3847  6.04  6.01  3.89       6.02
#> 10 Fair  D      0.91 SI1      64.4    60  3855  6.08  6.04  3.9        6.02
#> # ... with 312 more rows, and 8 more variables: y_estimate <dbl>,
#> #   z_estimate <dbl>, x_conf.low <dbl>, x_conf.high <dbl>, y_conf.low <dbl>,
#> #   y_conf.high <dbl>, z_conf.low <dbl>, z_conf.high <dbl>

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

【讨论】：

感谢您的解释。所以如果我是正确的，下面是过滤的结果？我想要的最终输出是减少 CI 中不存在的数据。我想减少这些数据并根据 BSTN 和 ASTN 进行汇总。与平均时间相比，有些值更大或更小，这就是我要过滤的原因
我想要的最终输出是消除不在上下边界的数据，然后用平均Ttime和TFtime汇总数据。
也许我更新的答案可以解决您的问题。我不确定。

【解决方案2】：

使用包dplyr（比plyr更新）你可以写如下。 “LB”和“UB”分别代表“Lower Bound”和“Upper Bound”。

library(dplyr)

A %>% 
  group_by(across(starts_with("BSTN") | starts_with("ASTN"))) %>% 
  summarise(
    count = n(),
    mean_Ttime = mean(Ttime),
    mean_TFtime = mean(TFtime),
    LB_Ttime = mean_Ttime - qnorm(0.975) * sd(Ttime) / sqrt(count),
    UB_Ttime = mean_Ttime + qnorm(0.975) * sd(Ttime) / sqrt(count),
    LB_TFtime = mean_TFtime - qnorm(0.975) * sd(TFtime) / sqrt(count),
    UB_TFtime = mean_TFtime + qnorm(0.975) * sd(TFtime) / sqrt(count)
  )

输出

# A tibble: 4 x 17
# Groups:   BSTN, BSTN2, BSTN3, BSTN4, BSTN5, ASTN, ASTN1, ASTN2, ASTN3 [4]
#    BSTN BSTN2 BSTN3 BSTN4 BSTN5  ASTN ASTN1 ASTN2 ASTN3 ASTN4 count mean_Ttime mean_TFtime LB_Ttime UB_Ttime LB_TFtime UB_TFtime
#   <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>      <dbl>       <dbl>    <dbl>    <dbl>     <dbl>     <dbl>
# 1  1001     0     0     0     0  1703     0     0     0     0     1       1005           0       NA       NA        NA        NA
# 2  1001  1703     0     0     0  1003  1703     0     0     0     1       2933         399       NA       NA        NA        NA
# 3  1001  1703     0     0     0  1006  1703     0     0     0     1       2284         399       NA       NA        NA        NA
# 4  1001  1703     0     0     0  1701  1703     0     0     0     3       1750         118     1750     1750       118       118

通过这个样本数据，我们获得了几个 NA，因为在这些情况下，该组的 count 为 1，但是当您拥有更大的数据集时，您很少会获得它们。

【讨论】：

谢谢！！抱歉回复晚了，那么最终输出行是否在置信区间边界内？
对不起，你到底是什么意思？我不明白你的问题
抱歉解释不好。那么应用此代码后的剩余数据是否在 0.975% 的边界内？我想消除具有置信区间的旅行时间和转移时间的数据
这些不再是您原来的行了。它们是原始行的摘要，按所有 BSTN 和 ASTN 变量分组。
@Yun Hynsoo：您是要删除低于/高于均值置信区间的观测值，还是只想削减前 2.5%数据？存在很大差异，对于后者，我们不需要计算均值的置信区间。如果您根据均值的 CI 删除观察值，您将丢失大量数据。