使用 dplyr/tidyverse 同时对多个变量进行多对 t 检验答案

【问题标题】：Multiple paired t-tests on multiple variables simultaneously using dplyr/tidyverse使用 dplyr/tidyverse 同时对多个变量进行多对 t 检验
【发布时间】：2019-03-08 18:12:32
【问题描述】：

假设这样的数据结构：

   ID testA_wave1 testA_wave2 testA_wave3 testB_wave1 testB_wave2 testB_wave3
1   1           3           2           3           6           5           3
2   2           4           4           4           3           6           6
3   3          10           2           1           4           4           4
4   4           5           3          12           2           7           4
5   5           5           3           9           2           4           2
6   6          10           0           2           6           6           5
7   7           6           8           4           6           8           3
8   8           1           5           4           5           6           0
9   9           3           2           7           8           4           4
10 10           4           9           5          11           8           8

我想要实现的是分别为每个测试计算配对 t 检验（在这种情况下意味着 testA 和 testB，但在现实生活中我有更多的测试）。我想这样做，我将给定测试的第一波与同一测试的所有其他后续波进行比较（在 testA 的情况下是 testA_wave1 与 testA_wave2 和 testA_wave1 与 testA_wave3）。

这样，我就可以实现了：

df %>%
 gather(variable, value, -ID) %>%
 mutate(wave_ID = paste0("wave", parse_number(variable)),
        variable = ifelse(grepl("testA", variable), "testA",
                     ifelse(grepl("testB", variable), "testB", NA_character_))) %>%
 group_by(wave_ID, variable) %>% 
 summarise(value = list(value)) %>% 
 spread(wave_ID, value) %>% 
 group_by(variable) %>% 
 mutate(p_value_w1w2 = t.test(unlist(wave1), unlist(wave2), paired = TRUE)$p.value,
        p_value_w1w3 = t.test(unlist(wave1), unlist(wave3), paired = TRUE)$p.value) %>%
 select(variable, matches("(p_value)"))

  variable p_value_w1w2 p_value_w1w3
  <chr>           <dbl>        <dbl>
1 testA           0.664        0.921
2 testB           0.146        0.418

但是，我希望看到不同/更优雅的解决方案，它们会产生类似的结果。我主要在寻找dplyr/tidyverse 的解决方案，但如果有完全不同的方法来实现它，我不反对。

样本数据：

set.seed(123)
df <- data.frame(ID = 1:20,
testA_wave1 = round(rnorm(20, 5, 3), 0),
testA_wave2 = round(rnorm(20, 5, 3), 0),
testA_wave3 = round(rnorm(20, 5, 3), 0),
testB_wave1 = round(rnorm(20, 5, 3), 0),
testB_wave2 = round(rnorm(20, 5, 3), 0),
testB_wave3 = round(rnorm(20, 5, 3), 0))

【问题讨论】：

testA_wave2 和 testA_wave3 怎么样？
@JasonAizkalns 我只对第一波与所有后续波的比较感兴趣。因此，不需要 testA_wave2 与 testA_wave3。
这更像是一个统计问题，但与编程相关 - 您正在运行大量配对 t 检验，所以您会针对多重比较进行调整吗？您可以进行更有效的纵向分析吗？
@Mike 我知道你的意思，但在这种情况下不需要更正。
关于@Mike 的评论（对于登陆这里但没有关注它的任何人）：stats.stackexchange.com/questions/88065/…

标签： r dplyr

【解决方案1】：

由于dplyr 0.8.0 我们可以使用group_split 将数据帧拆分为数据帧列表。

我们gather 数据帧并将其转换为长格式，然后将separate 列的名称（key）转换为不同的列（test 和wave）。然后我们使用group_split 将数据框拆分为基于test 列的列表。对于列表中的每个数据帧，我们将spread 转换为宽格式，然后计算t.test 值并使用map_dfr 将它们绑定到一个数据帧中。

library(tidyverse)

df %>%
  gather(key, value, -ID) %>%
  separate(key, c("test", "wave")) %>%
  group_split(test) %>% #Previously we had to do split(.$test) here
  map_dfr(. %>%
          spread(wave, value) %>%
          summarise(test = first(test),
                    p_value_w1w2 = t.test(wave1, wave2, paired = TRUE)$p.value, 
                    p_value_w1w3 = t.test(wave1, wave3, paired = TRUE)$p.value))


# A tibble: 2 x 3
#  test  p_value_w1w2 p_value_w1w3
#  <chr>        <dbl>        <dbl>
#1 testA        0.664        0.921
#2 testB        0.146        0.418

我们手动执行上面的 t 检验，因为只有 2 个值需要计算。如果有更多数量的 wave... 列，那么这可能会变得很麻烦。在这种情况下，我们可以这样做

df %>%
   gather(key, value, -ID) %>%
   separate(key, c("test", "wave")) %>%
   group_split(test) %>% 
   map_dfr(function(data) 
              data %>%
                   spread(wave, value) %>%
                   summarise_at(vars(setdiff(unique(data$wave), "wave1")), 
                   function(x) t.test(.$wave1, x, paired = TRUE)$p.value) %>%
                   mutate(test = first(data$test)))

#  wave2 wave3 test 
#  <dbl> <dbl> <chr>
#1 0.664 0.921 testA
#2 0.146 0.418 testB

这里它将对每个带有“wave1”列的“wave..”列执行 t 检验。

由于您也对其他解决方案持开放态度，因此这里尝试使用纯基础 R 解决方案

sapply(split.default(df[-1], sub("_.*", "", names(df[-1]))), function(x) 
 c(p_value_w1w2 = t.test(x[[1]], x[[2]],paired = TRUE)$p.value, 
   p_value_w1w3 = t.test(x[[1]], x[[3]],paired = TRUE)$p.value))


#                 testA     testB
#p_value_w1w2 0.6642769 0.1456059
#p_value_w1w3 0.9209554 0.4184603

我们根据test* 拆分列并创建数据框列表，并将t.test 应用于每个数据框的不同列组合。

【讨论】：

很好，学到了一些东西group_split。但是，需要相当多的硬编码，每次测试应该有 2 波以上。
@thothal 你是对的，但由于它在 OP 的示例中是硬编码的，我认为他们只有 2 个值要计算，但无论如何最好有一个可扩展的解决方案。谢谢，我已经更新了适用于许多“wave”列的答案。 :)

【解决方案2】：

2022 年 3 月 16 日更新

tidyverse 已经发展，这个解决方案也应该如此。

首先我做一个简单的假设：如果我们设计了实验，那么我们就知道这些组是什么以及我们跟随它们经过了多少波。如果我们不知道，那么我们可以从列名中提取此信息。见下文。

library("broom")
library("tidyverse")

tests <- c("A", "B")
waves <- 3

comparisons <-
  list(
    test = tests,
    first = 1,
    later = seq(2, waves)
  ) %>%
  cross_df()
comparisons
#> # A tibble: 4 × 3
#>   test  first later
#>   <chr> <dbl> <int>
#> 1 A         1     2
#> 2 B         1     2
#> 3 A         1     3
#> 4 B         1     3

将数据从宽格式转换为长格式。

data <- df %>%
  pivot_longer(
    -ID,
    names_to = "test_wave"
  ) %>%
  extract(
    test_wave, c("test", "wave"),
    regex = "test(.+)_wave(.+)",
    convert = TRUE
  )

然后将我们想要进行的比较与我们收集的数据配对。我添加了许多重命名语句以使代码更具可读性，但这并不是绝对必要的。

comparisons %>%
  inner_join(
    data,
    by = c("test", "first" = "wave")
  ) %>%
  rename(
    value.first = value
  ) %>%
  inner_join(
    data,
    by = c("test", "later" = "wave", "ID")
  ) %>%
  rename(
    value.later = value
  ) %>%
  group_by(
    test, first, later
  ) %>%
  group_modify(
    ~ tidy(t.test(.x$value.first, .x$value.later, paired = TRUE))
  ) %>%
  ungroup() %>%
  pivot_wider(
    id_cols = test,
    names_from = later,
    names_glue = "wave1_vs_wave{later}",
    values_from = p.value
  )
#> # A tibble: 2 × 3
#>   test  wave1_vs_wave2 wave1_vs_wave3
#>   <chr>          <dbl>          <dbl>
#> 1 A              0.664          0.921
#> 2 B              0.146          0.418

附录：从列名中提取测试名称和波数。

design <- df %>%
  select(starts_with("test")) %>%
  colnames() %>%
  str_match("test(.+)_wave(.+)")
tests <- unique(design[, 2])
waves <- max(as.integer(design[, 3]))

^{由reprex package 创建于 2022-03-16 (v2.0.1)}

旧解决方案

这是一种方法，使用purrr 相当多。

library("tidyverse")

set.seed(123)
df <- tibble(
  ID = 1:20,
  testA_wave1 = round(rnorm(20, 5, 3), 0),
  testA_wave2 = round(rnorm(20, 5, 3), 0),
  testA_wave3 = round(rnorm(20, 5, 3), 0),
  testB_wave1 = round(rnorm(20, 5, 3), 0),
  testB_wave2 = round(rnorm(20, 5, 3), 0),
  testB_wave3 = round(rnorm(20, 5, 3), 0)
)

pvalues <- df %>%
  # From wide tibble to long tibble
  gather(test, value, -ID) %>%
  separate(test, c("test", "wave")) %>%
  # Not stricly necessary; will order the waves alphabetically instead
  mutate(wave = parse_number(wave)) %>%
  inner_join(., ., by = c("ID", "test")) %>%
  # If there are two waves w1 and w2,
  # we end up with pairs (w1, w1), (w1, w2), (w2, w1) and (w2, w2),
  # so filter out to keep the pairing (w1, w2) only
  filter(wave.x == 1, wave.x < wave.y) %>%
  nest(ID, value.x, value.y) %>%
  mutate(pvalue = data %>%
           # Perform the test
           map(~t.test(.$value.x, .$value.y, paired = TRUE)) %>%
           map(broom::tidy) %>%
           # Also not strictly necessary; you might want to keep all
           # information about the test: estimate, statistic, etc.
           map_dbl(pluck, "p.value"))
pvalues
#> # A tibble: 4 x 5
#>   test  wave.x wave.y data              pvalue
#>   <chr>  <dbl>  <dbl> <list>             <dbl>
#> 1 testA      1      2 <tibble [20 x 3]>  0.664
#> 2 testA      1      3 <tibble [20 x 3]>  0.921
#> 3 testB      1      2 <tibble [20 x 3]>  0.146
#> 4 testB      1      3 <tibble [20 x 3]>  0.418

pvalues %>%
  # Drop the data in order to pivot the table
  select(- data) %>%
  unite("waves", wave.x, wave.y, sep = ":") %>%
  spread(waves, pvalue)
#> # A tibble: 2 x 3
#>   test  `1:2` `1:3`
#>   <chr> <dbl> <dbl>
#> 1 testA 0.664 0.921
#> 2 testB 0.146 0.418

^{由reprex package (v0.2.1) 于 2019 年 3 月 8 日创建}

【讨论】：

感谢您的解决方案，但是，它也提供了 wave2 与 wave3 的冗余比较。我想避免这样的步骤。
@tmfmnk 只需添加另一个过滤条件：wave.x == 1 应该覆盖它...所以filter(wave.x < wave.y, wave.x == 1)
是的。更改了过滤条件以比较第一波和以后的波。

【解决方案3】：

抛出data.table 解决方案：

library(stringr)
library(data.table)
library(magrittr) ## for the pipe operator

dt_sol <- function(df) {
  ## create patterns for the melt operation:
  ## all columns from the same wave should go in one column
  grps <- str_extract(names(df)[-1], 
                      "[0-9]+$") %>%
    unique() %>%
    paste0("wave", ., "$")
  grp_names <- sub("\\$", "", grps)
  ## melt the data table: all test*_wave_i data go into column wave_i
  df.m <- melt(df, 
               measure = patterns(grps),
               value.name = grp_names,
               variable.name = "test")
  ## define the names for the new column, we want to extract estimate and p.value
  new_cols <- c(outer(c("p.value", "estimate"), 
                      grp_names[-1],
                      paste, sep = "_"))
  ## use lapply on .SD which equals to all wave_i columns but the first one
  ## return estimate and p.value
  df.m[, 
       setNames(unlist(lapply(.SD, 
                              function(col) {
                                t.test(wave1, col, paired = TRUE)[c("p.value", "estimate")]
                              }), recursive = FALSE), new_cols),
       test, ## group by each test
       .SDcols = grp_names[-1]] 
}
dt <- copy(df)
setDT(dt)
dt_sol(dt)
#    test p.value_wave2 estimate_wave2 p.value_wave3 estimate_wave3
# 1:    1     0.6642769           0.40     0.9209554           -0.1
# 2:    2     0.1456059          -1.45     0.4184603            0.7

基准测试

将data.table 解决方案与tidyverse 解决方案进行比较，我们得到data.tablesolution 的3 倍速度提升：

dp_sol <- function(df) {
  df %>%
    gather(test, value, -ID) %>%
    separate(test, c("test", "wave")) %>%
    inner_join(., ., by = c("ID", "test")) %>%
    filter(wave.x == 1, wave.x < wave.y) %>%
    nest(ID, value.x, value.y) %>%
    mutate(pvalue = data %>%
             map(~t.test(.$value.x, .$value.y, paired = TRUE)) %>%
             map(broom::tidy) %>%
             map_dbl(pluck, "p.value"))
}

library(microbenchmark)

microbenchmark(dplyr = dp_sol(df),
               data.table = dt_sol(dt))


# Unit: milliseconds
#        expr      min       lq     mean   median       uq       max neval cld
#       dplyr 6.119273 6.897456 7.639569 7.348364 7.996607 14.938182   100   b
#  data.table 1.902547 2.307395 2.790910 2.758789 3.133091  4.923153   100  a

输入稍大：

make_df <- function(nr_tests = 2,
                    nr_waves = 3,
                    n_per_wave = 20) {
  mat <- cbind(seq(1, n_per_wave),
               matrix(round(rnorm(nr_tests * nr_waves * n_per_wave), 0),
                      nrow = n_per_wave))
  c_names <- c(outer(1:nr_waves, 1:nr_tests, function(w, t) glue::glue("test{t}_wave{w}")))
  colnames(mat) <- c("ID", c_names)
  as.data.frame(mat)
}

df2 <- make_df(100, 100, 10)
dt2 <- copy(df2)
setDT(dt2)

microbenchmark(dplyr = dp_sol(df2),
               data.table = dt_sol(dt2)

# Unit: seconds
#        expr      min       lq     mean   median       uq      max neval cld
#       dplyr 3.469837 3.669819 3.877548 3.821475 3.984518 5.268596   100   b
#  data.table 1.018939 1.126244 1.193548 1.173175 1.252855 1.743075   100  a

【讨论】：

【解决方案4】：

使用所有组合而不用替换：

只为testA组：

comb <- arrangements::combinations(names(df)[grep("testA",names(df))], k = 2,n =  3,replace = F )

tTest <- function(x, data = df){ 
  ttest <- t.test(x =data[x[1]] , y = data[x[2]])
  return(data.frame(var1 = x[1],
                    var2 = x[2],
                    t = ttest[["statistic"]][["t"]],
                    pvalue = ttest[["p.value"]]))
}

result <- apply(comb, 1, tTest, data = df)

结果：

dplyr::bind_rows(result)
         var1        var2          t    pvalue
1 testA_wave1 testA_wave2  0.5009236 0.6193176
2 testA_wave1 testA_wave3 -0.6426433 0.5243146
3 testA_wave2 testA_wave3 -1.1564854 0.2547069

对于所有组：

comb <- arrangements::combinations(x = names(df)[-1], k = 2,n =  6, replace = F )
result <- apply(comb, 1, tTest, data = df)

结果：

dplyr::bind_rows(result)

         var1        var2          t    pvalue
1  testA_wave1 testA_wave2  0.5009236 0.6193176
2  testA_wave1 testA_wave3 -0.6426433 0.5243146
3  testA_wave1 testB_wave1  0.4199215 0.6769510
4  testA_wave1 testB_wave2 -0.3447992 0.7321465
5  testA_wave1 testB_wave3  0.0000000 1.0000000
6  testA_wave2 testA_wave3 -1.1564854 0.2547069
7  testA_wave2 testB_wave1 -0.1070172 0.9153442
8  testA_wave2 testB_wave2 -0.8516264 0.3997630
9  testA_wave2 testB_wave3 -0.5640491 0.5762010
10 testA_wave3 testB_wave1  1.1068781 0.2754186
11 testA_wave3 testB_wave2  0.2966237 0.7683692
12 testA_wave3 testB_wave3  0.7211103 0.4755291
13 testB_wave1 testB_wave2 -0.7874100 0.4360152
14 testB_wave1 testB_wave3 -0.4791735 0.6346043
15 testB_wave2 testB_wave3  0.3865414 0.7013933

【讨论】：

【解决方案5】：

再加入另一个更简洁的data.table 解决方案，我们将数据融合为长格式：

setDT(df)
x = melt(df[,-1])[, tname := sub('_.+','',variable)][, wave := sub('.+_','',variable)]  

x[wave != 'wave1', .(p.value = 
   t.test(x[tname==test & wave == 'wave1', value], value, paired = TRUE)$p.value), 
  by = .(test=tname,wave)]
#     test  wave   p.value
# 1: testA wave2 0.6642769
# 2: testA wave3 0.9209554
# 3: testB wave2 0.1456059
# 4: testB wave3 0.4184603

【讨论】：