Calcular intervalos de confianza para una regresión logística

15

Estoy usando una regresión logística binomial para identificar si la exposición has_xo has_yimpacto tiene la probabilidad de que un usuario haga clic en algo. Mi modelo es el siguiente:

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

Este es el resultado de mi modelo:

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

Como cada coeficiente es significativo, usando este modelo puedo decir cuál es el valor de cualquiera de estas combinaciones usando el siguiente enfoque:

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

No entiendo cómo puedo informar sobre el estándar. Error de la predicción.

  1. ¿Solo necesito usar 1.96SE ? ¿O necesito convertir el SE usando un enfoque descrito aquí ?

  2. Si quiero entender el error estándar para ambas variables, ¿cómo lo consideraría?

A diferencia de esta pregunta , estoy interesado en comprender cuáles son los límites superior e inferior del error en un porcentaje. Por ejemplo, de mi predicción muestra un valor de 37% para True,True¿puedo calcular que esto es +/0.3 para un 95%CI ? (0.3% elegido para ilustrar mi punto)

celenius
fuente
Duplicado: stats.stackexchange.com/questions/5304/…
kjetil b halvorsen
@kjetilbhalvorsen ¿está seguro de que es un duplicado ya que el OP parece querer un intervalo de predicción pero parece estar funcionando en la escala OR en lugar de la escala logarítmica que puede ser la raíz del problema?
mdewey
2
Si desea evaluar qué tan buena predice una regresión logística, generalmente se usan diferentes medidas que predicción + SE. Una medida de evaluación popular es la curva ROC con las AUC respectivas
adibender el
1
¿Podría esto ser de alguna ayuda? stackoverflow.com/questions/47414842/…
Xavier Bourret Sicotte

Respuestas:

24

Su pregunta puede provenir del hecho de que se trata de Odds Ratios y Probabilidades, lo cual es confuso al principio. Dado que el modelo logístico es una transformación no lineal de la computación , los intervalos de confianza no son tan sencillos.βTx

Antecedentes

Recordemos que para el modelo de regresión logística

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

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

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

Considere el caso en el que tiene un aumento de una unidad en la variable , es decir, x 1 + 1 , entonces las nuevas probabilidades sonx1x1+1

Odds(Y=1)=eα+β1(x1+1)+β2x2=eα+β1x1+β1+β2x2
  • Por lo tanto, la Odds Ratio (OR) es

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

  • Riesgo relativo o (razón de probabilidad) = eα+β1x1+β1+β2x21+eα+β1x1+β1+β2x2eα+β1x1+β2x21+eα+β1x1+β2x2

Coeficientes de interpretación

¿Cómo interpretaría el valor del coeficiente ? Suponiendo que todo lo demás permanece fijo:βj

  • Por cada unidad de aumento en la relación log-odds aumenta en β j .xjβj
  • Por cada unidad de aumento en la razón de posibilidades aumenta en e β j .xjeβj
  • For every increase of xj from k to k+Δ the odds ratio increases by eβjΔ
  • If the coefficient is negative, then an increase in xj leads to a decrease in the odds ratio.

Confidence intervals for a single parameter βj

Do I just need to use 1.96SE? Or do I need to convert the SE using an approach described here?

Since the parameter βj is estimated using Maxiumum Likelihood Estimation, MLE theory tells us that it is asymptotically normal and hence we can use the large sample Wald confidence interval to get the usual

βj±zSE(βj)

Which gives a confidence interval on the log-odds ratio. Using the invariance property of the MLE allows us to exponentiate to get

eβj±zSE(βj)

which is a confidence interval on the odds ratio. Note that these intervals are for a single parameter only.

If I want to understand the standard-error for both variables how would I consider that?

If you include several parameters you can use the Bonferroni procedure, otherwise for all parameters you can use the confidence interval for probability estimates

Bonferroni procedure for several parameters

If g parameters are to be estimated with family confidence coefficient of approximately 1α, the joint Bonferroni confidence limits are

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

Confidence intervals for probability estimates

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