Obliczanie przedziałów ufności dla regresji logistycznej

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.

Czy muszę po prostu użyć ? Czy też muszę konwertować za pomocą opisanego tutaj podejścia ? $1.96*SE$ $SE$
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) $+/- 0.3%$ $95\% CI$

— celenius
źródło

Duplikat: stats.stackexchange.com/questions/5304/…

— kjetil b halvorsen

Możliwy duplikat Dlaczego istnieje różnica między ręcznym obliczeniem regresji logistycznej 95% przedziału ufności a użyciem funkcji confint () w R?

— 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

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

Czy to może pomóc? stackoverflow.com/questions/47414842/…

— Xavier Bourret Sicotte

Odpowiedzi:

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. $\beta^Tx$

tło

Przypomnijmy, że dla modelu regresji logistycznej

Prawdopodobieństwo z : $(Y = 1)$ $p = \frac{e^{\alpha + \beta_1x_1 + \beta_2 x_2}}{1 + e^{ \alpha + \beta_1x_1 + \beta_2 x_2}}$
Szanse na : $(Y = 1)$ $\left( \frac{p}{1-p}\right) = e^{\alpha + \beta_1x_1 + \beta_2 x_2}$
Log Szanse na : $(Y = 1)$ $\log \left( \frac{p}{1-p}\right) = \alpha + \beta_1x_1 + \beta_2 x_2$

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

Odds (Y = 1) = e^{α + β_{1} (x_{1} + 1) + β_{2} x_{2}} = e^{α + β_{1} x_{1} + β_{1} + β_{2} x_{2}}

$\text{Odds}(Y = 1) = e^{\alpha + \beta_1(x_1 + 1) + \beta_2x_2} = e^{\alpha + \beta_1 x_1 + \beta_1 + \beta_2x_2 }$

Iloraz szans (OR) jest zatem

\frac{Odds (x_{1} + 1)}{Odds (x_{1})} = \frac{e^{α + β_{1} (x_{1} + 1) + β_{2} x_{2}}}{e^{α + β_{1} x_{1} + β_{2} x_{2}}} = e^{β_{1}}

$\frac{\text{Odds}(x_1 + 1)}{\text{Odds}(x_1)} = \frac{e^{\alpha + \beta_1(x_1 + 1) + \beta_2x_2} }{e^{\alpha + \beta_1 x_1 + \beta_2x_2}} = e^{\beta_1}$

Log Iloraz szans = $\beta_1$
Ryzyko względne lub (współczynnik prawdopodobieństwa) = $\frac{ \frac{e^{\alpha + \beta_1x_1 + \beta_1 + \beta_2 x_2}}{1 + e^{ \alpha + \beta_1x_1 + \beta_1 + \beta_2 x_2}}}{ \frac{e^{\alpha + \beta_1x_1 + \beta_2 x_2}}{1 + e^{ \alpha + \beta_1x_1 + \beta_2 x_2}}}$

Współczynniki interpretacyjne

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

Z każdym wzrostem jednostki współczynnik logarytmiczny zwiększa się o . $x_j$ $\beta_j$
Z każdym wzrostem jednostki iloraz szans wzrasta o . $x_j$ $e^{\beta_j}$
Z każdym wzrostem od do współczynnik szans wzrasta o $x_j$ $k$ $k + \Delta$ $e^{\beta_j \Delta}$
Jeśli współczynnik jest ujemny, wówczas wzrost prowadzi do zmniejszenia ilorazu szans. $x_j$

Przedziały ufności dla jednego parametru $\beta_j$

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

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 $\beta_j$

β_{j} \pm z^{*} S E (β_{j})

$\beta_j \pm z^* SE(\beta_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} \pm z^{*} S E (β_{j})}

$e^{\beta_j \pm z^* SE(\beta_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ą $g$ $1 - \alpha$

β_{g} \pm z_{(1 - \frac{α}{2 g})} S E (β_{g})

$\beta_g \pm z_{(1 - \frac{\alpha}{2g})}SE(\beta_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 aby $p$ $Pr(p_{L} \leq p \leq p_{U}) = .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) $x^T\beta$
Zastosuj transformację monotoniczną do punktów końcowych aby uzyskać prawdopodobieństwa. $F(x^T\beta)$

Ponieważ jest monotoniczną transformacją $Pr(x^T\beta) = F(x^T\beta)$ $x^T\beta$

[P r (x^{T} β)_{L} \leq P r (x^{T} β) \leq P r (x^{T} β)_{U}] = [F (x^{T} β)_{L} \leq F (x^{T} β) \leq F (x^{T} β)_{U}]

$[Pr(x^T\beta)_L \leq Pr(x^T\beta) \leq Pr(x^T\beta)_U] = [F(x^T\beta)_L \leq F(x^T\beta) \leq F(x^T\beta)_U]$

Konkretnie oznacza to obliczenie a następnie zastosowanie transformacji logit do wyniku w celu uzyskania dolnej i górnej granicy: $\beta^Tx \pm z^* SE(\beta^Tx)$

[\frac{e^{x^{T} β - z^{*} S E (x^{T} β)}}{1 + e^{x^{T} β - z^{*} S E (x^{T} β)}}, \frac{e^{x^{T} β + z^{*} S E (x^{T} β)}}{1 + e^{x^{T} β + z^{*} S E (x^{T} β)}},]

$[\frac{e^{x^T\beta - z^* SE(x^T\beta)}}{1 + e^{x^T\beta - z^* SE(x^T\beta)}}, \frac{e^{x^T\beta + z^* SE(x^T\beta)}}{1 + e^{x^T\beta + z^* SE(x^T\beta)}},]$

Oszacowaną przybliżoną wariancję można obliczyć za pomocą macierzy kowariancji współczynników regresji za pomocą $x^T\beta$

V a r (x^{T} β) = x^{T} Σ x

$Var(x^T\beta) = x^T \Sigma 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,

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

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