# Visualizing effects of a categorical explanatory variable in a regression

Recently, I’ve been working on two problems that might be related to semiotic issues in predictive modeling (i.e. instead of a standard regression table, how can we plot coefficient values in a regression model). To be more specific, I have a variable of interest $Y$ that is observed for several individuals $i$, with explanatory variables $\mathbf{x}_i$, year $t$, in a specific region $z_i\in\{A,B,C,D,E\}$. Suppose that we have a simple (standard) linear model (forget about time here) $$y_i=\beta_0+\beta_1x_{1,i}+\cdots+\beta_kx_{k,i}+\sum_j \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i$$

Let us forget the temporal effect to focus on the spatial effect today. And consider some simulated dataset. There will be only one (continuous) explanatory variable. And I will generate correlated covariates, just to be more realistic.

n=1000 library(mnormt) r=.5 Sigma=matrix(c(1,r,r,1), 2, 2) set.seed(1) X=rmnorm(n,c(0,0),Sigma) X1=cut(X[,1],c(-100,quantile(X[,1],c(.1,.4,.7,.85)), 100),labels=LETTERS[1:5]) X2=X[,2] Y=5+X[,1]-X[,2]+rnorm(n)/2 db=data.frame(Y,X1,X2)

Here we have $$y_i=\beta_0+\beta_1x_{1,i}+\sum_{j\in\{A,B,C,D,E\}} \alpha_j \mathbf{1}(z_i\in j)+\varepsilon_i$$ The goal here is to get to graph to visualize the vector $\hat\alpha=(\hat\alpha_A,\cdots,\hat\alpha_E)$. Let us run the linear regression

reg1=lm(Y~X1+X2,data=db) idx=which(substr(names(reg1$coefficients), 1,2)=="X1") v1=reg1$coefficients[idx] names(v1)=LETTERS[2:5] barplot(v1,col=rgb(0,0,1,.4))

Note that it is possible to add some sort of “confidence interval” to discuss significance (or to avoid to spend hours discussing differences in bar heights that are not significantly different)

library(Hmisc) sv1=summary(reg1)$coefficients[idx,2] (bp1=barplot(v1,ylim=range(c(0,v1+2*sv1)))) errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE) My main concern here is the “reference” that is considered. Should $A$ be the reference? Why not $B$ db$X1=relevel(db$X1,"B") reg1=lm(Y~X1+X2,data=db) idx=which(substr(names(reg1$coefficients),1,2)=="X1") v1=reg1$coefficients[idx] names(v1)=LETTERS[c(1,3:5)] library(Hmisc) sv1=summary(reg1)$coefficients[idx,2] (bp1=barplot(v1) errbar(bp1[,1],v1,v1-2*sv1,v1+2*sv1,add=TRUE)

Why not the smallest one? Why not the largest one?… What if there is no simple way to choose. Furthermore, let us get back to the original point, which is that there might be some temporal aspects. More precisely, we can have $\hat\alpha^{(t)}=(\hat\alpha_A^{(t)},\cdots,\hat\alpha_E^{(t)})$. If we have also $\hat\alpha^{(t+1)}$ and we get another plot, how do we interpret it. If for $E$ the bar is taller, it means that relative to $A$, the difference has increased. I have the feeling that the interpretation is more complicated because we do not see, on that graph, changes in $\hat\alpha^{(t)}_A$.

Let us try something else. First, let us get back to the original setting

db$X1=relevel(db$X1,"A")

Consider here the regression without the intercept, so that all values remain

reg1=lm(Y~0+X1+X2,data=db) idx=which(substr(names(reg1$coefficients),1,2)=="X1") v1=reg2$coefficients[idx] names(v1)=LETTERS[1:5] barplot(v1)

It can be hard to read, especially if $Y$ takes (very) large values, and you think that barplots should start at 0. But still, having those 5 values is nice. Why not rescale that graph?

A natural idea my be to consider the case where no spatial component is considered, and to look at the difference with that reference.

reg1=lm(Y~1+X2,data=db) reg2=lm(Y~0+X1+X2,data=db) idx=which(substr(names(reg2$coefficients),1,2)=="X1") v1=reg2$coefficients[idx] v2=v1-reg1$coefficients["(Intercept)"] barplot(v2,col=rgb(0,0,1,.4)) sv2=summary(reg2)$coefficients[idx,2] (bp2=barplot(v2,ylim=range(c(v2-2*sv2,v2+2*sv2)))) errbar(bp2[,1],v2,v2-2*sv2,v2+2*sv2,add=TRUE)

I like that graph, I should admit it. Now, I still have some remaining questions. For instance, can we insure that when only the intercept is considered, the value of $\hat\beta_0$ is somewhere between $\hat\beta_A,\cdots,\hat\beta_E$? Is it possible that $\hat\beta_A-\hat\beta_0,\cdots,\hat\beta_E-\hat\beta_0$ are all positive? In that case, I would find that hard to interpret.

Actually, if I really want values that can be seen as compared to some average, why not consider a (weighted) average of $\hat\beta_A,\cdots,\hat\beta_E$? (weights being here proportion in each class, in each region)

# Logistic regression and categorical covariates

A short post to get back – for my nonlife insurance course – on the interpretation of the output of a regression when there is a categorical covariate. Consider the following dataset

> db = read.table("http://freakonometrics.free.fr/db.txt",header=TRUE,sep=";")
> attach(db)
> tail(db)
Y       X1       X2 X3
995  1 4.801836 20.82947  A
996  1 9.867854 24.39920  C
997  1 5.390730 21.25119  D
998  1 6.556160 20.79811  D
999  1 4.710276 21.15373  A
1000 1 6.631786 19.38083  A

Let us run a logistic regression on that dataset

> reg = glm(Y~X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.45885    1.04646  -4.261 2.04e-05 ***
X1           0.51664    0.11178   4.622 3.80e-06 ***
X2           0.21008    0.07247   2.899 0.003745 **
X3B          1.74496    0.49952   3.493 0.000477 ***
X3C         -0.03470    0.35691  -0.097 0.922543
X3D          0.08004    0.34916   0.229 0.818672
X3E          2.21966    0.56475   3.930 8.48e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 552.64  on 999  degrees of freedom
Residual deviance: 397.69  on 993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

Here, the reference is modality $A$. Which means that for someone with characteristics $(X_1,X_2,X_3=A)$, we predict the following probability

$p=H(\widehat\beta_0+\widehat\beta_1 X_1+\widehat\beta_2 X_2)$

where $H(\cdot)$ denotes the cumulative distribution function of the logistic distribution

$H(x)=\frac{e^x}{1+e^x}$

For someone with characteristics $(X_1,X_2,X_3=B)$, we predict the following probability

$p=H(\widehat\beta_0+\widehat\beta_1 X_1+\widehat\beta_2 X_2+\widehat\beta_3^{\ (B)})$

For someone with characteristics $(X_1,X_2,X_3=C)$, we predict the following probability

$p=H(\widehat\beta_0+\widehat\beta_1 X_1+\widehat\beta_2 X_2+\widehat\beta_3^{\ (C)})$

(etc.) Here, if we accept $H_0:\beta_3^{\ (C)}=0$ (against $H_1:\beta_3^{\ (C)}\neq0$), it means that modality $C$ cannot be considerd as different from $A$.

A natural idea can be to change the reference modality, and to look at the $p$-values. If we consider the following loop, we get

> M = matrix(NA,5,5)
> rownames(M)=colnames(M)=LETTERS[1:5]
> for(k in 1:5){
+ db$X3 = relevel(X3,LETTERS[k]) + reg = glm(Y~X1+X2+X3,family=binomial,data=db) + M[levels(db$X3)[-1],k] = summary(reg)$coefficients[4:7,4] + } > M A B C D E A NA 0.0004771853 9.225428e-01 0.8186723647 8.482647e-05 B 4.771853e-04 NA 4.841204e-04 0.0009474491 4.743636e-01 C 9.225428e-01 0.0004841204 NA 0.7506242347 9.194193e-05 D 8.186724e-01 0.0009474491 7.506242e-01 NA 1.730589e-04 E 8.482647e-05 0.4743636442 9.194193e-05 0.0001730589 NA and if we simply want to know if the $p$-value exceeds – or not – 5%, we get the following, > M.TF = M>.05 > M.TF A B C D E A NA FALSE TRUE TRUE FALSE B FALSE NA FALSE FALSE TRUE C TRUE FALSE NA TRUE FALSE D TRUE FALSE TRUE NA FALSE E FALSE TRUE FALSE FALSE NA The first column is obtained when $A$ is the reference, and then, we see which parameter should be considered as null. The interpretation is the following: • $C$ and $D$ are not different from $A$ • $E$ is not different from $B$ • $A$ and $D$ are not different from $C$ • $A$ and $C$ are not different from $D$ • $B$ is not different from $E$ Note that we only have, here, some kind of intuition. So, let us run a more formal test. Let us consider the following regression (we remove the intercept to get a model easier to understand) > library(car) > db$X3=relevel(X3,"A")
> reg=glm(Y~0+X1+X2+X3,family=binomial,data=db)
> summary(reg)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
X1   0.51664    0.11178   4.622 3.80e-06 ***
X2   0.21008    0.07247   2.899  0.00374 **
X3A -4.45885    1.04646  -4.261 2.04e-05 ***
X3E -2.23919    1.06666  -2.099  0.03580 *
X3D -4.37881    1.04887  -4.175 2.98e-05 ***
X3C -4.49355    1.06266  -4.229 2.35e-05 ***
X3B -2.71389    1.07274  -2.530  0.01141 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1386.29  on 1000  degrees of freedom
Residual deviance:  397.69  on  993  degrees of freedom
AIC: 411.69

Number of Fisher Scoring iterations: 7

It is possible to use Fisher test to test if some coefficients are equal, or not (more generally if some linear constraints are satisfied)

> linearHypothesis(reg,c("X3A=X3C","X3A=X3D","X3B=X3E"))
Linear hypothesis test

Hypothesis:
X3A - X3C = 0
X3A - X3D = 0
- X3E  + X3B = 0

Model 1: restricted model
Model 2: Y ~ 0 + X1 + X2 + X3

Res.Df Df  Chisq Pr(>Chisq)
1    996
2    993  3 0.6191      0.892

Here, we clearly accept the assumption that the first three factors are equal, as well as the last two. What is the next step? Well, if we believe that there are mainly two categories, $\{A,C,D\}$ and $\{B,E\}$, let us create that factor,

> X3bis=rep(NA,length(X3))
> X3bis[X3%in%c("A","C","D")]="ACD"
> X3bis[X3%in%c("B","E")]="BE"
> db$X3bis=as.factor(X3bis) > reg=glm(Y~X1+X2+X3bis,family=binomial,data=db) > summary(reg) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.39439 1.02791 -4.275 1.91e-05 *** X1 0.51378 0.11138 4.613 3.97e-06 *** X2 0.20807 0.07234 2.876 0.00402 ** X3bisBE 1.94905 0.36852 5.289 1.23e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 552.64 on 999 degrees of freedom Residual deviance: 398.31 on 996 degrees of freedom AIC: 406.31 Number of Fisher Scoring iterations: 7 Here, all the categories are significant. So we do have a proper model. # Une région géographique n’est pas une variable continue En relisant les devoirs maisons, je me suis rendu compte que certains avaient tenté de regrouper les régions (géographiques) par régions homogènes. Sauf que les régions étaient codées par un numéro (selon la codification officielle). Par exemple, dans une des bases, nous avions des assurés dans 4 zones géographiques, à savoir la région 82 (région Rhône-Alpes en rouge) la région 54 (région Poitou-Charentes en vert) la région 73 (région Midi-Pyrénées en bleu) et enfin la région 41 (région Lorraine en mauve). > unique(baseFREQ$region)
[1] 82 54 73 41

Une idée intéressante pour regrouper les régions pouvait être d’utiliser les arbres. Les régions étant des couleurs (on le voit bien sur la carte) et pas des variables quantitatives, il est normal de travailler sur des facteurs. D’ailleurs le code pour faire la carte est le suivant,

> library(maps)
>  france<-map(database="france")
>  dpt=c("Ain","Ardeche","Drome","Isere","Loire ","Rhone",
+ "Savoie","Haute-Savoie","Charente","Charente-Maritime",
+ "Deux-Sevres","Vienne","Ariege","Aveyron","Haute-Garonne",
+ "Gers","Lot","Hautes-Pyrenees","Tarn","Tarn-et-Garonne",
+ "Meurthe-et-Moselle","Meuse","Moselle","Vosges")
>  couleur=c(rep(2,8),rep(3,4),rep(4,8),rep(6,4))
>  match=match.map(france,dpt)
>  color=couleur[match]
>  map(database="france", fill=TRUE, col=color)

L’arbre sur les régions en tant que facteurs donne le découpage suivant

>  baseFREQ$fregion=as.factor(baseFREQ$region)
>  ARBRE1=tree(nombre~fregion,data=baseFREQ,split="gini")
>  plot(ARBRE1)
>  text(ARBRE1)

Bon, R a la mauvaise idée de recoder les classes (mais il garde l’ordre, i.e. a correspond à la région 41, b à 54, c à 73 et d à 82). Visuellement, on retient qu’il est possible de considérer deux grandes régions, AC (i.e. 41 et 73) et BD (i.e. 54 et 72). L’intérêt des arbres sur des variables qualitatives, des facteurs, c’est que tous les regroupements sont possibles. En revanche, si on fait un arbre sur la région qui est lue en tant que nombre (quantitatif), on obtient

>  ARBRE2=tree(nombre~region,data=baseFREQ,split="gini")
>  plot(ARBRE2)
>  text(ARBRE2)

Il est alors impossible de regrouper dans une même classe deux régions séparées par un nombre, i.e. on ne peut regrouper 41 et 82 dans la même classe. R suggère de distinguer peut être trois régions, à savoir 82 (à droite), puis 73 (au centre) et enfin de mettre éventuellement 41 et 54 ensemble. Ce qui n’est pas la stratégie optimale quand on regroupe des facteurs.

# Régression sur des variables catégorielles

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,
> 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,
> 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),
> 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),
> 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),
> 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),
> 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),
> 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),
> 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

> library(tree)
> arbre=tree(nbre~zone+offset(log(exposition)),
+ data=base,split="gini")
> plot(arbre)
> text(arbre)