Régression logistique: test anova-chi-carré vs signification des coefficients (anova () vs summary () en R)

35

J'ai un modèle logistique GLM avec 8 variables. J'ai effectué un test du chi-carré dans R anova(glm.model,test='Chisq')et 2 des variables se révèlent être prédictives lorsqu'elles sont ordonnées en haut du test et pas tellement lorsqu'elles sont ordonnées en bas. La summary(glm.model)donne à penser que leurs coefficients ne sont pas significatifs (p-valeur élevée). Dans ce cas, il semble que les variables ne sont pas significatives.

Je voulais demander quel est le meilleur test de la signification des variables - la signification du coefficient dans le résumé du modèle ou le test du Khi-deux de anova(). Aussi - quand l'un ou l'autre est-il meilleur que l'autre?

Je suppose que c’est une question générale, mais nous apprécierons tous les indicateurs.

StreetHawk
la source
4
Ceci est analogue à la distinction entre les sommes des carrés de type I et de type III pour tester les coefficients dans les modèles linéaires. Cela peut vous aider de lire ma réponse ici: comment interpréter une ANOVA séquentielle de type I et une MANOVA .
gung - Réintégrer Monica

Réponses:

61

En plus de la réponse de @ gung, je vais essayer de donner un exemple de ce que la anovafonction teste réellement. J'espère que cela vous permettra de décider quels tests sont appropriés pour les hypothèses que vous souhaitez tester.

yX1X2X3my.mod <- glm(y~x1+x2+x3, family="binomial")anova(my.mod, test="Chisq")

  1. glm(y~1, family="binomial") contre. glm(y~x1, family="binomial")
  2. glm(y~x1, family="binomial") contre. glm(y~x1+x2, family="binomial")
  3. glm(y~x1+x2, family="binomial") contre. glm(y~x1+x2+x3, family="binomial")

Ainsi, il compare séquentiellement le plus petit modèle avec le prochain modèle plus complexe en ajoutant une variable à chaque étape. Chacune de ces comparaisons est effectuée via un test du rapport de vraisemblance (test LR; voir l'exemple ci-dessous). À ma connaissance, ces hypothèses sont rarement intéressantes, mais vous devez en décider.

Voici un exemple dans R:

mydata      <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank)

my.mod <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
summary(my.mod)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
   ---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 

# The sequential analysis
anova(my.mod, test="Chisq")

Terms added sequentially (first to last)    

     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                   399     499.98              
gre   1  13.9204       398     486.06 0.0001907 ***
gpa   1   5.7122       397     480.34 0.0168478 *  
rank  3  21.8265       394     458.52 7.088e-05 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

# We can make the comparisons by hand (adding a variable in each step)

  # model only the intercept
mod1 <- glm(admit ~ 1,                data = mydata, family = "binomial") 
  # model with intercept + gre
mod2 <- glm(admit ~ gre,              data = mydata, family = "binomial") 
  # model with intercept + gre + gpa
mod3 <- glm(admit ~ gre + gpa,        data = mydata, family = "binomial") 
  # model containing all variables (full model)
mod4 <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial") 

anova(mod1, mod2, test="LRT")

Model 1: admit ~ 1
Model 2: admit ~ gre
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       399     499.98                          
2       398     486.06  1    13.92 0.0001907 ***

anova(mod2, mod3, test="LRT")

Model 1: admit ~ gre
Model 2: admit ~ gre + gpa
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
1       398     486.06                       
2       397     480.34  1   5.7122  0.01685 *

anova(mod3, mod4, test="LRT")

Model 1: admit ~ gre + gpa
Model 2: admit ~ gre + gpa + rank
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       397     480.34                          
2       394     458.52  3   21.826 7.088e-05 ***

psummary(my.mod)

  • Pour coefficient de x1: glm(y~x2+x3, family="binomial")vs. glm(y~x1+x2+x3, family="binomial")
  • Pour coefficient de x2: glm(y~x1+x3, family="binomial")vs.glm(y~x1+x2+x3, family="binomial")
  • Pour coefficient de x3: glm(y~x1+x2, family="binomial")vs.glm(y~x1+x2+x3, family="binomial")

Donc, chaque coefficient par rapport au modèle complet contenant tous les coefficients. Les tests de Wald sont une approximation du test du rapport de vraisemblance. Nous pourrions également faire les tests du rapport de vraisemblance (test LR). Voici comment:

mod1.2 <- glm(admit ~ gre + gpa,  data = mydata, family = "binomial")
mod2.2 <- glm(admit ~ gre + rank, data = mydata, family = "binomial")
mod3.2 <- glm(admit ~ gpa + rank, data = mydata, family = "binomial")

anova(mod1.2, my.mod, test="LRT") # joint LR test for rank

Model 1: admit ~ gre + gpa
Model 2: admit ~ gre + gpa + rank
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1       397     480.34                          
2       394     458.52  3   21.826 7.088e-05 ***

anova(mod2.2, my.mod, test="LRT") # LR test for gpa

Model 1: admit ~ gre + rank
Model 2: admit ~ gre + gpa + rank
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
1       395     464.53                       
2       394     458.52  1   6.0143  0.01419 *

anova(mod3.2, my.mod, test="LRT") # LR test for gre

Model 1: admit ~ gpa + rank
Model 2: admit ~ gre + gpa + rank
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
1       395     462.88                       
2       394     458.52  1   4.3578  0.03684 *

le psummary(my.mod)

rankanova(my.mod, test="Chisq")rankanova(mod1.2, my.mod, test="Chisq")p7.088dix-5 . C'est à chaque fois la comparaison entre le modèle sans ranket le modèle qui le contient.

COOLSerdash
la source
1
+1, c'est une bonne explication complète. Un petit point: je pense que lorsque test="Chisq"vous n'exécutez pas de test du rapport de probabilité, vous devez définir test="LRT"pour cela, voir ? Anova.glm .
gung - Réintégrer Monica
6
@gung Merci pour le compliment. test="LRT"et test="Chisq"sont synonymes (il est dit sur la page que vous avez lié).
COOLSerdash
2
Pas de problème, mais je pense que c'est un bon point. test="LRT"est mieux car il est tout de suite clair que c’est un test du rapport de vraisemblance. Je l'ai changé. Merci.
COOLSerdash
4
+1 Je suis impressionné par les progrès rapides que vous avez réalisés ici en seulement un mois et par votre capacité à fournir une explication claire et bien travaillée. Merci pour vos efforts!
whuber
1
Très bonne réponse. Puis-je demander comment les valeurs-p ( 7.088e-05, 0.01419, 00.03684) doivent être interprétées?
TheSimpliFire le