Tag Archives: Simpson

INF7100, les données

23/05/2020 Arthur Charpentier Leave a comment

La première partie de mon intervention sur la science des données, dans le cadre du cours INF7100 portera sur les données (et la distinction entre les données d’observations et les données d’expérience). Le plan est le suivant

111: Observation vs. Expérimentation pdf video (19:48)
112: Données Observationnelles et Biais pdf video (22:39)
121: Paradoxe de Simpson pdf video (22:24)
122: Recherche de Contrefactuels pdf video (25:27)
123: A/B Testing et Renforcement pdf video (12:13)
131: Incertitude (1) pdf video (39:01)
132: Incertitude (2) pdf video (41:45)

(oui je suis passé sur youtube pour ce cours… mon compte viméo a malheureusement été suspendu… et toutes les capsules de mon cours précédent ne sont plus visibles)

STT5100

STT5100, quiz (régression logistique)

03/04/2020 Arthur Charpentier Leave a comment

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(&gt;|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(&gt;|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(&gt;|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(&gt;|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(&gt;|z|)    
(Intercept)   0.58205    0.06899   8.436   &lt;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   &lt;2e-16 ***
DeptD        -1.29461    0.10582 -12.234   &lt;2e-16 ***
DeptE        -1.73931    0.12611 -13.792   &lt;2e-16 ***
DeptF        -3.30648    0.16998 -19.452   &lt;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")
&gt; 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(&gt;|z|)    
(Intercept)   0.58205    0.06899   8.436   &lt;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   &lt;2e-16 ***
DeptD        -1.29461    0.10582 -12.234   &lt;2e-16 ***
DeptE        -1.73931    0.12611 -13.792   &lt;2e-16 ***
DeptF        -3.30648    0.16998 -19.452   &lt;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 !

Decision Theory, Economics, Maths, Microeconomics, Probability, Risk

From Simpson’s paradox to pies

15/02/2013 Arthur Charpentier 3 Comments

Today, I wanted to publish a post on economics, and decision theory. And probability too… Those who do follow my blog should know that I am a big fan of Simpson’s paradox. I also love to mention it in my
econometric classes. It does raise important questions, that I do relate to multicolinearity, and interepretations of regression models, with multiple (negatively correlated) explanatory variables. This paradox has amazing pedogological virtues. I did mention it several times on this blog (I should probably mention that I discovered this paradox via Marco Scarsini, who did learn me a lot of things, in decision theory and in probability). For those who do not know this paradox, here is an example that Marco gave in one of his talk, a few years ago. Consider the following statistics, when healthy people entered in some hospital

hospital	total	survivors	deaths	survival rate
hospital A	600	590	10	98%
hospital B	900	870	30	97%

while, when sick people entered in the same hospitals

hospital	total	survivors	deaths	survival rate
hospital A	400	210	190	53%
hospital B	100	30	70	30%

Somehow, whatever your health situation, you should choose hospital A. Now, if we agregate

hospital	total	survivors	deaths	survival rate
hospital A	1000	800	200	80%
hospital B	1000	900	100	90%

i.e. without any doubts, people should choose hospital B.

Actually, Simpson’s paradox is called Simpson’s paradox because Colin Blyth named it that way in 1972, in his paper entitled on Simpson’s paradox and the sure-thing principle (an economic article in a statistical journal), that can be downloaded from http://www.stat.cmu.edu/~fienberg/…. He found this paradox in a paper published in 1951 by Edward Simpson, even if other papers actually did mention it earlier. The most popular application is probably admission at Berckley’s graduate studiesprograms, and sex bias, see Bickel, Hammel & O’Connell (1975), that can be downloaded from http://www.unc.edu/~nielsen/…. I also mentioned a geometric interpretation of this paradox a few years ago on my blog, which is so simple to understand that the paradox is no longer a paradox actually, since on the example above, we had

and

while

With symbolic notations, one can have at the same time

and

with also

as shown on the graph below

There should be connection between Simpson’s paradox and the ecological fallacy (which is an issue I recently discovered and that I found extremely interesting, related again to difficulties of interpreting
regressions). But that’s another story. My point today is that Colin Blyth did mention another nice paradox, that is related, this time, to stochastic orderings. The idea is the following. Consider the three spinners drawn below (imagine some arrows in those circles)

spinner A: no matter where the arrow stops, the gain is 3,
spinner B: 56% chances to gain 2, 22% chances to gain 4, and 22% chances to gain 6,
spinner C: 51% chances to gain 1, 49% chances to gain 5.

Instead of spinners, it is also possible to consider three different lotteries,

You play against a friend, you pick a spinner, while the friend picks another. Everyone flick his arrow, the highest number wins (no matter the difference). Let us compute the odds. First case, A against B, from
A’s perspective

	B-2	B-4	B-6
A-3	56% +1 win	22% -1 lose	22% -3 lose

In that case, A has 56% chance of beating B. Second case, A against C, from A’s perspective,

	C-1	C-5
A-3	51% +1 win	49% -2 lose

In that case, A has 51% chance of beating C. Third (an final) case, B against C, from B’s perspective. Assuming independence between the spinners, joint probabilities can easily be computed,

	C-1	C-5
B-2	28.56% +1 win	27.44% -3 lose
B-4	11.22% +3 win	10.78% -1 lose
B-6	11.22% +5 win	10.78% +1 win

In that case, B has 61.78% chance of beating C. So, if we try to summarize,

A is the best choice, since it beats both with – always – more than 50% chance,
C is the worst choice, since it is beaten by both with – always – more than 50% chance,

Now, assume that you play not against one friend, but two friends. An everyone picks a different spinner. Let
us compute the odds, one more time. First case, A against B and C, from A’s perspective

	B-2 C-1	B-2 C-5	B-4 C-1	B-4 C-5	B-6 C-1	B-6 C-5
A-3	28.56% +1 win	27.44% -2 lose	11.22% -1 lose	10.78% -1 lose	11.22% -3 lose	10.78% -3 lose

In that case, A has 28.56% chance of beating B and C. Second case, B against A and C, from B’s perspective,

	A-3 C-1	A-3 C-5
B-2	28.56% -1 lose	27.44% -2 lose
B-4	11.22% +1 win	10.78% -1 lose
B-6	11.22% +3 win	10.78% +1 win

In that case, B has 33.22% chance of beating A and B.Third (an final) case, C against A, from C’s perspective,

	A-3 B-2	A-3 B-4	A-3 B-6
C-1	28.56% -2 lose	11.22% -3 lose	11.22% -5 lose
C-5	27.44% +2 win	10.78% +1 win	10.78% -1 lose

In that case, C has 38.22% chance of beating A and B. So, if we try to summarize, this time

C is the best choice, since has (strictly) more than 1/3 chances to win, which the highest probability
A is the worst choice, since has (strictly) less than 1/3 chances to win, which the lowest probability

Odd isn’t it ? Now, is there an interpretation of that paradox ? Yes, Martin Gardner, in his paper on induction and probability, mentioned the case of drug testing. The value we had with the spinner is the health level, rated from 1 to 6. Thus, taking drug A, you always get an average health level of 3. With drug C, on the other hand, you get either very sick (level 1) or very well (level 5). Consider now a doctor who wants to maximize the patient’s chance of being well. If only pills A and C are available, then the doctor should choose A. This is what we’ve seen in the first part. Assume that now a company delivers a third pill, called drug B. Then the doctor should find C more interesting…. Odd, isn’t it ?

Colin Blyth gave a more amusing application. Assume that you like to go to the restaurant, and you like get a dessert there. Dessert A – the apple pie – is the average one, with a standard level, that you rank 3 (on a scale from 1 to 6). Dessert C – the cheese cake – can either be awfull (ranked 1) or delicious (ranked 5). You’d better go for the apple pie if you want to maximize the probability of not being disappointed (i.e. maximizing your “best chance” according to Colin Blyth, but I guess it can be interpreted as regret minimization too). Now assume that dessert B – the blueberry pie – is available (with ranks given by the spinner). Then you should go for the cheese cake. I let you imagine the discussion that you can have, then, with your favorite waitress

– Hi Mr Freakonometrics, do you want a piece of apple pie ? (yes, actually she also comes frequently on my blog, and knows me from my pseudo…)

– Probably. But actually, I was wondering if you did have your blueberry pie today ?

– Yes, in fact we do….

– Great, in that case, I’ll go for the cheese cake.

She’ll probably think that I am freak… so I hope she’ll come and read my post, to understand that, actually, it does make a lot of sense to go for what was supposed to be my worst case.

Maths

Paradoxe de Simpson expliqué par un dessin

10/06/2010 Arthur Charpentier

Il y a plusieurs mois, j’avais fait un billet sur le paradoxe de Simpson (ici). Pour reprendre un exemple simple, supposons que l’on a deux hôpitaux, affichant les statistiques suivantes, pour les personnes “saines”,

hôpital	total	survivants	décès	taux de survie
hopital A	600	590	10	98%	x
hopital B	900	870	30	97%

et pour les personnes “malades” (on suppose vraiment que le critère est identique pour les deux hôpitaux)

hôpital	total	survivants	décès	taux de survie
hopital A	400	210	190	53%	x
hopital B	100	30	70	30%

Bref, peu importe son état de santé, on a toujours intérêt à choisir l’hôpital A. D’où le paradoxe. En revanche, si on combine, on obtient “globalement”,

hôpital	total	survivants	décès	taux de survie
hopital A	1000	800	200	80%
hopital B	1000	900	100	90%	x

on note que l’on a intérêt à choisir l’hôpital B. Étonnant non ? (l’exemple vient de Marco Scarsini, et je n’ai jamais trouvé plus élégant).

J’évoque souvent ce paradoxe en économétrie, et c’est ce que j’avais fait dans mon précédant billet (ici). En fait, Roger Nelsen propose un petit exemple simple permettant de comprendre comment on peut obtenir un tel paradoxe. L’idée est de noter que

et pourtant

Cela n’a a priori rien d’exceptionnel… On peut parfaitement avoir

tout en ayant,

C’est ce que l’on peut visualiser ci-dessous,

Merci Roger….

"sendo l'intento mio scrivere cosa utile a chi la intende…"