Ich führe Post-hoc-Tests an einem linearen Mischeffektmodell in R
( lme4
Paket) durch. Ich verwende multcomp
Paket ( glht()
Funktion), um die Post-Hoc-Tests durchzuführen.
Mein experimenteller Entwurf besteht aus wiederholten Messungen mit einem zufälligen Blockeffekt. Die Modelle sind wie folgt spezifiziert:
mymod <- lmer(variable ~ treatment * time + (1|block), data = mydata, REML = TRUE)
Anstatt meine Daten hier anzuhängen, arbeite ich an den warpbreaks
im multcomp
Paket aufgerufenen Daten .
data <- warpbreaks
warpbreaks$rand <- NA
Ich habe eine zusätzliche Zufallsvariable hinzugefügt, um meinen "Block" -Effekt nachzuahmen:
warpbreaks$rand <- rep(c("foo", "bar", "bee"), nrow(warpbreaks)/3)
Dies ahmt mein Modell nach:
mod <- lmer(breaks ~ tension * wool + (1|rand), data = warpbreaks)
Mir ist das Beispiel in " Zusätzliche Multcomp-Beispiele - 2-Wege-Anova" bekannt. Dieses Beispiel führt Sie zu einem Vergleich der Spannungsniveaus innerhalb der Niveaus von wool
.
Was ist, wenn ich das Gegenteil tun möchte - die Ebenen wool
innerhalb der Ebenen von vergleichen tension
? (In meinem Fall würde dies einen Vergleich der Behandlungsniveaus (zwei - 0, 1) innerhalb der Zeitniveaus (drei - Juni, Juli, August) bedeuten.
Ich habe mir dazu den folgenden Code ausgedacht, aber er scheint nicht zu funktionieren (siehe Fehlermeldung unten).
Zunächst aus dem Beispiel (mit wool
und tension
vertauschten Stellen):
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 = ":")
Von hier nach unten mein eigener 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
quelle
time
einem numerischen Prädiktor ist. Ich vermute, Sie wollten es als Faktor.help("lsmeans", package = "lsmeans")
undvignette("using-lsmeans")
. Es gibt viele Dokumentationen und viele Beispiele.