Cara mengatur kontras khusus dengan lmer di R

9

Saya menggunakan lmer di R untuk memeriksa efek kondisi ( cond) pada beberapa hasil. Berikut adalah beberapa data yang dibuat, di mana s adalah pengidentifikasi subjek dan a, bdan ckondisi.

library("tidyr")
library("dplyr")
set.seed(123)
temp <- data.frame(s = paste0("S", 1:30), 
                   a = rnorm(30, -2, 1), 
                   b = rnorm(30, -3, 1), 
                   c = rnorm(30, -4, 1)) 

Saya ingin membandingkan

  1. level ake rata-rata level bdan cdan
  2. level bke level c.

Pertanyaan saya adalah, bagaimana cara mengatur kontras untuk melakukan ini sedemikian rupa sehingga intersep mencerminkan rata-rata dari tiga kondisi dan dua estimasi yang dihitung secara langsung mencerminkan perbedaan sebagaimana didefinisikan dalam 1. dan 2.?

Saya mencoba

c1 <- cbind(c(-0.5, 0.25, 0.25), c(0, -0.5, 0.5))
gather(temp, cond, result, a, b, c) %>%
  lmer(result ~ cond + (1|s), data = ., contrasts = list(cond = c1))

di mana cond2tampaknya OK, tetapi cond1tidak.

Mengikuti Bagaimana menafsirkan kontras khusus ini? , Saya mencoba menggunakan invers umum, sebagai gantinya, tetapi perkiraan ini juga tidak masuk akal.

c2 <- t(ginv(c1))
gather(temp, cond, result, a, b, c) %>%
  lmer(result ~ cond + (1|s), data = ., contrasts = list(cond = c2))

Saya mencoba kontras Helmert juga, tetapi cara masih tidak cocok.

gather(temp, cond, result, a, b, c) %>%
  mutate(cond = factor(cond, levels = c("c", "b", "a"))) %>%
  lmer(result ~ cond + (1|s), data = ., contrasts = list(cond = contr.helmert))

Apa cara yang benar untuk melakukan ini?

M4RT1NK4
sumber
Ini kedengarannya seperti kontras Helmert (c adalah tingkat pertama, lalu b, lalu a).
Michael M
Saya mencoba Helmert juga, tetapi jumlahnya bukan berarti saya mencari. Saya telah mengedit pertanyaan untuk memasukkan kontras Helmert, terima kasih.
M4RT1NK4

Jawaban:

13

Untuk langkah-langkah berikut, kita membutuhkan bingkai data dalam format panjang. Frame data datberisi variabel dependen result, yang kategoris prediktor cond(tingkatan: a, b, dan c), dan faktor acak s.

library(tidyr)
dat <- gather(temp, cond, result, a, b, c)

Berikut ini, saya akan menggambarkan dua pendekatan untuk membuat matriks kontras yang sesuai dengan kondisi yang ingin Anda bandingkan:

  1. ab+c2
  2. bc

Kontras khusus

Matriks matsesuai dengan perbedaan level.

mat <- rbind(c(1, -0.5, -0.5),     # a vs. (b + c) / 2
             c(0, 1, -1))          # b vs. c

Untuk membuat matriks kontras aktual, kami menghitung invers umum dengan ginv(dari MASS).

library(MASS)
cMat <- ginv(mat)
#            [,1]          [,2]
# [1,]  0.6666667 -7.130169e-17
# [2,] -0.3333333  5.000000e-01
# [3,] -0.3333333 -5.000000e-01

Matriks kontras ini cMatdapat digunakan di lmer.

library(lme4)
res <- lmer(result ~ cond + (1|s), data = dat, 
            contrasts = list(cond = cMat))
coef(summary(res))    
#              Estimate Std. Error    t value
# (Intercept) -2.948115  0.0946025 -31.163182
# cond1        1.351517  0.2006822   6.734612
# cond2        1.153918  0.2317279   4.979625

Seperti yang Anda lihat, perkiraan efek tetap sesuai dengan perbedaan yang ditentukan di atas. Selanjutnya, intersep mewakili rata-rata keseluruhan.

Helmert kontras dengan contr.helmert

Anda juga dapat menggunakan contr.helmertfungsi bawaan untuk membuat matriks kontras.

cHelmert <- contr.helmert(3)
#   [,1] [,2]
# 1   -1   -1
# 2    1   -1
# 3    0    2

Namun, pesanan tidak sesuai dengan yang Anda tentukan dalam pertanyaan. Karenanya, kita harus membalik urutan kolom dan baris. Kolom pertama berkorespondensi dengan bvs. adan yang kedua berkorespondensi dengan crata-rata dari bdan a.

cHelmert2 <- cHelmert[c(3:1), 2:1]
#   [,1] [,2]
# 3    2    0
# 2   -1    1
# 1   -1   -1

Bandingkan matriks kontras cHelmert2untuk cMat. Anda akan melihat bahwa kolom adalah versi skala dari matriks lain.

Hasilnya lmeradalah:

library(lme4)
res2 <- lmer(result ~ cond + (1|s), data = dat, 
             contrasts = list(cond = cHelmert2))
coef(summary(res2))    
#               Estimate Std. Error    t value
# (Intercept) -2.9481150 0.09460250 -31.163182
# cond1        0.4505056 0.06689407   6.734612
# cond2        0.5769590 0.11586393   4.979625

t

Sven Hohenstein
sumber
Terima kasih banyak! Hanya untuk memastikan saya mengerti ini sekarang - jika saya ingin membandingkan level pertama dengan sisa level dalam variabel 4 level, matakankah c(1, -1/3, -1/3, -1/3)? Jadi saya selalu mengatur angka seperti pada rumus (a + (b + c + d) / 3) dan kemudian ginvmenskala dengan tepat sehingga koefisien langsung mencerminkan perbedaannya. Dan ketika Anda mengubah urutan dalam contoh Helmert, itu hanya untuk mencocokkan pertanyaan? Kalau tidak, hasilnya harus sama, terlepas dari urutan kontrasnya, bukan?
M4RT1NK4
@ M4RT1NK4 Formula Anda dan kontras yang sesuai adalah benar. Urutan kolom baru saja diubah agar sesuai dengan urutan kolom dalam pertanyaan. Namun, urutan baris adalah penting, karena level pertama adalah level referensi. Dalam contoh Anda, level referensi adalah level ketiga.
Sven Hohenstein
@ SvenHohenstein Saya punya pertanyaan terkait berdasarkan jawaban ini, keberatan melihat-lihat? stats.stackexchange.com/questions/357781/...
mat