Tag Archives: regression

sur les transformations dans un modèle linéaire

Je voulais prendre 5 minutes pour reprendre une question posée par courriel, qui me permettra de poursuivre sur des choses évoquées mercredi dernier en cours. Je vais reformuler la question, mais en gros, cela disait: la méthode de Box-Cox (évoquée en début de semaine, ici) avait pour objet de choisir entre un modèle linéaire (que l’on étudié dans tous les sens depuis le début du cours) et un modèle log-linéaire (évoqué ici, par exemple). L’idée était d’associer le cas linéaire à  (dans la transformée de Box-Cox) et à  pour le cas multiplicatif (ou log-linéaire). Mais que se passe-t-il si la valeur optimale rejette ces deux cas, et est proche (disons) de  ?
Considérons le cas suivant (exemple que je ne me lasse d’utiliser)

> reglm=lm(dist~speed,data=cars)
> library(MASS)
> boxcox(reglm)

La valeur optimale est effectivement proche de . En posant , on aurait envie de considérer un modèle de la forme

(comme le suggère la transformation de Box-Cox). On peut faire la régression,

> regsqrt=lm(sqrt(dist)~speed,data=cars)
+ summary(regsqrt)

Call:
lm(formula = sqrt(dist) ~ speed, data = cars)

Residuals:
Min      1Q  Median      3Q     Max
-2.0684 -0.6983 -0.1799  0.5909  3.1534

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.27705    0.48444   2.636   0.0113 *
speed        0.32241    0.02978  10.825 1.77e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.102 on 48 degrees of freedom
Multiple R-squared: 0.7094,	Adjusted R-squared: 0.7034
F-statistic: 117.2 on 1 and 48 DF,  p-value: 1.773e-14

et effectivement, le résultat est concluant. Mais on ne peut pas en rester là, comme pour pour le passage au logarithme, notre variable d’intéret reste .
La prédiction sur le modèle transformé était

et

Comme avec le logarithme, on pourrait prendre comme prédiction pour  mais

Bref, on a un estimateur biaisé. Et là encore, comme on est sur des variables positives, l’inégalité de Jensen doit meme nous garantir que l’une des quantités domine toujours l’autre (je laisse les amateurs de convexité l’écrire dans le bon sens). Comment faire ?
La première solution est de noter que

Donc de manière très simple, on a notre prédiction (en prenant le carré, comme intuité, mais en rajoutant une terme – positif – lié à la variance).
Une seconde consiste à noter que Z (conditionnellement à la variable explicative) était gaussienne. Donc en prenant le carré (moyennant quelques changement d’échelle), on devrait tomber sur une loi du chi-deux qui est une loi que l’on connait bien (sinon je peux renvoyer ici).
Essayons de creuser un peu ces idées (en particulier, pour dériver des intervalles de confiance). Si on visualise la prédiction de ce modèle sur la racine carrée de la distance de freinage, on obtient

> plot(speed,sqrt(dist))
> x=seq(0,30,by=.2)
> distsqrtp=predict(regsqrt,newdata=
+ data.frame(speed=x),interval="prediction")
> polygon(c(x,rev(x)),c(distsqrtp[,2],
+ rev(distsqrtp[,3])),col="yellow",border=NA)
> lines(x,distsqrtp[,1],lwd=2,col="red")

Pour passer d’une loi normale à une loi du chi-deux, il faut retrancher la moyenne (pour centrer la variable) et diviser par l’écart-type (pour avoir une variance unitaire),

> s=summary(regsqrt)$sigma
> mu=predict(regsqrt)
> distsqrtp01=(sqrt(dist)-mu)/s

On a ainsi une variable qui, conditionnellement à la variable explicative est supposé suivre une loi ,

> plot(speed,distsqrtp01,ylim=c(-3,3))
> abline(h=qnorm(c(.025,.975),0,1),lty=2,col="red")

On prend ensuite le carré pour obtenir notre loi , on on trace la bande de confiance (il faut que la valeur à la deux soit faible),

> plot(speed,distsqrtp01^2)
> abline(h=qchisq(.95,df=1),lty=2,col="red")

On a maitenant nos intervalles de confiance, en utilisant cette loi du chi-deux, en inversant notre transformation,

> mu=predict(regsqrt,newdata=data.frame(speed=x))
> distsup=(mu+s*sqrt(qchisq(.95,df=1)))^2
> distinf=(mu-s*sqrt(qchisq(.95,df=1)))^2

que l’on peut visualiser simplement (sur les données de base, avec en ordonnée la distance)

> plot(cars)
> lines(x,distsup,lty=2,col="red")
> lines(x,distinf,lty=2,col="red")

Et la valeur prédite ? On va utiliser notre relation liant la variance, l’espérance du carré, et le carré de l’espérance,

> distesp=mu^2+s^2
> lines(x,distesp,lwd=2,col="red")

Nice, n’est-ce pas ?
Sinon, comme je le disais en cours, si la transformation optimale sur  est de prendre , peut-etre pourrait-on envisager de prendre  (de manière un peu duale) comme variable explicative. Encore une fois, c’est une idée, car rien ne le garantit. En effet

En particulier, le second modèle est un modèle linéaire classique, avec de l’homoscédasticité: on aura une dispersion uniforme autour de notre parabole. En revanche, avec le premier modèle, on a une espèce d’hétéroscédasticité (comme tenu du double produit), avec une variance du terme d’erreur qui croit avec  (car les coefficients sont positifs). Bref, en terme d’intervalles de confiance, on devrait avoir des choses assez différentes.
Regardons la régression de la distance (cette fois) sur la vitesse, et le carré de la vitesse,

> reglm=lm(dist~speed+I(speed^2),data=cars)
> distp=predict(reglm,newdata=
+ data.frame(speed=x),interval="prediction")

Si on visualise la prédiction, avec un intervalle de confiance, on obtient

> plot(cars)
> polygon(c(x,rev(x)),c(distp[,2],
+ rev(distp[,3])),col="yellow",border=NA)
> lines(x,distp[,1],lwd=2,col="red")

Autrement dit, on suppose que notre modèle est homoscédastique. Les régions de confiance entre les deux approches sont clairement différentes… meme si les prédictions sont presque superposées (oui oui, il y a deux courbes, une rouge et une bleue).

On Box-Cox transform in regression models

A few days ago, a former student of mine, David, came back to me about Box-Cox tests in linear models. It made me look more carefully at the test, and I do not understand what is computed, to be honest. Let us start with something simple, like a linear simple regression, i.e.

https://latex.codecogs.com/gif.latex?Y_i=\beta_0+\beta_1%20X_i+\varepsilon_i

Let us introduced – as suggested in Box & Cox (1964) – the following family of (power) transformations

https://latex.codecogs.com/gif.latex?Y_i^{(\lambda)}%20=%20\begin{cases}%20\dfrac{Y_i^\lambda-1}{\lambda}%20&\text{%20%20}%20(\lambda%20\neq%200)\\[8pt]%20\log{(Y_i)}%20%20&\text{%20}%20(\lambda%20=%200)%20\end{cases}

on the variable of interest. Then assume that

https://latex.codecogs.com/gif.latex?Y_i^{(\lambda)}=\beta_0+\beta_1%20X_i+\varepsilon_i

As mentioned in Chapter 14 of Davidson & MacKinnon (1993) – in French – the log-likelihood of this model (assuming that observations are independent, with distribution https://latex.codecogs.com/gif.latex?\varepsilon_i\sim\mathcal{N}(0,\sigma^2)) can be written

https://latex.codecogs.com/gif.latex?\log%20\mathcal{L}=-\frac{n}{2}\log(2\pi)-n\log\sigma%20\\-\frac{1}{2\sigma^2}\sum_{i=1}^n\left[Y_i^{(\lambda)}-(\beta_0+\beta_1%20X_i)\right]^2+(\lambda-1)\sum_{i=1}^n\log%20Y_i

We can then use profile-likelihood techniques (see here) to derive the optimal transformation.

This can be done in R extremely simply,

> library(MASS)
> boxcox(lm(dist~speed,data=cars),lambda=seq(0,1,by=.1))

we then get the following graph,

If we look at the code of the function, it is based on the QR decomposition of the https://latex.codecogs.com/gif.latex?\boldsymbol{X}matrix (since we assume that https://latex.codecogs.com/gif.latex?\boldsymbol{X} is a full-rank matrix). More precisely, https://latex.codecogs.com/gif.latex?\boldsymbol{X}=QR where https://latex.codecogs.com/gif.latex?\boldsymbol{X} is a https://latex.codecogs.com/gif.latex?n\times%202 matrix, https://latex.codecogs.com/gif.latex?Q is a https://latex.codecogs.com/gif.latex?n\times%202 orthonornal matrix, andhttps://latex.codecogs.com/gif.latex?R is a https://latex.codecogs.com/gif.latex?2\times2 upper triangle matrix. It might be convenient to use this matrix since, for instance, https://latex.codecogs.com/gif.latex?R\widehat{\boldsymbol{\beta}}=Q%27Y.  Thus, we do have an upper triangle system of equations.

> X=lm(dist~speed,data=cars)$qr

The code used to get the previous graph is (more or less) the following,

> g=function(x,lambda){
+ y=NA
+ if(lambda!=0){y=(x^lambda-1)/lambda}
+ if(lambda==0){y=log(x)}
+ return(y)} 
> n=nrow(cars)
> X=lm(dist~speed,data=cars)$qr
> Y=cars$dist
> logv=function(lambda){
+ -n/2*log(sum(qr.resid(X, g(Y,lambda)/
+ exp(mean(log(Y)))^(lambda-1))^2))}
> L=seq(0,1,by=.05)
> LV=Vectorize(logv)(L)
> points(L,LV,pch=19,cex=.85,col="red")

As we can see (with those red dots) we can reproduce the R graph. But it might not be consistent with other techniques (and functions described above). For instance, we can plot the profile likelihood function, https://latex.codecogs.com/gif.latex?\lambda\mapsto\log\mathcal{L}

> logv=function(lambda){
+ s=summary(lm(g(dist,lambda)~speed,
+ data=cars))$sigma
+ e=lm(g(dist,lambda)~speed,data=cars)$residuals
+ -n/2*log(2 * pi)-n*log(s)-.5/s^2*(sum(e^2))+
+ (lambda-1)*sum(log(Y))
+ }
> L=seq(0,1,by=.01)
> LV=Vectorize(logv)(L)
> plot(L,LV,type="l",ylab="")
> (maxf=optimize(logv,0:1,maximum=TRUE))
$maximum
[1] 0.430591

$objective
[1] -197.6966

> abline(v=maxf$maximum,lty=2)

The good point is that the optimal value of https://latex.codecogs.com/gif.latex?\lambda is the same as the one we got before. The only problem is that the https://latex.codecogs.com/gif.latex?y-axis has a different scale. And using profile likelihood techniques to derive a confidence interval will give us different results (with a larger confidence interval than the one given by the standard function),

> ic=maxf$objective-qchisq(.95,1)
> #install.packages("rootSolve")
> library(rootSolve)
> f=function(x)(logv(x)-ic)
> (lower=uniroot(f, c(0,maxf$maximum))$root)
[1] 0.1383507
> (upper=uniroot(f, c(maxf$maximum,1))$root)
[1] 0.780573
> segments(lower,ic,upper,ic,lwd=2,col="red")

Actually, it possible to rewrite the log-likelihood as

https://latex.codecogs.com/gif.latex?\mathcal{L}=\star-\frac{n}{2}\log\left[\sum_{i=1}^n\left(\frac{Y_i^{(\lambda)}-(\beta_0+\beta_1%20X_i)}{\dot{Y}^\lambda}\right)^2\right]

(let us just get rid of the constant), where

https://latex.codecogs.com/gif.latex?\dot{Y}=\exp\left[\frac{1}{n}\sum_{i=1}^n%20\log%20Y_i\right]

Here, it becomes

> logv=function(lambda){
+ e=lm(g(dist,lambda)~speed,data=cars)$residuals
+ elY=(exp(mean(log(Y))))
+ -n/2*log(sum((e/elY^lambda)^2))
+ }
>
> L=seq(0,1,by=.01)
> LV=Vectorize(logv)(L)
> plot(L,LV,type="l",ylab="")
> optimize(logv,0:1,maximum=TRUE)
$maximum
[1] 0.430591

$objective
[1] -47.73436

with again the same optimal value for https://latex.codecogs.com/gif.latex?\lambda, and the same confidence interval, since the function is the same, up to some additive constant.

So we have been able to derive the optimal transformation according to Box-Cox transformation, but so far, the confidence interval is not the same (it might come from the fact that here we substituted an estimator to the unknown parameter https://latex.codecogs.com/gif.latex?\sigma.

Why pictures are so important when modeling data?

(bis repetita) Consider the following regression summary,

Call:
lm(formula = y1 ~ x1)

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

Residual standard error: 1.237 on 9 degrees of freedom
Multiple R-squared: 0.6665,	Adjusted R-squared: 0.6295
F-statistic: 17.99 on 1 and 9 DF,  p-value: 0.00217
obtained from
> data(anscombe)
> reg1=lm(y1~x1,data=anscombe)
Can we say something if we look (only) at that output ? The intercept is significatively non-null, as well as the slope, the  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),
> reg2=lm(y2~x2,data=anscombe)
The output is here exactly the same as the one we had above
> summary(reg2)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Residual standard error: 0.003082 on 8 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 1.161e+06 on 1 and 8 DF,  p-value: < 2.2e-16
Finally consider
> reg4=lm(y4~x4,data=anscombe)
This time, there is an other kind of outlier, in ‘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

https://latex.codecogs.com/gif.latex?Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=\boldsymbol{X}_1\boldsymbol{\beta}_1+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon
Les équations normales sont ici (cf transparents du cours, ou équations écrites au tableau) https://latex.codecogs.com/gif.latex?\boldsymbol{X}%27\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\boldsymbol{X}%27Y, que l’on peut écrire sous la forme d’un système,

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\boldsymbol{X}_1%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_1%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_1%27Y%20\\%20\boldsymbol{X}_2%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1+\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
Le point intéressant est que si les variables dans le premier ensemble et les variables dans le second sont orthogonales, i.e. https://latex.codecogs.com/gif.latex?\boldsymbol{X}_1%27\boldsymbol{X}_2%20=\boldsymbol{0} et https://latex.codecogs.com/gif.latex?\boldsymbol{X}_2%27\boldsymbol{X}_1%20=\boldsymbol{0}, alors les termes diagonaux du système disparaissent, et on retrouve les équations normales des régressions simples, i.e.

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\boldsymbol{X}_1%27\boldsymbol{X}_1%20\widehat{\boldsymbol{\beta}}_1=\boldsymbol{X}_1%27Y%20\\%20\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
On a ici une conséquence du théorème de Frisch-Waugh: si les variables sont orthogonales, les estimateurs dans la régression multiple coïncident avec ceux obtenus avec des régressions simples.

Afin de mieux clarifier ce point, écrivons ce qui se passerait avec seulement deux variables. On va associer la constante à la première, i.e.

https://latex.codecogs.com/gif.latex?Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon=(1%20\%20\boldsymbol{X}_1)\begin{pmatrix}\boldsymbol{\beta}_0%20\\%20\boldsymbol{\beta}_1\end{pmatrix}+\boldsymbol{X}_2\boldsymbol{\beta}_2+\varepsilon
Dans ce cas, les équations normales sont

https://latex.codecogs.com/gif.latex?\left\{\begin{array}{l}\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}(1%20\%20\boldsymbol{X}_1)\begin{pmatrix}\widehat{\boldsymbol{\beta}}_0%20\\%20\widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\begin{pmatrix}1%20\\%20\boldsymbol{X}%27_1\end{pmatrix}Y%20\\%20\boldsymbol{X}_2%27(1\%20\boldsymbol{X}_1)%20\begin{pmatrix}\widehat{\boldsymbol{\beta}}_0%20\\%20\widehat{\boldsymbol{\beta}}_1\end{pmatrix}+\boldsymbol{X}_2%27\boldsymbol{X}_2%20\widehat{\boldsymbol{\beta}}_2=\boldsymbol{X}_2%27Y%20%20\end{array}\right.
Bref, il faut que https://latex.codecogs.com/gif.latex?\boldsymbol%20X_1 et https://latex.codecogs.com/gif.latex?\boldsymbol%20X_2 soient orthogonales pour que le terme sur la seconde diagonale disparaisse. C’est ce qu’on retrouve dans la plupart des livres d’économétrie, il suffit de chercher le théorème de Frisch-Waugh dans l’index (même si le théorème de Frisch-Waugh va beaucoup plus loin puisqu’il explique quels ajustement faire pour lier les estimateurs dans une régression simple, et dans le cas multiple).

Mais est-ce un résultat purement théorique ou peut-on l’invoquer dans un cas pratique ? Autrement dit, si deux variables semblent orthogonales, que se passent-il numériquement ? Si je simule des variables explicatives indépendantes, les estimateurs coïncident-ils (comme le suggère la théorie) ?

Considérons un modèle linéaire avec deux variables explicative que l’on supposera indépendantes, mais ayant toutes deux un impact sur la variable explicative. Le plus simple est de considérer un triple de variables telles que

https://latex.codecogs.com/gif.latex?\boldsymbol{Z}=\begin{pmatrix}%20Y%20\\%20X_1\\X_2\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix}%20\mu_Y%20\\%20\mu_{X_1}\\\mu_{X_2}\end{pmatrix},\begin{pmatrix}%201%20&%20\star%20&%20\star%20\\%20\star%20&%201&0\\\star&0&1\end{pmatrix}\right)
où https://latex.codecogs.com/gif.latex?\star désigne des valeurs non nulles. Par exemple,

> set.seed(1)
> n=100
> S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3)
> S
[,1] [,2] [,3]
[1,]  1.0  0.4  0.3
[2,]  0.4  1.0  0.0
[3,]  0.3  0.0  1.0
> M=c(1,3,6)
> Z=rmnorm(n,M,varcov=S)
> Y=Z[,1]; X1=Z[,2]; X2=Z[,3]
> base=data.frame(Y,X1,X2)
> cor(base)
Y          X1          X2
Y  1.0000000  0.37796580  0.27859381
X1 0.3779658  1.00000000 -0.05060539
X2 0.2785938 -0.05060539  1.00000000

On a simulé ici deux variables indépendantes, mais toutes deux corrélées avec la variable d’intérêt. On notera que numériquement

https://latex.codecogs.com/gif.latex?\frac{1}{n}\sum_{i=1}^n%20X_{1,i}X_{2,i}%20\neq%20\frac{1}{n}\sum_{i=1}^n%20X_{1,i}\cdot\frac{1}{n}\sum_{i=1}^n%20X_{2,i}

(autrement dit, en moyenne, on simulera des données indépendantes, mais dans ce cas précis, les données ne sont pas rigoureusement orthogonales).

> sum(X1*X2)/n
[1] 18.20194
> sum(X1)*sum(X2)/n^2
[1] 18.25097
Mais un test rapide nous pousse à accepter l’hypothèse d’indépendance entre les deux variables (sans trop d’hésitations).

> cor.test(X1,X2)

Pearsons product-moment correlation

data:  X1 and X2
t = -0.5016, df = 98, p-value = 0.6171
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2445923  0.1472766
sample estimates:
cor
-0.05060539

Si on regarde maintenant la régression,

> R=lm(Y~X1+X2,data=base)
> summary(R)

Call:
lm(formula = Y ~ X1 + X2, data = base)

Residuals:
Min       1Q   Median       3Q      Max
-1.87885 -0.57656  0.07957  0.52719  1.77407

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  4.71180    0.28246  16.682  < 2e-16 ***
X1           0.37232    0.08441   4.411 2.66e-05 ***
X2           0.25975    0.07755   3.350  0.00115 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7954 on 97 degrees of freedom
Multiple R-squared: 0.2317,	Adjusted R-squared: 0.2159
F-statistic: 14.63 on 2 and 97 DF,  p-value: 2.803e-06

Quel est le lien avec les régressions simples, sur chacune des variables ? Si on régresse sur la première, on obtient

> R1=lm(Y~X1,data=base)
> summary(R1)

Call:
lm(formula = Y ~ X1, data = base)

Residuals:
Min       1Q   Median       3Q      Max
-1.78973 -0.54060  0.04036  0.44138  1.94389

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.03165    0.27934  18.013  < 2e-16 ***
X1           0.35802    0.08859   4.041 0.000106 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8358 on 98 degrees of freedom
Multiple R-squared: 0.1429,	Adjusted R-squared: 0.1341
F-statistic: 16.33 on 1 and 98 DF,  p-value: 0.0001058

et si on régresse sur la seconde

> R2=lm(Y~X2,data=base)
> summary(R2)

Call:
lm(formula = Y ~ X2, data = base)

Residuals:
Min      1Q  Median      3Q     Max
-1.8652 -0.5507  0.1104  0.6109  2.1435

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.85053    0.12495  46.824  < 2e-16 ***
X2           0.24244    0.08443   2.872  0.00501 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.867 on 98 degrees of freedom
Multiple R-squared: 0.07761,	Adjusted R-squared: 0.0682
F-statistic: 8.246 on 1 and 98 DF,  p-value: 0.005005

Si les données ne sont pas rigoureusement orthogonales, alors https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 dans le modèle multiple est légèrement différent de https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 obtenu dans le modèle simple. Si on regarde la distribution de nos estimateurs (n’oublions pas qu’il y a une marge d’erreur), ils sont proches, ou au moins statistiquement proches (avec à gauche la distribution de https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 et à droite celle de https://latex.codecogs.com/gif.latex?\widehat{\beta}_2)

Si on génère 100 bases de données avec ce modèle (avec 500 observations à chaque fois) on obtient les estimateurs suivants (à gauche les https://latex.codecogs.com/gif.latex?\widehat{\beta}_1 et à droite les https://latex.codecogs.com/gif.latex?\widehat{\beta}_2, avec en abscisse le modèle multiple et en ordonnées, le modèle simple).

(la valeur centrale étant la vraie valeur, utilisée pour les simulations). On pourrait se demander si les valeurs éloignées du barycentre (ou de cette première bissectrice) correspondent à des cas où l’hypothèse d’indépendance entre les deux variables est rejetée ? On peut mettre des couleurs différentes en fonction de la https://latex.codecogs.com/gif.latex?p-value: plus la https://latex.codecogs.com/gif.latex?p-value est faible (rejet de l’hypothèse d’indépendance) plus la couleur sera vive (et rouge), et plus la https://latex.codecogs.com/gif.latex?p-value est forte (acceptation de l’hypothèse d’indépendance) plus la couleur sera terne (et jaune, voire blanche),

sur 1,000 simulations, on obtient

Le code est ici

> set.seed(1)
> library(mnormt)
> S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3)
> BX1=BX2=matrix(NA,1000,2)
> P=rep(NA,1000)
> for(s in 1:1000){
+   n=500
+   Z=rmnorm(n,c(6,3,1),varcov=S)
+   Y=Z[,1]; X1=Z[,2]; X2=Z[,3]
+   base=data.frame(Y,X1,X2)
+   R=lm(Y~X1+X2,data=base)
+   BX1[s,1]=coefficients(R)[2]
+   BX2[s,1]=coefficients(R)[3]
+   BX1[s,2]=coefficients(lm(Y~X1,data=base))[2]
+   BX2[s,2]=coefficients(lm(Y~X2,data=base))[2]
+   P[s]=cor.test(X1,X2)$p.value
+ }
> COL=(heat.colors(100))
> plot(BX1[,1],BX1[,2],pch=19,cex=.6,col=COL[trunc(100*P)])
> abline(a=0,b=1,lty=2)

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

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

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

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue")

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 – http://freakonometrics.blog.free.fr/public/perso5/logormal01.gif  – 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. http://freakonometrics.blog.free.fr/public/perso5/lognormal02.gif or http://freakonometrics.blog.free.fr/public/perso5/lognormal05.gif:

> 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 http://freakonometrics.blog.free.fr/public/perso5/lognormal03.gif or equivalently http://freakonometrics.blog.free.fr/public/perso5/lognormal04.gif.

> 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 http://freakonometrics.blog.free.fr/public/perso5/ols-002.gif était un estimateur Gaussien. En particulier, chacun des estimateurs est alors Gaussien, au sens où http://freakonometrics.blog.free.fr/public/perso5/ols--0003.gif, 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, http://freakonometrics.blog.free.fr/public/perso5/ols---ooo4.gif

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.

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

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 http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-01.gif of binary responses, i.e. http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-02.gif with http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-03.gif taking values 1 or 2. Assume that probability http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-04.gif 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 http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-06.gif such that http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-07.gif if http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-08.gif is lower than a given threshold, and 1 otherwise.
As in standard probit models, assume that

http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-09.gif

where we can assume that http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-10.gif is a Gaussian random vector. This assumption can be used to derive the likelihood of a sample http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-11.gif.

> 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

> (exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1)
[1] 0.5260951

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.

  1. Y[,1]==0)&(Y[,2]==1 []
  2. Y[,1]==1)&(Y[,2]==0 []
  3. Y[,1]==0)&(Y[,2]==0 []
  4. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==0 []
  5. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==0 []
  6. Y[,1]==0)&(Y[,2]==0)&(Y[,3]==1 []
  7. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==0 []
  8. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==1 []
  9. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==1 []
  10. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==1 []
  11. Y[,1]==0)&(Y[,2]==1 []
  12. Y[,1]==1)&(Y[,2]==0 []
  13. 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 http://freakonometrics.blog.free.fr/public/perso/wls08.png ce qui donne comme estimateur par moindres carrés pondérés

http://freakonometrics.blog.free.fr/public/perso/wls09.png

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 (), 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 ?