Calcul des intervalles de confiance pour une régression logistique

15

J'utilise une régression logistique binomiale pour identifier si l'exposition à has_xou has_yaffecte la probabilité qu'un utilisateur clique sur quelque chose. Mon modèle est le suivant:

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

Voici la sortie de mon modèle:

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

Comme chaque coefficient est significatif, en utilisant ce modèle, je peux dire quelle est la valeur de l'une de ces combinaisons en utilisant l'approche suivante:

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

Je ne comprends pas comment je peux faire rapport sur le Std. Erreur de prédiction.

  1. Dois-je simplement utiliser 1.96SE ? Ou dois-je convertir le SE utilisant une approche décrite ici ?

  2. Si je veux comprendre l'erreur standard pour les deux variables, comment puis-je considérer cela?

Contrairement à cette question , je suis intéressé à comprendre quelles sont les limites supérieure et inférieure de l'erreur en pourcentage. Par exemple, ma prédiction montre une valeur de 37% car True,Truepuis-je calculer que c'est pour un 95 % C I ? (0,3% choisi pour illustrer mon propos)+/0.395%CI

celenius
la source
Duplicate: stats.stackexchange.com/questions/5304/…
kjetil b halvorsen
@kjetilbhalvorsen êtes-vous sûr qu'il s'agit d'un doublon car l'OP semble vouloir un intervalle de prédiction mais semble fonctionner sur l'échelle OR plutôt que sur l'échelle logarithmique qui peut être à l'origine du problème?
mdewey
2
Si vous voulez évaluer la qualité d'une prédiction de régression logistique, on utilise généralement des mesures différentes de la prédiction + SE. Une mesure d'évaluation populaire est la courbe ROC avec l'AUC respective
adibender
1
Cela pourrait-il être utile? stackoverflow.com/questions/47414842/…
Xavier Bourret Sicotte

Réponses:

24

Votre question peut provenir du fait que vous traitez avec les rapports de cotes et les probabilités, ce qui prête à confusion au début. Puisque le modèle logistique est une transformation non linéaire de calculant les intervalles de confiance ne sont pas aussi simples.βTx

Contexte

Rappelons que pour le modèle de régression logistique

  • Probabilité de : p = e α + β 1 x 1 + β 2 x 2(Y=1)p=eα+β1x1+β2x21+eα+β1x1+β2x2

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

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

Considérons le cas où vous avez une augmentation d'une unité dans la variable , c'est-à-dire x 1 + 1 , alors les nouvelles cotes sontx1x1+1

Odds(Y=1)=eα+β1(x1+1)+β2x2=eα+β1x1+β1+β2x2
  • Les rapports de cotes (OR) sont donc

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

  • Risque relatif ou (rapport de probabilité) = eα+β1x1+β1+β2x21+eα+β1x1+β1+β2x2eα+β1x1+β2x21+eα+β1x1+β2x2

Interprétation des coefficients

Comment interpréteriez-vous la valeur du coefficient βj ? En supposant que tout le reste reste fixe:

  • Pour chaque augmentation unitaire de le log-odds ratio augmente de β jxjβj .
  • Pour chaque augmentation unitaire de le rapport de cotes augmente de e β jxjeβj .
  • Pour chaque augmentation de de k à k + Δ, le rapport de cotes augmente de e β j Δxjkk+ΔeβjΔ
  • Si le coefficient est négatif, une augmentation de entraîne une diminution du rapport de cotes.xj

Intervalles de confiance pour un seul paramètre βj

Dois-je simplement utiliser ? Ou dois-je convertir le SE en utilisant une approche décrite ici?1.96SE

Étant donné que le paramètre est estimé à l'aide de l'estimation du maximum de vraisemblance, la théorie MLE nous dit qu'il est asymptotiquement normal et que nous pouvons donc utiliser l' intervalle de confiance Wald à grand échantillon pour obtenir l'habituelβj

βj±zSE(βj)

Ce qui donne un intervalle de confiance sur le log-odds ratio. L'utilisation de la propriété d'invariance du MLE permet d'exponentier pour obtenir

eβj±zSE(βj)

qui est un intervalle de confiance sur le rapport de cotes. Notez que ces intervalles ne concernent qu'un seul paramètre.

Si je veux comprendre l'erreur standard pour les deux variables, comment puis-je considérer cela?

Si vous incluez plusieurs paramètres, vous pouvez utiliser la procédure de Bonferroni, sinon pour tous les paramètres, vous pouvez utiliser l'intervalle de confiance pour les estimations de probabilité

Procédure de Bonferroni pour plusieurs paramètres

g1α

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

Intervalles de confiance pour les estimations de probabilité

The logistic model outputs an estimation of the probability of observing a one and we aim to construct a frequentist interval around the true probability p such that Pr(pLppU)=.95

One approach called endpoint transformation does the following:

  • Compute the upper and lower bounds of the confidence interval for the linear combination xTβ (using the Wald CI)
  • Apply a monotonic transformation to the endpoints F(xTβ) to obtain the probabilities.

Since Pr(xTβ)=F(xTβ) is a monotonic transformation of xTβ

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

Concretely this means computing βTx±zSE(βTx) and then applying the logit transform to the result to get the lower and upper bounds:

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

The estimated approximate variance of xTβ can be calculated using the covariance matrix of the regression coefficients using

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
la source
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
la source
predict(fit, data.frame(Sex='male', Pclass='First'), type='response', se.fit=TRUE) will work.
Tony416