Dalam model multi-level, apa implikasi praktis dari memperkirakan dan tidak mengestimasi parameter korelasi efek acak?

Dalam model multi-level, apa implikasi praktis dan terkait interpretasi parameter estimasi efek korelasi acak dan tidak diperkirakan? Alasan praktis untuk menanyakan hal ini adalah bahwa dalam kerangka lmer dalam R, tidak ada metode yang diimplementasikan untuk memperkirakan nilai-p melalui teknik MCMC ketika estimasi dibuat dalam model korelasi antara parameter.

Misalnya, dengan melihat contoh ini (bagian yang dikutip di bawah), apa implikasi praktis dari M2 versus M3. Jelas, dalam satu kasus P5 tidak akan diestimasi dan yang lain akan.

Pertanyaan

1. Untuk alasan praktis (keinginan untuk mendapatkan nilai-p melalui teknik MCMC) orang mungkin ingin menyesuaikan model tanpa korelasi antara efek acak bahkan jika P5 secara substansial tidak nol. Jika seseorang melakukan ini, dan kemudian memperkirakan nilai-p melalui teknik MCMC, apakah hasilnya dapat diartikan? (Saya tahu @Ben Bolker sebelumnya menyebutkan bahwa "menggabungkan pengujian signifikansi dengan MCMC sedikit tidak jelas, secara statistik, meskipun saya memahami keinginan untuk melakukannya (mendapatkan interval kepercayaan lebih mendukung)" , jadi jika itu akan membuat Anda tidur lebih baik pada malam hari berpura-puralah aku mengatakan interval kepercayaan diri.)
2. Jika seseorang gagal memperkirakan P5, apakah itu sama dengan menyatakan bahwa itu adalah 0?
3. Jika P5 benar-benar bukan nol, maka dengan cara apa nilai estimasi P1-P4 terpengaruh?
4. Jika P5 benar-benar bukan nol, maka dengan cara apa estimasi kesalahan untuk P1-P4 terpengaruh?
5. Jika P5 benar-benar bukan nol, maka dalam hal apa interpretasi model gagal memasukkan P5 cacat?

Meminjam dari jawaban @Mike Lawrence (mereka yang lebih berpengetahuan daripada saya bebas untuk mengganti ini dengan notasi model penuh, saya tidak sepenuhnya yakin saya bisa melakukannya dengan kesetiaan yang wajar):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2)(Perkiraan P1 - P4)

M3: V1 ~ (1+V3|V2) + V3(Perkiraan P1-P5)

Parameter yang mungkin diperkirakan:

P1 : Pencegatan global

P2 : Penyadapan efek acak untuk V2 (yaitu untuk setiap tingkat V2, penyimpangan penyadapan tingkat itu dari penyadapan global)

P3 : Perkiraan global tunggal untuk efek (kemiringan) V3

P4 : Efek V3 dalam setiap level V2 (lebih khusus lagi, sejauh mana efek V3 dalam level tertentu menyimpang dari efek global V3), sementara menegakkan korelasi nol antara penyimpangan intersep dan penyimpangan efek V3 lintas level. dari V2.

P5 : Korelasi antara deviasi intersep dan deviasi V3 lintas level V2

Jawaban yang berasal dari simulasi yang cukup besar dan luas bersama dengan kode yang menyertai dalam R menggunakan lmer akan diterima.

Pertimbangkan data sleepstudy, termasuk dalam lme4. Bates membahas hal ini dalam buku daringnya tentang lme4. Dalam bab 3, ia mempertimbangkan dua model untuk data.

dan

Penelitian ini melibatkan 18 subjek, dipelajari selama 10 hari kurang tidur. Waktu reaksi dihitung pada awal dan pada hari-hari berikutnya. Ada efek yang jelas antara waktu reaksi dan durasi kurang tidur. Ada juga perbedaan yang signifikan antara subjek. Model A memungkinkan untuk kemungkinan interaksi antara mencegat acak dan efek lereng: bayangkan, katakanlah, bahwa orang dengan waktu reaksi yang buruk menderita lebih akut dari efek kurang tidur. Ini akan menyiratkan korelasi positif dalam efek acak.

Dalam contoh Bates, tidak ada korelasi yang jelas dari plot Lattice dan tidak ada perbedaan yang signifikan antara model. Namun, untuk menyelidiki pertanyaan yang diajukan di atas, saya memutuskan untuk mengambil nilai-nilai yang sesuai dari sleepstudy, meningkatkan korelasi dan melihat kinerja kedua model.

Seperti yang dapat Anda lihat dari gambar, waktu reaksi lama dikaitkan dengan kehilangan kinerja yang lebih besar. Korelasi yang digunakan untuk simulasi adalah 0,58

Saya mensimulasikan 1000 sampel, menggunakan metode simulasi di lme4, berdasarkan nilai-nilai yang pas dari data buatan saya. Saya cocok M0 dan Ma untuk masing-masing dan melihat hasilnya. Set data asli memiliki 180 pengamatan (10 untuk masing-masing 18 mata pelajaran), dan data yang disimulasikan memiliki struktur yang sama.

Intinya adalah bahwa ada sedikit perbedaan.

1. Parameter tetap memiliki nilai yang persis sama di bawah kedua model.
2. Efek acak sedikit berbeda. Ada 18 efek mencegat dan 18 kemiringan acak untuk setiap sampel yang disimulasikan. Untuk setiap sampel, efek-efek ini dipaksa untuk ditambahkan ke 0, yang berarti bahwa perbedaan rata-rata antara kedua model adalah (secara artifisial) 0. Tetapi varians dan kovarian berbeda. Median kovarians di bawah MA adalah 104, melawan 84 di bawah M0 (nilai aktual, 112). Varian dari lereng dan penyadapan lebih besar di bawah M0 dari MA, mungkin untuk mendapatkan ruang gerak tambahan yang dibutuhkan tanpa adanya parameter kovarian gratis.
3. Metode ANOVA untuk lmer memberikan statistik F untuk membandingkan model Slope dengan model dengan hanya intersep acak (tidak ada efek karena kurang tidur). Jelas, nilai ini sangat besar di bawah kedua model, tetapi biasanya (tetapi tidak selalu) lebih besar di bawah MA (rata-rata 62 vs rata-rata 55).
4. Perbedaan kovarians dan efek tetap berbeda.
5. Sekitar separuh waktu, ia tahu bahwa MA benar. Nilai p rata-rata untuk membandingkan M0 dengan MA adalah 0,0442. Meskipun terdapat korelasi yang berarti dan 180 pengamatan seimbang, model yang tepat hanya akan dipilih sekitar separuh waktu.
6. Nilai yang diprediksi berbeda di bawah kedua model, tetapi sangat sedikit. Perbedaan rata-rata antara prediksi adalah 0, dengan sd dari 2.7. Sd dari nilai yang diprediksi sendiri adalah 60,9

Jadi mengapa ini terjadi? @ung menduga, cukup, bahwa kegagalan untuk memasukkan kemungkinan korelasi memaksa efek acak menjadi tidak berkorelasi. Mungkin seharusnya; tetapi dalam implementasi ini, efek acak diizinkan untuk dikorelasikan, yang berarti bahwa data dapat menarik parameter ke arah yang benar, terlepas dari modelnya. Kesalahan dari model yang salah muncul dalam kemungkinan, itulah sebabnya Anda dapat (kadang-kadang) membedakan dua model di tingkat itu. Model efek campuran pada dasarnya sesuai dengan regresi linier untuk setiap subjek, dipengaruhi oleh apa yang menurut model seharusnya. Model yang salah memaksa kecocokan nilai yang kurang masuk akal daripada yang Anda dapatkan di bawah model yang tepat. Tetapi parameter, pada akhirnya, diatur oleh kecocokan dengan data aktual.

Ini kode saya yang agak kikuk. Idenya adalah agar sesuai dengan data studi tidur dan kemudian membangun set data simulasi dengan parameter yang sama, tetapi korelasi yang lebih besar untuk efek acak. Kumpulan data diumpankan ke simulasi.lmer () untuk mensimulasikan 1000 sampel, yang masing-masing sesuai dengan kedua cara. Setelah saya memasangkan objek yang pas, saya bisa mengeluarkan fitur yang berbeda dari fit dan membandingkannya, menggunakan uji-t, atau apa pun.

    # Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2) 
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01) 
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study. 
yhat <- XX.f %*%  fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10) 
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
    sim.data <- vector(mode="list",)
    tempReaction <- simulate(fm.sim, nsim=n.sim)
    tempdata <- model.frame(fm.sim)
    for (i in 1:n.sim){
        tempdata$Reaction <- tempReaction[,i]
			output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
			output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
			temp <- anova(output0,output1)
			pval <- temp$`Pr(>Chisq)`[2]
        sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
    }

Placidia telah memberikan jawaban menyeluruh menggunakan data simulasi berdasarkan pada sleepstudydataset. Berikut adalah jawaban lain (kurang ketat) yang juga menggunakan sleepstudydata.

Kita melihat bahwa seseorang dapat mempengaruhi perkiraan korelasi antara intersepsi acak dan kemiringan acak dengan "menggeser" variabel prediktor acak. Lihatlah hasil dari model fm1dan di fm2bawah ini:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

Dari output model, kita melihat bahwa korelasi varians acak telah berubah. Namun, lereng (tetap dan acak) tetap sama, seperti yang dilakukan estimasi varians residual. Estimasi intersep (tetap dan acak) berubah sebagai respons terhadap variabel yang bergeser.

De-berkorelasi kovariansi mencegat-lereng acak untuk LMM dibahas dalam catatan kuliah Dr. Jack Weiss di sini . Weiss mencatat bahwa mengurangi korelasi varians dalam mode ini kadang-kadang dapat membantu dengan konvergensi model, antara lain.

Contoh di atas memvariasikan korelasi acak (parameter "P5"). Mengatasi sebagian OP Q3, kita melihat dari output di atas bahwa:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected