¿Por qué hay una diferencia entre calcular manualmente un intervalo de confianza del 95% de regresión logística y usar la función confint () en R?

34

Queridos todos, he notado algo extraño que no puedo explicar, ¿verdad? En resumen: el enfoque manual para calcular un intervalo de confianza en un modelo de regresión logística y la función R confint()dan resultados diferentes.

He estado pasando por la regresión logística aplicada de Hosmer & Lemeshow (2a edición). En el tercer capítulo hay un ejemplo de cálculo de la razón de posibilidades y el intervalo de confianza del 95%. Usando R, puedo reproducir fácilmente el modelo:

Call:
glm(formula = dataset$CHD ~ as.factor(dataset$dich.age), family = "binomial")

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.734  -0.847  -0.847   0.709   1.549  

Coefficients:
                             Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -0.8408     0.2551  -3.296  0.00098 ***
as.factor(dataset$dich.age)1   2.0935     0.5285   3.961 7.46e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 136.66  on 99  degrees of freedom
Residual deviance: 117.96  on 98  degrees of freedom
AIC: 121.96

Number of Fisher Scoring iterations: 4

Sin embargo, cuando calculo los intervalos de confianza de los parámetros, obtengo un intervalo diferente al indicado en el texto:

> exp(confint(model))
Waiting for profiling to be done...
                                 2.5 %     97.5 %
(Intercept)                  0.2566283  0.7013384
as.factor(dataset$dich.age)1 3.0293727 24.7013080

Hosmer y Lemeshow sugieren la siguiente fórmula:

{mi}^{[{\hat{β}}_{1} \pm z_{1 - α / / 2} \times \hat{SE} ({\hat{β}}_{1})]}

$e^{[\hat\beta_1\pm z_{1-\alpha/2}\times\hat{\text{SE}}(\hat\beta_1)]}$

y calculan el intervalo de confianza para as.factor(dataset$dich.age)1ser (2.9, 22.9).

Esto parece sencillo de hacer en R:

# upper CI for beta
exp(summary(model)$coefficients[2,1]+1.96*summary(model)$coefficients[2,2])
# lower CI for beta
exp(summary(model)$coefficients[2,1]-1.96*summary(model)$coefficients[2,2])

da la misma respuesta que el libro.

Sin embargo, ¿alguna idea de por qué confint()parece dar resultados diferentes? He visto muchos ejemplos de personas que usan confint().

r regression logistic confidence-interval profile-likelihood correlation mcmc error mixture measurement data-augmentation r logistic goodness-of-fit r time-series exponential descriptive-statistics average expected-value data-visualization anova teaching hypothesis-testing multivariate-analysis r r mixed-model clustering categorical-data unsupervised-learning r logistic anova binomial estimation variance expected-value r r anova mixed-model multiple-comparisons repeated-measures project-management r poisson-distribution control-chart project-management regression residuals r distributions data-visualization r unbiased-estimator kurtosis expected-value regression spss meta-analysis r censoring regression classification data-mining mixture Andrés
fuente

1

¿Le importaría agregar la referencia bibliográfica exacta para Hosmer & Lemeshow? He estado buscando la sugerencia en sus publicaciones y libros durante bastante tiempo, pero aún no la he encontrado.

DavidR

36

Después de obtener los datos del sitio web adjunto , así es como lo haría:

chdage <- read.table("chdage.dat", header=F, col.names=c("id","age","chd"))
chdage$aged <- ifelse(chdage$age>=55, 1, 0)
mod.lr <- glm(chd ~ aged, data=chdage, family=binomial)
summary(mod.lr)

Los IC del 95% basados en la probabilidad de perfil se obtienen con

require(MASS)
exp(confint(mod.lr))

Esto suele ser el valor predeterminado si el MASSpaquete se carga automáticamente. En este caso, me sale

                2.5 %     97.5 %
(Intercept) 0.2566283  0.7013384
aged        3.0293727 24.7013080

Ahora, si quisiera comparar con el 95% de los CI de Wald (basados en la normalidad asintótica) como el que calculó a mano, usaría en su confint.default()lugar; esto produce

                2.5 %     97.5 %
(Intercept) 0.2616579  0.7111663
aged        2.8795652 22.8614705

Los IC de Wald son buenos en la mayoría de las situaciones, aunque el perfil basado en la probabilidad puede ser útil con estrategias de muestreo complejas. Si desea comprender la idea de cómo funcionan, aquí hay una breve descripción de los principios principales: Intervalos de confianza por el método de probabilidad de perfil, con aplicaciones en epidemiología veterinaria . También puede echar un vistazo al libro MASS de Venables y Ripley, §8.4, pp. 220-221.

chl
fuente

25

Seguimiento: los intervalos de confianza del perfil son más confiables (elegir el punto de corte apropiado para la probabilidad implica una suposición asintótica (muestra grande), pero esta es una suposición mucho más débil que la suposición de la superficie de probabilidad cuadrática subyacente a los intervalos de confianza de Wald). Hasta donde yo sé, no hay argumento para las estadísticas de Wald sobre los intervalos de confianza del perfil, excepto que las estadísticas de Wald son mucho más rápidas de calcular y pueden ser "lo suficientemente buenas" en muchas circunstancias (pero a veces fuera de lugar: busque Hauck- Efecto Donner).

Ben Bolker
fuente

2

Gracias por esto y por sugerir que busque el efecto Hauck-Donner. ¡El efecto no recibe mucho tratamiento en los libros de texto, pero parece bastante importante!

Andrew

18

Creo que si busca confint () en el archivo de ayuda, encontrará que el intervalo de confianza que se está construyendo es un intervalo de "perfil" en lugar de un intervalo de confianza de Wald (su fórmula de HL).

B_Miner
fuente

55

Ahh Eso responde la pregunta. Sin embargo, lleva al siguiente, ¿cuál es el preferido?

Andrew

¿Por qué hay una diferencia entre calcular manualmente un intervalo de confianza del 95% de regresión logística y usar la función confint () en R?

Respuestas: