1 Test de Wald

On commence par un jeu de données simulé : \(100\) observations conjointes de deux variables explicatives numériques x1 et x2, et une variable réponse y dichotomique \(0/1\)

load(url("https://irma.math.unistra.fr/~jberard/reg_log_sim_1")) 
str(donnees)
## 'data.frame':    100 obs. of  3 variables:
##  $ y : int  1 1 0 1 1 0 0 1 1 1 ...
##  $ x1: num  3.53 -2.21 -1.14 -2.03 1.21 ...
##  $ x2: num  3.3 -1.87 -1 -2.11 1.08 ...

On ajuste maintenant un modèle de régression logistique simple sur le jeu de données.

ajustement_100 <- glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), data = donnees)

La commande summary nous permet d’accéder aux résultats de tests de Wald de nullité des coefficients, les \(p-\)valeurs étant indiquées dans la colonne Pr(>|z|) pour chacun des coefficients.

Bien entendu, aucun coefficient n’est estimé à une valeur exactement égale à \(0\). La question que l’on se pose ici est de préciser dans quelle mesure cette valeur estimée non-nulle estimée est le signe d’un coefficient dont la vraie valeur est non-nulle, ou si cette valeur estimée non-nulle n’est que le résultat d’une fluctuation autour d’une vraie valeur égale à \(0\). On constate ici que :

  • l’hypothèse de nullité du coefficient associé au terme constant (Intercept) est assez nettement rejetée, avec une \(p-\)valeur à un peu moins de \(0.02\), ce qui lui vaut d’être assortie de deux étoiles **; l’écart à \(0\) de la valeur estimée est donc jugé “statistiquement significatif”;

  • l’hypothèse de nullité du coefficient associé à la variable explicative numérique x1 n’est pas rejetée (quel que soit le niveau “raisonnable” de risque de première espèce que l’on se fixe), avec une \(p-\)valeur d’environ \(0.89\); l’écart à \(0\) de la valeur estimée est donc jugé “statistiquement non-significatif”;

  • de même, l’hypothèse de nullité du coefficient associé à la variable explicative numérique x2 n’est pas rejetée (quel que soit le niveau “raisonnable” de risque de première espèce que l’on se fixe), avec une \(p-\)valeur d’environ \(0.27\).

summary(ajustement_100)
## 
## Call:
## glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), 
##     data = donnees)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   1.6838     0.5444   3.093  0.00198 **
## x1            0.1684     1.1828   0.142  0.88676   
## x2            1.5743     1.4144   1.113  0.26571   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 136.058  on 99  degrees of freedom
## Residual deviance:  47.701  on 97  degrees of freedom
## AIC: 53.701
## 
## Number of Fisher Scoring iterations: 7

Le jeu de données précédent est en fait extrait d’un jeu de données simulé plus vaste comportant \(10000\) observations.

load(url("https://irma.math.unistra.fr/~jberard/reg_log_sim_2")) 
str(donnees)
## 'data.frame':    10000 obs. of  3 variables:
##  $ y : int  1 1 0 1 1 0 0 1 1 1 ...
##  $ x1: num  3.53 -2.21 -1.14 -2.03 1.21 ...
##  $ x2: num  3.3 -1.87 -1 -2.11 1.08 ...

On procède de nouveau à l’ajustement d’un modèle de régression logistique simple, sur ce jeu de données plus étendu.

ajustement_10000 <- glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), data = donnees)
summary(ajustement_10000)
## 
## Call:
## glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), 
##     data = donnees)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.959080   0.033418  28.700   <2e-16 ***
## x1          0.977406   0.106204   9.203   <2e-16 ***
## x2          0.001149   0.115076   0.010    0.992    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 13360.5  on 9999  degrees of freedom
## Residual deviance:  6969.6  on 9997  degrees of freedom
## AIC: 6975.6
## 
## Number of Fisher Scoring iterations: 6

On constate que :

  • l’hypothèse de nullité du coefficient associé au terme constant est fortement rejetée, avec une \(p-\)valeur inférieure à \(2 \cdot 10^{-16}\), ce qui lui vaut d’être assortie de trois étoiles ***;

  • Contrairement à l’ajustement précédent, l’hypothèse de nullité du coefficient associé à la variable explicative numérique x1 est fortement rejetée, avec une \(p-\)valeur inférieure à \(2 \cdot 10^{-16}\);

  • l’hypothèse de nullité du coefficient associé à la variable explicative numérique x2 n’est pas rejetée (quel que soit le niveau “raisonnable” de risque de première espèce que l’on se fixe), avec une \(p-\)valeur d’environ \(0.99\).

Les deux jeux de données étudiés ici ont été simulés à l’aide du même modèle, à savoir un modèle de régression logistique pour la variable réponse \(y\), dont la composante systématique est de la forme \(\eta=1+x^{(1)}\). Le modèle de régression logistique est donc ici (par construction) valide, avec les vraies valeurs des coefficients données par \(\beta_0=1\), \(\beta_1=1\), \(\beta_2=0\).

Dans le premier cas (100 lignes), la non-nullité de \(\beta_0\) et la nullité de \(\beta_2\) sont correctement identifiées, mais pas la non-nullité de \(\beta_1\), du fait d’un nombre de données insuffisant. Dans le deuxième cas (10000 lignes au lieu de 100), la nullité/non-nullité des divers coefficients est correctement détectée, du fait de la taille plus grande du deuxième jeu de données (10000 lignes au lieu de 100), qui permet d’obtenir une estimation plus précise des coefficients, et en particulier une discrimination plus précise entre coefficients nuls et non-nuls.

On charge maintenant un troisième jeu de données, dont la structure est analogue aux deux précédents.

load(url("https://irma.math.unistra.fr/~jberard/reg_log_sim_3")) 
str(donnees)
## 'data.frame':    1000 obs. of  3 variables:
##  $ y : int  1 1 0 1 0 1 1 1 1 0 ...
##  $ x1: num  1.5 6.3 -2.17 6.76 -3.73 ...
##  $ x2: num  -2.0238 -4.3049 0.0014 -7.0384 1.5016 ...

On ajuste ensuite deux modèles de régression logistique : l’un comportant uniquement la variable explicative x2, et l’autre comportant les deux variables explicatives x1 et x2.

# Ajustement avec x2 toute seule
ajustement_1 <- glm(formula = y ~ x2, family = binomial(link = "logit"), data = donnees)
summary(ajustement_1)
## 
## Call:
## glm(formula = y ~ x2, family = binomial(link = "logit"), data = donnees)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.60266    0.08315   7.248 4.23e-13 ***
## x2          -0.59040    0.03946 -14.963  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1345.2  on 999  degrees of freedom
## Residual deviance:  940.9  on 998  degrees of freedom
## AIC: 944.9
## 
## Number of Fisher Scoring iterations: 5
# Ajustement avec x1 et x2
ajustement_2 <- glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), data = donnees)
summary(ajustement_2)
## 
## Call:
## glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), 
##     data = donnees)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.87531    0.10428   8.394   <2e-16 ***
## x1           1.08953    0.08708  12.513   <2e-16 ***
## x2           0.06773    0.06115   1.108    0.268    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1345.21  on 999  degrees of freedom
## Residual deviance:  692.67  on 997  degrees of freedom
## AIC: 698.67
## 
## Number of Fisher Scoring iterations: 6

Dans le premier cas (ajustement_1), on constate que l’hypothèse de nullité du coefficient associé à la variable x2 est fortement rejetée, avec une \(p-\)valeur inférieure à \(2 \cdot 10^{-16}\).

Dans le deuxième cas (ajustement_2), on constate que l’hypothèse de nullité du coefficient associé à la variable x2 n’est pas rejetée (quel que soit le niveau “raisonnable” de risque de première espèce que l’on se fixe), avec une \(p-\)valeur d’environ \(0.27\), tandis que l’hypothèse de nullité du coefficient associé à la variable x1 est fortement rejetée, avec une \(p-\)valeur inférieure à \(2 \cdot 10^{-16}\).

La valeur estimée des coefficients et leur significativité statistique éventuelle dépendent des variables incluses dans le modèle !!!

Ici, le modèle utilisé pour simuler les données est exactement le même que dans l’exemple précédent : modèle de régression logistique pour la variable réponse \(y\), dont la composante systématique est de la forme \(\eta=1+x^{(1)}\). Même si la valeur de x2 n’intervient aucunement dans la loi de la valeur simulée de y, le modèle dans lequel x2 est la seule variable explicative disponible, et la corrélation entre les valeurs de x1 et x2, expliquent le fait qu’un coefficient non-nul soit identifié pour x2. Lorsque les deux variables x1 et x2 sont utilisées comme variables explicatives, la valeur nulle de \(\beta_2\) et la valeur non-nulle de \(\beta_1\) sont correctement identifiées.

2 Test du rapport de vraisemblance (aka LRT : “Likelihood Ratio Test”)

On charge un jeu de données simulé comportant \(1000\) observations conjointes d’une variable réponse y dichotomique \(0/1\), d’une variable explicative numérique x1, et d’une variable explicative catégorielle x2 (à trois catégories : a, b, c)

 load(url("https://irma.math.unistra.fr/~jberard/reg_log_sim_4")) 
str(donnees)
## 'data.frame':    1000 obs. of  3 variables:
##  $ y : int  1 1 0 1 0 1 1 1 1 0 ...
##  $ x1: num  1.5 6.3 -2.17 6.76 -3.73 ...
##  $ x2: Factor w/ 3 levels "a","b","c": 3 3 2 3 1 3 3 2 3 1 ...

On compare ici deux ajustements d’un modèle de régression logistique : l’un comportant uniquement la variable x1, et l’autre comportant les deux variables x1 et x2. Le modèle comportant uniquement la variable x1 peut être vu comme une restriction du modèle comportant les deux variables x1 et x2 dans lequel on impose une valeur nulle aux coefficients associés à la variable x2.

ajustement_1 <- glm(formula = y ~ x1, family = binomial(link = "logit"), data = donnees) 
summary(ajustement_1)
## 
## Call:
## glm(formula = y ~ x1, family = binomial(link = "logit"), data = donnees)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.39329    0.11141   12.51   <2e-16 ***
## x1           0.84534    0.05514   15.33   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1255.24  on 999  degrees of freedom
## Residual deviance:  739.37  on 998  degrees of freedom
## AIC: 743.37
## 
## Number of Fisher Scoring iterations: 6
ajustement_2 <- glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), data = donnees)
summary(ajustement_2)
## 
## Call:
## glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), 
##     data = donnees)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.06044    0.20224  10.188  < 2e-16 ***
## x1           1.07027    0.07819  13.688  < 2e-16 ***
## x2b          0.11532    0.32704   0.353    0.724    
## x2c         -1.30227    0.27447  -4.745 2.09e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1255.24  on 999  degrees of freedom
## Residual deviance:  711.21  on 996  degrees of freedom
## AIC: 719.21
## 
## Number of Fisher Scoring iterations: 6

Pour effectuer le test du rapport de vraisemblance, on extrait les log-vraisemblances maximisées obtenues pour chacun des deux ajustements, à l’aide de la fonction logLik.

L1<-logLik(ajustement_1)
L1
## 'log Lik.' -369.6841 (df=2)
L2<-logLik(ajustement_2)
L2
## 'log Lik.' -355.6034 (df=4)

On calcule ensuite la différence des nombres de degré de liberté entre le “grand” modèle (avec les variables x1 et x2) et le “petit” modèle (avec uniquement la variable x1). On constate que cette différence est égale à \(2\), traduisant la contrainte selon laquelle les deux coefficients associés à la variable x2 dans le “grand” modèle sont contraints de valoir \(0\) dans le “petit” modèle.

diff_deg_lib <- attributes(L2)[["df"]] - attributes(L1)[["df"]]
diff_deg_lib
## [1] 2

On calcule ensuite la statistique du test du rapport de vraisemblance, c’est-à-dire le double de la différence des log-vraisemblances maximisées entre le “grand” modèle et le “petit” modèle.

LRT_stat <- 2 * (as.double(L2) - as.double(L1))
LRT_stat
## [1] 28.16134

La \(p-\)valeur du test est ensuite calculée à l’aide de la fonction de répartition de la loi du chi-deux à diff_deg_lib degrés de liberté.

LRT_pval <- 1 - pchisq(q = LRT_stat, df = diff_deg_lib) 
LRT_pval
## [1] 7.67083e-07

Le test du rapport de vraisemblance penche très clairement en faveur de l’inclusion de la variable x2 : l’hypothèse de nullité des coefficients associés à la variable x2 est nettement rejetée, avec une \(p-\)valeur d’environ \(7.7 \cdot 10^{-7}\).

Dans ce qui précède, on a décortiqué les différentes étapes conduisant au test du rapport de vraisemblance. On peut obtenir le même résultat de manière plus directe en utilisant la fonction anova (qui sera donc ici anova.glm), comme illustré ci-dessous.

anova(ajustement_1, ajustement_2, test = "LRT")

On utilise maintenant le test du rapport de vraisemblance dans un autre but : tester la pertinence d’une fusion des modalités a et b de la variable x2 dans le modèle utilisant les deux variables x1 et x2.

On effectue ceci de manière “manuelle”, en créant une nouvelle variable x2_fus dans laquelle ces deux modalités sont fusionnées en une seule modalité nommée ab, à l’aide de la commande revalue du paquet plyr.

donnees[["x2_fus"]] <- forcats::fct_recode(.f = donnees[["x2"]], "ab" = "a", "ab" = "b")

On ré-ajuste maintenant notre modèle de régression logistique en utilisant x2_fus au lieu de x2.

ajustement_3 <- glm(formula = y ~ x1 + x2_fus, family = binomial(link = "logit"), data = donnees)
summary(ajustement_3)
## 
## Call:
## glm(formula = y ~ x1 + x2_fus, family = binomial(link = "logit"), 
##     data = donnees)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.08830    0.18674   11.18  < 2e-16 ***
## x1           1.07517    0.07702   13.96  < 2e-16 ***
## x2_fusc     -1.33200    0.26167   -5.09 3.57e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1255.24  on 999  degrees of freedom
## Residual deviance:  711.33  on 997  degrees of freedom
## AIC: 717.33
## 
## Number of Fisher Scoring iterations: 6

On effectue maintenant un test du rapport de vraisemblance pour comparer les deux modèles ainsi ajustés. Le “grand” modèle est ici celui dans lequelle les modalités a et b de la variable x2 ne sont pas fusionnées, et le “petit” modèle celui dans lequel elles le sont, cette fusion pouvant se traduire par la contrainte d’égalité des coefficients associés aux modalités a et b de la variable x2.

anova(ajustement_3, ajustement_2, test = "LRT")

On constate que le test du rapport de vraisemblance penche ici très nettement en faveur de la fusion des modalités a et b. L’hypothèse d’égalité des coefficients associés aux modalités a et b n’est pas rejetée, avec une \(p-\)valeur d’environ \(0.72\).

Dans les deux exemples précédents, le modèle utilisé pour simuler les données est un modèle de régression logistique pour la variable réponse \(y\), dont la composante systématique est de la forme \(\eta=1+x^{(1)}+\mathbf{1}(x^{(2)} \in \{a,b\})\). Les deux tests effectués ont donc correctement identifié, respectivement, l’existence d’un effet pour x2, et l’absence de différence d’effet entre les modalités a et b de x2.

On charge maintenant un nouveau jeu de données, comportant \(1000\) observations conjointes d’une variable réponse y dichotomique \(0/1\), d’une variable explicative logique x1, et d’une variable explicative numérique x2.

load(url("https://irma.math.unistra.fr/~jberard/reg_log_sim_5")) 
str(donnees)
## 'data.frame':    1000 obs. of  3 variables:
##  $ y : int  1 1 1 0 1 1 1 1 0 1 ...
##  $ x1: Factor w/ 2 levels "FALSE","TRUE": 1 1 1 1 2 2 1 1 1 1 ...
##  $ x2: num  3.66 3.6 1.44 -1.41 5.4 ...

On ajuste cette fois deux modèles de régression logistique : l’un encodant un effet purement linéaire de la variable x1 sur la composante systématique (donc au moyen d’un seul coefficient), l’autre encodant l’effet de cette variable par un polynôme de degré \(2\), avec donc un possible effet quadratique.

ajustement_1 <- glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), data = donnees) 
ajustement_2 <- glm(formula = y ~ x1 + poly(x2, degree = 2), 
                    family = binomial(link = "logit"), 
                    data = donnees)

On effectue maintenant un test du rapport de vraisemblance pour comparer les deux modèles ainsi ajustés. Le “grand” modèle est ici celui dans lequelle un effet quadratique est ajusté pour x2, et le “petit” modèle celui dans lequel l’effet de x2 est purement linéaire, cette restriction pouvant se traduire par la nullité du coefficient associé au carré de x2 dans le modèle avec polynôme de degré \(2\).

anova(ajustement_1, ajustement_2, test = "LRT")

On constate que le test du rapport de vraisemblance penche ici très nettement en faveur de l’inclusion d’un effet quadratique. L’hypothèse de nullité du coefficient correspondant est rejetée, avec une \(p-\)valeur d’environ \(8 \cdot 10^{-6}\).

Dans l’exemple précédent, le modèle utilisé pour simuler les données est un modèle de régression logistique pour la variable réponse \(y\), dont la composante systématique est de la forme \(\eta=1+\mathbf{1}(x^{(1)}= TRUE)+x^{(2)}+0.008 \cdot (x^{(2)})^2\). Le test effectué a donc correctement identifié la présence d’un effet quadratique pour x2.

3 Test du score

On reprend simplement l’exemple précédent, pour illustrer le test du score, qui donne lieu à une conclusion similaire. On utilise de nouveau la fonction anova, mais avec l’argument test="Rao". (Voir également la fonction glmscoretest du paquet statmod, et la fonction test_score du paquet outilscoursglm.)

anova(ajustement_1, ajustement_2, test = "Rao")

4 Un exemple sur des données réelles

On commence par récupérer le jeu de données. Il s’agit d’un jeu de données d’assurance automobile, adapté du paquet CASdatasets, comportant les variables suivantes.

  • DriverAge : âge du conducteur, 5 catégories d’âge
  • CarAge : âge du véhicule, 2 catégories d’âge
  • Density : densité de la zone, 5 catégories
  • Brand : marque du véhicule, 2 catégories
  • Power : puissance du véhicule, 3 catégories
  • Gas : type de motorisation, 2 catégories
  • ClaimNb : nombre de sinistres observés
  • Exposition : en années-police

Les données ont été groupées par valeurs communes des variables explicatives (toutes catégorielles), si bien que l’on a en tout \(495\) lignes.

load(url("https://irma.math.unistra.fr/~jberard/freMTPLfreq_groupe"))
str(donnees)
## 'data.frame':    495 obs. of  8 variables:
##  $ DriverAge: Factor w/ 5 levels "(17,22]","(22,26]",..: 1 2 3 4 5 1 3 4 5 1 ...
##  $ CarAge   : Factor w/ 2 levels "[0,15]","(15,Inf]": 1 1 1 1 1 2 2 2 2 1 ...
##  $ Density  : Factor w/ 5 levels "[0,40]","(40,200]",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ Brand    : Factor w/ 2 levels "other","J": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Power    : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Gas      : Factor w/ 2 levels "Diesel","Regular": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ClaimNb  : num  5 8 21 43 1 0 1 2 0 11 ...
##  $ Exposure : num  34.7 66.2 402 519.4 23.6 ...

On ajuste maintenant un modèle log-Poisson pour le nombre de sinistres, en utiliant le logarithme de l’exposition comme terme d’offset.

ajustement <-
glm(formula = ClaimNb ~ DriverAge + CarAge + Density + Brand + Power + Gas + offset(log(Exposure)),
    data = donnees, 
    family = poisson(link = "log"))
summary(ajustement)
## 
## Call:
## glm(formula = ClaimNb ~ DriverAge + CarAge + Density + Brand + 
##     Power + Gas + offset(log(Exposure)), family = poisson(link = "log"), 
##     data = donnees)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -1.92779    0.04550 -42.373  < 2e-16 ***
## DriverAge(22,26]     -0.61540    0.04610 -13.349  < 2e-16 ***
## DriverAge(26,42]     -1.08342    0.03650 -29.680  < 2e-16 ***
## DriverAge(42,74]     -1.08009    0.03563 -30.315  < 2e-16 ***
## DriverAge(74,Inf]    -1.10306    0.05190 -21.252  < 2e-16 ***
## CarAge(15,Inf]       -0.24524    0.03083  -7.955 1.80e-15 ***
## Density(40,200]       0.18350    0.02676   6.856 7.09e-12 ***
## Density(200,500]      0.31562    0.02968  10.633  < 2e-16 ***
## Density(500,4.5e+03]  0.53311    0.02606  20.456  < 2e-16 ***
## Density(4.5e+03,Inf]  0.66589    0.03582  18.591  < 2e-16 ***
## BrandJ               -0.17933    0.02479  -7.234 4.68e-13 ***
## Power2                0.08689    0.02344   3.707  0.00021 ***
## Power3                0.23979    0.03012   7.961 1.70e-15 ***
## GasRegular           -0.17951    0.01671 -10.743  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2357.05  on 494  degrees of freedom
## Residual deviance:  563.82  on 481  degrees of freedom
## AIC: 2168.8
## 
## Number of Fisher Scoring iterations: 4

On s’intéresse en particulier ici à la variable DriverAge.

On vérifie sur le tableau ci-dessus que
tous les coefficients estimés pour cette variable (à l’exception de celui associé à la modalité de référence, bien entendu) possèdent des valeurs dont la différence vis-à-vis de \(0\) est jugée statistiquement significative, d’où les *** pour chacun des coefficients.

Néanmoins, le tracé effectué ci-dessous (au moyen de la fonction trace_terme_univar du fichier d’outils) montre que les différences entre les trois coefficients estimés DriverAge(26,42], DriverAge(42,74],
DriverAge(74,Inf] sont faibles en comparaison de l’incertitude portant sur ces estimations (traduite par le tracé des intervalles de confiance à \(95\)% obtenus par approximation normale, pour chaque coefficient).

outilscoursglm::trace_terme_univar(ajustement = ajustement,
                                   nom_variable = "DriverAge",
                                   nom_terme = "DriverAge",
                                   trace_exposition = TRUE,
                                   nom_exposition = "Exposure", 
                                   trace_int_conf = TRUE)
## Warning in outilscoursglm::trace_terme_univar(ajustement = ajustement, nom_variable = "DriverAge", : Utilisation des données de ajustement$data

Par conséquent, la question se pose de savoir si les différences entre les valeurs estimées de ces trois coefficients traduisent une différence réelle, ou si elles ne sont que le reflet de fluctuations différentes autour d’une valeur identique.

On commence par créer une nouvelle variable fusionnant les trois dernières modalités de la variable DriverAge, et sans cette fusion.

# Fusion des trois dernières modalités de DriverAge
donnees[["DriverAge_fus"]] <- 
  forcats::fct_recode(.f = donnees[["DriverAge"]],
                      "(26,Inf]" = "(26,42]",
                      "(26,Inf]" = "(42,74]",
                      "(26,Inf]" = "(74,Inf]")

On ré-ajuste maintenant le modèle log-Poisson, en utilisant la variable ainsi fusionnée.

ajustement_fus <-
glm(formula = ClaimNb ~ DriverAge_fus + CarAge + Density + Brand + Power + Gas + offset(log(Exposure)),
    data = donnees, 
    family = poisson(link = "log"))

On peut maintenant juger de l’impact de la fusion en utilisant le test du rapport de vraisemblance.

anova(ajustement_fus, ajustement, test = "LRT")

Celui-ci penche visiblement en faveur de la fusion, avec une \(p-\)valeur à \(0.85\) environ. On n’a donc pas détecté de différence statistiquement significative entre les effets des trois dernières modalités de la variable DriverAge. (Bien entendu, la portée de cette conclusion dépend de la validité du modèle GLM utilisé ici, qui n’a pas été examinée. Tout ce que l’on peut dire à ce stade est que l’on n’a pas détecté de différence statistiquement significative dans le cadre du modèle utilisé.)

Pour faire bonne mesure, on applique également le test du score, qui, de manière rassurante, fournit une \(p-\)valeur très voisine, avec dans tous les cas une conclusion identique.

anova(ajustement_fus, ajustement, test = "Rao")

On peut également utiliser le test de Wald, qui nous permet d’illustrer l’utilisation plutôt commode de la fonction linearHypothesis du paquet car pour tester l’égalité entre coefficients. Le résultat obtenu est de nouveau très proche des deux précédents (en termes de \(p-\)valeur) et donne lieu à une conclusion identique.

library(car)
## Loading required package: carData
linearHypothesis(model = ajustement,
                 c("DriverAge(42,74]=DriverAge(26,42]","DriverAge(26,42]=DriverAge(74,Inf]"))