All posts by Arthur Charpentier

Arthur Charpentier, professor in Montréal, in Actuarial Science. Former professor-assistant at ENSAE Paristech, associate professor at Ecole Polytechnique and assistant professor in Economics at Université de Rennes 1.  Graduated from ENSAE, Master in Mathematical Economics (Paris Dauphine), PhD in Mathematics (KU Leuven), and Fellow of the French Institute of Actuaries.

On the “correlation” between a continuous and a categorical variable

Let us get back on the Titanic dataset,

loc_fichier = "http://freakonometrics.free.fr/titanic.RData"
download.file(loc_fichier, "titanic.RData")
load("titanic.RData")
base = base[!is.na(base$Age),]

On consider two variables, the age x (the continuous one) and the survivor indicator y (the qualitative one)

X = base$Age
Y = base$Survived

It looks like the age might be a valid explanatory variable in the logistic regression,

summary(glm(Survived~Age,data=base,family=binomial))
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.05672    0.17358  -0.327   0.7438  
Age         -0.01096    0.00533  -2.057   0.0397 *
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 964.52  on 713  degrees of freedom
Residual deviance: 960.23  on 712  degrees of freedom
AIC: 964.23

The significance test here has a p-value just below 4%. Actually, one can relate it with the value of the deviance (the null deviance and the residual deviance). Recall thatD=2\big(\log\mathcal{L}(\boldsymbol{y})-\log\mathcal{L}(\widehat{\boldsymbol{\mu}})\big)whileD_0=2\big(\log\mathcal{L}(\boldsymbol{y})-\log\mathcal{L}(\overline{y})\big)Under the assumption that x is worthless, D_0-D tends to a \chi^2 distribution with 1 degree of freedom. And we can compute the p-value dof that likelihood ratio test,

1-pchisq(964.52-960.23,1)
[1] 0.03833717

(which is consistent with a Gaussian test). But if we consider a nonlinear transformation

summary(glm(Survived~bs(Age),data=base,family=binomial))
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.8648     0.3460   2.500 0.012433 *  
bs(Age)1     -3.6772     1.0458  -3.516 0.000438 ***
bs(Age)2      1.7430     1.1068   1.575 0.115299    
bs(Age)3     -3.9251     1.4544  -2.699 0.006961 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 964.52  on 713  degrees of freedom
Residual deviance: 948.69  on 710  degrees of freedom

which seems to be “more significant”

1-pchisq(964.52-948.69,3)
[1] 0.001228712

So it looks like the variable x is interesting here.

To visualize the non-null correlation, one can consider the condition distribution of x given y=1, and compare it with the condition distribution of x given y=0,

ks.test(X[Y==0],X[Y==1])
 
	Two-sample Kolmogorov-Smirnov test
 
data:  X[Y == 0] and X[Y == 1]
D = 0.088777, p-value = 0.1324
alternative hypothesis: two-sided

i.e. with a p-value above 10%, the two distributions are not significatly different.

F0 = function(x) mean(X[Y==0]<=x)
F1 = function(x) mean(X[Y==1]<=x)
vx = seq(0,80,by=.1)
vy0 = Vectorize(F0)(vx)
vy1 = Vectorize(F1)(vx)
plot(vx,vy0,col="red",type="s")
lines(vx,vy1,col="blue",type="s")

An alternative is discretize variable x and to use Pearson’s independence test,

k=5
LV = quantile(X,(0:k)/k)
LV[1] = 0
Xc = cut(X,LV)
table(Xc,Y)
           Y
Xc           0  1
  (0,19]    85 79
  (19,25]   92 45
  (25,31.8] 77 50
  (31.8,41] 81 63
  (41,80]   89 53
chisq.test(table(Xc,Y))
 
	Pearson's Chi-squared test
 
data:  table(Xc, Y)
X-squared = 8.6155, df = 4, p-value = 0.07146

The p-value is here 7%, with five categories for the age. And actually, we can compare the p-value

pvalue = function(k=5){
LV = quantile(X,(0:k)/k)
LV[1] = 0
Xc = cut(X,LV)
chisq.test(table(Xc,Y))$p.value}
vk = 2:20
vp = Vectorize(pvalue)(vk)
plot(vk,vp,type="l")
abline(h=.05,col="red",lty=2)

which gives a p-value close to 5%, as soon as we have enough categories. In the slides of the course (STT5100), I claim that actually, the age is an important variable when trying to predict if a passenger survived. Test mentioned here are not as conclusive, nevertheless…

STT5100, quiz (régression logistique)

Je mets en ligne des éléments de correction pour le quiz donné lundi dernier, en lien avec la régression logistique.

  • Quiz 1 – régression logistique (sur données individuelles)

On disposait de la base suivante, avec 44 bières testées, et un jugement porté (bon, ou pas),

download.file("http://freakonometrics.free.fr/Beer.RData","Beer.RData")
load("Beer.RData")
str(base)
'data.frame':	44 obs. of  7 variables:
 $ Brewery    : Factor w/ 8 levels "Bauhaus","Bent Paddle",..: 1 1 1 1 2 2 2 2 2 4 ...
 $ Beer       : Factor w/ 44 levels "14* ESB","B-Side Pils",..: 44 38 43 36 19 42 5 1 7 11 ...
 $ Description: Factor w/ 37 levels "American Blonde Ale",..: 22 15 37 20 32 28 16 13 3 26 ...
 $ Style      : Factor w/ 3 levels "Ale","IPA","Lager": 3 3 3 2 1 3 2 1 1 1 ...
 $ ABV        : num  5.4 5 5.4 6.7 4.8 5 6.2 5.6 6 5.4 ...
 $ IBU        : int  48 28 55 70 48 38 68 32 34 45 ...
 $ Good       : int  0 0 0 0 0 0 0 0 0 1 ...

Le but est d’estimer la probabilité de trouver une bière bonne, à partir de trois caractéristiques: le degré d’alcool (ABV-Alcool By Volume), la mesure d’amertume (IBU-International Bitterness Unit) et le style (trois styles sont ici présent, les lager, les ale, et les IPA – qui techniquement sont des ale, mais peu importe),

newbeer = data.frame(ABV=6,IBU=78,Style="IPA")

Pour répondre, on se lance

reg = glm(Good~ABV+IBU+Style,data=base,family=binomial)
summary(reg)
 
Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -10.98183    4.71577  -2.329   0.0199 *
ABV            1.42984    0.79147   1.807   0.0708 .
IBU            0.04905    0.04005   1.225   0.2207  
StyleIPA      -2.42680    1.67669  -1.447   0.1478  
StyleLager   -16.69972 2009.50187  -0.008   0.9934

Si on calcule la probabilité à partir de ce modèle, on obtient

predict(reg,newdata=newbeer,type="response")
        1 
0.2682531

On peut tenter de simplifier un peu le modèle, en notant que le style semblant peu significatif,

reg = glm(Good~ABV+IBU,data=base,family=binomial)
summary(reg)
 
Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept) -9.307699   3.678245  -2.530   0.0114 *
ABV          1.308989   0.721717   1.814   0.0697 .
IBU          0.008248   0.025139   0.328   0.7428

et comme IBU n’est toujours pas significatif, on utilise

reg = glm(Good~ABV+IBU,data=base,family=binomial)
summary(reg)
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -9.7881     3.3959  -2.882  0.00395 **
ABV           1.4662     0.5481   2.675  0.00747 **

Cette fois, la variable restante est significative. On obtient la prévision suivante

predict(reg,newdata=newbeer,type="response")
        1 
0.2707103

très proche de l’ancienne, de l’ordre de 27%. On peut juste regarder s’il n’y aurait pas d’effet non-linéaire

library(gam)
reg = gam(Good~bs(ABV),data=base,family=binomial)
plot(reg,se=TRUE)

et non, on peut garder l’idée d’un modèle linéaire.

Par contre, on est un peu frustrés parce que la prévision qu’on ferait correspond à la fréquence de bières jugées bonnes dans la population

mean(base$Good)
[1] 0.2727273

que l’on aurait obtenue avec

reg = glm(Good~1,data=base,family=binomial)
summary(reg)
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -0.9808     0.3385  -2.898  0.00376 **
 
predict(reg,newdata=newbeer,type="response")
        1 
0.2727273

Si le taux d’alcool passe de 6% à 9%, cette probabilité passe de 27% à 97%

reg = glm(Good~ABV,data=base,family=binomial)
newbeer = data.frame(ABV=9,IBU=78,Style="IPA")
predict(reg,newdata=newbeer,type="response")
        1 
0.9679383

Si on reprend le modèle avec trop de variables (non-significatives), on obtiendrait

reg = glm(Good~ABV+IBU+Style,data=base,family=binomial)
predict(reg,newdata=newbeer,type="response")
        1 
0.9639461

dont on peut également demander la précision, en calculant l’écart-type de cette prévision

predict(reg,newdata=newbeer,type="response",se.fit=TRUE)
$fit
        1 
0.9639461 
 
$se.fit
        1 
0.0680774
  • Quiz 2 : régression logistique sur données agrégées

La seconde question était un peu plus technique, avec la base suivante, contenant des données agrégées

download.file("http://freakonometrics.free.fr/BerkeleyAdmissions.RData","BerkeleyAdmissions.RData")
load("BerkeleyAdmissions.RData")
str(database)
'data.frame': 24 obs. of  4 variables:
$ Admit : Factor w/ 2 levels "Admitted","Rejected": 1 2 1 2 1 2 1 2 1 2 ...
$ Gender: Factor w/ 2 levels "Male","Female": 1 1 2 2 1 1 2 2 1 1 ...
$ Dept  : Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 2 2 2 2 3 3 ...
$ Freq  : num  512 313 89 19 353 207 17 8 120 205 ...

Par exemple, pour une personne de sexe féminin, ayant postulé dans le département A, il y a eu 89 admissions, sur 108 dossiers reçus, soit 82.4% de chances d’être admise. La question était d’estimer la probabilité d’admission avec une régression logistique, sur les facteurs Gender et Dept (sans effet croisé) était utilisée ?

Pour commencer, j’avoue que ce format de données ne me plait pas trop : par genre et par département, on va rajouter deux variables : le nombre d’acceptation et le nombre de rejets

BA = database[database$Admit=="Admitted",]
BA=BA[,2:4]
names(BA)=c("Gender","Dept","NAdm")
BR = database[database$Admit=="Rejected",]
BR=BR[,2:4]
names(BR)=c("Gender","Dept","NRej")
B=merge(BA,BR)
B$taux = B$NAdm/(B$NAdm+B$NRej)
B
   Gender Dept NAdm NRej       taux
1  Female    A   89   19 0.82407407
2  Female    B   17    8 0.68000000
3  Female    C  202  391 0.34064081
4  Female    D  131  244 0.34933333
5  Female    E   94  299 0.23918575
6  Female    F   24  317 0.07038123
7    Male    A  512  313 0.62060606
8    Male    B  353  207 0.63035714
9    Male    C  120  205 0.36923077
10   Male    D  138  279 0.33093525
11   Male    E   53  138 0.27748691
12   Male    F   22  351 0.05898123

On retrouve, sur la première ligne, le 82.4% mentionné. Pour construite le modèle, on utilise la syntaxe suivante

reggroup= glm(cbind(NAdm,NRej)~Gender+Dept,data=B,family=binomial)
summary(reggroup)
 
Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.58205    0.06899   8.436   <2e-16 ***
GenderFemale  0.09987    0.08085   1.235    0.217    
DeptB        -0.04340    0.10984  -0.395    0.693    
DeptC        -1.26260    0.10663 -11.841   <2e-16 ***
DeptD        -1.29461    0.10582 -12.234   <2e-16 ***
DeptE        -1.73931    0.12611 -13.792   <2e-16 ***
DeptF        -3.30648    0.16998 -19.452   <2e-16 ***

On peut alors faire notre prévision de la probabilité d’être acceptée, pour une femme postulant dans le département A. On peut rajouter les probabilités prédite dans notre base,

B$predgroup=predict(reggroup,type="response")
> B
   Gender Dept NAdm NRej       taux  predgroup
1  Female    A   89   19 0.82407407 0.66416742
2  Female    B   17    8 0.68000000 0.65441963
3  Female    C  202  391 0.34064081 0.35877694
4  Female    D  131  244 0.34933333 0.35144696
5  Female    E   94  299 0.23918575 0.25780964
6  Female    F   24  317 0.07038123 0.06757450
7    Male    A  512  313 0.62060606 0.64153930
8    Male    B  353  207 0.63035714 0.63149912
9    Male    C  120  205 0.36923077 0.33613931
10   Male    D  138  279 0.33093525 0.32903451
11   Male    E   53  138 0.27748691 0.23916654
12   Male    F   22  351 0.05898123 0.06154717

Et on obtient une probabilité d’être acceptatée de 66%. Mais notons qu’on peut procéder autrement (on reviendra sur cette idée plus tard dans le cours): on peut recréer la base originale, qui a permis d’obtenir le tableau de contingence.

base = NULL
for(i in 1:nrow(B)){
     L = B[i,1:4]
     vL = data.frame(Gender=L[1],Dept=L[2],Y=c(rep(1,L[3]),rep(0,L[4])))
     base = rbind(base,vL)}
str(base)
'data.frame':	4526 obs. of  3 variables:
 $ Gender: Factor w/ 2 levels "Male","Female": 2 2 2 2 2 2 2 2 2 2 ...
 $ Dept  : Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Y     : num  1 1 1 1 1 1 1 1 1 1 ...

On peut alors faire la régression sur cette base, avec les 4526 individus, et rajouter les prévisions, dans la base

regindiv = glm(Y~Gender+Dept,data=base,family=binomial)
B$predindiv=predict(regindiv,type="response",newdata=B)
B
   Gender Dept NAdm NRej       taux  predgroup  predindiv
1  Female    A   89   19 0.82407407 0.66416742 0.66416742
2  Female    B   17    8 0.68000000 0.65441963 0.65441963
3  Female    C  202  391 0.34064081 0.35877694 0.35877694
4  Female    D  131  244 0.34933333 0.35144696 0.35144696
5  Female    E   94  299 0.23918575 0.25780964 0.25780964
6  Female    F   24  317 0.07038123 0.06757450 0.06757450
7    Male    A  512  313 0.62060606 0.64153930 0.64153930
8    Male    B  353  207 0.63035714 0.63149912 0.63149912
9    Male    C  120  205 0.36923077 0.33613931 0.33613931
10   Male    D  138  279 0.33093525 0.32903451 0.32903451
11   Male    E   53  138 0.27748691 0.23916654 0.23916654
12   Male    F   22  351 0.05898123 0.06154717 0.06154717

Les prévisions sont identiques, et les modèles aussi,

summary(regindiv)
 
Call:
glm(formula = Y ~ Gender + Dept, family = binomial, data = base)
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4773  -0.9306  -0.3741   0.9588   2.3613  
 
Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.58205    0.06899   8.436   <2e-16 ***
GenderFemale  0.09987    0.08085   1.235    0.217    
DeptB        -0.04340    0.10984  -0.395    0.693    
DeptC        -1.26260    0.10663 -11.841   <2e-16 ***
DeptD        -1.29461    0.10582 -12.234   <2e-16 ***
DeptE        -1.73931    0.12611 -13.792   <2e-16 ***
DeptF        -3.30648    0.16998 -19.452   <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: 6044.3  on 4525  degrees of freedom
Residual deviance: 5187.5  on 4519  degrees of freedom
AIC: 5201.5
 
Number of Fisher Scoring iterations: 5

Bref, la probabilité prédite par le modèle linéaire, sans effet croisé, serait de l’ordre de 66%. Petit moment culturel: cette dernière base est tirée d’un article relativement connu, Bickel, Hammel & O’Connell (1975), en lien avec le paradoxe de Simpson. Et je recommande largement la lecture !