Множественные сравнения с использованием perm.t.test для сгруппированных переменных

У меня есть данные эксперимента для анализа с помощью R, но у меня возникла проблема, и после нескольких дней поиска я не могу найти решение.

Мне нужно запустить множественные t-тесты перестановок и тесты Манна-Уитни для моих данных, сгруппированных для разных переменных.

Например, я должен сказать, есть ли различия в моей переменной ответа (exparat) между обработками (обработкой) в каждый экспериментальный день (t).

Вот как выглядит мой набор данных:

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

Этот процесс очень трудоемкий и требует много времени, и у меня есть длинный список переменных для анализа, для которых факторы различаются, и поэтому я должен писать строки и строки кода для каждой из них.

Есть ли способ сделать все это попроще?


person Valerio M    schedule 26.05.2021    source источник


Ответы (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)
person Ronak Shah    schedule 26.05.2021
comment
Тоже работает, очень полезно. - person Valerio M; 28.05.2021
comment
Большое вам спасибо за вашу помощь (могу я это написать?). Как ты думаешь, сможешь ли ты мне помочь и с этой другой похожей проблемой? stackoverflow .com / questions / 67735153 / - person Valerio M; 29.05.2021

Рассмотрим 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
person Parfait    schedule 26.05.2021
comment
Работает, очень полезно. - person Valerio M; 28.05.2021