Teste post-hoc no multcomp :: glht para modelos de efeitos mistos (lme4) com interações

10

Estou realizando testes post-hoc em um modelo linear de efeitos mistos em R( lme4pacote). Estou usando o multcomppacote ( glht()função) para executar os testes post-hoc.

Meu projeto experimental é de medidas repetidas, com um efeito de bloqueio aleatório. Os modelos são especificados como:

mymod <- lmer(variable ~ treatment * time + (1|block), data = mydata, REML = TRUE)

Em vez de anexar meus dados aqui, estou trabalhando com os dados chamados warpbreaksdentro do multcomppacote.

data <- warpbreaks
warpbreaks$rand <- NA

Eu adicionei uma variável aleatória extra para imitar meu efeito de "bloqueio":

warpbreaks$rand <- rep(c("foo", "bar", "bee"), nrow(warpbreaks)/3)

Isso imita meu modelo:

mod <- lmer(breaks ~ tension * wool + (1|rand), data = warpbreaks)

Estou ciente do exemplo em " Exemplos adicionais de multcomp - 2 vias Anova" Este exemplo leva você a uma comparação de níveis de tensão dentro dos níveis de wool.

E se eu quiser fazer o oposto - comparar os níveis de wooldentro dos níveis de tension? (No meu caso, isso seria comparar os níveis de tratamento (dois - 0, 1) dentro dos níveis de tempo (três - junho, julho, agosto).

Eu vim com o código a seguir, mas ele não parece funcionar (veja a mensagem de erro abaixo).

Primeiro, do exemplo (com woole tensionlugares trocados):

tmp <- expand.grid(wool = unique(warpbreaks$wool), tension = unique(warpbreaks$tension))
X <- model.matrix(~ tension * wool, data = tmp)
glht(mod, linfct = X)

Tukey <- contrMat(table(warpbreaks$wool), "Tukey")

K1 <- cbind(Tukey, matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)))
rownames(K1) <- paste(levels(warpbreaks$tension)[1], rownames(K1), sep = ":")

K2 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(warpbreaks$tension)[2], rownames(K2), sep = ":")

Daqui para baixo, meu próprio código:

K3 <- cbind(matrix(0, nrow = nrow(Tukey), ncol = ncol(Tukey)), Tukey)
rownames(K2) <- paste(levels(warpbreaks$tension)[3], rownames(K3), sep = ":")

K <- rbind(K1, K2, K3)
colnames(K) <- c(colnames(Tukey), colnames(Tukey))

> summary(glht(mod, linfct = K %*% X))
Error in summary(glht(mod, linfct = K %*% X)) : 
  error in evaluating the argument 'object' in selecting a method for function 'summary': Error in K %*% X : non-conformable arguments
mixed-model post-hoc lsmeans lme4-nlme Ashley Asmus
fonte

Respostas:

6

É muito mais fácil usar o pacote lsmeans

library(lsmeans)
lsmeans(mod, pairwise ~ tension | wool)
lsmeans(mod, pairwise ~ wool | tension)
Russ Lenth
fonte
Ótimo, funciona! Obrigado. Nota: este código só funcionou para os meus dados depois de alterar minha variável de repetição de valores numéricos (3 e 6) para valores alfabéticos (A e B).
Bem, isso importa muito! Porque é um modelo diferente com timeum preditor numérico. Eu suspeito que você queria isso como um fator.
Russ Lenth
Como posso generalizar para mais preditores? se, por exemplo, tenho três preditores, como isso funciona?
divirta-se
11
@ havefun Por favor, olhe help("lsmeans", package = "lsmeans")e vignette("using-lsmeans"). Há muita documentação e muitos exemplos.
Russ Lenth
11
Conte o número de comparações que você obtém com cada método, elas não são as mesmas. Leia também sobre os ajustes de vários testes. Quando você tem uma família maior de testes, os valores de P ajustados são diferentes dos de uma família menor. Quando você usa uma variável by, os ajustes são aplicados separadamente a cada conjunto.
precisa saber é o seguinte