Tests d’adéquation globale pour les modèles GLM
1 Exemple : GLM log-Poisson
1.1 Chargement des données et ajustement du modèle
On reprend ici (sans plus de commentaires) l’exemple déjà traité dans les illustrations précédentes.
library(outilscoursglm)load(url("https://irma.math.unistra.fr/~jberard/donnees_swautoins_ext_sim"))
donnees[["Zone"]] <- relevel(donnees[["Zone"]], ref = "4")
donnees[["Bonus"]] <- relevel(donnees[["Bonus"]], ref = "7")
donnees[["Make"]] <- relevel(donnees[["Make"]], ref = "7-9")
ajustement <- glm(Claims ~ Kilometres + Zone + Bonus + Make + offset(log(Insured)),
family = poisson(link = "log"),
data = donnees)1.2 Test d’adéquation basé sur la déviance, données non-groupées
On commence par appliquer l’approximation douteuse de la loi de la déviance par une loi du \(\chi^2\) à \(N-(d+1)\) degrés de liberté.
A cette fin, on extrait la valeur de \(N-(d+1)\) et de la déviance du modèle ajusté.
deg_lib_res <- df.residual(ajustement)
deg_lib_res## [1] 239421dev <- deviance(ajustement)
dev## [1] 220632.9On calcule maintenant la \(p-\)valeur du test obtenu avec l’approximation en question.
p_val_chi2 <- 1-pchisq(q = dev, df = deg_lib_res)
p_val_chi2## [1] 1Le résultat est en quelque sorte trop beau pour être vrai : interprété littéralement, l’hypothèse que le modèle est adéquat ne sera jamais rejetée, quel que soit le niveau de risque de première espèce que l’on se fixe. En fait, c’est un signe visible que l’approximation sur laquelle est basée ce test n’est pas valable ici.
Pour clarifier cela, et pouvoir disposer d’une procédure de test plus correcte, une possibilité est de calculer une approximation par boostrap paramétrique de la distribution de la déviance sous l’hypothèse que le modèle est adéquat.
Attention : même en se limitant à \(B=1000\) réplications bootstrap, le temps
de calcul est élevé. L’utilisation de fonctions s’exécutant plus
rapidement que la fonction glm du paquet
stats, par exemple la fonction speedglm du
paquet speedglm, pourrait être envisagée afin de réduire le
temps de calcul.
dev_boot <- param_boot_glm(ajustement, B = 1000, f = deviance)On peut alors comparer la loi du \(\chi^2\) utilisée précédemment (dont la densité est représentée par la courbe rouge), et la loi empirique obtenue avec les réplications bootstrap (représentée par un histogramme). Le résultat est sans appel : l’approximation est tellement inexacte qu’elle est inutilisable pour effectuer un test.
hist(dev_boot, xlim = c(218000, 245000), prob = TRUE)
abscisses <- seq(218000, 245000, length.out = 200)
lines(x = abscisses, y = dchisq(x = abscisses, df = deg_lib_res), col = "red")
En revanche, nous pouvons calculer une \(p-\)valeur approchée, a priori plus correcte, grâce aux réplications bootstrap.
p_val_boot <- 1 - ecdf(dev_boot)(dev)
p_val_boot## [1] 0.071La \(p-\)valeur ainsi estimée ne suggère pas de problème majeur d’adéquation du modèle. (Il conviendrait de prendre également en compte l’impact du nombre limité de réplications bootstrap sur la précision d’estimation de la \(p-\)valeur, ce que l’on néglige ici.)
Il existe d’autres approches (moins coûteuses en temps de calcul) pour effectuer un test formel d’adéquation à partir de la déviance. Ci-après, on utilise celle décrite dans [McCullagh, “The Conditional Distribution of Goodness-of-Fit Statistics for Discrete Data”, Journal of the American Statistical Association 81(393):104-107, March 1986]. Ce test repose sur une approximation normale de la loi conditionnelle de la déviance sachant la valeur de \(\sum_{i=1}^N x^{(j)}_i Y_i\), pour \(j=0,\ldots, d\).
dev_test_pois_McCullagh1986(ajustement)## $p.valeur
## [1] 0.08538501
##
## $stat_test
## [1] 1.720261Ici encore, la \(p-\)valeur obtenue ne suggère pas de problème majeur d’adéquation.
1.3 Test d’adéquation basé sur la déviance, données groupées
On effectue de nouveau les étapes précédentes, cette fois sur les données groupées par valeurs communes des variables explicatives. Même si, comme c’est le cas dans notre exemple, on dispose des données individuelles, il n’est pas inintéressant de tester également l’adéquation sur les données groupées. [A détailler…]
donnees_groupe <- aggregate(cbind(Claims, Insured) ~ Kilometres + Bonus + Zone + Make,
FUN = sum,
data = donnees)
donnees <- donnees_groupe
ajustement <- glm(Claims ~ Kilometres + Zone + Bonus + Make + offset(log(Insured)),
family = poisson(link = "log"),
data = donnees)On commence donc par le test basé sur l’approximation (a priori douteuse) de la loi de la déviance par une loi du \(\chi^2\) à \(N-(d+1)\) degrés de liberté.
deg_lib_res <- df.residual(ajustement)
deg_lib_res## [1] 1680dev <- deviance(ajustement)
dev## [1] 1683.414p_val_chi2 <- 1 - pchisq(q = dev, df = deg_lib_res)
p_val_chi2## [1] 0.4719564La \(p-\)valeur ainsi obtenue ne suggère pas de problème d’adéquation, mais nous avons vu que l’approximation de la loi de la déviance utilisée par ce test pouvait être complètement erronée. On s’attend néanmoins à une situation moins extrême ici, du fait que les données sont groupées.
Pour préciser cela, on repasse par des réplications bootstrap. (On effectue ici \(B=10000\) simulations, ce qui est praticable du fait que la taille du jeu de données a été considérablement réduite par le groupement.)
dev_boot <- param_boot_glm(ajustement, B = 10000, f = deviance)On constate que, cette fois, l’approximation par la loi du \(\chi^2\) n’est pas catastrophique, quoiqu’inexacte.
hist(dev_boot, prob = TRUE)
abscisses <- seq(1400, 2000, length.out = 200)
lines(x = abscisses, y = dchisq(x = abscisses, df = deg_lib_res), col = "red")
p_val_boot <- 1 - ecdf(dev_boot)(dev)
p_val_boot## [1] 0.3972La \(p-\)valeur approchée obtenue avec les réplications bootstrap diffère de celle obtenue avec l’approximation par la loi du \(\chi^2\), mais conduit à une conclusion identique : pas de problème d’adéquation détecté. Au vu de ce qui précède, on déconseille néanmoins l’utilisation du test basé sur l’approximation par une loi du \(\chi^2\) à \(N-(d+1)\) degrés de liberté.
On utilise de nouveau le test proposé par [McCullagh, 1986], qui ne suggère pas non plus de problème d’adéquation du modèle.
dev_test_pois_McCullagh1986(ajustement)## $p.valeur
## [1] 0.8846602
##
## $stat_test
## [1] -0.1450641Ici encore, la \(p-\)valeur obtenue ne suggère pas de problème majeur d’adéquation.
2 Exemple : régression logistique
2.1 Chargement des données et ajustement du modèle
On utilise ici un jeu de données et un modèle de régression logistique étudié en détail dans [Hosmer, Lemeshow, Sturdivant, “Applied logistic regression”, 2013].
load(url("https://irma.math.unistra.fr/~jberard/GLOW500"))
ajustement <- glm(FRACTURE ~
AGE +
HEIGHT +
PRIORFRAC +
MOMFRAC +
ARMASSIST +
RATERISK_supegal_3 +
AGE:PRIORFRAC +
MOMFRAC:ARMASSIST,
family = binomial(link = "logit"),
data = donnees)
summary(ajustement)##
## Call:
## glm(formula = FRACTURE ~ AGE + HEIGHT + PRIORFRAC + MOMFRAC +
## ARMASSIST + RATERISK_supegal_3 + AGE:PRIORFRAC + MOMFRAC:ARMASSIST,
## family = binomial(link = "logit"), data = donnees)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.71724 3.32173 0.517 0.605175
## AGE 0.05731 0.01650 3.473 0.000515 ***
## HEIGHT -0.04674 0.01833 -2.550 0.010781 *
## PRIORFRAC1 4.61229 1.88018 2.453 0.014163 *
## MOMFRAC1 1.24664 0.39296 3.172 0.001512 **
## ARMASSIST1 0.64416 0.25193 2.557 0.010561 *
## RATERISK_supegal_3TRUE 0.46897 0.24078 1.948 0.051449 .
## AGE:PRIORFRAC1 -0.05527 0.02593 -2.132 0.033044 *
## MOMFRAC1:ARMASSIST1 -1.28054 0.62299 -2.055 0.039834 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 562.34 on 499 degrees of freedom
## Residual deviance: 500.50 on 491 degrees of freedom
## AIC: 518.5
##
## Number of Fisher Scoring iterations: 42.2 Test d’adéquation basé sur la statistique de Pearson
On commence par appliquer l’approximation douteuse de la loi de la statistique de Pearson par une loi du \(\chi^2\) à \(N-(d+1)\) degrés de liberté.
deg_lib_res <- df.residual(ajustement)
deg_lib_res## [1] 491stat_P <- sum(residuals(ajustement, type = "pearson")^2)
stat_P## [1] 496.2274p_val_chi2 <- 1 - pchisq(q = stat_P, df = deg_lib_res)
p_val_chi2## [1] 0.4256416La \(p-\)valeur ainsi obtenue ne suggère pas de problème d’adéquation, mais, comme précédemment, on va utiliser l’approche du bootstrap paramétrique pour vérifier (ou infirmer) l’approximation utilisée pour effectuer ce test.
# Bootstrap paramétrique pour la statistique de Pearson
stat_P_boot <- param_boot_glm(ajustement,
B = 10000,
f = function(x) {sum(residuals(x,type="pearson")^2)})La comparaison graphique effectuée ci-après entre la densité de la loi du \(\chi^2\) à \(N-(d+1)\) degrés de liberté et la loi empirique obtenue à l’aide des réplications bootstrap illustre de nouveau le manque de fiabilité de cette approximation lorsque l’on n’est pas dans la situation où chaque ligne du jeu de données représente un volume d’exposition important.
hist(stat_P_boot, prob = TRUE, breaks = 50)
abscisses <- seq(0, 2 * deg_lib_res, length.out = 400)
lines(x = abscisses, y = dchisq(x = abscisses, df = deg_lib_res), col = "red")
La \(p-\)valeur (a priori plus correcte) obtenue par bootstrap paramétrique ne suggère pas de problème d’adéquation.
p_val_boot <- 1 - ecdf(stat_P_boot)(stat_P) # p-valeur par bootstrap paramétrique
p_val_boot## [1] 0.6212Une alternative à l’approche bootstrap est l’approximation proposée par [G. Osius, D. Rojek (1992) “Normal goodness-of-fit tests for multinomial models with large degrees-of-freedom.”, Journal of the American Statistical Association 87: 1145–1152.], qui fournit une valeur comparable à celle obtenue précédemment.
test_Osius_Rojek_1992(ajustement)## $p.valeur
## [1] 0.6518748
##
## $stat.test
## [1] 0.45115922.3 Test d’adéquation de Hosmer-Lemeshow
On illustre maintenant l’utilisation du test d’adéquation de Hosmer-Lemeshow, basé sur un groupement en classes par valeur croissante du risque modélisé.
test_HL <- test_Hosmer_Lemeshow(ajustement) # Test de Hosmer-Lemeshow avec 10 groupes (par défaut)
test_HL## $p.valeur
## [1] 0.5704067
##
## $stat.test
## [1] 6.690106Pour vérifier la qualité de l’approximation utilisée par ce test (qui, pour \(K\) groupes, utilise une loi du \(\chi^2\) à \(K-2\) degrés de liberté pour la statistique de test), on effectue des simulations bootstrap.
stat_HL_boot <- param_boot_glm(ajustement,
B = 10000,
f = function(x) {test_Hosmer_Lemeshow(x)$stat.test})On peut maintenant comparer la distribution obtenue avec une loi du \(\chi^2\) à \(10-2=8\) degrés de liberté. On constate une approximation très satisfaisante.
hist(stat_HL_boot, prob = TRUE, breaks = 50)
abscisses <- seq(0, 30, length.out = 400)
lines(x = abscisses, y = dchisq(x = abscisses, df = 8), col = "red")
La \(p-\)valeur obtenue à l’aide des réplications bootstrap est voisine de celle obtenue par l’approximation par une loi du \(\chi^2\) à \(K-2\) degrés de liberté.
p_val_boot <- 1 - ecdf(stat_HL_boot)(test_HL[["stat.test"]]) # p-valeur par bootstrap paramétrique
p_val_boot## [1] 0.5615