---
title: TP chapitre 5
author: Antoine Lejay
output: html_document
---
```{r}
data(LifeCycleSavings)
lcs <- LifeCycleSavings
help(LifeCycleSavings)
head(lcs)
```
Considérons un modèle
$y_i=\beta_0 + \beta_1 x_{i,1} + \dotsb + \beta_{k} x_{i,k} +\varepsilon_i$.
Nous créons 2 fonctions, l'un qui donne la somme des
carrés résiduels (`ssr`) et l'autre la somme
des carrés expliqués (`sse`), c'est-à-dire
$$
SSR = \|\mathbf{y}-\widehat{\mathbf{y}}\|^2
\text{ et }
SSE = \|\widehat{\mathbf{y}}-\bar{y}\|^2.
$$
La moyenne des échantillons est la même que la moyenne
de la prédiction.
Nous définissons aussi la somme des carrés totales
$$
SST = \|\mathbf{y}-\bar{y}\|^2,
$$
c'est-à-dire la norme des observations centrées.
Rappelons que par le théorème de Pythagore,
$$
\sum_{i=1}^n (y_i-\bar{y})^2
=
\sum_{i=1}^n (y_i-\widehat{y}_i)^2
+
\sum_{i=1}^n (\widehat{y}_i-\bar{y})^2
$$
soit
$$ SST = SSE + SSR $$
```{r}
ssr <- function(reg)
{
return(sum(residuals(reg)^2))
}
sse <- function(reg)
{
return(sum((fitted(reg)-mean(fitted(reg)))^2))
}
sst <- function(reg)
{
return(sum((reg$model[,1]-mean(reg$model[,1]))^2))
}
```
Étudions cela sur un modèle à une variable :
```{r}
reg1 <- lm( sr ~ pop15, data = lcs)
cat(sprintf("SSR = %.3g\tSSE = %.3g\tSSR + SSE = %.3g\tSST = %.3g",ssr(reg1),sse(reg1),ssr(reg1)+sse(reg1),sst(reg1)))
```
## La commande `aov` pour une régression avec une variable
Ces calculs sont fournis par la commande `aov` :
```{r}
aov(reg1)
```
L'identification des _degrés de libertés_ est important, car nous
savons que les lois de $SSR$, $SST$ et $SSE$ sont de type $\sigma^2\chi^2_{f}$,
où $f$ est le nombre de degrés de libertés et $\sigma^2$ est la variance
du bruit (inconnue).
La commande `summary` appliquée à `aov` donne une information supplémentaire :
```{r}
summary(aov(reg1))
```
Les valeurs dans la colonne `Mean Sq` correspondent à la somme des carrés
divisée par le nombre de degrés de liberté.
##
La statistique $F$ est donc
$$
F = \frac{SSE/1}{SSR/(n-2)}.
$$
Ici, SSE n'a qu'une degré de libertés, et les résidus $n-2$
car 2 paramètres sont identifiés et que SSE calcule la moyenne.
Une conséquence du théorème de Cochran est qu'elle suit
sous l'hypothèse nulle « $\beta_1=0$ » une loi de Fisher
à $(1,n-2)$ degrés de libertés.
Cette statistique, retournée par `aov`, suit correspond aussi
à celle donnée par `summary(lm(...))`.
Dans le cas de la régression simple, la $p$-value
est la même que la statistique $T$ sur « $\beta_1=0$ ».
## Régression multiple et modèles emboîtés
Nous pouvons utiliser l'idée de Fisher sur des modèles
de régression multiple avec $\beta_1,\dotsc,\beta_k$
coefficients.
Considérons l'hypothèse nulle $H_0$ « $\beta_{q+1} = \dotsb = \beta_k =0$ »,
c'est-à-dire que seuls les $q$ premiers coefficients sont non nuls.
Sous l'hypothèse $H_0$, nous désignerons le modèle par _modèle réduit_.
Sans cette hypothèse, ce sera le _modèle plein_.
La statistique F consiste alors en
$$
F=\frac{SSE(\text{plein})-SSE(\text{réduit})}{k-q}\times
\frac{n-k}{SSR(\text{plein})}.
$$
Ici $k-q$ est le nombre de degrés de liberté de diffère
entre les deux modèles, et $n-k$ le nombre de degrés
de libertés de résidus sous le modèle plein.
Lorsque $k=1$ et $q=0$, alors $SSE(\text{réduit})=0$
et donc nous avons la même statistique que précédemment.
Notons que $SST(\text{plein})=SST(\text{réduit})$ et que
$$
SSR(\text{plein})+SSE(\text{plein}) = SST(\text{plein})
\text{ et }
SSR(\text{réduit})+SSE(\text{réduit}) = SST(\text{réduit})
$$
Ainsi
$$
SSE(\text{plein})-SSE(\text{réduit})
=
SSR(\text{réduit})-SSR(\text{plein}).
$$
Donc la statistique F est aussi
$$
F=\frac{SSR(\text{réduit})-SSR(\text{plein})}{k-q}\times
\frac{n-k}{SSR(\text{plein})}.
$$
Clairement, $SSR(\text{réduit}) \geq SSR(\text{plein})$
et donc $F\geq 0$.
Considérons maintenant une régression avec deux variables
explicatives (donc trois variables la moyenne est aussi estimée) :
```{r}
reg2 <- lm( sr ~ pop15 + pop75, data = lcs)
```
La commande `aov` donne alors
```{r}
aov(reg2)
```
Nous pouvons comparer ces quantités avec les différentes somme des carrés
```{r}
sse(reg1)
sse(reg2)
ssr(reg2)
```
Ainsi, nous voyons que la commande `aov` nous donne une décomposition
des sommes de carrés et des degrés de libertés :
```{r}
cat(sprintf("%2.f + %.2f = %.2f\n",sse(reg1),sse(reg2)-sse(reg1),sse(reg2)))
cat(sprintf("%2.f + %.2f = %.2f\n",sse(reg2),ssr(reg2),sst(reg2)))
```
La somme de la somme des carrés correspond à la somme des carrés totale,
alors que celle des degrés de liberté.
Attention aux degrés de liberté et à la définition de la somme des
carrés. Si nous imposons que $\beta_0=0$, alors les résultats
sont différents car il ne faut pas retirer la moyenne et
les degrés de liberté changent.
```{r}
aov(lm(sr ~ 0 + pop15 + pop75, data=lcs))
```
La commande `summary(aov(...))` donne en plus des informations sur la F-statistique.
Elle est obtenue par la somme des carrées moyennes (c'est‐à-dire divisée
par le nombre de degrés de libertés) sur la somme des carrés moyennes
des résidus :
```{r}
summary(aov(reg2))
```
Ici, le modèle plein est avec 2 variables explicatives, comparons par exemple avec
les sorties de
```{r}
reg2$df.residual*sse(reg1)/ssr(reg2)
reg2$df.residual*(ssr(reg1)-ssr(reg2))/ssr(reg2)
```
__Attention__ la commande `aov` est sensible à l'ordre dans lequel
les arguments de la formule sont entrées, comparons avec
```{r}
summary(aov(lm( sr ~ pop75 + pop15, data=lcs)))
```
La commande `anova` donne des informations similaires en rajoutant
les variables les unes après les autres et permet de s'interroger
sur l'opportunité de rajouter des variables :
```{r}
reg0 <- lm(sr ~ 1, data = lcs)
anova(reg0,reg1,reg2)
```
# Identifions les variables potentiellement significatives
__Q1__
Traçons les variables les unes contre les autres.
```{r}
pairs(lcs)
```
__Q2__
Regardons en détail la régression avec le modèle plein.
```{r}
lm( sr ~ . , data= lcs)
```
```{r}
summary(reg_all<-lm( sr ~ . , data= lcs))
```
Les variables `pop15` et `ddpi` sont identifiées comme significatives.
Sur la base du test de Fisher, on ne peut pas rejeter l'idée
que tous les coefficients $\beta_1,\dotsc,\beta_4$ sont nuls.
```{r}
summary(reg<-lm( sr ~ pop15 + ddpi , data= lcs))
```
Ajoutons la variable `pop75`.
```{r}
summary(reg_inter<-lm( sr ~ pop15 + ddpi + pop75 , data= lcs))
```
Regardons par rappot aux variables non significatives :
```{r}
summary(reg_inter<-lm( sr ~ dpi + pop75 , data= lcs))
```
Nous voyons que la statistique F ne permet par de rejeter
que les coefficients soient nuls.
__Q4.__ Utilisons la commande `drop1`, qui permet de tester
tous les sous-modèles en prenant comme critère l'AIC
(Akaike Information Criterion).
```{r}
drop1(reg_all)
```
Nous en concluons que le modèle sans la variable `dpi`,
c'est-à-dire juste avec `pop15`, `pop75` et `ddpi`
est meilleur en terme d'AIC que celui avec toutes les variables.
Exécutons la commande `anova` entre les deux modèles.
```{r}
anova(reg,reg_all)
```
Cela
correspond à calculer la F-statistique qui est donnée par
```{r}
a<-anova(reg,reg_all)
F<-(45/2)*a$`Sum of Sq`[2]/a$RSS[2]
cat(sprintf("F-statistique = %.4f\np-value = %.2g",F,pf(F,2,45,lower.tail=FALSE)))
```
Nous ne rejetons pas l'hypothèse nulle $H_0$ « $\beta_3=\beta_4=0$ »,
c'est-à-dire que les deux coefficients correspondant à `ddpi` et `pop75`
soient nuls simultanément.
