Как получить компактный буквенный дисплей Тьюки из GLM с взаимодействиями

У меня есть набор данных, которые я проанализировал с помощью обобщенной линейной модели, которая имеет три категориальных фактора в трехстороннем взаимодействии (factorA, factorB, factorC) и четвертый непрерывный фактор (factorD), который просто добавляется в модель. Я пытаюсь получить набор групп букв Тьюки (т. Е. Компактное отображение букв) из модели, но не нашел способа успешно включить взаимодействие. Меня не интересует включение factorD, только трое во взаимодействии.

Я получил попарные сравнения с поправкой на Тьюки:

lsmeans(my.glm, factorA*factorB*factorC)

Но я не смог понять, как из этого сделать компактный дисплей букв. Это можно сделать с multcomp пакетом, но я мог найти способы сделать это только с основными эффектами этого пакета, а не с взаимодействиями.

Тогда я попробовал пакет agricolae, так как этот пост (https://stats.stackexchange.com/questions/31547/how-to-obtain-the-results-of-a-tukey-hsd-post-hoc-test-in-a-table-pting-groupe) обсуждает, что это должно сработать. Однако следование инструкциям в этом ответе привело к нефункциональному ответу HSD.test. В частности, я мог бы заставить работать тесты основных эффектов, например HSD.test(my.glm,"factorA") но я не мог заставить взаимодействие работать. Я пробовал это:

intxns<-with(my.data, interaction(factorA,factorB,factorC))
HSD.test(my.glm,"intxns",group=TRUE)

Но получить сообщение об ошибке, которое указывает, что функция HSD.test не распознала intxns как допустимый объект, это выглядит так (я также проверил объект intxns, он выглядит хорошо, и количество строк соответствует количеству остатков моего glm):

Name: inxtns
 factorA factorB factorC factorD

Я получаю ту же ошибку, если просто помещаю ерунду в поле фактора при вызове функции HSD.test. Я проверил объект inxtns, он выглядит хорошо, и количество строк соответствует количеству остатков. Примечания agricolae на самом деле не охватывают использование взаимодействий в HSD.test, но я предполагаю, что это может сработать.

Кто-нибудь знает, как заставить HSD.test работать с взаимодействиями? Или есть какая-то другая функция, с которой вам нужно работать для создания компактных буквенных дисплеев для glm с взаимодействиями?

Я работаю над этим уже несколько дней и не могу найти решение, надеюсь, я не упускаю что-то очевидное.

Спасибо!


r interaction glm
person DirtStats    schedule 31.07.2013    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Я не знаю, как вы указали свою модель glm, но для HSD.test она пытается сопоставить конкретное имя лечения с тем же именем, указанным в формуле glm, а также во фрейме данных. Вот почему будет работать ваш главный эффект factorA, но не трехстороннее взаимодействие. Я считаю, что для множественных сравнительных тестов взаимодействий проще всего генерировать взаимодействия отдельно и добавлять их во фрейм данных в качестве дополнительных столбцов. Затем модель glm может быть указана с использованием новых переменных, которые кодируют взаимодействие.

Например,

set.seed(42)
glm.dat <- data.frame(y = rnorm(1000), factorA = sample(letters[1:2], 
   size = 1000, replace = TRUE),
 factorB = sample(letters[1:2], size = 1000, replace = TRUE),
 factorC = sample(letters[1:2], size = 1000, replace = TRUE))

# Generate interactions explicitly and add them to the data.frame
glm.dat$factorAB <- with(glm.dat, interaction(factorA, factorB))
glm.dat$factorAC <- with(glm.dat, interaction(factorA, factorC))
glm.dat$factorBC <- with(glm.dat, interaction(factorB, factorC))
glm.dat$factorABC <- with(glm.dat, interaction(factorA, factorB, factorC))

# General linear model 
 glm.mod <- glm(y ~ factorA + factorB + factorC +  factorAB + factorAC + 
   factorBC + factorABC, family = 'gaussian', data = glm.dat) 

# Multiple comparison test

library(agricolae)
comp <- HSD.test(glm.mod, trt = "factorABC", group = TRUE)

давая

comp$groups giving

    trt        means M
 1 a.a.a  0.070052189 a
 2 a.b.b  0.035684571 a
 3 b.a.a  0.020517535 a
 4 b.b.b -0.008153257 a
 5 a.b.a -0.036136140 a
 6 a.a.b -0.078891136 a
 7 b.a.b -0.080845419 a
 8 b.b.a -0.115808772 a
person Community    schedule 01.08.2013
comment
Большое спасибо, Джим! Кажется, это работает для моей модели и данных. Однако теперь lsmeans не нравится модель - он зависает, пережевывая память и вычислительную мощность всякий раз, когда я пытаюсь запустить модель через нее, даже когда я не использую новые условия взаимодействия. И если я все же использую новые термины взаимодействия, это дает в результате кучу НС. Вы сталкивались с этой проблемой раньше? - person DirtStats; 02.08.2013
comment
@DirtStats - Я столкнулся с этой проблемой и с моими образцами данных, но я не знал, была ли это проблема с моими данными или была какая-то основная проблема с моим решением. Возможно, стоит пометить модераторов, чтобы переместить ваш вопрос в stats.stackexchange. - person ; 02.08.2013