Cara menentukan model efek campuran Bayesian dalam BUGS

Saya memposting ini di awal minggu lalu mencabut pertanyaan ketika saya menemukan sumber yang bagus, tidak ingin membuang waktu orang. Saya belum membuat banyak kemajuan, saya khawatir. Dalam mencoba menjadi warga negara yang baik di sini saya akan membuat masalah sejelas mungkin. Saya menduga akan ada beberapa peminat.

Saya memiliki kerangka data di RI yang ingin dianalisis dalam BUGS atau R. Ini dalam format panjang. Ini terdiri dari beberapa pengamatan pada 120 individu, dengan total 885 baris. Saya sedang memeriksa terjadinya hasil kategoris - tetapi itu tidak benar-benar relevan di sini. Pertanyaannya adalah tentang sesuatu yang lebih dalam.

Model yang saya gunakan sampai sekarang adalah

  mymodel<-gee(Category ~ Predictor 1 + Predictor 2..family=binomial(link="logit"),
  data=mydata, 
   id=Person)

dengan model marginal pada dasarnya akuntansi untuk pengelompokan pasien. Saya kemudian memeriksanya

 mymodel<-gee(Category ~ Predictor 1 + Predictor 2.. , family=binomial(link="logit"), 
  corstr = "AR-M", 
  data=mydata, id=Person)

untuk menjelaskan waktu pemesanan pengamatan pada masing-masing orang.

Ini tidak banyak berubah.

Kemudian saya mencoba memodelkan mereka menggunakan set perintah MCMCPack berikut:

 mymodel<-MCMCglmm(category~  Predictor1 + Predictor2..,  
 data=mydata, family=binomial(link="logit"))

Pemeriksaan output mendebarkan, menunjukkan signifikansi statistik untuk banyak prediktor. Saya memuji diri sendiri sebagai bayesian yang baru bertobat, sampai saya menyadari bahwa saya belum memperhitungkan tindakan berulang dalam pasien.

Saya mengerti bahwa saya harus bertanggung jawab untuk itu. Saya mengerti bahwa ini mungkin berarti menyesuaikan hyperprior untuk setiap individu - apakah itu benar? Bentuk apa yang akan diambil dalam BUGS?

Berikut model log reg dasar: (pujian untuk Kruschke, J., Indiana)

    model {
  for( i in 1 : nData ) {
    y[i] ~ dbern( mu[i] )
    mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))
  }
  b0 ~ dnorm( 0 , 1.0E-12 )
   for ( j in 1 : nPredictors ) {
    b[j] ~ dnorm( 0 , 1.0E-12 )
  }
}

Namun, tidak ada hyperprior di sini untuk individu. Inilah upaya terbaik saya sejauh ini pada desain dalam-individu, yang menghitung tindakan berulang dalam diri orang:

Inilah model Jackman untuk JAGS

1 model{
2 ## loop over data for likelihood
3 for(i in 1:n){
4  y[i] ~ dbern( mu[i] )
    mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))
6 }
7 sigma ˜ dunif(0,20) ## prior on standard deviation
8 tau <- pow(sigma,-2) ## convert to precision
9
10 ## hierarchical model for each state’s intercept & slope
11 for(p in 1:50){
12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal
13 }
14
15 ## means, hyper-parameters
16 for(q in 1:2){
17 mu[q] ˜ dnorm(0,.0016)

}

Inilah model bajingan-anak saya untuk BUGS

1 model{
2 ## loop over data for likelihood
3 for(i in 1:n){
4 mu.y[i] <- alpha + beta[s[i],1] + beta[s[i],2]*(j[i]-jbar)
5 demVote[i] ˜ dnorm(mu.y[i],tau)
6 }
7 sigma ˜ dunif(0,20) ## prior on standard deviation
8 tau <- pow(sigma,-2) ## convert to precision
9
10 ## hierarchical model for each state’s intercept & slope
11 for(p in 1:120){
12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal
13 }
14
15 ## means, hyper-parameters
16 for(q in 1:2){
17 mu[q] ˜ dnorm(0,.0016)
  }

Adakah yang bisa memberi tahu saya jika saya menuju ke arah yang benar. Pemahaman saya tentang ini berkembang, tetapi lambat. Mohon bersikap lembut. Saya seorang tenaga medis, bukan statistik! Saya telah menggunakan R cukup banyak, tetapi saya baru mengenal BUGS, dan baru untuk Bayes.

Terima kasih,

R

Anda hampir berada di sana. Hanya beberapa komentar - Anda tidak perlu menjadikan prior untuk beta[,1:2]parameter bersama MV normal; Anda dapat membuat sebelumnya sehingga beta[i,1]dan beta[i,2]independen, yang menyederhanakan hal-hal (misalnya, tidak ada kovarians sebelum perlu ditentukan.) Catatan bahwa melakukan hal itu tidak berarti mereka akan independen di posterior.

Komentar lain: Karena Anda memiliki istilah konstan - alpha- dalam regresi, komponen beta[,1]harus memiliki rata-rata nol di awal. Juga, Anda tidak memiliki prioritas alphadalam kode.

Berikut adalah model dengan istilah intersep dan kemiringan hierarkis; Saya sudah mencoba untuk menjaga prior dan notasi Anda jika memungkinkan, mengingat perubahan:

model {
  for(i in 1:n){
    mu.y[i] <- alpha + beta0[s[i]] + beta1[s[i]]*(j[i]-jbar)
    demVote[i] ~ dnorm(mu.y[i],tau)
  }

  alpha ~ dnorm(0, 0.001) ## prior on alpha; parameters just made up for illustration
  sigma ~ dunif(0,20) ## prior on standard deviation
  tau <- pow(sigma,-2) ## convert to precision

  ## hierarchical model for each state’s intercept & slope
  for (p in 1:120) {
     beta0[p] ~ dnorm(0, tau0)
     beta1[p] ~ dnorm(mu1, tau1)
  }

  ## Priors on hierarchical components; parameters just made up for illustration
  mu1 ~ dnorm(0, 0.001) 
  sigma0 ~ dunif(0,20)
  sigma1 ~ dunif(0,20)
  tau0 <- pow(sigma0,-2)
  tau1 <- pow(sigma1,-2)
}

Sumber daya yang sangat berguna untuk model hierarkis, termasuk beberapa "trik" untuk mempercepat konvergensi, adalah Gelman dan Hill .

(Agak terlambat dengan jawabannya, tetapi mungkin bermanfaat bagi beberapa penanya di masa depan.)