Przeprowadzam testy post-hoc na liniowym modelu mieszanych efektów w R
( lme4
pakiecie). Używam multcomp
pakietu ( glht()
funkcji) do przeprowadzania testów post-hoc.
Mój plan eksperymentalny to powtarzane pomiary z przypadkowym efektem blokowania. Modele są określone jako:
mymod <- lmer(variable ~ treatment * time + (1|block), data = mydata, REML = TRUE)
Zamiast załączać tutaj moje dane, pracuję nad danymi wywoływanymi warpbreaks
w multcomp
pakiecie.
data <- warpbreaks
warpbreaks$rand <- NA
Dodałem dodatkową zmienną losową, aby naśladować mój efekt „bloku”:
warpbreaks$rand <- rep(c("foo", "bar", "bee"), nrow(warpbreaks)/3)
To naśladuje mój model:
mod <- lmer(breaks ~ tension * wool + (1|rand), data = warpbreaks)
Znam ten przykład z „ Dodatkowych przykładów Multcomp - 2-drożna Anova”. Ten przykład prowadzi do porównania poziomów napięcia w obrębie poziomów wool
.
Co jeśli chcę zrobić coś przeciwnego - porównać poziomy w wool
obrębie poziomów tension
? (W moim przypadku byłoby to porównanie poziomów leczenia (dwa - 0, 1) w ramach poziomów czasu (trzy - czerwiec, lipiec, sierpień).
Wymyśliłem następujący kod, aby to zrobić, ale wydaje się, że nie działa (patrz komunikat o błędzie poniżej).
Po pierwsze, z przykładu (z wool
i tension
zamieniły się miejscami):
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 = ":")
Od tego momentu mój własny kod:
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
źródło
time
predyktorem numerycznym. Podejrzewam, że chciałeś to jako czynnik.help("lsmeans", package = "lsmeans")
ivignette("using-lsmeans")
. Istnieje wiele dokumentacji i wielu przykładów.