Tag Archives: Fisher

Une courte histoire des expériences randomisées

On a beaucoup parlé ces derniers mois des expériences randomisées. Cet automne pour commencer, lorsque le prix Nobel d’économie a été attribué à Esther Duflo, Michael Kremer et Abhijit Banerjee. Plus récemment, lors de la polémique, en France, où, dans une conférence donnée le 13 février 2020, intitulée contre la méthode, un médecin marseillais affirmait « moi je n’ai jamais fait d’essais randomisés […] faire ça sur des maladies infectieuses, ça n’a pas de sens ».

Les expériences randomisées

Quand on se demande l’effet d’un traitement (par rapport à un autre), la bonne méthode consiste à se poser la question « que ce serait-il passé si la personne avait choisi l’autre option ? » – ce que John Stuart Mill parlait de la « méthode des différences ». On parle ici d’inférence causale car on ne peut pas savoir ce qu’aurait donné l’autre option, qui n’a pas été retenue. L’effet causal d’une option, ou d’un traitement, n’est jamais observé car il s’agit de la différence potentielle entre deux alternatives. La méthode la plus simple consiste à prendre deux patients proches (à défaut d’être identiques), et de voir comment les deux répondent aux deux options, la proximité se mesurant en fonction de covariables (comme l’âge, le genre, etc), et les patients n’ayant pas choisi les options proposées. Il y a trois points importants dans cette approche : (1) la notion de groupe de contrôle, autrement dit il faut qu’il y ait des personnes « non traitées » (2) la notion d’aveuglement, ou d’absence d’information, autrement dit, les patients ne doivent pas savoir quel traitement ils ont eu (3) la notion d’affectation aléatoire, autrement dit, le traitement doit être choisi de manière totalement indépendante des patients. Ce sont des points qui opposent fondamentalement les expériences randomisées aux données observationnelles, comme les oppose par exemple Paul Rosenbaum dans Observation and Experiment. Si la méthode d’expérience randomisée a été formalisée par Ronald Fisher en 1934, nous allons revenir ici sur deux moments importants dans l’histoire de cette méthodologie. La première est l’utilisation de la randomisation pour faire des tests autour de 1890, en particulier avec les travaux de Peirce puis Richet, sur le spiritualisme et la transmission de pensée (qui ont très probablement influencé Fisher). Et la seconde est l’expérience qui a été menée en 1954 aux États-Unis, pour tester l’efficacité d’un vaccin contre la polio (et qui montre au passage que faire des expériences randomisées dans un contexte de maladies infectieuses a tout son sens).

Continue reading Une courte histoire des expériences randomisées

De la pratique de la régression

Depuis le début de la session, j’ai imposé une petite innovation, en donnant, environ une semaine sur deux, un petit exercice (obligatoire mais non noté) avant le cours, en vue de forcer à réfléchir (et de donner des éléments de réponse). Par exemple pour le premier cours, il fallait “prévoir” une valeur manquante, et le but était de montrer que, naturellement, on choisit la valeur moyenne.

Pour demain, j’avais posé un exercice un peu plus compliqué, sachant qu’on avait vu, lors du dernier cours, comme faire une régression linéaire, et qu’on avait fini en discutant les tests simples (en lien avec la significativité) et les tests multiples. Pour l’exercice, j’avais mis en ligne une petite base de données,

download.file("http://freakonometrics.free.fr/data3.RData","data3.RData")
load("data3.RData")
str(df)
'data.frame':	147 obs. of  10 variables:
 $ Y : num  11.72 15.91 14.19 11.15 8.31 ...
 $ X1: num  1.33 3.18 0.28 2.08 0.11 1.67 1.97 1.27 4.38 0.52 ...
 $ X2: num  3.66 3.75 3.32 2.68 4.97 2.98 4.56 1.78 2.83 6.36 ...
 $ X3: num  1.41 3.01 0.34 2.19 0.25 1.69 2.01 1.25 4.41 0.43 ...
 $ X4: num  -3.53 -4.46 -3.35 -7.54 -7.02 -2.53 -6.1 -5.99 -3.92 -5.84 ...
 $ X5: num  0.57 0.01 -0.7 1.62 -0.95 -1.37 1.18 -0.72 2.63 -1.63 ...
 $ X6: num  -0.82 1.2 3.03 -0.91 -1.6 1.77 1 -1.33 1.31 -0.7 ...
 $ X7: num  1.01 0.06 2.02 3.63 2.66 2.53 1.29 3.5 1.17 1.8 ...
 $ X8: num  8.31 8.52 9.78 7.34 7.26 ...
 $ X9: num  6.04 6.53 7.52 5.61 4.52 6.06 6.2 5.99 6.93 4.38 ...

Je vais mettre ici les questions que je posais, et donner des pistes de réflexions, non pas sur les réponses attendues (je n’attends rien de cet exercice à part une réflexion), mais sur la discussion que peut amener chacune des questions,

  • Faîtes un modèle linéaire pour expliquer Y en fonction des neuf variables explicatives. Combien de variables explicatives garderiez-vous ?

Commençons par faire une régression sur toutes les variables

summary(lm(Y~., data=df))
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.046133   2.193448   3.668 0.000349 ***
X1           0.342293   0.915894   0.374 0.709186    
X2          -0.040479   0.103409  -0.391 0.696073    
X3           1.683875   0.897278   1.877 0.062693 .  
X4          -0.009254   0.062382  -0.148 0.882295    
X5          -1.085367   0.113840  -9.534  < 2e-16 ***
X6           0.983207   0.111830   8.792 5.49e-15 ***
X7          -0.015646   0.087483  -0.179 0.858327    
X8           0.012165   0.094756   0.128 0.898033    
X9           0.172210   0.180605   0.954 0.342005    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.019 on 137 degrees of freedom
Multiple R-squared:  0.9228,	Adjusted R-squared:  0.9178 
F-statistic:   182 on 9 and 137 DF,  p-value: < 2.2e-16

On a, dans la sortie, une dizaine de tests de Student qui sont évoqués, correspondant au test de l’hypothèse H_0:\beta_j=0 (contre l’hypothèse alternative (bilatérale) H_1:\beta_j\neq 0), dans un modèle de la forme y_i=\beta_0+\beta_1x_{1,i}+\dots+\beta_9x_{9,i}+\varepsilon_iC’est ce qu’on appelle le test de significativité de la variable x_j (oui, comme on l’a vu en cours, on peut dire le test parce que les autres tests classiques – Fisher, ou Wald – sont équivalent – sauf qu’au lieu de regarder t on regarde la statistique t^2 – ce qui présente l’avantage de voir la statistique de test comme une forme de distance à l’hypothèse H_0 : si c’est trop grand, on rejette…). Avec un seuil d’acceptation de l’ordre de 5%, on nous dit que 6 variables ne sont pas significatives. Mais gardons bien en mémoire que le test de significativité de x_j est fait ici en supposant que toutes les autres variables restent dans le modèle. Autre chose: avec un seuil d’acceptation de l’ordre de 10%, une des variables (la troisième) peut être vue comme significative.

Faisons un test multiple, pour savoir si on peut supprimer 6 des 9 variables explicatives possibles (faisons le à la main, inutile d’aller chercher un package pour le faire)

reg1 = lm(formula = Y ~ ., data = df)
reg0 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg0)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.11031    0.16788   54.27   <2e-16 ***
X3           1.96784    0.07604   25.88   <2e-16 ***
X5          -0.99296    0.08751  -11.35   <2e-16 ***
X6           1.08761    0.05934   18.33   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: < 2.2e-16 anova(reg0,reg1) Analysis of Variance Table Model 1: Y ~ X3 + X5 + X6 Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 Res.Df RSS Df Sum of Sq F Pr(>F)
1    143 146.16                           
2    137 142.12  6    4.0338 0.6481 0.6916

Le test de Fisher nous dit qu’on peut accepter l’hypothèse que les 6 coefficients sont nuls – ici H_0:\beta_1=\beta_2=\beta_4=\beta_7=\beta_8=\beta_9=0(qui est un test multiple, contrairement au test de Student précédant). Mais il nous dit aussi, qu’individuellement, les trois variables restants semblent significatives. Donc j’aurais tendance à garder 3 variables explicatives (je ne parle pas de la constante : on garde toujours la constante – qui n’explique pas grand chose, sauf la valeur moyenne de y).

  • faites une prévision pour un individu dont on sait que X3=1, X5=1 et X6=8

(en réalité, la vraie question que j’ai posée contenait une typo ce qui la rendait vicieuse parce que ce n’est pas le modèle qu’on vient de construire… mais pour commencer, regardons cette question)

On vient de calibrer ce modèle, donc il suffit de faire une prévision,

predict(reg0, newdata=data.frame(X3=1, X5=1, X6=8))
       1 
18.78603

Mais ce n’était pas la vraie question…

  • faites une prévision pour un individu dont on sait que X1=3, X5=1 et X6=8

Je pense que pour répondre à cette question, il convient d’oublier tout ce qu’on vient de voir. On nous donne trois informations, et on va voir si on peut les exploiter. Autrement dit, on va commencer par regarder la régression sur ces trois variables,

reg156 = lm(formula = Y ~ X1+X5+X6, data = df)
summary(reg156)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.06688    0.17222   52.65   <2e-16 ***
X1           1.98342    0.07800   25.43   <2e-16 ***
X5          -1.01438    0.08947  -11.34   <2e-16 ***
X6           1.07824    0.06039   17.86   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.026 on 143 degrees of freedom
Multiple R-squared:  0.9183,	Adjusted R-squared:  0.9166 
F-statistic: 535.8 on 3 and 143 DF,  p-value: < 2.2e-16

Le modèle est ici bon, les trois variables étant significatives. Cela dit, je dis qu’il est “bon” mais on verra vendredi comment discuter davantage ce point… En tous cas, on peut tenter de faire une prévision, et on obtient

predict(reg156, newdata=data.frame(X1=3, X5=1, X6=8))
       1 
22.62867
  • faites une prévision pour un individu dont on sait que X1=3, X3=2, X5=1 et X6=8

Comme auparavant, on nous donne 4 informations… on va regarder le modèle avec les 4 variables

reg1356 = lm(formula = Y ~ X1+X3+X5+X6, data = df)
summary(reg1356)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.10458    0.17131  53.148   <2e-16 ***
X1           0.16362    0.89028   0.184    0.854    
X3           1.80663    0.88052   2.052    0.042 *  
X5          -0.99558    0.08895 -11.192   <2e-16 ***
X6           1.08649    0.05986  18.152   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.014 on 142 degrees of freedom
Multiple R-squared:  0.9207,	Adjusted R-squared:  0.9184 
F-statistic: 411.9 on 4 and 142 DF,  p-value: < 2.2e-16

Cette fois, une des variables n’est pas significative (la première) donc on devrait l’enlever: on fait alors la régression juste sur les trois autres variables

reg356 = lm(formula = Y ~ X3+X5+X6, data = df)
summary(reg356)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.11031    0.16788   54.27   <2e-16 ***
X3           1.96784    0.07604   25.88   <2e-16 ***
X5          -0.99296    0.08751  -11.35   <2e-16 ***
X6           1.08761    0.05934   18.33   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.011 on 143 degrees of freedom
Multiple R-squared:  0.9206,	Adjusted R-squared:  0.919 
F-statistic:   553 on 3 and 143 DF,  p-value: < 2.2e-16

Cette fois, le modèle est “bon” et on peut alors faire la prévision

predict(reg356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.75388

A titre d’information, si on avait fait la prévision sans enlever la variable non significative, on aurait obtenu la valeur suivante

predict(reg1356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8))
       1 
20.90501

qui est globalement assez proche. On aura l’occasion d’en reparler en cours, quand on se demandera s’il est plus grave d’avoir une valeur non-significative dans la régression, ou d’oublier une variable importante… Mais encore une fois, le but de ces petits exercices est d’appliquer ce qu’on a vu en cours, et d’introduire des questions auxquelles j’apporterai des réponses au prochain cours !

Combiner les modalités d’une variable factorielle

Un billet rapide pour reprendre un point qu’on a vu ce matin en cours STT5100 pour illustrer le test de Fisher. On va utiliser les données de prix d’appartements en Pologne (données pas mal utilisées dans mon ébauche de notes de cours)

library(DALEX)
data(apartments)
with(data = apartments, boxplot(m2.price ~ district))

On a envie de faire ici des regroupements de modalités (c’est d’ailleurs suggéré par la régression simple, 5 variables explicatives étant ici non significatives). Pour mieux voir, on peut réordonner les modalités

A = with(data = apartments, aggregate(m2.price,by=list(district),FUN=mean))
A = A[order(A$x),]
L = as.character(A$Group.1)
apartments$district = factor(apartments$district, level=L)
with(data = apartments, boxplot(m2.price ~ district))

On va prendre ici le district le moins cher comme référence,

reg=lm(m2.price ~ district, data=apartments)
> summary(reg)
 
Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          2968.36      58.02  51.160   <2e-16 ***
districtBielany        17.38      84.16   0.207    0.836    
districtPraga          26.45      85.12   0.311    0.756    
districtUrsynow        42.01      82.65   0.508    0.611    
districtBemowo         80.10      83.71   0.957    0.339    
districtUrsus         102.01      82.25   1.240    0.215    
districtZoliborz      829.59      83.94   9.884   <2e-16 ***
districtMokotow       887.10      81.86  10.837   <2e-16 ***
districtOchota        987.93      84.16  11.738   <2e-16 ***
districtSrodmiescie  2214.39      83.28  26.591   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 597.4 on 990 degrees of freedom
Multiple R-squared:  0.5698,	Adjusted R-squared:  0.5659 
F-statistic: 145.7 on 9 and 990 DF,  p-value: < 2.2e-16

On peut tester si les 5 premières modalités sont nulles, ce qui est un test multiple, et on va utiliser le test de Ficher :

library(car)
linearHypothesis(reg, c("districtBielany = 0", 
                        "districtPraga = 0",
                        "districtUrsynow = 0",
                        "districtBemowo = 0",
                        "districtUrsus = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F Pr(>F)
1    995 354051715                           
2    990 353269202  5    782513 0.4386 0.8217

La statistique de Fisher est faible, et avec une p-value de 82%. On peut tenter le diable, et rajouter encore une modalité

library(car)
linearHypothesis(reg, c("districtBielany = 0", 
                        "districtPraga = 0",
                        "districtUrsynow = 0",
                        "districtBemowo = 0",
                        "districtUrsus = 0",
                        "districtZoliborz = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    996 405455409                                  
2    990 353269202  6  52186207 24.374 < 2.2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Mais on a peut être été trop gourmand cette fois. On va regrouper les 6 premières modalités (et appeler A le regroupement de districts). Si on regarde les prix moyens, par districts, on obtient

levels(apartments$district) = c(rep("A",6),levels(apartments$district)[7:11])
with(data = apartments, boxplot(m2.price ~ district))

apartments$district = relevel(apartments$district,"Zoliborz")

On recommence, en mettant le district le moins cher comme référence, et on veut tester si les deux suivants ont des coefficients nuls dans la régression linéaire.

reg=lm(m2.price ~ district, data=apartments)
linearHypothesis(reg, c("districtMokotow = 0",
                        "districtOchota = 0"))
Linear hypothesis test
 
Model 1: restricted model
Model 2: m2.price ~ district
 
  Res.Df       RSS Df Sum of Sq      F Pr(>F)
1    997 355292524                           
2    995 354051715  2   1240809 1.7435 0.1754

Avec une p-value de 17%, on peut accepter de regrouper les trois modalités ensemble. On a alors trois groupes de districts, dont les noms sont A, B et C. On obtient les boîtes à moustaches suivantes

levels(apartments$district) = c("B","A",rep("B",2),"C")
apartments$district = relevel(apartments$district,"A")
with(data = apartments, boxplot(m2.price ~ district))

Je laisse les plus courageux vérifier, mais on a trois districts effectivement différents, et si le but est de prévoir le prix des logements, inutiles d’utiliser un découpage avec 10 modalités, un découpage avec 3 suffit !

Classification from scratch, linear discrimination 8/8

Eighth post of our series on classification from scratch. The latest one was on the SVM, and today, I want to get back on very old stuff, with here also a linear separation of the space, using Fisher’s linear discriminent analysis.

Bayes (naive) classifier

Consider the follwing naive classification rulem^\star(\mathbf{x})=\text{argmin}_y\{\mathbb{P}[Y=y\vert\mathbf{X}=\mathbf{x}]\}orm^\star(\mathbf{x})=\text{argmin}_y\left\{\frac{\mathbb{P}[\mathbf{X}=\mathbf{x}\vert Y=y]}{\mathbb{P}[\mathbf{X}=\mathbf{x}]}\right\}(where \mathbb{P}[\mathbf{X}=\mathbf{x}] is the density in the continuous case).

In the case where y takes two values, that will be standard \{0,1\} here, one can rewrite the later asm^\star(\mathbf{x})=\begin{cases}1\text{ if }\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})>\displaystyle{\frac{1}{2}}\\0\text{ otherwise}\end{cases}and the set\mathcal{D}_S =\left\{\mathbf{x},\mathbb{E}(Y\vert \mathbf{X}=\mathbf{x})=\frac{1}{2}\right\}is called the decision boundary.

Assume that\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})and\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})then explicit expressions can be derived.m^\star(\mathbf{x})=\begin{cases}1\text{ if }r_1^2< r_0^2+2\displaystyle{\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}+\log\frac{\vert\mathbf{\Sigma}_0\vert}{\vert\mathbf{\Sigma}_1\vert}}\\0\text{ otherwise}\end{cases}where r_y^2 is the Manalahobis distance, r_y^2 = [\mathbf{X}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[\mathbf{X}-\mathbf{\mu}_y]

Let \delta_ybe defined as\delta_y(\mathbf{x})=-\frac{1}{2}\log\vert\mathbf{\Sigma}_y\vert-\frac{1}{2}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]^{\text{{T}}}\mathbf{\Sigma}_y^{-1}[{\color{blue}{\mathbf{x}}}-\mathbf{\mu}_y]+\log\mathbb{P}(Y=y)the decision boundary of this classifier is \{\mathbf{x}\text{ such that }\delta_0(\mathbf{x})=\delta_1(\mathbf{x})\}which is quadratic in {\color{blue}{\mathbf{x}}}. This is the quadratic discriminant analysis. This can be visualized bellow.

The decision boundary is here

But that can’t be the linear discriminant analysis, right? I mean, the frontier is not linear… Actually, in Fisher’s seminal paper, it was assumed that \mathbf{\Sigma}_0=\mathbf{\Sigma}_1.

In that case, actually, \delta_y(\mathbf{x})={\color{blue}{\mathbf{x}}}^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y-\frac{1}{2}\mathbf{\mu}_y^{\text{T}}\mathbf{\Sigma}^{-1}\mathbf{\mu}_y+\log\mathbb{P}(Y=y) and the decision frontier is now linear in {\color{blue}{\mathbf{x}}}. This is the linear discriminant analysis. This can be visualized bellow

Here the two samples have the same variance matrix and the frontier is

Link with the logistic regression

Assume as previously that\mathbf{X}\vert Y=0\sim\mathcal{N}(\mathbf{\mu}_0,\mathbf{\Sigma})and\mathbf{X}\vert Y=1\sim\mathcal{N}(\mathbf{\mu}_1,\mathbf{\Sigma})then\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}is equal to \mathbf{x}^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_y]-\frac{1}{2}[\mathbf{\mu}_1-\mathbf{\mu}_0]^{\text{{T}}}\mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0]+\log\frac{\mathbb{P}(Y=1)}{\mathbb{P}(Y=0)}which is linear in \mathbf{x}\log\frac{\mathbb{P}(Y=1\vert \mathbf{X}=\mathbf{x})}{\mathbb{P}(Y=0\vert \mathbf{X}=\mathbf{x})}=\mathbf{x}^{\text{{T}}}\mathbf{\beta}Hence, when each groups have Gaussian distributions with identical variance matrix, then LDA and the logistic regression lead to the same classification rule.

Observe furthermore that the slope is proportional to \mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0], as stated in Fisher’s article. But to obtain such a relationship, he observe that the ratio of between and within variances (in the two groups) was\frac{\text{variance between}}{\text{variance within}}=\frac{[\mathbf{\omega}\mathbf{\mu}_1-\mathbf{\omega}\mathbf{\mu}_0]^2}{\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_1\mathbf{\omega}+\mathbf{\omega}^{\text{T}}\mathbf{\Sigma}_0\mathbf{\omega}}which is maximal when \mathbf{\omega} is proportional to \mathbf{\Sigma}^{-1}[\mathbf{\mu}_1-\mathbf{\mu}_0], when \mathbf{\Sigma}_0=\mathbf{\Sigma}_1.

Homebrew linear discriminant analysis

To compute vector \mathbf{\omega}

m0 = apply(myocarde[myocarde$PRONO=="0",1:7],2,mean)
m1 = apply(myocarde[myocarde$PRONO=="1",1:7],2,mean)
Sigma = var(myocarde[,1:7])
omega = solve(Sigma)%*%(m1-m0)
omega
                 [,1]
FRCAR -0.012909708542
INCAR  1.088582058796
INSYS -0.019390084344
PRDIA -0.025817110020
PAPUL  0.020441287970
PVENT -0.038298291091
REPUL -0.001371677757

For the constant – in the equation \omega^T\mathbf{x}+b=0 – if we have equiprobable probabilities, use

b = (t(m1)%*%solve(Sigma)%*%m1-t(m0)%*%solve(Sigma)%*%m0)/2

Application (on the small dataset)

In order to visualize what’s going on, consider the small dataset, with only two covariates,

x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85)
y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3)
z = c(1,1,1,1,1,0,0,1,0,0)
df = data.frame(x1=x,x2=y,y=as.factor(z))
m0 = apply(df[df$y=="0",1:2],2,mean)
m1 = apply(df[df$y=="1",1:2],2,mean)
Sigma = var(df[,1:2])
omega = solve(Sigma)%*%(m1-m0)
omega
         [,1]
x1 -2.640613174
x2  4.858705676


Using R regular function, we get

library(MASS)
fit_lda = lda(y ~x1+x2 , data=df)
fit_lda
 
Coefficients of linear discriminants:
            LD1
x1 -2.588389554
x2  4.762614663

which is the same coefficient as the one we got with our own code. For the constant, use

b = (t(m1)%*%solve(Sigma)%*%m1-t(m0)%*%solve(Sigma)%*%m0)/2

If we plot it, we get the red straight line

plot(df$x1,df$x2,pch=c(1,19)[1+(df$y=="1")])
abline(a=b/omega[2],b=-omega[1]/omega[2],col="red")


As we can see (with the blue points), our red line intersects the middle of the segment of the two barycenters

points(m0["x1"],m0["x2"],pch=4)
points(m1["x1"],m1["x2"],pch=4)
segments(m0["x1"],m0["x2"],m1["x1"],m1["x2"],col="blue")
points(.5*m0["x1"]+.5*m1["x1"],.5*m0["x2"]+.5*m1["x2"],col="blue",pch=19)

Of course, we can also use R function

predlda = function(x,y) predict(fit_lda, data.frame(x1=x,x2=y))$class==1
vv=outer(vu,vu,predlda)
contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5)


One can also consider the quadratic discriminent analysis since it might be difficult to argue that \mathbf{\Sigma}_0=\mathbf{\Sigma}_1

fit_qda = qda(y ~x1+x2 , data=df)

The separation curve is here

plot(df$x1,df$x2,pch=19,
col=c("blue","red")[1+(df$y=="1")])
predqda=function(x,y) predict(fit_qda, data.frame(x1=x,x2=y))$class==1
vv=outer(vu,vu,predlda)
contour(vu,vu,vv,add=TRUE,lwd=2,levels = .5)

Classification from scratch, logistic regression 1/8

Let us start today our series on classification from scratch

The logistic regression is based on the assumption that given covariates \mathbf{x}, Y has a Bernoulli distribution,Y|\mathbf{X}=\mathbf{x}\sim\mathcal{B}(p_{\mathbf{x}}),~~~~p_\mathbf{x}=\frac{\exp[\mathbf{x}^T\mathbf{\beta}]}{1+\exp[\mathbf{x}^T\mathbf{\beta}]}The goal is to estimate parameter \mathbf{\beta}.

Recall that the heuristics for the use of that function for the probability is that\log[\text{odds}(Y=1)]=\log\frac{\mathbb{P}[Y=1]}{\mathbb{P}[Y=0]}=\mathbf{x}^T\mathbf{\beta}

Maximimum of the (log)-likelihood function

The log-likelihood is here\log\mathcal{L} = \sum_{i=1}^n y_i\log p_i+(1-y_i)\log (1-p_i) where p_{i}=(1+\exp[-\mathbf{x}_i^T\mathbf{\beta}])^{-1}. Numerical techniques are based on (numerical) gradient descent to compute the maximum of the likelihood function. The (negative) log-likelihood is the following function

y = myocarde$PRONO
X = cbind(1,as.matrix(myocarde[,1:7]))
negLogLik = function(beta){
 -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta)))
 }

We use the minus sign since standard optimization routines compute minima, not maxima. Now, to find the minimum of that function, we need a starting point to initiate the algorithm

beta_init = lm(PRONO~.,data=myocarde)$coefficients

Why not start with the parameter of the OLS. Somehow, we might think that at least, sign should be ok for instance. Anyway, we need a starting point, and let us use that one.

logistic_opt = optim(par = beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9))

Here, we obtain

 logistic_opt$par
 (Intercept)        FRCAR        INCAR        INSYS    
 1.656926397  0.045234029 -2.119441743  0.204023835 
       PRDIA        PAPUL        PVENT        REPUL 
-0.102420095  0.165823647 -0.081047525 -0.005992238

Let us verify here that this output is valid. For instance, what if we change the value of the starting point (randomly)

simu = function(i){
logistic_opt_i = optim(par = rnorm(8,0,3)*beta_init, 
negLogLik, hessian=TRUE, method = "BFGS", 
control=list(abstol=1e-9))
logistic_opt_i$par[2:3]
}
v_beta = t(Vectorize(simu)(1:1000))
plot(v_beta)
par(mfrow=c(1,2))
hist(v_beta[,1],xlab=names(myocarde)[1])
hist(v_beta[,2],xlab=names(myocarde)[2])

Ooops. There is a problem here. Clearly, we cannot rely on numerical optimization here. We can think about using another optimization routine

library(optimx)
logit = function(mX, vBeta) {
  exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) 
}
logLikelihoodLogitStable = function(vBeta, mX, vY) {
  -sum(vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta))) + 
(1-vY)*(-log(1 + exp(mX %*% vBeta)))) 
}
likelihoodScore = function(vBeta, mX, vY) {
  return(t(mX) %*% (logit(mX, vBeta) - vY) )
}
optimLogitLBFGS = optimx(beta_init, logLikelihoodLogitStable, 
method = 'L-BFGS-B', gr = likelihoodScore, 
mX = X, vY = y, hessian=TRUE)

The optimum is here

attr(optimLogitLBFGS, "details")[[2]]
              [,1]
       0.066680272
FRCAR  0.003080542
INCAR  0.079031364
INSYS -0.001586194
PRDIA  0.040500697
PAPUL -0.041870705
PVENT -0.014162756
REPUL  0.195632244

Let’s be honest here, I do not feel confortable with those techniques. So, what happened here ?

Here, the technique we use is based on the following idea,\mathbf{\beta}_{new}=\mathbf{\beta}_{old} -\left(\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}\right)^{-1}\cdot \frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}The problem is that my computer does not know this first and second derivatives. So it will compute them using approximation techniques.

Actually, it is possible to use functions dedicated to such computation

library(numDeriv)
library(MASS)
logit = function(x){1/(1+exp(-x))}
logLik = function(beta, X, y){
 -sum(y*log(logit(X%*%beta)) + 
(1-y)*log(1-logit(X%*%beta)))
}
optim_second = function(beta, num_iter){
  LL = vector()
  for(i in 1:num_iter){
    grad = (t(X)%*%(logit(X%*%beta) - y)) 
    H = hessian(logLik, beta, method = "complex", X = X, y = y)
    beta = beta - ginv(H)%*%grad
    LL[i] = logLik(beta, X, y)
  }
  result = list(beta, H)
return(result)
}

With our OLS starting point, we obtain

opt0 = optim_second(beta_init,500)
opt0[[1]]
             [,1]
[1,]  0.951074420
[2,]  0.018860280
[3,]  0.275428978
[4,]  0.144803636
[5,] -0.058535606
[6,]  0.001182178
[7,] -0.108651776
[8,] -0.002940315

But if we try with another starting point

opt1 = optim_second(beta_init*runif(8),500)
opt1[[1]]
             [,1]
[1,]  0.052894794
[2,]  0.024718435
[3,]  0.167953661
[4,]  0.171662947
[5,] -0.057458066
[6,] -0.011361034
[7,] -0.107532114
[8,] -0.002679064

Clearly, some coefficients are rather close. But other aren’t. From my point of viezw, that is a major problem (keep in mind that we do not deal here with massive data ! There are only 7 explanatory variables, and only 71 observations).

Why not try to be clever, and use the analytical values of those derivatives ? Even if some people claim the oppositive, sometimes, it can actually be usefull to do the maths, instead of considering only numerical values.

Newton (or Fisher) Algorithm

If you open any Econometrics textbooks (one can also try to derive it), you will get \frac{\partial\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}}=\mathbf{X}^T(\mathbf{y}-\mathbf{p}_{old})
while\frac{\partial^2\log\mathcal{L}(\mathbf{\beta}_{old})}{\partial\mathbf{\beta}\partial\mathbf{\beta}^T}=-\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X}

Y=myocarde$PRONO
X=cbind(1,as.matrix(myocarde[,1:7]))
colnames(X)=c("Inter",names(myocarde[,1:7]))
 beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)
 for(s in 1:9){
   pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s]))
   gradient=t(X)%*%(Y-pi)
   omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))
   Hessian=-t(X)%*%omega%*%X
   beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)}

Observe that here, I use only ten iterations of the algorithm !

 beta[,8:10]
                [,1]          [,2]          [,3]
XInter -10.187641685 -10.187641696 -10.187641696
XFRCAR   0.138178119   0.138178119   0.138178119
XINCAR  -5.862429035  -5.862429037  -5.862429037
XINSYS   0.717084018   0.717084018   0.717084018
XPRDIA  -0.073668171  -0.073668171  -0.073668171
XPAPUL   0.016756506   0.016756506   0.016756506
XPVENT  -0.106776012  -0.106776012  -0.106776012
XREPUL  -0.003154187  -0.003154187  -0.003154187

The thing is that is seems to converge extremely fast. And it is rather robust ! Look at what we get if we change our starting point

beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8)
 for(s in 1:9){
   pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s]))
   gradient=t(X)%*%(Y-pi)
   omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))
   Hessian=-t(X)%*%omega%*%X
   beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)}
 beta[,8:10]
                [,1]          [,2]          [,3]
XInter -10.187641586 -10.187641696 -10.187641696
XFRCAR   0.138178118   0.138178119   0.138178119
XINCAR  -5.862429017  -5.862429037  -5.862429037
XINSYS   0.717084013   0.717084018   0.717084018
XPRDIA  -0.073668172  -0.073668171  -0.073668171
XPAPUL   0.016756508   0.016756506   0.016756506
XPVENT  -0.106776012  -0.106776012  -0.106776012
XREPUL  -0.003154187  -0.003154187  -0.003154187

Nice, isn’t it? Looks like we got our winner, don’t we? And one can use the inverse of the Hessian matrix to get standard deviations.

Weighted Least-Squares

Let us go one step further. We’ve seen that we want to compute something like\mathbf{\beta}_{new} =(\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Delta}_{old}\mathbf{z}(if we do substitute matrices in the analytical expressions) where \mathbf{z}=\mathbf{X}\mathbf{\beta}_{old}+\mathbf{\Delta}_{old}^{-1}[\mathbf{y}-\mathbf{p}_{old}]. But actually, that’s simply a standard least-square problem\mathbf{\beta}_{new} = \text{argmin}\left\lbrace(\mathbf{z}-\mathbf{X}\mathbf{\beta})^T\mathbf{\Delta}_{old}^{-1}(\mathbf{z}-\mathbf{X}\mathbf{\beta})\right\rbraceThe only problem here is that weights \mathbf{\Delta}_{old} are functions of unknown \mathbf{\beta}_{old}. But actually, if we keep iterating, we should be able to solve it : given the \mathbf{\beta} we got the weights, and with the weights, we can use weighted OLS to get an updated \mathbf{\beta}. That’s the idea of iteratively reweighted least squares.

The algorithm will be

df = myocarde
beta_init = lm(PRONO~.,data=df)$coefficients
X = cbind(1,as.matrix(myocarde[,1:7]))
beta = beta_init
for(s in 1:1000){
p = exp(X %*% beta) / (1+exp(X %*% beta))
omega = diag(nrow(df))
diag(omega) = (p*(1-p))
df$Z = X %*% beta + solve(omega) %*% (df$PRONO - p)
beta = lm(Z~.,data=df[,-8], weights=diag(omega))$coefficients
}

and the output is here

 beta
  (Intercept)         FRCAR         INCAR         INSYS         PRDIA 
-10.187641696   0.138178119  -5.862429037   0.717084018  -0.073668171 
        PAPUL         PVENT         REPUL 
  0.016756506  -0.106776012  -0.003154187

which is almost what we’ve obtained before. Nice isn’t it ? Actually, here we also have standard deviations of estimators

summary( lm(Z~.,data=df[,-8], weights=diag(omega)))
 
Coefficients:
              Estimate Std. Error t value Pr(&gt;|t|)
(Intercept) -10.187642  10.668138  -0.955    0.343
FRCAR         0.138178   0.102340   1.350    0.182
INCAR        -5.862429   6.052560  -0.969    0.336
INSYS         0.717084   0.503527   1.424    0.159
PRDIA        -0.073668   0.261549  -0.282    0.779
PAPUL         0.016757   0.306666   0.055    0.957
PVENT        -0.106776   0.099145  -1.077    0.286
REPUL        -0.003154   0.004386  -0.719    0.475

The standard glm function

Of course, it is possible to use an R built-in function to get our estimate

summary(glm(PRONO~.,data=myocarde,family=binomial(link = "logit")))
 
Coefficients:
              Estimate Std. Error z value Pr(&gt;|z|)
(Intercept) -10.187642  11.895227  -0.856    0.392
FRCAR         0.138178   0.114112   1.211    0.226
INCAR        -5.862429   6.748785  -0.869    0.385
INSYS         0.717084   0.561445   1.277    0.202
PRDIA        -0.073668   0.291636  -0.253    0.801
PAPUL         0.016757   0.341942   0.049    0.961
PVENT        -0.106776   0.110550  -0.966    0.334
REPUL        -0.003154   0.004891  -0.645    0.519

Application and visualisation

Let us visualize the prediction obtained from the logistic regression, on our second dataset

x = c(.4,.55,.65,.9,.1,.35,.5,.15,.2,.85)
y = c(.85,.95,.8,.87,.5,.55,.5,.2,.1,.3)
z = c(1,1,1,1,1,0,0,1,0,0)
df = data.frame(x1=x,x2=y,y=as.factor(z))
reg = glm(y~x1+x2,data=df,family=binomial(link = "logit"))
u = seq(0,1,length=101)
p = function(x,y) predict.glm(reg,newdata=data.frame(x1=x,x2=y),type="response")
v = outer(u,u,p)
image(u,u,v,xlab="Variable 1",ylab="Variable 2",col=clr10,breaks=(0:10)/10)
points(x,y,pch=19,cex=1.5,col="white")
points(x,y,pch=c(1,19)[1+z],cex=1.5)
contour(u,u,v,levels = .5,add=TRUE)


Here level curves – or iso-probabilities – are linear, so the space is divided in two (0 and 1, survival and death, white and black) by a straight line (or an hyperplane in higher dimension). Furthermore, since we have a linear model, if we change the cutoff (the threshold used to create the two classes), we obtain another straight line (or hyperplane) parallel to the first one.

Next time, we will introduce splines to smooth those continuous covariates… to be continued.

La guerre des étoiles : distinguer le signal du bruit

La grande difficulté dans la modélisation et la construction de modèles prédictifs est de réussir à distinguer le signal et le bruit (pour reprendre le titre du classique de Nate Silver). La réponse statistique est la notion de significativité, et la recherche des ‘étoiles’ dans les sorties de régression. Avec l’explosion du nombre de données, il est devenu crucial de faire cette distinction, de savoir quelles sont les interactions qui sont significatives.

Approches historiques de cette notion de « significativité »

Le débat sur la significativité est ancien, même si sa formulation s’est faite historiquement dans des termes assez vagues. Par exemple, dès 1710, le médecin et mathématicien John Arbuthnot s’était interrogé sur le ratio du nombre de naissance de garçons et de filles, se demandant si la différence était « statistiquement significative ». Pour être plus précis, en utilisant des statistiques sur près de 90 ans, il avait noté “There seems no more probable Cause to be assigned in Physics for this Equality of the Births, than that in our ’first Parents Seed there were at first formed an equal Number of both Sexes.John Arbuthnot pose le premier la question en terme probabilistes. Un siècle plus tard, Pierre-Simon de Laplace a présenté ce que l’on peut interpréter comme un « test de significativité », (là encore avec notre terminologie actuelle). Il avait en effet noté, en prenant des mesures sur des baromètres que les observations à 9 heures du matin et 4 heures de l’après-midi était différentes. Significativement différentes. Et là encore posé la question en terme probabiliste, en se demandant s’il est « extrêmement probable » qu’il y ait une différence entre les deux mesures. Il avait alors introduit, le premier, un test de comparaison entre des valeurs moyennes : ayant noté que la différence excédait plusieurs écart-types, ce qu’il jugeait alors significativement important, il en conclue que les séries sont significativement différentes. En 1885, Francis Edgeworth avait repris ces idées, en comparant la taille des criminels, et la taille des gens ordinaires. Mais il faudra attendre les travaux de William Gosset (plus connu sous le nom Student), de Karl Pearson et surtout de Ronald Fisher pour avoir une définition plus rigoureuse de la significativité.

La notion de significativité est cruciale dans la construction de modèles prédictifs. En assurance automobile, l’âge du conducteur est une variable « significative » quand il s’agit d’expliquer la fréquence de sinistres. Formellement, cela signifie que l’âge et la fréquence de sinistre sont corrélées, et que cette corrélation, notée R, est « significativement non-nulle ». En 1921, Ronald Fisher, en proposant la construction d’un intervalle de confiance, propose ainsi un test de significativité. Comme il l’écrit, « from these values, we obtain the difference 0.0471 ± [0.0142] which might well be regarded as significant ». Sur la Figure 1, la figure de droite correspond au cas où la valeur est significativement non nulle (car l’intervalle de confiance à 95% ne contient pas 0) alors que la figure de gauche correspond au cas où la valeur est non-significative (car l’intervalle de confiance à 95% contient 0). Ronald Fisher passera alors plusieurs années à formaliser et expliquer cette idée de significativité statistique.

Figure 1 : distribution théorique de l’estimateur de la corrélation, R, et intervalle de confiance.

Tests, prise de décision et erreurs

Une des contributions majeures de la statistique des années 1920 a été de formaliser la prise de décision. La Figure 2 montre ainsi le mécanisme binaire de la prise de décision, et les deux types d’erreur : rejeter à tort une hypothèse, ou accepter à tort une hypothèse.

Figure 2 : schéma de la prise de décision.

Déclarer qu’une variable, ou une différence, est significative peut conduire à deux types d’erreur : la déclarer comme non-significative, alors qu’elle l’était (prendre un signal pour du bruit) ; et la déclarer comme significative, alors qu’elle ne l’était (prendre du bruit pour un signal). Quand on construit un test, on va essayer de contrôler la probabilité de commettre de telles erreurs.

Mais de faibles taux d’erreur ne veulent pas forcément dire qu’un test est efficace, et le principal danger des tests médicaux est lié à une mauvaise interprétation de ces taux d’erreur. Supposons qu’une maladie touche une personne sur mille, dans une population (selon le test, i.e. 1 personne sur 1000 est ‘positive’). Supposons aussi qu’il soit relative « efficace » au sens où 90% des cas ‘positifs’ sont effectivement atteints, et le test est négatif dans 99% des cas quand est n’est pas touché par la maladie. Ces chiffres sont bien au-dessus de la plupart des principaux tests couramment utilisés. Si le test est effectué sur 10000 personnes, le test sera positif sur 10 personnes (en moyenne), parmi lesquelles 9 sont effectivement malade, mais 1 est saine. A côté, le test sera négatif pour 9990 personnes. Et parmi ces 9990 négatifs, 9890 sont effectivement saines, et une centaine est pourtant malade. Donc au final, sur les 109 personnes malades, 9 ont été détectées, mais 100 sont passés inaperçues, soit un peu plus de 90% des malades ! Ce test supposé efficace ne l’est peut-être pas tant que ça.

Figure 3 : calculs des taux d’erreurs

Le concept de « p-value » et le mythe des 5%

Le concept de p-value est lié justement à l’erreur de rejeter, à tort une hypothèse. Si Ronald Fisher en parle abondement, il faut toutefois attribuer la paternité du concept à son collègue, Karl Pearson. Ce dernier, pour définir la notion de « significativité », utilise la formulation suivante : «  P=.1227 or the odds are now only 8 to 1 against a system of deviations as improbable as or more improbable than this one ». Pour reprendre l’exemple de la corrélation, Pearson nous dit qu’il y a une chance sur 8 pour obtenir une valeur aussi improbable. Comme l’illustre la Figure 4, la p-value est la probabilité d’avoir une statistique aussi grande ou aussi petite que celle obtenue sur l’échantillon, si effectivement la corrélation était nulle. Ou encore p=P{|R| > r | H0],

Figure 4 : distribution théorique de l’estimateur de la corrélation sous l’hypothèse où cette dernière serait nulle et probabilité que |R| dépasse la valeur empirique, observée.

Sur la figure de gauche, la p-value correspondant à l’aire rouge est de l’ordre de 10% alors que sur la figure de droite l’aire est de 0.1%. En fait, si la p-value est inférieur à 5%, alors 0 n’est pas dans l’intervalle de confiance à 95% de l’estimateur de corrélation. Donner la p-value est alors suffisant pour juger de la significativité d’une statistique.

C’est dans le chapitre 4 de son livre que Ronald Fisher pose les bases de la pratique (toujours en vigueur aujourd’hui) des tests statistiques : « if p is between 10% and 90% there is certainly no reason to suspect the hypothesis tested. If it is below 2% it is strongly indicated that the hypothesis fails to account for the whole of facts ». Autrement dit, avec une p-value supérieure à 10%, on peut accepter notre hypothèse (souvent que la variable n’est pas significative) et avec une p-value inférieure à 2%, on va la rejeter (la variable sera alors significative, si on fait un test de nullité d’une corrélation). Mais entre les deux ? Ronald Fisher clôt sans vraiment s’en rendre compte le débat, en affirmant « we whall not often be astray if we draw a conventional line at 5% ». Le choix – aujourd’hui dogmatique – d’un seuil à 5% repose sur l’idée de rejeter à tort une hypothèse avec une chance sur 20, ou bien correspond au fait de s’éloigner de 2 écart-types de la moyenne d’une loi normale centrée-réduite (ce qui arrive avec une chance sur 22).

Cette règle des 5% – une chance sur 20 – est encore en vigueur aujourd’hui. C’est elle que l’on utilise dans tous les modèles économétrique, en regardant les « étoiles » associées à chaque des variables explicatives (3 étoiles si la p-value est inférieure à 0.1%, 1 étoile si elle est inférieure à 5%, et rien au-delà de 10%). Si la -value est une mesure continue de la distance entre l’hypothèse que l’on cherche à tester et les données, ces étoiles ont instauré des seuils, malheureusement supposé infranchissables. Avec cette méthode, comme le notait Guelman (2015), «it seems impressive to see multiple independent findings that are statistically significant but with enough effort it is possible to find statistical significance anywhere. »

Pour illustrer la difficulté de la prise de décision, Power (2014) prend un exemple tiré de la survie à des séances de chimiothérapie, pour guérir un cancer.

Figure 5: chimiothérapie, données Williams et al. (1987).

Dans un “research hospital” le taux de survie est de 86%, alors qu’il est de 81% dans un “non-research hospital”. Un test du chi-deux, mesurant l’indépendance entre le choix de l’hôpital et la survie, donne une statistique de test de 0.5635, correspondant à une p-value de 45.29%. Autrement dit, on devrait être indifférent à l’endroit où on va se faire traiter, puisque l’hôpital est non significativement corrélé à la survie. Mais lorsque la vie est en jeu, est-on prêt à dire qu’une différence de 5 points (86% vs. 81%) n’est pas “statistiquement significative” ? Même avec une p-value proche de 50% !

Le problème des tests multiples

Si cette méthode pour juger de la significativité d’une variable dans un modèle prédictif a eu beaucoup de succès, il convient d’insister sur ses limites quand on se retrouver face à des données massives. Si on cherche à expliquer une variable (comme la sinistralité) par 100 variables, possiblement indépendantes de notre variable d’intérêt, l’analyse précédant devrait nous pousser à retenir 5 variables comme « significativement » corrélées avec la sinistralité.

Figure 6 : Exemple de séries significativement (très) corrélées, données Vigen (2015)

Le problème des tests multiples est d’autant plus important en imagerie médicale. Sur la Figure 7, 100 échantillons Gaussiens sont simulés, de moyenne nulle, et on teste la nullité des moyennes (seules les p-values sont présentées). Imaginons des tests sur des pixels d’une image d’IRM par exemple. Sur une image obtenue sur un individu sain, si des tests sont effectués, pixel par pixel, 10% des points présenteront une p-value inférieure à 10%, 5% une p-value inférieure à 5%, etc.

Figure 7 : 100 tests de nullité de moyenne, sur 100 échantillons simulés (de moyenne nulle).

En l’occurrence, la vision Bayésienne de la prise de décision (décrite dans Greenland & Poole (2013)) peut s’avérer plus juste, en plus d’offrir une interprétation plus claire.

Les alternatives pour juger de la significativité

Les p-values sont un outil important, central, dans la construction de modèles prédictifs. Et la pratique des tests est souvent un exercice intéressant : si la p-value est petite, on se sent conforté dans l’idée que la variable est significative. Mais si elle est grande (disons excède 10%), on ose parfois se dire qu’on n’a pas eu de chance… Dans certains cas, on peut même s’offrir la chance de faire le test sur un autre échantillon et – avec un peu de chance – la p-value sera plus faible.

Comme le note Briggs (2013), les p-value “encourage magical thinking […] they focus attention on the unobservable“. Même avec un regard critique, l’utilisation des p-value est dangereuse. Et l’utiliser de manière automatique, dans un algorithme d’apprentissage, l’est encore plus. Sans bon sens, on verra des variables extrinsèques exotiques utilisées dans les modèles prédictifs (pour reprendre la terminologie de Cass & Shell (1983). A quand un assureur qui utiliserait la pointure de pieds, les résultats au brevet des collèges ou la couleur de la boite à lettres dans son tarif d’assurance auto ?

Références

Arbuthnot, J. (1710). An argument for Divine Providence, taken from the constant regularity observed in the births of both sexes. Royal Society’s Philosophical Transactions.

Berkson, J. (1942). Tests of significance considered as evidence, Hoover, K.D. & Siegle.

Briggs, W. (2013) Everything Wrong With P-Values Under One Roof. http://wmbriggs.com/…

Cass, D. & Shell, K. (1983). Do Sunspots Matter?. Journal of Political Economy 91 (21): 193–228.

Fisher, R. (1925). Statistical Methods for Research Workers. Oliver & Boyd.

Greenland S. & Poole C. (2013). Living with P values: resurrecting a Bayesian perspective on frequentist statistics. Epidemiology, 24: 62-8.

Guelman, A. (2015). P-values and statistical practice. andrewgelman.com/…

Hall, P. and Selinger, B. (1986). Statistical significance: balancing evidence against doubt. Australian & New Zealand Journal of Statistics, 28.

Powers, P.R. (2014). Acts of God and Man Ruminations on Risk and Insurance. ColumbiaUniversity Press.

Silver, N. (2012) The Signal and the Noise: Why So Many Predictions Fail — but Some Don’t, Penguin Press

Vigen, T (2015). Spurious Correlations http://tylervigen.com/…

Calcul(s) d’information de Fisher

La semaine passée, on avait fait quelques calculs pour obtenir l’information de Fisher pour des lois classiques. Je voulais juste remettre au propre les calculs pour les lois à plusieurs paramètres. Pour la loi Gamma,

la log-vraisemblance s’écrit

de telle sorte que

Ici, même pas besoin de prendre une espérance car la Hessienne est constante

Continue reading Calcul(s) d’information de Fisher

Box plot, Fisher’s style

In a recent issue of Significance, I discovered an interesting – and amuzing – figure, about some box & beard plot, in Dr Fisher’s casebook: Beard the statistician in his den.

In French, the box plot (introduced by John Tukey, not George Box, as discussed in a previous post) is popular under the name boîte à moustaches (box with a mustache, for a simple translation).

> set.seed(2)
> x=rnorm(500)
> boxplot(x,horizontal=TRUE,axes=FALSE)
> axis(1)

I was wondering if it was possible to reproduce that Fisher’s style of box plot, with a beard (Ronald Fisher was famous for his beard). Technically, it is not complicated,

> boxplot(x,horizontal=TRUE,xlim=c(-1,1.3),axes=FALSE)
> axis(1)
> Q=quantile(x,c(.25,.75))
> y=cut(x[(x>=Q[1])&(x<=Q[2])],seq(Q[1],Q[2],length=11))
> tb=table(y)
> u=seq(Q[1],Q[2],length=11)
> umid=(u[1:10]+u[2:11])/2
> for(i in 1:10) segments(umid[i],1-.2,umid[i],1-.2-tb[i]/20,lwd=3)

Yes, I kept the mustaches here, but ploting a simple box, one could easily draw a goatee, or a chin strip

> rect(Q[1],1-.2,Q[2],1+.2)
> segments(median(x),1-.2,median(x),1+.2,lwd=2)

The only problem, is that between the first and the third quartiles, the distribution is much more ‘uniform‘ than it might look on the figure in Significance. Unless, perhaps, if add more bars on the left of the first quartile, and on the right of the third one,

> du=diff(umid)[1]
> y=cut(x,seq(Q[1]-du*8,Q[2]+du*8,length=11+16))
> tb=table(y)
> u=seq(Q[1]-du*8,Q[2]+du*8,length=11+16)
> umid=(u[1:26]+u[2:27])/2
> for(i in 1:8) segments(umid[i],1,umid[i],1-.2-tb[i]/20,lwd=3)
> for(i in 19:26) segments(umid[i],1,umid[i],1-.2-tb[i]/20,lwd=3)

Of course, we can try to add a smile on that face, below the mustache,

> vu=seq(-2.5,2.5,by=.02)
> vv=dnorm(vu)
> lines(vu,1-vv*4,col="red",lty=2)
> d=density(x,bw=.1)
> lines(d$x,1-d$y*4,col="red")

So far, I am not convinced. I mean, it should be possible to add something on the box plot,  to add more information. But I still wonder what… maybe using colors?

De la significativité (statistique), suite

Suite à mon rapide billet sur l’article donnant les conclusions d’une étude des effets de la e-cigarette (et la comparaison avec le patch pour arrêter de fumer), j’avais entendu cette phrase “ça veut dire quoi ‘statistiquement non significatif’ ?“. Comme c’est une très bonne question, qui plus est relativement complexe, j’ai pris un peu de temps pour y répondre… Il ne s’agit pas de complexité mathématique, en s’entend car tous les cours de statistique abordent cette notion, mais plus d’une notion complexe à vulgariser. On va revenir un peu sur l’histoire des tests de significativité et sur les travaux de Ronald Fisher (pour comprendre l’origine des 5% utilisé comme seuil lors de l’utilisation de la p-value), mais il est possible d’aller voir directement les exemples.

  • Un peu d’histoire

Si on veut commencer par la préhistoire, les travaux les plus anciens qui posent la question de la significativité sont ceux de John Arbuthnot, en 1710. En prenant 82 années de données de natalité, à Londres, il avait observé que chacune des années, il y avait eu plus de naissances de garçons que de filles (évoqué dans Hacking (1975)).

Il s’était demandé, avec la terminologie d’il y a plus de 300 ans, si cette différence pouvait être jugée comme ‘statistiquement significative. Je pense que l’on peut utiliser le terme ‘statistiquement‘ car John Arbuthnot pose la question en termes probabilistes, très clairement,

Un siècle plus tard, Pierre-Simon de Laplace a présenté ce que l’on peut interpréter comme un ‘test de significativité‘, là encore avec notre terminologie actuelle. Il avait en effet noté, en prenant des mesures sur des baromètres que les observations à 9 heures du matin et 4 heures de l’après midi était différentes. Significativement différentes. Et là encore posé la question en terme probabiliste (que l’on verra formalisé bien plus tard), en se demandant s’il est ‘extrêmement probable‘ qu’il y ait une différence entre les deux mesures. Il avait alors fait un test de comparaison des valeurs moyennes, et noté que la différence excédait plusieurs écart-types, ce qu’il jugeait alors significativement important. A la fin du XIXème siècle, Francis Edgeworth reprendra cette idée (sur la comparaison de moyennes) lorsqu’il écrira en 1885,

où on retrouve les premières idées qui serviront de base aux travaux ultérieurs. Car pour bien comprendre l’origine de la notion de significativité en statistique, il faut relire les travaux de Ronald Fisher des années 20, mais aussi de Student (William Gosset) et Karl Pearson. Ronald Fisher s’interroge, dans son papier de 1921, sur la valeur d’un coefficient de corrélation (c’est dans ce papier qu’il introduit formelle sa transformation avec une tangente hyperbolique z = {1 \over 2}\ln\left({1+r \over 1-r}\right) = \text{arctanh}(r) pour faire un test de nullité d’un coefficient de corrélation), où il note

Aussi, dans le papier, une valeur (ici une corrélation) peut être jugée comme ‘significative‘ (c’est à dire ‘significativement non nulle‘). Encore une fois, l’idée n’est pas nouvelle, Edgeworth aurait dit ‘non accidentelle‘ mais Fisher pose ici la terminologie, en parlant de ‘significativité‘. L’année suivante, en 1922, il parle de l’idée de tester la significativité, dans le cadre d’un modele de régression (on est toujours dans le contexte de la corrélation) ou en 1924, on retrouve l’idée de tester la significativité, Mais c’est surtout dans son livre Statistical Methods for Research Workers que l’on va trouver une formalisation claire (et très illustrée) de cette idée, et surtout la connecter à la notion de p-value (on ne peut pas parler de significatité sans parler de p-value). Pour l’anecdote, dans l’ouvrage de Ronald Fisher, on y parle de p-value (oui, avec une majuscule). Cela dit, la paternité de la p-value est à attribuer davantage à Pearson, par exemple dans son papier datant de 1900, ou encore, sur un autre exemple avec en plus l’interprétation usuelle en terme de significativité (sans toutefois utiliser ce terme), puisque Pearson dit qu’on a 1 chance sur 8 d’avoir une valeur aussi improbable. Plus loin, il dit aussi qu’arriver avec 1 chance sur 70 peut être qualifier de très improbable (un long débat suivra sur le seuil à partir duquel, effectivement, un événement peut être qualifié d’improbable) Mais revenons au livre de Ronald Fisher, et plus spécifiquement au chapitre 4, qui propose précisément d’étudier la ‘significativité de la différence entre deux moyennes

On y retrouve le test de Pearson basé sur la statistique du chi-deux, comme test d’ajustement (goodness of fit) et d’indépendance. et la comparaison de la p-value avec (ce qui va correspondre aux mythiques) 5%. Il parle plus loin d’une ‘ligne conventionnelle‘ que l’on tracera à 5%, Il faut d’ailleurs un lien entre une p-value à 5% (correspondant à un événement qui survient avec 1 chance sur 20) et le fait de s’éloigner de 2 écart-types de la moyenne pour une loi normale (qui survient avec 1 chances sur 22) Tout au long du livre, on retrouve évoqué le chiffre mythique 5% qui restera à la postérité. Le plus intéressant pour apprendre comme juger la significativité d’une valeur, on peut regarder les exemples, comme l’exemple 11 l’exemple 12, l’exemple 27, ou encore l’exemple 32 (pour en prendre juste quelques uns) (on retrouve l’utilisation de 2 écart-types pour une loi normale). La discussion amorcée dans ce livre sera poursuivie dans son papier de 1926, qui commence avec Le papier est passionnant, car il pose les bases de ce que font tous les statisticiens depuis 80 ans (il n’est pas surprenant de voir cet article reproduit dans le second volume des Breakthroughs in Statistics), où le seuil de 5% est encore utilisé. Il ne cesse tout au long de l’article de justifier cette règle de ‘un sur vingt‘, tout en autorisant d’autres valeurs, comme 2%, ou 10% (Fisher semble moins orthodoxe que certains utilisateurs de la notion de significativité) Il faudra ensuite attendre les travaux de Jerzy Neyman et Egon Pearson. qui vont axiomatiser et formaliser les tests statistiques, mais ça dépasse (de loin) l’objectif de mon billet d’aujourd’hui. Donc on va repousser à plus tard la présentation de leurs apports. Pour illustrer tout ce qu’on vient de voir, considérons deux exemples simples.

  • Petit exemple, tester l’indépendance

Comme Ronald Fisher, commençons par ce qui me semble le plus simple car on n’a pas besoin de parler modèle probabiliste ici. On se demande si deux variables sont indépendantes, ou pas. Deux cas peuvent se présenter:

  • on a des variables continues, et on veut voir si elles sont non-corrélées (oui, pour commencer on va simplifier, et dire que si elles sont non-corrélées… ça va se traduire par une notion que l’on peut assimiler – rapidement – à une notion d’indépendance. Promis, d’ici 6 mois il y aura des dizaines de billets sur le sujet sur le blog puisque je redonne mon cours sur la dépendance cet hiver). Je laisse les plus motivés jouer avec les données qui traînent sur le blog pour voir si la taille et le poids sont corrélées, ou utiliser un exemple que j’avais l’habitude de présenter dans mes cours de modèles de régression, où on se demande si – chez les femmes – le tour de poitrine et le tour de taille sont corrélés.
  • on a des variables catégorielles, par exemple la couleur des yeux, et la couleur des cheveux. Et on se demande si les deux sont couleurs sont indépendantes, ou pas.

Regardons ce dernier exemple, avec un petit jeu de 25 observations.

> X=c(rep("Brun",8+5),rep("Blond",9+3))
> Y=c(rep("Noisette",8),rep("Bleu",5),rep("Noisette",3),rep("Bleu",9))

On peut résumer l’information dans un tableau de contingence,

> T=as.matrix(table(X,Y))
> T
       Y
X       Bleu Noisette
  Blond    9        3
  Brun     5        8

On va essayer de comparer ce tableau avec ce qu’on devait avoir si les deux variables étaient indépendantes (car c’est ce que l’on souhaite tester : on veut savoir si le tableau que l’on a est significativement différent de ce qu’on aurait avec des variables indépendantes). Il faut ici ressortir la définition de l’indépendance entre deux variables aléatoires, \mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B) = \mathbb{P}(X\in\mathcal A) \times\mathbb{P}(Y\in\mathcal B) pour tout \mathcal A et \mathcal B. Par exemple \mathcal A peut être avoir des cheveux blonds, et \mathcal B avoir des yeux bleus. Pour rappel, on avait en tout 25 personnes ayant ces deux caractéristiques en même temps, donc empiriquement, on avait observé une probabilité de 9/25 soit 36% (c’est le terme de gauche). Pour le terme de droite, on notera que \mathbb{P}(X\in\mathcal A), c’est la probabilité d’avoir les cheveux blonds. Empiriquement, on en avait 9+3 sur 25, soit 48%. Et \mathbb{P}(Y\in\mathcal B) c’est la probabilité d’avoir les yeux bleus. Or on avait 9+5 personnes avec ces caractéristiques, soit 56%. Si les deux variables étaient indépendantes, avec la formule précédente, la probabilité d’avoir à la fois les cheveux blonds et les yeux bleus serait 48% x 56%, soit 26.88%. A comparer aux 36% observés. En fait, au lieu de comparer des pourcentages, on peut aussi comparer des nombres: comme il y a 25 personnes en tout, 26.88% signifie que dans notre groupe, 6.72 personnes devraient avoir les cheveux blonds et les yeux bleus. On peut s’amuser à faire ces calculs pour tous les \mathcal A et \mathcal B.

> Tind=chisq.test(T)$expected
> Tind
       Y
X       Bleu Noisette
  Blond 6.72     5.28
  Brun  7.28     5.72

Visuellement, on peut mettre ça en forme,

(sur le graphique de droite, les traits noirs correspondent aux nombres sous hypothèse d’indépendance, et en couleur, les valeurs observées). Bon, maintenant, il faut réfléchir un peu… On se demande si nos deux tableaux sont proches ou pas. Ou, pour utiliser notre mot de la journée, significativement proches. Alors là, ça va devenir technique le temps d’un rapide paragraphe. Dans chaque case du tableau, on compte des gens. Sur nos n personnes, on peut s’attendre à en observer un nombre aléatoire, noté N_{\mathcal{A},\mathcal{B}}, qui suit une loi binomiale N_{\mathcal{A},\mathcal{B}}\sim\mathcal{B}(n,\mathbb{P}(X\in\mathcal A\text{ et }Y\in\mathcal B)) (c’est la définition de la loi binomiale). On notera que quand n est grand, la loi normale est alors proche de la loi Gaussienne. Dans chaque case, on a des lois normales (non indépendantes, car les nombres sont contraints, par exemple au total, il faut n personnes), et en prenant les carrés des lois normales (centrées et réduites), on va obtenir une loi du chi-deux (je renvoie vers un précédant billet pour les aspects techniques). Prendre le carré est intéressant, et naturel, car ça fait penser à une distance Euclidienne (usuelle). On va alors poser
Q=\sum_{\mathcal{A},\mathcal{B}} \frac{[n_{\mathcal{A},\mathcal{B}}-n_{\mathcal{A},\mathcal{B}}^\perp]^2}{n_{\mathcal{A},\mathcal{B}}^\perp}
n_{\mathcal{A},\mathcal{B}} est le nombre de personnes observées, de caractéristiques (jointes) \mathcal A et \mathcal B. et n^\perp_{\mathcal{A},\mathcal{B}} est ce qu’on aurait si les variables étaient indépendantes (calculé auparavant), soit

n^\perp_{\mathcal{A},\mathcal{B}}=n\times \frac{n_{\mathcal{A}}}{n}\times \frac{n_{\mathcal{B}}}{n}

avec les notations naturelles. Ici,

> sum((T-Tind)^2/Tind)

[1] 3.380994

Bon… on en fait quoi maintenant de ce 3.38 ? C’est là qu’on utilise le petit résultat évoqué auparavant, qui nous dit que, si n est grand, et que les variables sont effectivement indépendantes, alors Q va suivre une loi du chi-deux (en l’occurrence à 1 degré de liberté, pour tenir compte des diverses contraintes). On va alors se demander s’il est possible, ou vraisemblable, d’avoir 3.38 avec une loi du chi-deux.

> pchisq(3.380994,df=1)
[1] 0.9340477
> 1-pchisq(3.380994,df=1)
[1] 0.06595227

Aussi, une loi du chi-deux à un degré de liberté a 93.4% chances d’être inférieure à la valeur obtenue, ou encore

\mathbb{P}(Q>3.38) \approx 6.59\% avec Q\sim \chi^2(1)

C’est cette grandeur que l’on appelle p-value et que l’on va essayer de comparer au seuil désormais mythique de 5%,

  • si cette probabilité est inférieure à 5%, on va rejeter l’hypothèse d’indépendance
  • si cette probabilité est supérieure à 5%, on va retenir l’hypothèse d’indépendance (pour faire simple)
> chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 3.381, df = 1, p-value = 0.06595

Les plus attentifs devraient me faire remarquer que c’est vaseux mon histoire : tout ce que je raconte est valide si n est grand, et dire que 25 c’est grand… c’est limite ! Tout d’abord, je dirais que c’est ce qui est utilisé dans un paquet d’études en médecine (peut-être un peu plus, mais guère plus). Ensuite je dirais que, malheureusement, la plupart des résultats que l’on voit dans les cours de statistiques sont de ce genre, à savoir des résultats asymptotiques, valides seulement si n est grand (c’est pareil en économie, ou des tests asymptotiques sont parfois évoqués avec 25 années d’observations). Mais heureusement, en s’accrochant un peu, on peut s’affranchir de cette hypothèse. A condition de faire un peu de bootstrap. En fait, l’idée est très très simple. On peut utiliser des générateurs de nombres aléatoires pour générer des tableaux de contingences pour lesquels les variables sont indépendantes. Car générer des variables aléatoires indépendantes, c’est ce qui est le plus simple. En fait, ici, on va se contenter de mélanger les variables,

> s=1
> set.seed(s)
> Xs=sample(X)
> Ys=sample(Y)

En faisant une permutation, je m’assure d’avoir autant de personnes dont la couleur des cheveux est \mathcal A que dans la population initiale (je ne fais que permuter) et de personnes dont la couleur des yeux est \mathcal B. Par contre, comme je mélange indépendamment, je suis certain que les variables sont indépendantes,

> Ts=as.matrix(table(Xs,Ys))
> Ts
       Ys
Xs      Bleu Noisette
  Blond    7        5
  Brun     7        6

Générer un tableau de contingence, c’est juste pour illustrer, car l’idée est de générer mille (ou plus encore) tableaux de contingence obtenus en simulant des variables indépendantes, et de voir ce que cette statistique du chi-deux donne

> Q=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ Xs=sample(X)
+ Ys=sample(Y)
+ Ts=as.matrix(table(Xs,Ys))
+ Q[s]=sum((Ts-Tind)^2/Tind)}

Si on regarde les valeurs obtenues

>  table(Q) 
Q
0.0509906759906761  0.337162837162837   1.06560106560107   1.92411754911755 
               298                262                200                133 
  3.38099400599401   4.81185481185481    6.9971694971695   9.00037462537462 
                60                 33                  9                  3 
  11.9141275391275    14.489676989677 
                 1                  1

autrement dit, dans 29.8% des scénarios, on avait une distance (un statistique) valant 0.05, dans 26.2% des scénarios, on avait 0.33, et dans 20% des scénarios, on avait 1.06. Etc. Si on cumule on obtient

>  cumsum(table(Q))/10
0.0509906759906761  0.33162837162837   1.06560106560107    1.92411754911755 
3.38099400599401    4.81185481185481    6.9971694971695     9.00037462537462 
11.9141275391275     14.489676989677 

Pour rappel, on avait eu

>  sum((T-Tind)^2/Tind) 
[1] 3.380994

qui peut effectivement être obtenu avec des variables indépendantes. Mais seulement dans 6% des scénarios. En fait, dans 95.3% des scénarios, on a eu une valeur inférieure ou égale à celle obtenue sur notre échantillon. Et dans 4.7% des scénarios, on a strictement dépassé cette valeur. On pourrait parler de p-value obtenue par bootstrap, ou rééchantillonage. On voit que l’on flirte ici avec le seuil des 5%… Visusellement, on a

> plot(c(0,as.numeric(names(table(Q)))),c(0,
+ cumsum(table(Q))/10),type="s",xlab="",ylab="")
> abline(v=sum((T-Tind)^2/Tind),lty=2)
> lines(seq(0,15,by=.01),pchisq(seq(0,15,by=.01),1)*100,col="red")

Ici, on est bien embêté pour conclure quoi que ce soit… Ce qui est amusant, c’est qu’avec un échantillon plus grand (trois fois plus grand par exemple), et les mêmes pourcentages, il n’y a plus aucune ambiguïté,

> T=T*3
> T
       Y
X       Bleu Noisette
  Blond   27        9
  Brun    15       24
>  Tind=chisq.test(T)$expected
> Tind Y X Bleu Noisette Blond 20.16 15.84 Brun 21.84 17.16

que l’on peut encore visualiser ci-dessous,

Ici, on obtient une p-value de l’ordre de 0.1%,

>  chisq.test(T,correct=FALSE)

	Pearson's Chi-squared test

data:  T
X-squared = 10.143, df = 1, p-value = 0.001449

On peut aussi faire du bootstrap si on trouve que 75 observations, c’est trop faible, mais la conclusion est la même,

  • Petit exemple, sur la comparaison de moyenne de deux groupes

Cette fois, on a deux populations, et on se demande si les deux groupes sont ‘comparables‘. On peut penser à comparer le nombre de personnes guéries au sein de deux populations, une qui a eu un placebo, et l’autre un médicament. Cela dit, cet exemple simple est en fait compliqué car ‘guérie’ n’est peut-être pas une notion claire. Prenons plutôt des quantités que l’on peut mesurer. Disons la taille des gens, et demandons nous si dans un groupe, les gens sont significativement plus grands que dans l’autre. Je vais tricher un peu… je vais prendre deux populations que je sais être différentes, a priori. Un groupes d’hommes, et un groupes de femmes. On va noter \{x_1,\cdots,x_{n_X}\} et \{y_1,\cdots,y_{n_Y}\} ces deux échantillons, car il va falloir formaliser un minimum.

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> Davis[12,c(2,3)]=Davis[12,c(3,2)] 
> Davis=Davis[order(Davis$height),]
> attach(Davis)
> set.seed(1)
> X=sample(height[sex=="F"],size=6)
> Y=sample(height[sex=="M"],size=5)

On a ici 5 personnes dans le premier groupe et 6 personnes dans le second. Et on se demande si les tailles sont significativement différentes entre les deux groupes,

> (mx=mean(X))
[1] 165.5
> (sx=sd(X))
[1] 4.679744
> (my=mean(Y))
[1] 178
> (sy=sd(Y))
[1] 7.968689

On veut savoir si l’écart entre 165,5 cm et 178 cm est significativement différent de zéro (on pourrait aussi se demander s’il est significativement positif, mais ça serait une autre histoire). On va faire une hypothèse ici : on va supposer que nos échantillons sont obtenus comme des tirages de variables aléatoires, indépendantes, x_i=X_i(\omega) et y_j=Y_j(\omega). Le modèle qu’on va supposer est que la taille des gens suit une loi normale, conditionnellement à leur sexe,

\left\{\begin{array}{l}X_i\sim\mathcal{N}(\mu_X,\sigma_X^2) \\Y_j\sim\mathcal{N}(\mu_Y,\sigma_Y^2)\end{array}\right.

Si on veut visualiser nos données, on peu utiliser

> u=seq(150,200,by=.25)
> v1=dnorm(u,mx,sx)
> v2=dnorm(u,my,sy)
> plot(u,v1,col="red",type="l",ylim=c(0,.09),axes=FALSE)
> axis(1)
> lines(u,v2,col="blue")
> points(X,rep(.075,length(X)),col="red")
> points(Y,rep(.075,length(Y)),col="blue")
> abline(v=mx,lty=2,col="red")
> abline(v=my,lty=2,col="blue")

 

Le plus simple aurait été de supposer que les deux lois ont la même variance, mais quand on regarde les densités sur le graphique ci-dessus, il faut rester crédible. Bref, on va utiliser le test de Welch. L’idée est d’utiliser un petit exercice de probabilité. Soient \{X_1,\cdots,X_{n_X}\} et \{Y_1,\cdots,Y_{n_Y}\} des variables aléatoires, Gaussiennes (comme supposé juste avant de faire un dessin) et indépendantes. Posons alors
\overline{X}=\frac{1}{n_X}\sum_{i=1}^{n_X} X_i
et
s_X=\frac{1}{n_X-1}\sum_{i=1}^{n_X} [X_i-\overline{X}]^2
(avec les quantités similaires pour le second échantillon) alors si \mu_X=\mu_Y, et si on note
s = \sqrt{{s_X^2 \over n_X} + {s_Y^2 \over n_Y}}.
la variable aléatoire
T = {\overline{X} - \overline{Y} \over s[latex]<br /> va suivre une loi de Student (connue depuis les travaux de <a href="http://freakonometrics.hypotheses.org/1968">William Gosset</a>), avec</p> [latex display="true"]\frac{(s_X^2/n_X + s_Y^2/n_Y)^2}{(s_X^2/n_X)^2/(n_X-1) + (s_Y^2/n_Y)^2/(n_Y-1)}

degrés de liberté. On notera qu'il ne s'agit pas ici d'un résultat asymptotique. Par contre, il repose intégralement sur une hypothèse - forte - de normalité des tailles des individus. L'idée est alors la même qu'auparavant : on va calculer cette statistique sur notre échantillon,

> nx=length(X)
> ny=length(Y)
> s=sqrt(sx^2/nx+sy^2/ny)
> T=(mx-my)/s
> T
[1] -3.091371

et se demander s'il est vraisemblable d'avoir une telle valeur avec une loi de Student. Le nombre de degrés de liberté est ici

> d=(sx^2/nx+sy^2/ny)^2/((sx^2/nx)^2/(nx-1)+(sy^2/ny)^2/(ny-1))
> d
[1] 6.218682

On pourrait alors calculer
\mathbb{P}(T>t)T\sim\mathcal{S}td(d)

mais ça serait un test pour savoir si la différence est significativement positive, ce qui n'est pas ce qu'on avait dit que l'on cherchait : on veut ici significativement non-nulle.

\mathbb{P}(\vert T\vert> \vert t\vert)T\sim\mathcal{S}td(d)

Si on veut visualiser un peu tout ça, on peut utiliser

> u=seq(-4,4,by=.01)
> v=dt(u,df=d)
> plot(u,v,type="l")
> abline(v=T,col="red")
> u=seq(-4,T,length=100)
> polygon(c(u,rev(u)),c(dt(u,df=d),rep(0,100)),col="red",border=NA)

Pour calculer ces probabilités, on utilise

> pt(T,d)+(1-pt(-T,d))
[1] 0.02038577

et comme la loi de Student est symétrique par rapport à l’origine,

> pt(T,d)*2           
[1] 0.02038577

Aussi, la p-value est de l'ordre de 2%, ce qui est suffisamment petit pour dire que, oui, nos moyennes sont différentes. C'est d'ailleurs ce que renvoie la commande

> t.test(X,Y)

	Welch Two Sample t-test

data:  X and Y
t = -3.0914, df = 6.219, p-value = 0.02039
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -22.31039  -2.68961
sample estimates:
mean of x mean of y 
    165.5     178.0

On avait vu tout à l'heure que l'idée de la valeur de 5% était liée au scénario qui se répète avec 1 chance sur 20 : si on pouvait générer d'autres échantillons, indépendants, alors 1 fois sur 20 on devrait atteindre une statistique de test de l'ordre de celle observée. Pour générer des échantillons vraiment indépendants, on avait utilisé du bootstrap, en rééchantillonnant dans notre population de test. Ici, on a bien plus que nos 11 observations, on en a en réalité 200. On peut se demande ce qui se serait passé si on avait constitué d'autres échantillons,

> Diff=rep(NA,1000)
> for(s in 1:1000){
+ set.seed(s)
+ X=sample(height[sex=="F"],size=6)
+ Y=sample(height[sex=="M"],size=5) 
+ Diff[s]=(mean(X)-mean(Y))/(sqrt(sd(X)^2/nx+sd(Y)^2/ny))}
> hist(Diff,col="light blue",border="white")

Sur aucun échantillon, la moyenne du premier groupe n'a dépassé la moyenne du second. Alors qu'on aurait pu penser en avoir une vingtaine (2% des mille scénarios générés). Ceci vient peut-être du fait que l'hypothèse de normalité pour les tailles peut être trop restrictive...

  • Pour aller plus loin

Le billet aura été un peu long, mais on aura eu le temps de voir en détails d'où vient cette notion de significativité statistique, et ses liens avec la p-value. On aura même eu le temps de voir sur deux exemples l'importance des hypothèses sous-jacentes, en particulier si le test est asymptotique, ou s'il repose sur des hypothèses fortes de lois. Mais ce n'était qu'une courte mise en bouche. Pour aller plus loin, je suggère la lecture de Berkson, J. (1942). Tests of significance considered as evidence, Hoover, K.D. & Siegle, M.V. (2008) Sound and fury: McCloskey and significance testing in economics, Liberman, M. (2013). "Significance", in 1885 and today et le passionnant Hall, P. and Selinger, B. (1986). Statistical significance: balancing evidence against doubt. Et sur les p-value, Goodman, S. (2002) A Dirty Dozen:Twelve P-Value Misconceptions. Je devrais aussi citer Ziliak, S. & McCloskey, D. (2008). The cult of statistical significance, mais je ne l'ai pas encore lu. Par contre, Xi'an en parlait sur son blog. Et si j'étais courageux, je tenterais une lecture bayésienne de ce billet.... promis, dès que je trouve un peu de temps !

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 application du test de Fisher, partie 2

Poursuivons un peu la discussion d’hier. En fait, si on poursuit, on arrive à regrouper des régions ensemble…. mais sans réelle cohérence car des régions proches au niveau électorale peuvent être éloignées géographiquement.

Si on veut que les régions regroupées soient proches, il faut une base avec des informations spatiales, comme savoir si deux départements se touchent, ou pas. Pour cela, on utilise la base suivante, qui contient

  • le numéro du département i
  • le numéro du département j
  • une variable indicatrice indiquant si les département se touche
  • la distance (en km) entre les centroïdes des départements
> france=read.table(
+ "http://freakonometrics.blog.free.fr/public/
data/departements-france.csv",
+ header=TRUE,sep=";")

On va nettoyer un peu, en enlevant les départements d’outre-mer,

> france0=france; i200=c(201,202)
> france0=france0[-c(which(france0$depi%in%i200),
+ which(france0$depj%in%i200)),]
> france0$ij=100*france0$depi+france0$depj

La base ressemble à ça,

> head(france0)
depi depj Cij Dij   ij
1   77   75   0  51 7775
2   78   75   0  39 7875
3   91   75   0  40 9175
4   92   75   1  10 9275
5   93   75   1  12 9375
6   94   75   1  13 9475

Plusieurs stratégies sont possible: regarde des départements voisins et voir si on peut les regrouper. Ou regarder des départements proches électoralement, et voir s’il sont proches géographiquement. Par exemple, si on ordonne les coefficients obtenus lors de la régression,

> reg=lm(X.VoixExpH~0+Dpt,weights=
+ Votants,data=sbelection)
> COEF=coefficients(summary(reg))
> dbcoef=data.frame(COEF[,1:3],Dpt=rownames(COEF))
> rownames(dbcoef)=rank(COEF[,1])
> COEFtrie=dbcoef[as.character(1:94),]
> COEFtrie[1:6,]
Estimate Std..Error   t.value    Dpt
1 35.69903  0.3246129 109.97414 DptD06
2 36.46840  0.3244806 112.39010 DptD67
3 36.67594  0.3896297  94.13025 DptD68
4 37.39283  0.3202112 116.77552 DptD83
5 39.92023  0.3921762 101.79156 DptD74
6 42.61811  0.6182921  68.92877 DptD10

On peut se demander si les quatre premiers départements pourraient être regroupés ensemble. Une des modalités sert de modalité de référence, et on regarde les trois autres par rapport à cette dernière

> sbelection$Dpt=relevel(sbelection$Dpt,"D67")
> reg0=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> COEF0=coefficients(summary(reg0))
> COEF0[c(1,7,67,82),]
Estimate Std. Error     t value   Pr(>|t|)
(Intercept) 36.4684029  0.3244806 112.3900971 0.00000000
DptD06      -0.7693736  0.4589784  -1.6762740 0.09369320
DptD68       0.2075400  0.5070493   0.4093093 0.68231515
DptD83       0.9244243  0.4558759   2.0277980 0.04258819

On fait ensuite une analyse de la variance,

> sbelection$DptGroupe=sbelection$Dpt
> sblction$DptGroupe[sblction$Dpt%in%c("D06","D68","D83")]="D67"
> model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
> model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F   Pr(>F)
1  36114 2281072960
2  36111 2280191044  3    881917 4.6556 0.002954 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

On rejette ici l’hypothèse que les quatre départements se comportent de la même manière. Mais si on tente de regrouper les trois premières, cette fois, on accepte l’hypothèse

> sbelection$DptGroupe=sbelection$Dpt
> sblction$DptGroupe[sblction$Dpt%in%c("D06","D68")]="D67"
> model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
> model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36113 2280476723
2  36111 2280191044  2    285679 2.2621 0.1041

Sauf que sur les trois départements, seulement deux se touchent,

> france0[france0$ij%in%c(6768,667,668),]
depi depj Cij Dij   ij
3510    6   67   0 533  667
3552   67   68   1  95 6768
3605    6   68   0 439  668

On peut essayer de parcourir les sorties de la régression, et regarder si on a des voisins. Pour chaque département, on regarde si le voisin est un voisin géographique,

> departement=as.numeric(substr(COEFtrie$Dpt,5,6))
> deptij=departement[1:(nrow(COEFtrie)-1)]*100+
> departement[(2):(nrow(COEFtrie))]
> voisins=france0[france0$ij%in%deptij,]
> voisinage=tapply(voisins$Cij,factor(voisins$ij),max)
> COEFtrie$voisin=NA
> v=voisinage[as.character(deptij)]
> COEFtrie$voisin[1:(nrow(COEFtrie)-1)]=
> v[is.na(v)==FALSE]
> COEFtrie[1:10,]
Estimate Std..Error   t.value    Dpt voisin
1  35.69903  0.3246129 109.97414 DptD06      0
2  36.46840  0.3244806 112.39010 DptD67      1
3  36.67594  0.3896297  94.13025 DptD68      0
4  37.39283  0.3202112 116.77552 DptD83      0
5  39.92023  0.3921762 101.79156 DptD74      0
6  42.61811  0.6182921  68.92877 DptD10      0
7  42.77423  0.4397101  97.27826 DptD01      0
8  43.53395  0.4448212  97.86843 DptD84      0
9  44.38759  0.3897553 113.88579 DptD85      0
10 44.65361  0.4577710  97.54575 DptD51      1

On note qu’il y a plusieurs voisins. On peut ensuite faire un test de Fisher. S’il est positif, on fusionne les classes. Pour la régression, mais aussi dans la base des distances et des voisinages,

>  if(COEFtrie$voisin[i]==1){
+  sbelection$Dpt[sbelection$Dpt==substr(
+ COEFtrie$Dpt[i+vsn],4,6)]=substr(COEFtrie$Dpt[i],4,6)
+  france0$depi[france0$depi==as.numeric(
+  substr(COEFtrie$Dpt[i+vsn],5,6))]=as.numeric(
+ substr(COEFtrie$Dpt[i],5,6))
+  france0$depi[france0$depj==as.numeric(
+ substr(COEFtrie$Dpt[i+vsn],5,6))]=as.numeric(
+ substr(COEFtrie$Dpt[i],5,6))

Dans la base des départements, il y aura plusieurs paires. Si A et B ont fusionné, il suffit que C touche A ou B pour que C soit voisin des régions fusionnées. C’est pour cela que l’on utilise

> voisinage=tapply(voisins$Cij,factor(voisins$ij),max)

pour savoir si les départements sont voisins. On peut alors répéter plusieurs fois l’opération.

Par exemple, on note que 28, 45 et 89 peuvent fusionner,

> sbelection$Dpt=relevel(sbelection$Dpt,"D45")
+ reg10=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
+ COEF10=coefficients(summary(reg10))
+ COEF10[c(1,28,88),]
Estimate Std. Error     t value  Pr(>|t|)
(Intercept) 45.9303788  0.4170114 110.1417762 0.0000000
DptD28       0.5643212  0.6594028   0.8558065 0.3921105
DptD89       0.9423212  0.7016026   1.3430982 0.1792486
+ sbelection$DptGroupe=sbelection$Dpt
+ sbelection$DptGroupe[sbelection$Dpt%in%
+ c("D28","D89")]="D45"
+ model0=lm(X.VoixExpH~DptGroupe,weights=Votants,data=sbelection)
+ model1=lm(X.VoixExpH~Dpt,weights=Votants,data=sbelection)
+ anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36113 2280314007
2  36111 2280191044  2    122963 0.9737 0.3777

Géographiquement, ce sont effectivement des voisins

Mais peu de départements voisins peuvent effectivement être fusionnés.

Une application du test de Fisher, partie 1

Un petit (?) billet pour revenir sur un point évoqué mardi dernier en cours, sur l’utilisation du test de Fisher (et des techniques d’analyse de variance). Je vais juste évoquer rapidement un vieux billet sur l’histoire du 5% utilisé lorsque l’on calcule la -value (ce qui complètera la réponse faite en classe mardi lorsque la question a été posée). Mardi dernier, j’avais évoqué rapidement le test d’une contrainte linéaire sur les paramètres. On se demande ici si les deux modèles suivants sont vraiment différents, entre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-2.png

et sa version contrainte

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-1.png

On cherche alors à tester l’hypothèse

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-3.png

contre l’hypothèse alternative

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-4.png

Pour ça, on peut utiliser ce qu’a proposé Ronald Fisher, à savoir d’utiliser la somme des carrés des résidus (qui correspondent – à un facteur multiplicatif près – à l’estimateur de la variance du bruit), et il a posé

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-9.png

Je parlais tout à l’heure de choix de modèle… on notera que l’on en est pas très loin car on peut réécrire cette expression en faisant apparaître le  des régressions, contraintes et noncontraintes.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-10.png

Fisher a montré que sous l’hypothèse H0,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fisher-11.png

i.e. la loi de la statistique de Fisher est précisément une loi de Fisher, Pour aller plus loin sur l’analyse de la variance (ANOVA, analysis of variance), je renvois àFaraway (2005).

J’avais dit oralement qu’une application intéressant du test de nullité de plusieurs coefficient dans une régression était son utilisation pour des variables catégorielles, et pour faire du regroupement de modalité.

Considérons un cas d’école. On dispose des résultats des votes dans plusieurs villes, en France, lors des dernières élections présidentielles (ces bases sont en ligne sur http://www.data.gouv.fr/), opposant François Hollande (avec un H) et Nicolas Sarkozy (avec un S),

> election=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/
205cf26f15a43974fdb7d3b25397613d.csv",
+ sep=";",header=TRUE,dec=",")
> head(election)
Code.du.dpt Inscrits Votants VoixH X.VoixExpH VoixS X.VoixExpS
1           1      592     506   195      41.05   280      58.95
2           1      215     175    82      52.23    75      47.77
3           1     8208    6447  3009      50.41  2960      49.59
4           1     1152     949   340      38.16   551      61.84
5           1      105      83    35      44.30    44      55.70
6           1     1702    1461   630      46.46   726      53.54

On va nettoyer un peu la base

> sbelection=election[election$Code.du.département
+ %in%as.character(1:95),]

car il y a des départements d’outre mers et la Corse, qui a une codification particulière. Bref, on retient les 95 départements métropolitain, moins le 20ème qui est la Corse. On va commencer par créer une variable de département qui soit un facteur qualitatif. Pour l’instant, le département est un numéro, entre 1 et 95 mais ce n’est aucunement un nombre !

> sbelection$Dpt=(paste("D",
> sbelection$Code.du.département,sep=""))
> I=as.numeric(as.character(sbelection$Code.du.département))<10
> sbelection$Dpt[I]=(paste("D0",
+ sbelection$Code.du.département,sep=""))
> sbelection$Dpt=factor(sbelection$Dpt)

On peut alors faire une régression, où on peut essayer de modéliser le taux de vote pour François Hollande (par exemple) en fonction de variables explicatives, et en particulier, on veut voir s’il n’y a pas un effet spatial. Bref, on va régresser sur le département. Et d’ailleurs que sur le département, inutile de mettre d’autre variable pour illustrer.

> reg=lm(X.VoixExpH~0+Dpt,weights=
+ Votants,data=sbelection)

Par contre, je fais une régression pondérée, car les villes (et villages) ne font pas du tout la même taille. C’est exactement la même idée que de faire une moyenne pondérée pour avoir le taux moyen dans une région (on rediscutera des régressions pondérées en fin de cours).

> COEF=coefficients(summary(reg))

La sortie ressemble à

> COEF
Estimate Std. Error   t value Pr(>|t|)
DptD01 42.77423  0.4397101  97.27826        0
DptD02 52.39423  0.4572016 114.59765        0
DptD03 56.89434  0.5468760 104.03518        0
DptD04 51.06544  0.7807748  65.40355        0
DptD05 50.70525  0.8403979  60.33481        0
DptD06 35.69903  0.3246129 109.97414        0
DptD07 53.44185  0.5564489  96.04089        0
DptD08 51.86319  0.6360901  81.53434        0
DptD09 64.67228  0.8013389  80.70528        0
DptD10 42.61811  0.6182921  68.92877        0

[…]

DptD89 46.87270  0.5642231  83.07476        0
DptD90 50.47561  0.9063735  55.68964        0
DptD91 53.41748  0.3172310 168.38670        0
DptD92 49.54305  0.2850597 173.79888        0
DptD93 65.29062  0.3338954 195.54217        0
DptD94 56.48855  0.3196643 176.71211        0
DptD95 53.89382  0.3344373 161.14778        0

On peut d’ailleurs faire une prédiction, par département (en prenant comme prédiction celle obtenue pour la première ville du département)

> chg=which(diff(as.numeric(as.character(
+ sbelection$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)

On peut ensuite essayer de faire un petite dessin pour visualiser cette prédiction,

> chg=which(diff(as.numeric(as.character(
+ sbelection$Code.du.département)))!=0)
> newelection=sbelection[1+c(0,chg),]
> PredDept=predict(reg,newdata=newelection)
> nomdpt=read.table(
+ "http://freakonometrics.blog.free.fr/
public/data/departements.txt",sep=";")
> nomdpt$nom=tolower(nomdpt$V2)
> dpt=as.character(c(paste("0",
+ 1:9,sep=""),c(10:19,21:95)))
> noms=as.character(nomdpt[nomdpt$V3%in%dpt,2])
> couleurliste=rev(rainbow(n=100))
> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur=couleurliste[trunc((PredDept-35)*3)]

C’est un peu long (et compliqué) car il faut identifier le bon département sur la carte R (qui est codée en lettres). D’ailleurs la fonction suivante ne marche pas

> match=match.map(francemap,dpt)
> francemap$names[which(is.na(match.map(
+ francemap,dpt[35]))==FALSE)]
[1] "Indre-et-Loire" "Indre"

car en l’occurrence, seule l’Indre aurait du être obtenue. On va faire la correspondance à la main,

> match=rep(NA,length(francemap$names))
> for(i in 1:length(francemap$names)){
+ if(sum(dpt==francemap$names[i])>0)
+ {match[i]=which(dpt==francemap$names[i])}
+ }

Une fois les départements identifiés, on colorie,

> color=couleur[match]
> map(database="france", fill=TRUE, col=color)

C’est joli, mais comme toujours c’est un peu compliqué de faire des cartes…

A titre d’information, les couleurs correspondent aux valeurs suivantes (avec en abscisse, le taux de vote obtenu par François Hollande)

Bref, on a beaucoup de modalité car il y a 95 département. Mais on devrait pouvoir regrouper. Si on ne sait pas trop par où commencer, prenons un département au hasard, disons le 35 (i.e. Rennes),

> k=34

On peut ensuite regarder les estimations des coefficients, ou mieux, les intervalles de confiances des différents estimateurs, et voir ceux qui chevauche celui de Rennes (on va prendre des intervalles de confiance à 70% histoire d’avoir des chances que les deux coefficients soient vraiment proches)

> alpha=.7
> CI.COEF=confint(reg,level=alpha)

On peut visualiser ces intervalles, avec une bande horizontale pour Rennes,

> plot(c(1:19,21:95),COEF[,1],ylim=c(35,
+ 60),pch=3,cex=.5)
> segments(c(1:19,21:95),CI.COEF[,1],
+ c(1:19,21:95),CI.COEF[,2])
> lk=which((CI.COEF[,1]<CI.COEF[k,2])&
+ (CI.COEF[,2]>CI.COEF[k,1]))
> segments(c(1:19,21:95)[lk],CI.COEF[lk,1],
+ c(1:19,21:95)[lk],CI.COEF[lk,2],col="red",lwd=2)

et en rouge, les intervalles qui la chevauchent

Ça y est, on a une liste potentielle de six départements à regrouper… On va prendre comme modalité de référence le département qui nous intéresse, et voir si les coefficients des autres départements sont (ou pas) significatifs.

> sbelection$Dpt=relevel(sbelection$Dpt,"D35")
> reg35=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)
> COEF35=coefficients(summary(reg35))
> COEF35[c(1,lk[lk<k]+1,lk[lk>k]),]
Estimate Std. Error     t value  Pr(>|t|)
(Intercept) 55.66226653  0.3275407 169.9399801 0.0000000
DptD11       0.59043789  0.6272106   0.9413710 0.3465210
DptD32       0.96845567  0.7855154   1.2328920 0.2176241
DptD36      -0.03007518  0.7414655  -0.0405618 0.9676455
DptD62       0.56238714  0.4247897   1.3239189 0.1855384
DptD75      -0.06226653  0.4098013  -0.1519432 0.8792326
DptD81      -0.11970183  0.6067040  -0.1972986 0.8435950

On a confirmation que ces six départements sont non significativement (individuellement) différents de celui de référence (et donc on accepte l’hypothèse de nullité du coefficient). Mais peut-on pour autant les regrouper? C’est là que l’on peut faire un test de Fisher. Pour cela on fait deux régressions. La première est celle avec un regroupement de départements en une modalité unique (les sept sont désormais ensemble)

> sbelection$DptGroupe=sbelection$Dpt
> sbelection$DptGroupe[sbelection$Dpt%
+ in%substr(rownames(COEF[lk,]),4,6)]=paste("D",k+1,sep="")
> model0=lm(X.VoixExpH~DptGroupe,
+ weights=Votants,data=sbelection)

et la seconde est celle que l’ont vient de faire (qui sera notre hypothèse alternative, i.e. il ne faut pas regrouper les départements ensemble car ils sont différents – ou au moins l’un est différent)

> model1=lm(X.VoixExpH~Dpt,
+ weights=Votants,data=sbelection)

On peut alors faire le test de Fisher (en continuant à corriger de l’effet taille de la population, et en faisant des sommes pondérées)

> SCR0=sum(sbelection$Votants*residuals(model0)^2)
> SCR1=sum(sbelection$Votants*residuals(model1)^2)
> dl0=model0$df.residual
> dl1=model1$df.residual
> (F=(SCR0-SCR1)/(dl0-dl1)/SCR1*dl1)
[1] 0.917397

On obtient la statistique de Fisher, et le plus simple est alors de calculer la -value

> 1-pf(F,dl0-dl1,dl1)
[1] 0.4809411

On a fait (sans le savoir) une analyse de la variance, qui peut se faire plus simplement (?) en utilisant la fonction suivante

> anova(model0,model1)
Analysis of Variance Table

Model 1: X.VoixExpH ~ DptGroupe
Model 2: X.VoixExpH ~ Dpt
Res.Df        RSS Df Sum of Sq      F Pr(>F)
1  36117 2280538612
2  36111 2280191044  6    347568 0.9174 0.4809

Et on accepte ici l’hypothèse que ces sept départements peuvent être regroupé en un seul.

Enfin, si statistiquement, il y a du sens à regrouper les modalités, on peut se demander si ce regroupement de régions géographique à du sens,

> library(maps)
> francemap<-map(database="france")
> dpt=noms
> couleur="red"
> match=match.map(francemap,dpt)
> color=couleur[match]
> map(dataLa première  estbase="france", fill=TRUE, col=color)

Les régions que l’on souhaite associer à Rennes sont fort éloignées, et les regrouper n’a peut-être pas de sens… (à suivre donc).

Fisher-Tippett theorem with an historical perspective

A couple of weeks ago, Rafael asked me if I had something on the history of extreme value theory. Since I will get back to fundamental results about extremes in my course, I promised I will write down a short post on all that issue.

To start from the beginning, in 1928, Ronald Fisher and Leonard Tippett formulated the three types of limiting distributions for the maximum term of a random sample (Fisher & Tippett (1928)). The problem was to characterize function http://freakonometrics.hypotheses.org/files/2015/12/ext-01.gif such that

http://freakonometrics.hypotheses.org/files/2015/12/ext-2.gif

where http://freakonometrics.hypotheses.org/files/2015/12/ext-3.gif where http://freakonometrics.hypotheses.org/files/2015/12/ext-4.gif‘s are i.i.d. with cumulative distribution function http://freakonometrics.hypotheses.org/files/2015/12/ext-5.gif. They had supporting arguments, but no (rigorous) proof. Nevertheless, the obtained that the only possible types for G were

http://freakonometrics.hypotheses.org/files/2015/12/ext-6.gif

i.e. Fréchet type (Pareto-type tails), or

http://freakonometrics.hypotheses.org/files/2015/12/ext-7.gif

i.e. Weibull type (bounded distribution type), or

http://freakonometrics.hypotheses.org/files/2015/12/ext-8.gif

i.e. Gumbel type (exponential-type tails). Emil Gumbel has been intensively using the so-called Gumbel distribution on river flows, since (as he explained in 1958), “it seems that the rivers know the theory. It only remains to convince the engineers of the validity of this analysis“.
Independently of that work (published in 1928), Maurice Fréchet considered in 1927 (in Sur la loi de probabilité de l’écart maximum) possible limits of

http://freakonometrics.hypotheses.org/files/2015/12/ext-9.gif

and obtained only http://freakonometrics.hypotheses.org/files/2015/12/ext-10.gif as possible limit. Richard von Mises gave in 1936 sufficient, but not necessary conditions for their (max) domain of attraction, i.e. characterization of function http://freakonometrics.hypotheses.org/files/2015/12/ext-11.gif such that the maxima converges to some specific function http://freakonometrics.hypotheses.org/files/2015/12/ext-01.gif (von Mises (1936)). E.g. he noticed that a sufficient condition on http://freakonometrics.hypotheses.org/files/2015/12/ext-11.gifto be in the (max) domain of attraction of the Gumbel distribution is that

http://freakonometrics.hypotheses.org/files/2015/12/ext-13.gif

Then in 1943, Boris Gnedenko gave a complete characterization of those three types, with a complete characterization for two of them (heavy tails, i.e. Fréchet type and bounded support, i.e. Weibull) but his necessary and sufficient condition was based on a function that was not explicitly defined (see Gnedenko (1943)). Laurens de Haan in the 70’s derived checkable condition for Gumbel’s type.
Boris Gnedenko proved (in Section 4 of his paper) that F is the (max) domain of attraction of http://freakonometrics.hypotheses.org/files/2015/12/ext-10.gif if and only if http://freakonometrics.hypotheses.org/files/2015/12/ext-16.gif is regularly varying at infinity, with index http://freakonometrics.hypotheses.org/files/2015/12/ext-17.gif (even if the term “regular variation” was not mentioned in the paper). Similar results were derived to characterize functions in the (max) domain of attraction of Weibull. For the (max) domain of attraction of http://freakonometrics.hypotheses.org/files/2015/12/ext-18.gif, Boris Gnedenko obtained that a necessary and sufficient condition was that there exists a function http://freakonometrics.hypotheses.org/files/2015/12/ext-19.gif such http://freakonometrics.hypotheses.org/files/2015/12/ext-19.gif goes to 0 at infinity and

http://freakonometrics.hypotheses.org/files/2015/12/ext-20.gif

Several papers have discussed what function http://freakonometrics.hypotheses.org/files/2015/12/ext-19.gif could be e.g. David Mejzler in 1949 (in Russian, but see also his 1965 paper), and Laurens de Hann in 1970 and 1971 (following the dramatic flood in the Netherlands in 1953, researchers in the Netherlands have focuses on dikes, and extreme value applications).

Mejzler’s idea was to work on quantiles, and not on the cumulative distribution function. I.e. define

http://freakonometrics.hypotheses.org/files/2015/12/ext-21.gif

Then a necessary and sufficient condition for F to be in the (max) domain of attraction of http://freakonometrics.hypotheses.org/files/2015/12/ext-18.gif is that

http://freakonometrics.hypotheses.org/files/2015/12/ext-23.gif

Laurens de Haan proved in 1971 that function http://freakonometrics.hypotheses.org/files/2015/12/ext-19.gif can be – in general – given by

http://freakonometrics.hypotheses.org/files/2015/12/ext-25.gif

And in 1976, Laurens de Haan obtained a three-type convergence working on quantile function http://freakonometrics.hypotheses.org/files/2015/12/ext-26.gif (with a much shorter proof).
There have been many many papers extending Fisher-Tippett’s theorem, e.g. on non-independent sequences, like exchangeable ones (in a paper by Simeon Berman in 1962, or on stationary Gaussian sequences in 1964).

Fisher-Tippett theorem and limiting distribution for the maximum

Tomorrow, we will discuss Fisher-Tippett theorem. The idea is that there are only three possible limiting distributions for normalized versions of the maxima of i.i.d. samples https://freakonometrics.hypotheses.org/files/2018/02/max-00.gif. For bounded distribution, consider e.g. the uniform distribution on the unit interval, i.e. https://freakonometrics.hypotheses.org/files/2018/02/max-09.gif on the unit interval. Let https://freakonometrics.hypotheses.org/files/2018/02/max-10.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-11.gif. Then, for all https://freakonometrics.hypotheses.org/files/2018/02/max-12.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-13.gif,

https://freakonometrics.hypotheses.org/files/2018/02/max-14.gif

i.e. the limiting distribution of the maximum is Weibull’s.

set.seed(1)
s=1000000
n=100
M=matrix(runif(s),n,s/n)
V=apply(M,2,max)
bn=1
an=1/n
U=(V-bn)/an
hist(U,probability=TRUE,,col="light green",
xlim=c(-7,1),main="",breaks=seq(-20,10,by=.25))
u=seq(-10,0,by=.1)
v=exp(u)
lines(u,v,lwd=3,col="red")

For heavy tailed distribution, or Pareto-type tails, consider Pareto samples, with distribution function https://freakonometrics.hypotheses.org/files/2018/02/max-05.gif. Let https://freakonometrics.hypotheses.org/files/2018/02/max-06.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-07.gif, then

https://freakonometrics.hypotheses.org/files/2018/02/max-08.gif

which means that the limiting distribution is Fréchet’s.

set.seed(1)
s=1000000
n=100
M=matrix((runif(s))^(-1/2),n,s/n)
V=apply(M,2,max)
bn=0
an=n^(1/2)
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(0,7),main="",breaks=seq(0,max(U)+1,by=.25))
u=seq(0,10,by=.1)
v=dfrechet(u,shape=2)
lines(u,v,lwd=3,col="red")

For light tailed distribution, or exponential tails, consider e.g. a sample of exponentially distribution variates, with common distribution function https://freakonometrics.hypotheses.org/files/2018/02/max-01.gif. Let https://freakonometrics.hypotheses.org/files/2018/02/max-02.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-03.gif, then

https://freakonometrics.hypotheses.org/files/2018/02/max-04.gif

i.e. the limiting distribution for the maximum is Gumbel’s distribution.

library(evd)
set.seed(1)
s=1000000
n=100
M=matrix(rexp(s,1),n,s/n)
V=apply(M,2,max)
(bn=qexp(1-1/n))
log(n)
an=1
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Consider now a Gaussian https://freakonometrics.hypotheses.org/files/2018/02/max-17.gif sample. We can use the following approximation of the cumulative distribution function (based on l’Hopital’s rule)

https://freakonometrics.hypotheses.org/files/2018/02/max-15.gif

as https://freakonometrics.hypotheses.org/files/2018/02/max-16.gif. Let https://freakonometrics.hypotheses.org/files/2018/02/max-18.gif and https://freakonometrics.hypotheses.org/files/2018/02/max-19.gif. Then we can get

https://freakonometrics.hypotheses.org/files/2018/02/max-20.gif

as https://freakonometrics.hypotheses.org/files/2018/02/max-21.gif. I.e. the limiting distribution of the maximum of a Gaussian sample is Gumbel’s. But what we do not see here is that for a Gaussian sample, the convergence is extremely slow, i.e., with 100 observations, we are still far away from Gumbel distribution,

and it is only slightly better with 1,000 observations,

set.seed(1)
s=10000000
n=1000
M=matrix(rnorm(s,0,1),n,s/n)
V=apply(M,2,max)
(bn=qnorm(1-1/n,0,1))
an=1/bn
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Even worst, consider lognormal observations. In that case, recall that if we consider (increasing) transformation of variates, we are in the same domain of attraction. Hence, since https://freakonometrics.hypotheses.org/files/2018/02/max-22.gif, if

https://freakonometrics.hypotheses.org/files/2018/02/max-23.gif

then

https://freakonometrics.hypotheses.org/files/2018/02/max-24.gif

i.e. using Taylor’s approximation on the right term,

https://freakonometrics.hypotheses.org/files/2018/02/max-25.gif

This gives us normalizing coefficients we should use here.

set.seed(1)
s=10000000
n=1000
M=matrix(rlnorm(s,0,1),n,s/n)
V=apply(M,2,max)
bn=exp(qnorm(1-1/n,0,1))
an=exp(qnorm(1-1/n,0,1))/(qnorm(1-1/n,0,1))
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,40,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Confidence interval for predictions with GLMs

Consider a (simple) Poisson regression http://freakonometrics.hypotheses.org/files/2016/11/poiss01.gif. Given a sample http://freakonometrics.hypotheses.org/files/2016/11/poiss02.gif where http://freakonometrics.hypotheses.org/files/2016/11/poiss03.gif, the goal is to derive a 95% confidence interval for http://freakonometrics.hypotheses.org/files/2016/11/poiss04.gif given http://freakonometrics.hypotheses.org/files/2016/11/poiss05.gif, where http://freakonometrics.hypotheses.org/files/2016/11/poiss04.gif is the prediction. Hence, we want to derive a confidence interval for the prediction, not the potential observation, i.e. the dot on the graph below

> r=glm(dist~speed,data=cars,family=poisson)
> P=predict(r,type="response",
+ newdata=data.frame(speed=seq(-1,35,by=.2)))
> plot(cars,xlim=c(0,31),ylim=c(0,170))
> abline(v=30,lty=2)
> lines(seq(-1,35,by=.2),P,lwd=2,col="red")
> P0=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> points(30,P1$fit,pch=4,lwd=3)

i.e.

Let http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif denote the maximum likelihood estimator of http://freakonometrics.hypotheses.org/files/2016/11/poiss07.gif. Then
http://freakonometrics.hypotheses.org/files/2016/11/poiss40.gif
where http://freakonometrics.hypotheses.org/files/2016/11/poiss101.gif is Fisher information of http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif (from standard maximum likelihood theory). Recall that
http://freakonometrics.hypotheses.org/files/2016/11/poiss13.gif
where computation of those values is based on the following calculations
http://freakonometrics.blog.fre<br /><br /> e.fr/public/latex/poiss21.gif
In the case of the log-Poisson regression
http://freakonometrics.hypotheses.org/files/2016/11/poiss36.gif
Let us get back to our initial problem.

  • confidence interval for the linear combination

A first idea to get a confidence interval for http://freakonometrics.hypotheses.org/files/2016/11/poiss49.gif is to get a confidence interval for http://freakonometrics.hypotheses.org/files/2016/11/poiss100.gif (by taking exponential values of bounds, since the exponential is a monotone function). Asymptotically, we know that
http://freakonometrics.hypotheses.org/files/2016/11/poiss40.gif

thus, an approximation for the variance matrix of http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif will be based on http://freakonometrics.hypotheses.org/files/2016/11/poiss45.gif, obtained by plugging estimators of the parameters.
Then, since http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif as an asymptotic multivariate distribution, any linear combination of the parameters will also be normal, i.e.
http://freakonometrics.hypotheses.org/files/2016/11/poiss47.gif has a normal distribution, centered on http://freakonometrics.hypotheses.org/files/2016/11/poiss49.gif, with variance http://freakonometrics.hypotheses.org/files/2016/11/poiss102.gif where http://freakonometrics.hypotheses.org/files/2016/11/Poiss110.gif is the variance of http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif. All those quantities can be easily computed. First, we can get the variance of the estimators

> i1=sum(predict(reg,type="response"))
> i2=sum(cars$speed*predict(reg,type="response"))
> i3=sum(cars$speed^2*predict(reg,type="response"))
> I=matrix(c(i1,i2,i2,i3),2,2)
> V=solve(I)

Hence, if we compare with the output of the regression,

> summary(reg)$cov.unscaled
(Intercept)         speed
(Intercept)  0.0066870446 -3.474479e-04
speed       -0.0003474479  1.940302e-05
> V
[,1]          [,2]
[1,]  0.0066871228 -3.474515e-04
[2,] -0.0003474515  1.940318e-05

Based on those values, it is easy to derive the standard deviation for the linear combination,

> x=30
> P2=predict(r,type="link",se.fit=TRUE,
+ newdata=data.frame(speed=x))
> P2
$fit
1
5.046034

$se.fit
[1] 0.05747075

$residual.scale
[1] 1

> sqrt(V[1,1]+2*x*V[2,1]+x^2*V[2,2])
[1] 0.05747084
> sqrt(t(c(1,x))%*%V%*%c(1,x))
[,1]
[1,] 0.05747084

And once we have the standard deviation, and normality (at least asymptotically), confidence intervals are derived, and then, taking the exponential of the bounds, we get confidence interval

> segments(30,exp(P2$fit-1.96*P2$se.fit),
+ 30,exp(P2$fit+1.96*P2$se.fit),col="blue",lwd=3)

Based on that technique, confidence intervals are no longer centered on the prediction. But who cares ?

  • delta method

Actually, those who like to use “more or less” expressions for confidence intervals will not like non centered intervals. So, an alternative is to use the delta method. Instead of writing (again) something on the theory, we can use a package which computes that method,

> estmean=t(c(1,x))%*%coef(reg)
> var=t(c(1,x))%*%summary(reg)$cov.unscaled%*%c(1,x)
> library(msm)
> deltamethod (~ exp(x1), estmean, var)
[1] 8.931232
> P1=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> P1
$fit
1
155.4048

$se.fit
1
8.931232

$residual.scale
[1] 1

The delta method gives us (asymptotic) normality, so once we have a standard deviation, we get the confidence interval.

> segments(30,P1$fit-1.96*P1$se.fit,30,
+ P1$fit+1.96*P1$se.fit,col="blue",lwd=3)

Note that those quantities – obtained with two different approaches – are rather close here

> exp(P2$fit-1.96*P2$se.fit)
1
138.8495
> P1$fit-1.96*P1$se.fit
1
137.8996
> exp(P2$fit+1.96*P2$se.fit)
1
173.9341
> P1$fit+1.96*P1$se.fit
1
172.9101
  • bootstrap techniques

And a third method (but far from what I expect to teach on that course) is to use bootstrap techniques to about those results based on asymptotic normality (we have only 50 observations). The idea is to sample from out dataset, and to run a log-Poisson regression on those new samples, and to repeat a lot of time,