Regresi umum tertimbang dalam BUGS, JAGS

10

Di Rkita dapat "bobot sebelumnya" glmregresi melalui parameter bobot . Sebagai contoh:

glm.D93 <- glm(counts ~ outcome + treatment, family = poisson(), weights=w)

Bagaimana ini bisa dicapai dalam model JAGSatau BUGS?

Saya menemukan beberapa makalah yang membahas hal ini, tetapi tidak ada yang memberikan contoh. Saya tertarik terutama ke dalam contoh-contoh Poisson dan regresi logistik.

pengguna28937
sumber
+1 pertanyaan yang sangat bagus! Saya menanyakan ini kepada beberapa spesialis bayesian dan mereka hanya mengatakan bahwa dalam beberapa kasus (bobot menurut kovariat kategori), Anda dapat menghitung distribusi posterior parameter (s) untuk setiap kategori dan kemudian menggabungkannya dalam rata-rata tertimbang. Mereka tidak memberi saya solusi umum, jadi saya akan sangat tertarik jika ada atau tidak!
Penasaran

Jawaban:

7

Mungkin terlambat ... tapi,

Harap dicatat 2 hal:

  • Menambahkan poin data tidak disarankan karena akan mengubah derajat kebebasan. Perkiraan efek tetap dapat diestimasi dengan baik, tetapi semua kesimpulan harus dihindari dengan model seperti itu. Sulit untuk "membiarkan data berbicara" jika Anda mengubahnya.
  • Tentu saja itu hanya bekerja dengan bobot bilangan bulat (Anda tidak dapat menduplikasi 0,5 titik data), yang bukan apa yang dilakukan dalam kebanyakan regresi berbobot (lm). Secara umum, penimbangan dibuat berdasarkan variabilitas lokal yang diperkirakan dari ulangan (mis. 1 / d atau 1 / d ^ 2 pada 'x' yang diberikan) atau berdasarkan tinggi respons (mis. 1 / Y atau 1 / Y ^ 2, pada diberikan 'x').

Di Jags, Bugs, Stan, proc MCMC, atau di Bayesian secara umum, kemungkinannya tidak berbeda dari pada frequentm lm atau glm (atau model apa pun), itu sama saja !! Cukup buat "bobot" kolom baru untuk respons Anda, dan tulis kemungkinannya sebagai

y [i] ~ dnorm (mu [i], tau / weight [i])

Atau poisson tertimbang:

y [i] ~ dpois (lambda [i] * weight [i])

Kode Bugs / Jags ini hanya untuk trik. Anda akan mendapatkan semuanya dengan benar. Jangan lupa untuk terus mengalikan posterior tau dengan bobot, misalnya saat membuat prediksi dan interval kepercayaan / prediksi.

Pierre Lebrun
sumber
Jika kita melakukannya seperti yang dinyatakan, kita mengubah mean dan varians. Mengapa bukan y [i] * berat [i] ~ dpois (lambda [i] * berat [i])? Itu hanya akan menyesuaikan varians. Masalahnya di sini adalah y [i] * berat [i] mungkin bertipe real.
user28937
memang, regresi berbobot memang berubah berarti (karena penimbangan mengarahkan regresi untuk lebih dekat ke titik-titik yang memiliki banyak bobot!) dan varians sekarang merupakan fungsi dari bobot (karenanya bukan model homoskedastik). Varians (atau presisi) tau tidak memiliki makna lagi, tetapi tau / berat [i] dapat diartikan persis seperti ketepatan model (untuk "x") yang diberikan. Saya tidak akan menyarankan penggandaan data (y) dengan bobot ... berharap jika ini adalah sesuatu yang ingin Anda lakukan, tapi saya tidak mengerti model Anda dalam hal ini ...
Pierre Lebrun
Saya setuju dengan Anda itu tidak mengubah rata-rata dalam contoh normal: y [i] ~ dnorm (mu [i], tau / weight [i]), tetapi itu terjadi pada yang kedua, karena lambda [i] * weight [ i] menjadi lambda "baru" untuk dpois dan ini tidak akan cocok dengan y [i] lagi. Saya harus memperbaiki diri sendiri seharusnya: ty [i] * exp (berat [i]) ~ dpois (lambda [i] * weight [i]). Gagasan dengan perkalian dalam kasus poisson adalah bahwa kita ingin menyesuaikan varians, tetapi juga menyesuaikan mean, jadi bukankah kita perlu mengoreksi mean?
user28937
Jika Anda perlu menyesuaikan varians secara independen, mungkin model Binomial Negatif mungkin berguna daripada Poisson? Itu menambahkan varians inflasi / parameter deflasi ke Poisson. Kecuali itu sangat mirip.
Pierre Lebrun
Pierre ide bagus. Saya juga memikirkan tentang representasi berkelanjutan dari distribusi poisson yang didefinisikan dalam slide 6/12 di linkd
user28937
4

Pertama, ada baiknya menunjukkan bahwa glmtidak melakukan regresi bayesian. Parameter 'bobot' pada dasarnya adalah kependekan dari "proporsi pengamatan," yang dapat diganti dengan pengambilan sampel dataset Anda secara tepat. Sebagai contoh:

x=1:10
y=jitter(10*x)
w=sample(x,10)

augmented.x=NULL
augmented.y=NULL    
for(i in 1:length(x)){
    augmented.x=c(augmented.x, rep(x[i],w[i]))
    augmented.y=c(augmented.y, rep(y[i],w[i]))
}

# These are both basically the same thing
m.1=lm(y~x, weights=w)
m.2=lm(augmented.y~augmented.x)

Jadi, untuk menambah bobot pada poin dalam JAGS atau BUGS, Anda dapat menambah dataset Anda dengan cara yang sama seperti di atas.

David Marx
sumber
2
ini tidak akan berhasil, bobot adalah bilangan real yang biasa, bukan bilangan bulat
Curious
Ini tidak menghalangi Anda untuk mendekati mereka dengan bilangan bulat. Solusi saya tidak sempurna, tetapi kira-kira berfungsi. Misalnya, dengan bobot (1/3, 2/3, 1), Anda bisa meng-upample kelas dua dengan faktor dua dan kelas ketiga dengan faktor tiga.
David Marx
0

Sudah mencoba menambahkan komentar di atas, tetapi perwakilan saya terlalu rendah.

Sebaiknya

y[i] ~ dnorm(mu[i], tau / weight[i])

tidak mungkin

y[i] ~ dnorm(mu[i], tau * weight[i])

dalam JAGS? Saya menjalankan beberapa tes yang membandingkan hasil dari metode ini di JAGS dengan hasil dari regresi tertimbang via lm () dan hanya dapat menemukan yang sesuai dengan yang terakhir. Berikut ini contoh sederhana:

aggregated <- 
  data.frame(x=1:5) %>%
  mutate( y = round(2 * x + 2 + rnorm(length(x)) ),
          freq = as.numeric(table(sample(1:5, 100, 
                 replace=TRUE, prob=c(.3, .4, .5, .4, .3)))))
x <- aggregated$x
y <- aggregated$y
weight <- aggregated$freq
N <- length(y)

# via lm()
lm(y ~ x, data = aggregated, weight = freq)

dan dibandingkan dengan

lin_wt_mod <- function() {

  for (i in 1:N) {
    y[i] ~ dnorm(mu[i], tau*weight[i])
    mu[i] <- beta[1] + beta[2] * x[i]
  }

  for(j in 1:2){
    beta[j] ~ dnorm(0,0.0001)
  }

  tau   ~ dgamma(0.001, 0.001)
  sigma     <- 1/sqrt(tau)
}

dat <- list("N","x","y","weight")
params <- c("beta","tau","sigma")

library(R2jags)
fit_wt_lm1 <- jags.parallel(data = dat, parameters.to.save = params,
              model.file = lin_wt_mod, n.iter = 3000, n.burnin = 1000)
fit_wt_lm1$BUGSoutput$summary
metrikrn
sumber
Terlepas dari reputasi, komentar tidak boleh diberikan sebagai jawaban.
Michael R. Chernick