Category Archives: Regression

Régression, variables explicatives et géométrie

La régression (comme tout calcul d’espérance conditionnelle) est un problème de projections

Un théorème intéressant est le théorème dit de Frisch-Waugh, permettant de comprendre la différence fondamentale entre un modèle de régression multiple, et les modèles de régression simple. On reviendra sur ce point en cours en évoquant rapidement le paradoxe de Simpson. La formulation est la suivante: on veut régresser http://freakonometrics.hypotheses.org/files/2016/05/FW15.gif sur http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW17.gif, deux ensembles (a priori disjoints) de variables explicatives,

http://freakonometrics.hypotheses.org/files/2016/05/FW02.gif

Les équations normales (associées au problème de minimisation de la somme des carrés des erreurs) sont

http://freakonometrics.hypotheses.org/files/2016/05/FW11.gif

de telle sorte qu’à l’optimum on peut relier les estimateurs des deux jeux de paramètres par

http://freakonometrics.hypotheses.org/files/2016/05/FW12.gif

La première partie correspond à la régression de http://freakonometrics.hypotheses.org/files/2016/05/FW15.gif sur http://freakonometrics.hypotheses.org/files/2016/05/FW13.gif, mais il reste un second terme dès lors que http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW17.gif ne sont pas orthogonales. Notons http://freakonometrics.hypotheses.org/files/2016/05/FW04.gifla matrice de projection (orthogonale) sur http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif

http://freakonometrics.hypotheses.org/files/2016/05/FW06.gif

Alors on peut écrire

http://freakonometrics.hypotheses.org/files/2016/05/FW03.gif

En posant

http://freakonometrics.hypotheses.org/files/2016/05/FW07.gif et http://freakonometrics.hypotheses.org/files/2016/05/FW08.gif,

on retrouve un modèle linéaire classique (à condition de travailler sur la projection des variables sur le sous-espace engendré par http://freakonometrics.hypotheses.org/files/2016/05/FW16.gif , i.e. en transformant les variables),

http://freakonometrics.hypotheses.org/files/2016/05/FW10.gif

Pour aller plus loin sur la géométrie des moindres carrées, et sur le théorème de Frisch-Waugh, je peux renvoyer à des notes de cours, et à quelques transparents

Le graphique ci-dessous correspond au cas où les variables explicatives sont orthogonales, et dans ce cas, la régression multiple est équivalente à deux régressions simples

http://freakonometrics.hypotheses.org/files/2016/05/FW1.gif

Le graphique ci-dessous au cas non orthogonal,

http://freakonometrics.hypotheses.org/files/2016/05/FW2.gif

Le code qui permet de vérifier ces histoires de projections successives est relativement simple. Tout d’abord on importe les données, et on regarde le modèle global,

> chicago=read.table(
+ "http://freakonometrics.free.fr/chicago.txt",
+ header=TRUE,sep=";")
> Y=chicago$Fire
> X1=chicago$X_1
> X2=chicago$X_2
> X3=chicago$X_3
> base=data.frame(Y,X1,X2,X3)
> tail(base)
Y    X1 X2     X3
42  4.8 0.152 19 13.323
43 10.4 0.408 25 12.960
44 15.6 0.578 28 11.260
45  7.0 0.114  3 10.080
46  7.1 0.492 23 11.428
47  4.9 0.466 27 13.731
> regression=lm(Y~X1+X2+X3)
> summary(regression)

Call:
lm(formula = Y ~ X1 + X2 + X3)

Residuals:
Min     1Q Median     3Q    Max
-9.737 -4.565 -1.479  3.751 16.079

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.07525    6.19447   3.564 0.000910 ***
X1          -0.62764    5.28130  -0.119 0.905953
X2           0.22378    0.06161   3.632 0.000744 ***
X3          -1.55059    0.38195  -4.060 0.000204 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.527 on 43 degrees of freedom
Multiple R-squared: 0.4417,	Adjusted R-squared: 0.4027
F-statistic: 11.34 on 3 and 43 DF,  p-value: 1.314e-05

> n=length(Y)
> X=matrix(c(rep(1,n),X1,X2,X3),n,4)
> X[1:5,]
[,1]  [,2] [,3]   [,4]
[1,]    1 0.604   29 11.744
[2,]    1 0.765   44  9.323
[3,]    1 0.735   36  9.948
[4,]    1 0.669   37 10.656
[5,]    1 0.814   53  9.730

Ensuite, on va projeter seulement sur les deux premières variables (et la constante)

> FWX1=X[,1:3]
> regression12=lm(Y~X1+X2)
> solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y
[,1]
[1,]  0.08069764
[2,] 11.56913900
[3,]  0.15108490
> summary(regression12)$coefficients
Estimate Std. Error    t value   Pr(>|t|)
(Intercept)  0.08069764 3.49182154 0.02311047 0.98166664
X1          11.56913900 5.05014782 2.29085156 0.02681947
X2           0.15108490 0.06854408 2.20420043 0.03278652
>
>
> FWX2=X[,4]
> H1=FWX1%*%solve(t(FWX1)%*%FWX1)%*%t(FWX1)
> M1=diag(rep(1,n))-H1
> FWX2s=M1%*%FWX2
> FWYs =M1%*%Y
> (beta2=solve(t(FWX2s)%*%FWX2s)%*%t(FWX2s)%*%FWYs)
[,1]
[1,] -1.550594
> summary(regression)$coefficients[4]
[1] -1.550594
>
> (beta1=solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y-
+        solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%FWX2%*%beta2)
[,1]
[1,] 22.0752495
[2,] -0.6276442
[3,]  0.2237765
> summary(regression)$coefficients[1:3]
[1] 22.0752495 -0.6276442  0.2237765

On retrouve bien l’estimateur du dernier paramètre à l’aide du théorème de Frisch-Waugh, puis en utilisant les équations normales, on en déduit les trois premiers.

Normalité des estimateurs dans une régression

Cette semaine en cours, on a évoqué la normalité des estimateurs dans une régression linéaire. En supposant la normalité des résidus, on a vu en cours que http://freakonometrics.blog.free.fr/public/perso5/ols-002.gif était un estimateur Gaussien. En particulier, chacun des estimateurs est alors Gaussien, au sens où http://freakonometrics.blog.free.fr/public/perso5/ols--0003.gif, ce qui peut se visualiser sur le graphique suivant (la constante est en abscisse, et la pente en ordonnée), avec un intervalle de confiance à 95%,

En fait, la variance étant inconnue (mais pouvant être estimée), si on remplace la variance des résidus par leur estimateur, la loi de l’estimateur est une loi de Student. De même pour l’estimateur de la pente de la droite de régression, http://freakonometrics.blog.free.fr/public/perso5/ols---ooo4.gif

Le fait que le couple soit également Gaussien fait que l’on peut construire non plus un intervalle de confiance (on n’est plus en dimension 1) mais une ellipse de confiance. On retient une forme elliptique car c’est la plus petite région dans laquelle on se trouvera avec une probabilité de 95% (comme discuté dans un vieux billet).

Mais pour mieux comprendre cette notion d’ellipse de confiance, le plus simple est d’avoir recours à du rééchantillonnage dans la base de données. En effet, comme on ne dispose que d’un jeu de données, on a un estimateur. Et la discussion sur la loi de notre estimateur est purement théorique. Pour visualiser la loi de notre estimateur, il faudrait des centaines, voire des milliers de bases similaires. Que l’on n’a pas. La solution est alors de tirer au hasard des points de notre échantillon (avec remise), i.e. de faire du bootstrap,

set.seed(1)
COEF=matrix(NA,10000,2)
for(s in 1:nrow(COEF)){
I=sample(1:nrow(cars),nrow(cars),replace=TRUE)
COEF[s,]=lm(dist~speed,data=cars[I,])$coefficients
}

Si on regarde ce qui se passe, tirage après tirage, on génère un paquet d’échantillon, et pour chaque échantillon, on ajuste une droite de régression.

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

Les distributions – sur tous ces échantillons – de la constante et de la pente semblent effectivement Gaussiennes,

hist(COEF[,1],col="light blue",prob=TRUE)
u=seq(min(COEF[,1]),max(COEF[,1]),length=500)
v=dnorm(u,mean(COEF[,1]),sd(COEF[,1]))
lines(u,v,lwd=3,col="red")
hist(COEF[,2],col="light blue",prob=TRUE)
u=seq(min(COEF[,2]),max(COEF[,2]),length=500)
v=dnorm(u,mean(COEF[,2]),sd(COEF[,2]))
lines(u,v,lwd=3,col="red")

pour la distribution de l’estimateur de la pente, alors que pour l’estimateur de la constante, on a la distribution suivante,

Si on regarde maintenant la loi jointe,

on retrouve une forme elliptique pour le nuage des points (et une forte corrélation négative entre les deux estimateurs). Mais si on creuse un peu

library(ellipse)
reg=lm(dist~speed,data=cars)
e=ellipse(reg)
plot(e,type="l",lwd=2)
polygon(e,col="light blue")
points(COEF,cex=.5)

l’ellipse ne coïncide pas (tout à fait) avec celle obtenu théoriquement. Ce qui laisse à penser que l’hypothèse de normalité des résidus est peut-être à revoir..,

Infidelity and econometrics

On http://www.bakadesuyo.com, there was recently an interesting discussion about infidelity, the key question being “at what ages are men and women most likely to have affairs?” The discussion is based on some graphs, e.g.

The source is a paper by Donald Cox. Based on a sample of 36 men and 22 women 3,432 respondent (NHSLS dataset) . And to be honest, I have been surprised by the shape of the curves. Especially for men… In order to compare, it is possible to use another dataset that can be found in R,

> library(Ecdat)
> data(Fair)
> tail(Fair)
sex age   ym child religious education occupation
596   male  47 15.0   yes         3        16          4
597   male  22  1.5   yes         1        12          2
598 female  32 10.0   yes         2        18          5
599   male  32 10.0   yes         2        17          6
600   male  22  7.0   yes         3        18          6
601 female  32 15.0   yes         3        14          1
rate nbaffairs
596    2         7
597    5         1
598    4         7
599    5         2
600    2         2
601    5         1

with 601 observations (from Fair (1977)). It is possible to run a Poisson regression to describe the number of affairs in the past year. E.g for men

> library(splines)
> regM=glm(nbaffairs~bs(age),family=poisson,
+ data=Fair[Fair$sex=="male",])
> a=seq(20,60)
> N=predict(regM,newdata=data.frame(age=a),type="response")
> plot(a,N,type="l",lwd=2,col="red")

or for women,

> regF=glm(nbaffairs~bs(age),family=poisson,
+ data=Fair[Fair$sex=="female",])
> N=predict(regF,newdata=data.frame(age=a),type="response")
> plot(a,N,type="l",lwd=2,col="red",lty=2)

On that (larger) dataset, we obtain curves that are more intuitive… But maybe the Poisson distribution is not an appropriate model. For instance, having no affairs do not mean that the person did not want to… So perhaps, a more interesting model would be a Poisson model with a zero-inflation, i.e. some people are honest and do not want to have affairs (and appear as 0), while some do want to have some affairs, and the number of affairs is Poisson distributed (and can take the value 0). If we focus on people wo do not want to have affairs, the model (and the prediction) is the following, where we plot the probability of not being interested in having an affair,

> library(pscl)
> regM0=zeroinfl(nbaffairs~bs(age)|bs(age),family=poisson,
+ link="logit",data=Fair[Fair$sex=="male",])
> N0=predict(regM0,newdata=data.frame(age=a),type="zero")
> plot(a,N0,type="l",lwd=2,col="blue")

For those willing to have an affair, here is the parameter of the Poisson distribution of the number of affairs,

> Nc=predict(regM0,newdata=data.frame(age=a),type="count")
> plot(a,Nc,type="l",lwd=2,col="purple")

The same can be done for women, with the probability of no-willing to have an affair,

and to Poisson rate for women willing to have an affair,

If we focus on people willing to have an affair, the curves are the following,

i.e. men below 40 have more interested, but after 40, the probability drops, while women are still more and more likely to be willing to have an affair. On the other hand, young women having affairs might be less, but they usually have much more affairs than men…

Sélection de variables versus sélection de modalités

En cours, nous avions évoqué (très rapidement) la sélection automatique de variables. La méthode la plus simple est une méthode stepwise, basé sur un critère de type AIC, ou BIC. Considérons la base suivante,

>  N = base$nbre
>  E = base$exposition
>  X1 = base$carburant
>  X2 = cut(base$agevehicule,c(0,3,10,101),
+ right=FALSE)
>  X3 = cut(base$ageconducteur,c(0,22,45,101),
+ right=FALSE)
>  X4 = as.factor(base$zone)
>  X5 = as.factor(base$puissance)
>  X6 = as.factor(base$region)
>  X7 = as.factor(base$marque)
>  base1=data.frame(N,E,X1,X2,X3,X4,X5,X6,X7)

Une méthode stepwise (backward) donne ici

> reg1=glm(N~X1+X2+X3+X4+X5+X6+X7+offset(log(E)),
+ family="poisson",data=base1)
> step(reg1)
Start:  AIC=20492.67
N ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + offset(log(E))

Df Deviance   AIC
- X5   11    15316 20482
- X3    2    15305 20490
<none>       15304 20493
- X2    2    15314 20499
- X1    1    15319 20506
- X7   10    15343 20511
- X4    5    15398 20576
- X6   14    15569 20729

Step:  AIC=20482.35
N ~ X1 + X2 + X3 + X4 + X6 + X7 + offset(log(E))

Df Deviance   AIC
- X3    2    15317 20479
<none>       15316 20482
- X2    2    15326 20488
- X1    1    15334 20498
- X7   10    15359 20505
- X4    5    15410 20566
- X6   14    15579 20717

Step:  AIC=20479.33
N ~ X1 + X2 + X4 + X6 + X7 + offset(log(E))

Df Deviance   AIC
<none>       15317 20479
- X2    2    15327 20485
- X1    1    15334 20495
- X7   10    15360 20502
- X4    5    15410 20563
- X6   14    15620 20754

Call:  glm(formula = N ~ X1 + X2 + X4 + X6 + X7 
       + offset(log(E)),
       family = "poisson",
data = base1)

Coefficients:
(Intercept)          X1E     X2[3,10)   X2[10,101)          X4B
-1.0588454   -0.1653822    0.0266763   -0.1135451   -0.0004047
X4C          X4D          X4E          X4F          X60
0.1497622    0.3748811    0.5052894    0.4292016   -0.3590838
X61          X62          X63          X64          X65
-0.9300641   -1.0278887   -1.1818218   -1.0971797   -0.9459414
X66          X67          X68          X69         X610
-1.3690795   -1.1425678   -1.5309402   -1.3883549   -1.4603624
X611         X612         X613          X72          X73
-1.6763206   -1.3974092   -1.4864404    0.0246113    0.1144990
X74          X75          X76         X710         X711
-0.0932555    0.1635397   -0.1478095    0.2502030    0.1967970
X712         X713         X714
-0.2420215    0.2161411   -0.1963162

Degrees of Freedom: 49999 Total (i.e. Null);  49967 Residual
Null Deviance:	    15810
Residual Deviance: 15320 	AIC: 20480

Autrement dit, on supprime la troisième (âge du conducteur principal, par classes arbitraires) et la cinquième variable (puissance du véhicule) en gardant toutes les autres. Mais ici, si une variable n’a pas été retenue, c’est que globalement, elle n’apportait pas beaucoup d’information. Il serait toutefois possible de garder une information partielle, en gardant éventuellement certaines modalités. L’idée est de disjoncter la base, en créant des variables indicatrices par modalités. La base sera beaucoup plus grosse, et la sélection prendra alors beaucoup plus de temps,

> base2=data.frame(model.matrix( ~ 0+X1+X2+X3+X4+X5+X6+X7,
+ data=base1))
> base2$E=base1$E
> base2$N=base1$N
> reg2=glm(N~.-E+offset(log(E)),family="poisson",
+ data=base2)
>  step(reg2)
Start:  AIC=20492.67
N ~ (X1D + X1E + X2.3.10. + X2.10.101. + X3.22.45. + X3.45.101.
X4B + X4C + X4D + X4E + X4F + X55 + X56 + X57 + X58 + X59 +
X510 + X511 + X512 + X513 + X514 + X515 + X60 + X61 + X62 +
X63 + X64 + X65 + X66 + X67 + X68 + X69 + X610 + X611 + X612 +
X613 + X72 + X73 + X74 + X75 + X76 + X710 + X711 + X712 +
X713 + X714 + E) - E + offset(log(E))

Step:  AIC=20492.67
N ~ X1D + X2.3.10. + X2.10.101. + X3.22.45. + X3.45.101. + X4B
X4C + X4D + X4E + X4F + X55 + X56 + X57 + X58 + X59 + X510 +
X511 + X512 + X513 + X514 + X515 + X60 + X61 + X62 + X63 +
X64 + X65 + X66 + X67 + X68 + X69 + X610 + X611 + X612 +
X613 + X72 + X73 + X74 + X75 + X76 + X710 + X711 + X712 +
X713 + X714 + offset(log(E))

Df Deviance   AIC
- X4B         1    15304 20491
- X58         1    15304 20491
- X511        1    15304 20491
- X2.3.10.    1    15304 20491
- X72         1    15304 20491
- X513        1    15304 20491
- X512        1    15304 20491
- X515        1    15304 20491
- X74         1    15305 20491
- X3.45.101.  1    15305 20491
- X714        1    15305 20491
- X55         1    15305 20492
- X3.22.45.   1    15305 20492
- X711        1    15306 20492
- X76         1    15306 20492
- X59         1    15306 20492
<none>             15304 20493
- X514        1    15306 20493
- X713        1    15306 20493
- X73         1    15307 20493
- X56         1    15307 20493
- X710        1    15307 20494
- X75         1    15308 20494
- X2.10.101.  1    15308 20495
- X57         1    15309 20495
- X4C         1    15310 20496
- X510        1    15310 20496
- X60         1    15312 20498
- X4F         1    15314 20500
- X712        1    15316 20503
- X1D         1    15319 20506
- X4D         1    15337 20524
- X61         1    15345 20532
- X65         1    15350 20536
- X62         1    15352 20538
- X64         1    15359 20545
- X4E         1    15362 20549
- X63         1    15366 20553
- X67         1    15370 20556
- X612        1    15381 20568
- X69         1    15382 20569
- X66         1    15387 20574
- X610        1    15389 20576
- X68         1    15393 20580
- X611        1    15406 20592
- X613        1    15451 20637

Step:  AIC=20490.67
N ~ X1D + X2.3.10. + X2.10.101. + X3.22.45. + X3.45.101. + X4C
X4D + X4E + X4F + X55 + X56 + X57 + X58 + X59 + X510 + X511 +
X512 + X513 + X514 + X515 + X60 + X61 + X62 + X63 + X64 +
X65 + X66 + X67 + X68 + X69 + X610 + X611 + X612 + X613 +
X72 + X73 + X74 + X75 + X76 + X710 + X711 + X712 + X713 +
X714 + offset(log(E))

etc etc… et si on va directement à la fin,

Step:  AIC=20469.18
N ~ X1D + X2.10.101. + X4C + X4D + X4E + X4F + X57 + X510 + X60
X61 + X62 + X63 + X64 + X65 + X66 + X67 + X68 + X69 + X610 +
X611 + X612 + X613 + X73 + X75 + X76 + X710 + X712 + X713 +
offset(log(E))

Df Deviance   AIC
<none>             15315 20469
- X76         1    15317 20470
- X713        1    15317 20470
- X73         1    15317 20470
- X57         1    15318 20470
- X75         1    15318 20471
- X710        1    15319 20471
- X510        1    15319 20471
- X4C         1    15322 20474
- X60         1    15322 20475
- X2.10.101.  1    15325 20478
- X4F         1    15325 20478
- X1D         1    15333 20485
- X712        1    15338 20490
- X61         1    15356 20508
- X4D         1    15359 20511
- X62         1    15363 20515
- X65         1    15363 20515
- X64         1    15371 20524
- X63         1    15378 20530
- X67         1    15383 20536
- X4E         1    15390 20543
- X612        1    15394 20547
- X69         1    15396 20548
- X66         1    15400 20553
- X610        1    15403 20555
- X68         1    15407 20559
- X611        1    15419 20572
- X613        1    15467 20619

Call:  glm(formula = N ~ X1D + X2.10.101. + X4C + X4D + X4E + X4F
X57 + X510 + X60 + X61 + X62 + X63 + X64 + X65 + X66 + X67 +
X68 + X69 + X610 + X611 + X612 + X613 + X73 + X75 + X76 +
X710 + X712 + X713 + offset(log(E)), family = "poisson",
data = base2)

Coefficients:
(Intercept)          X1D   X2.10.101.          X4C          X4D
-1.20880      0.16886     -0.13808      0.14888      0.37539
X4E          X4F          X57         X510          X60
0.50458      0.42768      0.08381      0.18722     -0.36509
X61          X62          X63          X64          X65
-0.93836     -1.03471     -1.18803     -1.10217     -0.95624
X66          X67          X68          X69         X610
-1.37463     -1.15391     -1.54213     -1.40188     -1.47217
X611         X612         X613          X73          X75
-1.68559     -1.40582     -1.49700      0.10874      0.15022
X76         X710         X712         X713
-0.15183      0.21948     -0.27400      0.19565

Degrees of Freedom: 49999 Total (i.e. Null);  49971 Residual
Null Deviance:	    15810
Residual Deviance: 15310 	AIC: 20470

Si la troisième variable (âge du conducteur principal, par classes arbitraires) disparait assez vite, en revanche, une information sur la cinquième (la puissance) est gardée car certaines modalités semblent être informative sur la fréquence d’accidents. En revanche, on notera qui si on fait un arbre, la troisième variable était toujours clairement significative, ce qui peut nous conforter dans l’idée de faire de la sélection de variables sur les modalités.

> library(tree)
> TREE= tree(N~X1+X2+X3+X4+X5+X6+X7+offset(log(E)),split="gini",
+ mincut = 2500,data=base1)
> plot(TREE)
> text(TREE,cex=.9)

Confidence interval for predictions with GLMs

Consider a (simple) Poisson regression http://freakonometrics.hypotheses.org/files/2016/11/poiss01.gif. Given a sample http://freakonometrics.hypotheses.org/files/2016/11/poiss02.gif where http://freakonometrics.hypotheses.org/files/2016/11/poiss03.gif, the goal is to derive a 95% confidence interval for http://freakonometrics.hypotheses.org/files/2016/11/poiss04.gif given http://freakonometrics.hypotheses.org/files/2016/11/poiss05.gif, where http://freakonometrics.hypotheses.org/files/2016/11/poiss04.gif is the prediction. Hence, we want to derive a confidence interval for the prediction, not the potential observation, i.e. the dot on the graph below

> r=glm(dist~speed,data=cars,family=poisson)
> P=predict(r,type="response",
+ newdata=data.frame(speed=seq(-1,35,by=.2)))
> plot(cars,xlim=c(0,31),ylim=c(0,170))
> abline(v=30,lty=2)
> lines(seq(-1,35,by=.2),P,lwd=2,col="red")
> P0=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> points(30,P1$fit,pch=4,lwd=3)

i.e.

Let http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif denote the maximum likelihood estimator of http://freakonometrics.hypotheses.org/files/2016/11/poiss07.gif. Then
http://freakonometrics.hypotheses.org/files/2016/11/poiss40.gif
where http://freakonometrics.hypotheses.org/files/2016/11/poiss101.gif is Fisher information of http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif (from standard maximum likelihood theory). Recall that
http://freakonometrics.hypotheses.org/files/2016/11/poiss13.gif
where computation of those values is based on the following calculations
http://freakonometrics.blog.fre<br /><br /> e.fr/public/latex/poiss21.gif
In the case of the log-Poisson regression
http://freakonometrics.hypotheses.org/files/2016/11/poiss36.gif
Let us get back to our initial problem.

  • confidence interval for the linear combination

A first idea to get a confidence interval for http://freakonometrics.hypotheses.org/files/2016/11/poiss49.gif is to get a confidence interval for http://freakonometrics.hypotheses.org/files/2016/11/poiss100.gif (by taking exponential values of bounds, since the exponential is a monotone function). Asymptotically, we know that
http://freakonometrics.hypotheses.org/files/2016/11/poiss40.gif

thus, an approximation for the variance matrix of http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif will be based on http://freakonometrics.hypotheses.org/files/2016/11/poiss45.gif, obtained by plugging estimators of the parameters.
Then, since http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif as an asymptotic multivariate distribution, any linear combination of the parameters will also be normal, i.e.
http://freakonometrics.hypotheses.org/files/2016/11/poiss47.gif has a normal distribution, centered on http://freakonometrics.hypotheses.org/files/2016/11/poiss49.gif, with variance http://freakonometrics.hypotheses.org/files/2016/11/poiss102.gif where http://freakonometrics.hypotheses.org/files/2016/11/Poiss110.gif is the variance of http://freakonometrics.hypotheses.org/files/2016/11/poiss06.gif. All those quantities can be easily computed. First, we can get the variance of the estimators

> i1=sum(predict(reg,type="response"))
> i2=sum(cars$speed*predict(reg,type="response"))
> i3=sum(cars$speed^2*predict(reg,type="response"))
> I=matrix(c(i1,i2,i2,i3),2,2)
> V=solve(I)

Hence, if we compare with the output of the regression,

> summary(reg)$cov.unscaled
(Intercept)         speed
(Intercept)  0.0066870446 -3.474479e-04
speed       -0.0003474479  1.940302e-05
> V
[,1]          [,2]
[1,]  0.0066871228 -3.474515e-04
[2,] -0.0003474515  1.940318e-05

Based on those values, it is easy to derive the standard deviation for the linear combination,

> x=30
> P2=predict(r,type="link",se.fit=TRUE,
+ newdata=data.frame(speed=x))
> P2
$fit
1
5.046034

$se.fit
[1] 0.05747075

$residual.scale
[1] 1

> sqrt(V[1,1]+2*x*V[2,1]+x^2*V[2,2])
[1] 0.05747084
> sqrt(t(c(1,x))%*%V%*%c(1,x))
[,1]
[1,] 0.05747084

And once we have the standard deviation, and normality (at least asymptotically), confidence intervals are derived, and then, taking the exponential of the bounds, we get confidence interval

> segments(30,exp(P2$fit-1.96*P2$se.fit),
+ 30,exp(P2$fit+1.96*P2$se.fit),col="blue",lwd=3)

Based on that technique, confidence intervals are no longer centered on the prediction. But who cares ?

  • delta method

Actually, those who like to use “more or less” expressions for confidence intervals will not like non centered intervals. So, an alternative is to use the delta method. Instead of writing (again) something on the theory, we can use a package which computes that method,

> estmean=t(c(1,x))%*%coef(reg)
> var=t(c(1,x))%*%summary(reg)$cov.unscaled%*%c(1,x)
> library(msm)
> deltamethod (~ exp(x1), estmean, var)
[1] 8.931232
> P1=predict(r,type="response",se.fit=TRUE,
+ newdata=data.frame(speed=30))
> P1
$fit
1
155.4048

$se.fit
1
8.931232

$residual.scale
[1] 1

The delta method gives us (asymptotic) normality, so once we have a standard deviation, we get the confidence interval.

> segments(30,P1$fit-1.96*P1$se.fit,30,
+ P1$fit+1.96*P1$se.fit,col="blue",lwd=3)

Note that those quantities – obtained with two different approaches – are rather close here

> exp(P2$fit-1.96*P2$se.fit)
1
138.8495
> P1$fit-1.96*P1$se.fit
1
137.8996
> exp(P2$fit+1.96*P2$se.fit)
1
173.9341
> P1$fit+1.96*P1$se.fit
1
172.9101
  • bootstrap techniques

And a third method (but far from what I expect to teach on that course) is to use bootstrap techniques to about those results based on asymptotic normality (we have only 50 observations). The idea is to sample from out dataset, and to run a log-Poisson regression on those new samples, and to repeat a lot of time,

Régression sur des variables catégorielles

Petit complément par rapport au cours de mardi. On avait évoqué tout d’abord la lecture des sorties lorsque l’on régresse sur des variables catégorielles (des facteurs). Commençons par supprimer la constante de la régression

> reg0=glm(nbre~0+zone,offset=log(exposition),data=base, 
+ family=poisson(link="log"))
> summary(reg0)

Call:
glm(formula = nbre ~ 0 + zone, family = poisson(link = "log"), 
    data = base, offset = log(exposition))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5717  -0.3968  -0.2996  -0.1547  12.6722  

Coefficients:
      Estimate Std. Error z value Pr(>|z|)    
zoneB -2.54187    0.06287  -40.43   <2e-16 ***
zoneA -2.54912    0.05285  -48.23   <2e-16 ***
zoneC -2.38525    0.03753  -63.56   <2e-16 ***
zoneD -2.13454    0.03878  -55.05   <2e-16 ***
zoneE -2.00204    0.03965  -50.49   <2e-16 ***
zoneF -2.06932    0.11547  -17.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 50966  on 50000  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(reg0,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
        1         2         3         4         5 
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On voit que toutes les modalités sont présentes, et toutes sont significatives. Si on régresse sur la constante, il faudra supprimer une modalité pour rendre le modèle identifiable. On peut forcer pour que la modalité de référence soit la seconde,

> base$zone=relevel(base$zone,"B")
> regB=glm(nbre~zone,offset=log(exposition),data=base,
+ family=poisson(link="log"))
> summary(regB)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54187    0.06287 -40.431  < 2e-16 ***
zoneA       -0.00725    0.08213  -0.088 0.929661
zoneC        0.15662    0.07322   2.139 0.032432 *
zoneD        0.40733    0.07387   5.514 3.50e-08 ***
zoneE        0.53983    0.07433   7.263 3.80e-13 ***
zoneF        0.47255    0.13148   3.594 0.000325 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

> predict(regB,newdata=data.frame(
+ zone=c("A","B","C","D","E"),exposition=rep(1,5)))
1         2         3         4         5
-2.549120 -2.541870 -2.385253 -2.134543 -2.002044

On notera que les prédictions ne changent pas. On peut aussi choisir la première comme modalité de référence,

> base$zone=relevel(base$zone,"A")
> reg=glm(nbre~zone,offset=log(exposition),
> data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zone, family = poisson(link = "log"), 
data = base,
offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3968  -0.2996  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54912    0.05285 -48.232  < 2e-16 ***
zoneB        0.00725    0.08213   0.088 0.929661
zoneC        0.16387    0.06482   2.528 0.011471 *
zoneD        0.41458    0.06555   6.324 2.54e-10 ***
zoneE        0.54708    0.06607   8.280  < 2e-16 ***
zoneF        0.47980    0.12699   3.778 0.000158 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49994  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Le fait que la seconde modalité ne soit pas significative se lit par rapport à la modalité de référence (en l’occurrence la première): non significatif signifie alors non significativement différente. Autrement dit, on peut regrouper les modalités en une seule.

> base$zonesimple=base$zone
> base$zonesimple[base$zone%in%c("A","B")]="A"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.54612    0.04046 -62.937  < 2e-16 ***
zonesimpleC  0.16087    0.05518   2.915  0.00355 **
zonesimpleD  0.41158    0.05604   7.345 2.06e-13 ***
zonesimpleE  0.54408    0.05665   9.605  < 2e-16 ***
zonesimpleF  0.47681    0.12235   3.897 9.74e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

On note qu’avec ce regroupement, les autres modalités sont sensiblement différentes. On peut aussi retenir la troisième comme modalité de référence

> base$zonesimple=relevel(base$zonesimple,"C")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.38525    0.03753 -63.557  < 2e-16 ***
zonesimpleA -0.16087    0.05518  -2.915  0.00355 **
zonesimpleD  0.25071    0.05396   4.646 3.39e-06 ***
zonesimpleE  0.38321    0.05460   7.019 2.24e-12 ***
zonesimpleF  0.31593    0.12142   2.602  0.00927 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Comme toutes les modalités semblent significatives, on peut tenter de prendre comme modalité de référence une des dernières (dont les estimations des coefficients donnent des résultats très proches)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5717  -0.3959  -0.2989  -0.1547  12.6722

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.06932    0.11547 -17.921  < 2e-16 ***
zonesimpleC -0.31593    0.12142  -2.602  0.00927 **
zonesimpleA -0.47681    0.12235  -3.897 9.74e-05 ***
zonesimpleD -0.06522    0.12181  -0.535  0.59232
zonesimpleE  0.06727    0.12209   0.551  0.58161
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15692  on 49995  degrees of freedom
AIC: 20798

Number of Fisher Scoring iterations: 6

Au vue de cette dernière sortie, on peut tenter de fusionner toutes les dernières classes ensembles

> base$zonesimple[base$zone%in%c("D","E","F")]="F"
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> summary(reg)

Call:
glm(formula = nbre ~ zonesimple, family = poisson(link = "log"),
data = base, offset = log(exposition))

Deviance Residuals:
Min       1Q   Median       3Q      Max
-0.5660  -0.3959  -0.3004  -0.1547  12.5929

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.07182    0.02696 -76.853  < 2e-16 ***
zonesimpleC -0.31344    0.04621  -6.783 1.18e-11 ***
zonesimpleA -0.47431    0.04861  -9.757  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 15809  on 49999  degrees of freedom
Residual deviance: 15698  on 49997  degrees of freedom
AIC: 20800

Number of Fisher Scoring iterations: 6

Bon, formellement, regrouper deux modalités (i.e. décréter que deux variables sont non significatives simultanément) demande un peu plus qu’un test de Student, ou que deux tests de Student…. Si on remonte un peu en arrière, on aurait pu faire un test multiple avant de fusionner les trois modalités (un test de type Fisher, ou une ANOVA)

> base$zonesimple=relevel(base$zonesimple,"F")
> reg=glm(nbre~zonesimple,offset=log(exposition),
+ data=base,family=poisson(link="log"))
> library(car)
> linearHypothesis(reg,c("zonesimpleD=0","zonesimpleE=0"))
Linear hypothesis test

Hypothesis:
zonesimpleD = 0
zonesimpleE = 0

Model 1: restricted model
Model 2: nbre ~ zonesimple

Res.Df Df  Chisq Pr(>Chisq)
1  49997
2  49995  2 5.7073    0.05763 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Manifestement, on peut accepter l’hypothèse que ces trois catégories n’en font qu’une. La zone géographique peut alors se découper en trois grandes zones, et pas en six. On notera que cela correspond à ce que propose un arbre de régression

> library(tree)
> arbre=tree(nbre~zone+offset(log(exposition)),
+ data=base,split="gini")
> plot(arbre)
> text(arbre)

ACT2040, introduction aux modèles linéaires généralisés

On commencera ce mardi les GLM, après avoir introduit les lois exponentielles (qui ont du être revues en démonstration vendredi dernier). La notation utilisée sera que la loi (densité ou fonction de probabilité) de http://freakonometrics.blog.free.fr/public/perso4/Yi-ltx.gif sera de la forme

http://freakonometrics.blog.free.fr/public/perso4/loi-exponentielle.gif

Pour un complément plus exhaustif, je renvoie au chapitre en ligne.

  • Le modèle linéaire (Gaussien)

Le modèle de base est le modèle Gaussien que l’on avait revu au dernier cours,

> X=c(1,2,3,4)
> Y=c(1,2,5,6)
> base=data.frame(X,Y)
> reg1=lm(Y~1+X,data=base)
> nbase=data.frame(X=seq(0,5,by=.1))
> Y1=predict(reg1,newdata=nbase)

Pour une prédiction (unique), on obtient la prédiction suivante

Le code pour une telle représentation est le suivant

> plot(X,Y,pch=3,cex=1.5,lwd=2,xlab="",ylab="")
> lines(nbase$X,Y1,col="red",lwd=2)
> u=2
> mu=predict(reg1)[2]
> sigma=summary(reg1)$sigma
> y=seq(0,7,.05)
> loi=dnorm(y,mu,sigma)
> segments(u,y,loi+u,y,col="light green")
> lines(loi+u,y)
> abline(v=u,lty=2)
> points(X[2],Y[2],pch=3,cex=1.5,lwd=2)
> points(X[2],predict(reg1)[2],pch=19,col="red")
> arrows(u-.2,qnorm(.05,mu,sigma),
+ u-.2,qnorm(.95,mu,sigma),length=0.1,code=3,col="blue")

On peut multiplier les prédictions, en se basant sur l’hypothèse d’homoscédasticité (la variance sera alors constante)

Mais on peut aller plus loin

  • Le modèle linéaire généralisé

Plusieurs modèles peuvent etre estimés, en changeant la loi de la variable à expliquer, et la fonction lien,

> reg2=glm(Y~1+X,data=base,family=poisson(link="identity"))
> Y2=predict(reg2,newdata=nbase,type="response")
> reg3=glm(Y~1+X,data=base,family=poisson(link="log"))
> Y3=predict(reg3,newdata=nbase,type="response")
> reg4=glm(Y~1+X,data=base,family=gaussian(link="log"))
> Y4=predict(reg4,newdata=nbase,type="response")
> sigma=sqrt(summary(reg4)$dispersion)

Pour le modèle Poissonnien avec un lien identité, on obtient

On obtient ainsi une variance qui augmente avec la prédiction,

Pour une régression de Poisson avec un lien logarithmique,

i.e. pour nos quatre prédictions

On peut comparer avec une prédiction d’un modèle Gaussien avec un lien logarithmique,

i.e. pour les quatre prédictions

C’est pas ma faute: c’est l’autre là, devant !

Les vacances sont l’occasion de se poser des questions que l’on n’a pas l’occasion de se poser d’ordinaire. Par exemple quand je m’étais demandé pourquoi, dans un petit village perdu, je continuais de fermer la porte de la voiture a clé (ici), il y a maintenant plusieurs mois.

Mais surtout, les longs trajets en voiture sont une opportunité pour discuter avec le copilote (qui en l’occurrence était la mère de mes enfants), et d’essayer de la convaincre du bien-fondé de certaines de nos décisions. Après une petite heure de route, la semaine dernière, je l’ai entendue m’interpeller “ralentis,  il y un radar et tu est à 115 !”. Bon, je dois avouer que l’argument “non non, je respecte la vitesse, mais de ton point de vue tu dois avoir une vision déformée et biaisée du tableau de bord car si tu étais pile en face, tu verrais que je suis un peu en dessous de 105 km/h” me marche plus…. car elle avait raison, j’étais bien à 115 km/h. Mais instinctivement, ma réponse a été la même que des millions d’automobilistes qui doivent faire face aux foudres de leur copilote: “oui, mais c’est pas ma faute”… Sans y réfléchir, cette excuse plus-que-centenaire, que je sortais a mes parents après avoir cassé la porte d’entrée après avoir donné un coup de pied dedans, ou que me ressortent mes enfants après lancé leur avion Playmobil par la fenêtre, est sortie toute seule de ma bouche. Et au moment de la sortir, j’ai compris immédiatement dans quel pétrin je me mettais… surtout que dans les yeux de ma copilote, je sentais pétiller le goût de la victoire (ou plutôt du je-sens-qu’on-va-rigoler-un-peu-avec-une-histoire-a-dormir-debout). “ah oui ? mais c’est la faute de qui alors… ?”.

Bon, alors la, on commence à visualiser en une fraction de seconde la liste des boucs émissaires possibles (le voisin, les enfants, mes parents, la belle-mère, la petite vieille du coin de la rue, la police, la pluie , le soleil, la fatigue, internet, la société, Nicolas Sarkozy, les normes IFRS, que sais-je encore….), et je ne sais pourquoi, j’ai retenu “ben c’est de la faute de l’autre, là, devant” (encore un 75 tiens !).
C’est vrai quoi, si ça se trouve, avec un effet d’entraînement, il m’aspire et sans que je fasse rien, il me fait gagner 10 km/h. Bon, mes compétences en mécanique des fluides sont limitées. Donc j’aurai du mal a justifier cette théorie (un peu fumeuse, j’en conviens). J’ai personnellement opté pour “quand on conduit. on ne peut pas passer son temps le nez sur le tableau de bord… donc on regarde très occasionnellement… Et la plupart du temps, on se positionne par rapport aux autres, et en particulier la voiture juste devant… on essaye de maintenir une distance à peu près constante, et donc on roule a la même vitesse. Si j’allais trop vite, c’est que lui allait trop vite !”. Et toc… pas mal ma théorie, non ?

Bon, maintenant va falloir le vérifier.

Mais ça tombe bien, car il y a quelque temps, j’avais travaillé sur des données de passage de voitures, dans la banlieue parisienne (ici).
Dans un premier temps, on peut faire un petit test d’indépendance, dichotomique (je dépasse la vitesse autorisée, versus la voiture devant dépasse la vitesse autorisée). Plaçons nous un moment ou ça ne bouchonne pas (disons entre 15 et 16 heures, un jour de semaine) et regardons la fréquence de véhicules en excès de vitesse, ainsi que la même fréquence pour les véhicules immédiatement derrières,

> M
     [,1] [,2]
[1,]  879  104
[2,]  103   36

avec par lignes le véhicule qui suit, et par colonne celui qui précède. Ici, 103 véhicules dépassaient la vitesse autorisée sans que le véhicule qui précédait ne dépasse les 110 km/h. Un test d’indépendance du chi-deux donne ici,

> chisq.test(M)
 
	Pearson's Chi-squared test with Yates' continuity correction
 
data:  M 
X-squared = 24.7862, df = 1, p-value = 6.405e-07
 
 
> chisq.test(M)$expected
         [,1]      [,2]
[1,] 860.3440 122.65597
[2,] 121.6560  17.34403

Bref, il semble que l’on rejette l’hypothèse d’indépendance.

On peut d’ailleurs le faire à plusieurs heures de la journée, et étudier l’évolution de la p-value du test d’indépendance,

Moralité, on rejette indépendance entre les vitesses de véhicules qui se suivent aux heures normales (i.e. en dehors des heures de pointe, et la nuit).

On peut aussi regarder en fonction de la vitesse du premier véhicule, la vitesse du véhicule juste derrière, avec une régression lissée,

Enfin, on peut aussi tenter une régression logistique, afin de quantifier la probabilité de faire un excès de vitesse, conditionnellement a la vitesse du véhicule qui le précède,

Si au lieu de le faire sur une journée, on regarde sur 38 jours (c’est ce dont je dispose dans ma base), on a

Autrement dit, la probabilité que je dépasse la vitesse autorisée (ici 110 km/h) dépend très fortement – et de manière sensiblement croissante – de la vitesse du véhicule qui me précède…

Moralité, mon argument pourrait être valide: si on va vite, c’est en partie parce que la voiture devant va vite, Mais si mon argumentaire me convainc (et pourrait – presque – convaincre ma femme), on verra s’il convainc la police. Car s’ils me contactent pour m’avoir flashé a une vitesse excessive, j’aurais envie de leur sortir ma théorie…. En plus ma théorie est basée sur une étude de corrélation (je fais des régressions toutes simples), et pas de causalité… ma conclusion pourrait alors parfaitement être inversée, ou complétée: je pourrais tout aussi légitimement expliquer à ma copilote “ben si j’allais trop vite, c’est de la faute de l’autre, la, derrière… il me collait, alors j’allais vite”.
Bref, il semble donc que l’excuse du “c’est pas ma faute, c’est l’autre” a de beaux jours devant elle… surtout que les élections approchent en France, donc je pense qu’on n’a pas fini de l’entendre !

Multivariate probit regression using (direct) maximum likelihood estimators

Consider a random pair http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-01.gif of binary responses, i.e. http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-02.gif with http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-03.gif taking values 1 or 2. Assume that probability http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-04.gif can be function of some covariates.

  • The Gaussian vector latent structure

A standard model is based a latent Gaussian structure, i.e. there exists some random vector http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-06.gif such that http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-07.gif if http://freakonometrics.hypotheses.org/files/2015/12/biv-prob-08.gif is lower than a given threshold, and 1 otherwise.
As in standard probit models, assume that

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

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

> logV=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+ sigma=matrix(c(1,CORRELATION,CORRELATION,1),2,2)
+     a11=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
+ for(i in 2:nrow(z)){a11=c(a11,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
+     a10=pnorm(z[1,1],sd=sqrt(sigma[1,1]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a10=c(a10,pnorm(z[i,1],sd=sqrt(sigma[1,1]))-pmnorm(z[i,],varcov=sigma))}
+     a01=pnorm(z[1,2],sd=sqrt(sigma[2,2]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a01=c(a01,pnorm(z[i,2],sd=sqrt(sigma[2,2]))-pmnorm(z[i,],varcov=sigma))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     1*log(a01) +
+     2*log(a10) +
+     3*log(a00) )
+ }
> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par

(the code is a bit long since I had trouble working properly with matrices – or more precisely to vectorize my functions – so I used loops… I am sure it is possible to write a better code).
It is possible to generate samples (based on that specific model) to check that we can actually derive proper maximum likelihood estimators,

> library(mnormt)
> set.seed(1)
> n=1000
> r=0.5
> X1=runif(n)
> X2=rnorm(n)
> Y1S=1+5*X1
> Y2S=8-5*X1
> RES=rmnorm(n,mean=c(0,0),varcov=matrix(c(1,r,r,1),2,2))
> YS=cbind(Y1S,Y2S)+RES
> Y1=(YS[,1]>quantile(YS[,1],.5))*1
> Y2=(YS[,2]>quantile(YS[,2],.5))*1
> base=data.frame(i,Y1,Y2,X1,X2,YS)
> head(base)
  i Y1 Y2        X1          X2      Y1S      Y2S
1 1  0  0 0.2655087  0.07730312 3.177587 5.533884
2 2  0  0 0.3721239 -0.29686864 1.935307 5.089524
3 3  1  0 0.5728534 -1.18324224 4.757848 5.172584
4 4  1  0 0.9082078  0.01129269 4.600029 3.878225
5 5  0  1 0.2016819  0.99160104 2.547362 6.743714
6 6  1  0 0.8983897  1.59396745 5.309974 4.421523

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

>  reg1=glm(Y1~X1+X2,data=base,family=binomial)
>  reg2=glm(Y2~X1+X2,data=base,family=binomial)
> summary(reg1)
 
Call:
glm(formula = Y1 ~ X1 + X2, family = binomial, data = base)
 
Deviance Residuals:
Min        1Q    Median        3Q       Max
-2.90570  -0.50126  -0.00266   0.49162   2.78256
 
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.291725   0.267149 -16.065   <2e-16 
X1           8.656836   0.510153  16.969   <2e-16 ***
X2           0.007375   0.090530   0.081    0.935
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1386.29  on 999  degrees of freedom
Residual deviance:  726.48  on 997  degrees of freedom
AIC: 732.48

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

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

Here, the optimization yields,

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

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

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

> matrix(coefficients(REG)[c(1:2,4:7)],2,3)
          [,1]      [,2]       [,3]
[1,] -2.451699  4.908617 0.01599183
[2,]  2.241964 -4.539951 0.10615443

The correlation here is also the same

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

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

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

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

  • marginally: for instance if we use a logistic cumulative distribution function, then we will have a bivariate logit regression
  • in terms of dependence structure: it is possible to consider another copula than the gaussian one, e.g. Gumbel’s copula (also called the bivariate logistic copula), or Clayton’s

Here, the following code can be used to extend the model to non Gaussian structures,

> F=function(x,r){pmnorm(x,rep(0,length(x)),
+                 varcov=matrix(c(1,r,r,1),2,2))}
> Fx=function(x1){F(c(x1,1e40),0)}
> Fy=function(x2){Fx(x2)}
> 
> logVgen=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+     a11=F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a11=c(a11,F(z[i,],r=CORRELATION))}
+     a10=Fx(z[1,1])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a10=c(a10,Fx(z[i,1])-F(z[i,],r=CORRELATION))}
+     a01=Fy(z[1,2])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a01=c(a01,Fy(z[i,2])-F(z[i,],r=CORRELATION))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     11*log(a01) +
+     12*log(a10) +
+     13*log(a00) )
+ }
>
> beta0=c(0,1,1,1,1,1,1)
> (OPT=optim(fn=logVgen,par=beta0,method="BFGS")$par)
[1]  0.52613820 -2.45172059  2.24196154  4.90863292 -4.53994592  0.01600769
[7]  0.10614807
There were 23 warnings (use warnings() to see them)

E.g.

> library(copula)
> F=function(x,r){pcopula(pnorm(x),
               claytonCopula(2, r))}
> Fx=function(x1){F(c(x1,1e40),0)
}
> Fy=function(x2){Fx(x2)}
  • An application to school tests

Consider the following dataset,

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

with here maths versus writing, with girls in red and boys in blue, where variables here are

  female :
    0: male
    1: female
  race :
    1: hispanic
    2: asian
    3: african-amer
    4: white
  ses :
    1: low
    2: middle
    3: high
  schtyp : type of school
    1: public
    2: private
  prog : type of program
    1: general
    2: academic
    3: vocation
  read : reading score
  write : writing score
  math : math score
  science : science score
  socst : social studies score

We can try to understand correlation between math and writing skills. Covariates can be the sex of the child, and his reading skills. The question will then be: are good students in maths and writing simply students that can read well ?

Here the code is simply

> W=hsb2$write>=50
> M=hsb2$math>=50
> base=data.frame(Y1=W,Y2=M,
+             X1=hsb2$female,X2=hsb2$read)
>
> library(Zelig)
> REG= zelig(list(mu1=Y1~X1+X2,
+             mu2=Y2~X1+X2,
+     rho=~1),
+     model="bprobit",data=base)
> summary(REG)
 
Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
    rho = ~1), model = "bprobit", data = base)
 
Pearson Residuals:
                Min        1Q  Median      3Q    Max
probit(mu1) -4.7518 -0.502594 0.15038 0.53038 1.8592
probit(mu2) -3.4243 -0.653537 0.23673 0.67011 2.6072
rhobit(rho) -4.9821  0.010481 0.13500 0.40776 2.9171
 
Coefficients:
                  Value Std. Error  t value
(Intercept):1 -5.484711   0.787101 -6.96825
(Intercept):2 -4.061384   0.633781 -6.40818
(Intercept):3  1.332187   0.322175  4.13497
X1:1           1.125924   0.233550  4.82092
X1:2           0.167258   0.202498  0.82598
X2:1           0.103997   0.014662  7.09286
X2:2           0.082739   0.012026  6.88017
 
Number of linear predictors:  3
 
Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
 
Dispersion Parameter for binom2.rho family:   1
 
Residual Deviance: 364.51 on 593 degrees of freedom
 
Log-likelihood: -182.255 on 593 degrees of freedom
 
Number of Iterations: 3
> (exp(summary(REG)@coef3[3])-1)/(exp(
summary(REG)@coef3[3])+1)
[1] 0.5824045

with a remaining correlation among residuals of 0.58. So with only the sex of the student, and his or her reading skill, we cannot explain the correlation between maths and writing skills. With our previous code, we have here

> beta0=c((exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1),
+      summary(REG)@coef3[c(1:2,4:7),1])
> beta0
              (Intercept):1 (Intercept):2          X1:1          X1:2
0.58240446   -5.48471133   -4.06138412    1.12592427    0.16725842
X2:1          X2:2
0.10399668    0.08273879
> (OPT=optim(fn=logV,par=beta0,method="BFGS")$par)
(Intercept):1 (Intercept):2          X1:1          X1:2
0.5824045    -5.4847113    -4.0613841     1.1259243     0.1672584
X2:1          X2:2
0.1039967     0.0827388

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

> (OPT=optim(fn=logV,par=beta0/2,method="BFGS")$par)
              (Intercept):1 (Intercept):2          X1:1          X1:2
   0.58233360   -5.49428984   -4.06839571    1.12696594    0.16760347
         X2:1          X2:2
   0.10417767    0.08287409
There were 12 warnings (use warnings() to see them)

So once again, it is possible to optimize numerically a likelihood function, and it works.

  1. Y[,1]==0)&(Y[,2]==1 []
  2. Y[,1]==1)&(Y[,2]==0 []
  3. Y[,1]==0)&(Y[,2]==0 []
  4. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==0 []
  5. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==0 []
  6. Y[,1]==0)&(Y[,2]==0)&(Y[,3]==1 []
  7. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==0 []
  8. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==1 []
  9. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==1 []
  10. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==1 []
  11. Y[,1]==0)&(Y[,2]==1 []
  12. Y[,1]==1)&(Y[,2]==0 []
  13. Y[,1]==0)&(Y[,2]==0 []

Too large datasets for regression ? What about subsampling….

recently, a classmate working in an insurance company told me he had too large datasets to run simple regressions (GLM, which involves optimization issues), and that they were thinking of a reward for the one who will write the best R-code (at least the fastest). My first idea was to use subsampling techniques, saying that 10 regressions on 100,000 observations can take less time than a regression on 1,000,000 observations. And perhaps provide also better results…

  • Time to run a regression, as a function of the number of observations

Here, I generate a dataset as follows

http://freakonometrics.hypotheses.org/files/2016/11/largesamp01.png

and we fit

http://freakonometrics.hypotheses.org/files/2016/11/largesamp02.png

where http://freakonometrics.hypotheses.org/files/2016/11/largesamp03.png is a spline function (just to make it as general as possible, since in insurance ratemaking, we include continuous variates that do not influence claims frequency linearly in the score). Yes, there might be also useless variables, including one of them which is strongly correlated with one that has an impact in the regression. The code to generate the dataset is simply

> n=10000
> X1=rexp(n)
> X2=sample(c("A","B","C"),size=n,replace=TRUE)
> X3=runif(n)
> Z=rmnorm(n,c(0,0),matrix(c(1,0.8,.8,1),2,2))
> X4=Z[,1]
> X5=Z[,2]
> X6=X1^2
> E=runif(n)
> lambda=.2*X5-4*dbeta(X3,2,5)+X1+
+1*(X2=="A")-2*(X2=="B")-5*(X2=="C")
> Y=rpois(n,exp(lambda))
> base=data.frame(Y,X1,X2,X3,X4,X5,X6,E)

We would like the study the time it takes to run a regression, as a function of the size (i.e. the number of lines http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png) of the dataset.

> system.time( glm(Y~bs(X1)+X2+X3+X4+
+ X5+X6+offset(log(E)),family=poisson,
+ data=base) )
utilisateur     système      écoulé
0.25        0.00        0.25

Here, the time I look at is the last one. But so far, it was rather simple, but it is not the best model I can get. Let us use a stepwise (backward) variable selection,

> system.time( step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=poisson,
+ data=base)) )
Start:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))
Step:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + offset(log(E))
Df Deviance    AIC
<none>        2236.0 2882.1
- X5      1   2240.1 2884.2
- X4      1   2244.1 2888.2
- X3      1   4783.2 5427.3
- X2      2   5311.4 5953.5
- bs(X1)  3   6273.7 6913.8
utilisateur     système      écoulé
1.82        0.03        1.86

Finally, from the first regression, we have points in black (based on 200 simulated datasets), and with a stepwise procedure, we have the points in red.

i.e. it might look linear (proportional), but if it was linear, then on a log-log scale, we should have also straigh lines, with slope 1,

Actually, it looks like a convex function.

The interpretation of that convexity might lead to misinterpretation. On the graph below on the left, on a dataset two times bigger than the previous one (black point) will be less than two times longer to run, while on the right, it will be more than two timess longer,

Convexity can simply be interpreted as “too large datasets take time, and too small too…”. Which is a first step: it should be interesting, in some cases, to run several regressions on smaller datasets….

  • Running 100 regressions on 100 lines, or running 1 regression on 10,000 lines ?

Here, we have datasets with http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=200,000 lines. The questions is how long will it take if we subdived into http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png subsamples (of equal size), and run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regressions ?

> nk=trunc(n/k)rep(1:k,each=nk); nt=nk*k
> base=data.frame(Y[1:nt],X1[1:nt],
+ X2[1:nt],X3[1:nt],X4[1:nt],X5[1:nt],
+ X6[1:nt],E[1:nt],classe)
> system.time( for(j in 1:k){
+  glm(Y~bs(X1)+X2+X3+X4+X5+
+ X6+offset(log(E)),family=poisson
+ ,data=base,subset=classe==j) })
utilisateur     système      écoulé
1.31        0.00        1.31
> system.time( for(j in 1:k){
+      step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=
+ poisson,data=base,subset=classe==j)) })
Start:  AIC=183.97
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))

[…]

  Df Deviance    AIC
<none>        117.15 213.04
- X2      2   250.15 342.04
- X3      1   251.00 344.89
- X4      1   420.63 514.53
- bs(X1)  3   626.84 716.74
utilisateur     système      écoulé
11.97        0.03       12.31

On the graph below, we have the time (y-axis, here on a log scale) it took to run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regression on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png, as function of http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png (x-axis), including the time it took to run the regression on a dataset of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png which is the concentration of dots on the left (i.e. http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png=1), both on the 6 regressors – in black – and with a strepwise procedure – in red. One has to keep in mind that I did not remove the printing option in the stepwise procedure, so it might be difficult to compare the two clouds (black vs. red). Nevertheless, we clearly see that if we run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regression on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png, when http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png is not too large, i.e. less than 10 or 15, it is not longer than the regression on http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=200,000 lines.

So here we see that running 100 regressions on 2,000 lines is longer than running 1 regression on 200,000 lines… But maybe we are not comparing things that are actually comparable: what if it takes a bit longer, but we strongely improve the quality of our estimators ?

  • What about the quality of the output ?

Here, we consider only one dataset, with http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=100,000 lines (just to make it run a bit faster). And http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png=20 subsets. Recall that the generated dataset is from

http://freakonometrics.hypotheses.org/files/2016/11/largesamp01.png

and we fit

http://freakonometrics.hypotheses.org/files/2016/11/largesamp02.png

Here, we plot here http://freakonometrics.hypotheses.org/files/2016/11/largesamp07.png and a confidence interval, defined as

http://freakonometrics.hypotheses.org/files/2016/11/largesamp08.png

The lightblue segment is the initial estimator, while the blue one is obtained from the stepwise procedure. The grey area represent the estimation on the overall sample, while the http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png segments on the right are the http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png estimators (each on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png).

We can see that we have much more volatility on those http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png estimators, but the average (horizontal doted lines) are not so bad… The true value (i.e. the one used to generate the dataset is the dotter black horizontal line).
And if we repeat that on 1,000 simulated dataset, we obtaind the following distribution for http://freakonometrics.hypotheses.org/files/2016/11/largesamp07.png (blue line), so we have an unbiased estimator of our parameter (the verticular line being here the true value), here including a stepwise procedure,

But if we add the the red curve is the average of the http://freakonometrics.hypotheses.org/files/2016/11/largesamp09.png the previous one being now the clear blue line in the back, we see that taking average of estimators on subsamples is not bad at all, on the contrary,

and for those who think that the stepwise procedure is a mistake, here is what we get without it,

So what we can see is that running 20 regressions can take (a little) more time (from what we’ve seen earlier) than running only one on the whole dataset…. but it provides better estimates. So the tradeoff is not that simple, and maybe running several regressions on huge datasets can be a proper alternative.

Qui peut m’aider à comprendre les sorties de SAS ?

Je m’étais promis que j’évoquerais une bizarrerie rencontrée avec SAS lors d’une formation…. Écrire ce billet permettra à ceux qui auraient des éléments d’explication de poster un commentaire.
Pour cela, comparons une régression logistique faite avec deux outils différents, sous SAS,

  • avec la procédure logistique

Le code pour faire une régression logistique ressemble à ça

PROC LOGISTIC DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / selection=stepwise;
RUN; QUIT;

ce qui donne la sortie suivante (je passe l’introduction pour insister sur les coefficients)

                                 The LOGISTIC Procedure

                   Analyse des estimations de la vraisemblance maximum
                                                     Erreur         Khi 2
  Paramètre                    DF    Estimation         std       de Wald    Pr > Khi 2

  Intercept                     1        1.7833      0.0676      696.9022        <.0001
  sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
  age_soc      1_AGESOC_-60     1        0.4578      0.0667       47.1020        <.0001
  fract_paiemt Annuel           1        0.6021      0.0997       36.4862        <.0001
  fract_paiemt Mensuel          1       -0.5410      0.0842       41.2342        <.0001
  • avec la procédure genmod (car la régression logistique est un glm)

On peut faire exactement la même chose (théoriquement) en ajustement un modèle GLM,

PROC GENMOD DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / dist = binomial;
RUN;

et la sortie ressemble à ça

                                  The GENMOD Procedure
                      Analyse des résultats estimés de paramètres

                                                  Erreur      Wald 95Limites
Paramètre                     DF   Estimation   standard      de confiance %       Khi 2
Intercept                      1       1.5073     0.1501     1.2131     1.8014    100.85
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000       .
age_soc        1_AGESOC_-60    1       0.9156     0.1334     0.6542     1.1771     47.10
age_soc        Z_AGESOC_+60    0       0.0000     0.0000     0.0000     0.0000       .
fract_paiemt   Annuel          1       0.6634     0.1770     0.3165     1.0104     14.05
fract_paiemt   Mensuel         1      -0.4798     0.1510    -0.7759    -0.1838     10.09
fract_paiemt   Semestriel      0       0.0000     0.0000     0.0000     0.0000       .
Scale                          0       1.0000     0.0000     1.0000     1.0000
  • comparaison des deux sorties

Si on regarde l’impact du sexe par exemple, dans la première sortie on peut lire

sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
alors que dans la seconde sortie, on a
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000

On dira ce qu’on veut, mais moi je trouve cette différence troublante…. Dans la seconde sortie, le coefficient vaut le double de l’autre….
Alors SAS semble s’y retrouver car si on lui demande d’afficher le score prédit pour un individu au hasard (le premier de la base par exemple), les prédictions sont très proches,

                           fract_                                 proba1_       proba1_
  Obs  sexe_soc   age_soc  paiemt      SPOCAM  proba1_logit        
1      Homme          71  Annuel         0      0.10242637    0.10241302

Si quelqu’un sait interpréter ce qui est fait avec cette procédure logistique (car R donne la même chose que la sortie GLM), je suis preneur…..


Vitesse et comportement dangereux au volant (partie 3)

Après plusieurs billets sur les accidents de la route (ici et ), parlons un peu de vitesse des véhicules. J’ai pu récupérer des données via l’inrets, d’observations près de Saclay en région parisienne, sur une 2×2 limitée à 110 km/h.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso3/.n118-saclay2_m.jpg

J’ai ainsi le passage de plusieurs millions de véhicules, sur quelques jours, incluant tous les véhicules. Mais si j’ai signé des papiers quant à la confidentialité des données, je dois pouvoir mettre des dessins construits à partir de ces données…

  • A quelle vitesse roulent les conducteurs parisiens ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, qui poussent les gens à ralentir. La courbe en noir est la vitesse moyenne des véhicules qui passent sur ce tronçon, et la courbe en bleu, la vitesse des 5% des véhicules les plus rapides,

Ces pics sont légèrement atténués le week end,

avec ici un trafic très dense tout l’après midi, entre 15 et 20 heures.

  • Quelle proportion de véhicules dépasse la limite autorisée ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, avec en gras la proportion des véhicules qui dépasse 110 km/h, et en trait fin, la proportion qui dépasse 115 km/h (soit 5% de plus),

avec des choses finalement assez proches le week end,

Autrement dit, aux heures de pointe, peu de monde dépasse la vitesse autorisée, mais vers midi, ou surtout la nuit, environ un véhicule sur quatre dépasse la vitesse autorisée (mais seulement un sur six la dépasse de plus de 5%).

  • Vitesse et respect des distances de sécurité

La base est très complète, puisque j’ai tous les véhicules qui sont passés près de Saclay, avec l’heure d’observation, la vitesse, ce qui permet d’obtenir la distance entre deux véhicules qui se suivent (en supposant que le premier véhicule maintienne sa vitesse constante).
Rappelons qu’il existe un lien théorique entre la vitesse d’un véhicule et la distance de freinage (expliqué ici ou ). Pour faire simple, la distance de freinage https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein02.png est la somme d’une distance parcourue avant de réagir https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein04.png, et d’une distance de freinage proprement dite https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein03.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein01.png

Si on suppose que la personne met une seconde à réagir, la distance de réaction est simplement la vitesse https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein05.png(exprimée en m/h), i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein06.png

Quant au freinage à proprement parler, rappelons que les équations qui caractérisent le mouvement d’un objet sont de la forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein10.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein11.png
https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein12.png

Aussi, en mettant la deuxième au carré, et en combinant avec la troisième, on a

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein13.png

avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein15.png (on cherche à arrêter le véhicule) et https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein14.png, aussi, avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein16.png, on obtient quelque chose qui peut s’écrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein17.png

à condition d’avoir des vitesses en km/h. Aussi en agglomération, il faut 23 m pour s’arrêter. Et sur notre route à Saclay, pour une personne qui roule à 110 km/h, il faut compter 110 m (ce qui correspond au message véhiculé par la sécurité routière).
Sur nos données, on peut représenter en rouge les individus qui ne respectent pas la distance de sécurité.

La courbe bleue est la distance moyenne des véhicule en fonction de la distance qui les sépare du véhicule qui les précède,
Aussi, en faisant une régression logistique (lissée), on obtient que 60% des conducteurs qui roulent entre 50 et 110 km/h ne respectent pas les distances de sécurité.

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 de R

> lm(dist~speed,cars)
Coefficients:
(Intercept)        speed 
    -17.579        3.932

Et on peut aussi vérifier qu’il minimise 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.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-1.png

En effet, sinon

https://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

https://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.

https://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é,
https://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.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-3.png
https://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 https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-5.png converge en probabilité vers https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-6.png, parfois noté

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-7.png

si

https://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.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-9.png

si l’espérance des variables est finie, alors

https://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,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-11.png

On parlera alors d’estimateur convergent. On peut aussi montrer que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-12.png

Mais formellement,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/latex-reponse-mail-13.png

si

https://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.

Delta method avec R, ou transformation des paramètres

Un petit billet aujourd’hui sur une question qui m’a été posée (et qui m’a permis de découvrir des nouvelles fonctions sous R). On suppose que l’on a un modèle économétrique de la forme

que l’on peut encore écrire

Supposons que ce soit la formulation en terme de ces derniers coefficients qui nous intéresse, et que l’on cherche un intervalle de confiance.
La Delta Method nous dit que si on dispose d’un estimateur asymptotiquement gaussien (mais c’est le cas pour la plupart des estimateurs que l’on construit en pratique), tout transformation bijective de cet estimateur reste asymptotiquement Gaussien. Bref, on va pouvoir utiliser ces résultats pour en déduire un intervalle de confiance (asymptotique) pour le nouveau paramètre. En effet, l’estimateur par moindre carrés vérifie la normalité asymptotique (sous quelques hypothèses)

ce qui donne avec les notations précédentes

Autrement dit, on peut facilement en déduire un intervalle de confiance (asymptotique). Pour cela, on a trois méthodes (au moins) avec R. Pour la suite, je vais considerer la régression sur la base de données cars (comme toujours), et on s’intéresse à l’estimation de

 

(ne me demandez pas pourquoi, il me fallait un exemple numérique, j’ai pris un truc nonlinéaire pour qu’au moins ça puisse avoir un vague intérêt).

  • En faisant un peu de calcul formel

Pour les calculs symboliques, il y a la fonction D dans le package de base,

> f <- expression(sin(cos(x + y^2)))
> ( Dfx <- D(f, "x") )
-(cos(cos(x + y^2)) * sin(x + y^2))

ce qui pourra alors être utilisé pour faire les calculs de la variance asymptotique.

> x=c(0,1,2); y=1
> eval(D.sc)
[1] -0.7216061 -0.8316919 -0.0774320

On peut également utiliser cette fonction pour faire plein de jolis dessins,

> x <-seq(-5, 5, length=100)
> fx=x^2
> Df <-D(expression(x^2), "x")
> Dfx<-eval(Df)
> plot(x, fx, type = "l", las = 1,ylim=c(-2,2))
> for (i in 1:length(x)) {
+ abline(a = fx[i] - Dfx[i] * x[i], b = Dfx[i],
+ col = rainbow(100)[i])
+ }

Mais bon, c’était juste pour rigoler un peu… On peut utiliser cette fonction pour calculer l’écart type de l’estimateur de la grandeur qui nous intéresse.

> reg1 = lm(dist~speed, data = cars)
> estmean = coef(reg1)
> estvar = summary(reg1)$cov.unscaled * summary(reg1)$sigma^2
> f = function(x,y){
+ -y/(x+y)
+ }
> X=f(estmean[1],estmean[2])
> g=expression(-y/(x+y))
> x<-as.numeric(estmean[1]); y<-as.numeric(estmean[2])
> V=eval(c(eval(D(g, "x")),eval(D(g, "y"))))
> sqrt(t(V) %*% estvar %*% V)
          [,1]
[1,] 0.1063255

Voilà pour la première méthode. Petite remarque: pour faire le calcul symbolique de la dérivée, j’utilise expression et non pas function.

  • En utilisant la fonction deltamethod de la librarie msm

Ici on commence comme avant, sauf que la phase de calcul est immédiat…

> reg2 <- lm(dist~speed, data = cars)
> estmean = coef(reg2)
> estvar = summary(reg2)$cov.unscaled * summary(reg2)$sigma^2
> deltamethod (~ -x2 / (x1 + x2), estmean, estvar) 
[1] 0.1063255

et d’ailleurs, ça peut s’appliquer pour un paramètre en dimension 2

> deltamethod (list(~ -x2 / (x1 + x2), ~ x1+x2), estmean, estvar) 
[1] 0.1063255 6.3664368

En revanche, on obtient des estimateur des écart-types, et pas de la matrice de variance-covariance des estimateurs. Autrement dit on ne peut pas en déduire l’ellipse de confiance, seulement des intervalles de confiance marginaux (en utilisant la normalité asymptotique).

  • En utilisant la fonction delta.method de librarie alr3

Pour finir, une dernière fonction. La seule chose qui change c’est la synthaxe. Mais la sortie est plus élégante,

> reg3 = lm(dist~speed, data = cars)
> delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255

En fait, l’écart-type indiqué ici est précisément celui qui est utilisé pour écrire l’intervalle de confiance. Dans cet exemple, l’intervalle de confiance approché (car on suppose que la normalité – en théorie seulement asymptotique – est vérifiée) et estimé (car les paramètres sont inconnus et estimés) est de la forme suivant, pour un seuil de confiance à 95%

> DM=delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255 
> DM$estimate+qt(.975,df=nrow(cars)-1)*DM$se
[1] 0.5018277
> DM$estimate+qt(.025,df=nrow(cars)-1)*DM$se
[1] 0.07448936

C’est à dire qu’on a un intervalle de confiance de la forme [0.074;0.502].