Несколько парных t-тестов для нескольких переменных одновременно с использованием dplyr / tidyverse

Question

Предположим структуру данных, подобную этой:

   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_wave1 против testA_wave2 и testA_wave1 против testA_wave3 в случае testA).

Таким образом, я смог достичь этого:

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))

Answer 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

Здесь он будет выполнять t-тест для каждого столбца "wave .." со столбцом "wave1".

---

Поскольку вы также открыты для других решений, здесь сделана попытка использовать чисто базовое решение 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 к различным комбинациям столбцов для каждого фрейма данных.

Answer 2

Вот один из способов сделать это, совсем немного используя 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

^{Создано в 2019-03-08 пакетом Представить (v0.2.1)}

Answer 3

Использование всех комбинаций без замены:

Только для 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

Answer 4

Чтобы добавить 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

---

Benchmark

Сравнивая решение data.table с решением tidyverse, мы получаем трехкратное увеличение скорости по сравнению с решением data.table:

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

Answer 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