Mengapa kita harus menggunakan REML (bukan ML) untuk memilih di antara model var-covar bersarang?

Berbagai deskripsi tentang pemilihan model pada efek acak dari Linear Mixed Models memerintahkan untuk menggunakan REML. Saya tahu perbedaan antara REML dan ML pada tingkat tertentu, tetapi saya tidak mengerti mengapa REML harus digunakan karena ML bias. Misalnya, apakah salah menjalankan LRT pada parameter varians dari model distribusi normal menggunakan ML (lihat kode di bawah)? Saya tidak mengerti mengapa lebih penting untuk tidak memihak daripada menjadi ML, dalam pemilihan model. Saya pikir jawabannya adalah "karena pemilihan model bekerja lebih baik dengan REML daripada dengan ML" tetapi saya ingin tahu sedikit lebih banyak dari itu. Saya tidak membaca derivasi LRT dan AIC (saya tidak cukup baik untuk memahaminya secara menyeluruh), tetapi jika REML secara eksplisit digunakan dalam derivasi, hanya mengetahui bahwa itu akan cukup memadai (misalnya,

n <- 100
a <- 10
b <- 1
alpha <- 5
beta <- 1
x <- runif(n,0,10)
y <- rnorm(n,a+b*x,alpha+beta*x)

loglik1 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
  -sum(dnorm(y,a+b*x,alpha,log=T))
}

loglik2 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
   beta <- p[4]
  -sum(dnorm(y,a+b*x,alpha+beta*x,log=T))
}

m1 <- optim(c(a,b,alpha),loglik1,x=x,y=y)$value
m2 <- optim(c(a,b,alpha,beta),loglik2,x=x,y=y)$value
D <- 2*(m1-m2)
1-pchisq(D,df=1) # p-value

mixed-model maximum-likelihood unbiased-estimator likelihood-ratio reml berdalih
sumber

Tentang REML dan AIC, Anda harus melihat pertanyaan ini .

Elvis

Jawaban:

Jawaban yang sangat singkat: REML adalah ML, jadi tes yang didasarkan pada REML benar. Karena estimasi parameter varians dengan REML lebih baik, wajar untuk menggunakannya.

Mengapa REML ML? Pertimbangkan misalnya model dengan , , dan adalah vektor efek tetap, adalah vektor efek acak, dan . Kemungkinan Terbatas dapat diperoleh dengan mempertimbangkan kontras untuk "menghapus" efek tetap. Lebih tepatnya, biarkan , sedemikian sehingga dan (yaitu, kolom

Y = X β + Z u + e

$Y = X\beta + Zu + e \def\R{\mathbb{R}}$

X \in R^{n \times p}

$X\in\R^{n\times p}$

Z \in R^{n \times q}

$Z\in\R^{n\times q}$

β \in R^{p}

$\beta \in \R^p$

u \sim N (0, τ I_{q})

$u \sim \mathcal N(0, \tau I_q)$

e \sim N (0, σ^{2} I_{n})

$e \sim \mathcal N(0, \sigma^2 I_n)$

n - p

$n-p$

C \in R^{(n - p) \times n}

$C \in \R^{(n-p)\times n}$

C X = 0

$CX = 0$

C C^{'} = I_{n - p}

$CC' = I_{n-p}$

C^{'}

$C'$ adalah basis ortonormal dari ruang vektor ortognal ke ruang yang dihasilkan oleh kolom ); kemudian dengan , dan kemungkinan untuk diberikan adalah Kemungkinan Terbatas.

X

$X$

C Y = C Z u + ϵ

$CY = CZ u + \epsilon$

ϵ \sim N (0, σ^{2} I_{n - p})

$\epsilon \sim \mathcal N(0, \sigma^2 I_{n-p})$

τ, σ^{2}

$\tau, \sigma^2$

C Y

$CY$

Elvis
sumber

Jawaban yang bagus (+1), apakah saya benar mengatakan bahwa matriks tergantung pada model untuk rata-rata? Jadi Anda hanya bisa membandingkan perkiraan REML untuk matriks sama ?

C

$C$

C

$C$

Ya, tergantung pada (saya akan mengedit jawaban dalam satu menit untuk membuatnya jelas), jadi model bersarang Anda harus memiliki variabel yang sama dengan efek tetap.

C

$C$

X

$X$

Elvis

REML adalah tidak seorang ML! The ML didefinisikan unik untuk model probabilitas yang diberikan tetapi REML tergantung pada fixed-efek parameterisasi. Lihat misalnya komentar ini oleh Doug Bates (juga banyak yang bersejarah tentang model campuran R-SIG).

Livius

@Livius Saya pikir jawaban saya menyatakan dengan cukup jelas bagaimana kemungkinan terbatas dibangun. Ini adalah kemungkinan, hanya saja bukan kemungkinan yang diberikan

diamati dalam model yang ditulis dalam persamaan yang ditampilkan pertama, tetapi diberikan vektor yang diproyeksikan

dalam model yang ditulis dalam persamaan yang ditampilkan kedua. REML adalah ML yang diperoleh dari kemungkinan ini.

Y

$Y$

C Y

$CY$

Elvis

Saya pikir itu adalah pokok dari protes DBates mengenai masalah ini: ini adalah model yang berbeda, dan ini adalah model yang sulit untuk diperbandingkan karena model dan parameterisasi saling terkait. Jadi Anda tidak komputasi yang ML untuk model asli Anda, tetapi yang ML untuk model yang berbeda yang timbul dari suatu parameterisasi tertentu model asli Anda. Oleh karena itu model yang dipasang REML dengan struktur efek tetap bersarang bukan lagi model bersarang (seperti yang Anda sebutkan di atas). Tetapi model yang sesuai dengan ML masih bersarang, karena Anda memaksimalkan kemungkinan pada model yang ditentukan.

Livius

Tes rasio kemungkinan adalah uji hipotesis statistik yang didasarkan pada rasio dua kemungkinan. Properti mereka ditautkan ke estimasi kemungkinan maksimum (MLE). (lihat mis. Estimasi Kemungkinan Maksimum (MLE) dalam istilah awam ).

Dalam kasus Anda (lihat pertanyaan) Anda ingin '' memilih '' di antara dua model var-covar yang bersarang, katakanlah Anda ingin memilih antara model di mana var-covar adalah dan model di mana var-covar berada mana yang kedua (model sederhana) adalah kasus khusus dari yang pertama (yang umum). $\Sigma_g$ $\Sigma_s$

Tes ini didasarkan pada kemungkinan rasio . Dimana dan adalah penduga kemungkinan maksimum. $LR=-2 (log(\mathcal{L}_s(\hat{\Sigma}_s)) - log(\mathcal{L}_g(\hat{\Sigma}_g) )$ $\hat{\Sigma}_s$ $\hat{\Sigma}_g$

Statistik adalah, asimtotik (!) . $LR$ $\chi^2$

Estimator kemungkinan maksimum diketahui konsisten, namun, dalam banyak kasus mereka bias. Ini adalah kasus untuk estimator MLE untuk dan , dapat menunjukkan bahwa mereka bias. Ini karena mereka dihitung menggunakan mean yang diperoleh dari data, sedemikian rupa sehingga penyebaran di sekitar 'perkiraan rata-rata' ini lebih kecil daripada penyebaran di sekitar rata-rata sebenarnya (lihat misalnya Penjelasan intuitif untuk membagi dengan saat menghitung standar deviasi ? ) $\hat{\Sigma}_s$ $\hat{\Sigma}_g$ $n-1$

Statistik di atas adalah dalam sampel yang besar, ini hanya karena fakta bahwa, dalam sampel dan konvergen ke nilai-nilai mereka yang sebenarnya (MLE konsisten). (Catatan: di tautan di atas, untuk sampel yang sangat besar, membaginya dengan n atau dengan (n-1), tidak ada bedanya) $LR$ $\chi^2$ $\hat{\Sigma}_s$ $\hat{\Sigma}_g$

Untuk sampel yang lebih kecil, MLE memperkirakan dari dan akan menjadi bias dan karena itu distribusi akan menyimpang dari , sedangkan perkiraan REML akan memberikan perkiraan berisi untuk dan , jadi jika Anda menggunakan , untuk pemilihan model var-covar, estimasi REML maka untuk sampel yang lebih kecil akan lebih baik didekati dengan . $\hat{\Sigma}_s$ $\hat{\Sigma}_g$ $LR$ $\chi^2$ $\Sigma_s$ $\Sigma_g$ $LR$ $\chi^2$

Perhatikan bahwa REML hanya boleh digunakan untuk memilih di antara struktur bersarang model var-covar dengan rata-rata yang sama, untuk model dengan cara yang berbeda, REML tidak sesuai, untuk model dengan cara yang berbeda kita harus menggunakan ML.

sumber

Pernyataan "Statistik LR adalah, asimptotik (!) Χ2" tidak benar dalam kasus ini. Ini karena jika

bersarang di

, maka

berada di batas

. Dalam hal ini, distribusi

tidak berlaku. Misalnya, lihat di sini

Σ_{s}

$\Sigma_s$

Σ_{g}

$\Sigma_g$

Σ_{s}

$\Sigma_s$

Σ_{g}

$\Sigma_g$

χ^{2}

$\chi^2$

Cliff AB

@Cliff AB, inilah yang dijelaskan di bawah pernyataan itu dan itulah alasan Anda harus menggunakan REML.

-4

Saya punya jawaban yang lebih berkaitan dengan akal sehat daripada dengan Statistik. Jika Anda melihat PROC MIXED di SAS, estimasi dapat dilakukan dengan enam metode:

http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect008.htm

tapi REML adalah default. Mengapa? Rupanya, pengalaman praktis menunjukkan ia memiliki kinerja terbaik (misalnya, peluang terkecil dari masalah konvergensi). Oleh karena itu, jika tujuan Anda dapat dicapai dengan REML, maka masuk akal untuk menggunakan REML sebagai lawan dari lima metode lainnya.

James
sumber

Itu harus dengan 'teori sampel besar' dan bias estimasi MLE, lihat jawaban saya.

"Ini default di SAS" bukan jawaban yang dapat diterima untuk pertanyaan "mengapa" di situs ini.

Paul

p-nilai untuk model campuran yang disediakan oleh SAS secara default tidak tersedia dengan desain di perpustakaan lme4 untuk R karena tidak dapat dipercaya ( stat.ethz.ch/pipermail/r-help/2006-May/094765.html ). Jadi "SAS default" bahkan bisa salah.

Tim