R：对分组变量的每个成对组合进行 t 检验，对 ID 变量中的每个元素进行答案

【问题标题】：R: t-test for each pairwise combination of a grouping variable, done for every element in an ID variableR：对分组变量的每个成对组合进行 t 检验，对 ID 变量中的每个元素进行
【发布时间】：2020-02-10 12:40:42
【问题描述】：

跟进this question，我正在尝试增加一层难度。

我有一个看起来像这样的data.frame：

> set.seed(123)
> mydf <- data.frame(Marker=rep(c('M1','M2'),each=15),
+                    Patient=rep(rep(c('P1','P2','P3'),each=5),2),
+                    Value=sample(1:1000, 30, replace = F))
> mydf
   Marker Patient Value
1      M1      P1   288
2      M1      P1   788
3      M1      P1   409
4      M1      P1   881
5      M1      P1   937
6      M1      P2    46
7      M1      P2   525
8      M1      P2   887
9      M1      P2   548
10     M1      P2   453
11     M1      P3   948
12     M1      P3   449
13     M1      P3   670
14     M1      P3   566
15     M1      P3   102
16     M2      P1   993
17     M2      P1   243
18     M2      P1    42
19     M2      P1   323
20     M2      P1   996
21     M2      P2   872
22     M2      P2   679
23     M2      P2   627
24     M2      P2   972
25     M2      P2   640
26     M2      P3   691
27     M2      P3   530
28     M2      P3   579
29     M2      P3   282
30     M2      P3   143

我想要做的是，在 Marker 基础上（我的 ID 变量）为每个 Patient 组合（我的分组变量）运行 t.test。

根据对上述相关问题的一个答案，我知道如何一次为一个 Marker 执行此操作。

我可以对mydf 进行子集化并执行以下操作：

> params_list <- utils::combn(levels(mydf$Patient), 2, FUN = list)
> mydf0 <- subset(mydf, Marker=="M1")
> model_t <- purrr::map(.x = params_list, 
+                       .f = ~ t.test(formula = Value ~ Patient, 
+                       data = subset(mydf0, Patient %in% .x)))
> t_pvals <- purrr::map_dbl(.x = model_t, .f  = "p.value")
> names(t_pvals) <- purrr::map_chr(.x = params_list, .f = ~ paste0(.x, collapse = "-vs-"))
> t_pvals
 P1-vs-P2  P1-vs-P3  P2-vs-P3 
0.3945742 0.5678729 0.7820905

现在我想以一种优雅的方式为mydf 中的所有标记做这件事，我选择了data.table。

我尝试以下操作，但无法重现 Marker M1 的上述 pvalue 结果。

> group1 <- unlist(lapply(params_list, '[', 1))
> group2 <- unlist(lapply(params_list, '[', 2))
> mydt <- data.table::data.table(mydf)
> results_df <- as.data.frame(mydt[, list(group1= unlist(lapply(params_list, '[', 1)),
+                                         group2= unlist(lapply(params_list, '[', 2)),
+                                         pvalue= purrr::map_dbl(.x = purrr::map(.x = params_list,
+                                                 .f = ~ stats::t.test(formula = Value ~ Patient, paired=FALSE,
+                                                 data = subset(mydt, Patient %in% .x))), .f  = "p.value") ),
+                                  by=list(Marker=Marker)])
> results_df
  Marker group1 group2    pvalue
1     M1     P1     P2 0.8092365
2     M1     P1     P3 0.5156313
3     M1     P2     P3 0.2879954
4     M2     P1     P2 0.8092365
5     M2     P1     P3 0.5156313
6     M2     P2     P3 0.2879954

results_df 的结构正是我想要的，但是 pvalues 显然是错误的。 M1与上面测试中的不同，M1和M2相同，表示相同的数据子集在这两种情况下都使用。

我认为我应该在 subset 命令中为每个 Marker 设置子集，所以我这样做了：

> markers_list <- as.list(levels(mydf$Marker))
> mydt <- data.table::data.table(mydf)
> results_df <- as.data.frame(mydt[, list(group1= unlist(lapply(params_list, '[', 1)),
+                                         group2= unlist(lapply(params_list, '[', 2)),
+                                         pvalue= purrr::map_dbl(.x = purrr::map(.x = params_list, .y = markers_list,
+                                                 .f = ~ stats::t.test(formula = Value ~ Patient, paired=FALSE,
+                                                 data = subset(mydt, Patient %in% .x & Marker==.y))), .f  = "p.value") ),
+                                  by=list(Marker=Marker)])
> results_df
  Marker group1 group2    pvalue
1     M1     P1     P2 0.7337355
2     M1     P1     P3 0.6930669
3     M1     P2     P3 0.3788015
4     M2     P1     P2 0.7337355
5     M2     P1     P3 0.6930669
6     M2     P2     P3 0.3788015

我以为就是这样，但我仍然得到不正确的 pvalues，并且 M1 和 M2 相同（两者仍使用相同的数据子集）...

所以现在我一无所知......我在这里做错了什么？该怎么做？

谢谢！

【问题讨论】：

标签： r dataframe testing data.table purrr

【解决方案1】：

这是一个data.table-解决方案

我无法重现您的示例数据，因此我读取了使用 data.table::fread() 提供的值。

您还可以在现有的 mydf 上使用 data.table::setDT(mydf) 将其转换为 data.table。

样本数据

library(data.table)
#setDT(mydf)   
mydf <- fread("   Marker Patient Value
      M1      P1   288
      M1      P1   788
      M1      P1   409
      M1      P1   881
      M1      P1   937
      M1      P2    46
      M1      P2   525
      M1      P2   887
      M1      P2   548
     M1      P2   453
     M1      P3   948
     M1      P3   449
     M1      P3   670
     M1      P3   566
     M1      P3   102
     M2      P1   993
     M2      P1   243
     M2      P1    42
     M2      P1   323
     M2      P1   996
     M2      P2   872
     M2      P2   679
     M2      P2   627
     M2      P2   972
     M2      P2   640
     M2      P3   691
     M2      P3   530
     M2      P3   579
     M2      P3   282
   M2      P3   143")

代码

我在代码中添加了一些简短的解释，并在 cmets 中添加了中间/临时结果。但它变得比代码更多的是注释；-)...
不管怎样，我们走吧……

mydf[, 
     #suppress immediate output using {}
     {
     # find all unique combinations of 2 patients (by Marker, see last line)
     # For Marker == "M1", this looks like:
      #    V1 V2
      # 1: P1 P2
      # 2: P1 P3
      # 3: P2 P3
     patientcomb <- data.table( t( combn( unique( Patient ), 2 ) ) )
     #set column names for V1 and V2 of patientcomb, for better readable code
     names( patientcomb ) <- c( "group1", "group2" )
     #now, using the temporarily created patientcomb-data.table...
     patientcomb[,
                 #... perform the t.test(), using the Values from mydf, 
                 #  where the patients match group1/group1
                 #remember, we are still grouped by Marker
                 data.table( p.value = t.test( Value[Patient == group1], 
                                               Value[Patient == group2])$p.value), 
                 #group by group1 and group2
                 by = .(group1, group2) ]
     # for Marker == M1, this looks like:
      #    group1 group2   p.value
      # 1:     P1     P2 0.3945742
      # 2:     P1     P3 0.5678729
      # 3:     P2     P3 0.7820905
     # for Marker == M2, this looks like:
      #    group1 group2   p.value
      # 1:     P1     P2 0.3098955
      # 2:     P1     P3 0.7505371
      # 3:     P2     P3 0.0372944
     }, 
    #main grouping by Marker
    by = .(Marker) ]

输出

似乎匹配所需的输出

#    Marker group1 group2   p.value
# 1:     M1     P1     P2 0.3945742
# 2:     M1     P1     P3 0.5678729
# 3:     M1     P2     P3 0.7820905
# 4:     M2     P1     P2 0.3098955
# 5:     M2     P1     P3 0.7505371
# 6:     M2     P2     P3 0.0372944

【讨论】：

【解决方案2】：

data.table 中的另一个选项：

mydf[, rbindlist(combn(split(Value, Patient), 2L, 
        function(x) c(as.list(names(x)), .(t.test(x[[1]], x[[2]])$p.value)), simplify=FALSE))
    , Marker]

输出：

   Marker V1 V2        V3
1:     M1 P1 P2 0.3945742
2:     M1 P1 P3 0.5678729
3:     M1 P2 P3 0.7820905
4:     M2 P1 P2 0.3098955
5:     M2 P1 P3 0.7505371
6:     M2 P2 P3 0.0372944

数据：

library(data.table)
mydf <- fread("
Marker Patient Value
M1      P1   288
M1      P1   788
M1      P1   409
M1      P1   881
M1      P1   937
M1      P2    46
M1      P2   525
M1      P2   887
M1      P2   548
M1      P2   453
M1      P3   948
M1      P3   449
M1      P3   670
M1      P3   566
M1      P3   102
M2      P1   993
M2      P1   243
M2      P1    42
M2      P1   323
M2      P1   996
M2      P2   872
M2      P2   679
M2      P2   627
M2      P2   972
M2      P2   640
M2      P3   691
M2      P3   530
M2      P3   579
M2      P3   282
M2      P3   143")

【讨论】：

我更喜欢这段代码的简单性，因为我尝试这样做的方式更像是......只是几个问题：1）2L 在这里是什么意思？与2 有什么区别？ - 2) 我注意到如果我这样做 group=as.list(names(x))... 我会得到 group1 和 group2 而不是 V1 和 V2 我想要...但是我将如何更改后缀？如果我想说group_A 和group_B 怎么办？非常感谢！
是的，2L 是整数 2。对于第二个 qn，你 acn 使用 setNames(as.list(names(x)), c("group_A","group_B"))

【解决方案3】：

这是一种tidyverse 方法：

library(tidyverse)

get_p_value <- function(df) {
   map_df(params_list,  ~{
     tibble(Marker = df[[1]][1], group1 = .x[1], group2 = .x[2], 
       pvalue =  t.test(df$Value[df$Patient == .x[1]], 
                        df$Value[df$Patient == .x[2]])$p.value)
      })
}

mydf %>% group_split(Marker) %>% map_df(get_p_value)

# A tibble: 6 x 4
#  Marker group1 group2 pvalue
#  <fct>  <chr>  <chr>   <dbl>
#1 M1     P1     P2     0.395 
#2 M1     P1     P3     0.568 
#3 M1     P2     P3     0.782 
#4 M2     P1     P2     0.310 
#5 M2     P1     P3     0.751 
#6 M2     P2     P3     0.0373

params_list 来自 OP。

params_list <- utils::combn(levels(mydf$Patient), 2, FUN = list)

【讨论】：

【解决方案4】：

对由Marker 分组的数据使用pairwise.t.test() 似乎是解决此问题的更好方法，并且无需显式生成Patient 组合。

library(dplyr)
library(tidyr)

mydf %>%
  group_by(Marker) %>%
  summarise(x = list(pairwise.t.test(Value, Patient, p.adjust.method = "none", pool.sd = FALSE)$p.value %>% as.data.frame.table(responseName = "p.value"))) %>%
  unnest(x) %>%
  filter(!is.na(p.value))

# A tibble: 6 x 4
  Marker Var1  Var2  p.value
  <fct>  <fct> <fct>   <dbl>
1 M1     P2    P1     0.395 
2 M1     P3    P1     0.568 
3 M1     P3    P2     0.782 
4 M2     P2    P1     0.310 
5 M2     P3    P1     0.751 
6 M2     P3    P2     0.0373

针对您的评论，还有一个成对版本的 wilcox 测试：

mydf %>%
  group_by(Marker) %>%
  summarise(x = list(pairwise.wilcox.test(Value, Patient, p.adjust.method = "none")$p.value %>% as.data.frame.table(responseName = "p.value"))) %>%
  unnest(x) %>%
  filter(!is.na(p.value))

# A tibble: 6 x 4
  Marker Var1  Var2  p.value
  <fct>  <fct> <fct>   <dbl>
1 M1     P2    P1     0.690 
2 M1     P3    P1     0.841 
3 M1     P3    P2     0.690 
4 M2     P2    P1     0.690 
5 M2     P3    P1     1     
6 M2     P3    P2     0.0556

【讨论】：

好吧，在我的现实生活数据中，大多数情况下它实际上是 wilcox.test（我只是在 MWE 中使用了 t.test 以便更容易理解）......但实际上是否存在非参数等效于pairwise.t.test()?
@DaniCee- wilcox 测试的成对版本确实有一个功能 - 它的工作方式相同。我在一个示例中进行了编辑。
我有一个关于pairwise.t.test 的一般性问题，看看你能不能帮助我... 我通常使用 Levene 的测试来测试等方差；此外，有时数据来自相同的患者，所以我认为它是配对的。因此，如果每个标记和分组变量级别的方差相等，我将运行 t.test 和 var.equal=TRUE 和 paired=TRUE（正确吗？）......但是，pairwise.t.test 不允许同时使用 @987654334 @ 和 paired 为 TRUE
@DaniCee - 不，这是不正确的。配对 t 检验等效于一个样本 t 检验，真实均值为零。在t.test(..., paired = TRUE) 中，同时设置var.equal = TRUE 会导致这个参数被默默地忽略，因为单一的方差源不能被合并。无论var.equal 是TRUE 还是FALSE，当paired = TRUE 时，您都会看到完全相同的结果。我认为pairwise.t.test() 可以通过防止这些论点同时为真来更好地处理这个问题。
哦哦哦好吧！再有一个理由来切换。非常感谢！