Bagaimana saya bisa memodelkan proporsi dengan BUGS / JAGS / STAN?

10

Saya mencoba membangun sebuah model di mana responsnya adalah proporsi (sebenarnya itu adalah bagian dari suara yang didapat partai di daerah pemilihan). Distribusi tidak normal, jadi saya memutuskan untuk memodelkannya dengan distribusi beta. Saya juga punya beberapa prediktor.

Namun, saya tidak tahu bagaimana menulisnya dalam BUGS / JAGS / STAN (JAGS akan menjadi pilihan terbaik saya, tetapi itu tidak terlalu penting). Masalah saya adalah saya membuat sejumlah parameter oleh prediktor, tapi lalu apa yang bisa saya lakukan dengannya?

Kode akan menjadi seperti ini (dalam sintaks JAGS), tapi saya tidak tahu bagaimana "menghubungkan" parameter y_hatdan y.

for (i in 1:n) {
 y[i] ~ dbeta(alpha, beta)

 y_hat[i] <- a + b * x[i]
}

( y_hathanya produk-silang dari parameter dan prediktor, maka hubungan deterministik. adan bmerupakan koefisien yang saya coba perkirakan, xmenjadi prediktor).

Terima kasih atas saran Anda!

Joël
sumber
Apa itu, b, y_hat? Anda harus mendefinisikan model Anda dengan jelas. Omong-omong sintaks BUGS dekat dengan sintaksis matematika. Jadi, jika Anda tahu cara menulis model Anda dalam bahasa matematika maka hampir semua pekerjaan dilakukan.
Stéphane Laurent
Stéphane, terima kasih. Saya mengedit pertanyaan untuk mendefinisikan a, b, y_hat. Saya juga tidak tahu jawabannya secara matematis, kalau tidak jawabannya akan jauh lebih mudah ;-)
Joël
Saya curiga saya bisa membangun fakta bahwa E (y) = alpha / (alpha + beta), tapi saya tidak bisa benar-benar mencari tahu bagaimana tepatnya.
Joël

Jawaban:

19

Pendekatan regresi beta adalah reparameterisasi dalam hal dan ϕ . Di mana μ akan setara dengan apa yang Anda prediksi. Dalam parameterisasi ini Anda akan memiliki α = μ × ϕ dan β = ( 1 - μ ) × ϕ . Kemudian Anda dapat memodelkan μ sebagai logit dari kombinasi linear. ϕ dapat memiliki sebelumnya sendiri (harus lebih besar dari 0), atau dapat dimodelkan pada kovariat juga (pilih fungsi tautan untuk tetap lebih besar dari 0, seperti eksponensial).μϕμα=μ×ϕβ=(1μ)×ϕμϕ

Mungkin sesuatu seperti:

for(i in 1:n) {
  y[i] ~ dbeta(alpha[i], beta[i])
  alpha[i] <- mu[i] * phi
  beta[i]  <- (1-mu[i]) * phi
  logit(mu[i]) <- a + b*x[i]
}
phi ~ dgamma(.1,.1)
a ~ dnorm(0,.001)
b ~ dnorm(0,.001)
Greg Snow
sumber
Terima kasih, ini sangat membantu! Saya mencoba menyesuaikan model dengan saran Anda.
Joël
Namun, ketika saya menjalankan model, saya mendapatkan kesalahan seperti: "Kesalahan dalam simpul y [6283] Nilai induk tidak valid". Adakah yang tahu apa yang terjadi di sini?
Joël
@ Joël, berapa nilai y [6283]? Sudahkah Anda memastikan bahwa nilai-nilai alfa dan beta dibatasi untuk nilai legal? Saya berharap bahwa sesuatu mungkin sudah ke 0 atau di bawah dan itu memberikan kesalahan.
Greg Snow
Tidak, saya memeriksa itu, semua nilai y saya benar-benar lebih unggul dari 0 (dan lebih rendah dari 1). Mungkin priors saya berbenturan dengan nilai-nilai empiris y di beberapa titik? Tapi saya tidak tahu bagaimana memeriksanya, dan masalah saya tampaknya masuk akal - setidaknya bagi saya!
Joël
1
@colin, saya tidak tahu JAGS dengan baik, jadi ini mungkin lebih baik ditanyakan di forum khusus untuk JAGS. Atau coba di alat yang berbeda, saya menemukan bahwa saya suka Stan for Bayes hari ini.
Greg Snow
18

Greg Snow memberikan jawaban yang bagus. Untuk kelengkapan, berikut ini ekuivalen dalam sintaks Stan. Meskipun Stan memiliki distribusi beta yang dapat Anda gunakan, lebih cepat untuk mengerjakan sendiri logaritma kepadatan beta karena konstanta log(y)dan log(1-y)dapat dihitung sekali pada permulaan (daripada setiap kali y ~ beta(alpha,beta)akan dipanggil). Dengan menambah lp__variabel yang dipesan (lihat di bawah), Anda dapat menjumlahkan logaritma densitas beta di atas pengamatan dalam sampel Anda. Saya menggunakan label "gamma" untuk vektor parameter dalam prediktor linier.

data {
  int<lower=1> N;
  int<lower=1> K;
  real<lower=0,upper=1> y[N];
  matrix[N,K] X;
}
transformed data {
  real log_y[N];
  real log_1my[N];
  for (i in 1:N) {
    log_y[i] <- log(y[i]);
    log_1my[i] <- log1m(y[i]);
  }
}
parameters {
  vector[K] gamma;
  real<lower=0> phi;
}
model {
  vector[N] Xgamma;
  real mu;
  real alpha_m1;
  real beta_m1;
  Xgamma <- X * gamma;
  for (i in 1:N) {
    mu <- inv_logit(Xgamma[i]);
    alpha_m1 <- mu * phi - 1.0;
    beta_m1 <- (1.0 - mu) * phi - 1.0;
    lp__ <- lp__ - lbeta(alpha,beta) + alpha_m1 * log_y[i] + 
                                        beta_m1 * log_1my[i];
  }
  // optional priors on gamma and phi here
}
Ben Goodrich
sumber
Ben terima kasih! Sangat berguna untuk memiliki sintaks Stan juga.
Joël
Stan v2 memiliki pernyataan pengambilan sampel "beta_proportion" yang saya yakin meniadakan perlunya memanipulasi secara langsung "lp__"
THK