# From multinomial regression to binary classification on some Siamese data

There are two kinds of people in the world: people who think there are two kinds of people in the world and people who don’t

(borrowed from Menand (2018)). Because things are always simpler when we face only binary choice, aren’t they? But consider here the case were multiple options are possible, and let us see if we cannot get back to simpler binary choices. Consider a collection of observations $(y_i,\boldsymbol{x}_i)$ where $y_i$ is some categorical variable, $y_i\in\mathcal{A}$ where $\mathcal{A}=\lbrace A_1,\cdots,A_\kappa \rbrace$, with $\kappa$ possible categories. Let $\mathcal{I}_k=\lbrace i:y_i\in A_k \rbrace$.

In a classical multinomial logistic regression, suppose that $A_1$ is the reference, then $$\mathbb{P}[Y=A_j|\boldsymbol{X}=\boldsymbol{x}]=\frac{\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_j]}{1+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_2]+\cdots+\exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k]}$$With a lot a categories, and a small number of observations, inference can be complicated, and non-robust.

• the Siamese dataset

The name Siamese I use, here, comes from Siamese Networks. Or sort of… As we say in French, it is an « histoire de l’homme qui a vu l’homme qui a vu l’ours » (story of the man who saw the man who saw the bear). A few years ago, a student tried to explain to me the idea of Siamese Networks and this is what I understood. I might be completely wrong, but the idea I got from it did make sense, in my mind at least. That is the story of that blog post…

The idea of the siamese algorithm will be to consider all pairs of observations, $(y_i,\boldsymbol{x}_i)$ and $(y_j,\boldsymbol{x}_j)$ :

1. $\tilde y_{i,j}=\boldsymbol{1}(y_i=y_j)$ indicating if individuals $i$ and $j$ are in the same category
2. $\tilde{\boldsymbol{x}}_{i,j}$ is a collection of $p-1$ variables,
• $\tilde {x}_{k:i,j}={x}_{k:i}-{x}_{k:j}$ if $x_k$ is continuous, or $\tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}|$ (we can use another metric, e.g. $\tilde {x}_{k:i,j}=|{x}_{k:i}-{x}_{k:j}|^2$, and this is why I decided to use some GAM model in the logistic regression on the Siamese dataset)
• $\tilde {x}_{k:i,j}=({x}_{k:i},{x}_{k:j})\in\mathcal{X}_k\times\mathcal{X}_k$ if $x_k$ is a categorical variables (taking values in the set $\mathcal{X}_k$), or $\tilde {x}_{k:i,j}=\boldsymbol{1}({x}_{k:i}\neq{x}_{k:j})\in\{0,1\}$

The original dataset was a $n\times p$ matrix, and (if there are no categorical variable), it becomes a $n(n-1)/2\times p$ matrix. The key point is that if the original variable $y_i$ was multinomial, $y_{i,j}$ is now binomial. For instance, if our initial dataset was the following, with two covariates, one continuous and one categorical its siamese counterpart is the following • Classification step

On the dataset $(\tilde y_{i,j},\tilde {x}_{k:i,j})_{i,j}$, fit a logistic regression, $$\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}]=\frac{\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}^\top\boldsymbol{\beta}]}$$(or any classification model – CART, random forest, etc). But that is the easy part (unless $n$ is large, because the siamese dataset has (roughly) $n^2/2$ rows). The difficult task is the prediction

• Prediction step

Consider a new input variable $\boldsymbol{x}_{\cdot}$, and define its siamese version, $\tilde{\boldsymbol{x}}_{\cdot}=(\tilde{\boldsymbol{x}}_{\cdot,j})_j$, i.e. a database with $n$ rows. Then compute
$$p_{\cdot,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{\cdot,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{\cdot,j}^\top\boldsymbol{\beta}]}$$ where $p_{\cdot,j}$ is the probability that $(y_j,\boldsymbol{x}_{j})$ and $(y_{\cdot},\boldsymbol{x}_{\cdot})$ are in the same category, as well as $$p_{i,j}=\mathbb{P}[\tilde Y|\tilde{\boldsymbol{X}}=\tilde{\boldsymbol{x}}_{i,j}]=\frac{\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}{1+\exp[\tilde{\boldsymbol{x}}_{i,j}^\top\boldsymbol{\beta}]}$$
Let $\boldsymbol{p}_{\cdot}=(p_{\cdot,j})$, and similarly $\boldsymbol{p}_{i}=(p_{i,j})$. Then several techniques can be used to predict $y_{\cdot}$.

1. $\widehat{y}_{\cdot}=y_{j^\star}$ where ${j^\star}=\underset{j=1,\cdots,n}{\text{argmax}}\{p_{\cdot,j}\}$: the predicted class is the one of the observation the most likely to be other same class
2. $\widehat{y}_{\cdot}=y_{j^\star}$ where ${j^\star}=\underset{\ell=1,\cdots,k}{\text{argmax}}\{\overline{p}_{\ell}\}$, where$$\overline{p}_{\ell} = \frac{1}{n_{\ell}}\sum_{j\in\mathcal{I}_s} \boldsymbol{1}(y_j =y_{\ell}),\text{ where }\mathcal{I}_s=\lbrace i:p_{\cdot,i}>s\rbrace$$consider only probabilities sufficiently high to be considered, and predict the most important class (majority rule)
3. $\widehat{y}_{\cdot}=y_{j^\star}$ where $${j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{\theta_i\}$$ where $$\theta_i=\cos(\boldsymbol{p}_{\cdot},\boldsymbol{p}_{i})=\displaystyle{\frac{\boldsymbol{p}_{\cdot}\cdot\boldsymbol{p}_{i}}{\|\boldsymbol{p}_{\cdot}\|\|\boldsymbol{p}_{i}\|}}$$
4. $\widehat{y}_{\cdot}=y_{j^\star}$ where $${j^\star}=\underset{i=1,\cdots,n}{\text{argmax}}\{KL_{\cdot|i}\}$$ and $KL_{\cdot|i}=\displaystyle{\sum_{j=1}^n p_{\cdot,j}\log\frac{p_{\cdot,j}}{p_{i,j}}}$ (but one can select another metric)
5. $\widehat{y}_{\cdot}=y_{j^\star}$ where ${j^\star}=\underset{j\in\mathcal{J}}{\text{argmax}}\{p_{\cdot,j}\}$ and $\mathcal{J}$ is a sample of $k$ observations, chosen randomly, one in each group (one-shot procedure): the predicted class is the one of the observation the most likely to be other same class

Heuristically, it can be related to some $k$ nearest neighbors strategy: we give the attribute that most neighbors have. The total distance is a weighted sum of the componentwise distances (for the logistic regression).

• Simulation study

In order to test that technique, let us generate some multinomial model where $y$ has 10 possible labels, with 6 (independent) covariates $x_1,\cdots,x_6$, and $\mathbb{P}[Y=A_k|\boldsymbol{X}=\boldsymbol{x}]\propto \exp[\boldsymbol{x}^\top\boldsymbol{\beta}_k]$ (where coefficients $\boldsymbol{\beta}_k$ where generated randomly) for $k\in\{1,2,\cdots,10\}$ (there were 10 categories) and with $n=700$ observations.

n=700 X1=rnorm(n) X2=rnorm(n) X3=rnorm(n) X4=rnorm(n) X5=rnorm(n) X6=rnorm(n) X=cbind(1,X1,X2,X3,sqrt(abs(X4)),X5*X1,X6) k = 10 PARAM = matrix(rnorm(k*6),k,6) PARAM[,1]=PARAM[,1]-1 PARAM=cbind(PARAM,0) P=matrix(NA,n,k-1) for(j in 1:(k-1)) P[,j] = X %*% (PARAM[j,])+rnorm(n) P=cbind(P,0) S=apply(exp(P),1,sum) Pb = exp(P)/S tirage = function(i){ sample(1:10,size=1,prob = Pb[i,]) } Y = LETTERS[Vectorize(tirage)(1:n)] dbase = data.frame(Y=as.factor(Y),X1,X2,X3,X4,X5)

In the paragraph previously, I suggested to take the most likely one. Being wrong means that it was not the first choice. But perhaps being the second or the third choice is not that bad, actually. So in my simulations, I look at the proportion of predictions were our prediction is the good one (top 1), if the true one is either the most likely or the second most likely (top 2), or in the top 3. That will be on my $x$-axis. I draw some line, but we simply have three points (top 1, top 2 and top 3). I compute the proportion of good prediction, using cross-validation techniques (10-fold). The black lines are one of the methods described above. The red one is the standard multinomial model (with a logistic link function). For the Siamese model, I tried several models. I tried a logistic regression, and some smooth version (GAM) on top and a classification tree, on the left, as well as some random forest on the right, below. It looks like the multinomial approach performs always better than any Siamese one… and to be honest, I am disappointed.

Here is the code I did use when I considered a logistic regression on the Siamese dataset,

set.seed(1) kfold = sample(rep(1:10,n/10))   KFOLDglm = function(i){ i_test=which(kfold==i) i_calibration=which(kfold!=i) y=credit[i_calibration,"Y"] tirage = function(){ v=c(sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1)) names(v)=levels(y) return(v) }   LogisticModel &lt;- multinom(Y ~ ., data = credit[i_calibration,], trace=FALSE)   comparaisonx = function(base,x=base[1,]){ mix_base = base for(j in 1:ncol(base)){ xj = as.numeric(x[j])-base[,j] mix_base[,j] = (xj) } mix_base } comparaisony = function(base,y=base){ as.factor(base == y) } creditx = credit[,-which(names(credit) == "Y")] nc=length(i_calibration) B=comparaisonx(base = creditx[i_calibration[2:nc],],x=creditx[i_calibration,]) B$Y=comparaisony(base = credit[i_calibration[2:nc],"Y"],y=credit[i_calibration,"Y"]) for(i in 2:(nc-1)){ B0=comparaisonx(base = creditx[i_calibration[(i+1):nc],],x=creditx[i_calibration[i],]) B0$Y=comparaisony(base = credit[i_calibration[(i+1):nc],"Y"],y=credit[i_calibration[i],"Y"]) B=rbind(B,B0) } credit_mix = B   OneShotLogisticModel &lt;- glm(Y ~ ., data = credit_mix, family=binomial) A_ref = table(credit[i_calibration,"Y"])/length(i_calibration)   vect_oneshot = function(i){ B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,]) predict(OneShotLogisticModel,type="response",newdata=B2) }   prediction_oneshot = function(i,type=1){ B2=comparaisonx(base = creditx[i_calibration,],x=creditx[i,]) p=predict(OneShotLogisticModel,type="response",newdata=B2) y=credit[i_calibration,"Y"] base = data.frame(p,y) base = base[rev(order(base$p)),] if(type==1){T = table(base$y[1:11]) return(names(which.max(T)))} if(type==2){return(base$y)} if(type==3){A=table(base$y[1:10])/10 T=A/A_ref return(names(which.max(T)))} if(type==4){ costheta = rep(NA,length(i_calibration)) for(j in 1:length(i_calibration)){ vecteur_proba = vect_oneshot(i_calibration[j]) costheta[j] = sum(vecteur_proba*p)/(sqrt(sum(vecteur_proba^2))*sqrt(sum(p^2))) } return(y[which.max(costheta)])} if(type==5){ kl = rep(NA,length(i_calibration)) for(j in 1:length(i_calibration)){ vecteur_proba = vect_oneshot(i_calibration[j]) kl[j] = sum(p*log(vecteur_proba/p)) } return(y[which.max(as.vector(kl))])} if(type==6){ ## one shot : tirer au hasard un de chaque, et dire lequel est plus credible !   y=credit[i_calibration,"Y"] tirage = function(){ v=c(sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1), sample(i_calibration[y==levels(y)],size=1)) names(v)=levels(y) return(v) } pd=rep(NA,101) for(ix in 1:101){ ids = tirage() B2=comparaisonx(base = creditx[ids,],x=creditx[i,]) p=predict(OneShotLogisticModel,type="response",newdata=B2) pd[ix]=levels(y)[which.max(p)] } levels(y)[which.max(table(pd))] } } PRED0=as.character(credit[i_test,"Y"]) PRED1=as.character(predict(LogisticModel,type = "class", newdata=credit[i_test,])) PRED21=as.character(Vectorize(function(i) prediction_oneshot(i,type=1)) (i_test)) PRED22=as.character(Vectorize(function(i) prediction_oneshot(i,type=2))(i_test)) PRED23=as.character(Vectorize(function(i) prediction_oneshot(i,type=3))(i_test)) PRED24=as.character(Vectorize(function(i) prediction_oneshot(i,type=4))(i_test)) PRED25=as.character(Vectorize(function(i) prediction_oneshot(i,type=5))(i_test)) PRED26=as.character(Vectorize(function(i) prediction_oneshot(i,type=6))(i_test)) B=data.frame(PRED0,PRED1,PRED21,PRED22,PRED23,PRED24,PRED25,PRED26) B }   s=1/100;setTxtProgressBar(pb, s*2) for(i in 2:10){ PREDICTION = rbind(PREDICTION,KFOLDglm(i)) s=s+1/100;setTxtProgressBar(pb, s*2)} for(j in 1:5) PREDICTION[,j]=as.character(PREDICTION[,j]) L=list() v=mean(PREDICTION[,1]!=PREDICTION[,2]) names(v)="logistic" L[["logistic"]]=v v=c(mean(PREDICTION[,1]!=PREDICTION[,3]), mean(PREDICTION[,1]!=PREDICTION[,4]), mean(PREDICTION[,1]!=PREDICTION[,5]), mean(PREDICTION[,1]!=PREDICTION[,6]), mean(PREDICTION[,1]!=PREDICTION[,7]), mean(PREDICTION[,1]!=PREDICTION[,8])) names(v)=c("top10","max","10norm","cos","KL","OS") L[["glm"]]=v

Cette fin de semaine, je mettais en ligne un court billet du deuxième effet kiss-cool où je rappelais que quand on fait une régression sur plusieurs variables explicatives corrélées, ce n’est pas équivalent à faire plusieurs régressions simple. C’est ce que disait le théorème de Frisch-Waugh, publié il y a près de 90 ans. Mais je glissais qu’il était possible de faire des régressions en cascade, sans aller beaucoup plus loin. Il faut relire l’article de Michael Lovell qui donner une interprétation géométrique de ce dessin. Mais programmons le, pour le voir… Considérons une base avec 3 variables explicatives, correspondant à des nombres d’incendies à Chicago

chicago=read.table("http://freakonometrics.free.fr/chicago.txt", header=TRUE,sep=";")

La régression multiple s’écrit ici$$y_i=\beta_0+\beta_1 x_{1,i}+\beta_2 x_{2,i}+\beta_3 x_{3,i}+\varepsilon_i$$On peut estimer ces paramètres par moindre carrés,

reg=lm(Fire~.,data=chicago) summary(reg)   Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 22.07525 6.19447 3.564 0.000910 *** X_1 -0.62764 5.28130 -0.119 0.905953 X_2 0.22378 0.06161 3.632 0.000744 *** X_3 -1.55059 0.38195 -4.060 0.000204 ***

et on a alors la prévision suivante
$$\widehat{y}_i=\widehat{\beta}_0+\widehat{\beta}_1 x_{1,i}+\widehat{\beta}_2 x_{2,i}+\widehat{\beta}_3 x_{3,i}$$En fait, je peux même visualiser cette prévision.

BC = cbind(x0=coefficients(reg), x1=coefficients(reg)*chicago[,2], x2=coefficients(reg)*chicago[,3], x3=coefficients(reg)*chicago[,4], yp=predict(reg),y=chicago$Fire) colr = 1:4 dessin_fire = function(i,BC=BC){ plot(0:3,0:3,xlim=c(0,35),col="white",ylim=c(-1,3),xlab="",ylab="",axes=FALSE) abline(v=BC[i,1],lty=2) rect(BC[i,1],3,BC[i,1]+BC[i,2],2,col=colr,border=NA) arrows(BC[i,1],2.5,BC[i,1]+BC[i,2],2.5,lwd=2,length=.1,col="white") rect(BC[i,1]+BC[i,2],2,BC[i,1]+BC[i,2]+BC[i,3],1,col=colr,border=NA) arrows(BC[i,1]+BC[i,2],1.5,BC[i,1]+BC[i,2]+BC[i,3],1.5,lwd=2,length=.1,col="white") rect(BC[i,1]+BC[i,2]+BC[i,3],1,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],0,col=colr,border=NA) segments(BC[i,5],0,BC[i,5],1,lwd=3) arrows(BC[i,1]+BC[i,2]+BC[i,3],.5,BC[i,1]+BC[i,2]+BC[i,3]+BC[i,4],.5,lwd=2,length=.1,col="white") abline(v=BC[i,1],lty=2) rect(BC[i,5],-1,BC[i,6],0,col=colr,density=20,border=NA) segments(BC[i,6],0,BC[i,6],-1,lwd=3,col=colr) text(26,2.5,expression(X),pos=4,col=colr) text(26,1.5,expression(X),pos=4,col=colr) text(26,.5,expression(X),pos=4,col=colr) text(26,-.5,expression(epsilon),pos=4,col=colr) axis(1)} Par exemple, pour la douzième observation de ma base dessin_fire(12) on obtient On lit de haut en bas : on commence par constante, vers 22. Puis on rajoute trois fois rien, parce que $x_1$ était non-significative dans notre régression. Puis on rajoute un petit quelque chose à cause de $x_2$, pour passer à 25, et $x_3$ nous fait plonger de plus de 20 points, pour finir vers 3.67.Le trait noir est la prévision que l’on obtient, à partir des trois variables. En bleu, on peut même voir l’erreur. On peut noter que $x_1$ est non-significative, tout en étant relativement corrélée avec $y$ cor(chicago[,1:2]) Fire X_1 Fire 1.0000000 0.3773486 X_1 0.3773486 1.0000000 Car si $x_1$ est non-significative, c’est parce que cette variable est très corrélée avec une autre variable explicative. En fait, on peut faire une série de régressions en cascade, en commençant par $x_1$, et plus spécifiquement, on peut écrire$$\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_1 x_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,2}+\widehat{b}_2 \tilde{x}_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}$$où le premier terme $(1)$ est obtenu en régression simplement $y$ sur $x_1$ (oui, on fait juste une régression simple), $$y_i=b_0+b_1 x_{1,i}+\eta_i$$puis on apporte une petite correction, pour tenir compte de ce que nous apprend $x_2$ une fois la première régression effectuée. C’est exactement ce que raconte le théorème de Frish-Waugh : on va projeter $y$ sur $x_1^{\perp}$ (car c’est justement ce que $x_1$ n’explique pas), et $x_2$ sur $x_1^{\perp}$ et faire la régression de ces deux projections, soit, pour $(2)$ $$\Pi_{x_1^{\perp}}y_i=b_{0,2}+b_2 \Pi_{x_1^{\perp}}x_{2,i}+\eta_{2,i}$$puis pour $(3)$ $$\Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i}$$En terme informatique, cela donne reg1=lm(Fire~X_1,data=chicago) reg2=lm(residuals(lm(Fire~X_1,data=chicago))~residuals(lm(X_2~X_1,data=chicago))) reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago))) soit ici BC123=cbind(x0=coefficients(reg1), x1=coefficients(reg1)*chicago[,2], x2=coefficients(reg2)+coefficients(reg2)*residuals(lm(X_2~X_1,data=chicago)), x3=coefficients(reg3)+coefficients(reg3)*residuals(lm(X_3~X_1+X_2,data=chicago)), yp=predict(reg),y=chicago$Fire)

En adaptant la fonction de dessin précédante on obtient

dessin_fire_123(12) Autrement dit, on commence avec une régression simple, en rouge: on part d’une constante vers 3.5, puis on augemente notre prévision, en tenant compte de $x_1$. Puis on corrige avec $x_2$ (ou plutôt la partie de $x_2$ non expliquée par $x_1$). sur ce qui n’avait pas été expliqué auparavant, ce qui pousse à revoir un peu à la baisse notre prévision, puis on tient compte de $x_3$ (ou là encore, sa projection sur l’orthogonal de $x_1$ et $x_2$). On note que la prévision est très exactement la même qu’auparavant.

Mais on peut aller plus loin… pourquoi commencer par $x_1$? On peut aussi commencer par $x_2$. On considère alors$$\widehat{y}_i=\underbrace{\widehat{b}_0+\widehat{b}_2 x_{2,i}}_{(2)}+\underbrace{\widehat{b}_{0,1}+\widehat{b}_1 \tilde{x}_{1,i}}_{(1)}+\underbrace{\widehat{b}_{0,3}+\widehat{b}_3 \tilde{x}_{3,i}}_{(3)}$$où le premier terme $(2)$ est obtenu en régression simplement $y$ sur $x_2$$$y_i=b_0+b_2 x_{2,i}+\eta_i$$puis on apporte une petite correction, pour tenir compte de ce que nous apprend $x_1$ une fois la première régression effectuée, $(1)$ $$\Pi_{x_2^{\perp}}y_i=b_{0,1}+b_1\Pi_{x_2^{\perp}}x_{1,i}+\eta_{1,i}$$puis pour $(3)$ $$\Pi_{(x_1,x_2)^{\perp}}y_i=b_{0,3}+b_3 \Pi_{(x_1,x_2)^{\perp}}x_{3,i}+\eta_{2,i}$$ Soit

reg1=lm(Fire~X_2,data=chicago) reg2=lm(residuals(lm(Fire~X_2,data=chicago))~residuals(lm(X_1~X_2,data=chicago))) reg3=lm(residuals(lm(Fire~X_1+X_2,data=chicago))~residuals(lm(X_3~X_1+X_2,data=chicago))) BC213=cbind(x0=coefficients(reg1), x1=coefficients(reg1)*chicago[,3], x2=coefficients(reg2)+coefficients(reg2)*residuals(lm(X_1~X_2,data=chicago)), x3=coefficients(reg3)+coefficients(reg3)*residuals(lm(X_3~X_1+X_2,data=chicago)), yp=predict(reg),y=chicago$Fire) et visuellement, on obtient dessin_fire_213(12) Cette fois, la constante est un peu plus élevée, et on commence par faire une régression seulement avec $x_2$, pour arriver un peu en dessous de 10. Puis on corrige. On note que cette seconde correction nous ramène… très exactement comme dans le cas précédant. Mais si on pense en projection successives, on ne doit pas être surpris. On peut récapituler avec une autre prévision, soit avec le modèle de régresion multiple (mais là encore, $x_1$ n’explique pas grand chose) dessin_fire(15) soit en commencant par une régression simple sur $x_1$ puis en faisant des régression en cascade, dessin_fire_123(15) soit, alternativement, en commencant par une régression simple sur $x_2$ dessin_fire_213(15) Comme annoncé, cette trois approches sont équivalentes, et donnent très exactement la même prévision. # Regression discontinuity model for TV series In September, we are usually happy to see our favorite TV series back on air… Or not? Because, admit it, if we are happy to see those characters back, most of the time, we are disappointed, too. So why not look at the data, to confirm this feeling? Nazareno Andrade shared some nice codes to get IMDB ratings in a nice csv file (you can either use the large csv file, or run your own codes) download.file("https://github.com/nazareno/imdb-series/raw/master/data/series_from_imdb.csv", destfile="series_from_imdb.csv") base = read.csv("series_from_imdb.csv") It is a large dataset, with more than 64,000 episodes of almost 890 TV series, str(base) 'data.frame': 64018 obs. of 18 variables:$ series_name: Factor w/ 889 levels "'Allo 'Allo!",..: 137 137 137 137 137 137 137 137 137 137 ... $episode : Factor w/ 54090 levels "-30-","¡Viva los muertos!",..: 32314 7446 16 7176 17748 9562 1379 36218 17845 5553 ...$ series_ep : int 1 2 3 4 5 6 7 8 9 10 ... $season : int 1 1 1 1 1 1 1 2 2 2 ...$ season_ep : int 1 2 3 4 5 6 7 1 2 3 ... $user_rating: num 8.9 8.7 8.7 8.2 8.3 9.2 8.8 8.7 9.2 8.3 ... Just pick a TV series, for instance Dan Harmon’s Community, sbase = base[base$series_name=="Community",]

We can plot the evolution of the rating over the 110 episodes.

sbase=sbase[!duplicated(sbase[,c(1,2,4,5)]),] sbase$series_ep=1:nrow(sbase) ()since there could be some problem with the data (such as duplicates, let us clean it quickly) plot(sbase$series_ep,sbase$UserRating,xlab=sbase$series_name) idx=c(0,which(diff(sbase$season)!=0),nrow(sbase)) abline(v=idx+.5,lty=2,col=colr) a = unique(sbase$season) for(u in a){ ssbase = sbase[sbase$season==u,] reg = lm(UserRating~series_ep,data=ssbase) lines(ssbase$series_ep,predict(reg),col=colr,lwd=2) } The vertical lines are here to visualize the seasons. On issue is that the lenght can vary with time. Consider Linwood Boomer’s Malcom in The Middle,

sbase = base[base$series_name=="Malcolm in the Middle",] or Craig Thomas and Carter Bays’s How I Met Your Mother, sbase = base[base$series_name=="How I Met Your Mother",] On those two, the evolution is rather stable. Look at AMC’s The Walking Dead,

sbase = base[base$series_name=="The Walking Dead",] Now, look at Howard Gordon and Alex Gansa’s Homeland, sbase = base[base$series_name=="Homeland",] There is an issue here with the last episode of season4, “Long Time Coming“, that has a very poor rating. If we remove that point, we get the thin line. Note that the regression line is always increasing. For Michael Hirst’s Vickings, we have

sbase = base[base$series_name=="Vicking",] If we look more carefully on the previous graph, for five seasons (out of six), we have a positive slope. Well, to be honest, it is not significantly positive most of the time, but still. Out of 80 shows, and a total of 583 seasons, the slope is postive 75% of the time (433) and negative 25% of the time (150). BASE = NULL L80 = unique(base$series_name) for(j in 1:length(L)){ sbase=base[base$series_name==L[j],] sbase=sbase[!duplicated(sbase[,c(1,2,4,5)]),] sbase=sbase[sbase$season&gt;0,] sbase$series_ep=1:nrow(sbase) a=unique(sbase$season) a=a[!is.na(a)] for(u in a){ ssbase=sbase[sbase$season==u,] reg=lm(UserRating~series_ep,data=ssbase) pente = NA if((!is.na(coefficients(reg)))&amp;(!is.na((summary(reg)$coefficients[2,4])))){ if((summary(reg)$coefficients[2,4]&lt;.05)&amp;(coefficients(reg)&gt;0)) pente="positive" if((summary(reg)$coefficients[2,4]&lt;.05)&amp;(coefficients(reg)&lt;0)) pente="negative" sdf=data.frame(nom=sbase$series_name,season=u,slope=coefficients(reg),inf=confint(reg)[2,1],sup=confint(reg)[2,2],signe=pente) BASE=rbind(BASE,sdf)} }} str(BASE) 'data.frame': 583 obs. of 6 variables:$ nom : Factor w/ 80 levels "Friends","Game of Thrones",..: 1 1 1 1 1 1 1 1 1 1 ...   mean(BASE$slope&gt;0)  0.7427101 table(BASE$signe) negative positive 15 144

Most of the time, the slope is not significant. To be more specific, 72% of the time, the slope is not significant. But when it is, 90% of the time, it is positive (144 seasons). Let us look at other TV series, for instance Joel Surnow and Robert Cochran’s 24,

sbase = base[base$series_name=="24",] Álex Pina’s La Casa de Papel, sbase = base[base$series_name=="La Casa de Papel",] Steven Knight’s Peaky Blinders,

sbase = base[base$series_name=="Peaky Blinders",] or David Simon’s The Wire, sbase = base[base$series_name=="The Wire",] The slope is increasing over almost all seasons. But a major drawback is that when we get back to our show, for a new season, we usually get disapointed. More specifically, we can quantify the difference in red below that can be estimated using

sbase12 = sbase[sbase$season%in%c(a[ij],a[ij+1]),] seuil = sbase12$series_ep[which(diff(sbase12$season)!=0)]+.5 s = function(x) (x-seuil)*(x&gt;seuil) reg = lm(UserRating~series_ep+s(series_ep)+I(series_ep&gt;seuil),data=sbase12) Here we have summary(reg) Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 8.45000 0.16338 51.719 2e-16 *** series_ep 0.10000 0.03235 3.091 0.008598 ** s(series_ep) 0.02000 0.04218 0.474 0.643291 I(series_ep)TRUE. -1.01778 0.20486 -4.968 0.000257 *** so the drop of 1 point (out of 10) cannot be claimed as being significant. That is the idea of regression discontinuity. If we loop again over all our series, we have 485 pairs of consecutive seasons. As expected, in 75% of the casse, from season $t-1$ to season $t$, we observe a negative rupture. As previously, in 70% of the cases, it is not significat (with linear models before and after), and when it is significant, it is negative in 96% of the cases ! But an alternative can be to use nonparametric models, on both sides. To illustrate, consider David Benioff and D. B. Weiss’s Game of Thrones, sbase = base[base$series_name=="Game of Thrones",] but let us remove the last season (no spoiler here, but clearly not worst watching) Consider for instance the drop between season 1 and season 2,

library(rdd) sbase12=sbase[sbase$season%in%c(1,2),] lmr=RDestimate(UserRating~series_ep,data=sbase12,cutpoint=mean(range(sbase12$series_ep))) plot(lmr) This is very consistent with what we observed with our linear regressions actually,

seuil=10.5 s = function(x) (x-seuil)*(x&gt;seuil) reg = lm(UserRating~series_ep+s(series_ep)+I(series_ep&gt;seuil),data=sbase12) summary(reg)   Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 8.70000 0.15458 56.281 2e-16 *** series_ep 0.07273 0.02491 2.919 0.01003 * s(series_ep) 0.01455 0.03523 0.413 0.68520 I(series_ep)TRUE -0.94000 0.20316 -4.627 0.00028 *** Here, the drop of one point is significant…

So, your favorite show had an outstanding finale ? and you can’t wait to watch the new season… Well, statistically, it’s very likely that you will be disapointed by the first episode of the forthcoming season…

# Du deuxième effet kiss-cool (régression multiple, scoring et évaluation)

Lorsque j’étais petit (il y a fort longtemps, à une époque où je regardais pas mal la télévision) il y avait une publicité pour les pastilles kiss cool,

Et quand je présente la régression multiple à mes étudiants, je ne peux m’empêcher d’y penser… Mais avant d’aller plus loin sur le parallèle, faisons un peu de mathématiques.

Les techniques de régression permettent d’avoir des jolis théorèmes, souvent d’une portée incroyablement générale (moyennant quelques petites hypothèses techniques). Par exemple le théorème de Frisch-Waugh, en régression multiple, dont j’ai déjà parlé dans des vieux billets. Un des corolaires est que lorsque les variables explicatives dans un modèle de régression sont orthogonales, la régression multiple correspond à une collection de régressions simples (autrement dit, les estimateurs par moindres carrés coïncident). Formellement, si on considère le modèle$$y_i=\beta_0+\beta_1x_{1,i}+\beta_2x_{2,i}+\varepsilon_i$$(avec les hypothèses usuelles des modèles de régression) alors, si les variables $x_1$ et $x_2$ sont non-corrélées, $\widehat{\beta}_1$ coïncide avec $\widehat{b}_1$ dans le modèle$$y_i=b_0+b_1x_{1,i}+\eta_i$$On peut faire une petite simulation pour confirmer (sinon, bien entendu, on peut regarder la démonstration qui se trouve dans tous les livres d’économétrie, qui est d’ailleurs un simple résultat d’algèbre linéaire, ou de géométrie, avec des projections successives sur des sous-espaces orthogonaux – même si c’est à Michael Lovell que l’on doit l’approche géométrique).

library(mnormt) r = 0 S = matrix(c(1,r,r,1),2,2) n = 1000 set.seed(1) X = rmnorm(n,c(0,0),S) E = rnorm(n,0,.3) Y = 2+X[,1]-2*X[,2]+E base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2]))

Petite note technique: je vais centrer les variables, histoire de ne pas avoir à garder la constante dans mon modèle (qui va compliquer les notations, et potentiellement embrouiller un peu le billet). La constante est nulle ici, on le voit,

reg = lm(Y~X1+X2,data=base) summary(reg)   Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 7.449e-18 9.777e-03 0.0 1 X1 1.012e+00 9.171e-03 110.3 2e-16 *** X2 -1.988e+00 9.719e-03 -204.6 2e-16 ***

Bref, je régresse sans constante

reg = lm(Y~0+X1+X2,data=base) summary(reg)   Coefficients: Estimate Std. Error t value Pr(|t|) X1 1.011520 0.009166 110.3 2e-16 *** X2 -1.988321 0.009714 -204.7 2e-16 ***

Maintenant, on va regarder les deux régressions simple

reg1 = lm(Y~0+X1,data=base) summary(reg1)   Coefficients: Estimate Std. Error t value Pr(|t|) X1 1.01300 0.06006 16.86 2e-16 ***

quand on régresse juste sur la première variable, et pour la seconde, on obtient

reg2 = lm(Y~0+X2,data=base) summary(reg2)   Coefficients: Estimate Std. Error t value Pr(|t|) X2 -1.98916 0.03528 -56.39 2e-16 ***

Autrement dit, nos estimateurs sont très proches (si on avait laissé la constante, ils coïncideraient).

Maintenant, le gros soucis est que ce résultat n’est plus valide lorsque les variables explicatives sont corrélées. Le théorème de Frisch-Waugh nous explique comment ces estimateurs divergent, mais le point ici est que si les variables sont corrélées, utiliser les régressions simples donne deux estimateurs biaisés des vrais paramètres (du modèle multiple). Recommençons l’exercice précédant

r=.9 S=matrix(c(1,r,r,1),2,2) set.seed(1) X=rmnorm(n,c(0,0),S) Y = 2+X[,1]-2*X[,2]+E base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2])) reg = lm(Y~X1+X2,data=base) summary(reg) reg = lm(Y~0+X1+X2,data=base) summary(reg)   Coefficients: Estimate Std. Error t value Pr(|t|) X1 0.98740 0.02205 44.79 2e-16 *** X2 -1.97321 0.02229 -88.54 2e-16 ***

(on retrouve des valeurs proches de celles utilisées pour simuler nos données, donc tout va bien). En revanche, pour les régressions simples, on obtient

reg1 = lm(Y~0+X1,data=base) summary(reg1)   Coefficients: Estimate Std. Error t value Pr(|t|) X1 -0.78784 0.02726 -28.9 2e-16 ***

et

reg2 = lm(Y~0+X2,data=base) summary(reg2)   Coefficients: Estimate Std. Error t value Pr(|t|) X2 -1.06543 0.01607 -66.31 2e-16 ***

Autrement dit, $\widehat{b}_1\neq \widehat{\beta}_1$ (et pareil pour le second). Ce qui signifie que si on construit une prévision à partir de ce modèle, $\widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}$, on sera potentiellement très loin de la “bonne” prévision $\widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i}$ (qui sera sans biais, etc, je renvoie ici vers n’importe quel cours d’économétrie). On peut le voir sur un dessin,

Yp=reg1$coefficients*base$X1+reg2$coefficients*base$X2 plot(base$Y,predict(reg),ylim=range(Yp),col=rgb(0,0,1,.5),cex=.7,xlab="observé",ylab="prédit") abline(a=0,b=1,lty=2) points(base$Y,Yp,col=rgb(1,0,0,.5),cex=.7) abline(lm(predict(reg)~base$Y),col="blue") abline(lm(Yp~base$Y),col="red") En bleu, on a les prévisions avec le modèle linéaire multiple, et en rouge, en faisant deux régressions indépendantes… Si on regarde sur la droite, le modèle rouge sur-valorise, ou disons sur-estime largement, alors qu’il sous-valorise à gauche. La différence entre les deux droites s’interprète ici comme un biais. Sur le graphique ci-dessous, on peut visualiser la distribution des $\widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}$, en rouge, et les $\widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i}$ , en bleu. Cet excès de dispersion, de variance, qu’on observe sur les points rouges, j’interprète ça comme de la polarisation

plot(density(Yp),col="red")lwd=2) lines(density(predict(reg)),col="blue",lwd=2) En fait, ce que raconte le théorème de Frisch-Waugh (et je renvoie à mon précédant billet pour plus de détails), c’est qu’on a le droit de faire plusieurs régressions, mais en cascade ! et surtout pas indépendamment : je peux expliquer $y$ avec la première variable $x_1$, et ensuite régresser le résidu (ce qu’on n’a pas pu expliquer) sur la seconde $x_2$. Cette méthode donnera la même prévision que le modèle multiple.

On peut aller un peu plus loin: on peut jouer sur la valeur de la corrélation, pour mesure l’écart entre les deux prévisions (ici je prévois pour une observation au hasard, la 78ème)

comp=function(r){ S=matrix(c(1,r,r,1),2,2) set.seed(1) X=rmnorm(n,c(0,0),S) Y = 2+X[,1]-2*X[,2]+E base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2])) reg = lm(Y~0+X1+X2,data=base) reg1 = lm(Y~0+X1,data=base) reg2 = lm(Y~0+X2,data=base) y1=predict(reg) y2=reg1$coefficients*base$X1+reg2$coefficients*base$X2 c(y1,y2,(y2-y1)/y1)} vR=seq(0,.98,by=.02) vc=Vectorize(comp)(vR) plot(vR,vc[3,]*100,ylab="Différence relative (%)",xlab="Corrélation",type="l") On prédit ici pour une observation avec un large $y_i$, ce qui correspondait à la partie de droite du graphique précédant (avec les points rouges et bleus).

C’est ce que j’appelle le deuxième effet kiss-cool. Quand on est dans le premier cas, avec les variables indépendantes (corrélation nulle, c’est à dire à gauche sur ma figure), on explique ce qu’on peut avec la première variable, et on rajoute l’impact de la seconde. Et $\widetilde{y}_i\sim\widehat{y}_i$. Le soucis ici est que je n’ai pas le droit de considérer deux modèles indépendants lors que les variables sont très corrélées. Une partie de l’explication fournie par la seconde variable était déjà inclue dans la première. Par exemple avec deux variables très (positivement) corrélées, la prévision qu’on obtient en ajoutant les deux effets estimés indépendamment avec deux régressions simple $\widetilde{y}_i=\widehat{b}_1x_{1,i}+\widehat{b}_2x_{2,i}$, on sur-estime de 50% à 70% la “vraie” prévision $\widehat{y}_i=\widehat{\beta}_1x_{1,i}+\widehat{\beta}_2x_{2,i}$.

Tous les chercheurs savent savent ça… on parle ici de résultats du tout premier cours de modèles linéaires. Et malgré tout, en pratique, on continue à utiliser cette seconde méthode. Un exemple bien connu est celui de l’évaluation (des étudiants, des chercheurs, peu importe). Par exemple, quand on évalue un dossier de financement pour un chercheur, on nous demande de mettre un score

• pour les publications scientifiques (nombre, qualité, etc)
• pour l’encadrement d’étudiants (nombre, niveau, etc)
• pour la qualité de l’environnement (prestige du labo, etc)
• etc

Et à la fin, on somme tout. Mais on le voit, ces variables sont très très corrélées: si vous êtes dans un labo prestigieux, vous attirez beaucoup de candidatures d’étudiants (et des bons), et avoir beaucoup d’étudiants va permettre d’avoir plus de publications (si on ajoute son nom comme co-auteur). Bref, on est typiquement dans un modèle à double (voire triple) effet kiss-cool. Quelqu’un dans un bon labo aura un bon score sur le troisième item, mais aussi un bon sur le nombre d’étudiants, et aussi sur les publications. Ajouter ces scores est stupide, car on a une spirale infernale (les bons sont sur-évalués, et les moins bon, sous-évalués), c’est ce que racontait mon premier dessin, avec les points rouges et bleus. C’est un effet clivant de polarisation forte.

Si on voulait faire les choses proprement, ce que dit le le théorème de Frisch-Waugh, c’est que les scores devraient être attribués en corrigeant de la corrélation entre les variables

• on peut commencer par calculer un score pour la qualité de l’environnement (prestige du labo, etc)
• à environnement donné, calculer un score pour l’encadrement d’étudiants (nombre, niveau, etc)
• à environnement donné, et à encadrement d’étudiants donné, calculer un score pour les publications
• etc

C’est comme la situation que je voyais en France, où on pouvait avoir une variable qui tenait compte d’un prestige du chercheur (par exemple, être chercheur CNRS donnait un bonus) et une autre sur le dossier de publications. Sauf que les deux sont corrélés. Et la plupart des classements d’universités sont construits à partir de scores qui sont loin d’être indépendants.

Bref, tant que l’évaluation se fera en sommant des scores qui sont construits sur des critères souvent très corrélés, on polarise fortement la population.

Est-ce gênant? A priori oui. Car le message que cela envoie est qu’il existe deux classes, les bons, et les mauvais, alors qu’en réalité, le niveau est beaucoup plus homogène qu’il n’y paraît.  Un petit effet positif se retrouve démultiplié par le fait qu’il va se répercuter (positivement) sur plein d’autres variables. C’est mon effet kiss-cool. Mais on pourrait se dire que c’est un problème de distribution des notes finales. Si l’ordre est préservé, on pourrait se dire que ce n’est pas très grave. Malheureusement, ce n’est pas le cas.

Si on quitte un instant le cas de la corrélation très forte, les rangs des prédictions (c’est à dire les rangs des chercheurs une fois donnés les notes $\widehat{y}_i$ ou $\widetilde{y}_i$ ) sont moins corrélés si la corrélation sous-jacente est importante (mais pas trop)

comp=function(r){ S=matrix(c(1,r,r,1),2,2) set.seed(1) X=rmnorm(n,c(0,0),S) Y = 2+X[,1]-2*X[,2]+E base = data.frame(Y=Y-mean(Y),X1=X[,1]-mean(X[,1]),X2=X[,2]-mean(X[,2])) reg = lm(Y~0+X1+X2,data=base) reg1 = lm(Y~0+X1,data=base) reg2 = lm(Y~0+X2,data=base) y1=predict(reg) y2=reg1$coefficients*base$X1+reg2$coefficients*base$X2 cor(y1,y2,method="spearman")} vR=seq(0,.98,by=.02) vc=Vectorize(comp)(vR) plot(vR,vc,ylab="Corrélation de rangs",xlab="Corrélation",type="l") Autrement dit, si les variables $x_1$ et $x_2$ sont très peu corrélées, on a les mêmes rangs  (globalement). En revanche, si la corrélation entre les variables $x_1$ et $x_2$ augmente, le rang des $\widetilde{y}_i$ est de moins en moins cohérent avec celui entre les $\widehat{y}_i$  (qui devrait être celui que l’on recherche).

Bref, il serait temps de comprendre enfin sérieusement les conséquences de ce joli papier, publié il y a presque 90 ans # Quantile Regression (home made, part 2)

A few months ago, I posted a note with some home made codes for quantile regression… there was something odd on the output, but it was because there was a (small) mathematical problem in my equation. So since I should teach those tomorrow, let me fix them.

## Median

Consider a sample $\{y_1,\cdots,y_n\}$. To compute the median, solve$$\min_\mu \left\lbrace\sum_{i=1}^n|y_i-\mu|\right\rbrace$$which can be solved using linear programming techniques. More precisely, this problem is equivalent to$$\min_{\mu,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^na_i+b_i\right\rbrace$$with $a_i,b_i\geq 0$ and $y_i-\mu=a_i-b_i$, $\forall i=1,\cdots,n$. Heuristically, the idea is to write $y_i=\mu+\varepsilon_i$, and then define $a_i$‘s and $b_i$‘s so that $\varepsilon_i=a_i-b_i$ and $|\varepsilon_i|=a_i+b_i$, i.e. $$a_i=(\varepsilon_i)_+=\max\lbrace0,\varepsilon_i\rbrace=|\varepsilon|\cdot\boldsymbol{1}_{\varepsilon_i>0}$$and$$b_i=(-\varepsilon_i)_+=\max\lbrace0,-\varepsilon_i\rbrace=|\varepsilon|\cdot\boldsymbol{1}_{\varepsilon_i<0}$$denote respectively the positive and the negative parts.

Unfortunately (that was the error in my previous post), the expression of linear programs is$$\min_{\mathbf{z}}\left\lbrace\boldsymbol{c}^\top\mathbf{z}\right\rbrace\text{ s.t. }\boldsymbol{A}\mathbf{z}=\boldsymbol{b},\mathbf{z}\geq\boldsymbol{0}$$In the equation above, with the $a_i$‘s and $b_i$‘s, we’re not far away. Except that we have $\mu\in\mathbb{R}$, while it should be positive. So similarly, set $\mu=\mu^+-\mu^-$ where $\mu^+=(\mu)_+$ and $\mu^-=(-\mu)_+$.

Thus, let$$\mathbf{z}=\big(\mu^+;\mu^-;\boldsymbol{a},\boldsymbol{b}\big)^\top\in\mathbb{R}_+^{2n+2}$$and then write the constraint as $\boldsymbol{A}\mathbf{z}=\boldsymbol{b}$ with $\boldsymbol{b}=\boldsymbol{y}$ and $$\boldsymbol{A}=\big[\boldsymbol{1}_n;-\boldsymbol{1}_n;\mathbb{I}_n;-\mathbb{I}_n\big]$$And for the objective function$$\boldsymbol{c}=\big(\boldsymbol{0},\boldsymbol{1}_n,-\boldsymbol{1}_n\big)^\top\in\mathbb{R}_+^{2n+2}$$

To illustrate, consider a sample from a lognormal distribution,

n = 101 set.seed(1) y = rlnorm(n) median(y)  1.077415

For the optimization problem, use the matrix form, with $3n$ constraints, and $2n+1$ parameters,

So far so good…

## Quantile Regression

Consider the following dataset, with rents of flat, in a major German city, as function of the surface, the year of construction, etc.

base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE)

The linear program for the quantile regression is now$$\min_{\boldsymbol{\beta}^+,\boldsymbol{\beta}^-,\mathbf{a},\mathbf{b}}\left\lbrace\sum_{i=1}^n\tau a_i+(1-\tau)b_i\right\rbrace$$with $a_i,b_i\geq 0$ and $$y_i=\boldsymbol{x}^\top[\boldsymbol{\beta}^+-\boldsymbol{\beta}^-]+a_i-b_i$$$\forall i=1,\cdots,n$ and $\beta_j^+,\beta_j^-\geq 0$ $\forall j=0,\cdots,k$. So use here

require(lpSolve) tau = .3 n=nrow(base) X = cbind( 1, base$area) y = base$rent_euro K = ncol(X) N = nrow(X) A = cbind(X,-X,diag(N),-diag(N)) c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N)) b = base$rent_euro const_type = rep("=",N) r = lp("min",c,A,const_type,b) beta = r$sol[1:K] - r$sol[(1:K+K)] beta  148.946864 3.289674 Of course, we can use R function to fit that model library(quantreg) rq(rent_euro~area, tau=tau, data=base) Coefficients: (Intercept) area 148.946864 3.289674 Here again, it seems to work quite well. We can use a different probability level, of course, and get a plot plot(base$area,base$rent_euro,xlab=expression(paste("surface (",m^2,")")), ylab="rent (euros/month)",col=rgb(0,0,1,.4),cex=.5) sf=0:250 yr=r$solution[2*n+1]+r$solution[2*n+2]*sf lines(sf,yr,lwd=2,col="blue") tau = .9 r = lp("min",c,A,const_type,b) tail(r$solution,2)  121.815505 7.865536 yr=r$solution[2*n+1]+r$solution[2*n+2]*sf lines(sf,yr,lwd=2,col="blue") And we can adapt the later to multiple regressions, of course,

X = cbind(1,base$area,base$yearc) K = ncol(X) N = nrow(X) A = cbind(X,-X,diag(N),-diag(N)) c = c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N)) b = base$rent_euro const_type = rep("=",N) r = lp("min",c,A,const_type,b) beta = r$sol[1:K] - r$sol[(1:K+K)] beta  -5542.503252 3.978135 2.887234 to be compared with library(quantreg) rq(rent_euro~ area + yearc, tau=tau, data=base) Coefficients: (Intercept) area yearc -5542.503252 3.978135 2.887234 Degrees of freedom: 4571 total; 4568 residual # De la pratique de la régression Depuis le début de la session, j’ai imposé une petite innovation, en donnant, environ une semaine sur deux, un petit exercice (obligatoire mais non noté) avant le cours, en vue de forcer à réfléchir (et de donner des éléments de réponse). Par exemple pour le premier cours, il fallait “prévoir” une valeur manquante, et le but était de montrer que, naturellement, on choisit la valeur moyenne. Pour demain, j’avais posé un exercice un peu plus compliqué, sachant qu’on avait vu, lors du dernier cours, comme faire une régression linéaire, et qu’on avait fini en discutant les tests simples (en lien avec la significativité) et les tests multiples. Pour l’exercice, j’avais mis en ligne une petite base de données, download.file("http://freakonometrics.free.fr/data3.RData","data3.RData") load("data3.RData") str(df) 'data.frame': 147 obs. of 10 variables:$ Y : num 11.72 15.91 14.19 11.15 8.31 ... $X1: num 1.33 3.18 0.28 2.08 0.11 1.67 1.97 1.27 4.38 0.52 ...$ X2: num 3.66 3.75 3.32 2.68 4.97 2.98 4.56 1.78 2.83 6.36 ... $X3: num 1.41 3.01 0.34 2.19 0.25 1.69 2.01 1.25 4.41 0.43 ...$ X4: num -3.53 -4.46 -3.35 -7.54 -7.02 -2.53 -6.1 -5.99 -3.92 -5.84 ... $X5: num 0.57 0.01 -0.7 1.62 -0.95 -1.37 1.18 -0.72 2.63 -1.63 ...$ X6: num -0.82 1.2 3.03 -0.91 -1.6 1.77 1 -1.33 1.31 -0.7 ... $X7: num 1.01 0.06 2.02 3.63 2.66 2.53 1.29 3.5 1.17 1.8 ...$ X8: num 8.31 8.52 9.78 7.34 7.26 ... $X9: num 6.04 6.53 7.52 5.61 4.52 6.06 6.2 5.99 6.93 4.38 ... Je vais mettre ici les questions que je posais, et donner des pistes de réflexions, non pas sur les réponses attendues (je n’attends rien de cet exercice à part une réflexion), mais sur la discussion que peut amener chacune des questions, • Faîtes un modèle linéaire pour expliquer Y en fonction des neuf variables explicatives. Combien de variables explicatives garderiez-vous ? Commençons par faire une régression sur toutes les variables summary(lm(Y~., data=df)) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 8.046133 2.193448 3.668 0.000349 *** X1 0.342293 0.915894 0.374 0.709186 X2 -0.040479 0.103409 -0.391 0.696073 X3 1.683875 0.897278 1.877 0.062693 . X4 -0.009254 0.062382 -0.148 0.882295 X5 -1.085367 0.113840 -9.534 &lt; 2e-16 *** X6 0.983207 0.111830 8.792 5.49e-15 *** X7 -0.015646 0.087483 -0.179 0.858327 X8 0.012165 0.094756 0.128 0.898033 X9 0.172210 0.180605 0.954 0.342005 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.019 on 137 degrees of freedom Multiple R-squared: 0.9228, Adjusted R-squared: 0.9178 F-statistic: 182 on 9 and 137 DF, p-value: &lt; 2.2e-16 On a, dans la sortie, une dizaine de tests de Student qui sont évoqués, correspondant au test de l’hypothèse $H_0:\beta_j=0$ (contre l’hypothèse alternative (bilatérale) $H_1:\beta_j\neq 0$), dans un modèle de la forme $$y_i=\beta_0+\beta_1x_{1,i}+\dots+\beta_9x_{9,i}+\varepsilon_i$$C’est ce qu’on appelle le test de significativité de la variable $x_j$ (oui, comme on l’a vu en cours, on peut dire le test parce que les autres tests classiques – Fisher, ou Wald – sont équivalent – sauf qu’au lieu de regarder $t$ on regarde la statistique $t^2$ – ce qui présente l’avantage de voir la statistique de test comme une forme de distance à l’hypothèse $H_0$ : si c’est trop grand, on rejette…). Avec un seuil d’acceptation de l’ordre de 5%, on nous dit que 6 variables ne sont pas significatives. Mais gardons bien en mémoire que le test de significativité de $x_j$ est fait ici en supposant que toutes les autres variables restent dans le modèle. Autre chose: avec un seuil d’acceptation de l’ordre de 10%, une des variables (la troisième) peut être vue comme significative. Faisons un test multiple, pour savoir si on peut supprimer 6 des 9 variables explicatives possibles (faisons le à la main, inutile d’aller chercher un package pour le faire) reg1 = lm(formula = Y ~ ., data = df) reg0 = lm(formula = Y ~ X3+X5+X6, data = df) summary(reg0) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 9.11031 0.16788 54.27 &lt;2e-16 *** X3 1.96784 0.07604 25.88 &lt;2e-16 *** X5 -0.99296 0.08751 -11.35 &lt;2e-16 *** X6 1.08761 0.05934 18.33 &lt;2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.011 on 143 degrees of freedom Multiple R-squared: 0.9206, Adjusted R-squared: 0.919 F-statistic: 553 on 3 and 143 DF, p-value: &lt; 2.2e-16 anova(reg0,reg1) Analysis of Variance Table Model 1: Y ~ X3 + X5 + X6 Model 2: Y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 Res.Df RSS Df Sum of Sq F Pr(&gt;F) 1 143 146.16 2 137 142.12 6 4.0338 0.6481 0.6916 Le test de Fisher nous dit qu’on peut accepter l’hypothèse que les 6 coefficients sont nuls – ici $$H_0:\beta_1=\beta_2=\beta_4=\beta_7=\beta_8=\beta_9=0$$(qui est un test multiple, contrairement au test de Student précédant). Mais il nous dit aussi, qu’individuellement, les trois variables restants semblent significatives. Donc j’aurais tendance à garder 3 variables explicatives (je ne parle pas de la constante : on garde toujours la constante – qui n’explique pas grand chose, sauf la valeur moyenne de $y$). • faites une prévision pour un individu dont on sait que X3=1, X5=1 et X6=8 (en réalité, la vraie question que j’ai posée contenait une typo ce qui la rendait vicieuse parce que ce n’est pas le modèle qu’on vient de construire… mais pour commencer, regardons cette question) On vient de calibrer ce modèle, donc il suffit de faire une prévision, predict(reg0, newdata=data.frame(X3=1, X5=1, X6=8)) 1 18.78603 Mais ce n’était pas la vraie question… • faites une prévision pour un individu dont on sait que X1=3, X5=1 et X6=8 Je pense que pour répondre à cette question, il convient d’oublier tout ce qu’on vient de voir. On nous donne trois informations, et on va voir si on peut les exploiter. Autrement dit, on va commencer par regarder la régression sur ces trois variables, reg156 = lm(formula = Y ~ X1+X5+X6, data = df) summary(reg156) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 9.06688 0.17222 52.65 &lt;2e-16 *** X1 1.98342 0.07800 25.43 &lt;2e-16 *** X5 -1.01438 0.08947 -11.34 &lt;2e-16 *** X6 1.07824 0.06039 17.86 &lt;2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.026 on 143 degrees of freedom Multiple R-squared: 0.9183, Adjusted R-squared: 0.9166 F-statistic: 535.8 on 3 and 143 DF, p-value: &lt; 2.2e-16 Le modèle est ici bon, les trois variables étant significatives. Cela dit, je dis qu’il est “bon” mais on verra vendredi comment discuter davantage ce point… En tous cas, on peut tenter de faire une prévision, et on obtient predict(reg156, newdata=data.frame(X1=3, X5=1, X6=8)) 1 22.62867 • faites une prévision pour un individu dont on sait que X1=3, X3=2, X5=1 et X6=8 Comme auparavant, on nous donne 4 informations… on va regarder le modèle avec les 4 variables reg1356 = lm(formula = Y ~ X1+X3+X5+X6, data = df) summary(reg1356) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 9.10458 0.17131 53.148 &lt;2e-16 *** X1 0.16362 0.89028 0.184 0.854 X3 1.80663 0.88052 2.052 0.042 * X5 -0.99558 0.08895 -11.192 &lt;2e-16 *** X6 1.08649 0.05986 18.152 &lt;2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.014 on 142 degrees of freedom Multiple R-squared: 0.9207, Adjusted R-squared: 0.9184 F-statistic: 411.9 on 4 and 142 DF, p-value: &lt; 2.2e-16 Cette fois, une des variables n’est pas significative (la première) donc on devrait l’enlever: on fait alors la régression juste sur les trois autres variables reg356 = lm(formula = Y ~ X3+X5+X6, data = df) summary(reg356) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 9.11031 0.16788 54.27 &lt;2e-16 *** X3 1.96784 0.07604 25.88 &lt;2e-16 *** X5 -0.99296 0.08751 -11.35 &lt;2e-16 *** X6 1.08761 0.05934 18.33 &lt;2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.011 on 143 degrees of freedom Multiple R-squared: 0.9206, Adjusted R-squared: 0.919 F-statistic: 553 on 3 and 143 DF, p-value: &lt; 2.2e-16 Cette fois, le modèle est “bon” et on peut alors faire la prévision predict(reg356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8)) 1 20.75388 A titre d’information, si on avait fait la prévision sans enlever la variable non significative, on aurait obtenu la valeur suivante predict(reg1356, newdata=data.frame(X1=3, X3=2, X5=1, X6=8)) 1 20.90501 qui est globalement assez proche. On aura l’occasion d’en reparler en cours, quand on se demandera s’il est plus grave d’avoir une valeur non-significative dans la régression, ou d’oublier une variable importante… Mais encore une fois, le but de ces petits exercices est d’appliquer ce qu’on a vu en cours, et d’introduire des questions auxquelles j’apporterai des réponses au prochain cours ! # De l’abus de notation dans les modèles de régression De manière un peu rituelle, je commence toujours mon cours de régression en revenant sur un point important de la statistique : les abus de notation ! Car tout le monde utilise les mêmes lettres (surtout les grecques) pour désigner des objets de nature différente. Dans la majorité des livres, on pourra nous dire sur la même page que $\widehat{\theta}=2.35$ et que $\text{Var}(\widehat{\theta})=1.07$, autrement dit $\widehat{\theta}$ peut désigner en même temps un nombre (dans le premier cas) et une variable aléatoire (dans le second). C’est pour le moins déroutant ! En fait, la raison est assez simple. La statistique commence toujours par un échantillon $\{y_1,y_2,\cdots,y_n\}$, des données, des chiffres. Si on reste là, on fait du descriptif. L’étape classique est ensuite de supposer que les observations $y_i$ sont des réalisations de variables aléatoires $Y_i$, qu’on supposera bien souvent indépendantes et identiquement distribuées. Et $\widehat{\theta}$ sera alors une statistique, c’est à dire une fonction de mes observations. Je peux alors définir $\widehat{\theta}=t(y_1,\cdots,y_n)$ comme étant la statistique observée sur mon échantillon, mais je peux aussi considérer $\widehat{\theta}=t(Y_1,\cdots,Y_n)$, qui est alors une variable aléatoire, mais avec la même notation. Si on voulait aider à comprendre, on utiliserait $\widehat{\Theta}$, mais bon, les choses sont ce qu’elles sont… Et en économétrie, ça devient rapidement un cauchemar quand on commence à parler des résidus… Autre particularité en statistique, c’est que si on distingue l’espérance et la moyenne (empirique), on a un seul mot pour parler de la variance, que ce soit pour une variable aléatoire, ou un vecteur de $\mathbb{R}^n$. On aura ainsi $$\mathbb{E}[Y]=\int y dF(y)$$et$$\overline{y}=\widehat{\mathbb{E}}[\boldsymbol{y}]=\frac{1}{n}\sum_{i=1} y_i$$alors que$$\text{Var}[Y]=\int [y-\mathbb{E}[Y]]^2 dF(y)$$et$$\widehat{\text{Var}}[\boldsymbol{y}]=\frac{1}{n}\sum_{i=1} (y_i-\overline{y})^2$$ Considérons un problème de régression maintenant, avec un modèle de la forme $y_i=\boldsymbol{x}_i^\top\boldsymbol{\beta}+\varepsilon_i$. Ici, $\varepsilon_i$ est un nombre réel, inconnu. Dans une écriture matricielle, on a $\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}$, où cette fois $\boldsymbol{\varepsilon}$ est un vecteur de $\mathbb{R}^n$ (et oui, je suis désolé mais ici $\boldsymbol{X}$ désigne la matrice des covariates, et non pas un vecteur aléatoire… je ferais un billet un jour pour parler du fait que parfois on dit que les $\boldsymbol{x}$ sont donnés et des fois – comme on conditionne suivant $\boldsymbol{X}$, autrement dit, on les voit comme aléatoires). On peut parfois faire une hypothèse quant à la distribution des résidus. Autrement dit, les $\varepsilon_i$ sont vues comme des réalisations de variables aléatoires $\varepsilon_i$, ainsi que $\boldsymbol{\varepsilon}$. On notera ainsi $\boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma})$. Ah oui, autre point juste pour perdre les élèves : $\text{Var}(\boldsymbol{\varepsilon})=\boldsymbol{\Sigma}$ alors que $\text{Var}(\varepsilon_i)=\sigma^2$… Bon, ici comme on suppose les observations indépendentes, et identiquement distribuées, on supposera que $\text{Var}(\boldsymbol{\varepsilon})=\boldsymbol{\Sigma}=\sigma^2\mathbb{I}$. Encore une fois, $\boldsymbol{\varepsilon}$ est (par définition) non observable. Par contre, on peut estimer ces résidus : à partir d’un estimateur $\widehat{\boldsymbol{\beta}}$ de $\boldsymbol{\beta}$, on peut définir $$\widehat{\boldsymbol{\varepsilon}}=\boldsymbol{y}-\widehat{\boldsymbol{y}}=\boldsymbol{y}-\boldsymbol{x}^\top\widehat{\boldsymbol{\beta}}$$Histoire de clarifier, je vais plutôt noter $\widehat{\boldsymbol{e}}$ ces résidus estimés, en utilisant l’estimateur par moindres carrés de $\boldsymbol{\beta}$. On peut noter que $\widehat{\boldsymbol{e}}=(\mathbb{I}-\boldsymbol{H})\boldsymbol{y}$ où classiquement $\boldsymbol{H}=\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top$ est la matrice de projection sur l’espace engendré par toutes les combinaisons linéaires des variables explicatives. Mais là encore, on peut voir le vecteur (numérique) $\widehat{\boldsymbol{e}}$ comme la réalisation d’une variable aléatoire $\widehat{\boldsymbol{E}}$. En particulier, $$\widehat{\boldsymbol{E}}=(\mathbb{I}-\boldsymbol{H})\boldsymbol{Y}=(\mathbb{I}-\boldsymbol{H})\boldsymbol{\varepsilon}$$$\boldsymbol{\varepsilon}$ est notre vecteur aléatoire, centré, de matrice de variance-covariance $\text{Var}(\boldsymbol{\varepsilon})=\sigma^2\mathbb{I}$. On peut alors en déduire que$$\mathbb{E}[\widehat{\boldsymbol{E}}]=(\mathbb{I}-\boldsymbol{H})\mathbb{E}[\boldsymbol{\varepsilon}]=\boldsymbol{0}$$et$$\text{Var}[\widehat{\boldsymbol{E}}]=(\mathbb{I}-\boldsymbol{H})\text{Var}[\boldsymbol{\varepsilon}](\mathbb{I}-\boldsymbol{H})^\top=\sigma^2(\mathbb{I}-\boldsymbol{H})$$(car $(\mathbb{I}-\boldsymbol{H})$ est idenpotent). Cette dernière relation est particulièrement importante, car on notera que $\text{Var}(\widehat{\boldsymbol{E}})\neq\sigma^2\mathbb{I}$. En particulier, si on prend un résidu estimé au hasard $\text{Var}(\widehat{E}_i)=\sigma^2(1-\boldsymbol{H}_{i,i})$ (on avait parlé longuement de $\boldsymbol{H}_{i,i}$ dans un billet récent, on leverage, en particulier on avait vu que $\boldsymbol{H}_{i,i}\in[0,1]$ (on avait discuté la borne inférieur, qui peut être améliorée, en fait $\boldsymbol{H}_{i,i}\in(0,1]$) de telle sorte que $\text{Var}(\widehat{E}_i)\leq\sigma^2$. Si on poursuit un peu, on peut regarder la somme des carrés estimés, et noter que$$\mathbb{E}\big[\sum_{i=1}^n \widehat{E}_i^2\big]=\mathbb{E}[\text{trace}( \widehat{\boldsymbol{E}}\widehat{\boldsymbol{E}}^\top)] =\text{trace}(\mathbb{E}[\text{trace}( \widehat{\boldsymbol{E}}\widehat{\boldsymbol{E}}^\top])$$i.e.$$\mathbb{E}\big[\sum_{i=1}^n \widehat{E}_i^2\big]=\sigma^2\text{trace}(\mathbb{I}-\boldsymbol{H})$$or $\text{trace}(\mathbb{I}-\boldsymbol{H})=n-p$ donc$$\widehat{\sigma}^2=\frac{1}{n-p}\sum_{i=1}^n \widehat{E}_i^2$$est un estimateur sans biais de $\sigma^2$. Et classiquement, on considèrera les résidus Studentisés$$\widehat{R}_i=\frac{\widehat{E}_i}{\widehat{\sigma}\sqrt{1-\boldsymbol{H}_{i,i}}}$$Si je voulais résumer un peu, on pourrait dire que$$\text{Var}(\boldsymbol{E})=\sigma^2\mathbb{I}$$$$\widehat{\text{Var}}(\boldsymbol{E})=\widehat{\sigma}^2\mathbb{I}$$$$\text{Var}(\widehat{\boldsymbol{E}})=\sigma^2(\mathbb{I}-\boldsymbol{H})$$$$\widehat{\text{Var}}(\widehat{\boldsymbol{E}})=\widehat{\sigma}^2(\mathbb{I}-\boldsymbol{H})$$En espérant que ça clarifie un peu…(?) # On leverage Last week, in our STT5100 (applied linear models) class, I’ve introduce the hat matrix, and the notion of leverage. In a classical regression model, $\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}$ (in a matrix form), the ordinary least square estimator of parameter $\boldsymbol{\beta}$ is $$\widehat{\boldsymbol{\beta}}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y}$$The prediction can then be written$$\widehat{\boldsymbol{y}}=\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\underbrace{\color{blue}{\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}}_{\color{blue}{\boldsymbol{H}}}\boldsymbol{y}$$where $\color{blue}{\boldsymbol{H}}$ is called the hat matrix. The matrix is idempotent, i.e. $$\boldsymbol{H}^2={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\textcolor{grey}{\boldsymbol{X}^\top{\boldsymbol{X}}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}}\boldsymbol{X}^\top}={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}=\boldsymbol{H}$$so it can be interpreted as a projection matrix. Furthermore, since$\boldsymbol{H}\boldsymbol{X}=\boldsymbol{X}$ (just do the maths), the projection is on a subspace that contains all the linear combinations of columns of $\boldsymbol{X}$. One can also observe that $\mathbb{I}-\boldsymbol{H}$ is also a projection matrix. And we can write$$\boldsymbol{y}=\underbrace{\boldsymbol{H}\boldsymbol{y}}_{\widehat{\boldsymbol{y}}}+\underbrace{(\mathbb{I}-\boldsymbol{H})\boldsymbol{y}}_{\widehat{\boldsymbol{\varepsilon}}}$$where $\widehat{\boldsymbol{y}}$ is the orthogonal projection of $\boldsymbol{y}$ on the (linear) space of linear combinations of columns of $\boldsymbol{X}$, and $\widehat{\boldsymbol{y}}\perp\widehat{\boldsymbol{\varepsilon}}$, which gives the classical interpretation of residuals, being unpredictible (at least with a linear model using variables $\boldsymbol{X}$). Let’s move a bit faster now (we’ve seen many other properties last week), and consider elements on the diagonal of matrix $\boldsymbol{H}$. Recall that we have so entry $\boldsymbol{H}_{i,i}$ is a measure of the influence of entry $\boldsymbol{y}_i$ on its prediction latex]\widehat{\boldsymbol{y}}_i[/latex]. We have seen that$$\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}(\boldsymbol{H})=\text{trace}(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top)$$which can be written$$\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}\boldsymbol{X}^\top(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1})=\text{trace}(\mathbb{I})=p$$where classically $p=k+1$, where $k$ is the number of explanatory variables. Further, since $\boldsymbol{H}$ is idempotent, we can write (from $\boldsymbol{H}=\boldsymbol{H}^2$) that$$\boldsymbol{H}_{i,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}\boldsymbol{H}_{j,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}^2$$One the one hand, since the second term is positive $\boldsymbol{H}_{i,i}\geq\boldsymbol{H}_{i,i}^2$, i.e. $1\geq\boldsymbol{H}_{i,i}$. And since both terms are positive, then $\boldsymbol{H}_{i,i}\in[0,1]$. And there was a question in the course on the sharpeness of the bounds. Using Anscombe’s dataset, we’ve seen that it was possible to get a leverage of 1. Using something rather similar df = data.frame(x = c(rep(1,10),6), y = c(1:10,8)) plot(df) we obtain model = lm(y~x,data=df) abline(model,col="red",lwd=2) H = lm.influence(model)$hat plot(1:11,H,type="h") The very last observation, the one one the right, is here extremely influencial : if we remove it, the model is completely different ! And here, we reach the upper bound, $\boldsymbol{H}_{11,11}=1$. Observe that all other points are equally influencial, and because on the constraint on the trace of the matrix, $\boldsymbol{H}_{i,i}=1/10$ when $i\in\{1,2,\cdots,10\}$.

Now, what about the lower bound ? In order to have some sort of “non-influencial” observations, consider the two following case.

• the case where one observation (below the first one) is such that $\widehat{\boldsymbol{y}}_{i}=\boldsymbol{y}_{i}$ (perfect prediction)
• the case where one observation (below the tenth one) is such that $\boldsymbol{x}_{i}=\overline{\boldsymbol{x}}$ and $\boldsymbol{y}_{i}=\overline{\boldsymbol{y}}$ (from the first order condition – or normal equation), the fitted regression line always go through point $(\overline{\boldsymbol{x}},\overline{\boldsymbol{y}})$

Let us move two observations from our dataset,

mean(c(4,rep(1,8),6))  1.8 df = data.frame(x = c(4,rep(1,8),6,1.8), y = c(predict(model,newdata=data.frame(x=4)), 2:9,8, predict(model,newdata=data.frame(x=1.8))))

We now have If we compute the leverages, we obtain

model = lm(y~x,data=df) H = lm.influence(model)$hat plot(1:11,H,type="h") so, for the first observation, its leverage actually increased (the blue part), and for the tenth one, we have the lowest influence, but it is not zero. Is it possible to reach zero ? Here, observe that for the tenth observation, $\boldsymbol{H}_{i,i}=1/n$. And actually, that’s the best we can do… We can prove that, in the case of a simple regression (as above)$$\boldsymbol{H}_{i,i}=\frac{1}{n}+\frac{(x_i-\overline{x})^2}{n\text{Var}(x)}$$which is minimum when $x_i=\overline{x}$, and then $\boldsymbol{H}_{i,i}=1/n$, otherwise $\boldsymbol{H}_{i,i}>1/n$. And this property is also valid in a multiple regression (as soon as an intercept is included in the regression – which should always be the case). To prove that result, let $\tilde{\boldsymbol{X}}$ denote the matrix of centered variables $\boldsymbol{X}$, then we can prove that $$\boldsymbol{H}_{i,i}=\frac{1}{n}+\big[\tilde{\boldsymbol{X}}(\tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}})^{-1}\tilde{\boldsymbol{X}}^\top\big]_{i,i}$$(which is basically a matrix version of the previous equation). I can maybe add another comment on Anscombe’s data. We’ve seen that on the right that we did reach 1. But I did not prove it. One way to prove it is actually to focus on the remaining $n-1$ points, on the left. Those have all the same $x$ values. We can prove that if $\boldsymbol{X}_{i_1}=\boldsymbol{X}_{i_2}$, then $$\boldsymbol{H}_{i_1,i_2}=\boldsymbol{X}_{i_1}^\top(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_{i_2}=\boldsymbol{H}_{i_1,i_1}$$hence, using the relationship obtained since the hat matrix is idempotent$$\boldsymbol{H}_{i_1,i_1}=2\boldsymbol{H}_{i_1,i_1}^2+\sum_{j\notin\{i_1,i_2\}}\boldsymbol{H}_{i_1,j}^2$$thus, we now have$$\boldsymbol{H}_{i_1,i_1}\big(1-2\boldsymbol{H}_{i_1,i_1}\big)>0$$i.e. $\boldsymbol{H}_{i_1,i_1}\in[0,1/2]$, where the upper bound becomes $1/(n-1)$ “duplicates”. So for $n-1$ $\boldsymbol{H}_{i,i}$‘s, we have values below $1/(n-1)$, the last one should be below $1$ and the sum has to be $k=2$ . So we have the value of the $n$ $\boldsymbol{H}_{i,i}$‘s. # Données Agrégées et Variables Compositionnelles Avec Enora Belz, nous venons de mettre en ligne une note méthodologique, Données Agrégées et Variables Compositionnelles, sur hal. La réforme du droit sur les données personnelles en Europe rend difficile l’accès aux données individuelles (même souvent non-nominatives), surtout quand on cherche des données jugées sensibles (et souvent, le revenu entre dans cette catégorie). Une solution souvent envisagée est la mise à disposition de données agrégées spatialement. Ces données posent toutefois deux problèmes techniques. Le premier est que les données catégorielles deviennent des compositions. Le second est lié au paradoxe écologique qui dit qu’il est dangereux d’inférer des relations économétriques individuelles à partir de données agrégées. Nous verrons ici comment travailler avec des données compositionnelles (pour éventuellement juste valider une approche classique de régression linéaire-plus simple à interpréter). Et nous évoquerons le second, mais qui reste malheureusement trop général pour pouvoir être traité de manière satisfaisante. # 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 n=100 df=data.frame(matrix(rnorm(n*n),n,n)) names(df)=c("Y",paste("X",1:99,sep="")) 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 reg=function(k){ frm=paste("Y~",paste("X",1:k,collapse="+",sep="")) model=lm(frm,data=df) summary(model)$adj.r.squared}

Let us see what’s going on…

plot(1:99,Vectorize(reg)(1:99)) (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…

reg=function(k){ frm=paste("Y~",paste("X",1:k,collapse="+",sep="")) model=lm(frm,data=df) summary(model)$r.squared} plot(1:99,Vectorize(reg)(1:99)) We can actually do the same with Akaike criteria $AIC_k$ and Schwarz (bayesian) criteria $BIC_k$. reg=function(k){ frm=paste("Y~",paste("X",1:k,collapse="+",sep="")) model=lm(frm,data=df) AIC(model)} plot(1:99,Vectorize(reg)(1:99)) 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 reg=function(k){ frm=paste("Y~",paste("X",1:k,collapse="+",sep="")) model=lm(frm,data=df) BIC(model)} plot(1:99,Vectorize(reg)(1:99)) 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 $i$th 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_i$$i.e. we do note have to estimate (at each round) $n$ models reg=function(k){ frm=paste("Y~",paste("X",1:k,collapse="+",sep="")) model=lm(frm,data=df) h=lm.influence(model)$hat/2 mean( (residuals(model)/1-h)^2 ))} plot(1:99,Vectorize(reg)(1:99)) 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$,

df=data.frame(matrix(rnorm(n*n),n,n)) df=df[,rev(order(abs(cor(df)[1,])))] names(df)=c("Y",paste("X",1:99,sep=""))}

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 2

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

## Geometric Properties of this Linear Model

Let’s define the scalar product in $\mathbb{R}^n$, $⟨\mathbf{a},\mathbf{b}⟩=\mathbf{a}^T\mathbf{b}$, and let’s note $\|\cdot\|$ the associated Euclidean standard, $\|\mathbf{a}\|=\sqrt{\mathbf{a}^T\mathbf{a}}$ (denoted $\|\cdot\|_{\ell_2}$ in the next post). Note $\mathcal{E}_X$ the space generated by all linear combinations of the $\mathbf{X}$ components (adding the constant). If the explanatory variables are linearly independent, $\mathbf{X}$ is a full (column) rank matrix and $\mathcal{E}_X$ is a space of dimension $p+1$. Let’s assume from now on that the variables $\mathbf{x}$  and $y$ are centered here. Note that no law hypothesis is made in this section, the geometric properties are derived from the properties of expectation and variance in the set of finite variance variables.

With this notation, it should be noted that the linear model is written $m(\mathbf{x})=⟨\mathbf{x},\beta⟩$. The space $H_z=\{\mathbf{x}\in\mathbb{R}^{p+1}:m(\mathbf{x})=z\}$ is a hyperplane (affine) that separates the space in two. Let’s define the orthogonal projection operator on $\mathcal{E}_X$, $\Pi_X =\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T$. Thus, the forecast that can be made for it is: $$\widehat{\mathbf{y}}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1} \mathbf{X}^T\mathbf{y}=\Pi_X\mathbf{y}$$. As, $\widehat{\varepsilon}=\mathbf{y}-\widehat{\mathbf{y}}=(\mathbb{I}-\Pi_X)\mathbf{y}=\Pi_{X^\perp}\mathbf{y}$, we note that $\widehat{\varepsilon}\perp\mathbf{x}$, which will be interpreted as meaning that residuals are a term of innovation, unpredictable in the sense that $\Pi_{X }\widehat{\varepsilon}=\mathbf{0}$. The Pythagorean theorem is written here: $$\Vert \mathbf{y} \Vert^2=\Vert \Pi_{ {X}}\mathbf{y} \Vert^2+\Vert \Pi_{ {X}^\perp}\mathbf{y} \Vert^2=\Vert \Pi_{ {X}}\mathbf{y}\Vert^2+\Vert \mathbf{y}-\Pi_{ {X}}\mathbf{y}\Vert^2=\Vert\widehat{\mathbf{y}}\Vert^2+\Vert\widehat{\mathbf{\varepsilon}}\Vert^2$$which is classically translated in terms of the sum of squares: $$\underbrace{\sum_{i=1}^n y_i^2}_{n\times\text{total variance}}=\underbrace{\sum_{i=1}^n \widehat{y}_i^2}_{n\times\text{explained variance}}+\underbrace{\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{n\times\text{residual variance}}$$The coefficient of determination, $R^2$, is then interpreted as the square of the cosine of the angle $\theta$ between $\mathbf{y}$ and $\Pi_X \mathbf{y}$ : $$R^2=\frac{\Vert \Pi_{{X}} \mathbf{y}\Vert^2}{\Vert \mathbf{y}\Vert^2}=1-\frac{\Vert \Pi_{ {X}^\perp} \mathbf{y}\Vert^2}{\Vert \mathbf {y}\Vert^2}=\cos^2(\theta)$$An important application was obtained by Frish & Waugh (1933), when the explanatory variables are divided into two groups, $\mathbf{X}=[\mathbf{X}_1 |\mathbf{X}_2]$, so that the regression becomes $y=\beta_0+\mathbf{X}_1 β_1+\mathbf{X}_2 β_2+\varepsilon$. Frish & Waugh (1933) showed that two successive projections could be considered. Indeed, if $\mathbf{y}_2^\star=\Pi_{X_1^\perp} \mathbf{y}$ and $X_2^\star=\Pi_{X_1^\perp}\mathbf{X}_2$, we can show that $$\widehat{\beta} _2=[{\mathbf{X}_2^\star}^T \mathbf{X}_2^\star]^{-1}{\mathbf{X}_2^\star}^T \mathbf{y}_2^\star$$ In other words, the overall estimate is equivalent to the combination of independent estimates of the two models if $\mathbf{X}_2^\star=\mathbf{X}_2$, i.e. $\mathbf{X}_2\in \mathcal{E}_{X_1}^\perp$, which can be noted $\mathbf{x}_1\perp\mathbf{x}_2$ We obtain here the Frisch-Waugh theorem which guarantees that if the explanatory variables between the two groups are orthogonal, then the overall estimate is equivalent to two independent regressions, on each of the sets of explanatory variables. This is a theorem of double projection, on orthogonal spaces. Many results and interpretations are obtained through geometric interpretations (fundamentally related to the links between conditional expectation and the orthogonal projection in space of variables of finite variance).

This geometric interpretation might help to get a better understanding of the problem of under-identification, i.e. the case where the real model would be $y_i=\beta_0+ \mathbf{x}_1^T \beta_1+\mathbf{x}_2^T \beta_2+\varepsilon_i$, but the estimated model is $y_i=b_0+\mathbf{x}_1^T \mathbf{b}_1+\eta_i$. The maximum likelihood estimator of $\mathbf{b}_1$ is $$\widehat{\mathbf{b}}_1=\mathbf {\beta}_1 + \underbrace{ (\mathbf {X}_1^T\mathbf {X}_1)^{-1} \mathbf {X}_1^T \mathbf {X}_{2} \mathbf{\beta}_2}_{\mathbf{\beta}_{12}}+\underbrace{(\mathbf{X}_1^{T}\mathbf{X}_1)^{-1} \mathbf{X}_1^T\varepsilon}_{\nu}$$so that $\mathbb{E}[\widehat{\mathbf{b}}_1]=\beta_1+\beta_{12}$, the bias ($\beta_{12}$) being null only in the case where $\mathbf{X}_1^T \mathbf{X}_2=\mathbf{0}$ (i. e. $\mathbf{X}_1\perp \mathbf{X}_2$): we find here a consequence of the Frisch-Waugh theorem.

On the other hand, over-identification corresponds to the case where the real model would be $y_i=\beta_0+\mathbf{x}_1^T \beta_1+\varepsilon_i$, but the estimated model is $y_i=b_0+ \mathbf{x}_1^T \mathbf{b} _1+\mathbf{x}_2^T \mathbf{b}_2+\eta_i$. In this case, the estimate is unbiased, in the sense that $\mathbb{E}[\widehat{\mathbf{b}}_1]=\beta_1$ but the estimator is not efficient. Later on, we will discuss an effective method for selecting variables (and avoid over-identification).

## From parametric to non-parametric

We can rewrite equation (4) in the form $\widehat{\mathbf{y}}=\Pi_X\mathbf{y}$ which helps us to see the forecast directly as a linear transformation of the observations. More generally, a linear predictor can be obtained by considering $m(\mathbf{x})=\mathbf{s}_{\mathbf{x}}^T \mathbf{y}$, where $\mathbf{s}_{\mathbf{x}}$ is a weight vector, which depends on $\mathbf{x}$, interpreted as a smoothing vector. Using the vectors $\mathbf{s}_{\mathbf{x}_i}$, calculated from the observations $\mathbf{x}_i$, we obtain a matrix $\mathbf{S}$ of size $n\times n$, and $\widehat{\mathbf{y}}=\mathbf{S}\mathbf{y}$. In the case of the linear regression described above, $\mathbf{s}_{\mathbf{x}}=\mathbf{X}[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{x}$, and in that case $\text{trace}(\mathbf{S})$ is the number of columns in the $\mathbf{X}$ matrix (the number of explanatory variables). In this context of more general linear predictors, $\text{trace}(\mathbf{S})$ is often seen as equivalent to the number of parameters (or complexity, or dimension, of the model), and $\nu=n-\text{trace}(\mathbf{S})$ is then the number of degrees of freedom (see Ruppert et al., 2003; Simonoff, 1996). The principle of parsimony says that we should minimize this dimension (the trace of the matrix $\mathbf{S}$) as much as possible. But in the general case, this dimension is more to obtain, explicitely.

The estimator introduced by Nadaraya (1964) and Watson (1964), in the case of a simple non-parametric regression, is also written in this form since$$\widehat{m}_h(x)=\mathbf{s}_{x}^T\mathbf{y}=\sum_{i=1}^n \mathbf{s}_{x,i}y_i$$where$$\mathbf{s}_{x,i}=\frac{K_h(x-x_i)}{K_h(x-x_1)+\cdots+K_h(x-x_n)}$$ where $K(\cdot)$ is a kernel function, which assigns a value that is lower the closer $x_i$ is to $x$, and $h>0$ is the bandwidth. The introduction of this metaparameter $h$ is an important issue, as it should be chosen wisely. Using asymptotic developments, we can show that if $X$ has density $f$, $$\text{biais}[\widehat{m}_h(x)]=\mathbb{E}[\widehat{m}_h(x)]-m(x)\sim {h^2}\left(\frac{C_1 }{2}m''(x)+C_2 m'(x)\frac{f'(x)}{f(x)}\right)$$and $$\displaystyle{{\text{Var}[\widehat{m}_h(x)]\sim\frac{C_3}{{nh}}\frac{\sigma(x)}{f(x)}}}$$for some constants that can be estimated (see Simonoff (1996) for a discussion). These two functions evolve inversely with $h$, as shown in Figure 1 (where the metaparameter on the $x$-axis is here, actually, $h^{-1}$). Keep in ming that we will see a similar graph in the context of machine learning models. Figure 1. Choice of meta-parameter and the Goldilocks problem: it must not be too large (otherwise there is too much variance), nor too small (otherwise there is too much bias).

The natural idea is then to try to minimize the mean square error, the MSE, defined as $bias[\widehat{m}_h (x)]^2+Var[\widehat{m}_h (x)]$, and them integrate over $x$, which gives an optimal value for $h$ of the form $h^\star=O(n^{-1/5})$, and reminds us of Silverman’s rule – see Silverman (1986). In larger dimensions, for continuous $\mathbf{x}$ variables, a multivariate kernel with matrix bandwidth $\mathbf{H}$ can be used, and$$\mathbb{E}[\widehat{m}_{\mathbf{H}}(\mathbf{x})]\sim m(\mathbf{x})+\frac{C_1}{2}\text{trace}\big(\mathbf{H}^Tm''(\mathbf{x})\mathbf{H}\big)+C_2\frac{m'(\boldsymbol{x})^T\mathbf{H}\mathbf{H}^T \nabla f(\mathbf{x})}{f(\mathbf{x})}$$while$$\text{Var}[\widehat{m}_{\mathbf{H}}(\mathbf{x})]\sim\frac{C_3}{n~\text{det}(\mathbf{H})}\frac{\sigma(\mathbf{x})}{f(\mathbf{x})}$$
If $\mathbf{H}$ is a diagonal matrix, with the same term $h$  on the diagonal, then $h^\star=O(n^{-1/(4+dim(\mathbf{x}))}$. However, in practice, there will be more interest in the integrated version of the quadratic error, $$MISE(\widehat{m}_{h})=\mathbb{E}[MSE(\widehat{m}_{h}(X))]=\int MSE(\widehat{m}_{h}(x))dF(x)$$and we can prove that $$MISE[\widehat{m}_h]\sim \overbrace{\frac{h^4}{4}\left(\int x^2k(x)dx\right)^2\int\big[m''(x)+2m'(x)\frac{f'(x)}{f(x)}\big]^2dx}^{\text{bias}^2} +\overbrace{\frac{\sigma^2}{nh}\int k^2(x)dx \cdot\int\frac{dx}{f(x)}}^{\text{variance}}$$as n→∞ and nh→∞. Here we find an asymptotic relationship that again recalls Silverman’s (1986) order of magnitude, $$h^\star =n^{-\frac{1}{5}}\left(\frac{C_1\int \frac{dx}{f(x)}}{C_2\int \big[m''(x)+2m'(x)\frac{f'(x)}{f(x)}\big]dx}\right)^{\frac{1}{5}}$$The main problem here, in practice, is that many of the terms in the expression above are unknown. Automatic learning offers computational techniques, when the econometrician used to searching for asymptotic (mathematical) properties.

To be continued (references mentioned above are online here)…

# Régression sur une variable qualitative et ANOVA

Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.

Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") Davis[12,c(2,3)]=Davis[12,c(3,2)] Davis=data.frame(Y=Davis$weight * 2.204622, X1=Davis$sex)

On voulait estimer le modèle $$y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_i$$On avait vu que l’on pouvait passer par l’écriture matricielle

 X=cbind(Davis$X1=='F',Davis$X1=='M') Y=Davis$Y car la matrice $\mathbf{X}^T\mathbf{X}$ est inversible (une fois que l’on enlève la constante)  solve(t(X)%*%X) [,1] [,2] [1,] 0.008928571 0.00000000 [2,] 0.000000000 0.01136364 et donc l’estimateur par moindres carrés est (classiquement)$$\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}$$  solve(t(X)%*%X) %*% (t(X)%*%Y) [,1] [1,] 125.4272 [2,] 167.3258 ce qui correspond effectivement à la sortie de R,  reg=lm(Y~0+X1,data=Davis) summary(reg) Coefficients: Estimate Std. Error t value Pr(&gt;|t|) X1F 125.427 1.960 64.00 &lt;2e-16 *** X1M 167.326 2.211 75.68 &lt;2e-16 *** Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes x=Y[X[,1]==1] y=Y[X[,2]==1] nx=length(x) ny=length(y) On avait vu en cours que les $\widehat{\mathbf{\beta}}$ avaient une interprétation très simple, puisque$$\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_i$$autrement dit $\widehat{{\beta}}_M$ est le poids moyen des hommes. Et en effet  mean(y)  167.3258 C’est finalement très naturel, ou intuitif. On peut maintenant s’interroger sur l’écart-type de l’estimateur de $\widehat{{\beta}}_M$. Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici  sqrt(var(y)/ny)  2.794391 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)  2.794391 car pour rappel$$\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}$$Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de $\mathbf{\beta}$ est proportionnel à $\sigma^2$, la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de $\sigma^2$  s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2)) sqrt(s2/ny)  2.210863 et en effet, on retombe sur la valeur donnée dans le tableau de régresion  sqrt(s2/nx)  1.959721 (pareil pour l’autre coefficient). On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).  X=cbind(1,Davis$X1=='M')

Là encore, le modèle devient identifiable, et on obtient ici

 solve(t(X)%*%X) %*% (t(X)%*%Y) [,1] [1,] 125.42724 [2,] 41.89855

On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence

mean(y)-mean(x)  41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis) summary(reg2)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 125.427 1.960 64.00 &lt;2e-16 *** X1M 41.899 2.954 14.18 &lt;2e-16 ***

Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)

 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)   Two Sample t-test   data: Y[X[, 1] == 1] and Y[X[, 2] == 1] t = -6.4475, df = 286, p-value = 4.826e-10 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -30.62603 -16.30035 sample estimates: mean of x mean of y 143.8626 167.3258

Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable $\mathbf{y}$), ça marche ! Posons

 Davis$Y=sample(Davis$Y)

ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici

 reg2=lm(Y~X1,data=Davis) summary(reg2)   Call: lm(formula = Y ~ X1, data = Davis)   Residuals: Min 1Q Median 3Q Max -57.458 -22.184 -5.512 17.809 118.912   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 143.4382 2.7820 51.56 &lt;2e-16 *** X1M 0.9645 4.1940 0.23 0.818 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 29.44 on 198 degrees of freedom Multiple R-squared: 0.000267, Adjusted R-squared: -0.004782 F-statistic: 0.05289 on 1 and 198 DF, p-value: 0.8183

autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient

 Y=Davis$Y t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE) Two Sample t-test data: Y[X[, 1] == 1] and Y[X[, 2] == 1] t = -0.22998, df = 198, p-value = 0.8183 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -9.235209 7.306165 sample estimates: mean of x mean of y 143.4382 144.4027 et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents. # Convex Regression Model This morning during the lecture on nonlinear regression, I mentioned (very) briefly the case of convex regression. Since I forgot to mention the codes in R, I will publish them here. Assume that $y_i=m(\mathbf{x}_i)+\varepsilon_i$ where $m:\mathbb{R}^d\rightarrow \mathbb{R}$ is some convex function. Then $m$ is convex if and only if $\forall\mathbf{x}_1,\mathbf{x}_2\in\mathbb{R}^d$, $\forall t\in[0,1]$, $$m(t\mathbf{x}_1+[1-t]\mathbf{x}_2) \leq tm(\mathbf{x}_1)+[1-t]m(\mathbf{x}_2)$$Hidreth (1954) proved that if$$m^\star=\underset{m \text{ convex}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-m(\mathbf{x_i})\big)^2\right\rbrace$$then $\mathbf{\theta}^\star=(m^\star(\mathbf{x_1}),\cdots,m^\star(\mathbf{x_n}))$ is unique. Let $\mathbf{y}=\mathbf{\theta}+\mathbf{\varepsilon}$, then $$\mathbf{\theta}^\star=\underset{\mathbf{\theta}\in \mathcal{K}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \big(y_i-\theta_i)\big)^2\right\rbrace$$where$$\mathcal{K}=\{\mathbf{\theta}\in\mathbb{R}^n:\exists m\text{ convex },m(\mathbf{x}_i)=\theta_i\}$$. I.e. $\mathbf{\theta}^\star$ is the projection of $\mathbf{y}$ onto the (closed) convex cone $\mathcal{K}$. The projection theorem gives existence and unicity. For convenience, in the application, we will consider the real-valued case, $m:\mathbb{R}\rightarrow \mathbb{R}$, i.e. $y_i=m(x_i)+\varepsilon_i$. Assume that observations are ordered $x_1\leq x_2\leq\cdots \leq x_n$. Here $$\mathcal{K}=\left\lbrace\mathbf{\theta}\in\mathbb{R}^n:\frac{\theta_2-\theta_1}{x_2-x_1}\leq \frac{\theta_3-\theta_2}{x_3-x_2}\leq \cdots \leq \frac{\theta_n-\theta_{n-1}}{x_n-x_{n-1}}\right\rbrace$$ Hence, quadratic program with $n-2$ linear constraints. $m^\star$ is a piecewise linear function (interpolation of consecutive pairs $(x_i,\theta_i^\star)$). If $m$ is differentiable, $m$ is convex if $$m(\mathbf{x})+ \nabla m(\mathbf{x})^{\text{T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})$$ More generally, if $m$ is convex, then there exists $\xi_{\mathbf{x}}\in\mathbb{R}^n$ such that $$m(\mathbf{x})+ \xi_{\mathbf{x}}^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y})$$ $\xi_{\mathbf{x}}$ is a subgradient of $m$ at ${\mathbf{x}}$. And then $$\partial m(\mathbf{x})=\big\lbrace m(\mathbf{x})+ \xi^{\text{ T}}\cdot[\mathbf{y}-\mathbf{x}] \leq m(\mathbf{y}),\forall \mathbf{y}\in\mathbb{R}^n\big\rbrace$$ Hence, $\mathbf{\theta}^\star$ is solution of $$\text{argmin}\big\lbrace\|\mathbf{y}-\mathbf{\theta}\|^2\big\rbrace$$$$\text{subject to }\theta_i+\xi_i^{\text{ T}}[\mathbf{x}_j-\mathbf{x}_i]\leq\mathbf{\theta}_j,~\forall i,j$$ and $\xi_1,\cdots,\xi_n\in\mathbb{R}^n$. Now, to do it for real, use cobs package for constrained (b)splines regression, library(cobs) To get a convex regression, use plot(cars) x = cars$speed y = cars$dist rc = conreg(x,y,convex=TRUE) lines(rc, col = 2) Here we can get the values of the knots rc Call: conreg(x = x, y = y, convex = TRUE) Convex regression: From 19 separated x-values, using 5 inner knots, 7, 8, 9, 20, 23. RSS = 1356; R^2 = 0.8766; needed (5,0) iterations and actually, if we use them in a linear-spline regression, we get the same output here reg = lm(dist~bs(speed,degree=1,knots=c(4,7,8,9,,20,23,25)),data=cars) u = seq(4,25,by=.1) v = predict(reg,newdata=data.frame(speed=u)) lines(u,v,col="green") Let us add vertical lines for the knots abline(v=c(4,7,8,9,20,23,25),col="grey",lty=2) # Parallelizing Linear Regression or Using Multiple Sources My previous post was explaining how mathematically it was possible to parallelize computation to estimate the parameters of a linear regression. More speficially, we have a matrix $\mathbf{X}$ which is $n\times k$ matrix and $\mathbf{y}$ a $n$-dimensional vector, and we want to compute $\widehat{\mathbf{\beta}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}$ by spliting the job. Instead of using the $n$ observations, we’ve seen that it was to possible to compute “something” using the first $n_1$ rows, then the next $n_2$ rows, etc. Then, finally, we “aggregate” the $m$ objects created to get our overall estimate. ## Parallelizing on multiple cores Let us see how it works from a computational point of view, to run each computation on a different core of the machine. Each core will see a slave, computing what we’ve seen in the previous post. Here, the data we use are y = cars$dist X = data.frame(1,cars$speed) k = ncol(X) On my laptop, I have three cores, so we will split it in $m=3$ chunks library(parallel) library(pbapply) ncl = detectCores()-1 cl = makeCluster(ncl) This is more or less what we will do: we have our dataset, and we split the jobs, We can then create lists containing elements that will be sent to each core, as Ewen suggested, chunk = function(x,n) split(x, cut(seq_along(x), n, labels = FALSE)) a_parcourir = chunk(seq_len(nrow(X)), ncl) for(i in 1:length(a_parcourir)) a_parcourir[[i]] = rep(i, length(a_parcourir[[i]])) Xlist = split(X, unlist(a_parcourir)) ylist = split(y, unlist(a_parcourir)) It is also possible to simplify the QR functions we will use compute_qr = function(x){ list(Q=qr.Q(qr(as.matrix(x))),R=qr.R(qr(as.matrix(x)))) } get_Vlist = function(j){ Q3 = QR1[[j]]$Q %*% Q2list[[j]] t(Q3) %*% ylist[[j]] } clusterExport(cl, c("compute_qr", "get_Vlist"), envir=environment())

Then, we can run our functions on each core. The first one is

 QR1 = parLapply(cl=cl,Xlist, compute_qr)

note that it is also possible to use

 QR1 = pblapply(Xlist, compute_qr, cl=cl)

which will include a progress bar (that can be nice when the database is rather large). Then use

 R1 = pblapply(QR1, function(x) x$R, cl=cl) %&gt;% do.call("rbind", .) Q1 = qr.Q(qr(as.matrix(R1))) R2 = qr.R(qr(as.matrix(R1))) Q2list = split.data.frame(Q1, rep(1:ncl, each=k)) clusterExport(cl, c("QR1", "Q2list", "ylist"), envir=environment()) Vlist = pblapply(1:length(QR1), get_Vlist, cl=cl) sumV = Reduce('+', Vlist) and finally the ouput is solve(R2) %*% sumV [,1] X1 -17.579095 X2 3.932409 which is what we were expecting… ## Using multiple sources In practice, it might also happen that various “servers” have the data, but we cannot get a copy. But it is possible to run some functions on their server, and get some output, that we can use afterwards. Datasets are supposed to be available somewhere. We can send a request, and get a matrix. Then we we aggregate all of them, and send another request. That’s what we will do here. Provider $j$ should run $f_1(\mathbf{X})$ on his part of the data, that function will return $R^{(1)}_j$. More precisely, to the first provider, send function1 = function(subX){ return(qr.R(qr(as.matrix(subX))))} R1 = function1(Xlist[]) and actually, send that function to all providers, and aggregate the output for(j in 2:m) R1 = rbind(R1,function1(Xlist[[j]])) The create on your side the following objects Q1 = qr.Q(qr(as.matrix(R1))) R2 = qr.R(qr(as.matrix(R1))) Q2list=list() for(j in 1:m) Q2list[[j]] = Q1[(j-1)*k+1:k,] Finally, contact one last time the providers, and send one of your objects function2=function(subX,suby,Q){ Q1=qr.Q(qr(as.matrix(subX))) Q2=Q return(t(Q1%*%Q2) %*% suby)} Provider $j$ should then run $f_2(\mathbf{X},\mathbf{y},Q_j^{(2)})$ on his part of the data, using also $Q_j^{(2)}$ as argument (that we obtained on own side) and that function will return $(\mathbf{Q}^{(2)}_j\mathbf{Q}^{(1)}_j)^{T}_j\mathbf{y}_j$. For instance, ask the first provider to run sumV = function2(Xlist[],ylist[], Q2list[]) and do the same with all providers for(j in 2:m) sumV = sumV+ function2(Xlist[[j]],ylist[[j]], Q2list[[j]]) solve(R2) %*% sumV [,1] X1 -17.579095 X2 3.932409 which is what we were expecting… # Linear Regression, with Map-Reduce Sometimes, with big data, matrices are too big to handle, and it is possible to use tricks to numerically still do the map. Map-Reduce is one of those. With several cores, it is possible to split the problem, to map on each machine, and then to agregate it back at the end. Consider the case of the linear regression, $\mathbf{y}=\mathbf{X}\mathbf{\beta}+\mathbf{\varepsilon}$ (with classical matrix notations). The OLS estimate of $\mathbf{\beta}$ is $\widehat{\mathbf{\beta}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}$. To illustrate, consider a not too big dataset, and run some regression. lm(dist~speed,data=cars)$coefficients (Intercept) speed -17.579095 3.932409 y=cars$dist X=cbind(1,cars$speed) solve(crossprod(X,X))%*%crossprod(X,y) [,1] [1,] -17.579095 [2,] 3.932409

How is this computed in R? Actually, it is based on the QR decomposition of $\mathbf{X}$, $\mathbf{X}=\mathbf{Q}\mathbf{R}$, where $\mathbf{Q}$ is an orthogonal matrix (ie $\mathbf{Q}^T\mathbf{Q}=\mathbb{I}$). Then $\widehat{\mathbf{\beta}}=[\mathbf{X}^T\mathbf{X}]^{-1}\mathbf{X}^T\mathbf{y}=\mathbf{R}^{-1}\mathbf{Q}^T\mathbf{y}$

solve(qr.R(qr(as.matrix(X)))) %*% t(qr.Q(qr(as.matrix(X)))) %*% y [,1] [1,] -17.579095 [2,] 3.932409

So far, so good, we get the same output. Now, what if we want to parallelise computations. Actually, it is possible.

Consider $m$ blocks

m = 5

and split vectors and matrices
$$\mathbf{y}=\left[\begin{matrix}\mathbf{y}_1\\\mathbf{y}_2\\\vdots \\\mathbf{y}_m\end{matrix}\right]$$ and $$\mathbf{X}=\left[\begin{matrix}\mathbf{X}_1\\\mathbf{X}_2\\\vdots\\\mathbf{X}_m\end{matrix}\right]=\left[\begin{matrix}\mathbf{Q}_1^{(1)}\mathbf{R}_1^{(1)}\\\mathbf{Q}_2^{(1)}\mathbf{R}_2^{(1)}\\\vdots \\\mathbf{Q}_m^{(1)}\mathbf{R}_m^{(1)}\end{matrix}\right]$$
To split vectors and matrices, use (eg)

Xlist = list() for(j in 1:m) Xlist[[j]] = X[(j-1)*10+1:10,] ylist = list() for(j in 1:m) ylist[[j]] = y[(j-1)*10+1:10]

and get small QR recomposition (per subset)

QR1 = list() for(j in 1:m) QR1[[j]] = list(Q=qr.Q(qr(as.matrix(Xlist[[j]]))),R=qr.R(qr(as.matrix(Xlist[[j]]))))

Consider the QR decomposition of $\mathbf{R}^{(1)}$ which is the first step of the reduce part$$\mathbf{R}^{(1)}=\left[\begin{matrix}\mathbf{R}_1^{(1)}\\\mathbf{R}_2^{(1)}\\\vdots \\\mathbf{R}_m^{(1)}\end{matrix}\right]=\mathbf{Q}^{(2)}\mathbf{R}^{(2)}$$where$$\mathbf{Q}^{(2)}=\left[\begin{matrix}\mathbf{Q}^{(2)}_1\\\mathbf{Q}^{(2)}_2\\\vdots\\\mathbf{Q}^{(2)}_m\end{matrix}\right]$$

R1 = QR1[]$R for(j in 2:m) R1 = rbind(R1,QR1[[j]]$R) Q1 = qr.Q(qr(as.matrix(R1))) R2 = qr.R(qr(as.matrix(R1))) Q2list=list() for(j in 1:m) Q2list[[j]] = Q1[(j-1)*2+1:2,]

Define – as step 2 of the reduce part$$\mathbf{Q}^{(3)}_j=\mathbf{Q}^{(2)}_j\mathbf{Q}^{(1)}_j$$
and$$\mathbf{V}_j=\mathbf{Q}^{(3)T}_j\mathbf{y}_j$$

Q3list = list() for(j in 1:m) Q3list[[j]] = QR1[[j]]\$Q %*% Q2list[[j]] Vlist = list() for(j in 1:m) Vlist[[j]] = t(Q3list[[j]]) %*% ylist[[j]]

and finally set – as the step 3 of the reduce part$$\widehat{\mathbf{\beta}}=[\mathbf{R}^{(2)}]^{-1}\sum_{j=1}^m\mathbf{V}_j$$

sumV = Vlist[] for(j in 2:m) sumV = sumV+Vlist[[j]] solve(R2) %*% sumV [,1] [1,] -17.579095 [2,] 3.932409

It looks like we’ve been able to parallelise our linear regression…