Por que existe uma diferença entre calcular manualmente um intervalo de confiança de 95% da regressão logística e usar a função confint () em R?

34

Caro pessoal, notei algo estranho que não sei explicar, não é? Em resumo: a abordagem manual para calcular um intervalo de confiança em um modelo de regressão logística e a função R confint()fornecem resultados diferentes.

Eu tenho passado pela regressão logística aplicada de Hosmer & Lemeshow (2ª edição). No terceiro capítulo, há um exemplo de cálculo da razão de chances e intervalo de confiança de 95%. Usando R, posso reproduzir facilmente o 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

No entanto, quando calculo os intervalos de confiança dos parâmetros, obtenho um intervalo diferente daquele indicado no 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 & Lemeshow sugerem a seguinte fórmula:

e^{[{\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)]}$

e calculam o intervalo de confiança para as.factor(dataset$dich.age)1ser (2,9, 22,9).

Isso parece simples de fazer no 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])

dá a mesma resposta que o livro.

No entanto, alguma opinião sobre o porquê confint()parece dar resultados diferentes? Eu já vi muitos exemplos de pessoas usando 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 Andrew
fonte

1

Você gostaria de adicionar a referência exata da literatura para Hosmer & Lemeshow? Eu tenho procurado a sugestão em suas publicações e livros há algum tempo, mas ainda não a encontrei.

9788209Data

36

Depois de buscar os dados no site que acompanha , eis como eu faria isso:

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)

Os ICs de 95% com base na probabilidade do perfil são obtidos com

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

Geralmente, esse é o padrão se o MASSpacote for carregado automaticamente. Nesse caso, eu recebo

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

Agora, se eu quisesse comparar com ICs Wald 95% (com base na normalidade assintótica) como o que você calculou manualmente, eu usaria confint.default(); isso produz

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

Os ICs de Wald são bons na maioria das situações, embora o perfil baseado em probabilidade possa ser útil com estratégias de amostragem complexas. Se você quiser entender a idéia de como eles funcionam, eis uma breve visão geral dos principais princípios: Intervalos de confiança pelo método de probabilidade de perfil, com aplicações em epidemiologia veterinária . Você também pode dar uma olhada no livro MASS de Venables e Ripley, §8.4, pp. 220-221.

chl
fonte

25

Acompanhamento: os intervalos de confiança do perfil são mais confiáveis (a escolha do ponto de corte apropriado para a probabilidade envolve uma suposição assintótica (amostra grande), mas essa é uma suposição muito mais fraca que a suposição da superfície da probabilidade quadrática subjacente aos intervalos de confiança de Wald). Até onde eu sei, não há argumento para as estatísticas de Wald sobre os intervalos de confiança do perfil, exceto que as estatísticas de Wald são muito mais rápidas de calcular e podem ser "boas o suficiente" em muitas circunstâncias (mas às vezes bem distantes: consulte a Efeito Donner).

Ben Bolker
fonte

2

Obrigado por isso e por sugerir que eu procure o efeito Hauck-Donner. O efeito não recebe muito tratamento nos livros didáticos, mas parece muito importante!

Andrew

18

Acredito que se você procurar no arquivo de ajuda por confint (), descobrirá que o intervalo de confiança que está sendo construído é um intervalo de "perfil" em vez de um intervalo de confiança de Wald (sua fórmula da HL).

B_Miner
fonte

5

Ahh Isso responde à pergunta. No entanto, leva ao próximo - qual é o preferido?

Andrew

Por que existe uma diferença entre calcular manualmente um intervalo de confiança de 95% da regressão logística e usar a função confint () em R?

Respostas: