Freakonometrics

To content | To menu | To search

Saturday, March 17 2012

De la normalité dans les régressions

Ce matin, Valérie m'a posé par mail une "question naïve" (ce sont ses termes). Et j'aime les "questions naïves" ! Tout d'abord parce que les personnes qui les posent montrent ainsi qu'elles ont vraiment envie de comprendre ! Et il n'y a rien de mieux pour flatter l'égo d'un enseignant que de lui demander d'expliquer des choses ! En plus, ce sont souvent des questions simples, mais qui traduisent derrières des choses plus fondamentales. Et donc que l'on peut souvent interpréter de manière plus complexe. 

Mais je m’égare. La question portait sur les liens entre "normalité des résidus" et "normalité de la variable explicative dans le modèle linéaire". Selon ses termes, "les hypothèses fondamentales sur les résidus peuvent être transférées sur la variable à expliquer" en particulier sur la normalité. Par exemple, Valérie voulait savoir si avoir des résidus Gaussiens rendait la variable explicative Gaussienne.

Pour faire les choses dans l'ordre, effectivement, il y a des liens très étroits entre les résidus et la variable à expliquer. C'est ce que dit Frees (2009) dans deux tableaux, lorsqu'il propose des hypothèses sur les résidus

et en face, des hypothèses sur la variable explicative,

avec quelque part la phrase suivante

Ça semble bien, et j'ai dit des choses proches en cours. Sauf que je ne suis pas (vraiment) d'accord avec l'expression

http://freakonometrics.blog.free.fr/public/perso5/Capture_d_ecran_2012-03-16_a_15.20.14.png

(ni vraiment avec d'autres ponts dans cette seconde table, mais la question de Valérie était sur la normalité). J'aurais plutôt envie de dire que c'est conditionnellement aux variables explicatives que la variable à expliquer est Gaussienne. Plus précisément, je traduirais la première hypothèse sous la forme

la troisième sous la forme

et enfin la cinquième sous la forme

Oui mais comme le note la seconde hypothèse, ça veut dire quoi de conditionner par une constante ? J'avoue que j'ai encore du mal à bien voir comment formaliser ce point, mais j'ai toujours du mal ce qui pourrait légitimer

http://freakonometrics.blog.free.fr/public/perso5/Capture_d_ecran_2012-03-16_a_15.20.14.png

Pour illustrer mon propos, construisons un modèle très simple, où la variable explicative (on va en prendre une seule) prend deux modalités,

n=1000
X=sample(c(0,1,1),size=n,replace=TRUE)
E=rnorm(n)
Y=-1+5*X+E

Les résidus sont Gaussien (je les ai choisis comme tels),

reg=lm(Y~X)
hist(residuals(reg),probability=TRUE,col="light blue") u=seq(-5,10,by=.025) v=dnorm(u,0,1) lines(u,v,lwd=2,col="red")

Par contre, dans ce cas, la loi de  n'est pas une loi Gaussienne: c'est un mélange de lois Gaussiennes,

hist(Y,probability=TRUE,col="light green")
u=seq(-5,10,by=.025)
v=dnorm(u,-1,1)/3+2*dnorm(u,5-1,1)/3
lines(u,v,lwd=2,col="red")

Bref, dire que la variable à expliquer est Gaussienne est une chose à éviter (ça restera vrai d'ailleurs avec les GLM, où, par exemple, dans une régression de Poisson: les comptages ne suivent pas des lois de Poisson). La loi conditionnelle en revanche est Gaussienne. Sauf qu'avec cette seconde hypothèse, il devient délicat de conditionner... A suivre probablement dans les commentaires si certains souhaitent réagir.

Thursday, February 23 2012

Visualization in regression analysis

Visualization is a key to success in regression analysis. This is one of the (many) reasons I am also suspicious when I read an article with a quantitative (econometric) analysis without any graph. Consider for instance the following dataset, obtained from http://data.worldbank.org/, with, for each country, the GDP per capita (in some common currency) and the infant mortality rate (deaths before the age of 5),
> library(gdata)
> XLS1=read.xls("http://api.worldbank.org/datafiles
/NY.GDP.PCAP.PP.CD_Indicator_MetaData_en_EXCEL.xls", sheet = 1)
> data1=XLS1[-(1:28),c("Country.Name","Country.Code","X2010")]
> names(data1)[3]="GDP"
> XLS2=read.xls("http://api.worldbank.org/datafiles
/SH.DYN.MORT_Indicator_MetaData_en_EXCEL.xls", sheet = 1)
> data2=XLS2[-(1:28),c("Country.Code","X2010")]
> names(data2)[2]="MORTALITY"
> data=merge(data1,data2)
> head(data)
Country.Code         Country.Name       GDP MORTALITY
1          ABW                Aruba        NA        NA
2          AFG          Afghanistan  1207.278     149.2
3          AGO               Angola  6119.930     160.5
4          ALB              Albania  8817.009      18.4
5          AND              Andorra        NA       3.8
6          ARE United Arab Emirates 47215.315       7.1

If we estimate a simple linear regression - http://freakonometrics.blog.free.fr/public/perso5/logormal01.gif  - we get
> regBB=lm(MORTALITY~GDP,data=data)
> summary(regBB)
 
Call:
lm(formula = MORTALITY ~ GDP, data = data)
 
Residuals:
Min     1Q Median     3Q    Max
-45.24 -29.58 -12.12  16.19 115.83
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 67.1008781  4.1577411  16.139  < 2e-16 ***
GDP         -0.0017887  0.0002161  -8.278 3.83e-14 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 39.99 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.2909,	Adjusted R-squared: 0.2867
F-statistic: 68.53 on 1 and 167 DF,  p-value: 3.834e-14 
We can look at the scatter plot, including the linear regression line, and some confidence bounds,
> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5)",cex=.5)
> text(data$GDP,data$MORTALITY,data$Country.Name,pos=3)
> x=seq(-10000,100000,length=101)
> y=predict(regBB,newdata=data.frame(GDP=x),
+ interval="prediction",level = 0.9)
> lines(x,y[,1],col="red")
> lines(x,y[,2],col="red",lty=2)
> lines(x,y[,3],col="red",lty=2)

We should be able to do a better job here. For instance, if we look at the Box-Cox profile likelihood,

> boxcox(regBB)

it looks like taking the logarithm of the mortality rate should be better, i.e. http://freakonometrics.blog.free.fr/public/perso5/lognormal02.gif or http://freakonometrics.blog.free.fr/public/perso5/lognormal05.gif:

> regLB=lm(log(MORTALITY)~GDP,data=data)
> summary(regLB)
 
Call:
lm(formula = log(MORTALITY) ~ GDP, data = data)
 
Residuals:
Min      1Q  Median      3Q     Max
-1.3035 -0.5837 -0.1138  0.5597  3.0583
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.989e+00  7.970e-02   50.05   <2e-16 ***
GDP         -6.487e-05  4.142e-06  -15.66   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.7666 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.5949,	Adjusted R-squared: 0.5925
F-statistic: 245.3 on 1 and 167 DF,  p-value: < 2.2e-16
 
> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5,log="y")
> text(data$GDP,data$MORTALITY,data$Country.Name)
> x=seq(300,100000,length=101)
> y=exp(predict(regLB,newdata=data.frame(GDP=x)))*
+ exp(summary(regLB)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLB,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLB)$sigma^2)
> lines(x,y,col="red",lty=2)
> y=qlnorm(.05, meanlog=predict(regLB,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLB)$sigma^2)
> lines(x,y,col="red",lty=2)

on the log scale or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita",
+ ylab="Mortality rate (under 5) log scale",cex=.5)

on the standard scale. Here we use quantiles of the log-normal distribution to derive confidence intervals.

But why shouldn't we take also the logarithm of the GDP ? We can fit a model http://freakonometrics.blog.free.fr/public/perso5/lognormal03.gif or equivalently http://freakonometrics.blog.free.fr/public/perso5/lognormal04.gif.


> regLL=lm(log(MORTALITY)~log(GDP),data=data)
> summary(regLL)
 
Call:
lm(formula = log(MORTALITY) ~ log(GDP), data = data)
 
Residuals:
Min       1Q   Median       3Q      Max
-1.13200 -0.38326 -0.07127  0.26610  3.02212
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.50192    0.31556   33.28   <2e-16 ***
log(GDP)    -0.83496    0.03548  -23.54   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.5797 on 167 degrees of freedom
(47 observations deleted due to missingness)
Multiple R-squared: 0.7684,	Adjusted R-squared: 0.767
F-statistic:   554 on 1 and 167 DF,  p-value: < 2.2e-16
 
> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5,log="xy")
> text(data$GDP,data$MORTALITY,data$Country.Name)
> x=exp(seq(1,12,by=.1))
> y=exp(predict(regLL,newdata=data.frame(GDP=x)))*
+ exp(summary(regLL)$sigma^2/2)
> lines(x,y,col="red")
> y=qlnorm(.95, meanlog=predict(regLL,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLL)$sigma^2)
> lines(x,y,col="red",lty=2)
> y=qlnorm(.05, meanlog=predict(regLL,newdata=data.frame(GDP=x)),
+ sdlog=summary(regLL)$sigma^2)
> lines(x,y,col="red",lty=2)

on the log scales or

> plot(data$GDP,data$MORTALITY,xlab="GDP per capita ",
+ ylab="Mortality rate (under 5)",cex=.5)

on the standard scale. If we compare the last two predictions, we have

with in blue is the log model, and in red is the log-log model (I did not include the first one for obvious reasons).

Thursday, February 2 2012

On linear models with no constant and R2

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

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

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

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

The least-square program is here

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

i.e. the estimate of the slope is

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

Numerically, we obtain

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

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

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

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

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

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

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

i.e.

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

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

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

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

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

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

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

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

Wednesday, May 11 2011

Multivariate probit regression using (direct) maximum likelihood estimators

Consider a random pair http://freakonometrics.blog.free.fr/public/perso3/biv-prob-01.gif of binary responses, i.e. http://freakonometrics.blog.free.fr/public/perso3/biv-prob-02.gif with http://freakonometrics.blog.free.fr/public/perso3/biv-prob-03.gif taking values 1 or 2. Assume that probability http://freakonometrics.blog.free.fr/public/perso3/biv-prob-04.gif can be function of some covariates http://freakonometrics.blog.free.fr/public/perso3/biv-prob-05.gif.
  • The Gaussian vector latent structure
A standard model is based a latent Gaussian structure, i.e. there exists some random vector http://freakonometrics.blog.free.fr/public/perso3/biv-prob-06.gif such that http://freakonometrics.blog.free.fr/public/perso3/biv-prob-07.gif if http://freakonometrics.blog.free.fr/public/perso3/biv-prob-08.gif is lower than a given threshold, and 1 otherwise.
As in standard probit models, assume that
http://freakonometrics.blog.free.fr/public/perso3/biv-prob-09.gif
where we can assume that http://freakonometrics.blog.free.fr/public/perso3/biv-prob-10.gif is a Gaussian random vector. This assumption can be used to derive the likelihood of a sample http://freakonometrics.blog.free.fr/public/perso3/biv-prob-11.gif.
> logV=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+ sigma=matrix(c(1,CORRELATION,CORRELATION,1),2,2)
+ a11=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
+ for(i in 2:nrow(z)){a11=c(a11,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
+ a10=pnorm(z[1,1],sd=sqrt(sigma[1,1]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a10=c(a10,pnorm(z[i,1],sd=sqrt(sigma[1,1]))-pmnorm(z[i,],varcov=sigma))}
+ a01=pnorm(z[1,2],sd=sqrt(sigma[2,2]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a01=c(a01,pnorm(z[i,2],sd=sqrt(sigma[2,2]))-pmnorm(z[i,],varcov=sigma))}
+ a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+ ((Y[,1]==0)&(Y[,2]==1))*log(a01) +
+ ((Y[,1]==1)&(Y[,2]==0))*log(a10) +
+ ((Y[,1]==0)&(Y[,2]==0))*log(a00) )
+ }
> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par
(the code is a bit long since I had trouble working properly with matrices - or more precisely to vectorize my functions - so I used loops... I am sure it is possible to write a better code).
It is possible to generate samples (based on that specific model) to check that we can actually derive proper maximum likelihood estimators,
> library(mnormt)
> set.seed(1)
> n=1000
> r=0.5
> X1=runif(n)
> X2=rnorm(n)
> Y1S=1+5*X1
> Y2S=8-5*X1
> RES=rmnorm(n,mean=c(0,0),varcov=matrix(c(1,r,r,1),2,2))
> YS=cbind(Y1S,Y2S)+RES
> Y1=(YS[,1]>quantile(YS[,1],.5))*1
> Y2=(YS[,2]>quantile(YS[,2],.5))*1
> base=data.frame(i,Y1,Y2,X1,X2,YS)
> head(base)
i Y1 Y2 X1 X2 Y1S Y2S
1 1 0 0 0.2655087 0.07730312 3.177587 5.533884
2 2 0 0 0.3721239 -0.29686864 1.935307 5.089524
3 3 1 0 0.5728534 -1.18324224 4.757848 5.172584
4 4 1 0 0.9082078 0.01129269 4.600029 3.878225
5 5 0 1 0.2016819 0.99160104 2.547362 6.743714
6 6 1 0 0.8983897 1.59396745 5.309974 4.421523

(the two columns on the right are latent observations, that cannot be used since theoretically they are unobservable). Note that it is a simple regression, one of the component is here only to bring some noise. First of all, let us look at marginal probit regression

>  reg1=glm(Y1~X1+X2,data=base,family=binomial)
> reg2=glm(Y2~X1+X2,data=base,family=binomial)
> summary(reg1)
 
Call:
glm(formula = Y1 ~ X1 + X2, family = binomial, data = base)
 
Deviance Residuals:
Min 1Q Median 3Q Max
-2.90570 -0.50126 -0.00266 0.49162 2.78256
 
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.291725 0.267149 -16.065 <2e-16
X1 8.656836 0.510153 16.969 <2e-16 ***
X2 0.007375 0.090530 0.081 0.935
---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1386.29 on 999 degrees of freedom
Residual deviance: 726.48 on 997 degrees of freedom
AIC: 732.48

Number of Fisher Scoring iterations: 5
> summary(reg2)
Call:
glm(formula = Y2 ~ X1 + X2, family = binomial, data = base)
Deviance Residuals:
Min 1Q Median 3Q Max
-
2.74682 -0.51814 -0.00001 0.57969 2.58565
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(
Intercept) 3.91709 0.24399 16.054 <2e-16 ***
X1 -7.89703 0.46277 -17.065 <2e-16 ***
X2 0.18360 0.08758 2.096 0.036 *
---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1386.29 on 999 degrees of freedom
Residual deviance: 777.61 on 997 degrees of freedom
AIC: 783.61
Number of Fisher Scoring iterations: 5

Here, the optimization yields,

> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par
> OPT[1]
[1] 0.5261382
> matrix(OPT[2:7],2,3)
[,1] [,2] [,3]
[1,] -2.451721 4.908633 0.01600769
[2,] 2.241962 -4.539946 0.10614807

Note that the coefficients we have obtained are almost identical to the ones obtained with R standard procedure,

>  library(Zelig)
> REG= zelig(list(mu1=Y1~X1+X2,
+ mu2=Y2~X1+X2,
+ rho=~1),
+ model="bprobit",data=base)
> summary(REG)
 
Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
rho = ~1), model = "bprobit", data = base)
 
Pearson Residuals:
Min 1Q Median 3Q Max
probit(mu1) -10.5442 -0.377243 0.0041803 0.36709 8.60398
probit(mu2) -7.8547 -0.376888 0.0083715 0.42923 5.88264
rhobit(rho) -13.8322 -0.091502 -0.0080544 0.37218 0.85101
 
Coefficients:
Value Std. Error t value
(Intercept):1 -2.451699 0.135369 -18.11116
(Intercept):2 2.241964 0.125072 17.92536
(Intercept):3 1.169461 0.189771 6.16249
X1:1 4.908617 0.252683 19.42602
X1:2 -4.539951 0.233632 -19.43203
X2:1 0.015992 0.050443 0.31703
X2:2 0.106154 0.049092 2.16235
 
Number of linear predictors: 3
 
Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
 
Dispersion Parameter for binom2.rho family: 1
 
Residual Deviance: 1460.355 on 2993 degrees of freedom
 
Log-likelihood: -730.1774 on 2993 degrees of freedom
 
Number of Iterations: 3

> matrix(coefficients(REG)[c(1:2,4:7)],2,3)
[,1] [,2] [,3]
[1,] -2.451699 4.908617 0.01599183
[2,] 2.241964 -4.539951 0.10615443
The correlation here is also the same
> (exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1)
[1] 0.5260951

That procedure works well an can be extended to ordinal responses (not only binary ones, or to three dimensional problems,

logV=function(beta){
BETA=matrix(beta[4:(3+ncol(Y)*ncol(X))],ncol(Y),ncol(X))
z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]),X%*%(BETA[3,]))
r12=beta[1]
r23=beta[2]
r31=beta[3]
s1=s2=s3=1
sigma=matrix(c(s1^2,r12*s1*s2,r31*s1*s3,
r12*s1*s2,s2^2,r23*s2*s3,
r31*s1*s3,r23*s2*s3,s3^2),3,3)
sigma1=matrix(c(s2^2,r23*s2*s3,
r23*s2*s3,s3^2),2,2)
sigma2=matrix(c(s1^2,r31*s1*s3,
r31*s1*s3,s3^2),2,2)
sigma3=matrix(c(s1^2,r12*s1*s2,
r12*s1*s2,s2^2),2,2)
a111=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a111=c(a111,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a011=pmnorm(z[1,2:3],varcov=sigma1)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a011=c(a011,pmnorm(z[i,2:3],varcov=sigma1)-pmnorm(z[i,],varcov=sigma))}
a101=pmnorm(z[1,c(1,3)],varcov=sigma2)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a101=c(a101,pmnorm(z[i,c(1,3)],varcov=sigma2)-pmnorm(z[i,],varcov=sigma))}
a110=pmnorm(z[1,1:2],varcov=sigma3)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a110=c(a110,pmnorm(z[i,1:2],varcov=sigma3)-pmnorm(z[i,],varcov=sigma))}
a100=pnorm(z[1,1],sd=s1)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a100=c(a100,pnorm(z[i,1],sd=s1)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a010=pnorm(z[1,2],sd=s2)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(2,3)],varcov=sigma1)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a010=c(a010,pnorm(z[i,2],sd=s2)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(2,3)],varcov=sigma1)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a001=pnorm(z[1,3],sd=s3)-pmnorm(z[1,c(2,3)],varcov=sigma1)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a001=c(a001,pnorm(z[i,3],sd=s3)-pmnorm(z[i,c(2,3)],varcov=sigma1)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a000=1-a111-a011-a101-a110-a001-a010-a100
 
a111[a111<=0]=1e-50
a110[a110<=0]=1e-50
a101[a101<=0]=1e-50
a011[a011<=0]=1e-50
a100[a100<=0]=1e-50
a010[a010<=0]=1e-50
a001[a001<=0]=1e-50
a000[a000<=0]=1e-50
 
-sum(((Y[,1]==0)&(Y[,2]==0)&(Y[,3]==0))*log(a111) +
((Y[,1]==1)&(Y[,2]==0)&(Y[,3]==0))*log(a011) +
((Y[,1]==0)&(Y[,2]==1)&(Y[,3]==0))*log(a101) +
((Y[,1]==0)&(Y[,2]==0)&(Y[,3]==1))*log(a110) +
((Y[,1]==1)&(Y[,2]==1)&(Y[,3]==0))*log(a001) +
((Y[,1]==1)&(Y[,2]==0)&(Y[,3]==1))*log(a010) +
((Y[,1]==0)&(Y[,2]==1)&(Y[,3]==1))*log(a100) +
((Y[,1]==1)&(Y[,2]==1)&(Y[,3]==1))*log(a000) )
}

A strong assumption in that bivariate model is that residuals have a Gaussian structure. It is possible to change that assumption

  • marginally: for instance if we use a logistic cumulative distribution function, then we will have a bivariate logit regression
  • in terms of dependence structure: it is possible to consider another copula than the gaussian one, e.g. Gumbel's copula (also called the bivariate logistic copula), or Clayton's
Here, the following code can be used to extend the model to non Gaussian structures,
> F=function(x,r){pmnorm(x,rep(0,length(x)),
+ varcov=matrix(c(1,r,r,1),2,2))}
> Fx=function(x1){F(c(x1,1e40),0)}
> Fy=function(x2){Fx(x2)}
>
> logVgen=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+ a11=F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a11=c(a11,F(z[i,],r=CORRELATION))}
+ a10=Fx(z[1,1])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a10=c(a10,Fx(z[i,1])-F(z[i,],r=CORRELATION))}
+ a01=Fy(z[1,2])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a01=c(a01,Fy(z[i,2])-F(z[i,],r=CORRELATION))}
+ a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+ ((Y[,1]==0)&(Y[,2]==1))*log(a01) +
+ ((Y[,1]==1)&(Y[,2]==0))*log(a10) +
+ ((Y[,1]==0)&(Y[,2]==0))*log(a00) )
+ }
>
> beta0=c(0,1,1,1,1,1,1)
> (OPT=optim(fn=logVgen,par=beta0,method="BFGS")$par)
[1] 0.52613820 -2.45172059 2.24196154 4.90863292 -4.53994592 0.01600769
[7] 0.10614807
There were 23 warnings (use warnings() to see them)
E.g.
> library(copula)
> F=function(x,r){pcopula(pnorm(x),
claytonCopula(2, r))}
> Fx=function(x1){F(c(x1,1e40),0)}
> Fy=function(x2){Fx(x2)}
  • An application to school tests

Consider the following dataset,

hsb2=read.table("http://freakonometrics.free.fr/hsb2.csv",
header=TRUE, sep=",")
math_male=hsb2$math[female==0]
write_male=hsb2$write[female==0]
math_female=hsb2$math[female==1]
write_female=hsb2$write[female==1]
plot(math_female, write_female, type="p",
pch=19,col="red",xlab="maths",ylab="writing",cex=.8)
points(math_male, write_male, cex=1.2, col="blue")

with here maths versus writing, with girls in red and boys in blue, where variables here are
  female :
0: male
1: female
race :
1: hispanic
2: asian
3: african-amer
4: white
ses :
1: low
2: middle
3: high
schtyp : type of school
1: public
2: private
prog : type of program
1: general
2: academic
3: vocation
read : reading score
write : writing score
math : math score
science : science score
socst : social studies score
We can try to understand correlation between math and writing skills. Covariates can be the sex of the child, and his reading skills. The question will then be: are good students in maths and writing simply students that can read well ?

Here the code is simply
> W=hsb2$write>=50
> M=hsb2$math>=50
> base=data.frame(Y1=W,Y2=M,
+ X1=hsb2$female,X2=hsb2$read)
>
> library(Zelig)
> REG= zelig(list(mu1=Y1~X1+X2,
+ mu2=Y2~X1+X2,
+ rho=~1),
+ model="bprobit",data=base)
> summary(REG)
 
Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
rho = ~1), model = "bprobit", data = base)
 
Pearson Residuals:
Min 1Q Median 3Q Max
probit(mu1) -4.7518 -0.502594 0.15038 0.53038 1.8592
probit(mu2) -3.4243 -0.653537 0.23673 0.67011 2.6072
rhobit(rho) -4.9821 0.010481 0.13500 0.40776 2.9171
 
Coefficients:
Value Std. Error t value
(Intercept):1 -5.484711 0.787101 -6.96825
(Intercept):2 -4.061384 0.633781 -6.40818
(Intercept):3 1.332187 0.322175 4.13497
X1:1 1.125924 0.233550 4.82092
X1:2 0.167258 0.202498 0.82598
X2:1 0.103997 0.014662 7.09286
X2:2 0.082739 0.012026 6.88017
 
Number of linear predictors: 3
 
Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
 
Dispersion Parameter for binom2.rho family: 1
 
Residual Deviance: 364.51 on 593 degrees of freedom
 
Log-likelihood: -182.255 on 593 degrees of freedom
 
Number of Iterations: 3
> (exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1)
[1] 0.5824045
with a remaining correlation among residuals of 0.58. So with only the sex of the student, and his or her reading skill, we cannot explain the correlation between maths and writing skills. With our previous code, we have here
> beta0=c((exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1),
+ summary(REG)@coef3[c(1:2,4:7),1])
> beta0
(Intercept):1 (Intercept):2 X1:1 X1:2
0.58240446 -5.48471133 -4.06138412 1.12592427 0.16725842
X2:1 X2:2
0.10399668 0.08273879
> (OPT=optim(fn=logV,par=beta0,method="BFGS")$par)
(Intercept):1 (Intercept):2 X1:1 X1:2
0.5824045 -5.4847113 -4.0613841 1.1259243 0.1672584
X2:1 X2:2
0.1039967 0.0827388

i.e. we obtain (almost) exactly the same estimators. But here I have used as starting values for the optimization procedure the estimators given by R. If we change them, hopefully we have a robust maximum likelihood estimator,

> (OPT=optim(fn=logV,par=beta0/2,method="BFGS")$par)
(Intercept):1 (Intercept):2 X1:1 X1:2
0.58233360 -5.49428984 -4.06839571 1.12696594 0.16760347
X2:1 X2:2
0.10417767 0.08287409
There were 12 warnings (use warnings() to see them)
So once again, it is possible to optimize numerically a likelihood function, and it works.

Tuesday, December 15 2009

Des quantiles aux régressions quantiles

Un petit billet rapide pour insister (davantage) sur l'importance des quantiles, et des applications possibles des régressions quantiles.

  • De la médiane aux quantiles
J'avais fait un petit billet (ici) pour me moquer des journalistes qui ne savent pas ce qu'est une médiane... mais j'ai parfois l'impression qu'ils ne sont pas les seuls à comprendre l'importance de ces grandeurs en modélisation. Les quantiles sont fondamentaux en tant qu'outils de mesure de risque (ici ou ). Ils sont aussi très présents dès lors que l'on s'intéresse aux richesses, et aux inégalités (ici par exemple). De même que l'on passe de la moyenne à l'espérance conditionnelle, on peut essayer de passer d'un quantile à un quantile conditionnel, et c'est l'intérêt de la régression quantile...
  • La régression quantile sur données climatiques
Vendredi dernier, James Elsner m'a présenté une jolie application des régressions quantiles sur les données climatiques (qu'il avait publié dans Nature). On peut récupérer les données sur sa page ouebe.
> StormMax=read.table("http://garnet.fsu.edu/~jelsner/extspace/extremedatasince1899.csv",header=TRUE,sep=",")
> StormMaxBasin=subset(StormMax,Region=="Basin"); StormMaxBasin=subset(StormMaxBasin,Yr>1977)
> attach(StormMaxBasin)
> boxplot(Wmax~as.factor(Yr),ylim=c(35,175),
+ xlab="Year",ylab="Intensity (kt)",col="light blue")
On peut regarder rapidement les données,

et éventuellement regarder l'allure d'une  régression linéaire (pour détecter une éventeullement tendance - linéaire - temporellle).

On peut aussi regarder des boxplot par annéeé de survenance,
> boxplot(Wmax~as.factor(Yr),ylim=c(35,175),
+ xlab="Year",ylab="Intensity (kt)",col="light blue")

On peut ainsi récupérer par années les quantiles élevés,
> x=boxplot(Wmax~as.factor(Yr),plot=F)
et tenter une régression sur les quantiles à 90%
>  xx=1:29
>  points(xx,x$stats[5,],pch=19,col="purple")
>  abline(lm(x$stats[5,]~xx),col="purple")

ou sur les quartiles à 75%
> points(xx,x$stats[4,],pch=19,col="purple")
> abline(lm(x$stats[4,]~xx),col="purple")

Bref, il semble que la tendance soit beaucoup plus marquée sur les évènements extrêmes que sur le reste... Et ça peut se vérifier en faisant plusieurs régressions quantiles, en regardant la significativité de la pente pour différents niveaux de quantiles.
> model=rq(Wmax~Yr,tau=seq(0.2,0.8,0.1))
> plot(summary(model,alpha=0.05,se="iid"),
+ parm=2,pch=19,cex=1.2,mar=c(5,5,4,2)+0.1,
+ ylab="Trend (kt/yr)",xlab="Quantile")

On note alors que la pente d'une régression quantile est strictement positive pour des quantiles élevés, supérieurs à 75%.

On peut d'ailleurs affiner un peu en prenant une séquence plus fine,

Pour la régression quantile locale, on peut utiliser la library(MASS)
> plot(Yr,Wmax,ylim=c(35,175),
+ xlab="Year",ylab="Intensity (kt)",col="blue")
> library(MASS)
> fit <- lprq(Yr,Wmax, h = 2, tau = 0.9)
> lines(fit$xx, fit$fv, col="red",lwd=2)
> fit <- lprq(Yr,Wmax, h = 4, tau = 0.9)
> lines(fit$xx, fit$fv, col="red")

ou encore la library(splines),
> plot(Yr,Wmax,ylim=c(35,175),
+ xlab="Year",ylab="Intensity (kt)",col="blue")
> library(splines)
> fit <- rq(Wmax~ bs(Yr, df = 6), tau = .8)
> X <- model.matrix(Wmax~ bs(Yr, df = 6))
> lines(Yr,X %*% fit$coef,col="red",lwd=2)
> fit <- rq(Wmax~ bs(Yr, df = 12), tau = .8)
> X <- model.matrix(Wmax~ bs(Yr, df = 12))
> lines(Yr,X %*% fit$coef,col="red")

Mais on peut faire la même chose, ou presque sur des vitesse de véhicules.... mais ça sera pour un autre billet dès que Ben aura fini de jouer avec les données.

Saturday, December 5 2009

De la surdispersion des nombres

J'avais déjà fait un billet (ici) pour rappeler que la loi de Poisson était très attirante, mais souvent peu adapté sur des données trop hétérogènes. En faisant des régressions, on essaye d'expliquer l'hétérogénéité avec des combinaisons linéaires de variables explicatives, autant que faire se peut. Et parfois, ce n'est pas suffisant, il reste de l'hétérogénéité résiduelle.

  • Effets fixes et effets aléatoires en économétrie
Quand on cherche é expliquer une variable (dite endogène, Y) par d'autres (dites exogènes), il y a celles dont on peut facilement disposer et qui sont observables (notées X), et celles qui sont cachées, et dont l'effet n'en est pas moins réel (notées Z). 
Formellement, on peut définir une vraisemblance sur l'ensemble de variables exogènes
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-01.png
et celle définie seulement sur les variables observables
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-02.png
Notons que
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-03.png
En actuariat, on s'intéresse à des fréquences de sinistres, supposées suivre des lois de Poisson (du moins si l'on arrive é constituer des classes suffisamment homogènes). Suppsons que
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-04.png
alors
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-05.png
où 
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-06.png
On a alors un modéle dit à effets fixes, au sens où
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-07.png
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-08.png
La vraisemblance se calcule alors par intégration sur les effets aléatoires. La variance de cet effet aléatoire va induire de l'hétérogénéité. Le cas le plus classique est obtenu lorsque l'on suppose que U suit une loi Gamma, auquel cas la loi de Y conditionnelle à X est une loi binomiale négative.
Si l'on suppose que http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-09.png, alors http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-10.png et
 http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-11.png
Alors, en posant http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-12.png on en déduit que
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-13.png
et donc
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-14.png
En revanche, si on calcule la loi nonconditionnelle, alors
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-15.png
alors que
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-16.png
soit
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-17.png
  • Tester la surdispersion
Un test graphique peut être facilement construit, en représentant la variance empirique dans un sous groupe de mêmes caractéristiques par rapport à l'espérance, par sous-classe. Dans un modèle d'équidispersion, les données seraient alignés sur une droite de pente 1.
> BASE=read.table("D:\ -data\\base.txt",header=TRUE,sep=";")
> I=is.na(BASE$NBSIN_DV_08)
> X1=BASE$ZONIER_08[I==FALSE]
> Y=BASE$NBSIN__08[I==FALSE]
> reglmqp=glm(Y~as.factor(X1)+0,family=quasipoisson)
> reglmp =glm(Y~as.factor(X1)+0,family=poisson)
> summary(reglmqp)
Call:
glm(formula = Y ~ as.factor(X1) + 0, family = quasipoisson)
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
as.factor(X1)1   -3.74950    0.73230  -5.120 3.06e-07 ***
as.factor(X1)2   -3.51155    0.73230  -4.795 1.63e-06 ***
as.factor(X1)3   -3.56311    0.31225 -11.411  < 2e-16 ***
[...]
as.factor(X1)27  -3.43071    0.03534 -97.091  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 1.072523)
On note qu'en ajustant une loi quasiPoisson, R nous dit qu'il y a un paramètre de surdispersion de l'ordre de 1,07 (correspond à un effet de surdispersion). Il s'agit de la courbe bleue sur le graphique ci-dessus. La courbe verte correspond à la régression (par moindres carrés) des variances empiriques sur les moyennes empiriques (sans constante,  avec une pente de l'ordre de 1,033).
> M = tapply(Y, as.factor(X1), mean)
> V = tapply(Y, as.factor(X1), sd)^2
> plot(M,V")
> abline(0,1,col="red")
> abline(0, 1.076690  ,col="blue")
> abline(lm(V~0+M),col="green")

Pour un test plus formel, il convient de spécifier davantage l'hypothèse alternative. On supposera par exemple que
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-18.png
qui correspond tout simplement à un modèle avec effet aléatoire, et http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-19.png est simplement la variance de l'effet aléatoire.
On cherche alors à tester
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-20.png contre http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-21.png
On peut alors considérer la statistique de test suivante
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-22.png
Cette statistique doit suivre, sous H0 une loi normale centrée réduite. Parmil les autres statistiques classiques, on pourra aussi utiliser
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/disper-23.png
Sous R, ces tests ont été implémentés dans la library(AER), via la fonction dispersiontest().
> library(AER)
> dispersiontest(reglmp)
        Overdispersion test
data:  reglmp
z = 5.9049, p-value = 1.764e-09
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion
  1.076690
> dispersiontest(reglmp,trafo = 2)
        Overdispersion test
data:  reglmp
z = 6.0101, p-value = 9.27e-10
alternative hypothesis: true alpha is greater than 0
sample estimates:
   alpha
2.566681
Ces deux tests sont détaillés dans les deux livres Cameron et Trivedi (1998, 2005).

Un autre exemple pourrait être le "nombre d'heures, en moyenne, passées à regarder la télévision" (via une des enquêtes de l'insee, ici). On se doute que la "consommation" de télévision (pour poursuivre dans la syntaxe de Patrick Le Lay) est liée à l'âge, et a priori pas de manière linéaire. On peut donc tenter un modèle gam,

> p1=read.dbf("D:\ -data\\hdv1.dbf")
> p2=read.dbf("D:\ -data\\hdv2.dbf")
> p3=read.dbf("D:\ -data\\hdv3.dbf")
> p4=read.dbf("D:\ -data\\hdv4.dbf")
> base=data.frame(p1,p2,p3,p4)
> reg3=gam(LT1FREQ~s(AGEE),
+ family=quasipoisson)
> summary(reg3)
Call: gam(formula = LT1FREQ ~ s(AGEE),
+ family = quasipoisson, data = base)
(Dispersion Parameter for quasipoisson family taken to be 2.5857)
Ce qui confirme la forte surdispersion observée sur le graphique moyenne-variance (par âge pris comme facteur) ci-dessus.
  • Quid de la sous dispersion ?
Il existe des cas où les nombres sont sous dispersés. Regardons par exemple le nombre de redoublements de classes par des enfants scolarisés.
On se doute que l'âge de l'enfant intervient. On peut alors régresser suivant l'âge, ou bien régresser sur une constante, et prendre l'âge (ou plutôt le logarithme de l'âge si l'on considère une fonction de lien logarithmique) en variable offset. A partir de la base fournie par l'insee (ici), on obtient
>  library(foreign)
> parent1=read.dbf("D:\ -data\\parent1.dbf")
> parent2=read.dbf("D:\ -data\\parent2.dbf")
> parent3=read.dbf("D:\ -data\\parent3.dbf")
> base=data.frame(parent1,parent2,parent3)
> base0=data.frame(Y=base$NBREDOUB,
+ E=base$AGEA)
> reg1=glm(Y~E,data=base0,family=quasipoisson)
> summary(reg1)
Call:
glm(formula = Y ~ E, family = quasipoisson, data = base0)
Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -3.474161   0.089512  -38.81   <2e-16 ***
E            0.158290   0.005087   31.12   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.8086682)
ou dans le cas de la régression où l'âge est déclaré en variable offset
> reg2=glm(Y~1,offset=log(E),data=base0,family=quasipoisson)
> summary(reg2)
Call:
glm(formula = Y ~ 1, family = quasipoisson, data = base0, offset = log(E))
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -3.65357    0.02421  -150.9   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.7639523)
Bref, les deux régression nous confortent dans l'idée d'une sous-dispersion (mais les tests ne sont pas implémentés pour ça).
Si on pense que l'effet de l'âge n'est pas spécialement linéaire, on peut aussi tenter un modèle gam.

> library(gam)
> reg3=gam(Y~s(E),data=base0,family=quasipoisson)
> summary(reg3)
Call: gam(formula = Y ~ s(E), family = quasipoisson, data = base0)
(Dispersion Parameter for quasipoisson family taken to be 0.8095)
Ce qui confirme cet effet de sousdispersion. Et là encore, un rapide test graphique confirme (plus ou moins) cette intuition. La courbe rouge  correspond là encore à l'équidispersion, la courbe bleue est de pente 0,81 (obtenue avec deux régressions) et la courbe mauve de pente 0,76.
Un autre exemple pourrait être le nombre de mariages d'une personne. Là aussi, ce nombre doit dépendre de l'âge de la personne interrogée.
> p1=read.dbf("D:\ -data\\hdv1.dbf")
> p2=read.dbf("D:\ -data\\hdv2.dbf")
> p3=read.dbf("D:\ -data\\hdv3.dbf")
> p4=read.dbf("D:\ -data\\hdv4.dbf")
> base=data.frame(p1,p2,p3,p4)
> reg1=glm(BNBMAR~AGEE,
+ family=quasipoisson)
> summary(reg1)
Call:
glm(formula = BNBMAR ~ AGEE, family = quasipoisson, data = base)
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.1685735  0.0230460  -50.71   <2e-16 ***
AGEE         0.0188916  0.0004085   46.25   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.3370578)
> reg2=glm(BNBMAR~1,offset=log(AGEE),data=base,family=quasipoisson)
> summary(reg2)
Call:
glm(formula = BNBMAR ~ 1, family = quasipoisson, data = base,
    offset = log(AGEE))
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.081337   0.006818  -598.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasipoisson family taken to be 0.3178957)
  • Trop de valeurs nulles ?
Il est possible, dans certains car qu'un problème vienne d'un trop grand nombre de valeurs nulles pour accepter une hypothèse de loi de Poisson. La library(pscl) propose une fonction zeroinfl() assez simple d'utilisation. Au niveau de la synthaxe, après le | figurent les variables expliquant le surnombre des valeurs nulles. On peut tester par exemple ce problème sur les redoublements (ou l'on peut suspecter que le "problème" vient du fait que trop peu de personnes ne redoublent pas, ou de manière duale, qu'il existe une spirale du redoublement)
> base=data.frame(parent1,parent2,parent3)
> base0=data.frame(Y=base$NBREDOUB,E=base$AGEA)
> Y=base0$Y;X1=base0$E
> library(pscl)
> reg4 <- zeroinfl(Y ~ E | E,
+ data = base0, dist = "poisson")
> summary(reg4)
Call:
zeroinfl(formula = Y ~ E | E, data = base0, dist = "poisson")
Count model coefficients (poisson with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.18824    0.19465 -11.242   <2e-16 ***
E            0.09077    0.01040   8.727   <2e-16 ***
Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  6.02989    0.60910   9.900   <2e-16 ***
E           -0.65280    0.07936  -8.225   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of iterations in BFGS optimization: 25
Log-likelihood: -2466 on 4 Df
Rapidement, notons que tous les coefficients sont significatifs, ce qui nous donne envie d'accepter l'hypothèse que la sous dispersion serait expliquée par un trop fort nombre de valeurs nulles.

Friday, October 30 2009

Qu'est-ce que la modéliation (économétrique, ou statistique) (1)

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/lunette.PNGAu vu des questions qui m'ont été posées en cours d'économétrie, j'ai décidé de faire un petit billet afin d'expliciter ce que j'ai pu raconter dans les premiers slides de remise à niveau de probabilité et de statistiques, en Master 1, où j'essayais de faire le lien entre le modèle probabiliste et l'analyse statistique, à partir du schéma (forcément simpliste) suivant

http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.formel-modelisation-slides_m.jpg
Pour mieux comprendre ce qu'est la modélisation, inspirons nous du papier de Joshua Angrist et Victor Lavy, "using Maimonides' rule to estimate the effect the effect of class size on scholastic achivement" (publié que le QJE en 1999, en ligne ici ou pour des transparents).
  • Partie 1: la problématique
La question que l'on se pose est simplement "quel est l'impact de la taille d'une classe sur les résultats scolaires ? ".
Beaucoup de monde (je pense aux parents ou aux enseignants, pas forcément aux ministres) pense que cet impact est réel, et négatif, à savoir que plus la taille et grande, moins bon sont les résultats, et inversement. Empiriquement, certaines études ont toutefois trouvé un effet positif, car les meilleurs élèves sont parfois regroupés dans les classes les plus chargées.
Bref, on souhaite répondre à cette question en faisant une étude empirique. La taille de la classe est une donnée facile à utiliser et à quantifier (encore que... le nombre sera celui au début du printemps), mais les résultats désignent une notion plus floue. Heureusement, certains pays diposent d'un test national (même s'il est parfois biaisé car les crédits des écoles étant liés à ces tests, certaines écoles semblent reforcer le bachotage).
  • Partie 2: les données (et le contexte)
Les données mises à disposition par Joshua Angrist et Victor Lavy sont celles d'écoles publiques en Israël, où deux tests ont été mis en oeuvre: un test national en fin d'étude (en mathématiques et en lecture) pour les élèves de CM1 et CM2 (4th et 5th grade) en 1991, et un test national en 1992 pour les élèves de CE2 (3rd grade). Nous n'avons pas des données individuelles (élève par élève) mais par classe. Afin de corriger d'éventuels effets exogènes (car on se doute que les résutats ne sont pas uniquement liés à la taille des classes), les auteurs utilisent des mesures de l'environnement social de l'école.
Voilà un peu pour les données, qui se trouvent ici et là. Sur le contexte, une information importante: l'étude se fait sur les écoles publiques israéliennes, où l'on croit dur comme fer à un effet négatif. En effet, depuis 1969, ces écoles appliquent (strictement) une règle talmudique énoncée par le rabin Maionide (alias Moshe ben Maimon, הרב משה בן מיימון, ou أبو عمران موسى بن ميمون بن عبد الله القرطبي الإسرائيلي ou Moussa ibn Maimoun ibn Abdallah al-Kourtoubi al-Israili dans la version arabe, plutôt connu pour ses réflexions en médecine ou en santé publique) au 12ème siècle, à savoir "dès qu'il y a plus de vingt cinq élèves dans une classe, il faut un assistant; dès qu'il y a plus de quarante élèves, il faut deux professeurs". 
Une fois formulée la question, et une fois à notre disposition un jeu de données pour faire une étude, on peut mettre en oeuvre un test économétrique (ou statistique).
  • Partie 3: la formalisation du problème
On commencer par formaliser le problème. On note
http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-1.png
le nombre d'élèves dans la classe i au printemps (ce qui sera noté comme la taille de la classe). On retiendra une indexation http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-2.png par école (school) et par classe (class). On note aussi
http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-3.png
la note moyenne obtenue dans la classe lors des tests nationaux. On a aussi ces variables explicatives
http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-4.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-X1.png désigne le nombre d'élèves dans l'école pour chacun des niveaux, http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-X1.png un identifiant pour la ville, http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-X1.png une variable indiquant le niveau socio-économique des élèves, et http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-X1.png un identifiant éthnique (juif ou arabe) et religieux (associé aux écoles). On notera que, compte tenu que compte tenu de la règle de Maionide,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-5.png           (1)
où la fonction http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-10.png désigne la partie entière. 
  • Partie 4: l'étude économétrique, les moindres carrés
L'étude économétrique est une phase d'estimation des paramètres du modèle formel. Mais auparavant, il est toujours intéressant (car instructif) de faire un peu de statistiques descriptives, en regardant quelques graphiques.
On peut déjà se demander si la règle de Maionide est effectivement mise en oeuvre dans les établissements scolaires, c'est à dire regarder si l'équation (1) est valide,
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.ecole-stats-1_m.jpg
Sur ce graphique, la courbe en trait pointillé correspond à la loi théorique, c'est à dire la partie de droite dans l'équation (1), et le trait plein correspond à la version empirique, c'est à dire la partie de gauche. Dans un second temps, regardons sur le même type de graphique comment évolue la note moyenne dans la classe (avec toujours en pointillé le nombre théorique d'élèves dans la casse).
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.ecole-stats-2_m.jpg
On note ici que les résultats aux tests sont générallement plus élevé dans les classes avec plus d'élèves dans l'école. Il faut aussi noter que les très grosses écoles sont plutôt dans des grandes villes riches, alors que les petites écoles sont des "development town" beaucoup plus pauvres (et à l'extérieur des villes). 
Bref, on voudrait pouvoir comprendre un peu mieux ce qui se passe... Pour cela , on (ou plutôt les auteurs) spécifie un modèle économétrique sous la forme  suivante
http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-7.png
Autrement dit dans le premier terme http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-8.png on va mettre tout ce qui, dans les aspects socio-économiques - pourrait expliquer les résultats aux tests. Le terme du milieu http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-9.png correspond à ce qui nous intéresse, à savoir l'impact de la taille des classes les tests. Enfin, le terme de droite http://perso.univ-rennes1.fr/arthur.charpentier/latex/model-latex-10.png correspond à une moyenne (par classe, mais aussi par école, afin de corriger d'un éventuel effet "très bonne école") mais aussi un bruit car, par définition d'un modèle, la modélisation est forcément imparfaite: il y aura toujours des effets que l'on ne pourra pas expliquer avec un modèle simple (car on cherche un modèle simple, ou plutôt parcimonieux, comme je l'évoquai ici).
Encore une fois, un modèle économétrique mesure une corrélation (éventuelle) entre des variables explicatives et la variable que l'on essaye d'expliquer, mais l'interprétation causale est fallacieuse (je reviendrais sur ce point dans un billet prochainement).
Mais dans un premier temps, on peut s'intéresser à une interprétation par moindres carrés du modèle, voire une simplification du modèle, où l'on régresse uniquement les notes sur la taille de la classe. Bon, on a accès aux données (ici pour une version csv, ou pour l'ensemble des données).Contentons nous déjà d'un niveau (4e), et regardons respectivement en lecture et en mathématique l'impact de la taille de la classe. On commence par un petit nettoyage de la base de données,
> base5=read.table("http://perso.univ-rennes1.fr/arthur.charpentier/Final4.csv",sep=";",header=TRUE)
> base=base5
> base$avgverb=base$avgverb-(base$avgverb>100)*100
> base$avgmath=base$avgmath-(base$avgmath>100)*100
> base$func1= base$c_size/(trunc((base$c_size-1)/40)+1)
> base$func2= base$cohsize/(trunc(base$cohsize/40)+1)
> base$verbsize[base$avgverb==NA]=0
> base$verbsize[base$passverb==NA]=0
> base$mathsize[base$avgmath==NA]=0
> base$mathsize[base$passmath==NA]=0
> base0=base[(base$classize>1)&(base$classize<45)&(base$c_size>5)&(base$c_leom==1)&(base$c_pik<3),]

Enuiste, on peut tenter de faire une régression
> reg=lm(avgverb~classize,data=base0)
> summary(reg)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 67.76309    0.78686  86.118   <2e-16 ***
classize     0.22119    0.02568   8.614   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.55 on 2016 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared: 0.0355,     Adjusted R-squared: 0.03502
F-statistic: 74.21 on 1 and 2016 DF,  p-value: < 2.2e-16

On note effectivement que la taille de la classe a un effet sur le nombre d'élèves, et que cet effet est croissant. On peut visualiser ça graphiquement
> taille=seq(5,45,by=.1)
> note=predict(reg,newdata=data.frame(classize=taille),interval="confidence")
> plot(base0$classize,base0$avgverb,cex=.6)
> lines(taille,note[,2],lty=2,col="red")

On peut aussi faire une régression locale pour détecter d'éventuelles nonlinéairités,
> regL=loess(avgverb~classize,data=base0,span=.5)
> note=predict(regL,newdata=data.frame(classize=taille),interval="confidence")
> plot(base0$classize,base0$avgverb,cex=.6)
> lines(taille,note,col="blue")

Si on regarde pour la note en maths, on a un effet assez proche,
> reg=lm(avgmath~classize,data=base0)
> summary(reg)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 57.36364    0.98724   58.10   <2e-16 ***
classize     0.33055    0.03222   10.26   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 9.472 on 2016 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared: 0.04963,    Adjusted R-squared: 0.04916
F-statistic: 105.3 on 1 and 2016 DF,  p-value: < 2.2e-16

Mais on se doute que quelque chose ne va pas, et que ce que l'on mesure ici doit être lié à autre chose.... et il faudra aller plus loin (ce qui fera l'objet d'un autre billet).

Thursday, October 1 2009

De l'interprétation d'un effet nonlinéaire en régression

Je vais poursuivre un peu le TD de statistique de l'actuariat 2, où nous avions parlé tarification et régression Poissonnienne.

  • Analyser un effet nonlinéaire 
Dans un premier temps, nous avions regardé un modèle simple, où l'on voulait voir si l'âge du conducteur expliquait la fréquence de sinistres. Bref, on peut utiliser une régression log-Poisson, qui donne la prédiction suivante, en fonction de l'âge,

Mais l'âge est une variable un peu particulière car c'est une fausse variable continue: seule les valeurs entières apparaissent. Comme on chercher à calculer
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-1.png
notons qu'il est possible de considérer l'estimateur empirique naturel de cette grandeur, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-2.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-glm-nb-3.png
Graphiquement, on obtient la courbe suivante,

Le côté ératique de la série à droite et à à gauche vient de la faible représentation dans la base, comme en atteste l'histogramme de l'âge du conducteur,

Je peux en profiter pour revenir à un commentaire posté ici sur les modèles GAM. En effet, au lieu de regarder la fréquence en fonction de l'âge, on peut s'intéresser à la la différence par rapport à la fréquence globale.
> m=mean(base$nbsin)
> Y=predict(REG,newdata=data.frame(age=X,expo=E),type = "response")
> lines(X,Y/m-1,col="red")
On obtient alors le graphique suivant (en superposant à l'ajustement obtenu par GAM)

Notons que si l'on ajuste un modèle log-Poisson uniquement sur la tranche d'âge [50-70], on est très proche de ce que propose le modèle GAM, d'où l'interprétation de l'ajustement GAM comme un ajustement local.


  • Changer de modèle (et les comparer)
On a choisi un modèle log-Poisson, mais on se doute que d'autres modèles seraient possible. En particulier dans notre base, les assurés ont eu 0 ou 1 sinistre l'année passée. Autrement dit, des lois binomiales peuvent être ajustée,
> REG0 = glm(nbsin~age,base,
+  family=poisson(link="log"),
+  offset=expo)
> REG1a = glm(nbsin~age,base,
+  family=binomial(link = "logit"),
+  offset=expo)
> REG1b = glm(nbsin~age,base,
+  family=binomial(link = "probit"),
+  offset=expo)
> REG1c = glm(nbsin~age,base,
+  family=gaussian(link = "identity"),
+  offset=expo)

Le critère AIC est souvent invoqué pour juger de l'adéquation d'un modèle. Il vaut ici
> AIC(REG0)
[1] 36470.7
> AIC(REG1a)
[1] 36044.45
> AIC(REG1b)
[1] 37304.02
> AIC(REG1c)
[1] 32549.83

Mais on peut également regarder le critère BIC,
> AIC(REG0, k = log(nrow(base)))
[1] 36488.81
> AIC(REG1a, k = log(nrow(base)))
[1] 36062.57
> AIC(REG1b, k = log(nrow(base)))
[1] 37322.14
> AIC(REG1c, k = log(nrow(base)))
[1] 32577
Autrement dit, en utilisant le règle "the smaller, the better" on retiendrait le dernier modèle qui est un modèle Gaussien, manifestement assez mauvais: sur le graphique ci-dessous, la régression log-Poisson est en mauve, presque confondu avec la régression logistique1, en bleu, la régression probit en vert et le modèle linéaire (Gaussien) en rouge,   

Bref, ce dernier modèle n'est pas convainquant... Par contre, notons que si l'on compare le modèle log-Poisson et logit, ce dernier conduit à tarifer moins cher les âges extrêmes,

Autrement dit, au lieu de comparer des modèles actuariels sur la base d'outils statistiques, et uniquement de ces outils, il peut être instructif d'essayer de comprendre qui verra sa prime baisser avec tel ou tel modèle... Mais je reviendrais plus longuement, une autre fois, sur les histoires de choix de modèles.....

Propriétés des estimateurs dans une régression

Comme les rappels (je devrais plutôt dire "remise à niveau") ont été particulièrement rapides, je vais prendre un peu de temps pour revenir sur quelques éléments du cours d'économétrie (d'autant plus que j'ai reçu quelques questions par mail).

  • sur le calcul des estimateurs
Je vais reprendre les questions que j'ai reçu, ça sera plus simple, "je comprend pas très bien quelle méthode utilise l'opération lm(Y~X.....) sur R pour trouver les coef" ou encore "je regresse un modèle simple, je sors le summary, et j'obtiens un betaREG 3 étoiles. Toujours sur R, je compare ce betaREG avec le betaCalc = (X'X)-1 * X'Y, il est différent ! Comment expliquer cela ? Est-ce possible ?".
Commençons par la première question: la fonction lm() utilise la méthode dite des moindres carrés, ce qui revient à calculer
en adoptant une écriture matricielle, i.e. la solution est alors
Cet estimateur est celui qui minimise la somme des carrés des erreurs. Et on peut vérifier numériquement que cet estimateur est bien celui calculé par R,
>data(cars)
>X=cbind(rep(1,nrow(cars)),cars$speed); Y=cars$dist
> X[1:5,]
      [,1] [,2]
 [1,]    1    4
 [2,]    1    4
 [3,]    1    7
 [4,]    1    7
 [5,]    1    8
> solve(t(X)%*%X)
            [,1]         [,2]
[1,]  0.19310949 -0.011240876
[2,] -0.01124088  0.000729927
> t(X)%*%Y
      [,1]
[1,]  2149
[2,] 38482
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
           [,1]
[1,] -17.579095
[2,]   3.932409

ce qui correspond très précisément à ce que calcule la fonction lm(),
> lm(dist~speed,cars)
Coefficients:
(Intercept)        speed
    -17.579        3.932
Et on peut aussi vérifier qu'il minimse bien la somme des carrés des erreurs,
> b0=seq(-30,10,by=3);
> b1=seq(1,7,by=.5)
> SC=matrix(NA,length(b0),length(b1))
> for(i in 1:length(b0)){
> for(j in 1:length(b1)){
 + SC[i,j]=sum((cars$dist-(b0[i]+b1[j]*cars$speed))^2)
}}
> contour(b0,b1,SC)
et je peux certifier que la valeur obtenue est bien le minimum de cette fonction (en fait c'est un résultat théorique d'optimisation).
Ca se généralise bien entendu en dimension plus grande. Alors attention, cette formule permet de calculer un estimateur du vecteur des paramètres, mais la formule est bien sûre fausse si on raisonne composante par composante.
> X=cars$speed; Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
         [,1]
[1,] 2.909132

Mais ça c'est un résultat (trivial) d'algèbre linéaire. En effet, le fait d'avoir

ne signifie pas que l'on ait ce genre de formule composante par composante, i.e.
Dit de manière un peu formelle, le produit matriciel n'est pas un produit terme à terme. La traduction économétrique de cette idée est qu'une régression multiple n'est pas une succession de régressions simples. En effet, dans l'équation précédante, le dernier terme correspond à la régression sans la constante, comme on peut le voir ci-dessous
> lm(dist~0+speed,cars)
Coefficients:
speed 
2.909 
et le premier est tout simplement la moyenne empirique,
> lm(dist~1,cars)
Coefficients:
(Intercept) 
      42.98 
> mean(cars$dist)
[1] 42.98
> X=rep(1,nrow(cars)); Y=cars$dist
> (solve(t(X)%*%X))%*%(t(X)%*%Y)
      [,1]
[1,] 42.98
  • sur les propriétés des estimateurs
Une question était "doit-on vérifier les propriétés de l'estimateur (sans biais, efficace etc...) du summary de la Reg ? pour chaque régression ?"
Bon, la réponse est non car c'est impossible ! Pour répéter tranquillement ce que j'avais dit plusieurs fois oralement: quand on fait un modèle (théorique), on part d'hypothèses (par exemple ici les résidus sont centrées, de variance constante). Sous ces hypothèses on a des propriétés: ce sont des théorèmes, autrement dit de la théorie. Le plus connu étant le théorème de Gauss-Markov sur le modèle linéaire, qui garantie que l'estimateur par moindre carré est BLUE (best linear unbiased estimator).
Autrement dit (si je retraduit ce que dit ce théorème), si les hypothèses sont valides, alors la théorie nous garantie que les estimateurs vérifient des propriétés, dont celle d'être sans biais, par exemple.
On ne peut pas vérifier qu'un estimateur est sans biais, on peut juste vérifier ex post que les hypothèses sont valides (ou non), ce qui garantit (ou pas) l'absence de biais.
Mais comme on va le voir par la suite, le biais ce n'est pas forcément génant, ce qui est important, c'est surtout la convergence. Comme le notait Clive Granger  "if you can't even get a consistant estimator, you shouldn't be in this business".
  • et sur des points plus théoriques
sinon la personne qui avait des soucis numériques à calculer les coefficients me demandait si cela avait un "rapport avec les plim X'E/n et plim (X'X)-1/ n".
Pour revenir sur ce point, une hypothèse forte dans le modèle linéaire (de base) est que le bruit soit soit un vrai bruit, c'est à dire non corrélé avec la variable explicative, i.e.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-1.png
En effet, sinon
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-2.png
autrement dit l'estimateur par moindre carré est biasé. La solution la plus naturelle pour éviter ce "problème" est d'utiliser un instrument (ou une variable instrumentale).
Mais surtout, sur cet exemple, on voit qu'on ne pourra avoir la convergence de notre estimateur par moindres carrés que si la corrélation entre le bruit et la variable explicative tend vers 0 asymptotiquement.
Il existe en effet une différence fondamentale entre le biais et la consistance (ou la convergence vers la vraie valeur).
  • un estimateur sans biais est généralement convergent (et de manière consistante)
  • un estimateur biaisé peut ne pas converger
Prenons un petit exemple pour illustrer ces points. On supposera n1 < n2 < n3par la suite. Comme je l'ai expliqué auparavant, la distribution de l'estimateur est une propriété théorique qui ne peut se voir sur un jeu de données (à moins de faire des simulations mais ça sort du cadre de mon billet d'aujourd'hui). La figure ci-dessous correspond au cas sans biais, et convergent
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-sans-biais-convergent_m.jpg
La figure ci-dessous correspond au cas biaisé, mais convergent, et consistant, au sens où asymptotiquement, l'estimateur sera sans biais.
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-biais-convergent_m.jpg
Enfin, un dernier cas correspondant au cas biaisé, convergent, mais non consistant. Asymptotiquement l'estimateur sera toujours biaisé,
http://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.graph-estimateur-biais-non-convergent_m.jpg
Mais revenons un peu sur la formaliation sous-jacente au modèle linéaire. Les propriétés sur les résidus sont conditionnelles aux variables explicatives, en particulier.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-3.png
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-4.png
(c'est l'hypothèse d'homoscédasticité et d'indépendance des résidus). Bon, il faut aussi l'absence de relation linéaire entre les variables explicatives, ce qui se traduit parfois comme une absence de collinéarité.
On peut s'intéresser aux propriétés asymptotiques de l'estimateur. Pour cela, il faut peut être rappeler ce que signifie la convergence pour une variable aléatoire (en statistique mathématique, les estimateurs sont vus comme des variables aléatoires). On dira que http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-5.png converge en probabilité vers http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-6.png, parfois noté
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-7.png
si
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-8.png
Par exemple la loi des grands nombres garantie que, pour un échantillon i.i.d.
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-9.png
 si l'espérance des variables est finie, alors
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-10.png
Alors sous les hypothèses mentionnées auparavant, on peut montrer que, dans le modèle linéaire,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-11.png
On parlera alors d'estimateur convergent. On peut aussi montrer que
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-12.png
Mais formellement,
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-13.png
si
http://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-14.png
Voilà en gros quelques éléments pour répondre aux questions, mais tout ça sera repris en détails dans le cours d'économétrie 1 qui commencera très bientôt... et surtout d'économétrie 2 qui abordera les aspects dynamiques.

Monday, June 15 2009

Séries temporelles versus séries individuelles

J'ai reçu un mail l'autre jour que me demandait "pourquoi ne peut-on pas utiliser directement les méthodes d'économétrie sur des séries chronologiques ?". Au lieu de fournir une réponse théorique ou formelle, je préférerais reprendre l'exemple que j'utilisais en cours. Pour cela, je vais utiliser une série que j'utilise dans mon cours d'économétrie.
Supposons que l'on s'interroge sur une éventuellement corrélation entre le tour de poitrine et le tour de taille chez les personnes de sexe féminin. Pour faire une étude économétrique, on a à notre disposition les mensurations des demoiselles qui posent nues, tous les mois, dans Playboy, depuis 1953. Autrement dit, on fait l'erreur de prendre ces données temporelles (tous les mois depuis 1953) pour une série d'observations individuelles. Regardons dans notre échantillon (presque 660 demoiselles) le tour de poitrine  et le tour de taille.Sur l'ensemble de la période considérée, les deux mesures semblent statistiquement indépendantes (ou alors avec une corrélation très faible). Ce qui est douteux  si on y réfléchit deux minutes. Et effectivement, si on fait l'étude par sous-période, on a effectivement une dépendance positive entre les deux grandeurs.

Et c'est quoi l'interprétation ? En fait, l'erreur dans l'analyse du premier graph vient du fait que la morphologie desdites demoiselles a beaucoup changé dans le temps. Autrement dit, la dépendance est sûrement beaucoup plus forte que ce que laisse sous-entendre le graphique ci-dessus. Regardons par exemple l'évolution du ratio entre le tour de taille et le tour de poitrine, On note que les deux variables évoluent en sens inverse. Plus précisément, le tour de poitrine a très fortement diminué, alors que le tour de taille a plutôt augmenté

En fait, plus généralement, beaucoup d'autres paramètres ont changé dans le temps, par exemple l'âge (ce qui peut poser la question de savoir s'il y a réellement un changement dans la corpulence dans le temps, ou c'est simplement du au fait que les demoiselles sont plus agées),    Moralité ? Il ne faut pas jeter ses vieux Playboy, on peut toujours faire des études économétriques avec....  ou alors il faut faire attention quand on étudie des variables à voir s'il n'y a pas une variable cachée derrière (ici le temps).

Sinon, une autre représentation graphique est possible, en interpolant (par une parabole) les mensurations de poitrine, tailles et hanches (et en tenant compte de la taille de la playmate). On a alors l'animation suivante

http://freakonometrics.blog.free.fr/public/perso3/animation-pb2.gif

(avec en haut la taille, au milieu la taille, et en bas les hanches, afin d'avoir une représentation de la morphologie moyenne des playmates cette année là).

Wednesday, March 18 2009

Constante dans une régression

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