L’examen intra avait lieu ce matin. L’examen est en ligne ici, avec les annexes, en ligne là. Des éléments de correction seront mis en ligne très bientôt.
Tag Archives: regression
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
Can we say something if we look (only) at that output ? The intercept is significatively non-null, as well as the slope, the 
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 
‘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 
‘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 ‘s and the ‘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),
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
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
This time, there is an other kind of outlier, in ‘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,
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
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.
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
> 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
(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).
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 dans le modèle multiple est légèrement différent de
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
et à droite celle de
)
 |
 |
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 et à droite les
, 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 -value: plus la
-value est faible (rejet de l’hypothèse d’indépendance) plus la couleur sera vive (et rouge), et plus la
-value est forte (acceptation de l’hypothèse d’indépendance) plus la couleur sera terne (et jaune, voire blanche),
 |
 |
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 -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.
Examen final ACT6420
Chose promise, chose due: les deux bases de données qui seront utilisées pour l’examen sont désormais en ligne. 10 questions porteront sur l’analyse de la série temporelle suivante
serie=read.table("http://freakonometrics.free.fr/exam-ts.txt") X=serie$x XTS=serie=ts(X,start=c(1999,9),frequency=12)
(qui peut être modélisée soit comme un objet série temporel, soit comme un vecteur numérique).
Enfin, 10 questions porteront sur la base de données suivante
examen=read.table( "http://freakonometrics.blog.free.fr/public/ data/basket-exam-v2.csv", header=TRUE,sep=";")
(si l’importation ne marche pas, veillez à ce que l’adresse du fichier tienne sur une ligne). Cette base devait servir initialement pour l’examen intra. Elle contient des matchs de baskets universitaires (aux États-Unis). Les variables sont les suivantes
- EQUIPE le nom de l’équipe qui sert de référence
- HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l’EQUIPE perdait)
- WINNER vaut 1 si l’équipe a gagné le match, 0 sinon,
- TOTALPOINTSDIFF le nombre de points de différence à la fin du match
- SECONDHALFPOINTS est le nombre de points marqués par l’EQUIPE pendant la seconde mi-temps
- NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
- HOMETEAM vaut 0 si l’équipe joue à l’extérieur, 1 si elle joue à domicile
- HALFWINNER vaut 1 si l’EQUIPE gagnait (supérieur ou nul) à la mi-temps
- HALFLOOSER vaut 1 si l’EQUIPE perdait (inférieur ou nul) à la mi-temps
- SECONDHALFPOINTSOTHER est le nombre de points marqués par l’adversaire au cours de la seconde mi-temps
La variable d’intérêt sera le nombre de points d’écart à la fin du match.
Enfin, 20 questions supplémentaires porteront sur la compréhension générale du cours. Les calculatrices seront autorisées. Maintenant, pour tous ceux qui ont encore du mal avec R, je recommande
- “R pour les débutants”d’Emmanuel Paradis (PDF)
- “Brise Glace-R” d’Andrew Robinson er Arnaud Schloesing (PDF)
- “Introduction à la programmation en R” de Vincent Goulet (PDF)
(les questions d’analyse de sortie ne porteront pas sur la programmation en R, mais davantage sur l’analyse critique des sorties, et leur interprétation).
Short versus long papers, in academic journals
This Monday, during my talk on quantile regressions (at the Montreal R-meeting), we’ve seen how those nice graphs could be interpreted, with the evolution of the slope of the linear regression, as a function of the probability level. One illustration was on large hurricanes, from Elsner, Kossin & Jagger (2008). The other one was on birthweight, from Abrevaya (2001).
It is also to illustrate that technique to academic publication, e.g. the length of papers, over time. Actually, the data we can extract from Scopus are quite similar to the ones uses on hurricanes. For several journals, it is possible to look at the length of articles. Since Scopus is quite expensive ($60,000 per year for the campus, as far as remember, so I can imagine the penalty I might have to pay for sharing such a dataset)
base=read.table("/home/scopus.csv", header=TRUE,sep=",") pages=base$Page.end-base$Page.start year=base$Year
Again, a first idea can be to look at boxplots, and regression on (nonparametric) quantiles, here for Econometrica,
boxplot(pages~as.factor(year),col="light blue") Q=function(p=.9) as.vector(by(pages,as.factor(year), function(x) quantile(x,p))) u=1:16 points(u,Q(p),pch=19,col="blue") abline(lm(Q(p)~u,weights=table(year)),lwd=2,col="blue")
Consider now (as in the slides in the previous post) a quantile regression (instead of a regression on quantiles), for instance in the Annals of Probability,
library(quantreg) u=seq(.05,.95,by=.01) coefstd=function(u) summary(rq(pages~year, tau=u))$coefficients[,2] coefest=function(u) summary(rq(pages~year, tau=u))$coefficients[,1] CS=Vectorize(coefstd)(u) CE=Vectorize(coefest)(u) k=2 plot(u,CE[k,],ylim=c(min(CE[k,]-2*CS[k,]), max(CE[k,]+2*CS[k,]))) polygon(c(u,rev(u)),c(CE[k,]+1.96*CS[k,], rev(CE[k,]-1.96*CS[k,])), col="light green",border=NA) lines(u,CE[k,],lwd=2,col="red") abline(h=0)
We have the following slope, for the year, as a function of the probability level,
The slope is always positive, so size of papers is increasing with time, short and long papers. But the influence of time is much larger for long paper than short one: for short papers (lower decile) every year, the size keeps increasing, with one more page every three years. For long paper (upper decile), it is two more pages every three years.
If we look now at the Annals of Statistics, we have
and for the evolution of the slope of the quantile regression,
Again the impact is positive: papers are longer in 2010 than 15 years ago. But the trend is the reverse: short papers (lower decile) are much longer, almost one more page every year, with long paper increase only by one more page every two years… Initially, I want to run such a study on a much longer term, with quantile regressions and splines to see when there might have been a change, both in lower and upper tails. Unfortunately, as suggested by some colleagues, there might have been some changes in the format of the journal (columns, margins, fonts, etc). That’s a shame, because I rediscover nice short papers of 5-10 pages published 20 or 30 years ago. They are nice to read (and also potentially interesting for a post on the blog). 5 pages, that’s perfect, but 40 pages, that’s way too long. I wonder if I am the only one having this feeling, missing those short but extremely interesting papers….
Régression médiane et géométrie
Suite à mon exposé d’hier, Pierre me demandait d’argumenter un point que j’avais évoqué oralement sans le justifier: “ une régression médiane passe forcément par deux points du nuage “. Par exemple,
x=c(1,2,3) y=c(3,7,8) plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10)) library(quantreg) abline(rq(y~x,tau=.5),col="red")
Essayons de le justifier… de manière un peu heuristique. Commençons un peu au hasard, avec une droite qui passe “entre” les points (on admettra que la régression médiane sépare l’espace en deux, avec la moitié des points en dessous, la moitié au dessus, à un près pour des histoires de parité).
c’est joli, mais on peut sûrement faire mieux. En particulier, on va chercher ici à minimiser la somme des valeurs absolue des erreurs (c’est le principe de la régression médiane). Essayons de translater la courbe, vers le haut, ou vers le bas,
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10)) abline(a=-1,b=3.2,col="blue",lwd=2) for(i in seq(-2,3,by=.5)){ abline(a=i,b=3.2,col="blue",lty=2)}
Si on regarde ce que vaut la somme des valeurs absolues des erreurs, on a
d=function(h) sum(abs(y-(h+3.2*x))) H=seq(-4,6,by=.01) D=Vectorize(d)(H) plot(H,D,type="l")
Formellement, l’optimum est ici
> optimize(d,lower=-5,upper=5) $minimum [1] -0.2000091 $objective [1] 2.200009
Bref, retenons cette courbe, et notons qu’elle passe par un des points,
Et c’est assez normal. On commençait avec deux points au dessus, et un en dessous. Soit la somme des valeurs absolues des erreurs . Si on translate de
, on passe de
à
(tant que l’on a toujours deux points au dessus, car celui en dessous sera toujours en dessous). Si on translate de
,
(là encore tant qu’il n’y a qu’un point en dessous). Bref, on a intérêt à translater vers le haut. Si on dépasse le premier point rencontré, on se retrouve dans la situation inverse, avec un point au dessus et deux au dessous, et on a intérêt à redescendre. Bref, la translation optimale revient a s’arrêter dès qu’on croise un point. Autrement dit, la régression passe forcément par un point. Au moins.
Maintenant, essayons de faire pivoter la courbe autour de ce point, là encore afin de minimiser la somme des valeurs absolues des erreurs,
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10)) abline(a=optimize(d,lower=-5,upper=5)$minimum,b=3.2, col="blue",lwd=2) points(x[1],y[1],cex=1.8,lwd=2,col="blue") for(i in seq(-1,5,by=.25)){ abline(a=(y[1]-i*x[1]),b=i,col="blue",lty=2)}
La distance est alors, en fonction de la pente de la droite
d2=function(h) sum(abs(y-((y[1]-h*x[1])+h*x))) H=seq(-4,6,by=.01) D=Vectorize(d2)(H) plot(H,D,type="l")
Là encore on peut formaliser un peu
> optimize(d2,lower=-5,upper=5) $minimum [1] 2.500018 $objective [1] 1.500018
Et si on regarde cette dernière droite, on passe par deux points,
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10)) h=optimize(d2,lower=-5,upper=5)$minimum abline(a=(y[1]-h*x[1]), b=h,col="purple",lwd=2)
Pour comprendre pourquoi, comparons les deux cas,
- en faisant un pivot vers le bas, pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue
- en faisant un pivot vers le haut, là encore pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue, mais en sens inverse par rapport au cas précédent,
Et pour faire simple, dans le premier cas, le gain (sur la baisse) compense la perte, alors que dans le second cas, c’est le contraire. Bref, on a intérêt à pivoter vers le bas, ou vers le point à droite. Jusqu’à l’atteindre. Et optimalement, on passera par ce point. Bref, on passera par deux points. Et je laisse les plus courageux regarder avec plus de trois points, mais c’est toujours pareil…
Talk on quantiles at the R Montreal group
This afternoon, I will be giving a two-hour talk at McGill on quantiles, quantile regressions, confidence regions, bagplots and outliers. Before defining (properly) quantile regressions, we will mention regression on (local) quantiles, as on the graph below, on hurricanes,
In order to illustrate quantile regression, consider the following natality database,
base=read.table( "http://freakonometrics.free.fr/natality2005.txt", header=TRUE,sep=";")
We can use it produce those nice graphs we can find in several papers, modeling weight of newborns,
u=seq(.05,.95,by=.01) coefstd=function(u) summary(rq(WEIGHT~SEX+ SMOKER+WEIGHTGAIN+BIRTHRECORD+AGE+ BLACKM+ BLACKF+COLLEGE,data=base,tau=u))$coefficients[,2] coefest=function(u) summary(rq(WEIGHT~SEX+ SMOKER+WEIGHTGAIN+BIRTHRECORD+AGE+ BLACKM+ BLACKF+COLLEGE,data=base,tau=u))$coefficients[,1] CS=Vectorize(coefstd)(u) CE=Vectorize(coefest)(u)
 |
 |
 |
 |
 |
 |
 |
 |
The slides can be downloaded on the blog, as well as the R-code.
Visualization in regression analysis
Visualization is a key to success in regression analysis. This is one of the (many) reasons I am also suspicious when I read an article with a quantitative (econometric) analysis without any graph. Consider for instance the following dataset, obtained from http://data.worldbank.org/, with, for each country, the GDP per capita (in some common currency) and the infant mortality rate (deaths before the age of 5),
> library(gdata) > XLS1=read.xls("http://api.worldbank.org/datafiles /NY.GDP.PCAP.PP.CD_Indicator_MetaData_en_EXCEL.xls", sheet = 1) > data1=XLS1[-(1:28),c("Country.Name","Country.Code","X2010")] > names(data1)[3]="GDP" > XLS2=read.xls("http://api.worldbank.org/datafiles /SH.DYN.MORT_Indicator_MetaData_en_EXCEL.xls", sheet = 1) > data2=XLS2[-(1:28),c("Country.Code","X2010")] > names(data2)[2]="MORTALITY" > data=merge(data1,data2) > head(data) Country.Code Country.Name GDP MORTALITY 1 ABW Aruba NA NA 2 AFG Afghanistan 1207.278 149.2 3 AGO Angola 6119.930 160.5 4 ALB Albania 8817.009 18.4 5 AND Andorra NA 3.8 6 ARE United Arab Emirates 47215.315 7.1
If we estimate a simple linear regression – 
– we get
> regBB=lm(MORTALITY~GDP,data=data) > summary(regBB) Call: lm(formula = MORTALITY ~ GDP, data = data) Residuals: Min 1Q Median 3Q Max -45.24 -29.58 -12.12 16.19 115.83 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 67.1008781 4.1577411 16.139 < 2e-16 *** GDP -0.0017887 0.0002161 -8.278 3.83e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 39.99 on 167 degrees of freedom (47 observations deleted due to missingness) Multiple R-squared: 0.2909, Adjusted R-squared: 0.2867 F-statistic: 68.53 on 1 and 167 DF, p-value: 3.834e-14
We can look at the scatter plot, including the linear regression line, and some confidence bounds,
> plot(data$GDP,data$MORTALITY,xlab="GDP per capita", + ylab="Mortality rate (under 5)",cex=.5) > text(data$GDP,data$MORTALITY,data$Country.Name,pos=3) > x=seq(-10000,100000,length=101) > y=predict(regBB,newdata=data.frame(GDP=x), + interval="prediction",level = 0.9) > lines(x,y[,1],col="red") > lines(x,y[,2],col="red",lty=2) > lines(x,y[,3],col="red",lty=2)
We should be able to do a better job here. For instance, if we look at the Box-Cox profile likelihood,
> boxcox(regBB)
it looks like taking the logarithm of the mortality rate should be better, i.e. 
or 
:
> regLB=lm(log(MORTALITY)~GDP,data=data) > summary(regLB) Call: lm(formula = log(MORTALITY) ~ GDP, data = data) Residuals: Min 1Q Median 3Q Max -1.3035 -0.5837 -0.1138 0.5597 3.0583 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.989e+00 7.970e-02 50.05 <2e-16 *** GDP -6.487e-05 4.142e-06 -15.66 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7666 on 167 degrees of freedom (47 observations deleted due to missingness) Multiple R-squared: 0.5949, Adjusted R-squared: 0.5925 F-statistic: 245.3 on 1 and 167 DF, p-value: < 2.2e-16 > plot(data$GDP,data$MORTALITY,xlab="GDP per capita", + ylab="Mortality rate (under 5) log scale",cex=.5,log="y") > text(data$GDP,data$MORTALITY,data$Country.Name) > x=seq(300,100000,length=101) > y=exp(predict(regLB,newdata=data.frame(GDP=x)))* + exp(summary(regLB)$sigma^2/2) > lines(x,y,col="red") > y=qlnorm(.95, meanlog=predict(regLB,newdata=data.frame(GDP=x)), + sdlog=summary(regLB)$sigma^2) > lines(x,y,col="red",lty=2) > y=qlnorm(.05, meanlog=predict(regLB,newdata=data.frame(GDP=x)), + sdlog=summary(regLB)$sigma^2) > lines(x,y,col="red",lty=2)
on the log scale or
> plot(data$GDP,data$MORTALITY,xlab="GDP per capita", + ylab="Mortality rate (under 5) log scale",cex=.5)
on the standard scale. Here we use quantiles of the log-normal distribution to derive confidence intervals.
But why shouldn’t we take also the logarithm of the GDP ? We can fit a model 
or equivalently 
.
> regLL=lm(log(MORTALITY)~log(GDP),data=data) > summary(regLL) Call: lm(formula = log(MORTALITY) ~ log(GDP), data = data) Residuals: Min 1Q Median 3Q Max -1.13200 -0.38326 -0.07127 0.26610 3.02212 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.50192 0.31556 33.28 <2e-16 *** log(GDP) -0.83496 0.03548 -23.54 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.5797 on 167 degrees of freedom (47 observations deleted due to missingness) Multiple R-squared: 0.7684, Adjusted R-squared: 0.767 F-statistic: 554 on 1 and 167 DF, p-value: < 2.2e-16 > plot(data$GDP,data$MORTALITY,xlab="GDP per capita ", + ylab="Mortality rate (under 5)",cex=.5,log="xy") > text(data$GDP,data$MORTALITY,data$Country.Name) > x=exp(seq(1,12,by=.1)) > y=exp(predict(regLL,newdata=data.frame(GDP=x)))* + exp(summary(regLL)$sigma^2/2) > lines(x,y,col="red") > y=qlnorm(.95, meanlog=predict(regLL,newdata=data.frame(GDP=x)), + sdlog=summary(regLL)$sigma^2) > lines(x,y,col="red",lty=2) > y=qlnorm(.05, meanlog=predict(regLL,newdata=data.frame(GDP=x)), + sdlog=summary(regLL)$sigma^2) > lines(x,y,col="red",lty=2)
on the log scales or
> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ", + ylab="Mortality rate (under 5)",cex=.5)
on the standard scale. If we compare the last two predictions, we have
with in blue is the log model, and in red is the log-log model (I did not include the first one for obvious reasons).
Examen intra, modèles de régression
L’examen intra comportera deux exercices proches de ceux faits en démonstration, et une analyse de sorties de régression. Cette dernière partie portera sur la base suivante (ou peut être une extraction de cette base)
> examen=read.table( + "http://freakonometrics.blog.free.fr/public/ data/basket-exam-v2.csv", + header=TRUE,sep=";") > tail(examen) EQUIPE halfdiff winner 9917 MARYLAND -1 0 9918 RHODE ISLAND -15 1 9919 ROBERT MORRIS 4 1 9920 NORTH CAROLINA-ASHEVILLE 11 1 9921 BROWN -2 1 9922 MOUNT ST. MARYS 1 0 totalpointsdiff secondhalfpoints 9917 -10 26 9918 1 44 9919 6 43 9920 19 39 9921 6 48 9922 -9 26 nodiff_pts_1000_left hometeam halfwinner 9917 9 1 0 9918 18 1 0 9919 17 1 1 9920 20 1 1 9921 16 1 0 9922 13 1 1 halflooser secondhalfpointsother 9917 1 35 9918 1 28 9919 0 41 9920 0 31 9921 1 40 9922 0 36
qui contient des matchs de baskets universitaires.
- EQUIPE le nom de l’équipe qui sert de référence
- HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l’EQUIPE perdait)
- WINNER vaut 1 si l’équipe a gagné le match, 0 sinon,
- TOTALPOINTSDIFF le nombre de points de différence à la fin du match
- SECONDHALFPOINTS est le nombre de points marqués par l’EQUIPE pendant la seconde mi-temps
- NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
- HOMETEAM vaut 0 si l’équipe joue à l’extérieur, 1 si elle joue à domicile
- HALFWINNER vaut 1 si l’EQUIPE gagnait (supérieur ou nul) à la mi-temps
- HALFLOOSER vaut 1 si l’EQUIPE perdait (inférieur ou nul) à la mi-temps
- SECONDHALFPOINTSOTHER est le nombre de points marqués par l’adversaire au cours de la seconde mi-temps
Normalité des estimateurs dans une régression
Cette semaine en cours, on a évoqué la normalité des estimateurs dans une régression linéaire. En supposant la normalité des résidus, on a vu en cours que 
était un estimateur Gaussien. En particulier, chacun des estimateurs est alors Gaussien, au sens où 
, ce qui peut se visualiser sur le graphique suivant (la constante est en abscisse, et la pente en ordonnée), avec un intervalle de confiance à 95%,
En fait, la variance étant inconnue (mais pouvant être estimée), si on remplace la variance des résidus par leur estimateur, la loi de l’estimateur est une loi de Student. De même pour l’estimateur de la pente de la droite de régression, 
Le fait que le couple soit également Gaussien fait que l’on peut construire non plus un intervalle de confiance (on n’est plus en dimension 1) mais une ellipse de confiance. On retient une forme elliptique car c’est la plus petite région dans laquelle on se trouvera avec une probabilité de 95% (comme discuté dans un vieux billet).
Mais pour mieux comprendre cette notion d’ellipse de confiance, le plus simple est d’avoir recours à du rééchantillonnage dans la base de données. En effet, comme on ne dispose que d’un jeu de données, on a un estimateur. Et la discussion sur la loi de notre estimateur est purement théorique. Pour visualiser la loi de notre estimateur, il faudrait des centaines, voire des milliers de bases similaires. Que l’on n’a pas. La solution est alors de tirer au hasard des points de notre échantillon (avec remise), i.e. de faire du bootstrap,
set.seed(1) COEF=matrix(NA,10000,2) for(s in 1:nrow(COEF)){ I=sample(1:nrow(cars),nrow(cars),replace=TRUE) COEF[s,]=lm(dist~speed,data=cars[I,])$coefficients }
Si on regarde ce qui se passe, tirage après tirage, on génère un paquet d’échantillon, et pour chaque échantillon, on ajuste une droite de régression.
Les distributions – sur tous ces échantillons – de la constante et de la pente semblent effectivement Gaussiennes,
hist(COEF[,1],col="light blue",prob=TRUE) u=seq(min(COEF[,1]),max(COEF[,1]),length=500) v=dnorm(u,mean(COEF[,1]),sd(COEF[,1])) lines(u,v,lwd=3,col="red") hist(COEF[,2],col="light blue",prob=TRUE) u=seq(min(COEF[,2]),max(COEF[,2]),length=500) v=dnorm(u,mean(COEF[,2]),sd(COEF[,2])) lines(u,v,lwd=3,col="red")
pour la distribution de l’estimateur de la pente, alors que pour l’estimateur de la constante, on a la distribution suivante,
Si on regarde maintenant la loi jointe,
on retrouve une forme elliptique pour le nuage des points (et une forte corrélation négative entre les deux estimateurs). Mais si on creuse un peu
library(ellipse) reg=lm(dist~speed,data=cars) e=ellipse(reg) plot(e,type="l",lwd=2) polygon(e,col="light blue") points(COEF,cex=.5)
l’ellipse ne coïncide pas (tout à fait) avec celle obtenu théoriquement. Ce qui laisse à penser que l’hypothèse de normalité des résidus est peut-être à revoir..,
Multivariate probit regression using (direct) maximum likelihood estimators
Consider a random pair 
of binary responses, i.e. 
with 
taking values 1 or 2. Assume that probability 
can be function of some covariates.
- The Gaussian vector latent structure
A standard model is based a latent Gaussian structure, i.e. there exists some random vector 
such that 
if 
is lower than a given threshold, and 1 otherwise.
As in standard probit models, assume that
where we can assume that 
is a Gaussian random vector. This assumption can be used to derive the likelihood of a sample 
.
> logV=function(parameter){ + CORRELATION=parameter[1] + BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X)) + z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,])) + sigma=matrix(c(1,CORRELATION,CORRELATION,1),2,2) + a11=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma) + for(i in 2:nrow(z)){a11=c(a11,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))} + a10=pnorm(z[1,1],sd=sqrt(sigma[1,1]))-pmnorm(z[1,],varcov=sigma) + for(i in + 2:nrow(z)){a10=c(a10,pnorm(z[i,1],sd=sqrt(sigma[1,1]))-pmnorm(z[i,],varcov=sigma))} + a01=pnorm(z[1,2],sd=sqrt(sigma[2,2]))-pmnorm(z[1,],varcov=sigma) + for(i in + 2:nrow(z)){a01=c(a01,pnorm(z[i,2],sd=sqrt(sigma[2,2]))-pmnorm(z[i,],varcov=sigma))} + a00=1-a10-a01-a11 + -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) + + 1*log(a01) + + 2*log(a10) + + 3*log(a00) ) + } > OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par
(the code is a bit long since I had trouble working properly with matrices – or more precisely to vectorize my functions – so I used loops… I am sure it is possible to write a better code).
It is possible to generate samples (based on that specific model) to check that we can actually derive proper maximum likelihood estimators,
> library(mnormt) > set.seed(1) > n=1000 > r=0.5 > X1=runif(n) > X2=rnorm(n) > Y1S=1+5*X1 > Y2S=8-5*X1 > RES=rmnorm(n,mean=c(0,0),varcov=matrix(c(1,r,r,1),2,2)) > YS=cbind(Y1S,Y2S)+RES > Y1=(YS[,1]>quantile(YS[,1],.5))*1 > Y2=(YS[,2]>quantile(YS[,2],.5))*1 > base=data.frame(i,Y1,Y2,X1,X2,YS) > head(base) i Y1 Y2 X1 X2 Y1S Y2S 1 1 0 0 0.2655087 0.07730312 3.177587 5.533884 2 2 0 0 0.3721239 -0.29686864 1.935307 5.089524 3 3 1 0 0.5728534 -1.18324224 4.757848 5.172584 4 4 1 0 0.9082078 0.01129269 4.600029 3.878225 5 5 0 1 0.2016819 0.99160104 2.547362 6.743714 6 6 1 0 0.8983897 1.59396745 5.309974 4.421523
(the two columns on the right are latent observations, that cannot be used since theoretically they are unobservable). Note that it is a simple regression, one of the component is here only to bring some noise. First of all, let us look at marginal probit regression
> reg1=glm(Y1~X1+X2,data=base,family=binomial) > reg2=glm(Y2~X1+X2,data=base,family=binomial) > summary(reg1) Call: glm(formula = Y1 ~ X1 + X2, family = binomial, data = base) Deviance Residuals: Min 1Q Median 3Q Max -2.90570 -0.50126 -0.00266 0.49162 2.78256 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.291725 0.267149 -16.065 <2e-16 X1 8.656836 0.510153 16.969 <2e-16 *** X2 0.007375 0.090530 0.081 0.935 --- 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 999 degrees of freedom Residual deviance: 726.48 on 997 degrees of freedom AIC: 732.48 Number of Fisher Scoring iterations: 5 > summary(reg2) Call: glm(formula = Y2 ~ X1 + X2, family = binomial, data = base) Deviance Residuals: Min 1Q Median 3Q Max -2.74682 -0.51814 -0.00001 0.57969 2.58565 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.91709 0.24399 16.054 <2e-16 *** X1 -7.89703 0.46277 -17.065 <2e-16 *** X2 0.18360 0.08758 2.096 0.036 * --- 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 999 degrees of freedom Residual deviance: 777.61 on 997 degrees of freedom AIC: 783.61 Number of Fisher Scoring iterations: 5
Here, the optimization yields,
> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par > OPT[1] [1] 0.5261382 > matrix(OPT[2:7],2,3) [,1] [,2] [,3] [1,] -2.451721 4.908633 0.01600769 [2,] 2.241962 -4.539946 0.10614807
Note that the coefficients we have obtained are almost identical to the ones obtained with R standard procedure,
> library(Zelig) > REG= zelig(list(mu1=Y1~X1+X2, + mu2=Y2~X1+X2, + rho=~1), + model="bprobit",data=base) > summary(REG) Call: zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2, rho = ~1), model = "bprobit", data = base) Pearson Residuals: Min 1Q Median 3Q Max probit(mu1) -10.5442 -0.377243 0.0041803 0.36709 8.60398 probit(mu2) -7.8547 -0.376888 0.0083715 0.42923 5.88264 rhobit(rho) -13.8322 -0.091502 -0.0080544 0.37218 0.85101 Coefficients: Value Std. Error t value (Intercept):1 -2.451699 0.135369 -18.11116 (Intercept):2 2.241964 0.125072 17.92536 (Intercept):3 1.169461 0.189771 6.16249 X1:1 4.908617 0.252683 19.42602 X1:2 -4.539951 0.233632 -19.43203 X2:1 0.015992 0.050443 0.31703 X2:2 0.106154 0.049092 2.16235 Number of linear predictors: 3 Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho) &n bsp; Dispersion Parameter for binom2.rho family: 1 Residual Deviance: 1460.355 on 2993 degrees of freedom Log-likelihood: -730.1774 on 2993 degrees of freedom Number of Iterations: 3 > matrix(coefficients(REG)[c(1:2,4:7)],2,3) [,1] [,2] [,3] [1,] -2.451699 4.908617 0.01599183 [2,] 2.241964 -4.539951 0.10615443
The correlation here is also the same
That procedure works well an can be extended to ordinal responses (not only binary ones, or to three dimensional problems,
logV=function(beta){ BETA=matrix(beta[4:(3+ncol(Y)*ncol(X))],ncol(Y),ncol(X)) z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]),X%*%(BETA[3,])) r12=beta[1] r23=beta[2] r31=beta[3] s1=s2=s3=1 sigma=matrix(c(s1^2,r12*s1*s2,r31*s1*s3, r12*s1*s2,s2^2,r23*s2*s3, r31*s1*s3,r23*s2*s3,s3^2),3,3) sigma1=matrix(c(s2^2,r23*s2*s3, r23*s2*s3,s3^2),2,2) sigma2=matrix(c(s1^2,r31*s1*s3, r31*s1*s3,s3^2),2,2) sigma3=matrix(c(s1^2,r12*s1*s2, r12*s1*s2,s2^2),2,2) a111=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma) for(i in 2:nrow(z)){a111=c(a111,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))} a011=pmnorm(z[1,2:3],varcov=sigma1)-pmnorm(z[1,],varcov=sigma) for(i in 2:nrow(z)){a011=c(a011,pmnorm(z[i,2:3],varcov=sigma1)-pmnorm(z[i,],varcov=sigma))} a101=pmnorm(z[1,c(1,3)],varcov=sigma2)-pmnorm(z[1,],varcov=sigma) for(i in 2:nrow(z)){a101=c(a101,pmnorm(z[i,c(1,3)],varcov=sigma2)-pmnorm(z[i,],varcov=sigma))} a110=pmnorm(z[1,1:2],varcov=sigma3)-pmnorm(z[1,],varcov=sigma) for(i in 2:nrow(z)){a110=c(a110,pmnorm(z[i,1:2],varcov=sigma3)-pmnorm(z[i,],varcov=sigma))} a100=pnorm(z[1,1],sd=s1)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma) for(i in 2:nrow(z)){a100=c(a100,pnorm(z[i,1],sd=s1)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))} a010=pnorm(z[1,2],sd=s2)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(2,3)],varcov=sigma1)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma) for(i in 2:nrow(z)){a010=c(a010,pnorm(z[i,2],sd=s2)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(2,3)],varcov=sigma1)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))} a001=pnorm(z[1,3],sd=s3)-pmnorm(z[1,c(2,3)],varcov=sigma1)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma) for(i in 2:nrow(z)){a001=c(a001,pnorm(z[i,3],sd=s3)-pmnorm(z[i,c(2,3)],varcov=sigma1)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))} a000=1-a111-a011-a101-a110-a001-a010-a100 a111[a111<=0]=1e-50 a110[a110<=0]=1e-50 a101[a101<=0]=1e-50 a011[a011<=0]=1e-50 a100[a100<=0]=1e-50 a010[a010<=0]=1e-50 a001[a001<=0]=1e-50 a000[a000<=0]=1e-50 -sum(((Y[,1]==0)&(Y[,2]==0)&(Y[,3]==0))*log(a111) + 4*log(a011) + 5*log(a101) + 6*log(a110) + 7*log(a001) + 8*log(a010) + 9*log(a100) + 10*log(a000) ) }
A strong assumption in that bivariate model is that residuals have a Gaussian structure. It is possible to change that assumption
- marginally: for instance if we use a logistic cumulative distribution function, then we will have a bivariate logit regression
- in terms of dependence structure: it is possible to consider another copula than the gaussian one, e.g. Gumbel’s copula (also called the bivariate logistic copula), or Clayton’s
Here, the following code can be used to extend the model to non Gaussian structures,
> F=function(x,r){pmnorm(x,rep(0,length(x)), + varcov=matrix(c(1,r,r,1),2,2))} > Fx=function(x1){F(c(x1,1e40),0)} > Fy=function(x2){Fx(x2)} > > logVgen=function(parameter){ + CORRELATION=parameter[1] + BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X)) + z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,])) + a11=F(z[1,],r=CORRELATION) + for(i in 2:nrow(z)){a11=c(a11,F(z[i,],r=CORRELATION))} + a10=Fx(z[1,1])-F(z[1,],r=CORRELATION) + for(i in 2:nrow(z)){a10=c(a10,Fx(z[i,1])-F(z[i,],r=CORRELATION))} + a01=Fy(z[1,2])-F(z[1,],r=CORRELATION) + for(i in 2:nrow(z)){a01=c(a01,Fy(z[i,2])-F(z[i,],r=CORRELATION))} + a00=1-a10-a01-a11 + -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) + + 11*log(a01) + + 12*log(a10) + + 13*log(a00) ) + } > > beta0=c(0,1,1,1,1,1,1) > (OPT=optim(fn=logVgen,par=beta0,method="BFGS")$par) [1] 0.52613820 -2.45172059 2.24196154 4.90863292 -4.53994592 0.01600769 [7] 0.10614807 There were 23 warnings (use warnings() to see them)
E.g.
> library(copula) > F=function(x,r){pcopula(pnorm(x), claytonCopula(2, r))} > Fx=function(x1){F(c(x1,1e40),0) } > Fy=function(x2){Fx(x2)}
- An application to school tests
Consider the following dataset,
hsb2=read.table("http://freakonometrics.free.fr/hsb2.csv", header=TRUE, sep=",") math_male=hsb2$math[female==0] write_male=hsb2$write[female==0] math_female=hsb2$math[female==1] write_female=hsb2$write[female==1] plot(math_female, write_female, type="p", pch=19,col="red",xlab="maths",ylab="writing",cex=.8) points(math_male, write_male, cex=1.2, col="blue")
with here maths versus writing, with girls in red and boys in blue, where variables here are
female :
0: male
1: female
race :
1: hispanic
2: asian
3: african-amer
4: white
ses :
1: low
2: middle
3: high
schtyp : type of school
1: public
2: private
prog : type of program
1: general
2: academic
3: vocation
read : reading score
write : writing score
math : math score
science : science score
socst : social studies score
We can try to understand correlation between math and writing skills. Covariates can be the sex of the child, and his reading skills. The question will then be: are good students in maths and writing simply students that can read well ?
Here the code is simply
> W=hsb2$write>=50 > M=hsb2$math>=50 > base=data.frame(Y1=W,Y2=M, + X1=hsb2$female,X2=hsb2$read) > > library(Zelig) > REG= zelig(list(mu1=Y1~X1+X2, + mu2=Y2~X1+X2, + rho=~1), + model="bprobit",data=base) > summary(REG) Call: zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2, rho = ~1), model = "bprobit", data = base) Pearson Residuals: Min 1Q Median 3Q Max probit(mu1) -4.7518 -0.502594 0.15038 0.53038 1.8592 probit(mu2) -3.4243 -0.653537 0.23673 0.67011 2.6072 rhobit(rho) -4.9821 0.010481 0.13500 0.40776 2.9171 Coefficients: Value Std. Error t value (Intercept):1 -5.484711 0.787101 -6.96825 (Intercept):2 -4.061384 0.633781 -6.40818 (Intercept):3 1.332187 0.322175 4.13497 X1:1 1.125924 0.233550 4.82092 X1:2 0.167258 0.202498 0.82598 X2:1 0.103997 0.014662 7.09286 X2:2 0.082739 0.012026 6.88017 Number of linear predictors: 3 Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho) Dispersion Parameter for binom2.rho family: 1 Residual Deviance: 364.51 on 593 degrees of freedom Log-likelihood: -182.255 on 593 degrees of freedom Number of Iterations: 3 > (exp(summary(REG)@coef3[3])-1)/(exp( summary(REG)@coef3[3])+1) [1] 0.5824045
with a remaining correlation among residuals of 0.58. So with only the sex of the student, and his or her reading skill, we cannot explain the correlation between maths and writing skills. With our previous code, we have here
> beta0=c((exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1), + summary(REG)@coef3[c(1:2,4:7),1]) > beta0 (Intercept):1 (Intercept):2 X1:1 X1:2 0.58240446 -5.48471133 -4.06138412 1.12592427 0.16725842 X2:1 X2:2 0.10399668 0.08273879 > (OPT=optim(fn=logV,par=beta0,method="BFGS")$par) (Intercept):1 (Intercept):2 X1:1 X1:2 0.5824045 -5.4847113 -4.0613841 1.1259243 0.1672584 X2:1 X2:2 0.1039967 0.0827388
i.e. we obtain (almost) exactly the same estimators. But here I have used as starting values for the optimization procedure the estimators given by R. If we change them, hopefully we have a robust maximum likelihood estimator,
> (OPT=optim(fn=logV,par=beta0/2,method="BFGS")$par) (Intercept):1 (Intercept):2 X1:1 X1:2 0.58233360 -5.49428984 -4.06839571 1.12696594 0.16760347 X2:1 X2:2 0.10417767 0.08287409 There were 12 warnings (use warnings() to see them)
So once again, it is possible to optimize numerically a likelihood function, and it works.
- Y[,1]==0)&(Y[,2]==1 [↩]
- Y[,1]==1)&(Y[,2]==0 [↩]
- Y[,1]==0)&(Y[,2]==0 [↩]
- Y[,1]==1)&(Y[,2]==0)&(Y[,3]==0 [↩]
- Y[,1]==0)&(Y[,2]==1)&(Y[,3]==0 [↩]
- Y[,1]==0)&(Y[,2]==0)&(Y[,3]==1 [↩]
- Y[,1]==1)&(Y[,2]==1)&(Y[,3]==0 [↩]
- Y[,1]==1)&(Y[,2]==0)&(Y[,3]==1 [↩]
- Y[,1]==0)&(Y[,2]==1)&(Y[,3]==1 [↩]
- Y[,1]==1)&(Y[,2]==1)&(Y[,3]==1 [↩]
- Y[,1]==0)&(Y[,2]==1 [↩]
- Y[,1]==1)&(Y[,2]==0 [↩]
- Y[,1]==0)&(Y[,2]==0 [↩]
On a toujours besoin de petits poids (chez soi)
Après avoir mis en ligne le billet d’hier sur les manifestants, j’ai passé la journée à me demander “pourquoi ai-je mis des poids tenant compte d’un effet taille dans ma régression” ?. Sur le moment ça me semblait très naturel, et plus la journée avançait, plus je me mettais à douter…
Bon, pour rappel quand on fait une régression classique, sans poids, et sans constante car on veut un coefficient de proportionnalité, par “moindres carrés“, on résout
ce qui donne
La dernière écriture montre que l’on fait une moyenne pondérée des pentes pour chacune des observations… Graphiquement, ça donne
les traits pointillés permettent de montrer ou se trouve le point moyen, avec la moyenne sur toutes les villes pour les syndicats et pour la police, pour les deux dates. Les flèches indiquent les pentes moyennes, c’est à dire le nombre total de manifestants en France selon les syndicats divisé par le nombre total de manifestants en France selon la police. Bref, la droite de régression “classique” n’a pas pour pente le ratio moyen…
Une alternative usuelle est d’introduire des poids, et on cherche à résoudre
dont la solution a la forme
Hier j’avais proposé de prendre des poids correspondant à la taille de la ville (pour intégrer un effet médiatique – les médias ne parlant parfois que des nombres de manifestants dans les très grosses villes), mais on pourrait proposer, sans le même genre d’idée que de prendre 
(et qui est représenté ci-dessous)
Toutefois, l’idée la plus pertinente aurait dû être 
ce qui donne comme estimateur par moindres carrés pondérés
c’est à dire le ratio moyen national… Bref, on est colinéaire avec nos vecteurs tracés depuis le début
Autrement dit, il aurait été plus malin ici d’attribuer des petits poids aux grosses villes pour avoir un coefficient interprétable….
Mais au débat du débat (probablement légitime) sur le choix de la régression, ce qui est troublant c’est que plus les poids donnés aux grosses villes sont importants, plus la pente est importante… autrement dit, on modèle quadratique pourrait être plus adapté, ou au moins quelque chose de non linéaire, et de strictement croissant. Sur le graphique suivant, on voit l’estimation par moindres carrées locaux, et par splines,
ce qui pourrait nous conforter dans l’idée que “les grandes villes ont des plus grands trottoirs“(ramené à la taille de la chaussée !), ce qui peut surprendre un peu.. Si on regarde l’estimation des pentes obtenues par ces deux méthodes, on a, pour les régressions locales
(avec en trait horizontal en pointillés les pentes des flèches, i.e. les taux observés nationalement), et par splines
On retrouve l’idée que dans les petites villes, le comptage est plus précis (des deux côtés) avec un facteur de 2 et 3 pour les deux manifestations respectivement. En mars dernier, on avait ensuite un facteur croissant passant de 2 à 4 en fonction de la taille de la manifestation. Et pour la dernière manifestation en revanche, c’est un peu n’importe quoi (si je peux me permettre).
Bref, pour revenir au point qui m’intéressait au départ, on note qu’il ne semble pas y avoir de méthode scientifique permettant de passer du nombre estimé par la police à celui proposé par les syndicats, car sinon il y aurait une relative stabilité entre les manifestations… stabilité que l’on n’a pas (quelle que soit le type de régression considéré…).
Qui a eu l’idée d’élargir les trottoirs en France ?
J’étais tombé l’autre jour sur un joli dessin expliquant pourquoi les policiers et les syndicats n’estimaient jamais le même nombre de manifestants, avec une explication scientifique, et vaguement rationnelle, ici.
C’est effectivement joli… et on retrouve un rapport de 3 dont pas mal de monde semble parler, sauf que si c’était vrai, le rapport entre les estimations de la police et des syndicats devrait être stable d’une manifestation à l’autre, non ?
Si on compare le 19 mars dernier (ici) et le 16 octobre (là), on arrive aux chiffres suivants,
Les traits pointillés sont la droite de régression. La courbe en trait plein et la régression (passant par l’origine car on cherche un rapport de proportionnalité) en prenant en compte un effet taille, i.e. en introduisant des poids différents dans ma régression. J’ai pris la taille des villes, autrement dit, plus la ville est grosse, plus je donne un poids important au rapport… L’idée étant d’intégrer un effet médiatique, même si la pente n’est plus le ratio moyen… Avec cette dernière méthode, j’obtiens des pentes de 4,06 pour mars, et 6,12 pour octobre… ce qui est à un peu différents des rapports nationaux mentionnés sur le site du monde, i.e. (3 millions selon les syndicats et 1,2 selon la police, en mars, contre 2,5 millions contre 825 milliers l’autre jour). Notons qu’on peut vérifier rapidement que ce n’est pas un problème statistique, les valeurs sont sensiblement différentes, si l’on compare les distributions des estimateurs..
Bref, quelque chose semble avoir changé entre ces deux dates…. la taille des trottoirs ?
De l’interprétation d’un effet nonlinéaire en régression
Je vais poursuivre un peu le TD de statistique de l’actuariat 2, où nous avions parlé tarification et régression Poissonnienne.
- Analyser un effet nonlinéaire
Dans un premier temps, nous avions regardé un modèle simple, où l’on voulait voir si l’âge du conducteur expliquait la fréquence de sinistres. Bref, on peut utiliser une régression log-Poisson, qui donne la prédiction suivante, en fonction de l’âge,
Mais l’âge est une variable un peu particulière car c’est une fausse variable continue: seule les valeurs entières apparaissent. Comme on chercher à calculer
notons qu’il est possible de considérer l’estimateur empirique naturel de cette grandeur, i.e.
où
Graphiquement, on obtient la courbe suivante,
Le côté ératique de la série à droite et à à gauche vient de la faible représentation dans la base, comme en atteste l’histogramme de l’âge du conducteur,
Je peux en profiter pour revenir à un commentaire posté ici sur les modèles GAM. En effet, au lieu de regarder la fréquence en fonction de l’âge, on peut s’intéresser à la la différence par rapport à la fréquence globale.
> m=mean(base$nbsin) > Y=predict(REG,newdata=data.frame(age=X,expo=E),type = "response") > lines(X,Y/m-1,col="red")
On obtient alors le graphique suivant (en superposant à l’ajustement obtenu par GAM)
Notons que si l’on ajuste un modèle log-Poisson uniquement sur la tranche d’âge [50-70], on est très proche de ce que propose le modèle GAM, d’où l’interprétation de l’ajustement GAM comme un ajustement local.
- Changer de modèle (et les comparer)
On a choisi un modèle log-Poisson, mais on se doute que d’autres modèles seraient possible. En particulier dans notre base, les assurés ont eu 0 ou 1 sinistre l’année passée. Autrement dit, des lois binomiales peuvent être ajustée,
> REG0 = glm(nbsin~age,base, + family=poisson(link="log"), + offset=expo) > REG1a = glm(nbsin~age,base, + family=binomial(link = "logit"), + offset=expo) > REG1b = glm(nbsin~age,base, + family=binomial(link = "probit"), + offset=expo) > REG1c = glm(nbsin~age,base, + family=gaussian(link = "identity"), + offset=expo)
Le critère AIC est souvent invoqué pour juger de l’adéquation d’un modèle. Il vaut ici
> AIC(REG0) [1] 36470.7 > AIC(REG1a) [1] 36044.45 > AIC(REG1b) [1] 37304.02 > AIC(REG1c) [1] 32549.83
Mais on peut également regarder le critère BIC,
> AIC(REG0, k = log(nrow(base))) [1] 36488.81 > AIC(REG1a, k = log(nrow(base))) [1] 36062.57 > AIC(REG1b, k = log(nrow(base))) [1] 37322.14 > AIC(REG1c, k = log(nrow(base))) [1] 32577
Autrement dit, en utilisant le règle “the smaller, the better” on retiendrait le dernier modèle qui est un modèle Gaussien, manifestement assez mauvais: sur le graphique ci-dessous, la régression log-Poisson est en mauve, presque confondu avec la régression logistique1, en bleu, la régression probit en vert et le modèle linéaire (Gaussien) en rouge,
Bref, ce dernier modèle n’est pas convainquant… Par contre, notons que si l’on compare le modèle log-Poisson et logit, ce dernier conduit à tarifer moins cher les âges extrêmes,
Autrement dit, au lieu de comparer des modèles actuariels sur la base d’outils statistiques, et uniquement de ces outils, il peut être instructif d’essayer de comprendre qui verra sa prime baisser avec tel ou tel modèle… Mais je reviendrais plus longuement, une autre fois, sur les histoires de choix de modèles…..
Propriétés des estimateurs dans une régression
Comme les rappels (je devrais plutôt dire “remise à niveau“) ont été particulièrement rapides, je vais prendre un peu de temps pour revenir sur quelques éléments du cours d’économétrie (d’autant plus que j’ai reçu quelques questions par mail).
- sur le calcul des estimateurs
Je vais reprendre les questions que j’ai reçu, ça sera plus simple, “je comprend pas très bien quelle méthode utilise l’opération lm(Y~X…..) sur R pour trouver les coef” ou encore “je regresse un modèle simple, je sors le summary, et j’obtiens un betaREG 3 étoiles. Toujours sur R, je compare ce betaREG avec le betaCalc = (X’X)-1 * X’Y, il est différent ! Comment expliquer cela ? Est-ce possible ?“.
Commençons par la première question: la fonction lm() utilise la méthode dite des moindres carrés, ce qui revient à calculer
en adoptant une écriture matricielle, i.e. la solution est alors
Cet estimateur est celui qui minimise la somme des carrés des erreurs. Et on peut vérifier numériquement que cet estimateur est bien celui calculé par R,
>data(cars) >X=cbind(rep(1,nrow(cars)),cars$speed); Y=cars$dist > X[1:5,] [,1] [,2] [1,] 1 4 [2,] 1 4 [3,] 1 7 [4,] 1 7 [5,] 1 8 > solve(t(X)%*%X) [,1] [,2] [1,] 0.19310949 -0.011240876 [2,] -0.01124088 0.000729927 > t(X)%*%Y [,1] [1,] 2149 [2,] 38482 > (solve(t(X)%*%X))%*%(t(X)%*%Y) [,1] [1,] -17.579095 [2,] 3.932409
ce qui correspond très précisément à ce que calcule la fonction de R
> lm(dist~speed,cars) Coefficients: (Intercept) speed -17.579 3.932
Et on peut aussi vérifier qu’il minimise bien la somme des carrés des erreurs,
> b0=seq(-30,10,by=3);
> b1=seq(1,7,by=.5)
> SC=matrix(NA,length(b0),length(b1))
> for(i in 1:length(b0)){
> for(j in 1:length(b1)){
+ SC[i,j]=sum((cars$dist-(b0[i]+b1[j]*cars$speed))^2)
}}
> contour(b0,b1,SC)
et je peux certifier que la valeur obtenue est bien le minimum de cette fonction (en fait c’est un résultat théorique d’optimisation).
Ca se généralise bien entendu en dimension plus grande. Alors attention, cette formule permet de calculer un estimateur du vecteur des paramètres, mais la formule est bien sûre fausse si on raisonne composante par composante.
> X=cars$speed; Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
[,1]
[1,] 2.909132
Mais ça c’est un résultat (trivial) d’algèbre linéaire. En effet, le fait d’avoir
ne signifie pas que l’on ait ce genre de formule composante par composante, i.e.
Dit de manière un peu formelle, le produit matriciel n’est pas un produit terme à terme. La traduction économétrique de cette idée est qu’une régression multiple n’est pas une succession de régressions simples. En effet, dans l’équation précédante, le dernier terme correspond à la régression sans la constante, comme on peut le voir ci-dessous
> lm(dist~0+speed,cars) Coefficients: speed 2.909
et le premier est tout simplement la moyenne empirique,
> lm(dist~1,cars)
Coefficients:
(Intercept)
42.98
> mean(cars$dist)
[1] 42.98
> X=rep(1,nrow(cars)); Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
[,1]
[1,] 42.98
- sur les propriétés des estimateurs
Une question était “doit-on vérifier les propriétés de l’estimateur (sans biais, efficace etc…) du summary de la Reg ? pour chaque régression ?”
Bon, la réponse est non car c’est impossible ! Pour répéter tranquillement ce que j’avais dit plusieurs fois oralement: quand on fait un modèle (théorique), on part d’hypothèses (par exemple ici les résidus sont centrées, de variance constante). Sous ces hypothèses on a des propriétés: ce sont des théorèmes, autrement dit de la théorie. Le plus connu étant le théorème de Gauss-Markov sur le modèle linéaire, qui garantie que l’estimateur par moindre carré est BLUE (best linear unbiased estimator).
Autrement dit (si je retraduit ce que dit ce théorème), si les hypothèses sont valides, alors la théorie nous garantie que les estimateurs vérifient des propriétés, dont celle d’être sans biais, par exemple.
On ne
peut pas vérifier qu’un estimateur est sans biais, on peut juste vérifier ex post que les hypothèses sont valides (ou non), ce qui garantit (ou pas) l’absence de biais.
Mais comme on va le voir par la suite, le biais ce n’est pas forcément gênant, ce qui est important, c’est surtout la convergence. Comme le notait Clive Granger “if you can’t even get a consistant estimator, you shouldn’t be in this business“.
- et sur des points plus théoriques
sinon la personne qui avait des soucis numériques à calculer les coefficients me demandait si cela avait un “rapport avec les plim X’E/n et plim (X’X)-1/ n“.
Pour revenir sur ce point, une hypothèse forte dans le modèle linéaire (de base) est que le bruit soit soit un vrai bruit, c’est à dire non corrélé avec la variable explicative, i.e.
En effet, sinon
autrement dit l’estimateur par moindre carré est biasé. La solution la plus naturelle pour éviter ce “problème” est d’utiliser un instrument (ou une variable instrumentale).
Mais surtout, sur cet exemple, on voit qu’on ne pourra avoir la convergence de notre estimateur par moindres carrés que si la corrélation entre le bruit et la variable explicative tend vers 0 asymptotiquement.
Il existe en effet une différence fondamentale entre le biais et la consistance (ou la convergence vers la vraie valeur).
- un estimateur sans biais est généralement convergent (et de manière consistante)
- un estimateur biaisé peut ne pas converger
Prenons un petit exemple pour illustrer ces points. On supposera n1 < n2 < n3par la suite. Comme je l’ai expliqué auparavant, la distribution de l’estimateur est une propriété théorique qui ne peut se voir sur un jeu de données (à moins de faire des simulations mais ça sort du cadre de mon billet d’aujourd’hui). La figure ci-dessous correspond au cas sans biais, et convergent
La figure ci-dessous correspond au cas biaisé, mais convergent, et consistant, au sens où asymptotiquement, l’estimateur sera sans biais.
Enfin, un dernier cas correspondant au cas biaisé, convergent, mais non consistant. Asymptotiquement l’estimateur sera toujours biaisé,
Mais revenons un peu sur la formaliation sous-jacente au modèle linéaire. Les propriétés sur les résidus sont conditionnelles aux variables explicatives, en particulier.
(c’est l’hypothèse d’homoscédasticité et d’indépendance des résidus). Bon, il faut aussi l’absence de relation linéaire entre les variables explicatives, ce qui se traduit parfois comme une absence de collinéarité.
On peut s’intéresser aux propriétés asymptotiques de l’estimateur. Pour cela, il faut peut être rappeler ce que signifie la convergence pour une variable aléatoire (en statistique mathématique, les estimateurs sont vus comme des variables aléatoires). On dira que 
converge en probabilité vers 
, parfois noté
si
Par exemple la loi des grands nombres garantie que, pour un échantillon i.i.d.
si l’espérance des variables est finie, alors
Alors sous les hypothèses mentionnées auparavant, on peut montrer que, dans le modèle linéaire,
On parlera alors d’estimateur convergent. On peut aussi montrer que
Mais formellement,
si
Voilà en gros quelques éléments pour répondre aux questions, mais tout ça sera repris en détails dans le cours d’économétrie 1 qui commencera très bientôt… et surtout d’économétrie 2 qui abordera les aspects dynamiques.