R: тест Крускала-Уоллиса в цикле по указанным столбцам в кадре данных

Question

Я хотел бы запустить KW-тест для определенных числовых переменных из фрейма данных, используя одну группирующую переменную. Я бы предпочел сделать это в цикле, вместо того, чтобы печатать все тесты, так как в них много переменных (больше, чем в примере ниже).

Имитация данных:

library(dplyr)
set.seed(123)
Data <- tbl_df(
data.frame(
muttype = as.factor(rep(c("missense", "frameshift", "nonsense"), each = 80)),
ados.tsc   = runif(240, 0, 10),
ados.sa    = runif(240, 0, 10),
ados.rrb   = runif(240, 0, 10))
) %>%
group_by(muttype)
ados.sim <- as.data.frame(Data)

Следующий код прекрасно работает вне цикла.

kruskal.test(formula (paste((colnames(ados.sim)[2]), "~ muttype")), data = 
ados.sim)

Но это не внутри цикла:

for(i in names(ados.sim[,2:4])){  
ados.mtp <- kruskal.test(formula (paste((colnames(ados.sim)[i]), "~ muttype")), 
data = ados.sim)
}

Я получаю ошибку:

Ошибка в термине. Формула (формула, данные = данные): недопустимый термин в формуле модели

Кто-нибудь знает, как это решить? Очень ценится!

Answer 1

Попробуйте:

results <- list()
for(i in names(ados.sim[,2:4])){  
  results[[i]] <- kruskal.test(formula(paste(i, "~ muttype")), data = ados.sim)
}

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

Обратите внимание на следующее:

for(i in names(ados.sim[,2:4])){  
   print(i)
}
[1] "ados.tsc"
[1] "ados.sa"
[1] "ados.rrb"

То есть i уже дает вам имя столбца. Проблема в вашем коде заключалась в том, что вы пытались использовать его как целое число для поднабора, что превратило результат в NA.

for(i in names(ados.sim[,2:4])){  
   print(paste((colnames(ados.sim)[i]), "~ muttype"))
}
[1] "NA ~ muttype"
[1] "NA ~ muttype"
[1] "NA ~ muttype"

И просто для справки, все это также может быть сделано следующими двумя способами, которые я часто предпочитаю, поскольку это немного облегчает последующий анализ:

Сначала сохраните все тестовые объекты в кадре данных:

library(tidyr)
df <- ados.sim %>% gather(key, value, -muttype) %>% 
      group_by(key) %>% 
      do(test = kruskal.test(x= .$value, g = .$muttype))

Затем вы можете установить подкадр данных для получения результатов теста:

df[df$key == "ados.rrb",]$test
[[1]]

    Kruskal-Wallis rank sum test

data:  .$value and .$muttype
Kruskal-Wallis chi-squared = 2.2205, df = 2, p-value = 0.3295

Либо получить все результаты непосредственно в кадре данных, не сохраняя тестовые объекты:

library(broom)
df2 <- ados.sim %>% gather(key, value, -muttype) %>% 
       group_by(key) %>% 
       do(tidy(kruskal.test(x= .$value, g = .$muttype)))
df2
# A tibble: 3 x 5
# Groups:   key [3]
       key statistic   p.value parameter                       method
     <chr>     <dbl>     <dbl>     <int>                       <fctr>
1 ados.rrb 2.2205031 0.3294761         2 Kruskal-Wallis rank sum test
2  ados.sa 0.1319554 0.9361517         2 Kruskal-Wallis rank sum test
3 ados.tsc 0.3618102 0.8345146         2 Kruskal-Wallis rank sum test