L'examen intra avait lieu ce matin. L'examen est en ligne ici, avec les annexes, en ligne là. Des éléments de correction seront mis en ligne très bientôt.
Tag - régression
Wednesday, November 7 2012
Examen ACT6420 sur les modèles de régression
By arthur charpentier on Wednesday, November 7 2012, 12:03 - ACT6420-A2012
Thursday, October 25 2012
Régression sur des variables orthogonales
By arthur charpentier on Thursday, October 25 2012, 08:47 - ACT6420-A2012
Hier, on a évoqué en cours le théorème de Frisch-Waugh (évoqué sur ce blog il y a quelques mois). Derrière, il y a l'estimation de modèles sur des sous-ensembles de valeurs. Formellement, on a un modèle de la forme
Les équations normales sont ici (cf transparents du cours, ou équations écrites au tableau) , que l'on peut écrire sous la forme d'un système,
Le point intéressant est que si les variables dans le premier ensemble et les variables dans le second sont orthogonales, i.e. et
, alors les termes diagonaux du système disparaissent, et on retrouve les équations normales des régressions simples, i.e.
On a ici une conséquence du théorème de Frisch-Waugh: si les variables sont orthogonales, les estimateurs dans la régression multiple coïncident avec ceux obtenus avec des régressions simples.
Afin de mieux clarifier ce point, écrivons ce qui se passerait avec seulement deux variables. On va associer la constante à la première, i.e.
Dans ce cas, les équations normales sont
Bref, il faut que et
soient orthogonales pour que le terme sur la seconde diagonale disparaisse. C'est ce qu'on retrouve dans la plupart des livres d'économétrie, il suffit de chercher le théorème de Frisch-Waugh dans l'index (même si le théorème de Frisch-Waugh va beaucoup plus loin puisqu'il explique quels ajustement faire pour lier les estimateurs dans une régression simple, et dans le cas multiple).
Mais est-ce un résultat purement théorique ou peut-on l'invoquer dans un cas pratique ? Autrement dit, si deux variables semblent orthogonales, que se passent-il numériquement ? Si je simule des variables explicatives indépendantes, les estimateurs coïncident-ils (comme le suggère la théorie) ?
Considérons un modèle linéaire avec deux variables explicative que l'on supposera indépendantes, mais ayant toutes deux un impact sur la variable explicative. Le plus simple est de considérer un triple de variables telles que
où désigne des valeurs non nulles. Par exemple,
> set.seed(1) > n=100 > S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3) > S [,1] [,2] [,3] [1,] 1.0 0.4 0.3 [2,] 0.4 1.0 0.0 [3,] 0.3 0.0 1.0 > M=c(1,3,6) > Z=rmnorm(n,M,varcov=S) > Y=Z[,1]; X1=Z[,2]; X2=Z[,3] > base=data.frame(Y,X1,X2) > cor(base) Y X1 X2 Y 1.0000000 0.37796580 0.27859381 X1 0.3779658 1.00000000 -0.05060539 X2 0.2785938 -0.05060539 1.00000000On a simulé ici deux variables indépendantes, mais toutes deux corrélées avec la variable d'intérêt. On notera que numériquement
(autrement dit, en moyenne, on simulera des données indépendantes, mais dans ce cas précis, les données ne sont pas rigoureusement orthogonales).
Mais un test rapide nous pousse à accepter l'hypothèse d'indépendance entre les deux variables (sans trop d'hésitations).
dans le modèle multiple est légèrement différent de obtenu dans le modèle simple. Si on regarde la distribution de nos estimateurs (n'oublions pas qu'il y a une marge d'erreur), ils sont proches, ou au moins statistiquement proches (avec à gauche la distribution de et à droite celle de 
)
Si on génère 100 bases de données avec ce modèle (avec 500 observations à chaque fois) on obtient les estimateurs suivants (à gauche les et à droite les , avec en abscisse le modèle multiple et en ordonnées, le modèle simple).
(la valeur centrale étant la vraie valeur, utilisée pour les simulations). On pourrait se demander si les valeurs éloignées du barycentre (ou de cette première bissectrice) correspondent à des cas où l’hypothèse d'indépendance entre les deux variables est rejetée ? On peut mettre des couleurs différentes en fonction de la 
-value: plus la -value est faible (rejet de l'hypothèse d'indépendance) plus la couleur sera vive (et rouge), et plus la -value est forte (acceptation de l'hypothèse d'indépendance) plus la couleur sera terne (et jaune, voire blanche),
Le code est ici
Autrement dit le fait que l'estimation dans le modèle multiple et les modèles simples diffèrent beaucoup n'a rien a voir avec la -value d'un test d'indépendance entre les deux variables... étonnant non ? En revanche, si on rajoute une observation pour rendre le modèle parfaitement non-corrélé, par exemple,
> cor.test(X1,X2) Pearsons product-moment correlation data: X1 and X2 t = -0.5016, df = 98, p-value = 0.6171 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.2445923 0.1472766 sample estimates: cor -0.05060539Si on regarde maintenant la régression,
> R=lm(Y~X1+X2,data=base) > summary(R) Call: lm(formula = Y ~ X1 + X2, data = base) Residuals: Min 1Q Median 3Q Max -1.87885 -0.57656 0.07957 0.52719 1.77407 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.71180 0.28246 16.682 < 2e-16 *** X1 0.37232 0.08441 4.411 2.66e-05 *** X2 0.25975 0.07755 3.350 0.00115 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.7954 on 97 degrees of freedom Multiple R-squared: 0.2317, Adjusted R-squared: 0.2159 F-statistic: 14.63 on 2 and 97 DF, p-value: 2.803e-06Quel est le lien avec les régressions simples, sur chacune des variables ? Si on régresse sur la première, on obtient
> R1=lm(Y~X1,data=base) > summary(R1) Call: lm(formula = Y ~ X1, data = base) Residuals: Min 1Q Median 3Q Max -1.78973 -0.54060 0.04036 0.44138 1.94389 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.03165 0.27934 18.013 < 2e-16 *** X1 0.35802 0.08859 4.041 0.000106 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.8358 on 98 degrees of freedom Multiple R-squared: 0.1429, Adjusted R-squared: 0.1341 F-statistic: 16.33 on 1 and 98 DF, p-value: 0.0001058et si on régresse sur la seconde
> R2=lm(Y~X2,data=base) > summary(R2) Call: lm(formula = Y ~ X2, data = base) Residuals: Min 1Q Median 3Q Max -1.8652 -0.5507 0.1104 0.6109 2.1435 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.85053 0.12495 46.824 < 2e-16 *** X2 0.24244 0.08443 2.872 0.00501 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.867 on 98 degrees of freedom Multiple R-squared: 0.07761, Adjusted R-squared: 0.0682 F-statistic: 8.246 on 1 and 98 DF, p-value: 0.005005Si les données ne sont pas rigoureusement orthogonales, alors
|
|
|
|
|
|
> set.seed(1) > library(mnormt) > S=matrix(c(1,.4,.3,.4,1,0,.3,0,1),3,3) > BX1=BX2=matrix(NA,1000,2) > P=rep(NA,1000) > for(s in 1:1000){ + n=500 + Z=rmnorm(n,c(6,3,1),varcov=S) + Y=Z[,1]; X1=Z[,2]; X2=Z[,3] + base=data.frame(Y,X1,X2) + R=lm(Y~X1+X2,data=base) + BX1[s,1]=coefficients(R)[2] + BX2[s,1]=coefficients(R)[3] + BX1[s,2]=coefficients(lm(Y~X1,data=base))[2] + BX2[s,2]=coefficients(lm(Y~X2,data=base))[2] + P[s]=cor.test(X1,X2)$p.value + } > COL=(heat.colors(100)) > plot(BX1[,1],BX1[,2],pch=19,cex=.6,col=COL[trunc(100*P)]) > abline(a=0,b=1,lty=2)
> difference=sum(X1*X2)*(n+1)/sum(X1)-sum(X2) > X1=c(X1,0) > X2=c(X2,difference) > n=n+1 > sum(X1*X2)/n [1] 18.02172 > sum(X1)*sum(X2)/n^2 [1] 18.02172
On retrouve ici une corrélation parfaitement nulle,
> Y=c(Y,0) > base=data.frame(Y,X1,X2) > cor(base) Y X1 X2 Y 1.0000000 3.946091e-01 2.923412e-01 X1 0.3946091 1.000000e+00 -4.817519e-16 X2 0.2923412 -4.817519e-16 1.000000e+00
Si on regarde maintenant la régression multiple
> R=lm(Y~X1+X2,data=base) > coefficients(R) (Intercept) X1 X2 -1.5028358 0.3589776 0.2531358 et les deux régressions simples, > R1=lm(Y~X1,data=base) > coefficients(R1) (Intercept) X1 0.02847362 0.35897764 > R2=lm(Y~X2,data=base) > coefficients(R2) (Intercept) X2 -0.4334011 0.2531358
on obtient exactement les mêmes coefficient.
Moralité ? l'hypothèse d'orthogonalité entre les variables explicatives n'est pas à prendre à la légère: il faut une orthogonalité stricte pour que numériquement, la régression multiple soit équivalente aux régressions simples (comme le prédit la théorie). Mais de manière un peu surprenante, si on est très proche de l'indépendance, les valeurs numériques des estimateurs peut différer. Et ce, parfois de manière non négligeable, même si statistiquement, on valide l'hypothèse d'orthogonalité entre les variables explicatives.
Wednesday, October 3 2012
Modèles de régression, premier cours
By arthur charpentier on Wednesday, October 3 2012, 02:26 - ACT6420-A2012
Le premier cours aura lieu tout à l'heure, le plan de cours est en ligne. Compte tenu du nombre d'étudiants inscrits, je vais utiliser des transparents, plutôt que de continuer à faire cours au tableau. La première série de transparents est en ligne,
Mais compte tenu des délais (et de la rentrée un peu intense), je n'aurais pas eu le temps d'imprimer les transparents pour le cours. Désolé...
Tuesday, February 21 2012
Examen intra, modèles de régression
By arthur charpentier on Tuesday, February 21 2012, 23:32 - ACT6420-H2012
L'examen intra comportera deux exercices proches de ceux faits en démonstration, et une analyse de sorties de régression. Cette dernière partie portera sur la base suivante (ou peut être une extraction de cette base)
> examen=read.table( + "http://freakonometrics.blog.free.fr/public/
data/basket-exam-v2.csv", + header=TRUE,sep=";") > tail(examen) EQUIPE halfdiff winner 9917 MARYLAND -1 0 9918 RHODE ISLAND -15 1 9919 ROBERT MORRIS 4 1 9920 NORTH CAROLINA-ASHEVILLE 11 1 9921 BROWN -2 1 9922 MOUNT ST. MARYS 1 0 totalpointsdiff secondhalfpoints 9917 -10 26 9918 1 44 9919 6 43 9920 19 39 9921 6 48 9922 -9 26 nodiff_pts_1000_left hometeam halfwinner 9917 9 1 0 9918 18 1 0 9919 17 1 1 9920 20 1 1 9921 16 1 0 9922 13 1 1 halflooser secondhalfpointsother 9917 1 35 9918 1 28 9919 0 41 9920 0 31 9921 1 40 9922 0 36
qui contient des matchs de baskets universitaires.
- EQUIPE le nom de l'équipe qui sert de référence
- HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l'EQUIPE perdait)
- WINNER vaut 1 si l'équipe a gagné le match, 0 sinon,
- TOTALPOINTSDIFF le nombre de points de différence à la fin du match
- SECONDHALFPOINTS est le nombre de points marqués par l'EQUIPE pendant la seconde mi-temps
- NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
- HOMETEAM vaut 0 si l'équipe joue à l'extérieur, 1 si elle joue à domicile
- HALFWINNER vaut 1 si l'EQUIPE gagnait (supérieur ou nul) à la mi-temps
- HALFLOOSER vaut 1 si l'EQUIPE perdait (inférieur ou nul) à la mi-temps
- SECONDHALFPOINTSOTHER est le nombre de points marqués par l'adversaire au cours de la seconde mi-temps
Monday, February 6 2012
Qualité d'une régression
By arthur charpentier on Monday, February 6 2012, 21:47 - ACT6420-H2012
La semaine passée, on avait commencé à travailler sur les incendies dans les différents quartiers de Chicago,
> 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) + reg123=lm(Y~.,data=base) + summary(reg123) Call: lm(formula = Y ~ ., data = base) 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 + reg23=lm(Y~X2+X3,data=base) > summary(reg23) Call: lm(formula = Y ~ X2 + X3, data = base) Residuals: Min 1Q Median 3Q Max -9.789 -4.644 -1.472 3.743 15.983 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 21.49646 3.78470 5.680 9.95e-07 *** X2 0.22129 0.05731 3.862 0.000366 *** X3 -1.52477 0.31059 -4.909 1.30e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.453 on 44 degrees of freedom Multiple R-squared: 0.4415, Adjusted R-squared: 0.4161 F-statistic: 17.39 on 2 and 44 DF, p-value: 2.721e-06
Comme seules deux variables semblent représentatives, on pourrait être tenté de visualiser la régression,
> library(scatterplot3d) > s3d=scatterplot3d(X2,X3,Y,type="h")
On peut ensuite représenter le plan de régression,
> s3d$plane3d(reg23)
Sur la figure ci-dessous, on peut visualiser les observations, et leur projection sur le plan de régression (le trait étant le résidu)
L'interprétation des signes des paramètres se fait sur la figure suivante: on fixe une des variables (c'est la notion de "toutes choses étant égales par ailleurs") on regarde ce que donnerait une hausse ou une baisse de l'autre.
en fonction de 
,
(correspondant au 
ème résidus) versus 
. Si les résidus étaient gaussiens, les points devraient être alignés suivant la première bissectrice,> plot(reg23,which=2) > s=sd(E) > Es=sort(E/s) > n=nrow(base) > Qtheorique=qnorm((1:n-.5)/(n)) > points(Qtheorique,Es,col="blue")
Pour aller plus loin dans les diagnostics, rappelons que si on considère la matrice de projection 
, alors les résidus empiriques 
peuvent s'obtenir par la relation
de telle sorte que 
. On notera alors
> plot(reg23,which=3) > s=summary(reg23)$sigma > X=matrix(c(rep(1,n),X2,X3),n,3) > H=X%*%solve(t(X)%*%X)%*%t(X) > T=E/s/sqrt(1-diag(H)) > points(Yp,sqrt(abs(T)),col="blue")
Parmi les autres informations pertinentes, il y a la distance de Cook,
(on parle de distance de Cook car il s'agit d'une distance entre 
et 
, correspondant à l'estimation du paramètre en enlevant la ème observation),
On peut ensuite continuer les graphiques,
Thursday, February 2 2012
un mot (rapide) sur le R2 et la sélection de modèle
By arthur charpentier on Thursday, February 2 2012, 23:40 - ACT6420-H2012
J'évoquais tout à l'heure dans un rapide billet l'utilisation du 
dans une régression sans constante. On va donc évoquer ici uniquement le cas où une constante est prise en compte dans la régression pour parler de la sélection de modèle, et plus particulièrement le 
, souvent appelé -ajusté. La sélection de modèle sera l'objet du prochain cours. Classiquement, on défini le coefficient d'ajustement comme
Dans une régression avec une constante, le résidus empiriques sont centrés, et on peut interpréter les termes comme des variances, d'où l'écriture souvent utilisée
. Une autre propriété, dès lors que l'on a une constante dans le modèle est que le peut être relié à la corrélation (au sens de Pearson) entre la prédiction et l'observation,
C'est d'ailleurs sous cette forme qu'a été introduit, historiquement, ce coefficient (on peut renvoyer aux travaux de Karl Pearson, ou Udny Yule sur le sujet, lors des travaux sur l'hérédité à la fin du XIXème siècle). Mais ce coefficient, s'il pourrait être utilisé comme mesure de la qualité de l'ajustement, ne peut être utilisé pour faire de la sélection de modèle. En effet, en rajoutant des variables, quelles qu'elles soient, améliorera toujours le modèle. C'est ce que l'on arrive à faire ci-dessous en rajoutant juste du bruit dans la régression
> summary(lm(dist~speed,data=cars)) Call: lm(formula = dist ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -29.069 -9.525 -2.272 9.215 43.201 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.5791 6.7584 -2.601 0.0123 * speed 3.9324 0.4155 9.464 1.49e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.38 on 48 degrees of freedom Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438 F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12 > set.seed(1) > cars$bruit=rnorm(nrow(cars)) > summary(lm(dist~speed+bruit,data=cars)) Call: lm(formula = dist ~ speed + bruit, data = cars) Residuals: Min 1Q Median 3Q Max -31.869 -9.412 -3.049 7.410 43.769 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.905 6.757 -2.650 0.0109 * speed 3.935 0.415 9.482 1.73e-12 *** bruit 2.788 2.639 1.056 0.2963 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.36 on 47 degrees of freedom Multiple R-squared: 0.6592, Adjusted R-squared: 0.6447 F-statistic: 45.45 on 2 and 47 DF, p-value: 1.034e-11 > set.seed(1) > cars$bruit=matrix(rnorm(20*nrow(cars)),nrow(cars),20) > summary(lm(dist~.,data=cars)) Call: lm(formula = dist ~ ., data = cars) Residuals: Min 1Q Median 3Q Max -30.147 -5.883 -2.070 7.354 32.096 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -23.8799 8.4775 -2.817 0.00879 ** speed 4.3870 0.5740 7.643 2.52e-08 *** bruit1 -1.7487 4.1122 -0.425 0.67391 bruit2 -2.7356 3.0890 -0.886 0.38338 bruit3 4.2464 3.5200 1.206 0.23778 bruit4 -4.5160 3.2715 -1.380 0.17838 ... bruit19 1.5769 2.6615 0.592 0.55828 bruit20 1.5107 2.7799 0.543 0.59113 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.87 on 28 degrees of freedom Multiple R-squared: 0.7834, Adjusted R-squared: 0.621 F-statistic: 4.822 on 21 and 28 DF, p-value: 7.557e-05
Plus on rajoute des variables, plus on augmente le . Un ajustement a été proposé afin de viser un principe de parcimonie. L'ajustement en question peut se retrouver dans Fisher (1924),
même si beaucoup de monde semble l'attribuer à Theil (1961). L'idée derrière est que l'association faite entre les variances et les sommes des carrés est rapide, et peut-être fallacieuse, en pratique. En effet, on utilisait comme estimateurs des variances
et
(on peut avoir égalité si on a un modèle parfait). On notera aussi qui si l’ajustement est trop mauvais, en particulier
est négatif. Le est populaire car il permet de pénaliser le , en prenant en compte la complexité du modèle (ici le nombre de variables explicatives). Pour qu'un modèle soit meilleur avec une variable en plus, il faut que la variable que l'on rajoute apporte vraiment quelque chose. Pour comprendre plus concrètement ce qui se passe, considérons la régression sur 
et celle sur 
. On notera 
et 
les associés. Pesaran (1974) a noté qu'il existait une relation simple entre et 
est l'estimateur du coefficient ajouté dans la régression la plus complète, et 
la statistique de Student du test de significativité. Autrement dit, 
si et seulement si la nouvelle variable n'est pas suffisamment significative, en l’occurrence 
. Mais ce n'est pas le coefficient
(même avec l'ajustement) qui est le plus utilisé pour faire de la sélection de modèle, comme nous le verrons en cours la semaine prochaine (mais une pénalisation de la vraisemblance). A suivre donc...Monday, January 30 2012
calculs des sorties de régression
By arthur charpentier on Monday, January 30 2012, 14:05 - ACT6420-H2012
Avec un peu de calcul matriciel (et en utilisant les fonctions de répartition des lois usuelles) on peut calculer à la main toute les grandeurs obtenues lorsque l'on fait une régression linéaire. Utilisons la base de donnée avec le poids d'une personne (qui sera notre variable d'intérêt, et que l'on essayera d'inférer), la taille et le sexe de la personne,
> TP = read.table(
"http://freakonometrics.free.fr/taillepoids.txt",
header=TRUE,sep=";") > tail(TP) indice sexe poids1 taille1 poids2 taille2 195 195 F 62 164 61 161 196 196 M 74 175 71 175 197 197 M 83 180 80 180 198 198 M 81 175 NA NA 199 199 M 90 181 91 178 200 200 M 79 177 81 178 > Y=TP$poids1 > X1=TP$sexe > X2=TP$taille1 > base=data.frame(Y,X1,X2) > tail(base) Y X1 X2 195 62 F 164 196 74 M 175 197 83 M 180 198 81 M 175 199 90 M 181 200 79 M 177
La régression linéaire donne ici
> regression=lm(Y~X1+X2) > summary(regression) Call: lm(formula = Y ~ X1 + X2) Residuals: Min 1Q Median 3Q Max -24.718 -7.159 -1.472 5.487 74.726 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 109.11420 14.20513 7.681 7.22e-13 *** X1M 22.49801 2.08690 10.781 < 2e-16 *** X2 -0.31298 0.08649 -3.619 0.000376 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 11.81 on 197 degrees of freedom Multiple R-squared: 0.3937, Adjusted R-squared: 0.3875 F-statistic: 63.95 on 2 and 197 DF, p-value: < 2.2e-16
On peut commencer par (re)calculer les estimateurs par moindres carrés ordinaires, en notant que de la variable associée au sexe, seul le sexe masculin a été retenu ici.
> n=length(Y) > X=matrix(c(rep(1,n),(X1=="M")*1,X2),n,3) > X[1:5,] [,1] [,2] [,3] [1,] 1 1 182 [2,] 1 0 161 [3,] 1 0 161 [4,] 1 1 177 [5,] 1 0 157 > (beta=solve(t(X)%*%X)%*%t(X)%*%Y) [,1] [1,] 109.1141967 [2,] 22.4980107 [3,] -0.3129827 > summary(regression)$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 109.1141967 14.20512536 7.681326 7.222037e-13 X1M 22.4980107 2.08690147 10.780581 1.347214e-21 X2 -0.3129827 0.08648531 -3.618912 3.760466e-04
On retrouve ici est estimation indiquée. Notons que si la modalité de référence était le sexe féminin, nous aurions,
> Xbis=matrix(c(rep(1,n),(X1=="F")*1,X2),n,3) > Xbis[1:5,] [,1] [,2] [,3] [1,] 1 0 182 [2,] 1 1 161 [3,] 1 1 161 [4,] 1 0 177 [5,] 1 1 157 > solve(t(Xbis)%*% Xbis)%*%t(Xbis)%*%Y [,1] [1,] 131.6122073 [2,] -22.4980107 [3,] -0.3129827
Si on compare les deux constantes, entre ces deux modèles, on retrouve la valeur du coefficient associé au sexe de l'individu,
> (solve(t(X)%*%X)%*%t(X)%*%Y)[1]- + (solve(t(Xbis)%*% Xbis)%*%t(Xbis)%*%Y)[1] [1] -22.49801 > (solve(t(X)%*%X)%*%t(X)%*%Y)[2] [1] 22.49801
Continuons un peu... L'estimateur de la volatilité des résidus est ici
> residus=Y-X%*%beta > t(residus)%*%residus [,1] [1,] 27493.32 > sum(residus^2) [1] 27493.32 > > (sigma=sqrt(sum(residus^2)/(n-ncol(X)))) [1] 11.81355 > summary(regression)$sigma [1] 11.81355
On peut utiliser cette grandeur pour estimer la matrice de variance-covariance de notre estimateur,
> (varbeta=sigma^2*solve(t(X)%*%X)) [,1] [,2] [,3] [1,] 201.785586 16.2312631 -1.224735573 [2,] 16.231263 4.3551577 -0.106737634 [3,] -1.224736 -0.1067376 0.007479709 > vcov(regression) (Intercept) X1M X2 (Intercept) 201.785586 16.2312631 -1.224735573 X1M 16.231263 4.3551577 -0.106737634 X2 -1.224736 -0.1067376 0.007479709
Le terme diagonal est celui qui est utilisé pour calculer l'écart-type des estimateurs de chacun des coefficients.
> (sigbeta=sqrt(diag(varbeta))) [1] 14.20512536 2.08690147 0.08648531 > summary(regression)$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 109.1141967 14.20512536 7.681326 7.222037e-13 X1M 22.4980107 2.08690147 10.780581 1.347214e-21 X2 -0.3129827 0.08648531 -3.618912 3.760466e-04
En faisant une hypothèse de normalité des résidus (ou en supposant le nombre d'observations grand, i.e. suffisamment grand pour invoquer des théorèmes asymptotiques), on peut faire un test de significativité des variables explicatives,
> (beta/sigbeta) [,1] [1,] 7.681326 [2,] 10.780581 [3,] -3.618912 > 2*(1-pt(abs(beta/sigbeta),df=n-ncol(X))) [,1] [1,] 7.220891e-13 [2,] 0.000000e+00 [3,] 3.760466e-04
avec respectivement la statistique de Student, et la p-value associé au test de nullité. Pour juger de la qualité de l'ajustement, on peut calculer
> (R2=cor(Y,X%*%beta)^2) [,1] [1,] 0.3936724 > summary(regression)$r.squared [1] 0.3936724 > > (R2a=1-(n-1)/(n-(ncol(X)))*(1-R2)) [,1] [1,] 0.3875168 > summary(regression)$adj.r.squared [1] 0.3875168
et enfin, le test de significativité global (de Fisher) donne ici
> (F=R2/(ncol(X)-1)/((1-R2)/(n-ncol(X)))) [,1] [1,] 63.95343 > summary(regression)$fstatistic value numdf dendf 63.95343 2.00000 197.00000 > 1-pf(F,(ncol(X)-1),(n-ncol(X))) [,1] [1,] 0
Bref, avec un peu de calcul, on peut retrouver toutes les valeurs données sous R. Une fois les calculs compris, on pourra passer à la géométrie, et à l'interprétation...
Wednesday, January 18 2012
de la normalité des estimateurs dans une régression
By arthur charpentier on Wednesday, January 18 2012, 02:40 - ACT6420-H2012
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 
était un estimateur Gaussien. En particulier, chacun des estimateurs est alors Gaussien, au sens où 
, 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, 
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.
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..,
Friday, January 13 2012
ACT6420 Méthodes de prévision, introduction à R
By arthur charpentier on Friday, January 13 2012, 11:28 - ACT6420-H2012
Cet après midi, j'assurerais la démonstration, en salle informatique à 14h30, PK S1 590, au sous sol (salle L). On fera une introduction à R, afin de présenter les rudiments du langage. Je mets en ligne quelques lignes de code,
> X=cars$speed > Y=cars$dist > n=length(X) > (mx=mean(X)) [1] 15.4 > sum(X)/n [1] 15.4 > s=0 > for(i in 1:n){s=s+X[i]} > print(s/n) [1] 15.4 > > my=mean(Y) > (sx=sd(X)) [1] 5.287644 > sqrt(sum((X-mx)^2)/(n-1)) [1] 5.287644 > sy=sd(Y) > (r=cor(X,Y)) [1] 0.8068949 > sum((X-mx)*(Y-my))/((n-1)*sx*sy) [1] 0.8068949 > > X<=15 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [9] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [17] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [25] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE [33] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [41] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [49] FALSE FALSE > sum(X<=15) [1] 26 > X[X<=15] [1] 4 4 7 7 8 9 10 10 10 11 11 12 12 12 12 13 13 [18] 13 13 14 14 14 14 15 15 15 > mean(Y[X<=15]) [1] 27.5 > mean(Y[X>15]) [1] 59.75 > > (b1=r*sy/sx) [1] 3.932409 > (b0=my-b1*mx) [1] -17.57909
> plot(X,Y) > plot(X,Y,xlab="Vitesse de la voiture", + ylab="Distance de freinage",pch=4, + xlim=c(0,30),axes=FALSE) > axis(1);axis(2) > abline(a=b0,b=b1,col="red",lty=2) > Yh=b0+b1*X > points(X,Yh,pch=19,col="blue")
> (SCT=sum((Y-my)^2)) [1] 32538.98 > (SCR=sum((Y-Yh)^2)) [1] 11353.52 > (SCE=sum((Yh-my)^2)) [1] 21185.46 > SCR+SCE [1] 32538.98 > > (s=sqrt(sum((Y-Yh)^2)/(n-2))) [1] 15.37959 > > (sb1=s/(sx*sqrt(n-1))) [1] 0.4155128 > t=(b1-0)/sb1 > 2*(1-pt(t,df=n-2)) [1] 1.489919e-12 > > lm(Y~X) Call: lm(formula = Y ~ X) Coefficients: (Intercept) X -17.579 3.932 > regression=lm(dist~speed,data=cars) > names(regression) [1] "coefficients" "residuals" "effects" [4] "rank" "fitted.values" "assign" [7] "qr" "df.residual" "xlevels" [10] "call" "terms" "model" > regression$coefficient (Intercept) speed -17.579095 3.932409 > summary(regression) Call: lm(formula = dist ~ speed, data = cars) Residuals: Min 1Q Median 3Q Max -29.069 -9.525 -2.272 9.215 43.201 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.5791 6.7584 -2.601 0.0123 * speed 3.9324 0.4155 9.464 1.49e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 15.38 on 48 degrees of freedom Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438 F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12 > summary(regression)$sigma [1] 15.37959
Monday, October 10 2011
ACT2040, régression binomiale et arbres
By arthur charpentier on Monday, October 10 2011, 23:44 - ACT2040-2011
Mardi, avant de parler de la régression poissonienne, on va revenir un
peu sur les modèles binomiaux, et présenter les arbres de régression.
Pour la théorie, je peux renvoyer vers des notes de cours sur internet (ici,
là, ou là
encore) mais je présenterais les grandes idées au tableau. Pour le
code, c'est assez facile
Pour rappel, les trois critères sont données par les relations suivantes. Pour le critère du chi-deux,
où, classiquement
pour le découpage en deux classes, puis
pour le découpage en trois classes, etc. Pour l'entropie,
avec deux classes, ou, si on a trois classes,
Enfin le critère de Gini,
avec deux classes, ou avec trois,
> sinistre=read.table("http://freakonometrics.free.fr/sinistreACT2040.txt",
+ header=TRUE,sep=";")
> sinistres=sinistre[sinistre$garantie=="1RC",]
> contrat=read.table("http://freakonometrics.free.fr/contractACT2040.txt",
+ header=TRUE,sep=";")
> T=table(sinistres$nocontrat)
> T1=as.numeric(names(T))
> T2=as.numeric(T)
> nombre1 = data.frame(nocontrat=T1,nbre=T2)
> I = contrat$nocontrat%in%T1
> T1= contrat$nocontrat[I==FALSE]
> nombre2 = data.frame(nocontrat=T1,nbre=0)
> nombre=rbind(nombre1,nombre2)
> base = merge(contrat,nombre)
> Y = base$nbre>0
> X1 = base$ageconducteur
> library(tree)
> arbre=tree(Y~X1,split="gini")
> plot(arbre)
> text(arbre,cex=.8)
> arbre=tree(Y~X1,split="gini", mincut = 5000)
> plot(arbre)
> text(arbre,cex=.8)
> age=base$ageconducteur
> car=base$agevehicule
> region=base$zone
> arbre=tree(Y~age+car+region,split="gini", mincut = 2600)
> plot(arbre)
> text(arbre,cex=.8)
> u=sort(unique(base$ageconducteur))
> base$Y = base$nbre>0
> gini=entropie=chideux=rep(NA,length(u))
> for(i in 2:(length(u))){
+ I=base$ageconducteur<(u[i]-.5)
+ ft1=sum(base$Y[I==TRUE])/nrow(base)
+ ff1=sum(base$Y[I==FALSE])/nrow(base)
+ ft0=sum(1-base$Y[I==TRUE])/nrow(base)
+ ff0=sum(1-base$Y[I==FALSE])/nrow(base)
+ M=matrix(c(ft0,ff0,ft1,ff1),2,2)
+ ft=ft0+ft1
+ f0=ft0+ff0
+ Mind=matrix(c(ft*f0,f0*(1-ft),(1-f0)*ft,(1-f0)*(1-ft)),2,2)
+ Q=sum(nrow(B)*(M-Mind)^2/Mind)
+ entropie[i] = -((ft1+ft0)*(ft1/(ft1+ft0)*log(ft1/(ft1+ft0))+
+ ft0/(ft1+ft0)*log(ft0/(ft1+ft0))) +
+ (ff1+ff0)*(ff1/(ff1+ff0)*log(ff1/(ff1+ff0))+
+ ff0/(ff1+ff0)*log(ff0/(ff1+ff0)))
+ )
+
+ gini[i] = -((ft1+ft0)*(ft1/(ft1+ft0)*(1-ft1/(ft1+ft0))+
+ ft0/(ft1+ft0)*(1-ft0/(ft1+ft0))) +
+ (ff1+ff0)*(ff1/(ff1+ff0)*(1-ff1/(ff1+ff0))+
+ ff0/(ff1+ff0)*(1-ff0/(ff1+ff0))))
+
+ chideux[i] = Q
+ }
> plot(u-.5,entropie,type="b",pch=19,col="blue",
+ ylab="Entropie ou déviance",xlab="Age du conducteur")
> n=which.min(entropie)
> points(u[n]-.5,entropie[n],pch=19,col="red",cex=1.5)
Tuesday, October 19 2010
On a toujours besoin de petits poids (chez soi)
By arthur charpentier on Tuesday, October 19 2010, 10:15 - Statistics
Après avoir mis en ligne le billet d'hier sur les manifestants, j'ai passé la journée à me demander "pourquoi ai-je mis des poids tenant compte d'un effet taille dans ma régression" ?. Sur le moment ça me semblait très naturel, et plus la journée avançait, plus je me mettais à douter...
Bon, pour rappel quand on fait une régression classique, sans poids, et sans constante car on veut un coefficient de proportionnalité, par "moindres carrés", on résout
ce qui donne
La dernière écriture montre que l'on fait une moyenne pondérée des pentes pour chacune des observations... Graphiquement, ça donne
les traits pointillés permettent de montrer ou se trouve le point moyen, avec la moyenne sur toutes les villes pour les syndicats et pour la police, pour les deux dates. Les flèches indiquent les pentes moyennes, c'est à dire le nombre total de manifestants en France selon les syndicats divisé par le nombre total de manifestants en France selon la police. Bref, la droite de régression "classique" n'a pas pour pente le ratio moyen...
Une alternative usuelle est d'introduire des poids, et on cherche à résoudre
dont la solution a la forme
Hier j'avais proposé de prendre des poids correspondant à la taille de la ville (pour intégrer un effet médiatique - les médias ne parlant parfois que des nombres de manifestants dans les très grosses villes), mais on pourrait proposer, sans le même genre d'idée que de prendre 
(et qui est représenté ci-dessous)
Toutefois, l'idée la plus pertinente aurait dû être 
ce qui donne comme estimateur par moindres carrés pondérés
c'est à dire le ratio moyen national... Bref, on est colinéaire avec nos vecteurs tracés depuis le début
Autrement dit, il aurait été plus malin ici d'attribuer des petits poids aux grosses villes pour avoir un coefficient interprétable....
Mais au débat du débat (probablement légitime) sur le choix de la régression, ce qui est troublant c'est que plus les poids donnés aux grosses villes sont importants, plus la pente est importante... autrement dit, on modèle quadratique pourrait être plus adapté, ou au moins quelque chose de non linéaire, et de strictement croissant. Sur le graphique suivant, on voit l'estimation par moindres carrées locaux, et par splines,
ce qui pourrait nous conforter dans l'idée que "les grandes villes ont des plus grands trottoirs"(ramené à la taille de la chaussée !), ce qui peut surprendre un peu.. Si on regarde l'estimation des pentes obtenues par ces deux méthodes, on a, pour les régressions locales
(avec en trait horizontal en pointillés les pentes des flèches, i.e. les taux observés nationalement), et par splines
On retrouve l'idée que dans les petites villes, le comptage est plus précis (des deux côtés) avec un facteur de 2 et 3 pour les deux manifestations respectivement. En mars dernier, on avait ensuite un facteur croissant passant de 2 à 4 en fonction de la taille de la manifestation. Et pour la dernière manifestation en revanche, c'est un peu n'importe quoi (si je peux me permettre).
Bref, pour revenir au point qui m'intéressait au départ, on note qu'il ne semble pas y avoir de méthode scientifique permettant de passer du nombre estimé par la police à celui proposé par les syndicats, car sinon il y aurait une relative stabilité entre les manifestations... stabilité que l'on n'a pas (quelle que soit le type de régression considéré...).
Monday, October 18 2010
Qui a eu l'idée d'élargir les trottoirs en France ?
By arthur charpentier on Monday, October 18 2010, 11:28 - Statistics
J'étais tombé l'autre jour sur un joli dessin expliquant pourquoi les policiers et les syndicats n'estimaient jamais le même nombre de manifestants, avec une explication scientifique, et vaguement rationnelle, ici.
Si on compare le 19 mars dernier (ici) et le 16 octobre (là), on arrive aux chiffres suivants,
Thursday, June 3 2010
Qu'est ce que la modélisation économétrique (3) ?
By arthur charpentier on Thursday, June 3 2010, 15:00 - Statistics
Compte tenu du nombre de billets qui agitent la blogosphère en ce moment, je me suis senti obligé de reprendre la plume pour poursuivre le billet précédant (ici). Initialement, le but était de montrer aux élèves comment mener une étude économétrique et de répondre à une question simple (en l'occurence "la taille d'un classe a-t-elle un impact positif ou négatif sur les résultats scolaires ? "). Malheureusement, j'ai été rattrapé par l'actualité. Le premier billet montrait qu'en faisant un modèle de régression simple, l'effet semblait positif: plus la taille est grande, meilleure est la moyenne. Dans le second billet, nous avions noté que cet effet pouvait cacher quelquechose, comme des conditions socioéconomiques. Essayons donc d'aller un peu plus loin....
- Régression multiple
> summary(lm(avgverb~tipuach+classize,data=base0))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 80.26363 0.74472 107.777 <2e-16 ***
tipuach -0.34994 0.01077 -32.486 <2e-16 ***
classize -0.03146 0.02221 -1.416 0.157
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
L'impact de l'indice socio-économique est clairement négatif, et si l'effet taille a un impact négatif, ce dernier n'est pas significatif. Il va falloir essayer d'aller un peu plus loin....
- Approche nonlinéaire
> library(mgcv)
> reg=gam(avgverb~s(tipuach,classize),data=base0)
> reg=gam(avgmath~s(tipuach,classize),data=base0)
soit visuellement, les prédictions suivantes pour l'épreuve de lecture
- Prise en compte de la règle de Maimonides
En fait, si on défini le nombre théorique d'élèves dans les classes, avec cette méthode (ce que j'avais détaillé dans la section 3 ici), on obtient des résultats qui pourrait ressemblerà ce que nous avions observé sur le nombre réel d'élèves dans les classes,
> reg=lm(avgverb~func1,data=base0)
> summary(reg)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 70.62076 0.88004 80.25 < 2e-16 ***
func1 0.12159 0.02789 4.36 1.37e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
On retrouve un effet positif (significatif) qui se visualise sur le graphique ci-dessous,
> reg=lm(avgverb~func1+tipuach,data=base0)
> summary(reg)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.91057 0.78669 105.392 < 2e-16 ***
func1 -0.11160 0.02321 -4.809 1.63e-06 ***
tipuach -0.35943 0.01050 -34.229 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
où cette fois l'impact de la taille de la classe est clairement significatif, et négatif: à charactéristiques socio-économique identiques, si la taille des classes suivait la règle de Maimonides, la taille de la classe aurait un impact négatif sur les résultats scolaires....
Une solution est d'utiliser une régression par variables instrumentales où l'instrument permet précisément de prendre en compte cette discontinuité. L'idée est ici qu'il existe une variable
telle que
Damned, l'étude n'est donc toujours pas finie, il va falloir creuser encore davantage (surtout qu'il reste un paquet de variables dans cette base)... mais pour les impatients, la conclusion à laquelle arrivent Joshua Angrist et Victor Lavy est résumé dans le tableau ci-dessous,
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 - econometrics
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
- Partie 1: la problématique
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)
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
par école (school) et par classe (class). On note aussi
désigne le nombre d'élèves dans l'école
pour chacun des niveaux, 
un identifiant pour la ville, 
une
variable indiquant le niveau socio-économique des
élèves, et 
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,
(1)
désigne la partie entière. - Partie 4: l'étude économétrique, les moindres carrés
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,
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
on va mettre tout ce qui, dans les aspects socio-économiques -
pourrait expliquer les résultats aux tests. Le terme du milieu 
correspond à ce qui nous intéresse, à savoir
l'impact de la taille des classes les tests. Enfin, le terme de droite 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")
> 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")
> 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
Wednesday, October 7 2009
J'ai de la multicolinéarité... c'est grave docteur ?
By arthur charpentier on Wednesday, October 7 2009, 12:33 - économétrie 1 - M1 - 09/10
Un petit billet sur la multicolinéarité (que j'avais un peu abordé ici par exemple), pour répondre à la question de savoir si c'est "grave" (je reprends la phrase d'un mail que j'ai reçu récemment).
- Ca existe la multicolinéarité ou bien on ne voit ça qu'à la télé ?
Par exemple (je reprends un exemple que je connais bien) en tarification des contrats d'assurance automobile, on peut avoir plusieurs informations à notre disposition, par exemple l'âge du conducteur, et l'ancienneté du permis de conduire... Ces données ne sont pas colinéaires (ou comonotone) car il existe des personnes qui passent leur permis de conduire tard, mais on se doute qu'elles sont très corrélées.... Voici un petit exemple sur de vraies données (même si on ne voit pas forcément grand chose à priori car les données sont discrètes), avec en abscisse l'ancienneté du permis et en ordonnée l'âge du conducteur
- Formulation du problème pour apporter des éléments de réponse
On sait que la variable Y dépend de la variable
via un modèle que l'on supposera linéaire (pour faire
simple, mais ça ne change pas grand chose sur le fond). Mais
dans la base, on dispose de deux variables explicatives très
corrélées (voire plus) et 
. La question telle que je la comprends est la suivante: même si
les variables sont très corrélées, est-ce que l'on
est en mesure de ne garder que la bonne variable (en l'occurence ). On peut simuler le modèle suivant,
est un bruit indépendant des variables explicatives, suivant une loi normale centrée et réduite.C'est le vrai modèle. Mais on cherche à estimer un modèle de la forme
- Tester la significativité du second coefficient
Sur la figure ci-dessous, on regarde la distribution de la p-value du test de significativité de Student sur la seconde variable (la distribution est ici estimée par noyau, avec tous les problèmes techniques que cela pose1), avec une corrélation de 0,25 (en rouge), de 0,5 (en mauve) et de 0,75 (en bleu),
- Utiliser des outils de sélection de modèle
Si l'on compare les R2 des deux modèles, on est rassuré car le R2 est presque toujours2 plus grand pour le premier (le bon) que pour le second, commen en témoigne la distribution de la différence entre les deux, sur autant de simulations que tout à l'heure, mais plus la corrélation est grande, plus la différence entre les R2 diminue....
2 sur 10,000 simulations, je n'ai obtenu aucun cas pour lequel le R2 du second modèle excédait le R2 du premier si la corrélation était inférieur à 0,99, mais avec une corrélation de 0,999, cela pouvait survenir... Maintenant, des corrélations de cet ordre me paraissent surprenantes en pratique (au sens où les variables sont quasiment identiques)...
