Petit complément par rapport au cours de mardi. On avait évoqué tout d’abord la lecture des sorties lorsque l’on régresse sur des variables catégorielles (des facteurs). Commençons par supprimer la constante de la régression
> reg0=glm(nbre~0+zone,offset=log(exposition),data=base, + family=poisson(link="log")) > summary(reg0) Call: glm(formula = nbre ~ 0 + zone, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3968 -0.2996 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) zoneB -2.54187 0.06287 -40.43 <2e-16 *** zoneA -2.54912 0.05285 -48.23 <2e-16 *** zoneC -2.38525 0.03753 -63.56 <2e-16 *** zoneD -2.13454 0.03878 -55.05 <2e-16 *** zoneE -2.00204 0.03965 -50.49 <2e-16 *** zoneF -2.06932 0.11547 -17.92 <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: 50966 on 50000 degrees of freedom Residual deviance: 15692 on 49994 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6 > predict(reg0,newdata=data.frame( + zone=c("A","B","C","D","E"),exposition=rep(1,5))) 1 2 3 4 5 -2.549120 -2.541870 -2.385253 -2.134543 -2.002044
On voit que toutes les modalités sont présentes, et toutes sont significatives. Si on régresse sur la constante, il faudra supprimer une modalité pour rendre le modèle identifiable. On peut forcer pour que la modalité de référence soit la seconde,
> base$zone=relevel(base$zone,"B") > regB=glm(nbre~zone,offset=log(exposition),data=base, + family=poisson(link="log")) > summary(regB) Call: glm(formula = nbre ~ zone, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3968 -0.2996 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.54187 0.06287 -40.431 < 2e-16 *** zoneA -0.00725 0.08213 -0.088 0.929661 zoneC 0.15662 0.07322 2.139 0.032432 * zoneD 0.40733 0.07387 5.514 3.50e-08 *** zoneE 0.53983 0.07433 7.263 3.80e-13 *** zoneF 0.47255 0.13148 3.594 0.000325 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49994 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6 > predict(regB,newdata=data.frame( + zone=c("A","B","C","D","E"),exposition=rep(1,5))) 1 2 3 4 5 -2.549120 -2.541870 -2.385253 -2.134543 -2.002044
On notera que les prédictions ne changent pas. On peut aussi choisir la première comme modalité de référence,
> base$zone=relevel(base$zone,"A") > reg=glm(nbre~zone,offset=log(exposition), > data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zone, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3968 -0.2996 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.54912 0.05285 -48.232 < 2e-16 *** zoneB 0.00725 0.08213 0.088 0.929661 zoneC 0.16387 0.06482 2.528 0.011471 * zoneD 0.41458 0.06555 6.324 2.54e-10 *** zoneE 0.54708 0.06607 8.280 < 2e-16 *** zoneF 0.47980 0.12699 3.778 0.000158 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49994 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6
Le fait que la seconde modalité ne soit pas significative se lit par rapport à la modalité de référence (en l’occurrence la première): non significatif signifie alors non significativement différente. Autrement dit, on peut regrouper les modalités en une seule.
> base$zonesimple=base$zone > base$zonesimple[base$zone%in%c("A","B")]="A" > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3959 -0.2989 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.54612 0.04046 -62.937 < 2e-16 *** zonesimpleC 0.16087 0.05518 2.915 0.00355 ** zonesimpleD 0.41158 0.05604 7.345 2.06e-13 *** zonesimpleE 0.54408 0.05665 9.605 < 2e-16 *** zonesimpleF 0.47681 0.12235 3.897 9.74e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49995 degrees of freedom AIC: 20798 Number of Fisher Scoring iterations: 6
On note qu’avec ce regroupement, les autres modalités sont sensiblement différentes. On peut aussi retenir la troisième comme modalité de référence
> base$zonesimple=relevel(base$zonesimple,"C") > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3959 -0.2989 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.38525 0.03753 -63.557 < 2e-16 *** zonesimpleA -0.16087 0.05518 -2.915 0.00355 ** zonesimpleD 0.25071 0.05396 4.646 3.39e-06 *** zonesimpleE 0.38321 0.05460 7.019 2.24e-12 *** zonesimpleF 0.31593 0.12142 2.602 0.00927 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49995 degrees of freedom AIC: 20798 Number of Fisher Scoring iterations: 6
Comme toutes les modalités semblent significatives, on peut tenter de prendre comme modalité de référence une des dernières (dont les estimations des coefficients donnent des résultats très proches)
> base$zonesimple=relevel(base$zonesimple,"F") > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5717 -0.3959 -0.2989 -0.1547 12.6722 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.06932 0.11547 -17.921 < 2e-16 *** zonesimpleC -0.31593 0.12142 -2.602 0.00927 ** zonesimpleA -0.47681 0.12235 -3.897 9.74e-05 *** zonesimpleD -0.06522 0.12181 -0.535 0.59232 zonesimpleE 0.06727 0.12209 0.551 0.58161 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 15809 on 49999 degrees of freedom Residual deviance: 15692 on 49995 degrees of freedom AIC: 20798 Number of Fisher Scoring iterations: 6
Au vue de cette dernière sortie, on peut tenter de fusionner toutes les dernières classes ensembles
> base$zonesimple[base$zone%in%c("D","E","F")]="F" > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > summary(reg) Call: glm(formula = nbre ~ zonesimple, family = poisson(link = "log"), data = base, offset = log(exposition)) Deviance Residuals: Min 1Q Median 3Q Max -0.5660 -0.3959 -0.3004 -0.1547 12.5929 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.07182 0.02696 -76.853 < 2e-16 *** zonesimpleC -0.31344 0.04621 -6.783 1.18e-11 *** zonesimpleA -0.47431 0.04861 -9.757 < 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: 15809 on 49999 degrees of freedom Residual deviance: 15698 on 49997 degrees of freedom AIC: 20800 Number of Fisher Scoring iterations: 6
Bon, formellement, regrouper deux modalités (i.e. décréter que deux variables sont non significatives simultanément) demande un peu plus qu’un test de Student, ou que deux tests de Student…. Si on remonte un peu en arrière, on aurait pu faire un test multiple avant de fusionner les trois modalités (un test de type Fisher, ou une ANOVA)
> base$zonesimple=relevel(base$zonesimple,"F") > reg=glm(nbre~zonesimple,offset=log(exposition), + data=base,family=poisson(link="log")) > library(car) > linearHypothesis(reg,c("zonesimpleD=0","zonesimpleE=0")) Linear hypothesis test Hypothesis: zonesimpleD = 0 zonesimpleE = 0 Model 1: restricted model Model 2: nbre ~ zonesimple Res.Df Df Chisq Pr(>Chisq) 1 49997 2 49995 2 5.7073 0.05763 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Manifestement, on peut accepter l’hypothèse que ces trois catégories n’en font qu’une. La zone géographique peut alors se découper en trois grandes zones, et pas en six. On notera que cela correspond à ce que propose un arbre de régression
