Как выполнить Tukey HSD.test() в списке данных?

Question

Я хотел бы выполнить тест Tukey Post Hoc для списка данных. В качестве результата я хотел бы получить буквы, указывающие, какие группы значительно отличаются друг от друга. HSD.test() из пакета Agricola делает это, но я не могу понять, как применить это к нескольким переменным одновременно. Это сэкономит мне много времени, так как мой набор данных содержит много переменных.

Это часть моих данных:

category <- c(rep("young", 3), rep("Middle", 4), rep("old", 5))
fat <- c(1857.87, 1953.90, 1440.70, 1553.81, 1785.91, 1893.82, 1483.75, 1784.99, 2011.01, 2023.04, 2011.05, 1788.81)
BMI <- c(21.1, 23.2, 24.5, 25.6, 21.8, 18.0, 19.2, 20.1, 22.1, 25.0, 26.1, 25.1)
age <- c(25, 23, 27, 55, 58, 62, 45, 75, 80, 75, 83, 89)
df2 <- data.frame(fat, BMI, age, category)

Я знаю, как сделать anova и HSD.test() для одной переменной.

lm.fat <- (lm(fat ~ as.factor(category), data = df2))
anova(lm.fat)
require(agricolae)
HSD.test(lm.fat, "as.factor(category)", group = TRUE, console = TRUE)

Кроме того, я знаю, как применить anova ко всем моим переменным в наборе данных с помощью функции sapply():

an <- lapply(df, function(x) aov(x~category, data = df))
sapply(an, anova, simplify=FALSE)

Но я не знаю, как выполнить posthoc HSD.test() для этих результатов. Я попробовал это:

lapply(an, function(m) HSD.test((m), "as.factor(category)", group = TRUE, console = TRUE))

Name:  as.factor(category) 
 category 
Name:  as.factor(category) 
 category 
Name:  as.factor(category) 
 category 

$fat
NULL

$BMI
NULL

$age
NULL

Я получаю NULL в качестве результата, поэтому что-то пошло не так, но я не могу понять, что. Я попробовал также другую функцию теста Tukey post hoc, а именно: TukeyHSD().

lapply(an, function(m) TukeyHSD(aov(m)))
$fat
Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = m)

$category
                   diff       lwr      upr     p adj
old-Middle    244.45750 -110.2508 599.1658 0.1874162
young-Middle   71.50083 -332.3523 475.3540 0.8757638
young-old    -172.95667 -559.1142 213.2008 0.4554780


$BMI
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = m)

$category
                   diff       lwr      upr     p adj
old-Middle    2.5300000 -2.494240 7.554240 0.3781804
young-Middle  1.7833333 -3.937015 7.503682 0.6711525
young-old    -0.7466667 -6.216366 4.723033 0.9237118


$age
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = m)

$category
              diff       lwr       upr     p adj
old-Middle    25.4  14.49330  36.30670 0.0002928
young-Middle -30.0 -42.41783 -17.58217 0.0002219
young-old    -55.4 -67.27372 -43.52628 0.0000010

Это работает для моего набора данных, но эта функция не дает букв, которые я действительно хотел бы получить. Кто-нибудь знает, как я могу сделать то же самое с HSD.test(), чтобы я получил буквы? Спасибо!

Answer 1

Вы можете попробовать следующий код:

List <- names(df2)[1:3] # select just the variables

model1 <- lapply(List, function(x) {
      lm(substitute(i~category, list(i = as.name(x))), data = df2)})


lapply(model1, summary)

letters = lapply(model1, function(m) HSD.test((m), "category", group = TRUE, console = TRUE))

Если бы это было взаимодействие, это было бы:

tx <- with(df2, interaction(category1, category2))  # determining the factors

model2 <- lapply(List, function(x) {
  glm(substitute(i~tx, list(i = as.name(x))), data = df2)}) # using the factors already in "tx"

lapply(model2, summary)

letters = lapply(model2, function(m) HSD.test((m), "tx", alpha = 0.05, group = TRUE, console = TRUE))

С наилучшими пожеланиями