Teste o modelo de regressão logística usando desvio residual e graus de liberdade

Eu estava lendo esta página em Princeton.edu . Eles estão executando uma regressão logística (com R). Em algum momento, eles calculam a probabilidade de obter um desvio residual maior do que aquele que obtiveram em uma com graus de liberdade iguais aos graus de liberdade do modelo. Copiando e colando do site deles ... $\chi^2$

> glm( cbind(using,notUsing) ~ age + hiEduc + noMore, family=binomial)

Call:  glm(formula = cbind(using, notUsing) ~ age + hiEduc + noMore,      
     family = binomial) 

Coefficients:
(Intercept)     age25-29     age30-39     age40-49       hiEduc       noMore  
    -1.9662       0.3894       0.9086       1.1892       0.3250       0.8330  

Degrees of Freedom: 15 Total (i.e. Null);  10 Residual
Null Deviance:      165.8 
Residual Deviance: 29.92        AIC: 113.4

O desvio residual de 29,92 em 10 df é altamente significativo:

> 1-pchisq(29.92,10)
[1] 0.0008828339

então precisamos de um modelo melhor

Por que faz sentido calcular 1-pchisq(29.92,10)e por que uma baixa probabilidade indica que algo está errado com seu modelo?

regression hypothesis-testing logistic generalized-linear-model chi-squared Remi.b
fonte

Respostas:

Eles estão usando um teste de desvio mostrado abaixo:

D (y) = - 2 ℓ (\hat{β}; y) + 2 ℓ ({\hat{θ}}^{(s)}; y)

$D(y) = -2\ell(\hat\beta;y) + 2\ell(\hat\theta^{(s)};y)$

Aqui, o representa o modelo de interesse ajustado e representa o modelo saturado. A probabilidade de log para o modelo saturado é (mais frequentemente do que não) , portanto, você fica com o desvio residual do modelo que eles ajustaram ( ). Esse teste de desvio é aproximadamente qui-quadrado com graus de liberdade ( sendo as observações sendo o número de variáveis ajustadas). Tem e de modo que o teste será de aproximadamente $\hat β$ $\hatθ(s)$ $0$ $29.92$ $n-p$ $n$ $p$ $n=16$ $p=6$ $\chi^2_{10}$ . O ponto nulo do teste é que seu modelo ajustado se ajusta bem aos dados e não há desajuste - você não perdeu nenhuma fonte de variação. No teste acima, você rejeita o valor nulo e, como resultado, perdeu algo no modelo que ajustou. A razão para usar esse teste é que o modelo saturado ajustará os dados perfeitamente; portanto, se você não rejeitar o nulo entre o modelo ajustado e o modelo saturado, isso indica que você não perdeu grandes fontes de dados variação no seu modelo.

francium87d
fonte

Sua pergunta, como indicado, foi respondida por @ francium87d. A comparação do desvio residual com a distribuição qui-quadrado apropriada constitui um teste do modelo ajustado em relação ao modelo saturado e mostra, nesse caso, uma falta significativa de ajuste.

Ainda assim, pode ser útil examinar mais detalhadamente os dados e o modelo para entender melhor o que significa que o modelo não tem ajuste:

d = read.table(text=" age education wantsMore notUsing using 
   <25       low       yes       53     6
   <25       low        no       10     4
   <25      high       yes      212    52
   <25      high        no       50    10
 25-29       low       yes       60    14
 25-29       low        no       19    10
 25-29      high       yes      155    54
 25-29      high        no       65    27
 30-39       low       yes      112    33
 30-39       low        no       77    80
 30-39      high       yes      118    46
 30-39      high        no       68    78
 40-49       low       yes       35     6
 40-49       low        no       46    48
 40-49      high       yes        8     8
 40-49      high        no       12    31", header=TRUE, stringsAsFactors=FALSE)
d = d[order(d[,3],d[,2]), c(3,2,1,5,4)]

library(binom)
d$proportion = with(d, using/(using+notUsing))
d$sum        = with(d, using+notUsing)
bCI          = binom.confint(x=d$using, n=d$sum, methods="exact")

m     = glm(cbind(using,notUsing)~age+education+wantsMore, d, family=binomial)
preds = predict(m, new.data=d[,1:3], type="response")

windows()
  par(mar=c(5, 8, 4, 2))
  bp = barplot(d$proportion, horiz=T, xlim=c(0,1), xlab="proportion",
               main="Birth control usage")
  box()
  axis(side=2, at=bp, labels=paste(d[,1], d[,2], d[,3]), las=1)
  arrows(y0=bp, x0=bCI[,5], x1=bCI[,6], code=3, angle=90, length=.05)
  points(x=preds, y=bp, pch=15, col="red")

A figura representa a proporção observada de mulheres em cada conjunto de categorias que usam controle de natalidade, juntamente com o intervalo exato de confiança de 95%. As proporções previstas do modelo são sobrepostas em vermelho. Podemos ver que duas proporções previstas estão fora dos ICs de 95% e outras cinco estão no ou muito perto dos limites dos respectivos ICs. São sete dos dezesseis ( ) que estão fora do objetivo. Portanto, as previsões do modelo não correspondem muito bem aos dados observados. $44\%$

Como o modelo poderia se encaixar melhor? Talvez haja interações entre as variáveis relevantes. Vamos adicionar todas as interações bidirecionais e avaliar o ajuste:

m2 = glm(cbind(using,notUsing)~(age+education+wantsMore)^2, d, family=binomial)
summary(m2)
# ...
#     Null deviance: 165.7724  on 15  degrees of freedom
# Residual deviance:   2.4415  on  3  degrees of freedom
# AIC: 99.949
# 
# Number of Fisher Scoring iterations: 4
1-pchisq(2.4415, df=3)  # [1] 0.4859562
drop1(m2, test="LRT")
# Single term deletions
# 
# Model:
# cbind(using, notUsing) ~ (age + education + wantsMore)^2
#                     Df Deviance     AIC     LRT Pr(>Chi)  
# <none>                   2.4415  99.949                   
# age:education        3  10.8240 102.332  8.3826  0.03873 *
# age:wantsMore        3  13.7639 105.272 11.3224  0.01010 *
# education:wantsMore  1   5.7983 101.306  3.3568  0.06693 .

O valor de p para o teste de falta de ajuste para este modelo é agora . Mas realmente precisamos de todos esses termos extras de interação? O comando mostra os resultados dos testes de modelo aninhados sem eles. A interação entre e não é muito significativa, mas eu ficaria bem com ela no modelo de qualquer maneira. Então, vamos ver como as previsões desse modelo se comparam aos dados: $0.486$ drop1()educationwantsMore

Elas não são perfeitas, mas não devemos assumir que as proporções observadas são um reflexo perfeito do verdadeiro processo de geração de dados. Parece-me que eles estão pulando em torno da quantidade apropriada (mais corretamente que os dados estão pulando em torno das previsões, suponho).

- Reinstate Monica
fonte

Não acredito que a estatística de desvio residual tenha uma . Eu acho que é uma distribuição degenerada porque a teoria assintótica não se aplica quando os graus de liberdade aumentam na mesma velocidade que o tamanho da amostra. De qualquer forma, duvido que o teste tenha poder suficiente e encorajo testes direcionados, como testes de linearidade usando splines de regressão e testes de interação. $\chi^2$

Frank Harrell
fonte

Eu acho que porque neste caso todos os preditores são categóricos, o não. os graus de liberdade do modelo saturado não aumentariam com o tamanho da amostra; portanto, a abordagem assintótica faz sentido. O tamanho da amostra ainda é bastante pequeno.

Scortchi - Restabelece Monica

Não tenho certeza disso. O df dos parâmetros do modelo é fixo, mas o df do " " residual é menos isso.

χ^{2}

$\chi^2$

n

$n$

Frank Harrell

Nesse caso, os dados consistem em 1607 indivíduos em uma tabela de contingência e o teste está comparando um modelo de 6 parâmetros com o modelo saturado de 16 parâmetros (em vez de um modelo de 1607 parâmetros).

Scortchi - Restabelece Monica

Então não deve ser rotulado como residual .

χ^{2}

$\chi^2$

31716 Frank

Concordo que essa terminologia é infeliz: glmfornece um "desvio residual" diferente quando os dados são agrupados quando não são - e um "desvio nulo" diferente para esse assunto.

Scortchi - Restabelece Monica