Freakonometrics

To content | To menu | To search

Friday, September 7 2012

That damn R-squared !

Another post about the R-squared coefficient, and about why, after some years teaching econometrics, I still hate when students ask questions about it. Usually, it starts with "I have a _____ R-squared... isn't it too low ?" Please, feel free to fill in the blanks with your favorite (low) number. Say 0.2. To make it simple, there are different answers to that question:

  1. if you don't want to waste time understanding econometrics, I would say something like "Forget about the R-squared, it is useless" (perhaps also "please, think twice about taking that econometrics course")
  2. if you're ready to spend some time to get a better understanding on subtle concepts, I would say "I don't like the R-squared. I might be interesting in some rare cases (you can probably count them on the fingers of one finger), like comparing two models on the same dataset (even so, I would recommend the adjusted one). But usually, its values has no meaning. You can compare 0.2 and 0.3 (and prefer the 0.3 R-squared model, rather than the 0.2 R-squared one), but 0.2 means nothing". Well, not exactly, since it means something, but it is not a measure tjat tells you if you deal with a good or a bad model. Well, again, not exactly, but it is rather difficult to say where bad ends, and where good starts. Actually, it is exactly like the correlation coefficient (well, there is nothing mysterious here since the R-squared can be related to some correlation coefficient, as mentioned in class)
  3. if you want some more advanced advice, I would say "It's complicated..." (and perhaps also "Look in a textbook write by someone more clever than me, you can find hundreds of them in the library !")
  4. if you want me to act like people we've seen recently on TV (during electoral debate), "It's extremely interesting, but before answering your question, let me tell you a story..."
Perhaps that last strategy is the best one, and I should focus on the story. I mean, this is exactly why I have my blog: to tell (nice) stories. With graphs, and math formulas inside. First of all, consider a regression model

so that the R-squared is defined as

Let us generate datasets, and then run regressions, to see what's going on...
For instance, consider 20 observations, with one variable of interest, one explanatory variable, and some low variance noise (to start with)
> set.seed(1)
> n=20
> X=runif(n)
> E=rnorm(n)
> Y=2+5*X+E*.5
> base=data.frame(X,Y)
> reg=lm(Y~X,data=base)
> summary(reg)
 
Call:
lm(formula = Y ~ X, data = base)
 
Residuals:
Min       1Q   Median       3Q      Max
-1.15961 -0.17470  0.08719  0.29409  0.52719
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   2.4706     0.2297   10.76 2.87e-09 ***
X             4.2042     0.3697   11.37 1.19e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.461 on 18 degrees of freedom
Multiple R-squared: 0.8778,	Adjusted R-squared: 0.871
F-statistic: 129.3 on 1 and 18 DF,  p-value: 1.192e-09 
The R-squared is high (almost 0.9). What if the underlying model is exactly the same, but now, the noise has a much higher variance ?
> Y=2+5*X+E*4
> base=data.frame(X,Y)
> reg=lm(Y~X,data=base)
> summary(reg)
 
Call:
lm(formula = Y ~ X, data = base)
 
Residuals:
Min      1Q  Median      3Q     Max
-9.2769 -1.3976  0.6976  2.3527  4.2175
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    5.765      1.837   3.138  0.00569 **
X             -1.367      2.957  -0.462  0.64953
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 3.688 on 18 degrees of freedom
Multiple R-squared: 0.01173,	Adjusted R-squared: -0.04318
F-statistic: 0.2136 on 1 and 18 DF,  p-value: 0.6495 
Now, the R-squared is rather low (around 0.01). Thus, the quality of the regression depends clearly on the variance of the noise. The higher the variance, the lower the R-squared. And usually, there is not much you can do about it ! On the graph below, the noise is changing, from no-noise, to extremely noisy, with the least square regression in blue (and a confidence interval on the prediction)
If we compare with the graph below, one can observe that the quality of the fit depends on the sample size, with now 100 observations (instead of 20),
So far, nothing new (if you remember what was said in class). In those two cases, here is the evolution of the R-squared, as a function of the variance of the noise (more precisely, here, the standard deviation of the noise) 
> S=seq(0,4,by=.2)
> R2=rep(NA,length(S))
> for(s in 1:length(S)){
+ Y=2+5*X+E*S[s]
+ base=data.frame(X,Y)
+ reg=lm(Y~X,data=base)
+ R2[s]=summary(reg)$r.squared}
with 20 obervations in blue, 100 in black. One important point is that in econometrics, we rarely choose the number of observations. If we have only 100 observations, we have to deal with it. Similarly, if observations are quite noisy there is (usually) not much we can do about it. All the more if you don't have any explanatory variable left. Perhaps you migh play (or try to play) with nonlinear effect...

Nevertheless, it looks like some econometricians really care about the R-squared, and cannot imagine looking at a model if the R-squared is lower than - say - 0.4. It is always possible to reach that level ! you just have to add more covariates ! If you have some... And if you don't, it is always possible to use polynomials of a continuous variate. For instance, on the previous example,

> S=seq(1,25,by=1)
> R2=rep(NA,length(S))
> for(s in 1:length(S)){
+ reg=lm(Y~poly(X,degree=s),data=base)
+ R2[s]=summary(reg)$r.squared}
If we plot the R-squared as a function of the degree of the polynomial regression, we have the following graph. Once again, the higher the degree, the more covariates, and the more covariates, the higher the R-squared,
I.e. with 22 degrees, it is possible to reach a 0.4 R-squared. But it might be interesting to the prediction we have with that model,
So, was it worth adding so much polynomial parts ? I mean, 22 is quite a large power... Here, the linear regression was significant, but not great. So what ? The R-squared was small ? so what ? sometimes, there's not much you can do about it... When dealing with individual observations (so called micro-econometrics), the variable of interest might be extremely noisy, and there is not much you can do. So your R-squared can be small, but your regression is perhaps still better than doing nothing... The good thing with a low R-squared is perhaps that it will recall us that we have to remain modest when we build a model. And always be cautious about conclusion. At least, provide a confidence interval...

Thursday, February 2 2012

un mot (rapide) sur le R2 et la sélection de modèle

J'évoquais tout à l'heure dans un rapide billet l'utilisation du http://freakonometrics.blog.free.fr/public/perso5/R2.gif dans une régression sans constante. On va donc évoquer ici uniquement le cas où une constante est prise en compte dans la régression pour parler de la sélection de modèle, et plus particulièrement le http://freakonometrics.blog.free.fr/public/perso5/r2ajuste.gif, souvent appelé http://freakonometrics.blog.free.fr/public/perso5/R2.gif-ajusté. La sélection de modèle sera l'objet du prochain cours. Classiquement, on défini le coefficient d'ajustement comme

http://freakonometrics.blog.free.fr/public/perso5/olssc13.gif

Dans une régression avec une constante, le résidus empiriques sont centrés, et on peut interpréter les termes comme des variances, d'où l'écriture souvent utilisée

http://freakonometrics.blog.free.fr/public/perso5/olssc10.gif

(qui présente l'avantage d'offrir une interprétation simple au coefficient).  En utilisant la formule de décomposition de la variance, c'est à dire le théorème de Pythagore,

http://freakonometrics.blog.free.fr/public/perso5/olssc12.gif

on peut réécrire ce coefficient

http://freakonometrics.blog.free.fr/public/perso5/olssc15.gif

ou encore

http://freakonometrics.blog.free.fr/public/perso5/olssc50.gif

(qui correspond presque à la formule calculée sous R). On en déduit de manière assez simple que l'on aura toujours http://freakonometrics.blog.free.fr/public/perso5/olssc11.gif. Une autre propriété, dès lors que l'on a une constante dans le modèle est que le http://freakonometrics.blog.free.fr/public/perso5/R2.gif peut être relié à la corrélation (au sens de Pearson) entre la prédiction et l'observation,

http://freakonometrics.blog.free.fr/public/perso5/olssc16.gif

C'est d'ailleurs sous cette forme qu'a été introduit, historiquement, ce coefficient (on peut renvoyer aux travaux de Karl Pearson, ou Udny Yule sur le sujet, lors des travaux sur l'hérédité à la fin du XIXème siècle). Mais ce coefficient, s'il pourrait être utilisé comme mesure de la qualité de l'ajustement, ne peut être utilisé pour faire de la sélection de modèle. En effet, en rajoutant des variables, quelles qu'elles soient, améliorera toujours le modèle. C'est ce que l'on arrive à faire ci-dessous en rajoutant juste du bruit dans la régression

> summary(lm(dist~speed,data=cars))
 
Call:
lm(formula = dist ~ speed, data = cars)
 
Residuals:
Min      1Q  Median      3Q     Max
-29.069  -9.525  -2.272   9.215  43.201
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791     6.7584  -2.601   0.0123 *
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared: 0.6511,	Adjusted R-squared: 0.6438
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12
 
> set.seed(1)
> cars$bruit=rnorm(nrow(cars))
> summary(lm(dist~speed+bruit,data=cars))
 
Call:
lm(formula = dist ~ speed + bruit, data = cars)
 
Residuals:
Min      1Q  Median      3Q     Max
-31.869  -9.412  -3.049   7.410  43.769
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -17.905      6.757  -2.650   0.0109 *
speed          3.935      0.415   9.482 1.73e-12 ***
bruit          2.788      2.639   1.056   0.2963
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15.36 on 47 degrees of freedom
Multiple R-squared: 0.6592,	Adjusted R-squared: 0.6447
F-statistic: 45.45 on 2 and 47 DF,  p-value: 1.034e-11
 
> set.seed(1)
> cars$bruit=matrix(rnorm(20*nrow(cars)),nrow(cars),20)
> summary(lm(dist~.,data=cars))
 
Call:
lm(formula = dist ~ ., data = cars)
 
Residuals:
Min      1Q  Median      3Q     Max
-30.147  -5.883  -2.070   7.354  32.096
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -23.8799     8.4775  -2.817  0.00879 **
speed         4.3870     0.5740   7.643 2.52e-08 ***
bruit1       -1.7487     4.1122  -0.425  0.67391
bruit2       -2.7356     3.0890  -0.886  0.38338
bruit3        4.2464     3.5200   1.206  0.23778
bruit4       -4.5160     3.2715  -1.380  0.17838
...
bruit19       1.5769     2.6615   0.592  0.55828
bruit20       1.5107     2.7799   0.543  0.59113
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 15.87 on 28 degrees of freedom
Multiple R-squared: 0.7834,	Adjusted R-squared: 0.621
F-statistic: 4.822 on 21 and 28 DF,  p-value: 7.557e-05 

Plus on rajoute des variables, plus on augmente le http://freakonometrics.blog.free.fr/public/perso5/R2.gif. Un ajustement a été proposé afin de viser un principe de parcimonie. L'ajustement en question peut se retrouver dans Fisher (1924),

même si beaucoup de monde semble l'attribuer à Theil (1961). L'idée derrière est que l'association faite entre les variances et les sommes des carrés est rapide, et peut-être fallacieuse, en pratique. En effet, on utilisait comme estimateurs des variances

http://freakonometrics.blog.free.fr/public/perso5/olssc17.gif

et 

http://freakonometrics.blog.free.fr/public/perso5/olssc18.gif

Or ces estimateurs sont des estimateurs biaisés des variances. Il serait plus judicieux de considérer

http://freakonometrics.blog.free.fr/public/perso5/olssc22.gif

et

http://freakonometrics.blog.free.fr/public/perso5/olssc20.gif

Ceci pousse à poser

http://freakonometrics.blog.free.fr/public/perso5/olssc26.gif

qui est effectivement le coefficient de détermination ajusté. On peut noter que comme on régresse toujours sur quelque chose (au pire une variable), alors http://freakonometrics.blog.free.fr/public/perso5/olssc27.gif (on peut avoir égalité si on a un modèle parfait). On notera aussi qui si l’ajustement est trop mauvais, en particulier

http://freakonometrics.blog.free.fr/public/perso5/olssc28.gif

alors http://freakonometrics.blog.free.fr/public/perso5/r2ajuste.gif est négatif. Le http://freakonometrics.blog.free.fr/public/perso5/r2ajuste.gif est populaire car il permet de pénaliser le http://freakonometrics.blog.free.fr/public/perso5/R2.gif, en prenant en compte la complexité du modèle (ici le nombre de variables explicatives). Pour qu'un modèle soit meilleur avec une variable en plus, il faut que la variable que l'on rajoute apporte vraiment quelque chose. Pour comprendre plus concrètement ce qui se passe, considérons la régression sur http://freakonometrics.blog.free.fr/public/perso5/olssc31.gif et celle sur http://freakonometrics.blog.free.fr/public/perso5/olssc32.gif. On notera http://freakonometrics.blog.free.fr/public/perso5/olssc.gif et http://freakonometrics.blog.free.fr/public/perso5/olssc33.gif les http://freakonometrics.blog.free.fr/public/perso5/r2ajuste.gif associés. Pesaran (1974) a noté qu'il existait une relation simple entre  http://freakonometrics.blog.free.fr/public/perso5/olssc.gif et http://freakonometrics.blog.free.fr/public/perso5/olssc33.gif

http://freakonometrics.blog.free.fr/public/perso5/olssc36.gif

  oùhttp://freakonometrics.blog.free.fr/public/perso5/ols37.gif est l'estimateur du coefficient ajouté dans la régression la plus complète, et http://freakonometrics.blog.free.fr/public/perso5/olssc51.gif la statistique de Student du test de significativité. Autrement dit, http://freakonometrics.blog.free.fr/public/perso5/olssc38.gif si et seulement si la nouvelle variable n'est pas suffisamment significative, en l’occurrence http://freakonometrics.blog.free.fr/public/perso5/olssc39.gif.
Mais ce n'est pas le coefficient http://freakonometrics.blog.free.fr/public/perso5/r2ajuste.gif (même avec l'ajustement) qui est le plus utilisé pour faire de la sélection de modèle, comme nous le verrons en cours la semaine prochaine (mais une pénalisation de la vraisemblance). A suivre donc...

On linear models with no constant and R2

In econometrics course we always say to our students that "if you fit a linear model with no constant, then you might have trouble. For instance, you might have a negative R-squared". So I tried to find databases on the internet such that, when we compute a linear regression, we actually obtain a negative R squared. I have generated hundreds to random databases that should exhibit such a property, in R. With no success. Perhaps to be more specific, I should explain what might happen if we do not include a constant in a linear model. Consider the following dataset, where points are on a straight line, with a negative slope, far from the origin, symmetric with respect to the first diagonal.

> x=1:3
> y=3:1
> plot(x,y)

Points are on a straight line, so it is actually possible to get a perfect linear model. But only if we integrate a constant in our model. This is related to the fact that the correlation between our two variates is -1,

> cor(x,y)
[1] -1

The least-square program is here

http://freakonometrics.blog.free.fr/public/perso5/olssc01b.gif

i.e. the estimate of the slope is

http://freakonometrics.blog.free.fr/public/perso5/olcsc02.gif

Numerically, we obtain

> sum(x*y)/sum(x^2)
[1] 0.7142857

which is the actual slope on the illustration above. If we compute the sum of squares of errors (as a function of the slope), we have here

> ssr=function(b){sum((y-b*x)^2)}
> SSR=Vectorize(ssr)
> B=seq(-1,3,by=.1)
> plot(B,SSR(B),ylim=c(0,ssr(3)),cex=.6,type="b")

so the value we have computed is actually the minimum of the sum of squares of errors. But note that the sum of squares always exceeds the total sum of squares in red on the graph above

> ssr(b)
[1] 6.857143
> sum((y-mean(y))^2)
[1] 2

Recall that the total "coefficient of variation", denoted http://freakonometrics.blog.free.fr/public/perso5/R2.gif, is defined as

http://freakonometrics.blog.free.fr/public/perso5/olsnc04.gif

i.e.

> 1-ssr(b)/sum((y-mean(y))^2)
[1] -2.428571

which is negative. It is also sometimes defined as "the square of the sample correlation coefficient between the outcomes and their predicted values". Here it would be related to 

> cor(b*x,y)
[1] -1

so we would have a unit http://freakonometrics.blog.free.fr/public/perso5/R2.gif . So obviously, using the http://freakonometrics.blog.free.fr/public/perso5/R2.gif in a model without a constant would give odd results. But the weird part is that if we run that regression with R, we get

> summary(lm(y~0+x))
 
Call:
lm(formula = y ~ 0 + x)
 
Residuals:
1       2       3
2.2857  0.5714 -1.1429
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x   0.7143     0.4949   1.443    0.286
 
Residual standard error: 1.852 on 2 degrees of freedom
Multiple R-squared: 0.5102,	Adjusted R-squared: 0.2653
F-statistic: 2.083 on 1 and 2 DF,  p-value: 0.2857 

Here, the estimation is correct. But the http://freakonometrics.blog.free.fr/public/perso5/R2.gif we obtain tells us that the model is not that bad... So if anyone knows what R computes, I'd be glad to know. The value given by R (thanks Vincent for asking me to look carefully at the R source code) is obtained using Pythagoras's theorem to compute the total sum of square,

> sum((b*x)^2)/(sum((b*x)^2)+sum((y-b*x)^2))
[1] 0.5102041

So be careful, the http://freakonometrics.blog.free.fr/public/perso5/R2.gif might look good, but meaningless !

Monday, January 30 2012

calculs des sorties de régression

Avec un peu de calcul matriciel (et en utilisant les fonctions de répartition des lois usuelles) on peut calculer à la main toute les grandeurs obtenues lorsque l'on fait une régression linéaire. Utilisons la base de donnée avec le poids d'une personne (qui sera notre variable d'intérêt, et que l'on essayera d'inférer), la taille et le sexe de la personne,

> TP = read.table(
"http://freakonometrics.free.fr/taillepoids.txt",
header=TRUE,sep=";") > tail(TP) indice sexe poids1 taille1 poids2 taille2 195 195 F 62 164 61 161 196 196 M 74 175 71 175 197 197 M 83 180 80 180 198 198 M 81 175 NA NA 199 199 M 90 181 91 178 200 200 M 79 177 81 178 > Y=TP$poids1 > X1=TP$sexe > X2=TP$taille1 > base=data.frame(Y,X1,X2) > tail(base) Y X1 X2 195 62 F 164 196 74 M 175 197 83 M 180 198 81 M 175 199 90 M 181 200 79 M 177

La régression linéaire donne ici

> regression=lm(Y~X1+X2)
> summary(regression)
 
Call:
lm(formula = Y ~ X1 + X2)
 
Residuals:
Min      1Q  Median      3Q     Max
-24.718  -7.159  -1.472   5.487  74.726
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 109.11420   14.20513   7.681 7.22e-13 ***
X1M          22.49801    2.08690  10.781  < 2e-16 ***
X2           -0.31298    0.08649  -3.619 0.000376 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 11.81 on 197 degrees of freedom
Multiple R-squared: 0.3937,	Adjusted R-squared: 0.3875
F-statistic: 63.95 on 2 and 197 DF,  p-value: < 2.2e-16

On peut commencer par (re)calculer les estimateurs par moindres carrés ordinaires, en notant que de la variable associée au sexe, seul le sexe masculin a été retenu ici.

> n=length(Y)
> X=matrix(c(rep(1,n),(X1=="M")*1,X2),n,3)
> X[1:5,]
[,1] [,2] [,3]
[1,]    1    1  182
[2,]    1    0  161
[3,]    1    0  161
[4,]    1    1  177
[5,]    1    0  157
> (beta=solve(t(X)%*%X)%*%t(X)%*%Y)
[,1]
[1,] 109.1141967
[2,]  22.4980107
[3,]  -0.3129827
> summary(regression)$coefficients
Estimate  Std. Error   t value     Pr(>|t|)
(Intercept) 109.1141967 14.20512536  7.681326 7.222037e-13
X1M          22.4980107  2.08690147 10.780581 1.347214e-21
X2           -0.3129827  0.08648531 -3.618912 3.760466e-04

On retrouve ici est estimation indiquée. Notons que si la modalité de référence était le sexe féminin, nous aurions,

> Xbis=matrix(c(rep(1,n),(X1=="F")*1,X2),n,3)
> Xbis[1:5,]
[,1] [,2] [,3]
[1,]    1    0  182
[2,]    1    1  161
[3,]    1    1  161
[4,]    1    0  177
[5,]    1    1  157
> solve(t(Xbis)%*% Xbis)%*%t(Xbis)%*%Y
[,1]
[1,] 131.6122073
[2,] -22.4980107
[3,]  -0.3129827

Si on compare les deux constantes, entre ces deux modèles, on retrouve la valeur du coefficient associé au sexe de l'individu,

> (solve(t(X)%*%X)%*%t(X)%*%Y)[1]-
+ (solve(t(Xbis)%*% Xbis)%*%t(Xbis)%*%Y)[1]
[1] -22.49801
> (solve(t(X)%*%X)%*%t(X)%*%Y)[2]
[1] 22.49801

Continuons un peu... L'estimateur de la volatilité des résidus est ici

> residus=Y-X%*%beta
> t(residus)%*%residus
[,1]
[1,] 27493.32
> sum(residus^2)
[1] 27493.32
>
> (sigma=sqrt(sum(residus^2)/(n-ncol(X))))
[1] 11.81355
> summary(regression)$sigma
[1] 11.81355

On peut utiliser cette grandeur pour estimer la matrice de variance-covariance de notre estimateur,

> (varbeta=sigma^2*solve(t(X)%*%X))
[,1]       [,2]         [,3]
[1,] 201.785586 16.2312631 -1.224735573
[2,]  16.231263  4.3551577 -0.106737634
[3,]  -1.224736 -0.1067376  0.007479709
> vcov(regression)
(Intercept)        X1M           X2
(Intercept)  201.785586 16.2312631 -1.224735573
X1M           16.231263  4.3551577 -0.106737634
X2            -1.224736 -0.1067376  0.007479709

Le terme diagonal est celui qui est utilisé pour calculer l'écart-type des estimateurs de chacun des coefficients.

> (sigbeta=sqrt(diag(varbeta)))
[1] 14.20512536  2.08690147  0.08648531
> summary(regression)$coefficients
Estimate  Std. Error   t value     Pr(>|t|)
(Intercept) 109.1141967 14.20512536  7.681326 7.222037e-13
X1M          22.4980107  2.08690147 10.780581 1.347214e-21
X2           -0.3129827  0.08648531 -3.618912 3.760466e-04

En faisant une hypothèse de normalité des résidus (ou en supposant le nombre d'observations grand, i.e. suffisamment grand pour invoquer des théorèmes asymptotiques), on peut faire un test de significativité des variables explicatives,

> (beta/sigbeta)
[,1]
[1,]  7.681326
[2,] 10.780581
[3,] -3.618912
> 2*(1-pt(abs(beta/sigbeta),df=n-ncol(X)))
[,1]
[1,] 7.220891e-13
[2,] 0.000000e+00
[3,] 3.760466e-04

avec respectivement la statistique de Student, et la p-value associé au test de nullité. Pour juger de la qualité de l'ajustement, on peut calculer

> (R2=cor(Y,X%*%beta)^2)
[,1]
[1,] 0.3936724
> summary(regression)$r.squared
[1] 0.3936724
>
> (R2a=1-(n-1)/(n-(ncol(X)))*(1-R2))
[,1]
[1,] 0.3875168
> summary(regression)$adj.r.squared
[1] 0.3875168

et enfin, le test de significativité global (de Fisher) donne ici

> (F=R2/(ncol(X)-1)/((1-R2)/(n-ncol(X))))
[,1]
[1,] 63.95343
> summary(regression)$fstatistic
value     numdf     dendf
63.95343   2.00000 197.00000
> 1-pf(F,(ncol(X)-1),(n-ncol(X)))
[,1]
[1,]    0

Bref, avec un peu de calcul, on peut retrouver toutes les valeurs données sous R. Une fois les calculs compris, on pourra passer à la géométrie, et à l'interprétation...

Thursday, May 21 2009

Le prochain qui me parle du R2...

Depuis que j'ai repris la recherche, j'ai pris l'habitude de débrancher mon téléphone. Mais lundi matin, comme j'avais un coup de fil à passer, j'ai du le rebrancher. J'ai alors reçu un appel d'un élève qui m'a agacé (et auquel j'ai eu la lâcheté de raccrocher en expliquant que je faisais autre chose - ce qui était néanmoins partiellement vrai) mais qui est assez symptomatique, "bonjour, j'ai fait une régression et j'ai un R2 de 0,312 et je voulais savoir si c'était bien ?" (ou "si c'était beaucoup ?", ou "si c'était assez ?", ou un truc du genre car j'avoue que j'avais effectivement la tête ailleurs).
Pourtant je croyais que mon cours avait précisément réussi à éviter de parler du R2. Je vais donc reprendre ici quelques points. J'espérais trouver un site internet de l'association des statisticiens qui en ont marre du R2 (ou d'opposants à la sectes des adorateurs du R2, car si j'ignore qui en est le gourou, je trouve que cette secte regroupe beaucoup d'adorateurs) malheureusement je n'ai rien trouvé de concluant sur internet. Donc je vais m'y mettre... même je risque de faire un très long billet,
  • définition et interprétation du R2

Un "bon" modèle permettra d'obtenir des estimations "proches" des valeurs observées. Sur la représentation dans l’espace des variables (comme le montre la figure ci-dessus) la qualité peut être évaluée par l’angle θ. Cet angle est compris entre -90˚et 90˚. Un angle proche de -90˚ou de 90˚indique un modèle de mauvaise qualité. Alors qu'un angle proche de 0˚ correspond à un bon modèle. Le cosinus carré de θ est donc une mesure possible de la qualité du modèle et cette mesure varie entre 0 et 1.
En utilisant le théorème de Pythagore nous permet d'écrireou encoreoù SCT (respectivement SCE et SCR) désigne la somme des carrés totale (respectivement expliquée par le modèle et résiduelle). Le coefficient de détermination R2 est défini parc’est-à-dire que le R2 est souvent interprété comme le "pourcentage de la variabilité de la variable endogène expliquée par le modèle".
Si on reprend une page du livre de Pierre-André Cornillon et Éric Matzner-Løber, on a les définitions suivantes
  • que mesure le R2 ? quelques exemples...
Comme je le dis auparavant, le R2 nous dit si la courbe de régression (pour prendre un modèle plus général que le cas linéaire) est proche des observations, mais ne nous dit nullement si le modèle est "bon".
Pour faire simple, avec des observations individuelles, le R2 sera souvent assez "faible", alors qu'avec des séries temporelles (en particulier des séries intégrées), le R2 sera souvent "plus élevé". Dans ce dernier cas, l'explication heuristique est qu'expliquer une série croissante par une autre série croissante marche souvent "bien".Sur l'exemple ci-dessus, on obtient - en ajustant un modèle linéaire - respectivement
Multiple R-squared: 0.3372,
Multiple R-squared: 0.6414,
Sur l'exemple ci-dessus, on obtient - en ajustant un modèle polynôme de degré 4 - respectivement
Multiple R-squared: 0.4968,
Multiple R-squared: 0.987,

Sur l'exemple 3, l'ajustement est "bon" au sens où il serait difficile de faire mieux, c'est simplement qu'il est très bruité. Avec un R2 de 0,5, l'exemple 3 me paraît proposer un "meilleur" ajustement que l'exemple 2, alors que le R2 est de 0,65. Dans le premier cas, l'ajustement semble bon, mais très bruité, alors que pour le second, le modèle est plutôt mauvais, mais le bruit est relativement faible.
  • le R2 n'est qu'un coefficient de corrélation, ou presque (et donc possède tous les défauts d'un coefficient de corrélation)
Si on reprend le premier point, il est possible d'écrire le R2 (dans un modèle avec une unique variable explicative) comme le carré du coefficient de corrélation (au sens de Pearson) entre la variable endogène et la variable explicative (plus généralement, ça serait entre les observations et la prédiction). Regardons par exemple sur la base de données cars,
> data(cars)
> cor(cars)
speed dist
speed 1.0000000 0.8068949
dist 0.8068949 1.0000000

Si on regarde la valeur du R2 on peut écrire
> summary(lm(dist~speed,data=cars))$r.squared
[1] 0.6510794
> sqrt(summary(lm(dist~speed,data=cars))$r.squared)
[1] 0.8068949

J'avais dit dans un précédant billet (ici) tout le mal que je pense du coefficient de corrélation, qui n'est pas une mesure de dépendance (ou une mesure de concordance au sens où l'avait défini Marco Scarsini, ici). Pour cette raison, je trouve qu'un R2 qui vaut 0,312 nous apprend aussi peu de chose qu'un coefficient de corrélation de 0,217 ou 0,013.
  • mais le R2 n'est pas exactement le carré d'une quantité (et donc peut être négatif, par exemple)
Formellement, le R2 est le carré du cosinus d'un angle, si on reprend l'interprétation géométrique de la régression linéaire. Mais si on regarde la version empirique, on peut être (désagréablement) surpris... En particulier avec Excel. Ce point est discuté sur de nombreux forum (ici ou ), je ne vais donc pas en rajouter.
  • le R2 augmente en rajoutant des variables explicatives (ce qui le rend difficilement utilisable pour comparer des modèles)
Ce point est souvent connu par tout le monde, c'est pour cela qu'on présente toujours le R2 ajusté après avoir présenté le R2.
  • le R2 ne peut pas être utiliser pour comparer des modèles pour lesquels la variable endogène change
Mais ceci n'est pas qu'un problème avec le R2, car il est très difficile de trouver un critère de choix suffisamment robuste. J'avais déjà abordé ce point lors de mon cours d'assurance non-vie à l'ensae il y a quelques années. J'avais vu des élèves essayer de comparer deux modèles (pour tarifer), à savoir
  • une régression gamma, i.e.
  • une "régression lognormale", c'est à dire que
le soucis est qu'il est difficile de comparer des deux modèles. En particulier, tous les élèves qui se posaient cette question à l'aide d'un critère de type AIC concluaient toujours à la supériorité de la loi lognormale, ou pour être un peu pointilleux, sur un modèle gaussien après une transformation logarithmique.
Dans le cas des modèles Gaussiens, c'est aussi la conclusion à laquelle arrive Valérie Mignon dans son livre d'économétrie, à savoir que si l'on cherche à comparer
  • un modèle additif, i.e.
  • un modèle multiplicatif, i.e.
on retient "toujours" le modèle multiplicatif si on retient le modèle qui a le meilleur R2. Je n'ai pas de raison a priori d'aller dans ce sens, mais je retiens que juger de la supériorité d'un modèle à l'aide uniquement de ce critère me paraît vraiment douteux.
L'exemple ci-dessous représente les deux régressions, avec un modèle linéaire versus un modèle exponentiel (i.e. avec comme variable endogène le logarithme).
Dans ce cas précis,
> REG.LINEAIRE=lm(Y~X)
> REG.MULTIPLICATIF=lm(log(Y)~X)
> summary(REG.LINEAIRE)$r.squared
[1] 0.8688
> summary(REG.MULTIPLICATIF)$r.squared
[1] 0.8654
Afin de répondre à cet exemple précis, je renvoie au test de Box-Cox (que j'avais abordé ici) qui vise précisément à comparer ces deux modèles. Ici, on conclue à la supériorité d'un modèle multiplicatif, mais ça n'a rien à voir avec la valeur du R2.

  • le R2 peut être changé artificellement (en réécrivant le modèle)
Je reprendrais un exemple tiré d'une épreuve du concours interne d'administrateur de l'insee (qui doit avoir une vingtaine d'années). Considérons un modèle assez classique en économie, où on essaye de lier taux d'intérêt et taux d'inflation. Je renvois à mes slides du cours d'économétrie où l'exemple est présenté, et mis en œuvre, ici (slides 39 à 42).
  • le R2 peut être changé artificellement (en supprimant les points aberrants)

Ce point est présenté dans certains livres, où il est expliqué qu'en enlevant des points (souvent les premières années quand on travaille sur des séries macroéconomiques) on peut améliorer le R2. Effectivement, si on enlève les quatre points les plus éloignés de la prédiction, on change le modèle, et le R2 est alors amélioré,
> summary(REG.TOTAL)$r.squared
[1] 0.6510794
> summary(REG.SUBSET)$r.squared
[1] 0.7539631

Je croyais là aussi que mon cours essayait d'expliquer l'importance (et le traitement) de ces points aberrants.
Pour conclure, je reprendrais ce que j'essaye d'expliquer dans mon cours d'économétrie (ici ou là), à savoir qu'au lieu de se contenter de regarder un R2, il vaut mieux passer un peu de temps à regarder quelques graphiques de diagnostique. Et pour comparer des modèles, je préfère toujours les critères AIC ou BIC au R2 (ou même au R2 ajusté). Quant à la normalité des résidus... je pense que je reviendrais là dessus aussi dans un billet les jours à venir....

Wednesday, March 18 2009

Constante dans une régression

Je reviens ici sur une question qui m'avait été posée il y a quelques mois en TD sur "supprimer la constante dans une régression". J'avais expliqué qu'il fallait faire très attention sur l'interprétation d'une régression sans constante.
J'en profite d'ailleurs pour noter que ce point était évoqué par Emmeline sur le blog http://www.mafeco.fr/ (qui pointait d'ailleurs ici, tout en mentionnant que mon blog n'était "pas vraiment grand public"... mais j'assume).
Bref, pour expliquer un peu cette histoire de constante, je ferais comme Boris Vian qui prétendait que "tout a déjà été dit cent fois, et beaucoup mieux  que par moi", donc je me contenterais de reprendre un passage du Davidson et MacKinnon (qu'on peut retrouver ici).
Bref, la crainte principale est que les résidus ne soient pas centrés, mais ici c'est surtout des discussions sur l'utilisation (et donc l'interprétation) du R2 qui sont évoquées. Je l'ai suffisamment dit dans mon cours, je ne suis pas un grand fan de cette mesure. Tout d'abord car il s'agit simplement d'une mesure de corrélation et que donc elle est très dépendante des lois sous-jacente. Un R2 de 0,4 (comme une corrélation de 0,4 d'ailleurs) ne veut rien dire en soit.