有時候我們會想要分組計算一些統計量，所以會使用 group_by(GROUP_LABEL) 並結合 summarize(STAT = FUN(VARIABLE)) 來得到想要的分組後的統計量。例如這個 minimal example
library(dplyr)
library(tidyr)
df <- bind_cols(diversity_rank = sample(1:5, 1000, replace = TRUE),
dei_count_rank = sample(1:5, 1000, replace = TRUE),
metricA = rnorm(1000, 10, 5),
metricB = rnorm(1000, 10, 10),
metricC = rnorm(1000, 1, 5))
df %>% group_by(diversity_rank, dei_count_rank) %>%
summarize(meanA = mean(metricA),
meanB = mean(metricB),
meanC = mean(metricC))
我們知道樣本中有兩個 dimension 的 group label ，所以最多可能有 25 組。也就是最後得到的統計量（平均數）是這 25 組個別的平均。
如果我們要的操作稍微複雜一點，例如我想要得到固定 diversity_rank = 1，問 dei_count_rank 分別為 1 及 為 5 時的平均數之差，還有在統計上是否有顯著差異（t-test）？這個可以怎麼做比較好呢？然後如法炮製，問當 diversity_rank = 2時， dei_count_rank分別為 1 及 5 時的統計量。
想法上，我們可以先消掉一個 dimension，只看那些 diversity_rank == 1 的傢伙，然後再看 dei_count_rank %in% c(1,5)的傢伙。
df %>%
filter(diversity_rank == 1) %>%
filter(dei_count_rank %in% c(1,5))
# A tibble: 81 × 5
diversity_rank dei_count_rank metricA metricB metricC
<int> <int> <dbl> <dbl> <dbl>
1 1 1 8.38 14.8 4.18
2 1 1 13.6 11.0 2.60
3 1 5 14.9 3.96 4.41
4 1 1 4.43 4.66 9.91
5 1 1 10.2 -0.162 7.65
6 1 1 0.316 -2.56 7.71
7 1 1 15.3 11.4 -6.40
8 1 1 5.36 18.2 -3.69
9 1 1 1.66 1.55 4.28
10 1 5 20.3 15.9 -5.46
# … with 71 more rows
我們可能會很想要接著寫 group_by 然後summarize吧！不過這樣接在summarize 裡的 function 只會幫我們計算出 within group 的統計量，不是我想要做的對於 rank = 1 及 rank = 5的兩個組別的計算。如果這裡只有一個metric是我們關心的，那我們或許可以簡單地用 spread 來把不同的 dei_count_rank 建立為不同的 column （ from long data to wide data），但這裡麻煩的是我們有三個 metric 都想要計算 rank 1及5的 t-test，該怎麼辦可以用 dplyr 的 style 來完成呢？
此時其實我們依然可以照著我們原始的想法，繼續接續 group_by 跟 summarize 操作，只是巧妙地讓一個 column 的data type不再只是 dbl 或 int ，而是 list
df %>%
filter(diversity_rank == 1) %>%
filter(dei_count_rank %in% c(1,5)) %>%
group_by(dei_count_rank) %>%
summarize_at(.vars = c("metricA", "metricB", "metricC"), .funs = list) %>%
ungroup()
# A tibble: 2 × 4
dei_count_rank metricA metricB metricC
<int> <list> <list> <list>
1 1 <dbl [43]> <dbl [43]> <dbl [43]>
2 5 <dbl [38]> <dbl [38]> <dbl [38]>
既然有了 list column，實際上我們至此並沒有損失任何資訊，只是讓 tibble 的排列方式改變了一下。而且在外觀上，我們已經得到一個縮減後的 dataframe，有時更有利於去想像最終 output 的樣子。這樣我們就可以得到一個用於計算最後統計量的中間產物。
此時為了簡化 case ，我們先專注在只針對 metricA 這一變數（要泛化的話只要針對不同的變數改變此後的步驟，然後bind_cols即可），求 diversity_rank 在 1-5 時，dei_count_rank 為 1 及 為 5 之差還有 t-test p-value。
注意到我們的例子現在是這樣：
df %>%
filter(dei_count_rank %in% c(1,5)) %>%
group_by(diversity_rank, dei_count_rank) %>%
summarize(metricA = list(na.omit(metricA))) %>%
ungroup()
# A tibble: 10 × 3
diversity_rank dei_count_rank metricA
<int> <int> <list>
1 1 1 <dbl [43]>
2 1 5 <dbl [38]>
3 2 1 <dbl [35]>
4 2 5 <dbl [40]>
5 3 1 <dbl [43]>
6 3 5 <dbl [38]>
7 4 1 <dbl [42]>
8 4 5 <dbl [44]>
9 5 1 <dbl [37]>
10 5 5 <dbl [50]>
實際上每2個row為即將要進行比較的跨組資料。我們喜歡用 mutate跨 column 地比較資料，所以可以用 spread 把 dei_count_rank 到底是 1 還是 5 放到 column 上去（注意，這裡 dei_count_rank 實際上可以超過兩個組別，也是一樣的作法）
我們利用 spread(key, value) 後就得到：
df %>%
filter(dei_count_rank %in% c(1,5)) %>%
group_by(diversity_rank, dei_count_rank) %>%
summarize(metricA = list(na.omit(metricA))) %>%
ungroup() %>%
spread(key = dei_count_rank, value = metricA)
# A tibble: 5 × 3
diversity_rank `1` `5`
<int> <list> <list>
1 1 <dbl [43]> <dbl [38]>
2 2 <dbl [35]> <dbl [40]>
3 3 <dbl [43]> <dbl [38]>
4 4 <dbl [42]> <dbl [44]>
5 5 <dbl [37]> <dbl [50]>
接著其實就是組內的運算了，只是記得 unlist()
df %>%
filter(dei_count_rank %in% c(1,5)) %>%
group_by(diversity_rank, dei_count_rank) %>%
summarize(metricA = list(na.omit(metricA))) %>%
ungroup() %>%
spread(key = dei_count_rank, value = metricA) %>%
group_by(diversity_rank) %>%
mutate(diff = mean(unlist(`1`)) - mean(unlist(`5`)),
p.value = t.test(unlist(`1`), unlist(`5`))$p.value)
# A tibble: 5 × 5
# Groups: diversity_rank [5]
diversity_rank `1` `5` diff p.value
<int> <list> <list> <dbl> <dbl>
1 1 <dbl [43]> <dbl [38]> -0.643 0.553
2 2 <dbl [35]> <dbl [40]> 0.611 0.654
3 3 <dbl [43]> <dbl [38]> 0.741 0.564
4 4 <dbl [42]> <dbl [44]> -0.982 0.335
5 5 <dbl [37]> <dbl [50]> 1.60 0.147
參考：
https://sebastiansauer.github.io/multiple-t-tests-with-dplyr/