Tag Archives: factor

Regression on factors

Most of our intuitions about regression models come from the Gaussian standard linear model. One interesting feature is that, when we have a factor explanatory variable, the sum of predictions per class is the sum of observations of the endogeneous variable, per class. To be more specific, consider some factor variable https://latex.codecogs.com/gif.latex?x_1\in\{0,1\}, and a regression model

https://latex.codecogs.com/gif.latex?y_i=\beta_0+\beta_1%20\boldsymbol{1}(x_1=1)+\beta_2%20x_2+\varepsilon_i

Use ordinary least squares to fit that model

https://latex.codecogs.com/gif.latex?\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1%20\boldsymbol{1}(x_1=1)+\widehat{\beta}_2%20x_2

Then for all https://latex.codecogs.com/gif.latex?x\in\{0,1\}

https://latex.codecogs.com/gif.latex?\sum_{i:x_i=x}%20y_i%20=%20\sum_{i:x_i=x}%20\widehat{y}_i

> n=200
> X1=rep(0:1,each=n/2)
> set.seed(1)
> X2=runif(2*n)
> L=X1-X2
> B=data.frame(Y=rnorm(n,L),X1=as.factor(X1),X2=X2)
> pd=aggregate(x=B$Y,by=list(B$X1),mean)$x
> pd
[1] -0.4881735  0.5341301
> fit=lm(Y~X1+X2,data=B)
> B2=data.frame(x=B$X1,y=predict(fit))
> aggregate(x=B2$y,by=list(B2$x),mean)$x
[1] -0.4881735  0.5341301

Continue reading Regression on factors

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 . Which means that for someone with characteristics , we predict the following probability

where  denotes the cumulative distribution function of the logistic distribution

For someone with characteristics , we predict the following probability

For someone with characteristics , we predict the following probability

(etc.) Here, if we accept  (against ), it means that modality  cannot be considerd as different from .

A natural idea can be to change the reference modality, and to look at the -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 -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  is the reference, and then, we see which parameter should be considered as null. The interpretation is the following:

  •  and  are not different from 
  •  is not different from 
  •  and  are not different from 
  •  and  are not different from 
  •  is not different from 

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,  and , 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, 
+ 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

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