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) |
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 |
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 |
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 |
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 *** |
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) |
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
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 |
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 |
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
(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”).
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 |
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 |
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 *** |
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 |
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
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 |
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 |
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.