Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.
Davis=read.table( "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") Davis[12,c(2,3)]=Davis[12,c(3,2)] Davis=data.frame(Y=Davis$weight * 2.204622, X1=Davis$sex) |
On voulait estimer le modèle y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_iOn avait vu que l’on pouvait passer par l’écriture matricielle
X=cbind(Davis$X1=='F',Davis$X1=='M') Y=Davis$Y |
car la matrice \mathbf{X}^T\mathbf{X} est inversible (une fois que l’on enlève la constante)
solve(t(X)%*%X) [,1] [,2] [1,] 0.008928571 0.00000000 [2,] 0.000000000 0.01136364 |
et donc l’estimateur par moindres carrés est (classiquement)\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}
solve(t(X)%*%X) %*% (t(X)%*%Y) [,1] [1,] 125.4272 [2,] 167.3258 |
ce qui correspond effectivement à la sortie de R,
reg=lm(Y~0+X1,data=Davis) summary(reg) Coefficients: Estimate Std. Error t value Pr(>|t|) X1F 125.427 1.960 64.00 <2e-16 *** X1M 167.326 2.211 75.68 <2e-16 *** |
Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes
x=Y[X[,1]==1] y=Y[X[,2]==1] nx=length(x) ny=length(y) |
On avait vu en cours que les \widehat{\mathbf{\beta}} avaient une interprétation très simple, puisque\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_iautrement dit \widehat{{\beta}}_M est le poids moyen des hommes. Et en effet
mean(y) [1] 167.3258 |
C’est finalement très naturel, ou intuitif.
On peut maintenant s’interroger sur l’écart-type de l’estimateur de \widehat{{\beta}}_M . Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici
sqrt(var(y)/ny) [1] 2.794391 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny) [1] 2.794391 |
car pour rappel\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de \mathbf{\beta} est proportionnel à \sigma^2 , la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de \sigma^2
s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2)) sqrt(s2/ny) [1] 2.210863 |
et en effet, on retombe sur la valeur donnée dans le tableau de régresion
sqrt(s2/nx) [1] 1.959721 |
(pareil pour l’autre coefficient).
On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).
X=cbind(1,Davis$X1=='M') |
Là encore, le modèle devient identifiable, et on obtient ici
solve(t(X)%*%X) %*% (t(X)%*%Y) [,1] [1,] 125.42724 [2,] 41.89855 |
On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence
mean(y)-mean(x) [1] 41.89855 |
La sortie de régression devient ici
reg2=lm(Y~X1,data=Davis) summary(reg2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 125.427 1.960 64.00 <2e-16 *** X1M 41.899 2.954 14.18 <2e-16 *** |
Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)
t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE) Two Sample t-test data: Y[X[, 1] == 1] and Y[X[, 2] == 1] t = -6.4475, df = 286, p-value = 4.826e-10 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -30.62603 -16.30035 sample estimates: mean of x mean of y 143.8626 167.3258 |
Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable \mathbf{y} ), ça marche ! Posons
Davis$Y=sample(Davis$Y) |
ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici
reg2=lm(Y~X1,data=Davis) summary(reg2) Call: lm(formula = Y ~ X1, data = Davis) Residuals: Min 1Q Median 3Q Max -57.458 -22.184 -5.512 17.809 118.912 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 143.4382 2.7820 51.56 <2e-16 *** X1M 0.9645 4.1940 0.23 0.818 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 29.44 on 198 degrees of freedom Multiple R-squared: 0.000267, Adjusted R-squared: -0.004782 F-statistic: 0.05289 on 1 and 198 DF, p-value: 0.8183 |
autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient
Y=Davis$Y t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE) Two Sample t-test data: Y[X[, 1] == 1] and Y[X[, 2] == 1] t = -0.22998, df = 198, p-value = 0.8183 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -9.235209 7.306165 sample estimates: mean of x mean of y 143.4382 144.4027 |
et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents.
