使用 perm.t.test 对分组变量进行多重比较

Question

我有一些实验数据要用 R 分析，但我遇到了问题，经过几天的搜索后我找不到解决方案。

我需要对针对不同变量分组的数据进行运行多重排列 t 检验和 Mann-Whitney 检验。

例如，我必须说明在每个实验日 (t) 的处理 (treat) 之间我的响应变量 (exparat) 是否存在差异。

这是我的数据集的样子：

cham spg treat  t exparat
1 T2S2  A6     T T0   1e-04
2 T2S2  A7     T T0   1e-04
3 T2S2  A9     T T0   1e-04
4 T2S2 A10     T T0   1e-04
5 T3S2 A11     T T0   1e-04
6 C1S2 A17     C T0   1e-04

如果我必须使用参数 t 检验，我会使用 dplyr 管道和函数 group_by:

 stat.test <- data %>%
  group_by(t) %>%
  t_test(RespVar ~ treat)

但是我不能对排列 t 检验做同样的事情（我正在使用函数 perm.t.test，perm.t.test {RVAideMemoire}。所以我必须多次编写代码这样：

    dt1 = perm.t.test(exparat~treat,data = subset(data, t == "T0"), nperm=999)
    dt2 = perm.t.test(exparat~treat,data = subset(data, t == "T1"), nperm=999)
    dt3 = perm.t.test(exparat~treat,data = subset(data, t == "T2"), nperm=999)
    dt4 = perm.t.test(exparat~treat,data = subset(data, t == "T3"), nperm=999)
    dt5 = perm.t.test(exparat~treat,data = subset(data, t == "T4"), nperm=999)
    dt6 = perm.t.test(exparat~treat,data = subset(data, t == "T5"), nperm=999)
    dt7 = perm.t.test(exparat~treat,data = subset(data, t == "T6"), nperm=999)
    dt8 = perm.t.test(exparat~treat,data = subset(data, t == "T7"), nperm=999)
    dt9 = perm.t.test(exparat~treat,data = subset(data, t == "T8"), nperm=999)
    dt10 = perm.t.test(exparat~treat,data = subset(data, t == "T9"), nperm=999)
    dt11 = perm.t.test(exparat~treat,data = subset(data, t == "T10"), nperm=999)
    dt12 = perm.t.test(exparat~treat,data = subset(data, t == "T11"), nperm=999)
    dt12 = perm.t.test(exparat~treat,data = subset(data, t == "T12"), nperm=999)
    dt13 = perm.t.test(exparat~treat,data = subset(data, t == "T13"), nperm=999)
    dt14 = perm.t.test(exparat~treat,data = subset(data, t == "T14"), nperm=999)

和提取，并结合每个分析的结果进行 data.frame 导出：

t = dt1$statistic
perm =dt1$permutations
p_val = dt1$p.value
dt1 = data.frame(t, perm, p_val)
row.names(dt1) <- "dt1"

t = dt2$statistic
perm =dt2$permutations
p_val = dt2$p.value
dt2 = data.frame(t, perm, p_val)
row.names(dt2) <- "dt2"

t = dt3$statistic
perm =dt3$permutations
p_val = dt3$p.value
dt3 = data.frame(t, perm, p_val)
row.names(dt3) <- "dt3"

t = dt4$statistic
perm =dt4$permutations
p_val = dt4$p.value
dt4 = data.frame(t, perm, p_val)
row.names(dt4) <- "dt4"

t = dt5$statistic
perm =dt5$permutations
p_val = dt5$p.value
dt5 = data.frame(t, perm, p_val)
row.names(dt5) <- "dt5"

t = dt6$statistic
perm =dt6$permutations
p_val = dt6$p.value
dt6 = data.frame(t, perm, p_val)
row.names(dt6) <- "dt6"

t = dt7$statistic
perm =dt7$permutations
p_val = dt7$p.value
dt7 = data.frame(t, perm, p_val)
row.names(dt7) <- "dt7"

t = dt8$statistic
perm =dt8$permutations
p_val = dt8$p.value
dt8 = data.frame(t, perm, p_val)
row.names(dt8) <- "dt8"

t = dt9$statistic
perm =dt9$permutations
p_val = dt9$p.value
dt9 = data.frame(t, perm, p_val)
row.names(dt9) <- "dt9"

t = dt10$statistic
perm =dt10$permutations
p_val = dt10$p.value
dt10 = data.frame(t, perm, p_val)
row.names(dt10) <- "dt10"

t = dt11$statistic
perm =dt11$permutations
p_val = dt11$p.value
dt11 = data.frame(t, perm, p_val)
row.names(dt11) <- "dt11"

t = dt12$statistic
perm =dt12$permutations
p_val = dt12$p.value
dt12 = data.frame(t, perm, p_val)
row.names(dt12) <- "dt12"

t = dt13$statistic
perm =dt13$permutations
p_val = dt13$p.value
dt13 = data.frame(t, perm, p_val)
row.names(dt13) <- "dt13"

t = dt14$statistic
perm =dt14$permutations
p_val = dt14$p.value
dt14 = data.frame(t, perm, p_val)
row.names(dt14) <- "dt14"

expar2 = rbind(dt1, dt2, dt3, dt4, dt5,dt6,dt7,dt8,dt9,dt10,dt11,dt12,dt13,dt14)
expar2 = data.frame(expar2)

expar2$padj <- p.adjust(expar2$p_val,method="BH")
options(scipen=999)
expar2$padj <- round(expar2$padj,4)
expar2$p.value <- round(expar2$p.value,4)
expar2

这个过程非常费力和耗时，我有一长串变量要分析哪些因素是不同的，所以我为每个变量写了一行又一行的代码。

有没有更简单的方法来完成这一切？

Answer 1

考虑by对数据进行子集处理，处理每个子集，然后do.call+rbind对子集进行堆叠：

dt_list <- by(data, data$t, function(sub) {
   dt <- perm.t.test(exparat~treat, data=sub, nperm=999)
   df <- data.frame(
       t = dt$statistic, 
       perm = dt$permutations,
       p_val = dt$p.value,
       row.names = sub$t[[1]]
   )
   return(df)
})

expar2 <- do.call(rbind.data.frame, dt_list)

expar2 <- within(expar2, {
    padj <- round(p.adjust(p_val, method="BH"), 4)
    p.value <- round($p.value, 4)
})

expar2

Answer 2

在tidyverse中你可以试试下面的-

library(tidyverse)

data %>%
  group_by(t) %>%
  summarise(model = list(perm.t.test(exparat~treat,
                         data = cur_data(), nperm=999))) %>%
  mutate(res = map(model, broom::tidy)) %>%
  unnest(res)

使用 perm.t.test 对分组变量进行多重比较

Multiple comparisons using perm.t.test on grouped variables

r

loops

permutation

t-test