Avant de parler davantage du premier devoir, un mot pour expliquer pourquoi dans la première partie il est important d’éviter d’utiliser des données temporelles (le cours présente des outils pour modéliser des données individuelles). Considérons un exemple simple (évoqué voilà plusieurs mois ici même), où on cherche à prédire un tour de poitrine conditionnellement à un tour de taille (pour une personne de sexe féminin, cf. le mannequin ci-contre). L’idée intuitive est qu’il doit exister une relation positive, possiblement linéaire entre les deux mensurations. A partir de données, on devrait pouvoir quantifier cette relation. On dispose pour cela de données, en ligne sur

> mensurations<-read.table("
+ http://freakonometrics.free.fr/mensurations.csv",
+ header=TRUE,sep=";")
> head(mensurations[,-1])
POITRINE TAILLE HANCHE HAUTEUR POIDS      AGE
1       NA     NA     NA      NA    NA       NA
2       94     58     89     165    51 22.24230
3       NA     NA     NA      NA    NA       NA
4       94     58     89     165    51 22.40383
5       91     58     97     165    55       NA
6       NA     NA     NA      NA    NA       NA

Si on regarde rapidement la corrélation entre la taille et la poitrine, cette dernière est faible,

> Z=mensurations[,c("TAILLE","POITRINE")]
> I=apply(Z,1,function(x) sum(is.na(x)))
> cor(Z[I==0,])
TAILLE   POITRINE
TAILLE   1.00000000 0.09700002
POITRINE 0.09700002 1.00000000

… mais significativement non nulle,

> cor.test(Z[,1],Z[,2],,method = "pearson")

Pearsons product-moment correlation

data:  Z[, 1] and Z[, 2]
t = 2.4559, df = 635, p-value = 0.01432
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.0194634 0.1733769
sample estimates:
cor
0.09700002

Avec une -value de l’ordre de 1.5%, on va rejeter l’hypothèse  de nullité de la corrélation (c’est ici un test basé sur la transformation de Fisher). On peut aussi faire une régression simple,

> reg=lm(POITRINE~TAILLE,data=mensurations)
> summary(reg)

Call:
lm(formula = POITRINE ~ TAILLE, data = mensurations)

Residuals:
Min      1Q  Median      3Q     Max
-9.7253 -3.4424 -0.2091  1.5576 13.8579

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 82.67798    2.82953  29.220   <2e-16 ***
TAILLE       0.11663    0.04749   2.456   0.0143 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.133 on 635 degrees of freedom
(22 observations deleted due to missingness)
Multiple R-squared: 0.009409,	Adjusted R-squared: 0.007849
F-statistic: 6.031 on 1 and 635 DF,  p-value: 0.01432

(on notera que le test sur la corrélation, le test de significativité de la – seule – variable explicative, et le test de significativité globale donnent ici exactement la même -value de l’ordre de 1.5%). Là encore, on voit que le tour de taille a un effet significativement non nul… mais marginal. Graphiquement, la prédiction avec le tour de taille est le trait rouge,

La valeur moyenne est en trait pointillé. On peut voir que la différence est infime,

> base=data.frame(TAILLE=seq(50,70,by=5))
> base$Pred.POITRINE=predict(reg,newdata=base)
> base$Moy.POITRINE=mean(POITRINE,na.rm=TRUE)
> base
TAILLE Pred.POITRINE Moy.POITRINE
1     50      88.50935     89.61538
2     55      89.09248     89.61538
3     60      89.67562     89.61538
4     65      90.25876     89.61538
5     70      90.84189     89.61538

Le soucis est que cette base est en fait une base sur données individuelles. En fait, il s’agit des mensurations des demoiselles qui ont posé comme playmate dans la version U.S. du magasine Playboy. Si on se restreint sur des sous-période, par exemple les années 70, on a des prédictions relativement différentes,

> submensurations=mensurations[abs(DATE-1975)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> summary(subreg)

Call:
lm(formula = POITRINE ~ TAILLE, data = submensurations)

Residuals:
Min     1Q Median     3Q    Max
-8.655 -1.871 -0.655  1.750 13.155

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  66.1485     7.3839   8.958 6.12e-15 ***
TAILLE        0.4053     0.1230   3.296   0.0013 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.574 on 117 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.08496,	Adjusted R-squared: 0.07714
F-statistic: 10.86 on 1 and 117 DF,  p-value: 0.001299

> base$Pred70s.POITRINE=predict(subreg,newdata=base)
> base
TAILLE Pred.POITRINE Moy.POITRINE Pred70s.POITRINE
1     50      88.50935     89.61538         86.41274
2     55      89.09248     89.61538         88.43917
3     60      89.67562     89.61538         90.46559
4     65      90.25876     89.61538         92.49201
5     70      90.84189     89.61538         94.51844

que l’on peut visualiser ci-dessous,

L’explication est que la morphologie des femmes (ou en tous les cas de celles qui ont posé nues dans Playboy) a beaucoup changé en 60 ans, comme évoquéprécédemment. On peut le voir en regardant l’évolution de la pente, par date (5ans)

(le trait horizontal est la pente obtenue sur l’ensemble de la période). On peut d’ailleurs intégrer l’année dans notre modèle, en rajoutant le temps comme variable explicative dans la régression,

> reg=lm(POITRINE~TAILLE+DATE,data=mensurations)
> summary(reg)

Call:
lm(formula = POITRINE ~ TAILLE + DATE, data = mensurations)

Residuals:
Min      1Q  Median      3Q     Max
-9.7304 -2.2184 -0.5405  2.0742 12.1041

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 368.296266  18.008455  20.451  < 2e-16 ***
TAILLE        0.323319   0.042140   7.673 6.39e-14 ***
DATE         -0.150301   0.009393 -16.002  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.491 on 634 degrees of freedom
(22 observations deleted due to missingness)
Multiple R-squared: 0.2944,	Adjusted R-squared: 0.2922
F-statistic: 132.3 on 2 and 634 DF,  p-value: < 2.2e-16

Visuellement, on aurait la prédiction suivante,

Si on se place autour de 1970,

> submensurations=mensurations[abs(DATE-1970)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> predict(subreg,newdata=data.frame(TAILLE=c(55,70)))
1        2
89.14361 95.75476

alors qu’autour de 1990,

> submensurations=mensurations[abs(DATE-1990)<=5,]
> subreg=lm(POITRINE~TAILLE,data=submensurations)
> predict(subreg,newdata=data.frame(TAILLE=c(55,70)))
1        2
87.69853 92.70228

Notons qu’on fait ici des régressions locales (dans ces deux derniers calculs), ce qui n’est pas ce qui est représenté visuellement… Si on veut une meilleure modélisation locale, il faudrait quitter le monde des modèles linéaires pour aller un peu plus loin. On verra en fin de cours si on a le temps d’en parler

Cela dit, globalement, on avait un effet positif, mais il aurait été tout à fait possible d’avoir un effet (global négatif). C’est ce qui arrive dans le paradoxe de Simpson.