Bagaimana cara menghitung interval kepercayaan x-intersep dalam regresi linier?

9

Karena standard error dari regresi linier biasanya diberikan untuk variabel respon, saya bertanya-tanya bagaimana cara mendapatkan interval kepercayaan di arah lain - misalnya untuk x-intersep. Saya dapat memvisualisasikan apa yang mungkin terjadi, tetapi saya yakin pasti ada cara mudah untuk melakukan ini. Di bawah ini adalah contoh dalam R bagaimana memvisualisasikan ini:

set.seed(1)
x <- 1:10
a <- 20
b <- -2
y <- a + b*x + rnorm(length(x), mean=0, sd=1)

fit <- lm(y ~ x)
XINT <- -coef(fit)[1]/coef(fit)[2]

plot(y ~ x, xlim=c(0, XINT*1.1), ylim=c(-2,max(y)))
abline(h=0, lty=2, col=8); abline(fit, col=2)
points(XINT, 0, col=4, pch=4)
newdat <- data.frame(x=seq(-2,12,len=1000))

# CI
pred <- predict(fit, newdata=newdat, se.fit = TRUE) 
newdat$yplus <-pred$fit + 1.96*pred$se.fit 
newdat$yminus <-pred$fit - 1.96*pred$se.fit 
lines(yplus ~ x, newdat, col=2, lty=2)
lines(yminus ~ x, newdat, col=2, lty=2)

# approximate CI of XINT
lwr <- newdat$x[which.min((newdat$yminus-0)^2)]
upr <- newdat$x[which.min((newdat$yplus-0)^2)]
abline(v=c(lwr, upr), lty=3, col=4)

masukkan deskripsi gambar di sini

Marc di dalam kotak
sumber
1
Anda bisa bootstrap ini: library(boot); sims <- boot(data.frame(x, y), function(d, i) { fit <- lm(y ~ x, data = d[i,]) -coef(fit)[1]/coef(fit)[2] }, R = 1e4); points(quantile(sims$t, c(0.025, 0.975)), c(0, 0)). Untuk interval prediksi terbalik, file bantuan chemCal:::inverse.predictmemberikan referensi berikut yang mungkin juga membantu menurunkan CI: Massart, LM, Vandenginste, BGM, Buydens, LMC, De Jong, S., Lewi, PJ, Smeyers-Verbeke, J. (1997 ) Buku Pegangan Chemometrics dan Qualimetrics: Bagian A, hlm. 200
Roland
1
Apa yang Anda tunjukkan dalam grafik bukanlah CI untuk intersep. Anda menunjukkan titik-titik di mana garis keyakinan bawah dan atas dari prediksi melewati sumbu.
Roland
1
Seringkali dalam regresi linier seseorang memiliki model yang mengatakan sesuatu seperti ini: sehingga diperlakukan sebagai acak dan sebagai tetap. Itu dapat dibenarkan dengan mengatakan Anda mencari distribusi bersyarat yang diberikan . Dalam praktiknya jika Anda mengambil sampel baru, biasanya bukan hanya tetapi juga yang berubah, menunjukkan dalam beberapa keadaan mereka juga harus dianggap acak. Saya ingin tahu apakah ini sesuai dengan kepatutan
Yi=α+βxi+εiwhere ε1,εni.i.d. N(0,σ2),
YxxYx
Michael Hardy
1
@AdrienRenaud - Tampaknya bagi saya bahwa jawaban Anda terlalu sederhana mengingat aspek asimetris yang saya sebutkan, dan disorot oleh latihan bootstrap yang diilustrasikan Roland. Jika saya tidak bertanya terlalu banyak, mungkin Anda bisa memperluas pendekatan kemungkinan yang Anda sebutkan.
Marc di dalam kotak

Jawaban:

9

Bagaimana cara menghitung interval kepercayaan x-intersep dalam regresi linier?

Asumsi

  • Gunakan model regresi sederhana .yi=α+βxi+εi
  • Kesalahan memiliki distribusi normal yang tergantung pada regressorϵ|XN(0,σ2In)
  • Pas menggunakan kuadrat terkecil biasa

3 prosedur untuk menghitung interval kepercayaan pada x-intersep

Urutan pertama ekspansi Taylor

Model Anda dengan perkiraan standar deviasi dan pada dan parameter dan diperkirakan kovarians . Anda memecahkanY=aX+bσaσbabσab

aX+b=0X=ba.

Kemudian standar deviasi pada diberikan oleh:σXX

(σXX)2=(σbb)2+(σaa)22σabab.

MIB

Lihat kode dari Marc di kotak di Bagaimana cara menghitung interval kepercayaan x-intersep dalam regresi linier? .

CAPITANI-POLLASTRI

CAPITANI-POLLASTRI menyediakan Fungsi Distribusi Kumulatif dan Fungsi Kepadatan untuk rasio dua variabel acak Normal berkorelasi. Ini dapat digunakan untuk menghitung interval kepercayaan x-intersep dalam regresi linier. Prosedur ini memberikan (hampir) hasil yang identik dengan yang dari MIB.

Memang, menggunakan kuadrat terkecil biasa dan mengasumsikan normalitas kesalahan, (diverifikasi) dan berkorelasi (terverifikasi).β^N(β,σ2(XTX)1)β^

Prosedurnya adalah sebagai berikut:

  • dapatkan penaksir OLS untuk dan .ab
  • dapatkan matriks varians-kovarians dan ekstrak, .σa,σb,σab=ρσaσb
  • Asumsikan bahwa dan mengikuti distribusi Normal Berkorelasi Bivariat, . Kemudian fungsi kerapatan dan Fungsi Distribusi Kumulatif diberikan oleh CAPITANI-POLLASTRI.abN(a,b,σa,σb,ρ)xintercept=ba
  • Gunakan Fungsi Distribusi Kumulatif untuk menghitung kuantil yang diinginkan dan mengatur interval cofidence.xintercept=ba

Perbandingan 3 prosedur

Prosedur dibandingkan menggunakan konfigurasi data berikut:

  • x <- 1:10
  • a <- 20
  • b <- -2
  • y <- a + b * x + rnorm (panjang (x), rata-rata = 0, sd = 1)

10000 sampel berbeda dihasilkan dan dianalisis menggunakan 3 metode. Kode (R) yang digunakan untuk menghasilkan dan menganalisis dapat ditemukan di: https://github.com/adrienrenaud/stackExchange/blob/master/crossValidated/q221630/answer.ipynb

  • MIB dan CAPITANI-POLLASTRI memberikan hasil yang setara.
  • Urutan pertama ekspansi Taylor berbeda secara signifikan dari dua metode lainnya.
  • MIB dan CAPITANI-POLLASTRI mengalami kekurangan cakupan. 68% (95%) ci ditemukan mengandung nilai sebenarnya 63% (92%) saat itu.
  • Ekspansi Taylor urutan pertama mengalami over-coverage. 68% (95%) ci ditemukan mengandung nilai sebenarnya 87% (99%) saat itu.

Kesimpulan

Distribusi x-intersep bersifat asimetris. Ini membenarkan interval kepercayaan asimetris. MIB dan CAPITANI-POLLASTRI memberikan hasil yang setara. CAPITANI-POLLASTRI memiliki justifikasi teori yang bagus dan memberikan dasar bagi MIB. MIB dan CAPITANI-POLLASTRI menderita dari cakupan sedang dan dapat digunakan untuk mengatur interval kepercayaan.

Adrien Renaud
sumber
Terima kasih atas jawaban yang bagus ini. Apakah metode ini menyiratkan bahwa kesalahan standar dari intersep x simetris? Interval prediksi pada gambar saya menyiratkan bahwa ini bukan masalahnya, dan saya telah melihat referensi untuk ini di tempat lain.
Marc di dalam kotak
Ya, itu menyiratkan interval simetris. Jika Anda menginginkan yang asimetris, Anda dapat menggunakan kemungkinan profil memperlakukan parameter model Anda sebagai parameter gangguan. Tapi ini lebih banyak pekerjaan :)
Adrien Renaud
Bisakah Anda menjelaskan lebih detail bagaimana Anda mendapatkan ekspresi untuk ? (σX/X)2
@ fcop Ini adalah ekspansi Taylor. Silahkan lihat di en.wikipedia.org/wiki/Propagation_of_uncertainty
Adrien Renaud
2

Saya akan merekomendasikan bootstrap residunya:

library(boot)

set.seed(42)
sims <- boot(residuals(fit), function(r, i, d = data.frame(x, y), yhat = fitted(fit)) {

  d$y <- yhat + r[i]

  fitb <- lm(y ~ x, data = d)

  -coef(fitb)[1]/coef(fitb)[2]
}, R = 1e4)
lines(quantile(sims$t, c(0.025, 0.975)), c(0, 0), col = "blue")

plot yang dihasilkan

Apa yang Anda tunjukkan dalam grafik adalah titik-titik di mana batas bawah / atas dari pita kepercayaan prediksi melewati sumbu. Saya tidak berpikir ini adalah batas kepercayaan pencegatan, tapi mungkin mereka perkiraan kasar.

Roland
sumber
Hebat - ini sudah terlihat lebih masuk akal daripada contoh dari komentar Anda. Terima kasih lagi.
Marc di dalam kotak