Tag Archives: R2

Random thoughts on econometric models with (pure) random features

For my lectures on applied linear models, I wanted to illustrate the fact that the R^2 is never a good measure of the goodness of the model, since it’s quite easy to improve it. Consider the following dataset


with one variable of interest y, and 99 features x_j. All of them being (by construction) independent. And we have 100 observations… Consider here the regression on the first k features, and compute R_k^2 of that regression


Let us see what’s going on…


(actually, it’s not exactly what we have on the graph…. we have the average obtained over 1,000 samples randomly generated, with 90% confidence bands). Oberve that \mathbb{E}[R^2_k]=k/n, i.e. if we add some pure random noise, we keep increasing the R^2 (up to 1, actually).

Good news, as we’ve seen in the course, the adjusted R^2 – denoted \bar R^2-might help. Observe that \mathbb{E}[\barR^2_k]=0, so, in some sense, adding features does not help here…


We can actually do the same with Akaike criteria AIC_k and Schwarz (bayesian) criteria BIC_k.


For the AIC, the intitial increase makes sense : we should not prefer the model with 10 covariates, compared with nothing. The strange thing is the far right behavior : we prefer here 80 random noise features to none ! Which I find hard to interprete… For the BIC the code is simply


and here also, we have the same pattern, where we prefer a big model with juste pure noise to nothing…

A last one to conclude (or not) : what about the leave-one-out cross validation mean squared error ? More precisely, CV=\frac{1}{n}\sum_{i=1}\widehat{\varepsilon}^2_{-i}where \widehat{\varepsilon}^2_{-i}=y_i-\widehat{y}_{-i} where \widehat{y}_{-i} is the predicted value obtained with the model is estimated when the ith observation is deleted. One can prove that \widehat{\beta}_{-i}=\widehat{\beta}-(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{x}_i\hat\varepsilon_i(1-H_{i,i})^{-1}where H is the classical hat matrix, thus\widehat{\varepsilon}_{-i}=(1-H_{i,i})^{-1}\hat\varepsilon_ii.e. we do note have to estimate (at each round) n models

  mean( (residuals(model)/1-h)^2 ))}

Here, it make sense : adding noisy features yields overfit ! So the mean squared error is decreasing !

That’s all nice, but it might not be very realistic… Here, for my model with only one variable, I just pick one, at random…. In practice, we try to get the “best one”… So a more natural idea would be to order the variables according to their correlations with y,


and as before, we can plot the evolution of R^2_k as a function of k the number of features considered,

which is increasing, with a higher slope at the beginning… For the \bar R^2_k we might actually prefer a correlated noise to nothing (which makes sense actually). So here since we somehow chose our variables, \bar R^2_k seems to be always positive…

For the AIC_k here also, there is an improvement. Before coming back to the original situation (with about 80 features) and here also, we observe the drop on the far right part of the graph

The BIC_k might like the top three features, but soon, we have a deterioration…. even if here also, we have the drop at the far right (with more than 95 features… for 100 observations).

Finally, observe that here again, our (leave-one-out) cross-validation has not been mesled by our noisy variables : it is always decreasing !

So it seems that cross-validation techniques are more robust than the AIC and BIC (even if we mentioned in a previous post connexions between all those concepts) when we have a lot a noisy (non-relevent) features.

Probabilistic Foundations of Econometrics, part 4

This post is the fourth one of our series on the history and foundations of econometric and machine learning models. Part 3 is online here.

Goodness of Fit, and Model

In the Gaussian linear model, the determination coefficient – noted R^2 – is often used as a measure of fit quality. It is based on the variance decomposition formula \underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\bar{y})^2}_{\text{total variance}}=\underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{\text{residual variance}}+\underbrace{\frac{1}{n}\sum_{i=1}^n (\widehat{y}_i-\bar{y})^2}_{\text{explained variance}} The R^2 is defined as the ratio of explained variance and total variance, another interpretation of the coefficient that we had introduced from the geometry of the least squares R^2= \frac{\sum_{i=1}^n (y_i-\bar{y})^2-\sum_{i=1}^n (y_i-\widehat{y}_i)^2}{\sum_{i=1}^n (y_i-\bar{y})^2}The sums of the error squares in this writing can be rewritten as a log-likelihood. However, it should be remembered that, up to one additive constant (obtained with a saturated model) in generalized linear models, deviance is defined by {Deviance}(\widehat{\beta}) = -2\log[\mathcal{L}] which can also be noted Deviance(\widehat{\mathbf{y}}). A null deviance can be defined as the one obtained without using the explanatory variables \mathbf{x}, so that \widehat{y}_i=\overline{y}. It is then possible to define, in a more general context (with a non-Gaussian distribution for y)R^2=\frac{{Deviance}(\overline{y})-{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}=1-\frac{{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}However, this measure cannot be used to choose a model, if one wishes to have a relatively simple model in the end, because it increases artificially with the addition of explanatory variables without significant effect. We will then tend to prefer the adjusted R^2,\bar R^2 = {1-(1-R^{2})\cdot{n-1 \over n-p}} = R^{2}-\underbrace{(1-R^{2})\cdot{p-1 \over n-p}}_{\text{penalty}}where p is the number of parameters of the model. Measuring the quality of fit will penalize overly complex models.

This idea will be found in the Akaike criterion, where AIC=Deviance+2\cdot p or in the Schwarz criterion, BIC=Deviance+log(n)\cdot p. In large dimensions (typically p>\sqrt{n}), we will tend to use a corrected AIC, defined by AIC_c=Deviance+2⋅p⋅n/(n-p-1) .

These criterias are used in so-called “stepwise” methods, introducing the set methods. In the “forward” method, we start by regressing to the constant, then we add one variable at a time, retaining the one that lowers the AIC criterion the most, until adding a variable increases the AIC criterion of the model. In the “backward” method, we start by regressing on all variables, then we remove one variable at a time, removing the one that lowers the AIC criterion the most, until removing a variable increases the AIC criterion from the model.

Another justification for this notion of penalty (we will come back to this idea in machine learning) can be the following. Let us consider an estimator in the class of linear predictors, \mathcal{M}=\big\lbrace m:~m(\mathbf{x})=s_h(\mathbf{x})^T\mathbf{y} \text{ where }S=(s(\mathbf{x}_1),\cdots,s(\mathbf{x}_n))^T\text{ is some smoothing matrix}\big\rbrace and assume that y=m_0 (x)+\varepsilon, with \mathbb{E}[\varepsilon]=0 and Var[\varepsilon]=\sigma^2\mathbb{I}, so that m_0 (x)=\mathbb{E}[Y|X=x] . From a theoretical point of view, the quadratic risk, associated with an estimated model \widehat{m}, \mathbb{E}\big[(Y-\widehat{m}(\mathbf{X}))^2\big], is written\mathcal{R}(\widehat{m})=\underbrace{\mathbb{E}\big[(Y-m_0(\mathbf{X}))^2\big]}_{\text{error}}+\underbrace{\mathbb{E}\big[(m_0(\mathbf {X})-\mathbb{E}[\widehat{m}(\mathbf{X})])^2\big]}_{\text{bias}^2}+\underbrace{\mathbb{E}\big[(\mathbb{E}[\widehat{m}(\mathbf{X})]-\widehat{m}(\mathbf{X}))^2\big]}_{\text{variance}} if m_0 is the true model. The first term is sometimes called “Bayes error”, and does not depend on the estimator selected, \widehat{m}.

The empirical quadratic risk, associated with a model m, is here: \widehat{\mathcal{R}}_n(m)=\frac{1}{n}\sum_{i=1}^n (y_i-m(\mathbf{x}_i))^2 (by convention). We recognize here the mean square error, “mse”, which will more generally give the “risk” of the model m when using another loss function (as we will discuss later on). It should be noted that:\displaystyle{\mathbb{E}[\widehat{\mathcal{R}}_n(m)]=\frac{1}{n}\|m_0(\mathbf{x})-m(\mathbf{x})\|^2+\frac{1}{n}\mathbb{E}\big(\|{Y}-m_0(\mathbf{X})\|^2\big)} We can show that:n\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]=\mathbb{E}\big(\|Y-\widehat{m}(\mathbf{x})\|^2\big)=\|(\mathbb{I}-\mathbf{S})m_0\|^2+\sigma^2\|\mathbb{I}-\mathbf{S}\|^2so that the (real) risk of \widehat{m} is: {\mathcal{R}}_n(\widehat{m})=\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]+2\frac{\sigma^2}{n}\text{trace}(\boldsymbol{S})So, if \text{trace}(\boldsymbol{S})\geq0 (which is not a too strong assumption), the empirical risk underestimates the true risk of the estimator. Actually, we recognize here the number of degrees of freedom of the model, the right-hand term corresponding to Mallow’s C_p, introduced in Mallows (1973) using not deviance but R^2.

Statistical Tests

The most traditional test in econometrics is probably the significance test, corresponding to the nullity of a coefficient in a linear regression model. Formally, it is the test of H_0:\beta_k=0 against H_1:\beta_k\neq 0. The so-called Student test, based on the statistics t_k=\widehat{\beta}_k/se_{\widehat{β}_k}, allows to decide between the two alternatives, using the test p-value, defined by \mathbb{P}[|T|>|t_k|] avec T\overset{\mathcal{L}}{\sim} Std_\nu, where \nu is the number of degrees of freedom of the model (\nu=p+1 for the standard linear model). In large dimension, however, this statistic is of very limited interest, given a significant FDR (“False Discovery Ratio”). Classically, with a level of significance \alpha=0.05, 5% of the variables are falsely significant. Suppose that we have p=100 explanatory variables, but that 5 (only) are really significant. We can hope that these 5 variables will pass the Student test, but we can also expect that 5 additional variables (false positive test) will emerge. We will then have 10 variables perceived as significant, while only half are significant, i.e. an FDR ratio of 50%. In order to avoid this recurrent pitfall in multiple tests, it is natural to use the procedure of Benjamini & Hochberg (1995).

From a correlation to some causal effect

Econometric models are used to implement public policy evaluations. It is therefore essential to fully understand the underlying mechanisms in order to know which variables actually make it possible to act on a variable of interest. But then we move on to another important dimension of econometrics. Jerry Neyman was responsible for the first work on the identification of causal mechanisms, and then Rubin (1974) formalized the test, called the “Rubin causal model” in Holland (1986). The first approaches to the notion of causality in econometrics were based on the use of instrumental variables, models with discontinuity of regression, analysis of differences in differences, and natural or unnatural experiments. Causality is usually inferred by comparing the effect of a policy – or more generally of a treatment – with its counterfactual, ideally given by a random control group. The causal effect of the treatment is then defined as \Delta=y_1-y_0, i.e. the difference between what the situation would be with treatment (noted t=1) and without treatment (noted t=0). The concern is that only y=t\cdot y_1+(1-t)\cdot y_0 and t are observed. In other words, the causal effect of variable t  on t  is not observed (since only one of the two potential variables – y_0 or y_1  is observed for each individual), but it is also individual, and therefore a function of x-covariates. Generally, by making assumptions about the distribution of the triplet (Y_0,Y_1,T) , some parameters of the causal effect distribution become identifiable, based on the density of the observable variables (Y,T) . Classically, we will be interested in the moments of this distribution, in particular the average effect of treatment in the population, \mathbb{E}[\Delta] , or even just the average effect of treatment in the case of treatment \mathbb{E}[\Delta|T=1] . If the result (Y_0,Y_1) is independent of the processing access variable T, it can be shown that \mathbb{E}[\Delta]=\mathbb{E}[Y|T=1]- \mathbb{E} [Y|T=0]. But if this independence hypothesis is not verified, there is a selection bias, often associated with \mathbb{E}[Y_0|T=1]- \mathbb{E} [Y_0|T=0]. Rosenbaum & Rubin (1983) propose to use a propensity to be treated score, p(x)=\mathbb{P}[T=1|X=x] , noting that if variable Y_0\ is independent of access to treatment T conditionally to the explanatory variables X, then it is independent of T  conditionally to the score p(X) : it is sufficient to match them using their propensity score. Heckman et al (2003) thus proposes a kernel estimator on the propensity score, which simply provides an estimator of the effect of the treatment, provided that it is treated.

To be continued next time, we’ll introduce “machine learning techniques” (references mentioned above are online here)

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


i.e. the estimate of the slope is


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



> 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))

lm(formula = y ~ 0 + x)

1       2       3
2.2857  0.5714 -1.1429

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 !

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,
> 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….