Pengujian post-hoc dalam multcomp :: glht untuk model efek-campuran (lme4) dengan interaksi

10

Saya melakukan tes post-hoc pada model efek campuran linear dalam R( lme4paket). Saya menggunakan multcomppaket ( 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 warpbreaksdalam multcomppaket.

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 wooldalam 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 wooldan tensionbertukar 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
Ashley Asmus
sumber

Jawaban:

6

Jauh lebih mudah menggunakan paket lsmeans

library(lsmeans)
lsmeans(mod, pairwise ~ tension | wool)
lsmeans(mod, pairwise ~ wool | tension)
Russ Lenth
sumber
Bagus, itu berhasil! Terima kasih. Catatan: kode ini hanya berfungsi untuk data saya setelah mengubah variabel pengulangan dari nilai numerik (3 & 6) menjadi nilai alfabet (A & B).
Yah, itu sangat penting! Karena itu model yang berbeda dengan timesebagai prediktor numerik. Saya menduga Anda menginginkannya sebagai faktor.
Russ Lenth
Bagaimana saya bisa menggeneralisasi lebih banyak prediktor? jika misalnya saya punya 3 alat prediksi bagaimana cara kerjanya?
bersenang-senang
1
@havefun Silakan lihat help("lsmeans", package = "lsmeans")dan vignette("using-lsmeans"). Ada banyak dokumentasi dan banyak contoh.
Russ Lenth
1
Hitung jumlah perbandingan yang Anda peroleh dengan masing-masing metode, mereka tidak sama. Baca juga penyesuaian beberapa pengujian. Saat Anda memiliki kelompok tes yang lebih besar, nilai P yang disesuaikan berbeda dengan untuk kelompok yang lebih kecil. Ketika Anda menggunakan oleh variabel, penyesuaian diterapkan secara terpisah untuk setiap set.
Russ Lenth