Menghitung interval kepercayaan untuk regresi logistik

15

Saya menggunakan regresi logistik binomial untuk mengidentifikasi apakah paparan has_xatau has_ydampak kemungkinan bahwa pengguna akan mengklik sesuatu. Model saya adalah sebagai berikut:

fit = glm(formula = has_clicked ~ has_x + has_y, 
          data=df, 
          family = binomial())

Ini output dari model saya:

Call:
glm(formula = has_clicked ~ has_x + has_y, 
    family = binomial(), data = active_domains)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9869  -0.9719  -0.9500   1.3979   1.4233  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -0.504737   0.008847 -57.050  < 2e-16 ***
has_xTRUE -0.056986   0.010201  -5.586 2.32e-08 ***
has_yTRUE  0.038579   0.010202   3.781 0.000156 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 217119  on 164182  degrees of freedom
Residual deviance: 217074  on 164180  degrees of freedom
AIC: 217080

Number of Fisher Scoring iterations: 4

Karena masing-masing koefisien signifikan, dengan menggunakan model ini saya dapat memberi tahu apa nilai kombinasi ini menggunakan pendekatan berikut:

predict(fit, data.frame(has_x = T, has_y=T), type = "response")

Saya tidak mengerti bagaimana saya bisa melaporkan Std. Kesalahan prediksi.

  1. Apakah saya hanya perlu menggunakan ? Atau apakah saya perlu mengonversi menggunakan pendekatan yang dijelaskan di sini ?1.96SESE

  2. Jika saya ingin memahami kesalahan standar untuk kedua variabel, bagaimana saya mempertimbangkannya?

Tidak seperti pertanyaan ini , saya tertarik untuk memahami apa batas kesalahan atas dan bawah dalam persentase. Misalnya, prediksi saya menunjukkan nilai 37% untuk True,Truedapatkah saya menghitung bahwa ini untuk ? (0,3% dipilih untuk menggambarkan poin saya)95 % C I+/0.395%CI

celenius
sumber
@ kjetilbhalvorsen apakah Anda yakin ini duplikat karena OP tampaknya menginginkan interval prediksi tetapi tampaknya bekerja pada skala OR daripada skala log yang mungkin menjadi akar masalah?
mdewey
2
Jika Anda ingin mengevaluasi seberapa bagus regresi logistik, biasanya menggunakan ukuran yang berbeda dari prediksi + SE. Salah satu ukuran evaluasi populer adalah ROC-Curve dengan masing-masing AUC
adibender
1
Mungkinkah ini membantu? stackoverflow.com/questions/47414842/…
Xavier Bourret Sicotte

Jawaban:

24

Pertanyaan Anda mungkin berasal dari fakta bahwa Anda berurusan dengan Odds Ratios and Probabilities yang membingungkan pada awalnya. Karena model logistik adalah transformasi non linear dari menghitung interval kepercayaan tidak sesederhana.βTx

Latar Belakang

Ingat itu untuk model regresi logistik

  • Probabilitas dari : p = e α + β 1 x 1 + β 2 x 2(Y=1)p=eα+β1x1+β2x21+eα+β1x1+β2x2

  • Odds dari : ( p(Y=1)(p1p)=eα+β1x1+β2x2

  • Log Odds dari : log ( p(Y=1)log(p1p)=α+β1x1+β2x2

Pertimbangkan kasus di mana Anda memiliki kenaikan satu unit dalam variabel , yaitu x 1 + 1 , maka peluang baru adalahx1x1+1

Odds(Y=1)=eα+β1(x1+1)+β2x2=eα+β1x1+β1+β2x2
  • Odds Ratio (OR) oleh karena itu

Odds(x1+1)Odds(x1)=eα+β1(x1+1)+β2x2eα+β1x1+β2x2=eβ1
  • Log Odds Ratio = β1

  • Risiko relatif atau (rasio probabilitas) = eα+β1x1+β1+β2x21+eα+β1x1+β1+β2x2eα+β1x1+β2x21+eα+β1x1+β2x2

Menafsirkan koefisien

βj

  • xjβj
  • xjeβj
  • xjkk+ΔeβjΔ
  • xj

βj

1.96SE

βj

βj±zSE(βj)

eβj±zSE(βj)

yang merupakan interval kepercayaan pada rasio odds. Perhatikan bahwa interval ini hanya untuk parameter tunggal.

Jika saya ingin memahami kesalahan standar untuk kedua variabel, bagaimana saya mempertimbangkannya?

Jika Anda memasukkan beberapa parameter, Anda dapat menggunakan prosedur Bonferroni, jika tidak, untuk semua parameter Anda dapat menggunakan interval kepercayaan untuk estimasi probabilitas

Prosedur Bonferroni untuk beberapa parameter

g1α

βg±z(1α2g)SE(βg)

Interval kepercayaan untuk estimasi probabilitas

pPr(pLppU)=.95

Satu pendekatan yang disebut transformasi titik akhir melakukan hal berikut:

  • xTβ
  • F(xTβ)

Pr(xTβ)=F(xTβ)xTβ

[Pr(xTβ)LPr(xTβ)Pr(xTβ)U]=[F(xTβ)LF(xTβ)F(xTβ)U]

βTx±zSE(βTx)

[exTβzSE(xTβ)1+exTβzSE(xTβ),exTβ+zSE(xTβ)1+exTβ+zSE(xTβ),]

xTβ

Var(xTβ)=xTΣx

The advantage of this method is that the bounds cannot be outside the range (0,1)

There are several other approaches as well, using the delta method, bootstrapping etc.. which each have their own assumptions, advantages and limits.


Sources and info

My favorite book on this topic is "Applied Linear Statistical Models" by Kutner, Neter, Li, Chapter 14

Otherwise here are a few online sources:

Xavier Bourret Sicotte
sumber
Much of this is about CI for the coefficients which is a fine thing for the OP to know about but are we sure that is what he needs? You later section seems to me more relevant but perhaps the distinctions may be missed if read too quickly?
mdewey
2
Yes you are probably right - but understanding odds, log odds and probabilities for log regression is something I struggled with in the past - I hope this post summarises the topic well enough to such that it might help someone in the future. Perhaps I could answer the question more explicitly by providing a CI but we would need the covariance matrix
Xavier Bourret Sicotte
5

To get the 95% confidence interval of the prediction you can calculate on the logit scale and then convert those back to the probability scale 0-1. Here is an example using the titanic dataset.

library(titanic)
data("titanic_train")

titanic_train$Pclass = factor(titanic_train$Pclass, levels = c(1,2,3), labels = c('First','Second','Third'))

fit = glm(Survived ~ Sex + Pclass, data=titanic_train, family = binomial())

inverse_logit = function(x){
  exp(x)/(1+exp(x))
}

predicted = predict(fit, data.frame(Sex='male', Pclass='First'), type='link', se.fit=TRUE)

se_high = inverse_logit(predicted$fit + (predicted$se.fit*1.96))
se_low = inverse_logit(predicted$fit - (predicted$se.fit*1.96))
expected = inverse_logit(predicted$fit)

The mean and low/high 95% CI.

> expected
        1 
0.4146556 
> se_high
        1 
0.4960988 
> se_low
        1 
0.3376243 

And the output from just using type='response', which only gives the mean

predict(fit, data.frame(Sex='male', Pclass='First'), type='response')
        1 
0.4146556
Shawn
sumber
predict(fit, data.frame(Sex='male', Pclass='First'), type='response', se.fit=TRUE) will work.
Tony416