# Some general thoughts on Partial Dependence Plots with correlated covariates

The partial dependence plot is a nice tool to analyse the impact of some explanatory variables when using nonlinear models, such as a random forest, or some gradient boosting.The idea (in dimension 2), given a model $m(x_1,x_2)$ for $\mathbb{E}[Y|X_1=x_1,X_2=x_2]$. The partial dependence plot for variable $x_1$ is model $m$ is function $p_1$ defined as $x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)]$. This can be approximated, using some dataset using $$\widehat{p}_1(x_1)=\frac{1}{n}\sum_{i=1}^n m(x_1,x_{2,i})$$My concern here what the interpretation of that plot when there are some (strongly) correlated covariates. Let us generate some dataset to start with

n=1000 library(mnormt) r=.7 set.seed(1234) X = rmnorm(n,mean = c(0,0),varcov = matrix(c(1,r,r,1),2,2)) Y = 1+X[,1]-2*X[,2]+rnorm(n)/2 df = data.frame(Y=Y,X1=X[,1],X2=X[,2])

As we can see, the true model is here is $y_i=\beta_0+\beta_1 x_{1,i}+\beta_2x_{2,i}+\varepsilon_i$ where $\beta_1 =1$ but the two variables are positively correlated, and the second one has a strong negative impact. Note that here

reg = lm(Y~.,data=df) summary(reg)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 1.01414 0.01601 63.35 &lt;2e-16 *** X1 1.02268 0.02305 44.37 &lt;2e-16 *** X2 -2.03248 0.02342 -86.80 &lt;2e-16 ***

If we estimate a wrongly specified model $y_i=b_0+b_1 x_{1,i}+\eta_i$, we would get

reg1 = lm(Y~X1,data=df) summary(reg1)   Coefficients: Estimate Std. Error t value Pr(&gt;|t|) (Intercept) 1.03522 0.04680 22.121 &lt;2e-16 *** X1 -0.44148 0.04591 -9.616 &lt;2e-16 ***

Thus, on the proper model, $\widehat{\beta}_1\sim+1.02$ while $\widehat{b}_1\sim-0.44$  on the mispecified model.

Now, let us look at the parial dependence plot of the good model, using standard R dedicated packages,

library(pdp) pdp::partial(reg, pred.var = "X1", plot = TRUE, plot.engine = "ggplot2") which is the linear line $y=1+x$, that corresponds to $y=\beta_0+\beta_1x$.

library(DALEX) plot(DALEX::single_variable(DALEX::explain(reg, data=df),variable = "X1",type = "pdp")) which corresponds to the previous graph. Here, it is also possible to creaste our own function to compute that partial dependence plot,

pdp1 = function(x1){ nd = data.frame(X1=x1,X2=df$X2) mean(predict(reg,newdata=nd)) } that will be the straight line below (the dotted line is the theoretical one $y=1+x$, vx=seq(-3.5,3.5,length=101) vpdp1 = Vectorize(pdp1)(vx) plot(vx,vpdp1,type="l") abline(a=1,b=1,lty=2) which is very different from the univariate regression on $x_1$ abline(reg1,col="red") Actually, the later is very consistent with a local regression, only on $x_1$ library(locfit) lines(locfit(Y~X1,data=df),col="blue") Now, to get back to the definition of the partial dependence plot, $x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2}}[m(x_1,X_2)]$, in the context of correlated variable, I was wondering if it would not make more sense to consider some local version actually, something like $x_1\mapsto\mathbb{E}_{\mathbb{P}_{X_2|X_1}}[m(x_1,X_2)]$. My intuition was that, somehow, it did not make any sense to consider any $X_2$ while $X_1$ was fixed (and equal to $x_1$). But it would make more sense actually to look at more valid $X_2$‘s given the value of $X_1$. And a natural estimate could be some $k$ neareast-neighbors, i.e. $$\tilde{p}_1(x_1)=\frac{1}{k}\sum_{i\in\mathcal{V}_k(x)}^n m(x_1,x_{2,i})$$ where $\mathcal{V}_k(x)$ is the set of indices of the $k$ $x_i$‘s that are the closest to $x$, i.e. lpdp1 = function(x1){ nd = data.frame(X1=x1,X2=df$X2) idx = rank(abs(df$X1-x1)) mean(predict(reg,newdata=nd[idx&lt;50,])) } vlpdp1 = Vectorize(lpdp1)(vx) lines(vx,vlpdp1,col="darkgreen",lwd=2) Surprisingly (?), this local partial dependence plot gives a curve that corresponds to the simple regression… # 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.

# 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. # How long could it take to run a regression This afternoon, while I was discussing with Montserrat (aka @mguillen_estany) we were wondering how long it might take to run a regression model. More specifically, how long it might take if we use a Bayesian approach. My guess was that the time should probably be linear in , the number of observations. But I thought I would be good to check. Let us generate a big dataset, with one million rows, > n=1e6 > X=runif(n) > Y=2+5*X+rnorm(n) > B=data.frame(X,Y) Consider as a benchmark the standard linear regression, > lm_freq = function(n){ + idx = sample(1:1e6,size=n) + reg = lm(Y~X,data=B[idx,]) + summary(reg) + } Here the regression is a subset of smaller size. We can do the same with a Bayesian approach, using stan, > stan_lm =" + data { + int N; + vector[N] x; + vector[N] y; + } + parameters { + real alpha; + real beta; + real tau; + } + transformed parameters { + real sigma; + sigma <- 1 / sqrt(tau); + } + model{ + y ~ normal(alpha + beta * x, sigma); + alpha ~ normal(0, 10); + beta ~ normal(0, 10); + tau ~ gamma(0.001, 0.001); + } + " Define then the model > library(rstan) > system.time( stanmodel <<- stan_model(model_code = stan_lm)) utilisateur système écoulé 0.043 0.000 0.043 We want to see how long it might take to run a regression, > lm_bayes = function(n){ + idx = sample(1:1e6,size=n) + fit = sampling(stanmodel, + data = list(N=n, + x=X[idx], + y=Y[idx]), + iter = 1000, warmup=200) + summary(fit) + } We use the following package to see how long it takes > library(microbenchmark) > time_lm = function(n){ + M = microbenchmark(lm_freq(n), + lm_bayes(n),times=50) + return(apply( matrix(M$time,nrow=2),1,mean))
+ }

We can now compare the time it took with ten, one hundred, on thousand, and ten thousand observations,

> vN = c(10,100,1000,10000)
> T = Vectorize(time_lm)(vN)

we can then plot it

> plot(vN,T[2,]/1e6,log="xy",col="red",type="b",
+      xlab="Number of Observations",ylab="Time")
> lines(vN,T[1,]/1e6,col="blue",type="b") It looks like (if we forget about the very small sample) that the time it takes to run a regression is linear, with the two techniques (the frequentist and the Bayesian ones).

And actually, the same story olds for logistic regressions. Consider the following dataset

> n=1e6
> X=runif(n)
> S=-3+2*X+rnorm(n)
> Y=rbinom(n,size=1,prob=exp(S)/(1+exp(S)))
> B=data.frame(X,Y)

The frequentist version of the logistic regression is

> glm_freq = function(n){
+   idx = sample(1:1e6,size=n)
+   reg = glm(Y~X,data=B[idx,],family=binomial)
+   summary(reg)
+ }

and the Bayesian one, using stan,

> stan_glm = "
+ data {
+ int N;
+ vector[N] x;
+ int<lower=0,upper=1> y[N];
+ }
+ parameters {
+ real alpha;
+ real beta;
+ }
+ model {
+ alpha ~ normal(0, 10);
+ beta ~ normal(0, 10);
+ y ~ bernoulli_logit(alpha + beta * x);
+ }
+ "
> stanmodel = stan_model(model_code = stan_glm) )
> glm_bayes = function(n){
+   idx = sample(1:1e6,size=n)
+   fit = sampling(stanmodel,
+        data = list(N=n,
+        x = X[idx],
+        y = Y[idx]),
+        iter = 1000, warmup=200)
+   summary(fit)
+ }

Again, we can see how long it takes to run those regression models

> time_gl m= function(n){
+   M = microbenchmark(glm_freq(n),
+   glm_bayes(n),times=50)
+   return(apply( matrix(M$time,nrow=2),1,mean)) + } # Regression Models, It’s Not Only About Interpretation Yesterday, I did upload a post where I tried to show that “standard” regression models where not performing bad. At least if you include splines (multivariate splines) to take into accound joint effects, and nonlinearities. So far, I do not discuss the possible high number of features (but with boostrap procedures, it is possible to assess something related to variable importance, that people from machine learning like). But my post was not complete: I was simply plotting the prediction obtained by some model. And it “looked like” the regression was nice, but so were the random forrest, the -nearest neighbour and boosting algorithm. What if we compare those models on new data? # On Some Alternatives to Regression Models When you start discussing with people in machine learning, you quickly hear something like “forget your econometric models, your GLMs, I can easily find a machine learning ‘model’ that can beat yours”. I am usually very sceptical, especially when I hear “easily” or “always“. I have no problem about the fact that I use old econometric models, but I had the feeling that things aren’t that easy. I can understand that we might have problems when we do have a lot of features (I am still working on that, I’ll get back to this point soon), but I have the feeling that I can still capture interactions, and non-linearities with standard econometric models as well as any machine learning algorithm. Just to illustrate, consider the following ‘model where is (just to illustrate) > n <- 5000 > rtf <- function(x1, x2) { sin(x1+x2)/(x1+x2) } > xgrid <- seq(1,6,length=31) > ygrid <- seq(1,6,length=31) > zgrid <- outer(xgrid,ygrid,rtf) > persp(xgrid,ygrid,zgrid,theta=30, phi=30, + col="green", ticktype="detailed",shade=TRUE) # Modèles de prévision, fourre-tout Dans le dernier cours de modèle de prévision, la semaine passée, nous avions passé un peu de temps sur l’étude des points aberrants et des points influents. Tout est expliqué dans les slides (avec les codes) donc je ne reviendrais pas dessus. Je pourrais juste évoquer quelques lignes de codes utilisées pour voir l’impact d’une observation sur la régression, en enlevant l’observation de la base, et en regardant ce que ça donne sur la prévision > plot(cars) > text(cars[c(23,49),1],cars[c(23,49),2]+4,c(23,49)) > abline(lm(dist~speed,data=cars)) > abline(lm(dist~speed,data=cars[-23,]),col="red") > abline(lm(dist~speed,data=cars[-49,]),col="blue") On a ensuite vu les problèmes de choix de modèles, et les méthodes de pénalisation (pénalisation du R2, et de la log-vraisemblance). On a ensuite vu les méthodes automatiques de choix de modèle, > US=read.table("http://freakonometrics.free.fr/US.txt",sep=";") > reg=lm(Murder~.,data=US) > regs=step(reg) Start: AIC=63.01 Murder ~ Population + Income + Illiteracy + Life.Exp + HS.Grad + Frost + Area Df Sum of Sq RSS AIC - Income 1 0.236 128.27 61.105 - HS.Grad 1 0.973 129.01 61.392 <none> 128.03 63.013 - Area 1 7.514 135.55 63.865 - Illiteracy 1 8.299 136.33 64.154 - Frost 1 9.260 137.29 64.505 - Population 1 25.719 153.75 70.166 - Life.Exp 1 127.175 255.21 95.503 Step: AIC=61.11 Murder ~ Population + Illiteracy + Life.Exp + HS.Grad + Frost + Area Df Sum of Sq RSS AIC - HS.Grad 1 0.763 129.03 59.402 <none> 128.27 61.105 - Area 1 7.310 135.58 61.877 - Illiteracy 1 8.715 136.98 62.392 - Frost 1 9.345 137.61 62.621 - Population 1 27.142 155.41 68.702 - Life.Exp 1 127.500 255.77 93.613 Step: AIC=59.4 Murder ~ Population + Illiteracy + Life.Exp + Frost + Area Df Sum of Sq RSS AIC <none> 129.03 59.402 - Illiteracy 1 8.723 137.75 60.672 - Frost 1 11.030 140.06 61.503 - Area 1 15.937 144.97 63.225 - Population 1 26.415 155.45 66.714 - Life.Exp 1 140.391 269.42 94.213 Enfin, on a aussi vu le passage au logarithme sur la variable d’intérêt, avec un modèle log-normal, > plot(cars$speed,cars$dist) > regln=lm(log(dist)~speed,data=cars) > summary(regln) Call: lm(formula = log(dist) ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -1.46604 -0.20800 -0.01683 0.24080 1.01519 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.67612 0.19614 8.546 3.34e-11 *** speed 0.12077 0.01206 10.015 2.41e-13 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.4463 on 48 degrees of freedom Multiple R-squared: 0.6763, Adjusted R-squared: 0.6696 F-statistic: 100.3 on 1 and 48 DF, p-value: 2.413e-13 > beta=coefficients(regln) > for(x in seq(2,30,by=.1)){ + zhat=beta+beta*x + yhat=exp(zhat+1/2*0.4463^2) + points(x,exp(zhat),col="red") + points(x,yhat,col="blue") + } > abline(lm(dist~speed,data=cars)) Pour comparer les modèles (linéaire versus log-linéaire), on avait suggéré de comparer les résidus des deux modèles, > U=cars$dist - exp(predict(regln)+1/2*0.4463^2)

> sd(U)
 17.12949
> sd(residuals(lm(dist~speed,data=cars)))
 15.22184

Ce soir, on fini cette première partie du cours, sur les données individuelles.

# Régression linéaire, quelques codes

Un rapide billet pour mettre en ligne les codes utilisés la semaine passée, complétant les codes des transparents. On travaille toujours sur la même base, ou on cherche à prévoir une distance de freinage d’un véhicule, tenant compte de la vitesse du véhicule.

> plot(cars)
> reg=lm(dist~speed,data=cars)
> summary(reg)

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

Residuals:
Min      1Q  Median      3Q     Max
-29.069  -9.525  -2.272   9.215  43.201

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791     6.7584  -2.601   0.0123 *
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared:  0.6511,	Adjusted R-squared:  0.6438
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Pour faire plusieurs prévisions, à la main, on peut utiliser le code suivant (la boucle permet de faire des prévisions pour plusieurs valeurs)

> for(x in seq(3,30,by=.25)){
+ b0=coef(reg)
+ b1=coef(reg)
+ Yx=b0+b1*x
+ V=vcov(reg)
+ Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2
+ IC1=Yx+c(-1,+1)*1.96*sqrt(Vx)
+ s=summary(reg)$sigma + IC2=Yx+c(-1,+1)*1.96*s + points(x,Yx,pch=19,col="red") + points(c(x,x),IC1,pch=3,col="blue") + points(c(x,x),IC2,pch=3,col="purple")} On avait ensuite fait une régression linéaire sur une sous-base, avec 20 observations tirées au hasard > I=sample(1:50,size=20) > reg=lm(dist~speed,data=cars[I,]) Le but était de visualiser l’impact du nombre d’observation sur la qualité de la régression > summary(reg) Call: lm(formula = dist ~ speed, data = cars[I, ]) Residuals: Min 1Q Median 3Q Max -23.529 -7.998 -5.394 11.634 39.348 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -20.7408 9.4639 -2.192 0.0418 * speed 4.2247 0.6129 6.893 1.91e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 16.62 on 18 degrees of freedom Multiple R-squared: 0.7252, Adjusted R-squared: 0.71 F-statistic: 47.51 on 1 and 18 DF, p-value: 1.91e-06 > for(x in seq(3,30,by=.25)){ + b0=coef(reg) + b1=coef(reg) + Yx=b0+b1*x + V=vcov(reg) + Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2 + IC=Yx+c(-1,+1)*1.96*sqrt(Vx) + points(x,Yx,pch=19,col="purple") + points(c(x,x),IC,pch=3,col="green")} Notons qu’il est possible d’utiliser des fonctions de R pour faire des prévisions, avec des intervalles de confiance > predict(reg, + newdata=data.frame(speed=c(15,25)),interval= "confidence") fit lwr upr 1 42.62976 34.75450 50.50502 2 84.87677 68.92746 100.82607 > predict(reg, + newdata=data.frame(speed=15),interval= "prediction") fit lwr upr 1 42.62976 6.836077 78.42344 Quand on a plus d’une variable explicative, c’est plus compliqué de “visualiser” la régression > chicago=read.table("http://freakonometrics.free.fr/chicago.txt", + header=TRUE,sep=";") > Y=chicago$Fire
>  X1=chicago$X_1 > X2=chicago$X_2
>  X3=chicago$X_3 > base=data.frame(Y,X1,X2,X3) > plot(X2,X3) > reg=lm(Y~X2+X3,data=base) > y=function(x2,x3) predict(reg,newdata=data.frame(X2=x2,X3=x3)) > VX2=seq(0,80) > VX3=seq(5,25) > VY=outer(VX2,VX3,y) > image(VX2,VX3,VY) > contour(VX2,VX3,VY,add=TRUE) qui correspond à un plan de régression > persp(VX2,VX3,VY,theta=30,ticktype=detailed) On reviendra plus en détails sur ce point, mais il est possible de faire des régressions non linéaires assez facilement, à partir de ce modèle linéaire. On avait commencé par un modèle linéaire sur le logarithme de la distance > plot(cars$speed,log(cars$dist)) > reg1=lm(log(dist)~speed,data=cars) > abline(reg1,col="red") (on le verra, ce n’est pas fini, car on n’a pas ici de prévision sur la distance, juste sur son logarithme… mais promis, on en reparlera) ou sur la racine carrée > plot(cars$speed,sqrt(cars$dist)) > reg1=lm(sqrt(dist)~speed,data=cars) > abline(reg1,col="red") Au lieu de transformer la variable d’intérêt, on peut aussi transformer la variable explicative. On peut pendre des puissances, ou des fonctions simples, mais aussi mettre des ruptures. On avait commencé par une variable indicatrice, > plot(cars$speed,cars$dist) > s=10 > abline(v=s,col="green") > regs=lm(dist~speed+I(speed>s),data=cars) > summary(regs) Call: lm(formula = dist ~ speed + I(speed > s), data = cars) Residuals: Min 1Q Median 3Q Max -29.472 -9.559 -2.088 7.456 44.412 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.2964 6.7709 -2.555 0.0139 * speed 4.3140 0.5762 7.487 1.5e-09 *** I(speed > s)TRUE -7.5116 7.8511 -0.957 0.3436 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.39 on 47 degrees of freedom Multiple R-squared: 0.6577, Adjusted R-squared: 0.6432 F-statistic: 45.16 on 2 and 47 DF, p-value: 1.141e-11 Mais on peut aussi mettre des fonctions afin d’avoir un modèle linéaire par morceaux, tout en étant continu > plot(cars) > s=15 > abline(v=s,col="green") > positive=function(x) ifelse(x>0,x,0) > regs=lm(dist~speed+positive(speed-s),data=cars) > summary(regs) Call: lm(formula = dist ~ speed + positive(speed - s), data = cars) Residuals: Min 1Q Median 3Q Max -29.502 -9.513 -2.413 5.195 45.391 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -7.6519 10.6254 -0.720 0.47500 speed 3.0186 0.8627 3.499 0.00103 ** positive(speed - s) 1.7562 1.4551 1.207 0.23350 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.31 on 47 degrees of freedom Multiple R-squared: 0.6616, Adjusted R-squared: 0.6472 F-statistic: 45.94 on 2 and 47 DF, p-value: 8.761e-12 On a ici une rupture, mais on pourrait imaginer en avoir plusieurs > nouvellebase=data.frame(speed=5:25) > y=predict(regs,newdata=nouvellebase) > lines(5:25,y,col="red") > > plot(cars$speed,cars$dist) > s1=10 > s2=20 > abline(v=c(s1,s2),col="green") > positive=function(x) ifelse(x>0,x,0) > regs=lm(dist~speed+positive(speed-s1)+positive(speed-s2),data=cars) > summary(regs) Call: lm(formula = dist ~ speed + positive(speed - s1) + positive(speed - s2), data = cars) Residuals: Min 1Q Median 3Q Max -24.374 -9.475 -2.625 6.639 43.914 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -7.6305 16.2941 -0.468 0.6418 speed 3.0630 1.8238 1.679 0.0998 . positive(speed - s1) 0.2087 2.2453 0.093 0.9263 positive(speed - s2) 4.2812 2.2843 1.874 0.0673 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15 on 46 degrees of freedom Multiple R-squared: 0.6821, Adjusted R-squared: 0.6613 F-statistic: 32.89 on 3 and 46 DF, p-value: 1.643e-11 Comme vu en cours, le test de significativité des deux derniers coefficients ne veut pas dire que la pente est nulle, mais qu’elle est significativement différente de cette obtenue sur la zone de gauche (avant les deux seuils). # Some heuristics about local regression and kernel smoothing In a standard linear model, we assume that . Alternatives can be considered, when the linear assumption is too strong. • Polynomial regression A natural extension might be to assume some polynomial function, Again, in the standard linear model approach (with a conditional normal distribution using the GLM terminology), parameters can be obtained using least squares, where a regression of on is considered. Even if this polynomial model is not the real one, it might still be a good approximation for . Actually, from Stone-Weierstrass theorem, if is continuous on some interval, then there is a uniform approximation of by polynomial functions. Just to illustrate, consider the following (simulated) dataset set.seed(1) n=10 xr = seq(0,n,by=.1) yr = sin(xr/2)+rnorm(length(xr))/2 db = data.frame(x=xr,y=yr) plot(db) with the standard regression line reg = lm(y ~ x,data=db) abline(reg,col="red") Consider some polynomial regression. If the degree of the polynomial function is large enough, any kind of pattern can be obtained, reg=lm(y~poly(x,5),data=db) But if the degree is too large, then too many ‘oscillations’ are obtained, reg=lm(y~poly(x,25),data=db) and the estimation might be be seen as no longer robust: if we change one point, there might be important (local) changes plot(db) attach(db) lines(xr,predict(reg),col="red",lty=2) yrm=yr;yrm=yr-2 regm=lm(yrm~poly(xr,25)) lines(xr,predict(regm),col="red") • Local regression Actually, if our interest is to have locally a good approximation of , why not use a local regression? This can be done easily using a weighted regression, where, in the least square formulation, we consider (it is possible to consider weights in the GLM framework, but let’s keep that for another post). Two comments here: • here I consider a linear model, but any polynomial model can be considered. Even a constant one. In that case, the optimization problem is which can be solve explicitly, since • so far, nothing was mentioned about the weights. The idea is simple, here: if you can a good prediction at point , then should be proportional to some distance between and : if is too far from , then it should not have to much influence on the prediction. For instance, if we want to have a prediction at some point , consider . With this model, we remove observations too far away, Actually, here, it is the same as reg=lm(yr~xr,subset=which(abs(xr-x0)<1) A more general idea is to consider some kernel function that gives the shape of the weight function, and some bandwidth (usually denoted h) that gives the length of the neighborhood, so that This is actually the so-called Nadaraya-Watson estimator of function . In the previous case, we did consider a uniform kernel , with bandwith , But using this weight function, with a strong discontinuity may not be the best idea… Why not a Gaussian kernel, This can be done using fitloc0 = function(x0){ w=dnorm((xr-x0)) reg=lm(y~1,data=db,weights=w) return(predict(reg,newdata=data.frame(x=x0)))} On our dataset, we can plot ul=seq(0,10,by=.01) vl0=Vectorize(fitloc0)(ul) u0=seq(-2,7,by=.01) linearlocalconst=function(x0){ w=dnorm((xr-x0)) plot(db,cex=abs(w)*4) lines(ul,vl0,col="red") axis(3) axis(2) reg=lm(y~1,data=db,weights=w) u=seq(0,10,by=.02) v=predict(reg,newdata=data.frame(x=u)) lines(u,v,col="red",lwd=2) abline(v=c(0,x0,10),lty=2) } linearlocalconst(2) Here, we want a local regression at point 2. The horizonal line below is the regression (the size of the point is proportional to the wieght). The curve, in red, is the evolution of the local regression Let us use an animation to visualize the construction of the curve. One can use library(animate) but for some reasons, I cannot install the package easily on Linux. And it is not a big deal. We can still use a loop to generate some graphs vx0=seq(1,9,by=.1) vx0=c(vx0,rev(vx0)) graphloc=function(i){ name=paste("local-reg-",100+i,".png",sep="") png(name,600,400) linearlocalconst(vx0[i]) dev.off()} for(i in 1:length(vx0)) graphloc(i) and then, in a terminal, I simply use  convert -delay 25 /home/freak/local-reg-1*.png /home/freak/local-reg.gif Of course, it is possible to consider a linear model, locally, fitloc1 = function(x0){ w=dnorm((xr-x0)) reg=lm(y~poly(x,degree=1),data=db,weights=w) return(predict(reg,newdata=data.frame(x=x0)))} or even a quadratic (local) regression, fitloc2 = function(x0){ w=dnorm((xr-x0)) reg=lm(y~poly(x,degree=2),data=db,weights=w) return(predict(reg,newdata=data.frame(x=x0)))} Of course, we can change the bandwidth To conclude the technical part this post, observe that, in practise, we have to choose the shape of the weight function (the so-called kernel). But there are (simple) technique to select the “optimal” bandwidth h. The idea of cross validation is to consider where is the prediction obtained using a local regression technique, with bandwidth . And to get a more accurate (and optimal) bandwith is obtained using a model estimated on a sample where the ith observation was removed. But again, that is not the main point in this post, so let’s keep that for another one… Perhaps we can try on some real data? Inspired from a great post on http://f.briatte.org/teaching/ida/092_smoothing.html, by François Briatte, consider the Global Episode Opinion Survey, from some TV show, http://geos.tv/index.php/index?sid=189 , like Dexter. library(XML) library(downloader) file = "geos-tww.csv" html = htmlParse("http://www.geos.tv/index.php/list?sid=189&collection=all") html = xpathApply(html, "//table[@id='collectionTable']")[] data = readHTMLTable(html) data = data[,-3] names(data)=c("no",names(data)[-1]) data=data[-(61:64),] Let us reshape the dataset, data$no = 1:96
data$mu = as.numeric(substr(as.character(data$Mean), 0, 4))
data$se = sd(data$mu,na.rm=TRUE)/sqrt(as.numeric(as.character(data$Count))) data$season = 1 + (data$no - 1)%/%12 data$season = factor(data$season) plot(data$no,data$mu,ylim=c(6,10)) segments(data$no,data$mu-1.96*data$se,
data$no,data$mu+1.96*data$se,col="light blue") As done by François, we compute some kind of standard error, just to reflect uncertainty. But we won’t really use it. plot(data$no,data$mu,ylim=c(6,10)) abline(v=12*(0:8)+.5,lty=2) for(s in 1:8){reg=lm(mu~no,data=db,subset=season==s) lines((s-1)*12+1:12,predict(reg)[1:12],col="red") } Henre, we assume that all seasons should be considered as completely independent… which might not be a great assumption. db = data NW = ksmooth(db$no,db$mu,kernel = "normal",bandwidth=5) plot(data$no,data$mu) lines(NW,col="red") We can try to look the curve with a larger bandwidth. The problem is that there is a missing value, at the end. If we (arbitrarily) fill it, we can run a kernel regression, db$mu=7
NW = ksmooth(db$no,db$mu,kernel = "normal",bandwidth=12)
plot(data$no,data$mu,ylim=c(6,10))
lines(NW,col="red") # Non-observable vs. observable heterogeneity factor

This morning, in the ACT2040 class (on non-life insurance), we’ve discussed the difference between observable and non-observable heterogeneity in ratemaking (from an economic perspective). To illustrate that point (we will spend more time, later on, discussing observable and non-observable risk factors), we looked at the following simple example. Let  denote the height of a person. Consider the following dataset

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")

There is a small typo in the dataset, so let us make manual changes here

> Davis[12,c(2,3)]=Davis[12,c(3,2)]

Here, the variable of interest is the height of a given person,

> X=Davis$height  If we look at the histogram, we have > hist(X,col="light green", border="white",proba=TRUE,xlab="",main="") Can we assume that we have a Gaussian distribution ? Maybe not… Here, if we fit a Gaussian distribution, plot it, and add a kernel based estimator, we get > (param <- fitdistr(X,"normal")$estimate)
> f1 <- function(x) dnorm(x,param,param)
> x=seq(100,210,by=.2)
> lines(x,f1(x),lty=2,col="red")
> lines(density(X)) If you look at that black line, you might think of a mixture, i.e. something like

(using standard mixture notations). Mixture are obtained when we have a non-observable heterogeneity factor: with probability , we have a random variable  (call it type ), and with probability , a random variable  (call it type ). So far, nothing new. And we can fit such a mixture distribution, using e.g.

> library(mixtools)
> mix <- normalmixEM(X)
number of iterations= 335
> (param12 <- c(mix$lambda,mix$mu,mix$sigma))  0.4002202 178.4997298 165.2703616 6.3561363 5.9460023  If we plot that mixture of two Gaussian distributions, we get > f2 <- function(x){ param12*dnorm(x,param12,param12) + (1-param12)*dnorm(x,param12,param12) } > lines(x,f2(x),lwd=2, col="red") lines(density(X)) Not bad. Actually, we can try to maximize the likelihood with our own codes, > logdf <- function(x,parameter){ + p <- parameter + m1 <- parameter + s1 <- parameter + m2 <- parameter + s2 <- parameter + return(log(p*dnorm(x,m1,s1)+(1-p)*dnorm(x,m2,s2))) + } > logL <- function(parameter) -sum(logdf(X,parameter)) > Amat <- matrix(c(1,-1,0,0,0,0, + 0,0,0,0,1,0,0,0,0,0,0,0,0,1), 4, 5) > bvec <- c(0,-1,0,0) > constrOptim(c(.5,160,180,10,10), logL, NULL, ui = Amat, ci = bvec)$par

   0.5996263 165.2690084 178.4991624   5.9447675   6.3564746

Here, we include some constraints, to insurance that the probability belongs to the unit interval, and that the variance parameters remain positive. Note that we have something close to the previous output.

Let us try something a little bit more complex now. What if we assume that the underlying distributions have the same variance, namely

In that case, we have to use the previous code, and make small changes,

> logdf <- function(x,parameter){
+ p <- parameter
+ m1 <- parameter
+ s1 <- parameter
+ m2 <- parameter
+ s2 <- parameter
+ return(log(p*dnorm(x,m1,s1)+(1-p)*dnorm(x,m2,s2)))
+ }
> logL <- function(parameter) -sum(logdf(X,parameter))
> Amat <- matrix(c(1,-1,0,0,0,0,0,0,0,0,0,1), 3, 4)
> bvec <- c(0,-1,0)
> (param12c= constrOptim(c(.5,160,180,10), logL, NULL, ui = Amat, ci = bvec)$par)  0.6319105 165.6142824 179.0623954 6.1072614 This is what we can do if we cannot observe the heterogeneity factor. But wait… we actually have some information in the dataset. For instance, we have the sex of the person. Now, if we look at histograms of height per sex, and kernel based density estimator of the height, per sex, we have So, it looks like the height for male, and the height for female are different. Maybe we can use that variable, that was actually observed, to explain the heterogeneity in our sample. Formally, here, the idea is to consider a mixture, with an observable heterogeneity factor: the sex, We now have interpretation of what we used to call class  and  previously: male and female. And here, estimating parameters is quite simple, > (pM <- mean(sex=="M"))  0.44 > (paramF <- fitdistr(X[sex=="F"],"normal")$estimate)
mean         sd
164.714286   5.633808
>  (paramM <- fitdistr(X[sex=="M"],"normal")$estimate) mean sd 178.011364 6.404001 And if we plot the density, we have > f4 <- function(x) pM*dnorm(x,paramM,paramM)+(1-pM)*dnorm(x,paramF,paramF) > lines(x,f4(x),lwd=3,col="blue") What if, once again, we assume identical variance? Namely, the model becomes Then a natural idea to derive an estimator for the variance, based on previous computations, is to use The code is here > s=sqrt((sum((height[sex=="M"]-paramM)^2)+sum((height[sex=="F"]-paramF)^2))/(nrow(Davis)-2)) > s  6.015068 and again, it is possible to plot the associated density, > f5 <- function(x) pM*dnorm(x,paramM,s)+(1-pM)*dnorm(x,paramF,s) > lines(x,f5(x),lwd=3,col="blue") Now, if we think a little about what we’ve just done, it is simply a linear regression on a factor, the sex of the person, where . And indeed, if we run the code to estimate this linear model, > summary(lm(height~sex,data=Davis)) Call: lm(formula = height ~ sex, data = Davis) Residuals: Min 1Q Median 3Q Max -16.7143 -3.7143 -0.0114 4.2857 18.9886 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 164.7143 0.5684 289.80 <2e-16 *** sexM 13.2971 0.8569 15.52 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 6.015 on 198 degrees of freedom Multiple R-squared: 0.5488, Adjusted R-squared: 0.5465 F-statistic: 240.8 on 1 and 198 DF, p-value: < 2.2e-16 we get the same estimators for the means and the variance as the ones obtained previously. So, as mentioned this morning in class, if you have a non-observable heterogeneity factor, we can use a mixture model to fit a distribution, but if you can get a proxy of that factor, that is observable, then you can run a regression. But most of the time, that observable variable is just a proxy of a non-observable one… # Linear regression from a contingency table This morning, Benoit sent me an email, about an exercise he found in an econometric textbook, about linear regression. Consider the following dataset, Here, variable X denotes the income, and Y the expenses. The goal was to fit a linear regression (actually, in the email, it was mentioned that we should try to fit an heteroscedastic model, but let us skip this part). So Benoit’s question was more or less: how do you fit a linear regression from a contingency table? Usually, when I got an email on Saturday morning, I try to postpone. But the kids had their circus class, so I had some time to answer. And this did not look like a complex puzzle… Let us import this dataset in R, so that we can start playing > df=read.table("http://freakonometrics.free.fr/baseexo.csv",sep=";",header=TRUE) > M=as.matrix(df[,2:ncol(df)]) > M[is.na(M)]<-0 > M X14 X19 X21 X23 X25 X27 X29 X31 X33 X35 [1,] 74 13 7 1 0 0 0 0 0 0 [2,] 6 4 2 7 4 0 0 0 0 0 [3,] 2 3 2 2 4 0 0 0 0 0 [4,] 1 1 2 3 3 2 0 0 0 0 [5,] 2 0 1 3 2 0 6 0 0 0 [6,] 2 0 2 1 0 0 1 2 1 0 [7,] 0 0 0 2 0 0 1 1 3 0 [8,] 0 1 0 1 0 0 0 0 2 0 [9,] 0 0 0 0 1 1 0 1 0 1 The first idea I had was to use those counts as weights. Weighted least squares should be perfect. The dataset is built from this matrix, > W=as.vector(M) > x=df[,1] > X=rep(x,ncol(M)) > y=as.numeric(substr(names(df)[-1],2,3)) > Y=rep(y,each=nrow(M)) > base1=data.frame(X1=X,Y1=Y,W1=W) Here we have > head(base1,10) X1 Y1 W1 1 16 14 74 2 23 14 6 3 25 14 2 4 27 14 1 5 29 14 2 6 31 14 2 7 33 14 0 8 35 14 0 9 37 14 0 10 16 19 13 The regression is the following, > summary(reg1) Call: lm(formula = Y1 ~ X1, data = base1, weights = W1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.35569 2.03022 2.145 0.038 * X1 0.68263 0.09016 7.572 3.04e-09 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 7.892 on 40 degrees of freedom Multiple R-squared: 0.589, Adjusted R-squared: 0.5787 F-statistic: 57.33 on 1 and 40 DF, p-value: 3.038e-09 It looks like the output is the same as what Benoit found, so we should be happy. Now, I had a second thought. Why not create the implied dataset. Using replicates, we should be able to create the dataset that was used to get this contingency table, > vX=vY=rep(NA,sum(W)) > sumW=c(0,cumsum(W)) > for(i in 1:length(W)){ + if(W[i]>0){ + vX[(1+sumW[i]):sumW[i+1]]=X[i] + vY[(1+sumW[i]):sumW[i+1]]=Y[i] + }} > base2=data.frame(X2=vX,Y2=vY) Here, the dataset is much larger, and there is no weight, > tail(base2,10) X2 Y2 172 31 31 173 33 31 174 37 31 175 31 33 176 33 33 177 33 33 178 33 33 179 35 33 180 35 33 181 37 35 If we run a linear regression on this dataset, we obtain > summary(reg2) Call: lm(formula = Y2 ~ X2, data = base2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.35569 0.95972 4.538 1.04e-05 *** X2 0.68263 0.04262 16.017 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 3.731 on 179 degrees of freedom Multiple R-squared: 0.589, Adjusted R-squared: 0.5867 F-statistic: 256.5 on 1 and 179 DF, p-value: < 2.2e-16 If we compare the two regressions, we have > rbind(coefficients(summary(reg1)), + coefficients(summary(reg2))) Estimate Std. Error t value Pr(>|t|) (Intercept) 4.3556857 2.03021637 2.145429 3.804237e-02 X1 0.6826296 0.09015771 7.571506 3.038443e-09 (Intercept) 4.3556857 0.95972279 4.538483 1.036711e-05 X2 0.6826296 0.04261930 16.016913 2.115373e-36 The estimators are exactly the same (which does not surprise me), but standard deviation (ans significance levels) are quite different. And to be honest, I find that surprising. Which approach here is the most legitimate (since they are finally not equivalent)? # Why pictures are so important when modeling data? (bis repetita) Consider the following regression summary, Call: lm(formula = y1 ~ x1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0001 1.1247 2.667 0.02573 * x1 0.5001 0.1179 4.241 0.00217 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.237 on 9 degrees of freedom Multiple R-squared: 0.6665, Adjusted R-squared: 0.6295 F-statistic: 17.99 on 1 and 9 DF, p-value: 0.00217 obtained from > data(anscombe) > reg1=lm(y1~x1,data=anscombe) Can we say something if we look (only) at that output ? The intercept is significatively non-null, as well as the slope, the is large (66%). It looks like we do have a nice model here. And in a perfect world, we might hope that data are coming from this kind of dataset, But it might be possible to have completely different kinds of patterns. Actually, four differents sets of data are coming from Anscombe (1973). And that all those datasets are somehow equivalent: the ‘s have the same mean, and the same variance > apply(anscombe[,1:4],2,mean) x1 x2 x3 x4 9 9 9 9 > apply(anscombe[,1:4],2,var) x1 x2 x3 x4 11 11 11 11 and so are the ‘s > apply(anscombe[,5:8],2,mean) y1 y2 y3 y4 7.500909 7.500909 7.500000 7.500909 > apply(anscombe[,5:8],2,var) y1 y2 y3 y4 4.127269 4.127629 4.122620 4.123249 Further, observe also that the correlation between the ‘s and the ‘s is the same > cor(anscombe)[1:4,5:8] y1 y2 y3 y4 x1 0.8164205 0.8162365 0.8162867 -0.3140467 x2 0.8164205 0.8162365 0.8162867 -0.3140467 x3 0.8164205 0.8162365 0.8162867 -0.3140467 x4 -0.5290927 -0.7184365 -0.3446610 0.8165214 > diag(cor(anscombe)[1:4,5:8])  0.8164205 0.8162365 0.8162867 0.8165214 which yields the same regression line (intercept and slope) > cbind(coef(reg1),coef(reg2),coef(reg3),coef(reg4)) [,1] [,2] [,3] [,4] (Intercept) 3.0000909 3.000909 3.0024545 3.0017273 x1 0.5000909 0.500000 0.4997273 0.4999091 But there is more. Much more. For instance, we always have the standard deviation for residuals > c(summary(reg1)$sigma,summary(reg2)$sigma, + summary(reg3)$sigma,summary(reg4)$sigma)  1.236603 1.237214 1.236311 1.235695 Thus, all regressions here have the same R2 > c(summary(reg1)$r.squared,summary(reg2)$r.squared, + summary(reg3)$r.squared,summary(reg4)$r.squared)  0.6665425 0.6662420 0.6663240 0.6667073 Finally, Fisher’s F statistics is also (almost) the same. + c(summary(reg1)$fstatistic,summary(reg2)$fstatistic, + summary(reg3)$fstatistic,summary(reg4)$fstatistic) value value value value 17.98994 17.96565 17.97228 18.00329 Thus, with the following datasets, we have the same prediction (and the same confidence intervals). Consider for instance the second dataset (the first one being mentioned above), > reg2=lm(y2~x2,data=anscombe) The output is here exactly the same as the one we had above > summary(reg2) Call: lm(formula = y2 ~ x2, data = anscombe) Residuals: Min 1Q Median 3Q Max -1.9009 -0.7609 0.1291 0.9491 1.2691 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.001 1.125 2.667 0.02576 * x2 0.500 0.118 4.239 0.00218 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.237 on 9 degrees of freedom Multiple R-squared: 0.6662, Adjusted R-squared: 0.6292 F-statistic: 17.97 on 1 and 9 DF, p-value: 0.002179 Here, the perfect model is the one obtained with a quadratic regression. > reg2b=lm(y2~x2+I(x2^2),data=anscombe) > summary(reg2b) Call: lm(formula = y2 ~ x2 + I(x2^2), data = anscombe) Residuals: Min 1Q Median 3Q Max -0.0013287 -0.0011888 -0.0006294 0.0008741 0.0023776 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -5.9957343 0.0043299 -1385 <2e-16 *** x2 2.7808392 0.0010401 2674 <2e-16 *** I(x2^2) -0.1267133 0.0000571 -2219 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.001672 on 8 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 7.378e+06 on 2 and 8 DF, p-value: < 2.2e-16 Consider now the third one > reg3=lm(y3~x3,data=anscombe) i.e. > summary(reg3) Call: lm(formula = y3 ~ x3, data = anscombe) Residuals: Min 1Q Median 3Q Max -1.1586 -0.6146 -0.2303 0.1540 3.2411 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0025 1.1245 2.670 0.02562 * x3 0.4997 0.1179 4.239 0.00218 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.236 on 9 degrees of freedom Multiple R-squared: 0.6663, Adjusted R-squared: 0.6292 F-statistic: 17.97 on 1 and 9 DF, p-value: 0.002176 This time, the linear model could have been perfect. The problem is one outlier. If we remove it, we have > reg3b=lm(y3~x3,data=anscombe[-3,]) > summary(reg3b) Call: lm(formula = y3 ~ x3, data = anscombe[-3, ]) Residuals: Min 1Q Median 3Q Max -0.0041558 -0.0022240 0.0000649 0.0018182 0.0050649 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.0056494 0.0029242 1370 <2e-16 *** x3 0.3453896 0.0003206 1077 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.003082 on 8 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: 1 F-statistic: 1.161e+06 on 1 and 8 DF, p-value: < 2.2e-16 Finally consider > reg4=lm(y4~x4,data=anscombe) This time, there is an other kind of outlier, in ‘s, but again, the regression is exactly the same, > summary(reg4) Call: lm(formula = y4 ~ x4, data = anscombe) Residuals: Min 1Q Median 3Q Max -1.751 -0.831 0.000 0.809 1.839 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0017 1.1239 2.671 0.02559 * x4 0.4999 0.1178 4.243 0.00216 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.236 on 9 degrees of freedom Multiple R-squared: 0.6667, Adjusted R-squared: 0.6297 F-statistic: 18 on 1 and 9 DF, p-value: 0.002165 The graph is here So clearly, looking at the summary of a regression does not tell us anything… This is why we do spend some time on diagnostic, looking at graphs with the errors (the graphs above could be obtained only with one explanatory variable, while errors can be studied in any dimension): everything can be seen on thise graphs. E.g. for the first dataset, or the second one the third one or the fourth one, # Le passage au log dans les modèles linéaires Un billet rapide pour compléter et illustrer le passage au log dans un modèle linéaire (que l’on abordera cette semaine en cours). Le point de départ est le modèle linéaire, où on suppose que, conditionnellement à suit une loi normale. Pour rappel, si on a une loi normale, , alors et . Les intervalles de confiance à 90% et 95% sont symétriques par rapport à la moyenne (qui est aussi la médiane, soit dit en passant), Dans un modèle Gaussien avec homoscédasiticité, i.e. alors que . On a alors les bandes de confiance suivantes, pour un modèle de régression linéaire, Bon, maintenant, que se passe-t-il si on prend l’exponentiel ? Pour la loi normale, rappelons que l’on obtient une loi lognormale, i.e. , les deux paramètres étant liés à la loi normale sous jacente, car désormais alors que Graphiquement, on a la loi suivante, avec les intervalles de confiance à 90% et 95% représentés ci-dessous. Le point noir est alors que le point bleu est l’espérance de la loi lognormale. On notera que le quantile de la loi log-normale est l’exponentiel du quantile de la loi normale. En effet, si alors . En particulier, n’est pas la moyenne de , mais la médiane (puisque était la médiane de ). Mais il n’est pas rare de voir utilisé un intervalle de confiance de la forme qui est la forme classique de l’intervalle de confiance Gaussien (symétrique autour de la moyenne). Ici, on aurait les niveaux suivants Notons qu’il n’y a aucune raison ici d’avoir une probabilité d’être dans l’intervalle de confiance obtenu avec les quantiles de la loi normale. Maintenant, si on prend l’exponentiel d’un modèle linéaire (i.e. le logarithme de la variable d’intérêt est modélisé par un modèle linéaire) on a avec une variance (conditionnelle) qui dépend de la variable explicative Là encore, le plus naturel est d’utiliser comme bornes de l’intervalle de confiance des quantiles associés à la loi lognormale, mais il n’est pas rare de voir utilisé des intervalles de type Gaussiens, On perd là encore en interprétation car les bornes n’ont plus rien à voir avec les quantiles. # Visualization in regression analysis Visualization is a key to success in regression analysis. This is one of the (many) reasons I am also suspicious when I read an article with a quantitative (econometric) analysis without any graph. Consider for instance the following dataset, obtained from http://data.worldbank.org/, with, for each country, the GDP per capita (in some common currency) and the infant mortality rate (deaths before the age of 5), > library(gdata) > XLS1=read.xls("http://api.worldbank.org/datafiles /NY.GDP.PCAP.PP.CD_Indicator_MetaData_en_EXCEL.xls", sheet = 1) > data1=XLS1[-(1:28),c("Country.Name","Country.Code","X2010")] > names(data1)="GDP" > XLS2=read.xls("http://api.worldbank.org/datafiles /SH.DYN.MORT_Indicator_MetaData_en_EXCEL.xls", sheet = 1) > data2=XLS2[-(1:28),c("Country.Code","X2010")] > names(data2)="MORTALITY" > data=merge(data1,data2) > head(data) Country.Code Country.Name GDP MORTALITY 1 ABW Aruba NA NA 2 AFG Afghanistan 1207.278 149.2 3 AGO Angola 6119.930 160.5 4 ALB Albania 8817.009 18.4 5 AND Andorra NA 3.8 6 ARE United Arab Emirates 47215.315 7.1 If we estimate a simple linear regression – – we get > regBB=lm(MORTALITY~GDP,data=data) > summary(regBB) Call: lm(formula = MORTALITY ~ GDP, data = data) Residuals: Min 1Q Median 3Q Max -45.24 -29.58 -12.12 16.19 115.83 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 67.1008781 4.1577411 16.139 < 2e-16 *** GDP -0.0017887 0.0002161 -8.278 3.83e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 39.99 on 167 degrees of freedom (47 observations deleted due to missingness) Multiple R-squared: 0.2909, Adjusted R-squared: 0.2867 F-statistic: 68.53 on 1 and 167 DF, p-value: 3.834e-14 We can look at the scatter plot, including the linear regression line, and some confidence bounds, > plot(data$GDP,data$MORTALITY,xlab="GDP per capita", + ylab="Mortality rate (under 5)",cex=.5) > text(data$GDP,data$MORTALITY,data$Country.Name,pos=3)
> x=seq(-10000,100000,length=101)
> y=predict(regBB,newdata=data.frame(GDP=x),
+ interval="prediction",level = 0.9)
> lines(x,y[,1],col="red")
> lines(x,y[,2],col="red",lty=2)
> lines(x,y[,3],col="red",lty=2) We should be able to do a better job here. For instance, if we look at the Box-Cox profile likelihood,

> boxcox(regBB) it looks like taking the logarithm of the mortality rate should be better, i.e. or :

> regLB=lm(log(MORTALITY)~GDP,data=data)
> summary(regLB)

Call:
lm(formula = log(MORTALITY) ~ GDP, data = data)

Residuals:
Min      1Q  Median      3Q     Max
-1.3035 -0.5837 -0.1138  0.5597  3.0583

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.989e+00  7.970e-02   50.05   <2e-16 ***
GDP         -6.487e-05  4.142e-06  -15.66   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7666 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.5949,	Adjusted R-squared: 0.5925
F-statistic: 245.3 on 1 and 167 DF,  p-value: < 2.2e-16

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5,log="y")
> text(data$GDP,data$MORTALITY,data$Country.Name) > x=seq(300,100000,length=101) > y=exp(predict(regLB,newdata=data.frame(GDP=x)))* + exp(summary(regLB)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLB,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLB)$sigma^2) > lines(x,y,col="red",lty=2) > y=qlnorm(.05, meanlog=predict(regLB,newdata=data.frame(GDP=x)), + sdlog=summary(regLB)$sigma^2)
> lines(x,y,col="red",lty=2) on the log scale or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5) on the standard scale. Here we use quantiles of the log-normal distribution to derive confidence intervals.

But why shouldn’t we take also the logarithm of the GDP ? We can fit a model or equivalently .

> regLL=lm(log(MORTALITY)~log(GDP),data=data)
> summary(regLL)

Call:
lm(formula = log(MORTALITY) ~ log(GDP), data = data)

Residuals:
Min       1Q   Median       3Q      Max
-1.13200 -0.38326 -0.07127  0.26610  3.02212

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.50192    0.31556   33.28   <2e-16 ***
log(GDP)    -0.83496    0.03548  -23.54   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5797 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.7684,	Adjusted R-squared: 0.767
F-statistic:   554 on 1 and 167 DF,  p-value: < 2.2e-16

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5,log="xy")
> text(data$GDP,data$MORTALITY,data$Country.Name) > x=exp(seq(1,12,by=.1)) > y=exp(predict(regLL,newdata=data.frame(GDP=x)))* + exp(summary(regLL)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLL,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLL)$sigma^2) > lines(x,y,col="red",lty=2) > y=qlnorm(.05, meanlog=predict(regLL,newdata=data.frame(GDP=x)), + sdlog=summary(regLL)$sigma^2)
> lines(x,y,col="red",lty=2) on the log scales or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5) on the standard scale. If we compare the last two predictions, we have with in blue is the log model, and in red is the log-log model (I did not include the first one for obvious reasons).