Obliczanie przedziałów ufności dla regresji logistycznej

15

Korzystam z dwumianowej regresji logistycznej, aby określić, czy narażenie has_xlub ma has_ywpływ na prawdopodobieństwo kliknięcia przez użytkownika. Mój model jest następujący:

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

To wynik z mojego modelu:

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

Ponieważ każdy współczynnik jest znaczący, za pomocą tego modelu jestem w stanie powiedzieć, jaka jest wartość dowolnej z tych kombinacji, stosując następujące podejście:

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

Nie rozumiem, jak mogę raportować o Std. Błąd prognozy.

  1. Czy muszę po prostu użyć ? Czy też muszę konwertować za pomocą opisanego tutaj podejścia ?S E1.96SESE

  2. Jeśli chcę zrozumieć błąd standardowy dla obu zmiennych, jak bym to rozważyć?

W przeciwieństwie do tego pytania , jestem zainteresowany zrozumieniem, jakie są górne i dolne granice błędu. Na przykład moje przewidywania pokazują wartość 37% dla True,Trueczy mogę obliczyć, że jest to dla ? (Wybrano 0,3%, aby zilustrować mój punkt widzenia)95 % C I+/0.395%CI

celenius
źródło
Duplikat: stats.stackexchange.com/questions/5304/…
kjetil b halvorsen
@kjetilbhalvorsen, czy jesteś pewien, że jest to duplikat, ponieważ OP wydaje się chcieć interwału przewidywania, ale wydaje się, że działa w skali OR, a nie w skali dziennika, która może być źródłem problemu?
mdewey
2
Jeśli chcesz ocenić, jak dobra jest regresja logistyczna, zwykle stosuje się inne miary niż przewidywanie + SE. Jednym z popularnych mierników oceny jest krzywa ROC z odpowiednim AUC
adibender
1
Czy to może pomóc? stackoverflow.com/questions/47414842/…
Xavier Bourret Sicotte

Odpowiedzi:

24

Twoje pytanie może wynikać z faktu, że masz do czynienia z ilorazem szans i prawdopodobieństwami, co na początku jest mylące. Ponieważ model logistyczny jest nieliniową transformacją obliczeń przedziały ufności nie są tak proste.βTx

tło

Przypomnijmy, że dla modelu regresji logistycznej

  • Prawdopodobieństwo z :p = e α + β 1 x 1 + β 2 x 2(Y=1)p=eα+β1x1+β2x21+eα+β1x1+β2x2

  • Szanse na :( str(Y=1)(p1p)=eα+β1x1+β2x2

  • Log Szanse na :log ( str(Y=1)log(p1p)=α+β1x1+β2x2

Rozważ przypadek, w którym masz wzrost o jedną jednostkę w zmiennej , tj. , wtedy nowe szanse sąx1x1+1

Odds(Y=1)=eα+β1(x1+1)+β2x2=eα+β1x1+β1+β2x2
  • Iloraz szans (OR) jest zatem

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

  • Ryzyko względne lub (współczynnik prawdopodobieństwa) =eα+β1x1+β1+β2x21+eα+β1x1+β1+β2x2eα+β1x1+β2x21+eα+β1x1+β2x2

Współczynniki interpretacyjne

Jak interpretowałbyś wartość współczynnika ? Zakładając, że wszystko inne pozostaje nieruchome:βj

  • Z każdym wzrostem jednostki współczynnik logarytmiczny zwiększa się o .xjβj
  • Z każdym wzrostem jednostki iloraz szans wzrasta o .xjeβj
  • Z każdym wzrostem od do współczynnik szans wzrasta oxjkk+ΔeβjΔ
  • Jeśli współczynnik jest ujemny, wówczas wzrost prowadzi do zmniejszenia ilorazu szans.xj

Przedziały ufności dla jednego parametruβj

Czy muszę tylko użyć ? Czy też muszę konwertować SE za pomocą opisanego tutaj podejścia?1.96SE

Ponieważ parametr jest szacowany za pomocą oszacowania wiarygodności Maxiumum, teoria MLE mówi nam, że jest asymptotycznie normalna, a zatem możemy użyć dużego przedziału ufności Walda, aby uzyskać zwykłeβj

βj±zSE(βj)

Co daje przedział ufności w stosunku logarytmicznym. Użycie właściwości niezmienniczości MLE pozwala nam na wykładnik, aby uzyskać

eβj±zSE(βj)

który jest przedziałem ufności dla ilorazu szans. Pamiętaj, że te przedziały dotyczą tylko jednego parametru.

Jeśli chcę zrozumieć błąd standardowy dla obu zmiennych, jak bym to rozważyć?

Jeśli podasz kilka parametrów, możesz użyć procedury Bonferroniego, w przeciwnym razie dla wszystkich parametrów możesz użyć przedziału ufności do oszacowania prawdopodobieństwa

Procedura Bonferroniego dla kilku parametrów

Jeżeli parametry mają być oszacowane przy współczynniku ufności rodziny wynoszącym około , wspólne granice ufności Bonferroniego wynosząg1α

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

Przedziały ufności dla oszacowań prawdopodobieństwa

Model logistyczny generuje oszacowanie prawdopodobieństwa zaobserwowania jednego, a naszym celem jest zbudowanie interwału częstości wokół rzeczywistego prawdopodobieństwa tak abypPr(pLppU)=.95

Jedno podejście zwane transformacją punktu końcowego wykonuje następujące czynności:

  • Oblicz górne i dolne granice przedziału ufności dla kombinacji liniowej (używając Wald CI)xTβ
  • Zastosuj transformację monotoniczną do punktów końcowych aby uzyskać prawdopodobieństwa.F(xTβ)

Ponieważ jest monotoniczną transformacjąPr(xTβ)=F(xTβ)xTβ

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

Konkretnie oznacza to obliczenie a następnie zastosowanie transformacji logit do wyniku w celu uzyskania dolnej i górnej granicy:βTx±zSE(βTx)

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

Oszacowaną przybliżoną wariancję można obliczyć za pomocą macierzy kowariancji współczynników regresji za pomocąxTβ

Var(xTβ)=xTΣx

Zaletą tej metody jest to, że granice nie mogą znajdować się poza zakresem(0,1)

Istnieje również kilka innych podejść, przy użyciu metody delta, ładowania początkowego itp., Które mają swoje własne założenia, zalety i ograniczenia.


Źródła i informacje

Moją ulubioną książką na ten temat jest „Zastosowane liniowe modele statystyczne” Kutnera, Neter, Li, rozdział 14

W przeciwnym razie oto kilka źródeł online:

Xavier Bourret Sicotte
źródło
Wiele z tego dotyczy CI dla współczynników, co jest dobrą rzeczą dla OP, ale czy jesteśmy pewni, że tego właśnie potrzebuje? Późniejsze sekcje wydają mi się bardziej odpowiednie, ale być może rozróżnienia mogą zostać pominięte, jeśli zostaną przeczytane zbyt szybko?
mdewey,
2
Tak, prawdopodobnie masz rację - ale zrozumienie szans, logarytmów i prawdopodobieństwa regresji logów było czymś, z czym borykałem się w przeszłości - mam nadzieję, że ten post podsumuje ten temat wystarczająco dobrze, aby mógł komuś pomóc w przyszłości. Być może mógłbym bardziej precyzyjnie odpowiedzieć na to pytanie, podając CI, ale potrzebowalibyśmy macierzy kowariancji
Xavier Bourret Sicotte
5

Aby uzyskać 95% przedział ufności prognozy, możesz obliczyć na skali logit, a następnie przekonwertować je z powrotem na skalę prawdopodobieństwa 0-1. Oto przykład wykorzystujący tytaniczny zestaw danych.

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)

Średni i niski / wysoki 95% CI.

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

I wynik z samego użycia type='response', co daje tylko średnią

predict(fit, data.frame(Sex='male', Pclass='First'), type='response')
        1 
0.4146556
Shawn
źródło
predict(fit, data.frame(Sex='male', Pclass='First'), type='response', se.fit=TRUE)będzie działać.
Tony416