跟进 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
我想做的是在标记的基础上为每个Patient组合(我的分组变量)运行t.test (我的 ID 变量)。
根据对上述相关问题的一个回答，我知道如何一次对一个标记执行此操作。
我可以子集 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 中的所有 Markers 做这件事，我选择了 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 相同(相同的数据子集仍在用于两者)...
所以现在我一无所知......我在这里做错了什么？应该怎么做呢？
谢谢!
最佳答案
这是一个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")
代码
我在代码的注释中添加了一些简短的解释和中间/临时结果。但它已经变得比代码更多的注释 ;-)...
不管怎样，我们开始吧……
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
关于R:对分组变量的每个成对组合进行 t 检验，对 ID 变量中的每个元素进行，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/60150655/