# Why pictures are so important when modeling data?

(bis repetita) Consider the following regression summary,

Call:
lm(formula = y1 ~ x1)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0001     1.1247   2.667  0.02573 *
x1            0.5001     0.1179   4.241  0.00217 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared: 0.6665,	Adjusted R-squared: 0.6295
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217
obtained from
> data(anscombe)
> reg1=lm(y1~x1,data=anscombe)
Can we say something if we look (only) at that output ? The intercept is significatively non-null, as well as the slope, the $R^2$ is large (66%). It looks like we do have a nice model here. And in a perfect world, we might hope that data are coming from this kind of dataset,
But it might be possible to have completely different kinds of patterns. Actually, four differents sets of data are coming from Anscombe (1973). And that all those datasets are somehow equivalent: the $X$‘s have the same mean, and the same variance
> apply(anscombe[,1:4],2,mean)
x1 x2 x3 x4
9  9  9  9
> apply(anscombe[,1:4],2,var)
x1 x2 x3 x4
11 11 11 11
and so are the $Y$‘s
> apply(anscombe[,5:8],2,mean)
y1       y2       y3       y4
7.500909 7.500909 7.500000 7.500909
> apply(anscombe[,5:8],2,var)
y1       y2       y3       y4
4.127269 4.127629 4.122620 4.123249
Further, observe also that the correlation between the $X$‘s and the $Y$‘s is the same
> cor(anscombe)[1:4,5:8]
y1         y2         y3         y4
x1  0.8164205  0.8162365  0.8162867 -0.3140467
x2  0.8164205  0.8162365  0.8162867 -0.3140467
x3  0.8164205  0.8162365  0.8162867 -0.3140467
x4 -0.5290927 -0.7184365 -0.3446610  0.8165214
> diag(cor(anscombe)[1:4,5:8])
[1] 0.8164205 0.8162365 0.8162867 0.8165214
which yields the same regression line (intercept and slope)
> cbind(coef(reg1),coef(reg2),coef(reg3),coef(reg4))
[,1]     [,2]      [,3]      [,4]
(Intercept) 3.0000909 3.000909 3.0024545 3.0017273
x1          0.5000909 0.500000 0.4997273 0.4999091
But there is more. Much more. For instance, we always have the standard deviation for residuals
> c(summary(reg1)$sigma,summary(reg2)$sigma,
+ summary(reg3)$sigma,summary(reg4)$sigma)
[1] 1.236603 1.237214 1.236311 1.235695
Thus, all regressions here have the same R2
> c(summary(reg1)$r.squared,summary(reg2)$r.squared,
+ summary(reg3)$r.squared,summary(reg4)$r.squared)
[1] 0.6665425 0.6662420 0.6663240 0.6667073
Finally, Fisher’s F statistics is also (almost) the same.
+ c(summary(reg1)$fstatistic[1],summary(reg2)$fstatistic[1],
+ summary(reg3)$fstatistic[1],summary(reg4)$fstatistic[1])
value    value    value    value
17.98994 17.96565 17.97228 18.00329
Thus, with the following datasets, we have the same prediction (and the same confidence intervals). Consider for instance the second dataset (the first one being mentioned above),
> reg2=lm(y2~x2,data=anscombe)
The output is here exactly the same as the one we had above
> summary(reg2)

Call:
lm(formula = y2 ~ x2, data = anscombe)

Residuals:
Min      1Q  Median      3Q     Max
-1.9009 -0.7609  0.1291  0.9491  1.2691

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    3.001      1.125   2.667  0.02576 *
x2             0.500      0.118   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared: 0.6662,	Adjusted R-squared: 0.6292
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002179
Here, the perfect model is the one obtained with a quadratic regression.
> reg2b=lm(y2~x2+I(x2^2),data=anscombe)
> summary(reg2b)

Call:
lm(formula = y2 ~ x2 + I(x2^2), data = anscombe)

Residuals:
Min         1Q     Median         3Q        Max
-0.0013287 -0.0011888 -0.0006294  0.0008741  0.0023776

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.9957343  0.0043299   -1385   <2e-16 ***
x2           2.7808392  0.0010401    2674   <2e-16 ***
I(x2^2)     -0.1267133  0.0000571   -2219   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.001672 on 8 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 7.378e+06 on 2 and 8 DF,  p-value: < 2.2e-16
Consider now the third one
> reg3=lm(y3~x3,data=anscombe)
i.e.
> summary(reg3)

Call:
lm(formula = y3 ~ x3, data = anscombe)

Residuals:
Min      1Q  Median      3Q     Max
-1.1586 -0.6146 -0.2303  0.1540  3.2411

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0025     1.1245   2.670  0.02562 *
x3            0.4997     0.1179   4.239  0.00218 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared: 0.6663,	Adjusted R-squared: 0.6292
F-statistic: 17.97 on 1 and 9 DF,  p-value: 0.002176
This time, the linear model could have been perfect. The problem is one outlier. If we remove it, we have
> reg3b=lm(y3~x3,data=anscombe[-3,])
> summary(reg3b)

Call:
lm(formula = y3 ~ x3, data = anscombe[-3, ])

Residuals:
Min         1Q     Median         3Q        Max
-0.0041558 -0.0022240  0.0000649  0.0018182  0.0050649

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.0056494  0.0029242    1370   <2e-16 ***
x3          0.3453896  0.0003206    1077   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.003082 on 8 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 1.161e+06 on 1 and 8 DF,  p-value: < 2.2e-16
Finally consider
> reg4=lm(y4~x4,data=anscombe)
This time, there is an other kind of outlier, in $X$‘s, but again, the regression is exactly the same,
> summary(reg4)

Call:
lm(formula = y4 ~ x4, data = anscombe)

Residuals:
Min     1Q Median     3Q    Max
-1.751 -0.831  0.000  0.809  1.839

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   3.0017     1.1239   2.671  0.02559 *
x4            0.4999     0.1178   4.243  0.00216 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 9 degrees of freedom
Multiple R-squared: 0.6667,	Adjusted R-squared: 0.6297
F-statistic:    18 on 1 and 9 DF,  p-value: 0.002165
The graph is here
So clearly, looking at the summary of a regression does not tell us anything… This is why we do spend some time on diagnostic, looking at graphs with the errors (the graphs above could be obtained only with one explanatory variable, while errors can be studied in any dimension): everything can be seen on thise graphs. E.g. for the first dataset,
or the second one
the third one
or the fourth one,

# Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 2)

Continuons la discussion de ce matin, avec l’exemple des capitales américaines. On peut trouver une telle liste liste en ligne sur wikipedia,

> capital=read.table(
+ "http://freakonometrics.free.fr/US-capital.csv",
> capital=capital[-c(2,11),]

On va enlever Hawaï et l’Alaska, qui forcent à faire de trop grands détours. Dans TSP, sous R, il y a des listes de villes américaines, avec les coordonnées. En bricolant un peu, on va récupérer les coordonnées de la plupart des capitales,

> library(maps)
> library(maptools)
> library(sp)
> library(TSP)
> Lcap=paste(as.character(capital[,1]),
+ as.character(capital[,2]),sep=", ")
> M=as.matrix(USCA312)
> L=labels(M)[[1]]
> Lcap=Lcap[which(Lcap%in%L)]
> k=which(L%in%Lcap)
> M=M[k,]
> M=M[,k]
> listeUSA=TSP(M,L[k])

Et cette fois on est prêt,

> tour=solve_TSP(listeUSA, method = "nn")

on a lancé l’algorithme, et on peut faire un dessin,

> COORD=data.frame(USCA312_coords[k,])
> COORD=COORD[as.numeric(tour),]
> map('state',fill=TRUE,col="yellow")
> lines(COORD$coords.x1,COORD$coords.x2,
+ lwd=3,col="blue")
> points(COORD$coords.x1,COORD$coords.x2,
+ pch = 19, cex = 1.2, col = "red")

pas mal, n’est pas ? Bon, mais comme d’habitude, rien ne garantie que l’on ait trouvé le trajet optimal… En fait, le soucis est que si on relance l’algorithme, on retombe sur un autre trajet,

voire encore un autre si on lancer une nouvelle fois l’algorithme,

C’est pénible, n’est ce pas ? En fait, si on lance 1,000 fois l’algorithme on obtient la distribution (des distances optimales) suivante

autrement dit, on a encore beaucoup de variabilité, avec une distance optimale qui peut varier d’environ 10% entre deux boucles d’algorithmes. Cela dit, cela reste faible par rapport à un trajet au hasard,

Et d’ailleurs, si on compare la distance obtenue en tirant les villes au hasard, on arrive à un taux de l’ordre de 27% (qui peut être comparé à ce que nous avions obtenu en tirant au hasard dans le carré unité)

>  HAZARD=OPTIMAL=RATIO=rep(NA,500)
>  for(s in 1:500){
+   HAZARD[s]=tsp.longueur(M,sample(1:nrow(M)))
+   tour=solve_TSP(listeUSA, method = "nn")
+   OPTIMAL[s]=tsp.longueur(M,as.numeric(tour))
+   RATIO[s]=HAZARD[s]/
+            OPTIMAL[s]
+ }

Bref, compte tenu de la (grande) dispersion des résultats optimaux obtenu, on pourrait dire que cet algorithme est assez mauvais. Même s’il est très rapide. Un autre avantage est que l’on peut spécifier la ville de départ, par exemple. Ici, si on souhaite partir de Floride, on peut obtenir le trajet suivant

> tour=solve_TSP(listeUSA, method = "nn",
+ control=list(start=41))

Passionnant n’est-ce pas ? Et encore, ça ne va pas m’aider à trouver le parcours optimal pour les enfants, car pour les enfants, il faut spécifier le point de départ (voire le point d’arrivée), il faut en plus rester sur la route (voire traverser aux passages pour piétons). On peut aussi vouloir faire une pause pipi au milieu de la promenade…. Bref, on a pas mal de contraintes ! Mais tout est expliqué dans le livre fabuleux de William Cook, In Pursuit of the Traveling Salesman, qui reprend toute l’histoire des algorithmes les plus fins pour résoudre ce joli problème de maths…

# Quel trajet minimal pour récupérer le maximum de bonbons à l’Halloween ? (partie 1)

Bon, ben cette fois ça y est, c’est l’Halloween… et comme tous les ans, mon rôle devrait se limiter à creuser les citrouilles, décorer la maison… et attendre à la porte de la maison que les enfants des voisins viennent me dévaliser ma réserve de bonbons ! Je pourrais à la rigueur revendiquer un petit quelque chose dans le déguisement de mon fils (on a passé pas mal de temps à feuilleter les premiers tomes de Walking Deads, histoire de mieux saisir l’essence des personnages de morts-vivants, en plus de la Zombie Walk la fin de semaine passée, en centre-ville).

Mais cette année, pas de billet sur les zombies comme l’an passé, mais un vrai problème pratique: comment optimiser son trajet, quand on sait où sont les bonnes maisons, où on a de fortes chances d’avoir plein de bonbons (en particulier chez les parents des amis) ? C’est assez proche de l’idée de passer par un maximum de maisons, et ce à une heure raisonnable (histoire qu’il reste des bonbons à récupérer).

Pour tous ceux qui ont fait un peu de programmation dynamique, c’est un problème classique, connu sous le nom du “voyageur de commerce“. Ce problème est né en juillet 1954 avec “the shortest route for a traveling salesman” paru dans Newsweek ci-dessous (que l’on retrouve en détails dans le livre 50 Years of Integer Programming)

Enfin, disons qu’il a été popularisé à cette époque, car Euler évoquait déjà un problème similaire quand il travaillait sur le problème des sept ponts de Königsberg, dans “Solutio problematis ad geometriam situs pertinentis“,

(même si ce problème consiste à étudier l’existence d’un trajet sous des contraintes de passages par des nœuds, et pas vraiment de l’optimisation). On retrouve aussi un exemple en 1832, en Allemagne, dans un livre intitulé Der Handlungsreisende – wie er sein soll und was er zu thun hat, um Auftraege zu erhalten und eines gluecklichen Erfolgs in seinen Geschaeften gewiss zu sein – Von einem alten Commis-Voyageur où un trajet passant par 45 villes en Allemagne et en Suisse devait être optimisé

(le schéma ci-dessus montre celle décrite en 1832, pour un total de 1,285km, mais il est possible de faire une boucle de 1,248km).

En fait le problème a explosé en 1962 (à cette époque, il semble que les gens pouvaient encore penser que faire des maths pouvait être le fun), lors que Procter & Gamble ont offert un prix de 10,000$à celui qui trouverait le chemin le plus court pour que Toody et Muldoom, les conducteurs de la voiture 54 dans une série populaire à l’époque, passent par 33 villes, aux États-Unis. Cela dit, Carl Menger avait noté dès 1930 a noté que ce problème devrait etre difficile à résoudre. Facile à écrire, mais difficile à résoudre explicitement. Dans son livre, William Cooke rappelle qu’avec 33 villes, il y aurait 131,565,418,466,846,765,083,609,006,080,000,000 chemins possibles, dont il faudrait calculer la longueur. Et qu’avec l’ordinateur le plus puissant en 2009 (effectuant 1.5 million de milliards d’opérations à la seconde) il faudrait 28,000 milliards d’année pour mener à bien les calculs. Même si on essayait de tenir compte des progrès de l’informatique dans les mille prochaines années, ça risque de prendre du temps !.. En 1967, Jack Edmonds notait “I conjecture that there is no good algorithm for the traveling salesman problem” où “good” est une manière élégante pour dire un temps polynomial en $n$, le nombre de villes. Ce qui s’appelle les problèmes de la classe $\mathcal{P}$. Il semblerait que le problème du voyageur de commerce soit un problème $\mathcal{NP}$ et pas de la classe $\mathcal{P}$(cf le chapitre du livre de Sanjeev Arora et Boaz Barak, ou le chapitre 11 des notes de cours de Olivier Bournez, que j’ai pu découvrir grâce à @dmonniaux, ou formellement de la classe $\mathcal{NP}$ et pas $\mathcal{P}$, i.e. $\mathcal{NP}$ complet), voir aussi la discussion sur wikipedia), Et depuis, l’institut Clay a proposé$1,000,000 à celui qui prouverait que les problèmes $\mathcal{NP}$  ne sont pas de type $\mathcal{P}$ (ou qui prouveraient qu’ils le sont… laconjecture reste ouverte, semble-t-il). Bref, les enjeux financiers derrière ont considérablement augmenté (on pourra consulter Sipser (2007) pour la petite histoire)

Bon, et si l’histoire est intéressante, ce n’est pas une raison pour la lire en se croisant les bras… essayons de programmer un algorithme (simple) afin d’approcher une solution, car je ne pourrais me contenter de dire à mes enfants “c’est un problème $\mathcal{NP}$ complet, papa a capitulé“.

Commençons par tirer 5 points au hasard dans le carré unité,

>  n = 5
>  (x = matrix(runif(2*n),nr=n))
[,1]       [,2]
[1,] 0.4142621 0.35600725
[2,] 0.9477231 0.66938653
[3,] 0.2518656 0.19900873
[4,] 0.5314231 0.22954670
[5,] 0.5156522 0.09066741
et commençons par calculer la distance entre les points (i.e. une matrice, ou ici une matrice triangulaire inférieure)
>  (d = dist(x))
1         2         3         4
2 0.6186980
3 0.2258786 0.8399243
4 0.1723919 0.6056111 0.2812205
5 0.2840514 0.7222195 0.2851688 0.1397719
On va mettre ce triangle sous forme de matrice (pleine), ça sera plus simple à manipuler par la suite
>  (d = as.matrix(d))
1         2         3         4         5
1 0.0000000 0.6186980 0.2258786 0.1723919 0.2840514
2 0.6186980 0.0000000 0.8399243 0.6056111 0.7222195
3 0.2258786 0.8399243 0.0000000 0.2812205 0.2851688
4 0.1723919 0.6056111 0.2812205 0.0000000 0.1397719
5 0.2840514 0.7222195 0.2851688 0.1397719 0.0000000
Le principe, pour initialiser l’algorithme est de tirer au hasard un trajet, c’est à dire une permutation de nos points,
>  (o = sample(1:n))
[1] 2 4 1 5 3
Il faut alors calculer le 4 distances entre ces 5 villes (on ne prévoit pas de retour à la ville de départ, pour l’instant). Le plus simple est de jouer avec cette matrice de la réarranger,
>  d[o[1:(n-1)],][,o[2:n]]
4         1         5         3
2 0.6056111 0.6186980 0.7222195 0.8399243
4 0.0000000 0.1723919 0.1397719 0.2812205
1 0.1723919 0.0000000 0.2840514 0.2258786
5 0.1397719 0.2840514 0.0000000 0.2851688
et de prendre les 4 valeurs sur la diagonale
>  diag(d[o[1:(n-1)],][,o[2:n]])
[1] 0.6056111 0.1723919 0.2840514 0.2851688
Ce sont les 4 distances que l’on cherche, entre les 5 villes. La distance totale parcourue est alors
>  sum(diag(d[o[1:(n-1)],][,o[2:n]]))
[1] 1.347223
L’idée est ensuite de permuter deux des villes
>  (i=sample(1:n,2))
[1] 3 1
(ici la première et la troisième de la liste, mais on pourrait décider de fixer la première si on doit partir de là) et de voir si la nouvelle distance est plus courte, ou pas,
>  os=o
>  os[i]=o[rev(i)]
>  sum(diag(d[os[1:(n-1)],][,os[2:n]]))
[1] 1.785391
La distance étant plus longue, on ne va pas permuter. On pourrait tenter 1000 fois de permuter deux villes. Si on ne parvient pas à réduire la distance totale, c’est qu’on a atteint une valeur intéressante (à défaut de pouvoir prouver qu’elle est optimale). Et si on peut améliorer la distance, on permute. Le plus simple est donc de faire des petites fonctions. La première reprend le calcul de la distance totale, sur la diagonale de la matrice permutée, à matrice de distance donnée et à permutation donnée,
> tsp.longueur=function(matrice,ordres) {
+ n=length(ordres)
+ sum(diag(matrice[ordres[1:(n-1)],][,ordres[2:n]])) }
La seconde est juste une boucle,
> tsp.optimal=function(d,N=1000) {
+ d=as.matrix(d)
+ n=ncol(d)
+ o=sample(1:n)
+ v=tsp.longueur(d,o)
+ k=0
+ while(k < N) {
+   i<-sample(1:n,2)
+   os=o
+   os[i]=o[rev(i)]
+   w=tsp.longueur(d,os)
+   if(w < v) {
+     v=w
+     o=os
+     k=0} else {k=k+1}}
+ list(ordre=o,longueur=v)}
On peut se lancer maintenant, avec par exemple 30 points simulés, sur le carré unité,
>  set.seed(1)
>  n=30
>  x=matrix(runif(2*n),nr=n)
>  o=sample(1:n)
>  plot(x,xlim=0:1,ylim=0:1,xlab="",ylab="")
>  lines(x[o,],col='blue')
>  r=tsp.optimal(dist(x))
>  os=r$ordre > lines(x[os,],col='blue',lwd=1.5) avec la première trajectoire, à gauche, qui reste en trait fin sur le dessin de droite, Pas mal, non ? Bon, en fait le soucis est que si on joue plusieurs fois avec notre fonction, on obtient toujours un graphique différent, Autrement dit, je trouve (certes) des chemins pas trop long, mais je suis loin de trouver le chemin optimal. Et histoire de clotûrer le débat, notons qu’on peut faire un peu de simulation, afin de comparer la distance obtenue en faisant le trajet au hasard, et un trajet “optimal“, > RATIO=HASARD=OPTIMAL=matrix(NA,500,10) > for(m in 1:10){ + n=5*m + for(s in 1:500){ + x=matrix(runif(2*n),nr=n) + r=tsp.opt(dist(x)) + HAZARD[s,m]=tsp.longueur(as.matrix(dist(x)),1:n) + OPTIMAL[s,m]=r$longueur
+  RATIO[s,m]=HAZARD[s,m]/
+             OPTIMAL[s,m]
+ }
+}

Si on compare la longueur moyenne du trajet (sur 500 scénarios) avec un tour au hasard (en rouge) et ou tour optimisé (en bleu), on a

ce qui donne un ratio moyen (ou une moyenne des ratios, les deux sont représentés)

Autrement dit, avec une dizaine de points éparpillés (au hasard) dans le carré unité, un tour au hasard ne sera que deux fois plus long… étonnant non ?

Mais on doit pouvoir aller plus loin, parce que sous R, il y a des fonctions (et des packages) pour faire des algorithmes plus compliqués… à suivre avant la tournée des bonbons de demain soir !

# Le passage au log dans les modèles linéaires

Un billet rapide pour compléter et illustrer le passage au log dans un modèle linéaire (que l’on abordera cette semaine en cours). Le point de départ est le modèle linéaire, où on suppose que, conditionnellement à $X$$Y$ suit une loi normale. Pour rappel, si on a une loi normale, $Y\sim \mathcal{N}(\mu,\sigma^2)$, alors $\mathbb{E}(Y)=\mu$ et $\text{Var}(Y)=\sigma^2$. Les intervalles de confiance à 90% et 95% sont symétriques par rapport à la moyenne (qui est aussi la médiane, soit dit en passant),

Dans un modèle Gaussien avec homoscédasiticité, $Y|X\sim \mathcal{N}(\beta_0+\beta_1X,\sigma^2)$ i.e. $\mathbb{E}(Y|X)=\beta_0+\beta_1X$ alors que $\text{Var}(Y|X)=\sigma^2$. On a alors les bandes de confiance suivantes, pour un modèle de régression linéaire,

Bon, maintenant, que se passe-t-il si on prend l’exponentiel ? Pour la loi normale, rappelons que l’on obtient une loi lognormale, i.e. $\exp(Y)\sim LN(\mu,\sigma^2)$, les deux paramètres étant liés à la loi normale sous jacente, car désormais

$\mathbb{E}(\exp(Y))=\exp\left(\mu+\frac{\sigma^2}{2}\right)$

alors que

$\text{Var}(\exp(Y))=(\exp[\sigma^2]-1)\cdot\exp[2\mu+\sigma^2]$

Graphiquement, on a la loi suivante, avec les intervalles de confiance à 90% et 95% représentés ci-dessous. Le point noir est $\exp(\mu)$  alors que le point bleu est l’espérance de la loi lognormale.

On notera que le quantile de la loi log-normale est l’exponentiel du quantile de la loi normale. En effet, si $\mathbb{P}(Y\leq q)=\alpha$ alors $\mathbb{P}(\exp(Y)\leq \exp(q))=\alpha$. En particulier, $\exp(\mu)$ n’est pas la moyenne de $\exp(Y)$, mais la médiane (puisque $\mu$ était la médiane de $Y$).

Mais il n’est pas rare de voir utilisé un intervalle de confiance de la forme

$\mathbb(\exp(Y))\pm q_{1-\alpha/2}\sqrt{\text{Var}(\exp(Y))}$

qui est la forme classique de l’intervalle de confiance Gaussien (symétrique autour de la moyenne). Ici, on aurait les niveaux suivants

Notons qu’il n’y a aucune raison ici d’avoir une probabilité $1-\alpha$ d’être dans l’intervalle de confiance obtenu avec les quantiles $q_{1-\alpha/2}$ de la loi normale.

Maintenant, si on prend l’exponentiel d’un modèle linéaire (i.e. le logarithme de la variable d’intérêt est modélisé par un modèle linéaire) on a

$\mathbb{E}(\exp(Y)|X)=\exp\left(\beta_0+\beta_1X+\frac{\sigma^2}{2}\right)$

avec une variance (conditionnelle) qui dépend de la variable explicative

$\text{Var}(\exp(Y)|X)=(\exp[\sigma^2]-1)\cdot\exp[2(\beta_0+\beta_1X)+\sigma^2]$

Là encore, le plus naturel est d’utiliser comme bornes de l’intervalle de confiance des quantiles associés à la loi lognormale,

mais il n’est pas rare de voir utilisé des intervalles de type Gaussiens,

On perd là encore en interprétation car les bornes n’ont plus rien à voir avec les quantiles.

# ACT2121 premier examen

Cet après-midi avait lieu le premier examen d’actuariat 1, préparatoire pour l’examen P de la SOA. Comme pour l’examen P, il y avait des questions à choix multiples, en un temps (très limité). Le sujet est en ligne, ainsi que des éléments de corrections (bientôt) ici. Si vous n’êtes pas d’accord avec mes corrections, merci de me le dire rapidement, si possible avant que je n’attribue les notes…

# 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. # Somewhere else, part 18 A nice post (found a couple of days ago) on an insteresting topic, I should perhaps discuss more often on this blog Also, this week, a series of posts and articles about L’Aquila trial, Still, a lot of very interesting posts found here, and there, Also several documentaries, found, online Et comme toujours, quelques articles en francais, Did I miss something ? # 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 $p$-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 et sa version contrainte On cherche alors à tester l’hypothèse contre l’hypothèse alternative 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é 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 $R^2$ des régressions, contraintes et noncontraintes. Fisher a montré que sous l’hypothèse H0, 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 $p$-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). # Régression sur des variables orthogonales Hier, on a évoqué en cours le théorème de Frisch-Waugh (évoqué sur ce blog il y a quelques mois). Derrière, il y a l’estimation de modèles sur des sous-ensembles de valeurs. Formellement, on a un modèle de la forme $Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=\boldsymbol{X}_1\boldsymbol{\beta}_1+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon$ Les équations normales sont ici (cf transparents du cours, ou équations écrites au tableau) $\boldsymbol{X}'\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\boldsymbol{X}'Y$, que l’on peut écrire sous la forme d’un système, $\left\{\begin{array}{l}\boldsymbol{X}_1'\boldsymbol{X}_1 \widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_1'\boldsymbol{X}_2 \widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_1'Y \\ \boldsymbol{X}_2'\boldsymbol{X}_1 \widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_2'\boldsymbol{X}_2 \widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2'Y \end{array}\right.$ Le point intéressant est que si les variables dans le premier ensemble et les variables dans le second sont orthogonales, i.e. $\boldsymbol{X}_1'\boldsymbol{X}_2 =\boldsymbol{0}$ et $\boldsymbol{X}_2'\boldsymbol{X}_1 =\boldsymbol{0}$, alors les termes diagonaux du système disparaissent, et on retrouve les équations normales des régressions simples, i.e. $\left\{\begin{array}{l}\boldsymbol{X}_1'\boldsymbol{X}_1 \widehat{\boldsymbol{\beta}}_1=\boldsymbol{X}_1'Y \\ \boldsymbol{X}_2'\boldsymbol{X}_2 \widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2'Y \end{array}\right.$ On a ici une conséquence du théorème de Frisch-Waugh: si les variables sont orthogonales, les estimateurs dans la régression multiple coïncident avec ceux obtenus avec des régressions simples. Afin de mieux clarifier ce point, écrivons ce qui se passerait avec seulement deux variables. On va associer la constante à la première, i.e. $Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=(1 \ \boldsymbol{X}_1)\begin{pmatrix}\boldsymbol{\beta}_0 \\ \boldsymbol{\beta}_1\end{pmatrix}+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon$ Dans ce cas, les équations normales sont $\left\{\begin{array}{l}\begin{pmatrix}1 \\ \boldsymbol{X}'_1\end{pmatrix}(1 \ \boldsymbol{X}_1)\begin{pmatrix}\widehat{\boldsymbol{\beta}}_0 \\ \widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\begin{pmatrix}1 \\ \boldsymbol{X}'_1\end{pmatrix}\boldsymbol{X}_2 \widehat{\boldsymbol{\beta}}_2=\begin{pmatrix}1 \\ \boldsymbol{X}'_1\end{pmatrix}Y \\ \boldsymbol{X}_2'(1\ \boldsymbol{X}_1) \begin{pmatrix}\widehat{\boldsymbol{\beta}}_0 \\ \widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\boldsymbol{X}_2'\boldsymbol{X}_2 \widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2'Y \end{array}\right.$ Bref, il faut que $\boldsymbol X_1$ et $\boldsymbol X_2$ soient orthogonales pour que le terme sur la seconde diagonale disparaisse. C’est ce qu’on retrouve dans la plupart des livres d’économétrie, il suffit de chercher le théorème de Frisch-Waugh dans l’index (même si le théorème de Frisch-Waugh va beaucoup plus loin puisqu’il explique quels ajustement faire pour lier les estimateurs dans une régression simple, et dans le cas multiple). Mais est-ce un résultat purement théorique ou peut-on l’invoquer dans un cas pratique ? Autrement dit, si deux variables semblent orthogonales, que se passent-il numériquement ? Si je simule des variables explicatives indépendantes, les estimateurs coïncident-ils (comme le suggère la théorie) ? Considérons un modèle linéaire avec deux variables explicative que l’on supposera indépendantes, mais ayant toutes deux un impact sur la variable explicative. Le plus simple est de considérer un triple de variables telles que $\boldsymbol{Z}=\begin{pmatrix} Y \\ X_1\\X_2\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix} \mu_Y \\ \mu_{X_1}\\\mu_{X_2}\end{pmatrix},\begin{pmatrix} 1 & \star & \star \\ \star & 1&0\\\star&0&1\end{pmatrix}\right)$ où $\star$ désigne des valeurs non nulles. Par exemple, > set.seed(1) > n=100 > S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3) > S [,1] [,2] [,3] [1,] 1.0 0.4 0.3 [2,] 0.4 1.0 0.0 [3,] 0.3 0.0 1.0 > M=c(1,3,6) > Z=rmnorm(n,M,varcov=S) > Y=Z[,1]; X1=Z[,2]; X2=Z[,3] > base=data.frame(Y,X1,X2) > cor(base) Y X1 X2 Y 1.0000000 0.37796580 0.27859381 X1 0.3779658 1.00000000 -0.05060539 X2 0.2785938 -0.05060539 1.00000000 On a simulé ici deux variables indépendantes, mais toutes deux corrélées avec la variable d’intérêt. On notera que numériquement $\frac{1}{n}\sum_{i=1}^n X_{1,i}X_{2,i} \neq \frac{1}{n}\sum_{i=1}^n X_{1,i}\cdot\frac{1}{n}\sum_{i=1}^n X_{2,i}$ (autrement dit, en moyenne, on simulera des données indépendantes, mais dans ce cas précis, les données ne sont pas rigoureusement orthogonales). > sum(X1*X2)/n [1] 18.20194 > sum(X1)*sum(X2)/n^2 [1] 18.25097 Mais un test rapide nous pousse à accepter l’hypothèse d’indépendance entre les deux variables (sans trop d’hésitations). > cor.test(X1,X2) Pearsons product-moment correlation data: X1 and X2 t = -0.5016, df = 98, p-value = 0.6171 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.2445923 0.1472766 sample estimates: cor -0.05060539 Si on regarde maintenant la régression, > R=lm(Y~X1+X2,data=base) > summary(R) Call: lm(formula = Y ~ X1 + X2, data = base) Residuals: Min 1Q Median 3Q Max -1.87885 -0.57656 0.07957 0.52719 1.77407 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.71180 0.28246 16.682 < 2e-16 *** X1 0.37232 0.08441 4.411 2.66e-05 *** X2 0.25975 0.07755 3.350 0.00115 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7954 on 97 degrees of freedom Multiple R-squared: 0.2317, Adjusted R-squared: 0.2159 F-statistic: 14.63 on 2 and 97 DF, p-value: 2.803e-06 Quel est le lien avec les régressions simples, sur chacune des variables ? Si on régresse sur la première, on obtient > R1=lm(Y~X1,data=base) > summary(R1) Call: lm(formula = Y ~ X1, data = base) Residuals: Min 1Q Median 3Q Max -1.78973 -0.54060 0.04036 0.44138 1.94389 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.03165 0.27934 18.013 < 2e-16 *** X1 0.35802 0.08859 4.041 0.000106 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.8358 on 98 degrees of freedom Multiple R-squared: 0.1429, Adjusted R-squared: 0.1341 F-statistic: 16.33 on 1 and 98 DF, p-value: 0.0001058 et si on régresse sur la seconde > R2=lm(Y~X2,data=base) > summary(R2) Call: lm(formula = Y ~ X2, data = base) Residuals: Min 1Q Median 3Q Max -1.8652 -0.5507 0.1104 0.6109 2.1435 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.85053 0.12495 46.824 < 2e-16 *** X2 0.24244 0.08443 2.872 0.00501 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.867 on 98 degrees of freedom Multiple R-squared: 0.07761, Adjusted R-squared: 0.0682 F-statistic: 8.246 on 1 and 98 DF, p-value: 0.005005 Si les données ne sont pas rigoureusement orthogonales, alors $\widehat{\beta}_1$ dans le modèle multiple est légèrement différent de $\widehat{\beta}_1$ obtenu dans le modèle simple. Si on regarde la distribution de nos estimateurs (n’oublions pas qu’il y a une marge d’erreur), ils sont proches, ou au moins statistiquement proches (avec à gauche la distribution de $\widehat{\beta}_1$ et à droite celle de $\widehat{\beta}_2$) Si on génère 100 bases de données avec ce modèle (avec 500 observations à chaque fois) on obtient les estimateurs suivants (à gauche les $\widehat{\beta}_1$ et à droite les $\widehat{\beta}_2$, avec en abscisse le modèle multiple et en ordonnées, le modèle simple). (la valeur centrale étant la vraie valeur, utilisée pour les simulations). On pourrait se demander si les valeurs éloignées du barycentre (ou de cette première bissectrice) correspondent à des cas où l’hypothèse d’indépendance entre les deux variables est rejetée ? On peut mettre des couleurs différentes en fonction de la $p$-value: plus la $p$-value est faible (rejet de l’hypothèse d’indépendance) plus la couleur sera vive (et rouge), et plus la $p$-value est forte (acceptation de l’hypothèse d’indépendance) plus la couleur sera terne (et jaune, voire blanche), sur 1,000 simulations, on obtient Le code est ici > set.seed(1) > library(mnormt) > S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3) > BX1=BX2=matrix(NA,1000,2) > P=rep(NA,1000) > for(s in 1:1000){ + n=500 + Z=rmnorm(n,c(6,3,1),varcov=S) + Y=Z[,1]; X1=Z[,2]; X2=Z[,3] + base=data.frame(Y,X1,X2) + R=lm(Y~X1+X2,data=base) + BX1[s,1]=coefficients(R)[2] + BX2[s,1]=coefficients(R)[3] + BX1[s,2]=coefficients(lm(Y~X1,data=base))[2] + BX2[s,2]=coefficients(lm(Y~X2,data=base))[2] + P[s]=cor.test(X1,X2)$p.value
+ }
> COL=(heat.colors(100))
> plot(BX1[,1],BX1[,2],pch=19,cex=.6,col=COL[trunc(100*P)])
> abline(a=0,b=1,lty=2)

Autrement dit le fait que l’estimation dans le modèle multiple et les modèles simples diffèrent beaucoup n’a rien a voir avec la $p$-value d’un test d’indépendance entre les deux variables… étonnant non ? En revanche, si on rajoute une observation pour rendre le modèle parfaitement non-corrélé, par exemple,

> difference=sum(X1*X2)*(n+1)/sum(X1)-sum(X2)
> X1=c(X1,0)
> X2=c(X2,difference)
> n=n+1
> sum(X1*X2)/n
[1] 18.02172
> sum(X1)*sum(X2)/n^2
[1] 18.02172
On retrouve ici une corrélation parfaitement nulle,
> Y=c(Y,0)
>  base=data.frame(Y,X1,X2)
>  cor(base)
Y            X1            X2
Y  1.0000000  3.946091e-01  2.923412e-01
X1 0.3946091  1.000000e+00 -4.817519e-16
X2 0.2923412 -4.817519e-16  1.000000e+00
Si on regarde maintenant la régression multiple
>  R=lm(Y~X1+X2,data=base)
>  coefficients(R)
(Intercept)          X1          X2
-1.5028358   0.3589776   0.2531358
et les deux régressions simples,
>  R1=lm(Y~X1,data=base)
>  coefficients(R1)
(Intercept)          X1
0.02847362  0.35897764
>  R2=lm(Y~X2,data=base)
>  coefficients(R2)
(Intercept)          X2
-0.4334011   0.2531358
on obtient exactement les mêmes coefficient.
Moralité ? l’hypothèse d’orthogonalité entre les variables explicatives n’est pas à prendre à la légère: il faut une orthogonalité stricte pour que numériquement, la régression multiple soit équivalente aux régressions simples (comme le prédit la théorie). Mais de manière un peu surprenante, si on est très proche de l’indépendance, les valeurs numériques des estimateurs peut différer. Et ce, parfois de manière non négligeable, même si statistiquement, on valide l’hypothèse d’orthogonalité entre les variables explicatives.

# ACT2121, éléments de correction (2)

Un petit complément pour revenir sur un calcul fait hier, en classe, au tableau.

Une police d’assurance couvre des pertes X et Y qui sont des variables aléatoires continues fonction dont on donne la densité conjointes. Trouver l’espérance du remboursement s’il est de X + Y avec un déductible de 1.

Au tableau, je m’étais lancé dans des calculs un peu compliqué (notant que les variables étaient indépendantes, on pouvait utiliser la convolution pour trouver la loi de la somme facilement… mais on reviendra sur ce calcul plus tard). La solution la plus simple était la suivante: la loi de $(X,Y)$ admet pour densité $f_{X,Y}(x,y)=2x\boldsymbol{1}_{[0,1]\times[0,1]}(x,y)$. Le remboursement $R$ est la variable aléatoire

$R=h(X,Y)=(X+Y-1)\boldsymbol{1}_{]1,+\infty[}(X+Y)$

On peut alors faire le calcul de
$\mathbb{E}(R)=\int_{\mathbb{R}^2} h(x,y)f(x,y)dxdy$
i.e.
$\mathbb{E}(R)=\int_0^1 \int_{1-x}^1 [x+y-1]2x dy dx$
soit
$\int_0^1 2x\left[(x-1)y+\frac{y^2}{2}\right]_{1-x}^1 dx=\int_0^1 x^3dx=\frac{1}{4}$
Voilà pour la version simple (et finalement assez générale). Pour la version plus compliquée
$f_{X,Y}(x,y)=2x\boldsymbol{1}_{[0,1]\times[0,1]}(x,y)=\underbrace{2x \boldsymbol{1}_{[0,1]}(x)}\underbrace{ \boldsymbol{1}_{[0,1]}(y)}$
autrement dit, on a obtenu que les deux composantes étaient indépendantes, et que $Y$ est uniformément distribuée sur $[0,1]$. Aussi, pour la loi de la somme
$\mathbb{P}(S>s)=\int_{0}^1 \mathbb{P}(S>s|Y=y)f_Y(y)dy=\int_{0}^1 \mathbb{P}(X>s-y|Y=y)dy$
et par indépendance,

$\mathbb{P}(X>s-y|Y=y)=\mathbb{P}(X>s-y)=\int_{s-y}^1 f_X(x)dx$

i.e.
$\mathbb{P}(X>s-y|Y=y)= \int_{s-y}^1 2x\boldsymbol{1}_{[0,1]}(x) dx=\int_{\max\{s-y,0\}}^1 2xdx$
soit

$\mathbb{P}(X>s-y|Y=y)=\left[x^2\right]_{\max\{s-y,0\}}^1=1-(\max\{s-y,0\})^2$

Aussi

$\mathbb{P}(S>s)=\int_{0}^1 \big( 1-(\max\{s-y,0\})^2 \big)dy$

Première chose, on note que
$\mathbb{P}(S>1)=\int_{0}^1 \big( 1-(\max\{1-y,0\})^2 \big)dy=\int_{0}^1 \big( 1-(1-y)^2 \big)dy$
i.e.
$\mathbb{P}(S>1)=\left[y+\frac{(1-y)^3}{3}\right]_0^1=1+0-\frac{1}{3}=\frac{2}{3}$
J’avais alors suggéré d’utiliser (en utilisant la propriété que pour une variable positive, l’espérance est l’intégrale de la fonction de survie) que
$\mathbb{E}(R)=\int_0^1 \mathbb{P}(R>s)ds$
Or cette expression n’est valide que pour une variable strictement positive ! Or ici, il y a de forte chance (1/3) que le remboursement soit nul. Le plus simple est alors de faire des calculs classiques d’espérance. Notons que
$\mathbb{P}(S>1+s)=\left[y+\frac{(1+s-y)^3}{3} \right]_0^1=1+\frac{s^3}{3}-\frac{(s+1)^3}{3}$
donc un peu de calcul de dérivées permet d’écrire $f_S(s)=s(2-s)$ pour $s\in[1,2]$. On peut alors calculer l’espérance du remboursement
$\mathbb{E}(R)=\int_1^2(s-1)f_S(s)ds=\int_1^2(s-1)s(2-s)ds$
soit
$\mathbb{E}(R)=\int_1^2(s^3+3s^2-2s)ds=\left[\frac{s^4}{4}+s^3-s^2\right]_1^2=\frac{1}{4}$
Youpi, avec plein de calculs, on arrive au même résultat…

Faisons un peu de simulations pour vérifier,

> n=1000000
> set.seed(1)
> X=sqrt(runif(n))
> Y=runif(n)

On peut vérifier par exemple la densité de $S$ sur $[1,2]$,

> u=seq(1,2,by=.01)
> hist(X+Y,probability=TRUE,col="light blue")
> lines(u,u*(2-u),lwd=2,col="red")

Si on regarde maintenant la prime pure, on retrouve le résultat obtenu par deux calculs différents,

> R=(X+Y-1)*(X+Y>1)
> mean(R)
[1] 0.2497528

# Somewhere else, part 17

This week, the most important news is undoubtly

I also read a very interesting post

among many others

with several articles early this week on the Nobel price (in Economics)

avec plusieurs articles en francais

mais aussi

Did I miss something ?

# Actuariat 1, ACT2121, troisième cours

Pour la troisième série d’exercices, les exercices sont en ligne ici.

Je rappelle que lundi 29 octobre aura lieu le premier examen (qui compte pour 33,33% de la note finale, comme annoncé dans le plan de cours, en ligne ici). Il y aura 30 questions, en français, du même type que celles vues toutes les semaines. La correction se fera suivant le principe de l’examen P de la SOA

A candidate’s score will be based on the number of questions answered correctly. No credit will be given for omitted answers, and no credit will be lost for wrong answers; therefore, a candidate should answer all questions, even if the candidate needs to guess.

# Fractals and Kronecker product

A few years ago, I went to listen to Roger Nelsen who was giving a talk about copulas with fractal support. Roger is amazing when he gives a talk (I am also a huge fan of his books, and articles), and I really wanted to play with that concept (that he did publish later on, with Gregory Fredricks and José Antonio Rodriguez-Lallena). I did mention that idea in a paper, writen with Alessandro Juri, just to mention some cases where deriving fixed point theorems is not that simple (since the limit may not exist).

The idea in the initial article was to start with something quite simple, a the so-called transformation matrix, e.g.

$T=\frac{1}{8}\left(\begin{matrix}1& 0 & 1 \\ 0 & 4 & 0 \\ 1 & 0&1\end{matrix}\right)$
Here, in all areas with mass, we spread it uniformly (say), i.e. the support of $T(C^\perp)$ is the one below, i.e. $1/8$th of the mass is located in each corner, and $1/2$ is in the center. So if we spread the mass to have a copula (with uniform margin,)we have to consider squares on intervals $[0,1/4]$$[1/4,3/4]$ and $[3/4,1]$,

Then the idea, then, is to consider $T^2=\otimes^2T$, where  $\otimes^2T$ is the tensor product (also called Kronecker product) of $T$ with itself. Here, the support of $T^2(C^\perp)$ is

Then, consider $T^3=\otimes^3T$, where $\otimes^3T$ is the tensor product of $T$ with itself, three times. And the support of $T^3(C^\perp)$ is

Etc. Here, it is computationally extremely simple to do it, using this Kronecker product. Recall that if $\mathbf{A}=(a_{i,j})$, then

$\mathbf{A}\otimes\mathbf{B} = \begin{pmatrix} a_{11} \mathbf{B} & \cdots & a_{1n}\mathbf{B} \\ \vdots & \ddots & \vdots \\ a_{m1} \mathbf{B} & \cdots & a_{mn} \mathbf{B} \end{pmatrix}$

So, we need a transformation matrix: consider the following $4\times4$ matrix,

> k=4
> M=matrix(c(1,0,0,1,
+            0,1,1,0,
+            0,1,1,0,
+            1,0,0,1),k,k)
> M
[,1] [,2] [,3] [,4]
[1,]    1    0    0    1
[2,]    0    1    1    0
[3,]    0    1    1    0
[4,]    1    0    0    1

Once we have it, we just consider the Kronecker product of this matrix with itself, which yields a $4^2\times4^2$ matrix,

> N=kronecker(M,M)
> N[,1:4]
[,1]  [,2] [,3] [,4]
[1,]     1    0    0    1
[2,]     0    1    1    0
[3,]     0    1    1    0
[4,]     1    0    0    1
[5,]     0    0    0    0
[6,]     0    0    0    0
[7,]     0    0    0    0
[8,]     0    0    0    0
[9,]     0    0    0    0
[10,]    0    0    0    0
[11,]    0    0    0    0
[12,]    0    0    0    0
[13,]    1    0    0    1
[14,]    0    1    1    0
[15,]    0    1    1    0
[16,]    1    0    0    1

And then, we continue,

> for(s in 1:3){N=kronecker(N,M)}

After only a couple of loops, we have a $4^5\times4^5$ matrix. And we can plot it simply to visualize the support,

> image(N,col=c("white","blue"))

As we zoom in, we can visualize this fractal property,

# ACT2121, éléments de correction

La plupart des calculs pouvant se faire sans calculs (trop) complexes avec une calculatrice, je vais revenir sur un exercice (l’exercice 11 de la seconde série) pour proposer une correction.

Le petit Nestor collectionne les cartes de joueurs de Baseball dans les paquets de gommes à mâcher. Il y a en tout 20 cartes différentes (réparties aléatoirement, une par paquet). Combien de paquets de gommes Nestor devrait-il s’attendre à avoir à acheter pour obtenir la collection complète ?

La bonne stratégie était d’ordonner les cartes par ordre d’apparition, pour la première fois, et de modéliser le nombre de paquets entre deux premières apparitions, comme indiqué sur le dessin ci-dessous

Si $N$ désigne le nombre total de paquets à acheter, on note $M_i$ le nombre de paquets à acheter, entre l’apparition de la $i-1$ème nouvelle carte, et la $i$ème (avec la convention que $M_1$vaut 1, i.e. à l’achat du premier paquet, on a notre première carte). On notera que
$M_i\sim\mathcal{G}((20-(i-1))/20)$
et à partir de là, les calculs sont simples, puisque

$\mathbb{E}(N)=\mathbb{E}(M_1+M_2+\cdots+M_{20})$

soit (par linéarité de l’espérance)

$\mathbb{E}(N)=\mathbb{E}(M_1)+\mathbb{E}(M_2)+\cdots+\mathbb{E}(M_{20})$

i.e.

$\mathbb{E}(N)=\left( 1+ \frac{20}{19} + \frac{20}{18} + \frac{20}{17} + ... + \frac{20}{2}+\frac{20}{1} \right)$

ou encore

$\mathbb{E}(N)=20\left( \frac{1}{20}+ \frac{1}{19} + \frac{1}{18} + \frac{1}{17} + ... + \frac{1}{2}+1 \right)$

La bonne réponse était alors (il faut sommer les vingt termes)

> sum(20/(1:20)) [1] 71.95479

Mais cette sommation de vingt termes n’est pas triviale, aussi, en cours, j’avais suggéré que

$1+ \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... + \frac{1}{n} \sim \ln(n)$
> log(20)*20 [1] 59.91465

qui diffère du résultat numérique car il manquait la constante d’Euler

$\gamma = \lim_{n \rightarrow \infty } \left( 1+ \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... + \frac{1}{n} - \ln(n) \right)$

i.e.

$1+ \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... + \frac{1}{n} \sim \ln(n)+\gamma$

Numériquement, on obtient, en prenant pour $\gamma$ 0.57721 (cf Google)

> (log(20)+ 0.57721)*20 [1] 71.45885

# Somewhere else, part 16

A very interesting post discovered recently