1

I have some data from an experiment to analyse with R but I have a problem and after days of search I can't find a solution.

I need to run multiple permutation t-tests and Mann-Whitney tests on my data grouped for different variables.

For examples, I have to say if there are differences in my response variable (exparat) between treatments (treat) on each experimental day (t).

This is how my dataset looks like:

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

If I had to use a parametric t-test I would have used a dplyr pipe and the function group_by:

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

But I can't do the same thing for permutation t-tests (I'm using the function perm.t.test, perm.t.test {RVAideMemoire}. So I have to write the code several times in this way:

    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)

and the extract, and combine the results for each of the analyses to make a data.frame to export:

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

This process is very laborious and time consuming, and I have a long list of variables to analyse for which the factors are different and so I have write lines and lines of codes for each of them.

Is there a simpler way to do all this?

Behnam Hedayat
  • 837
  • 4
  • 18
Valen78
  • 115
  • 6

2 Answers2

0

Consider by for subsetting data, processing each subset, and then do.call + rbind to stack subsets:

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
Parfait
  • 104,375
  • 17
  • 94
  • 125
0

In tidyverse you can try the following -

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)
Ronak Shah
  • 377,200
  • 20
  • 156
  • 213
  • Thank you so much for all your help (can I write it?). Do you think you can help me also with this other similar problem? https://stackoverflow.com/questions/67735153/multiple-paired-permutation-t-tests-using-perm-t-test – Valen78 May 29 '21 at 01:01