econometrics - Freakonometrics

Miscellaneous › econometrics

Saturday, March 17 2012

De la normalité dans les régressions

By arthur charpentier on Saturday, March 17 2012, 00:01

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

(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

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.

no trackback

Thursday, February 23 2012

Visualization in regression analysis

By arthur charpentier on Thursday, February 23 2012, 10:52

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

no trackback

Thursday, February 2 2012

On linear models with no constant and R2

By arthur charpentier on Thursday, February 2 2012, 09:31

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

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

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

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

The least-square program is here

i.e. the estimate of the slope is

Numerically, we obtain

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

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

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

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

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

Recall that the total "coefficient of variation", denoted , is defined as

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 . So obviously, using the 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 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 might look good, but meaningless !

4 comments no trackback

Wednesday, May 11 2011

Multivariate probit regression using (direct) maximum likelihood estimators

By arthur charpentier on Wednesday, May 11 2011, 13:46

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

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.

2 comments no trackback

Tuesday, December 15 2009

Des quantiles aux régressions quantiles

By arthur charpentier on Tuesday, December 15 2009, 15:42

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 là). 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.

2 comments no trackback

Saturday, December 5 2009

De la surdispersion des nombres

By arthur charpentier on Saturday, December 5 2009, 23:08

overdispersion

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

où

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

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 H₀ 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,

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.

no trackback

Friday, October 30 2009

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

By arthur charpentier on Friday, October 30 2009, 17:25

Au 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 là 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$

où $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 là 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).

no trackback

Thursday, October 1 2009

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

By arthur charpentier on Thursday, October 1 2009, 11:10

regression

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$

où

$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 logistique¹, 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.....

no trackback

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

By arthur charpentier on Thursday, October 1 2009, 11:09

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 n₁ < n₂ < n₃par 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$

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

$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.

no trackback

Monday, June 15 2009

Séries temporelles versus séries individuelles

By arthur charpentier on Monday, June 15 2009, 21:06

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

(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à).

no trackback

Wednesday, March 18 2009

Constante dans une régression

By arthur charpentier on Wednesday, March 18 2009, 22:38

constante
R2

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 R² 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 R² de 0,4 (comme une corrélation de 0,4 d'ailleurs) ne veut rien dire en soit.

Freakonometrics

Miscellaneous › econometrics

De la normalité dans les régressions

Visualization in regression analysis

On linear models with no constant and R2

Multivariate probit regression using (direct) maximum likelihood estimators

Des quantiles aux régressions quantiles

De la surdispersion des nombres

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

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

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

Séries temporelles versus séries individuelles

Constante dans une régression

Presentation

Keep in touch

Pages

Categories

Plan du site

Search