Je vais reprendre ici le code tapé au tableau ce matin, mais sur une autre base de données. On va revenir sur des données qui ont donné le nom de méthodes de "régression" aux techniques que nous abordons en cours. Pour rappel, le mot régression apparaît dans le titre de l'article de Galton (1886) "regression towards mediocrity in hereditary stature". Mais on va plutôt travailler sur les données collectées dans Pearson (1903) sur la taille des pères et de leurs fils, à l'âge adulte (pour les personnes qui voudraient travailler sur la taille des mères et de leurs filles, je renvoie à Pearson & Lee (1903)).
> father.son=read.table( + "http://stat-www.berkeley.edu/users/juliab/141C/pearson.dat", + sep=" ") > attach(father.son)
On va chercher à prévoir la taille d'un fils, conditionnellement à la taille de son père,
> Y=father.son$V3 > X1=father.son$V2
Les estimateurs par moindres carrés sont obtenus en estimant d'abord la pente (à partir de la corrélation entre les deux variables)
> (r=cor(X1,Y)) [1] 0.5013383 > (sx=sd(X1)) [1] 2.744868 > (sy=sd(Y)) [1] 2.814702 > (b1=r*sy/sx) [1] 0.514093
et ensuite, on fixe la constante de manière à passer par le barycentre du nuage de points,
> (b0=mean(Y)-b1*mean(X1)) [1] 33.8866
Visuellement, on obtient
> plot(X1,Y) > abline(a=b0,b=b1,col="red") > abline(h=mean(Y),lty=2) > abline(v=mean(X1),lty=2)
(le barycentre étant ici à l'intersection des lignes horizontale et verticale, la droite bleue étant la première bissectrice). On retrouve des deux valeurs en faisant la régression sous R
> lm(Y~X1) Call: lm(formula = Y ~ X1) Coefficients: (Intercept) X1 33.8866 0.5141
Maintenant, on peut aussi faire un peu de calcul matriciel, pour estimer la droite de régression
> X=cbind(1,X1) > (B=solve(t(X)%*%X)%*%t(X)%*%Y) [,1] 33.886604 X1 0.514093
mais aussi la variance de la série des résidus,
> n=length(Y) > E=Y-X%*%B > (s2=sum(E^2)/(n-2)) [1] 5.936804
On peut alors aller un peu plus loin, en estimant la variance de nos estimateurs obtenus par moindres carrés,
> VB=s2*solve(t(X)%*%X) > sqrt(diag(VB)) X1 1.83235382 0.02704874
On retrouve ces valeurs sur la sortie détaillée de la régression,
> reg=lm(Y~X1) > summary(reg) Call: lm(formula = Y ~ X1) Residuals: Min 1Q Median 3Q Max -8.8772 -1.5144 -0.0079 1.6285 8.9685 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 33.88660 1.83235 18.49 <2e-16 *** X1 0.51409 0.02705 19.01 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.437 on 1076 degrees of freedom Multiple R-squared: 0.2513, Adjusted R-squared: 0.2506 F-statistic: 361.2 on 1 and 1076 DF, p-value: < 2.2e-16
Pour retrouver la valeur du ,
> sum((X%*%B-mean(Y))^2)/ + sum((Y-mean(Y))^2) [1] 0.2513401
Si on regarde la statistique du test de Student de significativité de la variable explicative (i.e. test pour savoir si la pente de la régression pourrait être nulle, ou pas)
> (T=(B[2]-0)/sqrt(VB[2,2])) [1] 19.00618
Tomber sur une telle valeur avec une loi de Student est rare, très rare... Pour s'en convaincre, on peut simuler quelques valeurs pour voir ce que peut être un tirage suivant une telle loi,
> rt(50,n-2) [1] 1.08379425 -0.17139018 0.36563361 -0.11331197 [5] -1.43101846 1.36197766 0.62777125 -0.86627679 [9] -0.69083285 1.40111081 0.19765215 -0.29186439 [13] 1.41921424 -0.38485241 -0.90524764 -0.58007929 [17] 0.20974133 0.93225320 1.30109930 0.67961650 [21] 0.45066183 0.19906039 -0.75011473 -1.50135529 [25] -0.20852248 -0.81572301 -0.38352095 -2.06752958 [29] -0.39492703 1.35945460 -0.36081938 -0.56273880 [33] -0.57315938 -0.87391435 -0.30124439 -0.55463941 [37] -0.08113229 -0.35297047 -0.81998094 -1.11087831 [41] 0.71881538 -0.08276505 1.14767768 1.85126104 [45] 0.14713420 -0.03935444 -0.92559603 -0.80283654 [49] -0.10799174 -2.16898630
Autrement dit, on n'atteint jamais de telles valeurs. On peut aussi regarder la densité de cette loi, pour voir qu'en moyenne, on est plutôt autour de 0, entre -2 et +2 dans 95% des cas,
> u=seq(-5,5,by=.1) > plot(u,dt(u,n-2),type="l") > lines(u,dnorm(u),col="red")
Si on essaye de calculer la -value, on voit qu'il est improbable d'obtenir une telle valeur avec une loi de Student, donc on va rejeter l'hypothèse que la pente puisse être nulle,
> 2*pt(-T,n-2) [1] 1.121268e-69
Cela dit, si on testait si la pente pouvait être 1/2, on aurait un résultat totalement différent,
> (T=(B[2]-.5)/sqrt(VB[2,2])) [1] 0.5210239 > 2*pt(-T,n-2) [1] 0.6024573
car avec une -value de l'ordre de 60%, on va accepter l'hypothèse de pente valant 1/2.
- we get
or 
:
or equivalently 
.
sera de la forme
