Interval kepercayaan untuk parameter regresi: Bayesian vs klasik

Diberikan dua array x dan y, keduanya panjang n, saya cocok dengan model y = a + b * x dan ingin menghitung interval kepercayaan 95% untuk lereng. Ini adalah (b - delta, b + delta) di mana b ditemukan dengan cara biasa dan

delta = qt(0.975,df=n-2)*se.slope

dan se.slope adalah kesalahan standar pada lereng. Salah satu cara untuk mendapatkan standard error dari lereng dari R adalah summary(lm(y~x))$coef[2,2].

Sekarang anggaplah saya menulis kemungkinan kemiringan yang diberikan x dan y, kalikan ini dengan "flat" sebelumnya dan gunakan teknik MCMC untuk menggambar sampel m dari distribusi posterior. Menetapkan

lims = quantile(m,c(0.025,0.975))

Pertanyaan saya: (lims[[2]]-lims[[1]])/2kira-kira sama dengan delta seperti yang didefinisikan di atas?

Tambahan Di bawah ini adalah model JAGS sederhana di mana keduanya tampak berbeda.

model {
 for (i in 1:N) {
  y[i] ~ dnorm(mu[i], tau)
  mu[i] <- a + b * x[i]
 }
 a ~ dnorm(0, .00001)
 b ~ dnorm(0, .00001)
 tau <- pow(sigma, -2)
 sigma ~ dunif(0, 100)
}

Saya menjalankan yang berikut ini di R:

N <- 10
x <- 1:10
y <- c(30.5,40.6,20.5,59.1,52.5,
       96.0,121.4,78.9,112.1,128.4)
lin <- lm(y~x)

#Calculate delta for a 95% confidence interval on the slope
delta.lm <- qt(0.975,df=N-2)*summary(lin)$coef[2,2]

library('rjags')
jags <- jags.model('example.bug', data = list('x' = x,'y' = y,'N' = N),
                   n.chains = 4,n.adapt = 100)
update(jags, 1000)
params <- jags.samples(jags,c('a', 'b', 'sigma'),7500)
lims <- quantile(params$b,c(0.025,0.975))
delta.bayes <- (lims[[2]]-lims[[1]])/2

cat("Classical confidence region: +/-",round(delta.lm, digits=4),"\n")
cat("Bayesian confidence region:  +/-",round(delta.bayes,digits=4),"\n")

Dan dapatkan:

Wilayah kepercayaan klasik: +/- 4.6939

Wilayah kepercayaan Bayesian: +/- 5.1605

Menjalankan kembali ini beberapa kali, wilayah kepercayaan Bayesian secara konsisten lebih luas daripada yang klasik. Jadi apakah ini karena prior yang saya pilih?

Masalahnya ada pada prior pada sigma. Coba pengaturan yang kurang informatif

tau ~ dgamma(1.0E-3,1.0E-3)
sigma <- pow(tau, -1/2)

dalam file jags Anda. Kemudian perbarui a bunch

update(10000)

ambil parameternya, dan rangkum jumlah minat Anda. Seharusnya cukup selaras dengan versi klasik.

Klarifikasi: Pembaruan hanya untuk memastikan Anda mendapatkan ke mana Anda pergi dengan pilihan apa pun sebelum Anda memutuskan, meskipun rantai untuk model seperti ini dengan prius yang tersebar dan nilai awal acak memang membutuhkan waktu lebih lama untuk bertemu. Dalam masalah nyata Anda akan memeriksa konvergensi sebelum merangkum apa pun, tetapi konvergensi bukanlah masalah utama dalam contoh Anda, saya tidak berpikir.

Jika Anda mengambil sampel dari posterior b | y dan menghitung lims (seperti yang Anda tentukan) itu harus sama dengan (b - delta, b + delta). Khususnya, jika Anda menghitung distribusi posterior b | y di bawah flat sebelumnya, sama dengan distribusi sampling klasik b.

Untuk detail lebih lanjut lihat: Gelman et al. (2003). Analisis Data Bayesian. CRC Tekan. Bagian 3.6

Edit:

Ringold, perilaku yang Anda amati konsisten dengan ide Bayesian. Interval Kredibel Bayesian (CI) umumnya lebih lebar daripada yang klasik. Dan alasannya adalah, seperti yang Anda tebak dengan tepat, hyperpriors memperhitungkan variabilitas karena parameter yang tidak diketahui.

Untuk skenario sederhana seperti ini (BUKAN DALAM UMUM):

Baysian CI> Empiris Bayesian CI> CI Klasik; > == lebih luas

Untuk model Gaussian linier, lebih baik menggunakan paket bayesm. Ini mengimplementasikan keluarga semi-konjugat dari prior, dan sebelumnya Jeffreys adalah batas kasus keluarga ini. Lihat contoh saya di bawah ini. Ini adalah simulasi klasik, tidak perlu menggunakan MCMC.

Saya tidak ingat apakah interval kredibilitas tentang parameter regresi persis sama dengan interval kepercayaan kuadrat biasa, tetapi dalam kasus apa pun mereka sangat dekat.

> # required package
> library(bayesm)
> # data
> age <- c(35,45,55,65,75)
> tension <- c(114,124,143,158,166)
> y <- tension
> # model matrix
> X <- model.matrix(tension~age)
> # prior parameters
> Theta0 <- c(0,0)
> A0 <- 0.0001*diag(2)
> nu0 <- 0
> sigam0sq <- 0
> # number of simulations
> n.sims <- 5000
> # run posterior simulations
> Data <- list(y=y,X=X)
> Prior <- list(betabar=Theta0, A=A0, nu=nu0, ssq=sigam0sq)
> Mcmc <- list(R=n.sims)
> bayesian.reg <- runireg(Data, Prior, Mcmc)
> beta.sims <- t(bayesian.reg$betadraw) # transpose of bayesian.reg$betadraw
> sigmasq.sims <- bayesian.reg$sigmasqdraw
> apply(beta.sims, 1, quantile, probs = c(0.025, 0.975))
[,1] [,2]
2.5% 53.33948 1.170794
97.5% 77.23371 1.585798
> # to be compared with: 
> frequentist.reg <- lm(tension~age)

Mengingat bahwa regresi linier sederhana secara analitik identik antara analisis klasik dan Bayesian dengan Jeffrey sebelumnya, keduanya analitik, tampaknya agak aneh untuk menggunakan metode numerik seperti MCMC untuk melakukan analisis Bayesian. MCMC hanyalah alat integrasi numerik, yang memungkinkan metode Bayesian untuk digunakan dalam masalah yang lebih rumit yang secara analitis tidak dapat dipecahkan, sama seperti Newton-Rhapson atau Fisher Scoring adalah metode numerik untuk menyelesaikan masalah klasik yang tidak dapat dipecahkan.

Distribusi posterior p (b | y) menggunakan p Jeffrey sebelumnya (a, b, s) sebanding dengan 1 / s (di mana s adalah standar deviasi kesalahan) adalah distribusi t siswa dengan lokasi b_ols, skala se_b_ols (" ols "untuk" kuadrat terkecil biasa "), dan n-2 derajat kebebasan. Tetapi distribusi sampling b_ols juga merupakan t siswa dengan lokasi b, skala se_b_ols, dan n-2 derajat kebebasan. Jadi mereka identik kecuali bahwa b dan b_ols telah ditukar, jadi ketika datang untuk menciptakan interval, "est + - bound" interval kepercayaan akan dibalik menjadi "est - + terikat" dalam interval kredibel.

Jadi interval kepercayaan dan interval kredibel identik secara analitis, dan tidak masalah metode mana yang digunakan (asalkan tidak ada informasi tambahan sebelumnya) - jadi ambil metode yang lebih murah secara komputasi (mis. Yang memiliki inversi matriks lebih sedikit). Apa hasil Anda dengan MCMC menunjukkan bahwa pendekatan tertentu yang digunakan dengan MCMC memberikan interval kredibel yang terlalu lebar dibandingkan dengan interval kredibel analitik yang tepat. Ini mungkin hal yang baik (meskipun kami ingin perkiraannya menjadi lebih baik) bahwa solusi Bayesian yang diperkirakan muncul lebih konservatif daripada solusi Bayesian yang tepat.