Pruebas post hoc en multcomp :: glht para modelos de efectos mixtos (lme4) con interacciones

Estoy realizando pruebas post-hoc en un modelo lineal de efectos mixtos en R( lme4paquete). Estoy usando el multcomppaquete ( glht()función) para realizar las pruebas post-hoc.

Mi diseño experimental es de medidas repetidas, con un efecto de bloqueo aleatorio. Los modelos se especifican como:

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

En lugar de adjuntar mis datos aquí, estoy trabajando con los datos llamados warpbreaksdentro del multcomppaquete.

data <- warpbreaks
warpbreaks$rand <- NA

He agregado una variable aleatoria adicional para imitar mi efecto de "bloqueo":

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

Esto imita mi modelo:

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

Soy consciente del ejemplo en los " Ejemplos adicionales de Multcomp- 2 Way Anova" Este ejemplo lo lleva a una comparación de niveles de tensión dentro de los niveles de wool.

¿Qué sucede si quiero hacer lo contrario: comparar los niveles wooldentro de los niveles de tension? (En mi caso, esto sería comparar los niveles de tratamiento (dos - 0, 1) dentro de los niveles de tiempo (tres - junio, julio, agosto).

Se me ocurrió el siguiente código para hacerlo, pero parece que no funciona (vea el mensaje de error a continuación).

Primero, del ejemplo (con wooly tensionlugares intercambiados):

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

De aquí a abajo, mi propio 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

Es mucho más fácil hacerlo usando el paquete lsmeans

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