ddply 总结置信区间？
ddply summarise with Confidence interval?
我正在尝试使用 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"s 和 "BSTN"s 分组的 Ttime 和 TFtime 的值。
对于"Ttime"和"TFtime"的平均值我想反映95%的置信区间。所以计算[=的平均值28=] 和 "TFtime"s 在 95% 边界内。 如果 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))
使用包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，但是当您拥有更大的数据集时，您很少会获得它们。
更新答案
我不确定，但我猜你真正想做的是过滤所有大于/小于 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-06-23 创建
旧答案
有了更好的示例数据，我们可以实现更具程序化的方法。例如，我使用 ggplot2s diamonds 数据集。在下面的代码中查看我的评论。
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-06-19 创建
如果您想根据均值的置信度 iIntervals 过滤观察结果，您可以按如下方式调整我上面的方法。请注意，这与过滤每个变量的顶部和底部 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 日创建
我正在尝试使用 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"s 和 "BSTN"s 分组的 Ttime 和 TFtime 的值。
对于"Ttime"和"TFtime"的平均值我想反映95%的置信区间。所以计算[=的平均值28=] 和 "TFtime"s 在 95% 边界内。 如果 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))
使用包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，但是当您拥有更大的数据集时，您很少会获得它们。
更新答案
我不确定，但我猜你真正想做的是过滤所有大于/小于 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-06-23 创建旧答案
有了更好的示例数据，我们可以实现更具程序化的方法。例如，我使用 ggplot2s diamonds 数据集。在下面的代码中查看我的评论。
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-06-19 创建如果您想根据均值的置信度 iIntervals 过滤观察结果，您可以按如下方式调整我上面的方法。请注意，这与过滤每个变量的顶部和底部 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 日创建