14

Misalkan saya memiliki peserta, yang masing-masing memberikan respons 20 kali, 10 dalam satu kondisi dan 10 di lain. Saya cocok dengan model efek campuran linear yang membandingkan di setiap kondisi. Berikut adalah contoh yang dapat direplikasi mensimulasikan situasi ini menggunakan paket di : $N$ $Y$ $Y$ lme4R

library(lme4)
fml <- "~ condition + (condition | participant_id)"
d <- expand.grid(participant_id=1:40, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))

set.seed(23432)
d <- cbind(d, simulate(formula(fml), 
                       newparams=list(beta=c(0, .5), 
                                      theta=c(.5, 0, 0), 
                                      sigma=1), 
                       family=gaussian, 
                       newdata=d))

m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)

Model ini mmenghasilkan dua efek tetap (intersep dan lereng untuk kondisi), dan tiga efek acak (intersep acak partisipan, lereng acak partisipan untuk kondisi, dan korelasi intersep-slope).

Saya ingin membandingkan secara statistik ukuran varians mencegat acak oleh-peserta di seluruh kelompok yang ditentukan oleh condition(yaitu, menghitung komponen varians yang disorot dengan warna merah secara terpisah dalam kondisi kontrol dan eksperimental, kemudian menguji apakah perbedaan dalam ukuran komponen tidak nol). Bagaimana saya melakukan ini (lebih disukai dalam R)?

BONUS

Katakanlah modelnya sedikit lebih rumit: Para peserta masing-masing mengalami 10 rangsangan masing-masing 20 kali, 10 dalam satu kondisi dan 10 di lain. Dengan demikian, ada dua set efek acak silang: efek acak untuk peserta dan efek acak untuk stimulus. Berikut ini contoh yang dapat direproduksi:

library(lme4)
fml <- "~ condition + (condition | participant_id) + (condition | stimulus_id)"
d <- expand.grid(participant_id=1:40, stimulus_id=1:10, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))

set.seed(23432)
d <- cbind(d, simulate(formula(fml), 
                       newparams=list(beta=c(0, .5), 
                                      theta=c(.5, 0, 0, .5, 0, 0), 
                                      sigma=1), 
                       family=gaussian, 
                       newdata=d))

m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)

Saya ingin membandingkan secara statistik besarnya varian intersep oleh peserta secara acak di seluruh kelompok yang ditentukan oleh condition. Bagaimana saya melakukan itu, dan apakah prosesnya berbeda dari yang ada dalam situasi yang dijelaskan di atas?

EDIT

Untuk sedikit lebih spesifik tentang apa yang saya cari, saya ingin tahu:

Apakah pertanyaannya, "apakah respons rata-rata bersyarat dalam setiap kondisi (yaitu, nilai intersepsi acak di setiap kondisi) secara substansial berbeda satu sama lain, di luar yang kami harapkan karena kesalahan pengambilan sampel" adalah pertanyaan yang terdefinisi dengan baik (yaitu, apakah pertanyaan ini bahkan secara teoritis jawab)? Jika tidak, mengapa tidak?
Jika jawaban atas pertanyaan (1) adalah ya, bagaimana saya akan menjawabnya? Saya lebih suka Rimplementasi, tetapi saya tidak terikat dengan lme4paket - misalnya, sepertinya OpenMxpaket tersebut memiliki kemampuan untuk mengakomodasi analisis multi-grup dan multi-level ( https: //openmx.ssri.psu. edu / openmx-features ), dan ini sepertinya jenis pertanyaan yang harus dijawab dalam kerangka SEM.

r mixed-model lme4-nlme Patrick S. Forscher
sumber

1

@ MarkWhite, saya telah memperbarui pertanyaan untuk menanggapi komentar Anda. Maksud saya, saya ingin membandingkan deviasi standar dari intersep peserta ketika mereka memberikan respons dalam kondisi kontrol vs ketika mereka memberikan respons dalam kondisi eksperimental. Saya ingin melakukan ini secara statistik, yaitu, menguji apakah perbedaan dalam standar deviasi dari penyadapan berbeda dari 0.

Patrick S. Forscher

2

Saya menulis jawaban, tetapi akan tidur karena saya tidak yakin itu sangat berguna. Pertanyaannya adalah bahwa saya pikir seseorang tidak dapat melakukan apa yang Anda minta. Efek acak dari intersep adalah varians dalam rata-rata peserta ketika mereka berada dalam kondisi kontrol. Jadi kita tidak bisa melihat perbedaannya untuk pengamatan dalam kondisi eksperimental. Penyadapan didefinisikan di tingkat orang, dan kondisinya ada di tingkat pengamatan. Jika Anda mencoba membandingkan varian antara kondisi, saya akan memikirkan model heteroscedastic kondisional.

Mark White

2

Saya sedang mengerjakan revisi & kirim ulang untuk makalah di mana saya memiliki peserta yang memberikan tanggapan terhadap set rangsangan. Setiap peserta terpapar pada berbagai kondisi dan setiap stimulus menerima respons dalam berbagai kondisi - dengan kata lain, penelitian saya meniru pengaturan yang saya jelaskan dalam deskripsi "BONUS" saya. Dalam salah satu grafik saya, tampak bahwa respons peserta rata-rata memiliki variabilitas yang lebih besar di salah satu kondisi daripada yang lain. Peninjau meminta saya untuk menguji apakah ini benar.

Patrick S. Forscher

2

Silakan lihat di sini stats.stackexchange.com/questions/322213 untuk cara mengatur model lme4 dengan berbagai parameter varians untuk setiap level variabel pengelompokan. Saya tidak yakin bagaimana melakukan tes hipotesis apakah dua parameter varians sama; secara pribadi, saya selalu lebih suka melakukan bootstrap pada subjek dan rangsangan untuk mendapatkan interval kepercayaan, atau mungkin untuk menyiapkan semacam tes hipotesis permutasi (berbasis resampling).

Amoeba berkata Reinstate Monica

3

Saya setuju dengan komentar @MarkWhite bahwa pertanyaan "adalah perbedaan intersepsi acak yang secara substansial berbeda satu sama lain ..." paling tidak jelas dan paling buruk tidak masuk akal, karena intersep perlu merujuk pada nilai-Y dalam satu kelompok tertentu ( grup menetapkan nilai 0), jadi membandingkan "intersep" di seluruh grup secara tegas tidak masuk akal. Saya pikir cara yang lebih baik untuk mengulangi pertanyaan Anda, seperti yang saya pahami, akan menjadi sesuatu seperti: "Apakah variasi tanggapan rata-rata kondisional peserta dalam kondisi A vs kondisi B tidak setara?"

Jake Westfall

6

Ada lebih dari satu cara untuk menguji hipotesis ini. Misalnya, prosedur yang diuraikan oleh @amoeba harus berfungsi. Tapi menurut saya cara paling sederhana dan paling bijaksana untuk mengujinya adalah menggunakan uji rasio kemungkinan lama yang baik membandingkan dua model bersarang. Satu-satunya bagian yang berpotensi sulit dari pendekatan ini adalah mengetahui cara mengatur pasangan model sehingga menjatuhkan satu parameter akan dengan bersih menguji hipotesis yang diinginkan dari varian yang tidak sama. Di bawah ini saya jelaskan bagaimana melakukannya.

Jawaban singkat

Beralih ke pengkodean kontras (jumlah ke nol) untuk variabel independen Anda dan kemudian lakukan uji rasio kemungkinan membandingkan model lengkap Anda dengan model yang memaksa korelasi antara lereng acak dan penyadapan acak menjadi 0:

# switch to numeric (not factor) contrast codes
d$contrast <- 2*(d$condition == 'experimental') - 1

# reduced model without correlation parameter
mod1 <- lmer(sim_1 ~ contrast + (contrast || participant_id), data=d)

# full model with correlation parameter
mod2 <- lmer(sim_1 ~ contrast + (contrast | participant_id), data=d)

# likelihood ratio test
anova(mod1, mod2)

Penjelasan / intuisi visual

Agar jawaban ini masuk akal, Anda perlu memiliki pemahaman intuitif tentang apa perbedaan nilai parameter korelasi untuk data yang diamati. Pertimbangkan garis regresi subjek-spesifik (berbeda-beda). Pada dasarnya, parameter korelasi mengontrol apakah garis regresi peserta "menyebar ke kanan" (korelasi positif) atau "menyebar ke kiri" (korelasi negatif) relatif terhadap titik , di mana X adalah kode kontras Anda yang independen variabel. Salah satu dari ini menyiratkan varians yang tidak sama dalam tanggapan rata-rata bersyarat peserta. Ini diilustrasikan di bawah ini: $X=0$

Dalam plot ini, kami mengabaikan beberapa pengamatan yang kami miliki untuk setiap subjek dalam setiap kondisi dan alih-alih hanya plot dua rata-rata subjek masing-masing, dengan garis yang menghubungkan mereka, mewakili kemiringan acak subjek itu. (Ini terdiri dari data dari 10 subjek hipotetis, bukan data yang diposting di OP.)

Di kolom di sebelah kiri, di mana ada korelasi kemiringan-intersepsi negatif yang kuat, garis regresi menyebar ke kiri relatif ke titik . Seperti yang dapat Anda lihat dengan jelas pada gambar, ini mengarah ke variasi yang lebih besar dalam rata-rata acak subjek dalam kondisi daripada dalam kondisi . $X=0$ $X=-1$ $X=1$

Kolom di sebelah kanan menunjukkan gambar cermin terbalik dari pola ini. Dalam hal ini ada perbedaan yang lebih besar dalam mean acak subjek dalam kondisi daripada dalam kondisi . $X=1$ $X=-1$

Kolom di tengah menunjukkan apa yang terjadi ketika lereng acak dan intersep acak tidak berkorelasi. Ini berarti bahwa garis regresi melebar ke kiri persis sebanyak yang melebar ke kanan, relatif ke titik . Ini menyiratkan bahwa varians rata-rata subjek dalam kedua kondisi adalah sama. $X=0$

Sangat penting di sini bahwa kami telah menggunakan skema pengkodean kontras jumlah-ke-nol, bukan kode dummy (yaitu, tidak menetapkan grup pada $X=0$ vs $X=1$ ). Hal ini hanya di bawah kontras skema pengkodean bahwa kita memiliki hubungan ini dimana varians adalah sama jika dan hanya jika lereng-intercept korelasi adalah 0. Gambar di bawah mencoba untuk membangun intuisi:

Apa yang ditunjukkan gambar ini adalah dataset yang sama persis di kedua kolom, tetapi dengan variabel independen yang dikodekan dua cara berbeda. Di kolom di sebelah kiri kami menggunakan kode kontras - ini persis situasi dari gambar pertama. Di kolom di sebelah kanan kami menggunakan kode dummy. Ini mengubah arti intersep - sekarang intersep mewakili respons prediksi subjek dalam kelompok kontrol. Panel bawah menunjukkan konsekuensi dari perubahan ini, yaitu, bahwa korelasi slope-intercept tidak lagi mendekati 0, meskipun datanya sama dalam arti yang dalam dan varians kondisional sama dalam kedua kasus. Jika ini masih tidak masuk akal, mempelajari jawaban saya sebelumnya di mana saya berbicara lebih banyak tentang fenomena ini dapat membantu.

Bukti

Biarkan $y_{ijk}$ menjadi respons $j$ dari subjek ke- $i$ dalam kondisi $k$ . (Kami hanya memiliki dua kondisi di sini, jadi $k$ hanya 1 atau 2.) Kemudian model campuran dapat ditulis

y_{i j k} = α_{i} + β_{i} x_{k} + e_{i j k},

$y_{ijk} = \alpha_i + \beta_ix_k + e_{ijk},$ mana

α_{i}

$\alpha_i$ adalah subjeknya ' penyadapan acak dan memiliki varian

σ_{α}^{2}

$\sigma^2_\alpha$ ,

β_{i}

$\beta_i$ adalah kemiringan acak subjek dan memiliki varian

σ_{β}^{2}

$\sigma^2_\beta$ ,

e_{i j k}

$e_{ijk}$ adalah istilah kesalahan tingkat observasi, dan

cov (α_{i}, β_{i}) = σ_{α β}

$\text{cov}(\alpha_i, \beta_i)=\sigma_{\alpha\beta}$ .

Kami ingin menunjukkan bahwa

var (α_{i} + β_{i} x_{1}) = var (α_{i} + β_{i} x_{2}) \Leftrightarrow σ_{α β} = 0.

$\text{var}(\alpha_i + \beta_ix_1) = \text{var}(\alpha_i + \beta_ix_2) \Leftrightarrow \sigma_{\alpha\beta}=0.$

Dimulai dengan sisi kiri dari implikasi ini, kita memiliki

\begin{aligned} var (α_{i} + β_{i} x_{1}) & = var (α_{i} + β_{i} x_{2}) \\ σ_{α}^{2} + x_{1}^{2} σ_{β}^{2} + 2 x_{1} σ_{α β} & = σ_{α}^{2} + x_{2}^{2} σ_{β}^{2} + 2 x_{2} σ_{α β} \\ σ_{β}^{2} (x_{1}^{2} - x_{2}^{2}) + 2 σ_{α β} (x_{1} - x_{2}) & = 0. \end{aligned}

$\begin{aligned} \text{var}(\alpha_i + \beta_ix_1) &= \text{var}(\alpha_i + \beta_ix_2) \\ \sigma^2_\alpha + x^2_1\sigma^2_\beta + 2x_1\sigma_{\alpha\beta} &= \sigma^2_\alpha + x^2_2\sigma^2_\beta + 2x_2\sigma_{\alpha\beta} \\ \sigma^2_\beta(x_1^2 - x_2^2) + 2\sigma_{\alpha\beta}(x_1 - x_2) &= 0. \end{aligned}$

Kode kontras jumlah-ke-nol menyiratkan bahwa $x_1 + x_2 = 0$ dan $x_1^2 = x_2^2 = x^2$ . Kemudian kita dapat mengurangi baris terakhir di atas menjadi

\begin{aligned} σ_{β}^{2} (x^{2} - x^{2}) + 2 σ_{α β} (x_{1} + x_{1}) & = 0 \\ σ_{α β} & = 0, \end{aligned}

$\begin{aligned} \sigma^2_\beta(x^2 - x^2) + 2\sigma_{\alpha\beta}(x_1 + x_1) &= 0 \\ \sigma_{\alpha\beta} &= 0, \end{aligned}$ yang ingin kami buktikan. (Untuk menetapkan arah lain dari implikasi, kita bisa mengikuti langkah-langkah yang sama ini secara terbalik.)

Untuk mengulangi, ini menunjukkan bahwa jika variabel independen adalah kode kontras (jumlah ke nol) , maka varian rata-rata acak subjek dalam setiap kondisi adalah sama jika dan hanya jika korelasi antara lereng acak dan penyadapan acak adalah 0. Kuncinya titik ambil dari semua ini adalah bahwa menguji hipotesis nol bahwa $\sigma_{\alpha\beta} = 0$ akan menguji hipotesis nol dari varian yang sama yang dijelaskan oleh OP.

Ini TIDAK berfungsi jika variabel independen, katakanlah, kode dummy. Secara khusus, jika kita memasukkan nilai $x_1=0$ dan $x_2=1$ ke dalam persamaan di atas, kita menemukan bahwa

var (α_{i}) = var (α_{i} + β_{i}) \Leftrightarrow σ_{α β} = - \frac{σ_{β}^{2}}{2} .

$\text{var}(\alpha_i) = \text{var}(\alpha_i + \beta_i) \Leftrightarrow \sigma_{\alpha\beta} = -\frac{\sigma^2_\beta}{2}.$

Jake Westfall
sumber

Ini sudah merupakan jawaban yang hebat, terima kasih! Saya pikir ini yang paling dekat dengan menjawab pertanyaan saya, jadi saya menerimanya dan memberi Anda hadiah (akan segera kedaluwarsa), tapi saya ingin melihat pembenaran aljabar jika Anda punya waktu dan energi untuk itu.

Patrick S. Forscher

1

@ PatrickS.Forscher Saya baru saja menambahkan bukti

Jake Westfall

1

@JakeWestfall Dalam contoh mainan saya, subjek membalik tanggapan dalam dua kondisi. Jika subjek memiliki respons

dalam kondisi A dan

dalam kondisi B, lalu apa nilai BLUP dari intersep acak untuk subjek ini saat kami menggunakan model? Saya pikir itu hanya bisa 0. Jika semua subjek memiliki BLUP sama dengan nol, maka varians mencegat acak juga nol. Jadi model ini tidak cocok dengan contoh mainan ini sama sekali. Sebaliknya, model didefinisikan di atas melalui akan memiliki dua BLUPs untuk setiap mata pelajaran, dan mereka dapat dengan mudah menjadi

dan

. Apakah saya melewatkan sesuatu di sini?

a

$a$

- a

$-a$ (1 | subject)dummy

a

$a$

- a

$-a$

Amuba mengatakan Reinstate Monica

1

Saya melihat sekarang bahwa Anda benar @amoeba, terima kasih telah menjelaskan. Saya akan mengedit jawaban saya sesuai.

Jake Westfall

1

@amoeba Anda benar bahwa kemungkinan BLUP dapat berkorelasi bahkan tanpa parameter korelasi dalam model. Tetapi saya percaya bahwa untuk tujuan pengujian prosedur masih berfungsi sebagaimana dimaksud (misalnya, ia memiliki tingkat kesalahan tipe 1 nominal) karena hanya model dengan parameter korelasi yang dapat memasukkan itu ke dalam fungsi kemungkinan dan dengan demikian "menerima kredit" untuk itu . Artinya, bahkan jika BLUP keluar berkorelasi dalam model yang lebih sederhana, masih seolah-olah efeknya tidak berkorelasi sejauh kemungkinan total yang bersangkutan, sehingga uji LR akan bekerja. Saya pikir :)

Jake Westfall

6

Anda dapat menguji signifikansi, parameter model, dengan bantuan perkiraan interval kepercayaan untuk fungsi paket lme4 confint.merMod.

bootstrap (lihat Interval Keyakinan misalnya dari bootstrap )

> confint(m, method="boot", nsim=500, oldNames= FALSE)
Computing bootstrap confidence intervals ...
                                                           2.5 %     97.5 %
sd_(Intercept)|participant_id                         0.32764600 0.64763277
cor_conditionexperimental.(Intercept)|participant_id -1.00000000 1.00000000
sd_conditionexperimental|participant_id               0.02249989 0.46871800
sigma                                                 0.97933979 1.08314696
(Intercept)                                          -0.29669088 0.06169473
conditionexperimental                                 0.26539992 0.60940435

profil kemungkinan (lihat misalnya Apa hubungan antara kemungkinan profil dan interval kepercayaan? )

> confint(m, method="profile", oldNames= FALSE)
Computing profile confidence intervals ...
                                                          2.5 %     97.5 %
sd_(Intercept)|participant_id                         0.3490878 0.66714551
cor_conditionexperimental.(Intercept)|participant_id -1.0000000 1.00000000
sd_conditionexperimental|participant_id               0.0000000 0.49076950
sigma                                                 0.9759407 1.08217870
(Intercept)                                          -0.2999380 0.07194055
conditionexperimental                                 0.2707319 0.60727448

Ada juga metode 'Wald'tetapi ini diterapkan untuk efek tetap saja.
Ada juga beberapa jenis jenis anova (rasio kemungkinan) dalam paket lmerTestyang dinamai ranova. Tapi sepertinya aku tidak masuk akal dari ini. Distribusi perbedaan dalam logLikelihood, ketika hipotesis nol (nol varians untuk efek acak) benar tidak didistribusikan secara chi-square (mungkin ketika jumlah peserta dan uji coba tinggi, uji rasio kemungkinan kemungkinan masuk akal).

Varians dalam kelompok tertentu

Untuk mendapatkan hasil varians dalam grup tertentu, Anda dapat melakukan reparameterisasi

# different model with alternative parameterization (and also correlation taken out) 
fml1 <- "~ condition + (0 + control + experimental || participant_id) "

Di mana kami menambahkan dua kolom ke kerangka data ^{(ini hanya diperlukan jika Anda ingin mengevaluasi 'kontrol' yang tidak berkorelasi dan 'eksperimental', fungsinya (0 + condition || participant_id)tidak akan mengarah pada evaluasi berbagai faktor dalam kondisi sebagai tidak berkorelasi)}

#adding extra columns for control and experimental
d <- cbind(d,as.numeric(d$condition=='control'))
d <- cbind(d,1-as.numeric(d$condition=='control'))
names(d)[c(4,5)] <- c("control","experimental")

Sekarang lmerakan memberikan perbedaan untuk kelompok yang berbeda

> m <- lmer(paste("sim_1 ", fml1), data=d)
> m
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: paste("sim_1 ", fml1)
   Data: d
REML criterion at convergence: 2408.186
Random effects:
 Groups           Name         Std.Dev.
 participant_id   control      0.4963  
 participant_id.1 experimental 0.4554  
 Residual                      1.0268  
Number of obs: 800, groups:  participant_id, 40
Fixed Effects:
          (Intercept)  conditionexperimental  
               -0.114                  0.439

Dan Anda dapat menerapkan metode profil untuk ini. Misalnya sekarang confint memberikan interval kepercayaan untuk kontrol dan varian ekserimental.

> confint(m, method="profile", oldNames= FALSE)
Computing profile confidence intervals ...
                                    2.5 %     97.5 %
sd_control|participant_id       0.3490873 0.66714568
sd_experimental|participant_id  0.3106425 0.61975534
sigma                           0.9759407 1.08217872
(Intercept)                    -0.2999382 0.07194076
conditionexperimental           0.1865125 0.69149396

Kesederhanaan

Anda dapat menggunakan fungsi kemungkinan untuk mendapatkan perbandingan yang lebih maju, tetapi ada banyak cara untuk membuat perkiraan di sepanjang jalan (mis. Anda bisa melakukan tes anova / lrt konservatif, tetapi apakah itu yang Anda inginkan?).

Pada titik ini membuat saya bertanya-tanya apa sebenarnya poin dari perbandingan (tidak begitu umum) antar varian. Saya bertanya-tanya apakah itu mulai menjadi terlalu canggih. Mengapa perbedaan antara varian bukan rasio antara varian (yang terkait dengan distribusi-F klasik)? Mengapa tidak melaporkan interval kepercayaan saja? Kita perlu mengambil langkah mundur, dan mengklarifikasi data dan cerita yang seharusnya diceritakan, sebelum masuk ke jalur lanjutan yang mungkin berlebihan dan tidak berhubungan dengan masalah statistik dan pertimbangan statistik yang sebenarnya menjadi topik utama.

Saya bertanya-tanya apakah seseorang harus melakukan lebih dari sekadar menyatakan interval kepercayaan (yang sebenarnya bisa mengatakan lebih dari tes hipotesis. Tes hipotesis memberikan jawaban ya tidak tetapi tidak ada informasi tentang sebaran populasi yang sebenarnya. Berikan data yang cukup yang Anda bisa membuat sedikit perbedaan untuk dilaporkan sebagai perbedaan signifikan). Untuk masuk lebih dalam ke masalah ini (untuk tujuan apa pun), saya percaya, membutuhkan pertanyaan penelitian yang lebih spesifik (didefinisikan secara sempit) untuk memandu mesin matematika untuk membuat penyederhanaan yang tepat (bahkan ketika perhitungan yang tepat mungkin dilakukan atau ketika itu bisa diperkirakan dengan simulasi / bootstrap, itupun dalam beberapa pengaturan masih memerlukan beberapa interpretasi yang sesuai). Bandingkan dengan tes eksak Fisher untuk menyelesaikan pertanyaan (khusus) (tentang tabel kontingensi),

Contoh sederhana

Untuk memberikan contoh kesederhanaan yang mungkin saya tunjukkan di bawah ini perbandingan (dengan simulasi) dengan penilaian sederhana dari perbedaan antara dua varian kelompok berdasarkan pada uji-F yang dilakukan dengan membandingkan perbedaan dalam tanggapan rata-rata individu dan dilakukan dengan membandingkan varians yang diturunkan dari model campuran.

$j$

{\hat{Y}}_{i, j} \sim N (μ_{j}, σ_{j}^{2} + \frac{σ_{ϵ}^{2}}{10})

$\hat{Y}_{i,j} \sim N(\mu_j, \sigma_j^2 + \frac{\sigma_{\epsilon}^2}{10})$

$\sigma_\epsilon$ $\sigma_{j}$ $j = \lbrace 1,2 \rbrace$

Anda dapat melihat ini dalam simulasi grafik di bawah ini selain untuk skor-F berdasarkan sampel berarti skor-F dihitung berdasarkan variasi yang diprediksi (atau jumlah kesalahan kuadrat) dari model.

^{$\sigma_{j=1} = \sigma_{j=2} = 0.5$ $\sigma_\epsilon=1$}

Anda dapat melihat bahwa ada beberapa perbedaan. Perbedaan ini mungkin karena fakta bahwa model linear efek campuran mendapatkan jumlah kesalahan kuadrat (untuk efek acak) dengan cara yang berbeda. Dan istilah kesalahan kuadrat ini tidak (lagi) dinyatakan dengan baik sebagai distribusi Chi-kuadrat sederhana, tetapi masih terkait erat dan mereka dapat diperkirakan.

$\sigma_{j=1} \neq \sigma_{j=2}$ $\hat{Y}_{i,j}$ $\sigma_j$ $\sigma_\epsilon$

^{$\sigma_{j=1} = 0.5$ $\sigma_{j=2} = 0.25$ $\sigma_\epsilon=1$}

Jadi model berdasarkan pada sarana sangat tepat. Tapi itu kurang kuat. Ini menunjukkan bahwa strategi yang tepat tergantung pada apa yang Anda inginkan / butuhkan.

^{Dalam contoh di atas ketika Anda menetapkan batas ekor kanan pada 2.1 dan 3.1, Anda mendapatkan sekitar 1% dari populasi dalam kasus varians yang sama (resp 103 dan 104 dari 10.000 kasus) tetapi dalam kasus varians yang tidak setara, batas-batas ini berbeda. banyak (memberikan 5334 dan 6716 kasus)}

kode:

set.seed(23432)

# different model with alternative parameterization (and also correlation taken out)
fml1 <- "~ condition + (0 + control + experimental || participant_id) "
fml <- "~ condition + (condition | participant_id)"

n <- 10000

theta_m <- matrix(rep(0,n*2),n)
theta_f <- matrix(rep(0,n*2),n)

# initial data frame later changed into d by adding a sixth sim_1 column
ds <- expand.grid(participant_id=1:40, trial_num=1:10)
ds <- rbind(cbind(ds, condition="control"), cbind(ds, condition="experimental"))
  #adding extra columns for control and experimental
  ds <- cbind(ds,as.numeric(ds$condition=='control'))
  ds <- cbind(ds,1-as.numeric(ds$condition=='control'))
  names(ds)[c(4,5)] <- c("control","experimental")

# defining variances for the population of individual means
stdevs <- c(0.5,0.5) # c(control,experimental)

pb <- txtProgressBar(title = "progress bar", min = 0,
                    max = n, style=3)
for (i in 1:n) {

  indv_means <- c(rep(0,40)+rnorm(40,0,stdevs[1]),rep(0.5,40)+rnorm(40,0,stdevs[2]))
  fill <- indv_means[d[,1]+d[,5]*40]+rnorm(80*10,0,sqrt(1)) #using a different way to make the data because the simulate is not creating independent data in the two groups 
  #fill <- suppressMessages(simulate(formula(fml), 
  #                     newparams=list(beta=c(0, .5), 
  #                                    theta=c(.5, 0, 0), 
  #                                    sigma=1), 
  #                     family=gaussian, 
  #                     newdata=ds))
  d <- cbind(ds, fill)
  names(d)[6] <- c("sim_1")


  m <- lmer(paste("sim_1 ", fml1), data=d)
  m
  theta_m[i,] <- m@theta^2

  imeans <- aggregate(d[, 6], list(d[,c(1)],d[,c(3)]), mean)
  theta_f[i,1] <- var(imeans[c(1:40),3])
  theta_f[i,2] <- var(imeans[c(41:80),3])

  setTxtProgressBar(pb, i)
}
close(pb)

p1 <- hist(theta_f[,1]/theta_f[,2], breaks = seq(0,6,0.06))       
fr <- theta_m[,1]/theta_m[,2]
fr <- fr[which(fr<30)]
p2 <- hist(fr, breaks = seq(0,30,0.06))



plot(-100,-100, xlim=c(0,6), ylim=c(0,800), 
     xlab="F-score", ylab = "counts [n out of 10 000]")
plot( p1, col=rgb(0,0,1,1/4), xlim=c(0,6), ylim=c(0,800), add=T)  # means based F-score
plot( p2, col=rgb(1,0,0,1/4), xlim=c(0,6), ylim=c(0,800), add=T)  # model based F-score
fr <- seq(0, 4, 0.01)
lines(fr,df(fr,39,39)*n*0.06,col=1)
legend(2, 800, c("means based F-score","mixed regression based F-score"), 
       fill=c(rgb(0,0,1,1/4),rgb(1,0,0,1/4)),box.col =NA, bg = NA)
legend(2, 760, c("F(39,39) distribution"), 
       lty=c(1),box.col = NA,bg = NA)
title(expression(paste(sigma[1]==0.5, " , ", sigma[2]==0.5, " and ", sigma[epsilon]==1)))

Sextus Empiricus
sumber

Itu berguna tetapi tampaknya tidak menjawab pertanyaan tentang bagaimana membandingkan varian dalam dua kondisi.

Amoeba berkata Reinstate Monica

@amoeba Saya menemukan bahwa jawaban ini memberikan inti masalah (tentang menguji komponen varians acak). Apa yang diinginkan OP sebenarnya sulit untuk dibaca di seluruh teks. Apa yang dimaksud dengan "varian intersepsi acak"? (jamak dalam kaitannya dengan mencegat membingungkan saya) Satu kasus yang mungkin adalah menggunakan model sim_1 ~ condition + (0 + condition | participant_id)"di mana Anda mendapatkan parameterisasi menjadi dua parameter (satu untuk setiap kelompok) daripada dua parameter satu untuk mencegat dan satu untuk efek (yang perlu digabungkan untuk grup).

Sextus Empiricus

Setiap subjek memiliki beberapa respons rata-rata dalam kondisi A dan beberapa respons rata-rata dalam kondisi B. Pertanyaannya adalah apakah varians antar subjek dalam A berbeda dari varians antar subjek dalam B.

Amuba mengatakan Reinstate Monica

Ini tidak menyelesaikan tugas yang diajukan dalam judul "Bandingkan komponen varians acak di seluruh level variabel pengelompokan". Saya perhatikan ada kesalahan ketik yang membingungkan di tubuh pertanyaan, yang telah saya perbaiki. Saya juga mencoba untuk lebih memperjelas kata-kata dari pertanyaan itu.

Patrick S. Forscher

Dimungkinkan untuk menjawab pertanyaan menggunakan car::linearHypothesisTest( math.furman.edu/~dcs/courses/math47/R/library/car/html/… ), yang memungkinkan pengguna untuk menguji hipotesis sewenang-wenang dengan model yang sesuai. Namun, saya harus menggunakan metode @ amoeba untuk mendapatkan kedua intersep acak dalam model model yang sama sehingga mereka dapat dibandingkan dengan fungsi ini. Saya juga sedikit tidak yakin dengan validitas metode ini.

Patrick S. Forscher

5

Salah satu cara yang relatif mudah adalah dengan menggunakan tes rasio kemungkinan melalui anovaseperti yang dijelaskan dalam lme4FAQ .

Kita mulai dengan model penuh di mana varians tidak dibatasi (yaitu, dua varians yang berbeda diperbolehkan) dan kemudian cocok dengan satu model dibatasi di mana kedua varians diasumsikan sama. Kami hanya membandingkannya dengan anova()(perhatikan bahwa saya mengatur REML = FALSEmeskipun REML = TRUEdengan anova(..., refit = FALSE)sepenuhnya layak ).

m_full <- lmer(sim_1 ~ condition + (condition | participant_id), data=d, REML = FALSE)
summary(m_full)$varcor
 # Groups         Name                  Std.Dev. Corr  
 # participant_id (Intercept)           0.48741        
 #                conditionexperimental 0.26468  -0.419
 # Residual                             1.02677     

m_red <- lmer(sim_1 ~ condition + (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

anova(m_full, m_red)
# Data: d
# Models:
# m_red: sim_1 ~ condition + (1 | participant_id)
# m_full: sim_1 ~ condition + (condition | participant_id)
#        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
# m_red   4 2396.6 2415.3 -1194.3   2388.6                         
# m_full  6 2398.7 2426.8 -1193.3   2386.7 1.9037      2      0.386

Namun, tes ini cenderung konservatif . Misalnya, FAQ mengatakan:

Perlu diingat bahwa uji hipotesis nol berbasis LRT adalah konservatif ketika nilai nol (seperti σ2 = 0) berada pada batas ruang yang layak; dalam kasus paling sederhana (varians efek acak tunggal), nilai-p kira-kira dua kali lebih besar dari seharusnya (Pinheiro dan Bates 2000).

Ada beberapa alternatif:

$\chi^2$
Simulasikan distribusi yang benar menggunakan RLRsim(seperti juga dijelaskan dalam FAQ).

Saya akan menunjukkan opsi kedua sebagai berikut:

library("RLRsim")
## reparametrize model so we can get one parameter that we want to be zero:
afex::set_sum_contrasts() ## warning, changes contrasts globally
d <- cbind(d, difference = model.matrix(~condition, d)[,"condition1"])

m_full2 <- lmer(sim_1 ~ condition + (difference | participant_id), data=d, REML = FALSE)
all.equal(deviance(m_full), deviance(m_full2))  ## both full models are identical

## however, we need the full model without correlation!
m_full2b <- lmer(sim_1 ~ condition + (1| participant_id) + 
                   (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_full2b)$varcor
 # Groups           Name        Std.Dev.
 # participant_id   (Intercept) 0.44837 
 # participant_id.1 difference  0.13234 
 # Residual                     1.02677 

## model that only has random effect to be tested
m_red <- update(m_full2b,  . ~ . - (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name       Std.Dev.
 # participant_id difference 0.083262
 # Residual                  1.125116

## Null model 
m_null <- update(m_full2b,  . ~ . - (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_null)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

exactRLRT(m_red, m_full2b, m_null)
# Using restricted likelihood evaluated at ML estimators.
# Refit with method="REML" for exact results.
# 
#   simulated finite sample distribution of RLRT.
#   
#   (p-value based on 10000 simulated values)
# 
# data:  
# RLRT = 1.9698, p-value = 0.0719

Seperti yang bisa kita lihat, hasilnya menunjukkan bahwa dengan REML = TRUEkita akan mendapatkan hasil yang tepat. Tapi ini dibiarkan sebagai latihan untuk pembaca.

Mengenai bonus, saya tidak yakin apakah RLRsimmemungkinkan pengujian simultan beberapa komponen, tetapi jika demikian, ini dapat dilakukan dengan cara yang sama.

Tanggapan terhadap komentar:

Jadi benar, maka, secara umum kemiringan acak $\theta_X$ memungkinkan intersepsi acak $\theta_0$ bervariasi di berbagai tingkatan $X$ ?

Saya tidak yakin pertanyaan ini dapat menerima jawaban yang masuk akal.

Penyadapan acak memungkinkan perbedaan istimewa pada tingkat keseluruhan untuk setiap tingkat faktor pengelompokan. Misalnya, jika variabel dependen adalah waktu respons, beberapa peserta lebih cepat dan beberapa lebih lambat.
Lereng acak memungkinkan setiap tingkat faktor pengelompokan efek istimewa dari faktor yang diperkirakan lereng acak. Misalnya, jika faktornya adalah kongruensi, maka beberapa peserta dapat memiliki efek kongruensi yang lebih tinggi daripada yang lain.

Jadi apakah random-slope memengaruhi intercept acak? Dalam beberapa hal ini mungkin masuk akal, karena mereka memungkinkan setiap tingkat faktor pengelompokan efek yang sama sekali istimewa untuk setiap kondisi. Pada akhirnya, kami memperkirakan dua parameter istimewa untuk dua kondisi. Namun, saya pikir perbedaan antara tingkat keseluruhan yang ditangkap oleh intersep dan efek khusus kondisi yang ditangkap oleh lereng acak adalah penting dan kemudian lereng acak tidak dapat benar-benar memengaruhi intersep acak. Namun, masih memungkinkan setiap tingkat faktor pengelompokan menjadi istimewa secara terpisah untuk setiap tingkat kondisi.

Namun demikian, pengujian saya masih melakukan apa yang diinginkan pertanyaan aslinya. Ini menguji apakah perbedaan varian antara kedua kondisi adalah nol. Jika nol, maka varians di kedua kondisi sama. Dengan kata lain, hanya jika tidak perlu untuk kemiringan acak adalah varians dalam kedua kondisi identik. Saya harap itu masuk akal.

Henrik
sumber

1

Anda menggunakan kontras pengobatan ( contr.treatment) yang mana kondisi kontrol adalah referensi (yaitu, yang dihitung intersep acak). Parametrization yang saya usulkan, saya menggunakan jumlah kontras (yaitu, contr.sum) dan intersep adalah mean besar. Saya merasa lebih masuk akal untuk menguji apakah perbedaannya nol ketika intersepnya adalah rata-rata utama alih-alih kondisi kontrol (tetapi penulisan yang dilakukan menunjukkan bahwa itu mungkin relatif tidak penting). Anda mungkin ingin membaca halaman 24 hingga 26 dari: singmann.org/download/publications/…

Henrik

1

Terima kasih! Namun, pertanyaan saya sedikit berbeda: (1) Jawaban Anda sepertinya menyiratkan bahwa pertanyaan saya berkurang menjadi "adalah kemiringan acak untuk kondisi yang berbeda dari 0". Apakah ini benar? (2) Jika jawaban untuk (1) adalah "ya", ini menyarankan interpretasi lain dari kemiringan acak untuk condition: ini memungkinkan intersepsi acak bervariasi di berbagai tingkat condition. Apakah ini benar?

Patrick S. Forscher

2

2 ¢ saya: Contoh tandingan @amoeba dengan prosedur yang diajukan Henrik adalah benar. Henrik hampir benar, tetapi ia membandingkan model yang salah. Model perbandingan bahwa pertanyaan jawabannya Patrick adalah perbandingan antara model Henrik disebut m_fullvs m_full2b. Yaitu: varian tanggapan rata-rata kondisional peserta dalam A vs B tidak sama jika jika korelasi intersep-slope acak adalah nol --- yang penting, di bawah parameterisasi pengkodean kontras jumlah ke nol . Pengujian varians kemiringan acak tidak diperlukan. Mencoba memikirkan bagaimana menjelaskan ini dengan ringkas ...

Jake Westfall

2

Ini sebenarnya bukan penjelasan yang tepat, tetapi mempelajari jawaban saya di sini mungkin sedikit menjelaskan masalah ini. Pada dasarnya, parameter korelasi mengontrol apakah garis regresi peserta "melebar ke kanan" (koreksi positif) atau "melebar ke kiri" (kor negatif). Salah satu dari ini menyiratkan varians yang tidak sama dalam tanggapan rata-rata bersyarat peserta. Penjumlahan jumlah-ke-nol kemudian memastikan bahwa kami sedang mencari korelasi pada titik yang tepat di X

Jake Westfall

2

Saya akan mempertimbangkan memposting jawaban dengan gambar jika saya dapat menemukan waktu ...

Jake Westfall

5

Model Anda

m = lmer(sim_1 ~ condition + (condition | participant_id), data=d)

sudah memungkinkan varians lintas subjek dalam kondisi kontrol berbeda dari varians lintas subjek dalam kondisi eksperimental. Ini dapat dibuat lebih eksplisit dengan parametriisasi ulang yang setara:

m = lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=d)

Matriks kovarian acak sekarang memiliki interpretasi yang lebih sederhana:

Random effects:
 Groups         Name                  Variance Std.Dev. Corr
 participant_id conditioncontrol      0.2464   0.4963       
                conditionexperimental 0.2074   0.4554   0.83

Di sini kedua varians tersebut adalah tepat dua varians yang Anda minati: varians [lintas subjek] dari respons rata-rata bersyarat dalam kondisi kontrol dan sama dalam kondisi eksperimental. Dalam dataset disimulasikan Anda, mereka adalah 0,25 dan 0,21. Perbedaannya diberikan oleh

delta = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]

dan sama dengan 0,039. Anda ingin menguji apakah secara signifikan berbeda dari nol.

EDIT: Saya menyadari bahwa tes permutasi yang saya jelaskan di bawah ini salah; itu tidak akan berfungsi sebagaimana dimaksud jika sarana dalam kondisi kontrol / eksperimen tidak sama (karena pengamatan tidak dapat ditukar dengan nol). Mungkin ide yang lebih baik untuk mem-bootstrap subyek (atau subjek / item dalam kasing Bonus) dan mendapatkan interval kepercayaan untuk delta.

Saya akan mencoba memperbaiki kode di bawah ini untuk melakukannya.

Saran berbasis permutasi asli (salah)

Saya sering menemukan bahwa seseorang dapat menyelamatkan diri dari banyak masalah dengan melakukan tes permutasi. Memang, dalam hal ini sangat mudah diatur. Mari ubah kontrol / kondisi eksperimental untuk setiap subjek secara terpisah; maka perbedaan dalam varian harus dihilangkan. Mengulangi ini berkali-kali akan menghasilkan distribusi nol untuk perbedaan.

(Saya tidak memprogram dalam R; semua orang jangan ragu untuk menulis ulang yang berikut dengan gaya R. yang lebih baik.)

set.seed(42)
nrep = 100
v = matrix(nrow=nrep, ncol=1)
for (i in 1:nrep)
{
   dp = d
   for (s in unique(d$participant_id)){             
     if (rbinom(1,1,.5)==1){
       dp[p$participant_id==s & d$condition=='control',]$condition = 'experimental'
       dp[p$participant_id==s & d$condition=='experimental',]$condition = 'control'
     }
   }
  m <- lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=dp)
  v[i,] = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]
}
pvalue = sum(abs(v) >= abs(delta)) / nrep

Menjalankan ini menghasilkan nilai-p $p=0.7$ . Seseorang dapat meningkat nrepmenjadi 1000 atau lebih.

Logika yang persis sama dapat diterapkan dalam kasus Bonus Anda.

amuba kata Reinstate Monica
sumber

Sangat menarik, terima kasih! Saya harus berpikir lebih banyak tentang mengapa reparameterisasi Anda berfungsi, karena ini tampaknya menjadi wawasan utama dari jawaban ini.

Patrick S. Forscher

Anehnya, nilai intersepsi per-grup dalam jawaban Anda tampaknya berbeda dari yang ada di jawaban @MartijnWeterings '.

Patrick S. Forscher

@ PatrickS.Forscher Itu karena dia, saya pikir, menghasilkan dataset yang berbeda. Saya dapat menggunakan sim_1 ~ 0 + condition + (0 + dummy(condition, "control") + dummy(condition, "experimental") | participant_id)formulasi dan mendapatkan hasil yang sama seperti pada jawaban saya.

Amoeba berkata Reinstate Monica

1

@ PatrickS.Forscher Tidak, saya menggunakan data yang dihasilkan oleh kode Anda (dengan seed Anda). Saya mengatur seed menjadi 42 hanya ketika melakukan pengujian permutasi. Martijn yang mengubah dataset, bukan saya.

Amoeba berkata Reinstate Monica

1

Proposal ini jelas masuk akal. Seperti yang saya pikir Anda sudah alami, menyiapkan tes permutasi untuk data bertingkat tidak sepenuhnya mudah. Pendekatan serupa yang akan sedikit lebih mudah untuk diterapkan adalah bootstrap parametrik, yang cukup sederhana untuk dilakukan dengan lme4 menggunakan metode simulate () dari objek lmer yang dipasang, yaitu, panggil simulasikan (m) berkali-kali untuk membangun bootstrap distribusi. Hanya ide untuk bermain-main.

Jake Westfall

0

Melihat masalah ini dari perspektif yang sedikit berbeda dan mulai dari bentuk "umum" dari model campuran linier, kita miliki

y_{i j k} = μ + α_{j} + d_{i j} + e_{i j k}, d_{i} \sim N (0, Σ), e_{i j k} \sim N (0, σ^{2})

$y_{ijk} = \mu + \alpha_j + d_{ij} + e_{ijk}, \quad d_{i} \sim N(0, \Sigma), \quad e_{ijk} \sim N(0, \sigma^2)$ where

α_{j}

$\alpha_j$ is the fixed effect of the

j

$j$ 'th condition and

d_{i} = (d_{i 1}, \dots, d_{i J})^{⊤}

$d_i = (d_{i1}, \ldots, d_{iJ})^\top$ is a random vector (some call it vector-valued random effect, I think) for the

i

$i$ 'th participant in the

j

$j$ 'th condition.
In your example we have two conditions

y_{i 1 k}

$y_{i1k}$ and

y_{i 2 k}

$y_{i2k}$ which I'll denote as

A

$A$ and

B

$B$ in what follows. So the covariance matrix of the two-dimensional random vector

d_{i}

$d_{i}$ is of the general form

Σ = [\begin{matrix} σ_{A}^{2} & σ_{A B} \\ σ_{A B} & σ_{B}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_A & \sigma_{AB}\\ \sigma_{AB} & \sigma^2_{B} \end{bmatrix}$

with non-negative $\sigma^2_A$ and $\sigma^2_{B}$ .

Let's first see how the re-parameterized version of $\Sigma$ looks when we use sum contrasts.

The variance of the intercept, which corresponds to the grand mean, is

σ_{1}^{2} := Var(grand mean) = Var (\frac{1}{2} \cdot (A + B)) = \frac{1}{4} \cdot (Var (A) + Var (B) + 2 \cdot Cov (A, B)) .

$\sigma^2_{1} := \text{Var(grand mean)} = \text{Var}(\frac{1}{2} \cdot (A+B)) = \frac{1}{4} \cdot (\text{Var}(A) + \text{Var}(B) + 2 \cdot \text{Cov}(A,B)).$

The variance of the contrast is

σ_{2}^{2} := Var(contrast) = Var (\frac{1}{2} \cdot (A - B)) = \frac{1}{4} \cdot (Var (A) + Var (B) - 2 \cdot Cov (A, B)) .

$\sigma^2_{2} := \text{Var(contrast)} = \text{Var}(\frac{1}{2} \cdot (A-B)) = \frac{1}{4} \cdot (\text{Var}(A) + \text{Var}(B) - 2 \cdot \text{Cov}(A,B)).$

And the covariance between the intercept and the contrast is

σ_{12} := Cov (grand mean, contrast) = Cov (\frac{1}{2} \cdot (A + B), \frac{1}{2} \cdot (A - B)) = \frac{1}{4} \cdot (Var (A) - Var (B)) .

$\sigma_{12} := \text{Cov}(\text{grand mean, contrast}) = \text{Cov}(\frac{1}{2} \cdot (A+B), \frac{1}{2} \cdot (A-B)) \\ = \frac{1}{4} \cdot (\text{Var}(A) - \text{Var}(B)).$

Thus, the re-parameterized $\Sigma$ is

Σ = [\begin{matrix} σ_{1}^{2} + σ_{2}^{2} + 2 σ_{12} & σ_{1}^{2} - σ_{2}^{2} \\ σ_{1}^{2} - σ_{2}^{2} & σ_{1}^{2} + σ_{2}^{2} - 2 σ_{12} \end{matrix}] = [\begin{matrix} σ_{A}^{2} & σ_{A B} \\ σ_{A B} & σ_{B}^{2} \end{matrix}] .

$\Sigma = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 + 2\sigma_{12} & \sigma^2_1 - \sigma^2_2\\ \sigma^2_1 - \sigma^2_2 & \sigma^2_1 + \sigma^2_2 - 2\sigma_{12} \end{bmatrix} = \begin{bmatrix}\sigma^2_A & \sigma_{AB}\\ \sigma_{AB} & \sigma^2_{B} \end{bmatrix}.$

$\Sigma$ can be decomposed into

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} σ_{2}^{2} & - σ_{2}^{2} \\ - σ_{2}^{2} & σ_{2}^{2} \end{matrix}] + 2 [\begin{matrix} σ_{12} & 0 \\ 0 & - σ_{12} \end{matrix}] .

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & -\sigma^2_2\\ -\sigma^2_2 & \sigma^2_2\end{bmatrix} + 2\begin{bmatrix}\sigma_{12} & 0\\ 0 & -\sigma_{12}\end{bmatrix}.$

Setting the covariance parameter $\sigma_{12}$ to zero we get

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} σ_{2}^{2} & - σ_{2}^{2} \\ - σ_{2}^{2} & σ_{2}^{2} \end{matrix}] = [\begin{matrix} σ_{1}^{2} + σ_{2}^{2} & σ_{1}^{2} - σ_{2}^{2} \\ σ_{1}^{2} - σ_{2}^{2} & σ_{1}^{2} + σ_{2}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & -\sigma^2_2\\ -\sigma^2_2 & \sigma^2_2\end{bmatrix} = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 & \sigma^2_1 - \sigma^2_2 \\ \sigma^2_1 - \sigma^2_2 & \sigma^2_1 + \sigma^2_2\end{bmatrix}$

which, as @Jake Westfall derived slightly differently, tests the hypothesis of equal variances when we compare a model without this covariance parameter to a model where the covariance parameter is still included/not set to zero.

Notably, introducing another crossed random grouping factor (such as stimuli) does not change the model comparison that has to be done, i.e., anova(mod1, mod2) (optionally with the argument refit = FALSE when you use REML estimation) where mod1 and mod2 are defined as @Jake Westfall did.

Taking out $\sigma_{12}$ and the variance component for the contrast $\sigma^2_2$ (what @Henrik suggests) results in

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix}$

which tests the hypothesis that the variances in the two conditions are equal and that they are equal to the (positive) covariance between the two conditions.

When we have two conditions, a model that fits a covariance matrix with two parameters in a (positive) compound symmetric structure can also be written as

# code snippet from Jake Westfall
d$contrast <- 2*(d$condition == 'experimental') - 1

# new model
mod3 <- lmer(sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id), 
             data = d, REML = FALSE)

or (using the categorical variable/factor condition)

mod4 <- lmer(sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id), 
             data = d, REML = FALSE)

with

Σ = [\begin{matrix} σ_{1}^{2} + σ_{2}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} + σ_{2}^{2} \end{matrix}] = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} σ_{2}^{2} & 0 \\ 0 & σ_{2}^{2} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1 + \sigma^2_2 \end{bmatrix} = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & 0\\ 0 & \sigma^2_2\end{bmatrix}$

where $\sigma^2_1$ and $\sigma^2_2$ are the variance parameters for the participant and the participant-condition-combination intercepts, respectively. Note that this $\Sigma$ has a non-negative covariance parameter.

Below we see that mod1, mod3, and mod4 yield equivalent fits:

# code snippet from Jake Westfall
d$contrast <- 2*(d$condition == 'experimental') - 1

mod1 <- lmer(sim_1 ~ contrast + (contrast || participant_id),
             data = d, REML = FALSE)

mod2 <- lmer(sim_1 ~ contrast + (contrast | participant_id),
             data = d, REML = FALSE)

# new models 
mod3 <- lmer(sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id), 
             data = d, REML = FALSE) 

mod4 <- lmer(sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id), 
             data = d, REML = FALSE)

anova(mod3, mod1)
# Data: d
# Models:
# mod3: sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id)
# mod1: sim_1 ~ contrast + ((1 | participant_id) + (0 + contrast | participant_id))
#      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
# mod3  5 2396.9 2420.3 -1193.5   2386.9                        
# mod1  5 2396.9 2420.3 -1193.5   2386.9     0      0          1

anova(mod4, mod3)
# Data: d
# Models:
# mod4: sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id)
# mod3: sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id)
#      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
# mod4  5 2396.9 2420.3 -1193.5   2386.9                        
# mod3  5 2396.9 2420.3 -1193.5   2386.9     0      0          1

With treatment contrasts (the default in R) the re-parameterized $\Sigma$ is

Σ = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} + σ_{12} \\ σ_{1}^{2} + σ_{12} & σ_{1}^{2} + σ_{2}^{2} + 2 σ_{12} \end{matrix}] = [\begin{matrix} σ_{1}^{2} & σ_{1}^{2} \\ σ_{1}^{2} & σ_{1}^{2} \end{matrix}] + [\begin{matrix} 0 & 0 \\ 0 & σ_{2}^{2} \end{matrix}] + [\begin{matrix} 0 & σ_{12} \\ σ_{12} & 2 σ_{12} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1 + \sigma_{12}\\ \sigma^2_1 + \sigma_{12} & \sigma^2_1 + \sigma^2_2 + 2\sigma_{12} \end{bmatrix} = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}0 & 0\\ 0 & \sigma^2_2\end{bmatrix} + \begin{bmatrix}0 & \sigma_{12}\\ \sigma_{12} & 2\sigma_{12}\end{bmatrix}$

where $\sigma^2_1$ is the variance parameter for the intercept (condition $A$ ), $\sigma^2_2$ the variance parameter for the contrast ( $A - B$ ), and $\sigma_{12}$ the corresponding covariance parameter.

We can see that neither setting $\sigma_{12}$ to zero nor setting $\sigma^2_2$ to zero tests (only) the hypothesis of equal variances.

However, as shown above, we can still use mod4 to test the hypothesis as changing the contrasts has no impact on the parameterization of $\Sigma$ for this model.

statmerkur
sumber

Model efek campuran: Bandingkan komponen varians acak di seluruh level variabel pengelompokan

Jawaban:

Jawaban singkat

Penjelasan / intuisi visual

Bukti