Tests post-hoc en multcomp :: glht pour les modèles à effets mixtes (lme4) avec interactions

J'effectue des tests post-hoc sur un modèle linéaire à effets mixtes dans R( lme4package). J'utilise multcomppackage ( glht()fonction) pour effectuer les tests post-hoc.

Mon plan expérimental est des mesures répétées, avec un effet de bloc aléatoire. Les modèles sont spécifiés comme:

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

Plutôt que de joindre mes données ici, je travaille à partir des données appelées warpbreaksdans le multcomppackage.

data <- warpbreaks
warpbreaks$rand <- NA

J'ai ajouté une variable aléatoire supplémentaire pour imiter mon effet "bloc":

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

Cela imite mon modèle:

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

Je connais l'exemple dans les " Exemples Multcomp supplémentaires - Anova à 2 voies" Cet exemple vous amène à une comparaison des niveaux de tension au sein des niveaux de wool.

Et si je veux faire le contraire - comparer les niveaux de wooldans les niveaux de tension? (Dans mon cas, ce serait comparer les niveaux de traitement (deux - 0, 1) dans les niveaux de temps (trois - juin, juillet, août).

J'ai trouvé le code suivant pour le faire, mais cela ne semble pas fonctionner (voir le message d'erreur ci-dessous).

Tout d'abord, à partir de l'exemple (avec woolet tensionlieux échangés):

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 = ":")

D'ici en bas, mon propre code:

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

C'est beaucoup plus facile à faire en utilisant le paquet lsmeans

library(lsmeans)
lsmeans(mod, pairwise ~ tension | wool)
lsmeans(mod, pairwise ~ wool | tension)