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?
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)