Saya melakukan tes post-hoc pada model efek campuran linear dalam R
( lme4
paket). Saya menggunakan multcomp
paket ( glht()
fungsi) untuk melakukan tes post-hoc.
Desain eksperimental saya adalah tindakan berulang, dengan efek blok acak. Model ditentukan sebagai:
mymod <- lmer(variable ~ treatment * time + (1|block), data = mydata, REML = TRUE)
Daripada melampirkan data saya di sini, saya bekerja dari data yang disebut warpbreaks
dalam multcomp
paket.
data <- warpbreaks
warpbreaks$rand <- NA
Saya telah menambahkan variabel acak tambahan untuk meniru efek "blok" saya:
warpbreaks$rand <- rep(c("foo", "bar", "bee"), nrow(warpbreaks)/3)
Ini meniru model saya:
mod <- lmer(breaks ~ tension * wool + (1|rand), data = warpbreaks)
Saya menyadari contoh dalam " Contoh Multcomp Tambahan- 2 Way Anova" Contoh ini membawa Anda ke perbandingan tingkat ketegangan dalam level wool
.
Bagaimana jika saya ingin melakukan yang sebaliknya - membandingkan level wool
dalam level tension
? (Dalam kasus saya, ini akan membandingkan tingkat pengobatan (dua - 0, 1) dalam tingkat waktu (tiga - Juni, Juli, Agustus).
Saya datang dengan kode berikut untuk melakukannya, tetapi tampaknya tidak berfungsi (lihat pesan kesalahan di bawah).
Pertama, dari contoh (dengan wool
dan tension
bertukar tempat):
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 = ":")
Dari sini ke bawah, kode saya sendiri:
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
sumber
time
sebagai prediktor numerik. Saya menduga Anda menginginkannya sebagai faktor.help("lsmeans", package = "lsmeans")
danvignette("using-lsmeans")
. Ada banyak dokumentasi dan banyak contoh.