Commit cc8efc91 authored by Isabelle Sanchez's avatar Isabelle Sanchez
Browse files

version OK sur poste - 00 quasi vide

parent 905453e2
Pipeline #41568 failed with stage
in 16 minutes and 17 seconds
......@@ -43,1360 +43,3 @@ Les données se trouvent dans le package `ade4`. Les variables qui nous intéres
On examine la relation entre une variable à expliquer quantitative $Y$ aléatoire et une variable explicative quantitative $x$, qui est non aléatoire et fixe.
On observe $n$ expériences indépendantes, soit $n$ paires d'observations $(x_i,y_i)_{i=1,\ldots,n}$.
<span style="color:Sienna">**Exemple:** On veut expliquer le nombre de caeca pyloriques `meris$caec` en fonction du diamètre orbital `morpho$DO`. On observe 306 truites. La première chose à faire est de représenter la relation entre des deux variables par un nuage de points.</span>
```{r reglinsimple1, message=FALSE, warning=FALSE,fig.cap="Nombre de caeca pyloriques en fonction du diamètre orbital."}
ggplot(lascaux2, aes(y = caec, x = DO)) + geom_point()
```
On suppose une relation linéaire entre $x$ et $Y$, soit:
$$Y_i=\mu+b x_i, \qquad i=1,\ldots,n.$$
Cependant cette relation n'est jamais exactement vérifiée pour toutes les observations, car aucune ligne droite ne peut passer par tous les points. En effet, il y a des fluctuations aléatoires dues à l'échantillonnage de $Y_i$ (si l'expérience est répétée, une autre valeur $y_i$ de $Y_i$ sera observée). on doit donc prendre en compte un bruit aléatoire et le modèle est défini par:
\begin{equation}
Y_i=\mu+b x_i+\epsilon_i, \qquad i=1,\ldots,n, \qquad \textrm{avec} \quad \epsilon_i\overset{i.i.d.}{\sim} \mathcal{N}(0,\sigma^2). (\#eq:reglinsimple)
\end{equation}
$Y_i$ aléatoire, $\mu$ la constante ou intercept du modèle, $b$ la pente ou paramètre de régression associé à $x$ et $\epsilon_i$ bruit aléatoire gaussien.
Cela revient donc à supposer $Y_i \sim \mathcal{N}(\mu + b x_i, \sigma^2)$, et les $Y_i$ indépendants entre eux. On a trois paramètres inconnus: $\mu$, $b$ et $\sigma^2$.
```{r, reglinsimple2graph,fig.width=8,fig.height=7,echo=FALSE,fig.cap="Illustration du modèle de régression linéaire simple."}
x <- c(3.591013, 2.000000, 3.200000, 3.744146, 4.380719, 0.882786)
y <- c(2.436851, 2.488175, 2.800000, 3.148497, 3.116439, 1.578209)
mod <- lm(y~x)
new <- data.frame(x = x)
pred <- predict(lm(y ~ x), new, se.fit = TRUE)$fit
temp <- data.frame(x,y)
temp2 <- seq(-2.5,0, length = 200)
xnorm2 <- x[2]+c(dnorm(temp2))*2
ynorm2 <- pred[2]-c(temp2)/6
xnorm3 <- x[3]+c(dnorm(temp2))*2
ynorm3 <- pred[3]-c(temp2)/6
temp2 <- seq(0,2.5, length = 200)
xnorm4 <- x[2]+c(dnorm(temp2))*2
ynorm4 <- pred[2]-c(temp2)/6
xnorm5 <- x[3]+c(dnorm(temp2))*2
ynorm5 <- pred[3]-c(temp2)/6
gaussian <- data.frame(temp2,xnorm2,ynorm2,xnorm3,ynorm3,xnorm4,ynorm4,xnorm5,ynorm5)
ggplot(temp, aes(y = y, x =x)) +
geom_point() +
scale_x_continuous(breaks=c(0,x[2],x[3],4.6), labels=c("",expression(x[2]),expression(x[3]),""), limits=c(0,4.6)) +
scale_y_continuous(breaks=c(1.4,y[2],pred[2],y[3],pred[3],3.2), labels=c("",expression(y[2]),expression(a+b*x[2]),expression(y[3]),expression(a+b*x[3]),""),limits=c(1.4,3.2)) +
theme(text = element_text(size=rel(5)), plot.title = element_blank(),axis.title.x = element_blank(),axis.title.y = element_blank()) +
geom_hline(yintercept=1.4) +
geom_vline(xintercept=0) +
stat_smooth(method="lm", se=FALSE, color="blue",fullrange=TRUE) +
geom_segment(aes(x = x[2], y = pred[2], xend = x[2], yend = y[2]),col="seagreen",linetype="dashed") +
geom_segment(aes(x = x[3], y = pred[3], xend = x[3], yend = y[3]),col="seagreen",linetype="dashed") +
geom_text(x=x[2]-0.4,y=2.36,label=expression(y[2]-(a+b*x[2])),color="seagreen",size=6) +
geom_text(x=x[3]-0.4,y=2.76,label=expression(y[3]-(a+b*x[3])),color="seagreen",size=6) +
geom_line(data=gaussian, aes(x=xnorm2,y=ynorm2), col="red") +
geom_line(data=gaussian, aes(x=xnorm3,y=ynorm3), col="red") +
geom_line(data=gaussian, aes(x=xnorm4,y=ynorm4), col="red") +
geom_line(data=gaussian, aes(x=xnorm5,y=ynorm5), col="red") +
geom_segment(aes(x = x[2], y = pred[2], xend = x[2]+dnorm(0)*2, yend = pred[2]),col="red",linetype="dashed") +
geom_segment(aes(x = x[3], y = pred[3], xend = x[3]+dnorm(0)*2, yend = pred[3]),col="red",linetype="dashed") +
geom_text(x=x[2]+0.05,y=1.8,label=expression(paste("distribution de ",Y[2]," en ",x[2])),color="red",size=5, hjust=0) +
geom_text(x=x[2]+0.05,y=1.73,label=expression(paste("mean ", a+b*x[2])),color="red",size=5, hjust=0) +
geom_text(x=x[2]+0.05,y=1.66,label=expression(paste("variance ", sigma^2)),color="red",size=5, hjust=0) +
geom_text(x=x[3]+0.05,y=2.28,label=expression(paste("distribution de ",Y[3]," en ",x[3])),color="red",size=5, hjust=0) +
geom_text(x=x[3]+0.05,y=2.21,label=expression(paste("mean ", a+b*x[3])),color="red",size=5, hjust=0) +
geom_text(x=x[3]+0.05,y=2.14,label=expression(paste("variance ", sigma^2)),color="red",size=5, hjust=0) +
geom_text(x=2.3,y=3.15,label="y = a+b*x",color="blue",size=10,hjust=0)
```
### Estimation {#estimation}
On estime les paramètres $\mu$ et $b$ par la méthode des moindres carrés ou par le maximum de vraisemblance (résultats identiques). De même pour $\sigma^2$ (mais résultats différents entre les moindres carrés et le maximum de vraisemblance). L'hypothèse de normalité des résidus est nécessaire pour la méthode du maximum de vraisemblance, mais ne l'est pas pour les moindres carrés.
\begin{eqnarray}
\nonumber \textrm{MCO:} \quad (\widehat{\mu},\widehat{b}) &=& \textrm{valeurs minimisant} \quad \sum_{i=1}^n \Big(Y_i-(a+bx_i)\Big)^2.\\
\nonumber \textrm{MV:} \quad (\widehat{\mu},\widehat{b}) &=& \textrm{valeurs maximisant la log-vraisemblance}.
\end{eqnarray}
Une fois que les estimateurs $(\widehat{a},\widehat{b})$ sont obtenus, on peut obtenir la droite de régression qui est donnée par:
\begin{displaymath}
y=\widehat{a}+\widehat{b}x.
\end{displaymath}
```{remark}
Les estimateurs $\widehat{a}$ and $\widehat{b}$ sont aléatoires car dépendent des $Y_i$ qui sont eux-mêmes aléatoires. En utilisant nos observations on obtient certaines réalisations de nos estimateurs, mais si l'expérience est répétée, alors de nouvelles paires d'observations seront obtenues et de nouvelles valeurs pour $\widehat{a}$ et $\widehat{b}$ seront obtenues.
```
<span style="color:Sienna">**Exemple:** La méthode utilisée par R est celle du maximum de vraisemblance.</span>
```{r reglinsimple2}
mod1 <- lm(caec ~ DO, data=lascaux2)
summary(mod1)
```
<span style="color:Sienna">On a $\widehat{a}=66.6307$, $\widehat{b}=-2.0982$ et $\widehat{\sigma}=9.963$.
La droite de régression est $y=66.6307-2.0982 DO$.</span>
### Interprétation
Pour une valeur donnée de $x$, la moyenne théorique de $Y$ est $$\mathbb{E}(Y \mid x) = a + b x$$ (le terme d'erreur est supposé centré).
##### L'intercept $\mu$
Il représente la moyenne théorique de $Y$ pour $x=0$. En effet, dans ce cas $\mathbb{E}(Y \mid x) = \mu + b \times 0 = \mu$.
Voir Figure \@ref(fig:interpretation).
Ce paramètre n'est pas toujours biologiquement pertinent. Car parfois $x=0$ ne veut rien dire d'un point de vue biologique.
##### La pente $b$
Imaginons que $x$ augmente de 1. Par exemple $x_2=x_1+1$. La moyenne théorique de $Y_1$ est $\mu + b x_1$, et celle de $Y_2$ est $\mu + b x_2 = \mu +b x_1 + b$. Voir Figure \@ref(fig:interpretation).
Donc, quand $x$ augmente de 1:
* Si $b>0$: $b$ représente l'augmentation dans la moyenne de $Y$.
* Si $b<0$: $b$ représente la diminution dans la moyenne de $Y$.
```{r, interpretation,fig.width=8,fig.height=7,echo=FALSE,fig.cap="Représentation de l'intercept $\\mu$ et de la pente $b$."}
x <- c(3.591013, 2.000000, 3.200000, 3.744146, 4.380719, 0.882786)
y <- c(2.436851, 2.488175, 2.800000, 3.148497, 3.116439, 1.578209)
mod <- lm(y~x)
new <- data.frame(x = x)
pred <- predict(lm(y ~ x), new, se.fit = TRUE)$fit
temp <- data.frame(x,y)
ggplot(temp, aes(y = y, x =x)) +
geom_point() +
scale_x_continuous(breaks=c(0,x[2],x[3],4.6), labels=c("",expression(x[2]),expression(paste(x[3],"=",x[2]+1)),""),limits=c(0,4.6)) +
scale_y_continuous(breaks=c(1.4,y[2],pred[2],y[3],pred[3],3.2), labels=c("",expression(y[2]),expression(mu+b*x[2]),expression(y[3]),expression(mu +b*x[3]),""),limits=c(1.4,3.2)) +
theme(text = element_text(size=rel(5)), plot.title = element_blank(),axis.title.x = element_blank(),axis.title.y = element_blank()) +
geom_hline(yintercept=1.4) +
geom_vline(xintercept=0) +
stat_smooth(method="lm", se=FALSE, color="blue",fullrange=TRUE) +
geom_segment(aes(x = x[3], y = pred[2], xend = x[3], yend = pred[3]),col="seagreen",arrow = arrow(length = unit(0.3, "cm"))) +
geom_segment(aes(x = x[3], y = pred[3], xend = x[3], yend = pred[2]),col="seagreen",arrow = arrow(length = unit(0.3, "cm"))) +
geom_segment(aes(x = x[2], y = pred[2], xend = x[3], yend = pred[2]),col="seagreen",arrow = arrow(length = unit(0.3, "cm"))) +
geom_segment(aes(x = x[3], y = pred[2], xend = x[2], yend = pred[2]),col="seagreen",arrow = arrow(length = unit(0.3, "cm"))) +
geom_text(x=x[3]+0.07,y=y[2]-0.05,label="b",color="seagreen",size=6) +
geom_text(x=2.6,y=pred[2]-0.05,label="1",color="seagreen",size=6) +
geom_text(x=2.3,y=3.15,label=expression(y = mu+b*x),color="blue",size=10,hjust=0)
```
<span style="color:Sienna">**Exemple:** </span>
* <span style="color:Sienna">L'ordonnée à l'origine vaut 66.6307, mais est peu pertinente, car on ne s'intéresse pas au nombre de caeca pyloriques d'une truite de diamètre orbital 0! </span>
* <span style="color:Sienna">La pente vaut -2.0982, donc quand le diamètre orbital augmente de 1mm, le nombre de caeca pyloriques diminue en moyenne de 2.0982. </span>
* <span style="color:Sienna">On peut également calculer par exemple la diminution théorique moyenne du nombre de caeca pyloriques quand le diamètre orbital augmente de 2mm:
\begin{displaymath}
\bigg(\mu + b \times (x+2)\bigg) - \bigg(\mu + b \times x\bigg) = b \times 2,
\end{displaymath}
qui est estimée par $2 \times -2.0982 = -4.1964$, ce qui correspond à une diminution moyenne de 4.1964 caeca pyloriques. </span>
```{r, reglinsimple3,fig.width=7,fig.height=5,fig.cap="Droite de régression estimée sur nos données."}
ggplot(lascaux2, aes(y = caec, x = DO)) +
geom_point() +
stat_smooth(method="lm", se=FALSE, color="red",fullrange=TRUE) +
labs(x ="Diamètre orbital", y = "Nombre de caeca pyloriques")
```
### Valeurs ajustées, résidus et prédiction
* La valeur ajustée de $Y_i$ est $\widehat{Y_i}=\widehat{\mu} + \widehat{b}x_i$.
* Le résidu associé à la $i^{\textrm{ème}}$ mesure est $\widehat{\epsilon_i}=Y_i-\widehat{Y_i}$.
* Soit $x^*$ une nouvelle valeur de $x$ qui n'a pas été observée. La prédiction associée est $\widehat{Y^*}=\widehat{\mu} + \widehat{b}x^*$.
```{r, reglinsimple4,fig.width=8,fig.height=7,echo=FALSE,fig.cap="Droite de régression, valeurs observées, résidus calculés et prédiction."}
x <- c(3.591013, 2.000000, 3.200000, 3.744146, 4.380719, 0.882786)
y <- c(2.436851, 2.488175, 2.800000, 3.148497, 3.116439, 1.578209)
mod <- lm(y~x)
new <- data.frame(x = x)
pred <- predict(lm(y ~ x), new, se.fit = TRUE)$fit
temp <- data.frame(x,y)
ggplot(temp, aes(y = y, x =x)) +
geom_point() +
scale_x_continuous(breaks=c(0,1.5,x[2],x[3],4.6), labels=c("","x*",expression(x[2]),expression(paste(x[3],"=",x[2]+1)),""),limits=c(0,4.6)) +
scale_y_continuous(breaks=c(1.4,mod$coef[1]+mod$coef[2]*1.5,pred[2],y[2],pred[3],y[3],3.2), labels=c("",expression(paste(hat(y),"*=",~hat(mu)+hat(b)*x,"*",sep="")),expression(hat(y[2])==hat(mu)+hat(b)*x[2]),expression(y[2]),expression(hat(y[3])==hat(mu)+hat(b)*x[3]),expression(y[3]),""),limits=c(1.4,3.2)) +
theme(text = element_text(size=rel(5)), plot.title = element_blank(),axis.title.x = element_blank(),axis.title.y = element_blank()) +
geom_hline(yintercept=1.4) +
geom_vline(xintercept=0) +
stat_smooth(method="lm", se=FALSE, color="blue",fullrange=TRUE) +
geom_segment(aes(x = x[2], y = y[2], xend = x[2], yend = pred[2]),col="seagreen",linetype="dashed") +
geom_segment(aes(x = x[3], y = y[3], xend = x[3], yend = pred[3]),col="seagreen",linetype="dashed") +
geom_text(x=x[2]-0.82,y=2.32,label=expression(paste(hat(epsilon)[2]," calculé",sep=" ")),color="seagreen",size=6,hjust=0) +
geom_text(x=x[3]-0.82,y=2.73,label=expression(paste(hat(epsilon)[3]," calculé",sep=" ")),color="seagreen",size=6,hjust=0) +
geom_text(x=2.3,y=3.1,label=expression(LS: y == hat(mu)+hat(b)*x),color="blue",size=8,hjust=0) +
geom_point(x=1.5,y=mod$coef[1]+mod$coef[2]*1.5, color="red", shape=3, size=5) +
geom_segment(aes(x = 0, y = mod$coef[1]+mod$coef[2]*1.5, xend = 1.5, yend = mod$coef[1]+mod$coef[2]*1.5),color="red",linetype="dotted") +
geom_segment(aes(x = 1.5, y = 1.4, xend = 1.5, yend = mod$coef[1]+mod$coef[2]*1.5),color="red",linetype="dotted")
```
<span style="color:Sienna">**Exemple:** </span>
* <span style="color:Sienna">La valeur ajustée pour un diamètre orbital de 10.01 vaut $66.6307 - 10.01 \times 2.0982 = 45.62772$ caeca pyloriques. De même la valeur ajustée pour un diamètre orbital de 9.99 vaut $45.66968$ caeca.</span>
* <span style="color:Sienna">La $47^{\textrm{ème}}$ truite observée a un diamètre orbital de 10.01 mm. Pour cette truite on observe 45 caeca. Le résidu associé vaut donc -0.63.</span>
* <span style="color:Sienna">La $46^{\textrm{ème}}$ truite observée a un diamètre orbital de 9.99 mm. Pour cette truite on observe 37 caeca. Le résidu associé vaut donc -8.67.</span>
* <span style="color:Sienna">La valeur prédite pour un diamètre orbital de 10 (jamais observé) vaut $66.6307 - 10 \times 2.0982 = 45.6487$ caeca pyloriques.</span>
### Inférence statistique
Grâce à l'hypothèse que les $\epsilon_i\overset{i.i.d.}{\sim} \mathcal{N}(0,\sigma^2)$, on obtient des distributions de probabilités pour les paramètres $\widehat{a}$, $\widehat{b}$, $\widehat{\sigma}^2$, les valeurs ajustées $\widehat{Y_i}$ et les prédictions $\widehat{Y^*}$. On peut alors effectuer des tests et obtenir des intervalles de confiance pour ces paramètres, valeurs ajustées et prédiction.
#### Loi de l'estimateur $\widehat{b}$
```{proposition}
Sous le modèle \@ref(eq:reglinsimple), la loi de l'estimateur $\widehat{b}$ est la suivante:
$$\frac{\widehat{b}-b}{\widehat{\sigma_{\widehat{b}}}} \sim \mathcal{S}t_{n-2}.$$
```
On a expliqué précédemment que l'estimateur $\widehat{b}$ est aléatoire. En d'autres termes, si l'expérience est répétée 1000 fois, nous obtiendrons 1000 estimations différentes de $b$, car les observations issues de chaque expérience seront toutes différentes (conséquence de l'aléatoire). La proposition précédente explique que les 1000 estimations différentes de $\frac{\widehat{b}-b}{\hat{\sigma_{\widehat{b}}}}$ suivront une distribution de Student. Si nous faisons un histogramme de ces 1000 estimations, celui-ci sera proche d'une $\mathcal{S}t_{n-2}$, voir la \@ref(fig:lawestimate).
```{r, lawestimate,fig.width=8,fig.height=6,echo=FALSE,fig.cap="Histogramme de 1000 estimations."}
student <- rt(1000,304)
X1 <- as.data.frame(student)
ggplot(X1, aes(x = student)) +
geom_histogram(aes(x = student, y=..density..),fill="lightblue", color="black",bins=40) +
xlim(-4,4) +
geom_vline(aes(xintercept=0), color="red", linetype="dashed", size=1) +
stat_function(fun = dt, args = list(df=304),size=1.5,color="blue") +
ggtitle(expression(paste("Histogram and distribution law for ",(hat(b)-b)/hat(sigma(hat(b))): St[304]))) +
ylab("") +
theme(plot.title=element_text(size = 20), text=element_text(size=rel(4)))
```
#### Test statistique pour $b$
Il est intéressant de tester si $b$ est significativement différent de 0. En effet, $b=0$ indique qu'il n'y a pas de relation linéaire entre $x$ et $Y$, tandis que $b \neq 0$ indique une relation linéaire entre $x$ et $Y$. En pratique on obtiendra toujours $\widehat{b} \neq 0$, mais cela ne veut pas dire que $b \neq 0$. Un test est alors nécessaire:
$$H_0: b=0 \qquad vs \qquad H_1: b \neq 0.$$
En utilisant la proposition précédente, la statistique de test sous $H_0$ est donc:
$$T = \frac{\widehat{b}}{\widehat{\sigma_{\widehat{b}}}} \underset{H_0}{\sim} \mathcal{T}_{n-2}.$$
La région critique se visualise sur la figure \@ref(fig:reglinsimple5).
```{r, reglinsimple5,fig.width=8,fig.height=5,echo=FALSE,fig.cap="Région critique du test de Student pour la pente $b$."}
x <- seq(-2, 12, length = 5000)
plot(c(-2,12), c(-0.06,0.25), type="n",xlab="", ylab="", axes=F)
temp2 <- seq(qnorm(0.975,mean=5,sd=2),12,length=500)
x2 <- c(temp2,qnorm(0.975,mean=5,sd=2))
y2 <- c(dnorm(temp2,mean=5,sd=2),0)
polygon(x2, y2, col="khaki1",border=NA)
temp3 <- seq(-2,qnorm(0.025,mean=5,sd=2),length=500)
x3 <- c(temp3,qnorm(0.025,mean=5,sd=2))
y3 <- c(dnorm(temp3,mean=5,sd=2),0)
polygon(x3, y3, col="khaki1",border=NA)
abline(h=0, lwd=2)
lines(x,dnorm(x,mean=5,sd=2),col="tomato3",lwd=2)
segments(qnorm(0.5,mean=5,sd=2),0,qnorm(0.5,mean=5,sd=2),dnorm(qnorm(0.5,mean=5,sd=2),mean=5,sd=2),col="gray30",lty=2)
segments(qnorm(0.975,mean=5,sd=2),0,qnorm(0.975,mean=5,sd=2),dnorm(qnorm(0.975,mean=5,sd=2),mean=5,sd=2),col="gray30",lwd=2)
segments(qnorm(0.025,mean=5,sd=2),0,qnorm(0.025,mean=5,sd=2),dnorm(qnorm(0.025,mean=5,sd=2),mean=5,sd=2),col="gray30",lwd=2)
text(9.5, 0.03,expression(alpha/2), col = "darkblue",font=1, cex=1.5, pos=4)
text(-1, 0.03,expression(alpha/2), col = "darkblue",font=1, cex=1.5, pos=4)
text(qnorm(0.5,mean=5,sd=2), -0.01, 0, col = "darkblue",font=1)
segments(qnorm(0.025,mean=5,sd=2)+0.05,-0.04,qnorm(0.975,mean=5,sd=2)-0.05,-0.04,col="tomato3",lwd=2)
segments(qnorm(0.975,mean=5,sd=2)+0.05,-0.04,14,-0.04,col="seagreen",lwd=2)
segments(-12,-0.04,qnorm(0.025,mean=5,sd=2)-0.05,-0.04,col="seagreen",lwd=2)
text(5,-0.06, "H0 acceptée", col = "tomato3",font=1, cex=1.4)
text(11,-0.06, "H0 rejetée", col = "seagreen",font=1, cex=1.4)
text(-1,-0.06, "H0 rejetée", col = "seagreen",font=1, cex=1.4)
text(2,0.21, "sous H0", col = "tomato3", font=1, cex=1.4)
text(2,0.23, "distribution de T", col = "tomato3", font=1, cex=1.4)
text(qnorm(0.975,mean=5,sd=2)-1,-0.02,expression(t[paste("n-2,",1-alpha/2)]), col = "darkblue",font=1, cex=1.5, pos=4)
text(qnorm(0.025,mean=5,sd=2)-1.1,-0.02,expression(t[paste("n-2,",alpha/2)]), col = "darkblue",font=1, cex=1.5, pos=4)
```
<span style="color:Sienna">**Exemple:** Nous obtenons $T_{obs}=-2.0982/0.3386=-6.196$.
La région critique est $]-\infty;-1.97] \cup [1.97,\infty[$.
$T_{obs}$ se trouve dans cette région, nous rejetons donc $H_0$: les données ne valident pas l'absence de relation linéaire entre le diamètre orbital et le nombre de caeca pyloriques.</span>
#### Intervalle de confiance pour $b$
```{proposition}
Grâce à la proposition prédente, nous obtenons l'intervalle de confiance de niveau $(1-\alpha)$ pour $b$:
$$ CI_{1-\alpha}(b) \quad = \quad [\widehat{b} - t_{n-2,1-\frac{\alpha}{2}} \times\widehat{\sigma_{\widehat{b}}} \quad,\quad \widehat{b} + t_{n-2,1-\frac{\alpha}{2}} \times \widehat{\sigma_{\widehat{b}}}].$$
```
<span style="color:Sienna">**Exemple:** On obtient les intervalles de confiance de niveau 95\% pour $\mu$ et $b$ avec R de la manière suivante:</span>
```{r reglinsimple6}
confint(mod1)
```
<span style="color:Sienna">L'interprétation pour $b$ est la suivante: si nous répétons l'expérience de nombreuses fois (pêcher 306 truites), dans 95\% des cas le ``vrai'' $b$ sera compris dans l'intervalle de confiance associé à l'expérience. Ici l'intervalle associé à notre expérience est [-2.764582,-1.431855] (l'IC varie entre les différentes expériences possibles, il est aléatoire, tandis que $b$ est fixe, il existe une vraie valeur de $b$ que nous souhaitons estimer).</span>
#### Intervalle de prédiction et hyperboles de prédiction
Pour une nouvelle valeur $x^*$, la prédiction associée $\widehat{Y^*}$ est aléatoire. En effet, si une nouvelle expérience est menée, de nouvelles estimations des paramètres $\widehat{a}$ et $\widehat{b}$ seront obtenues, donc une nouvelle valeur de prédiction $\widehat{Y^*}$ sera obtenue. Pour avoir une idée de comment la prédiction $\widehat{Y^*}$ varie, on utilise un intervalle de prédiction.
```{proposition}
La loi de distribution pour l'erreur de prédiction $Y^* - \widehat{Y^*}$ est:
$$\frac{Y^*-\widehat{Y^*}}{\widehat{\sigma}\sqrt{1+\frac1n + \frac{x^*-\overline{x}}{\sum_{i=1}^n(x_i-\overline{x})^2}}} \sim \mathcal{T}_{n-2}.$$
Un intervalle de prédiction pour $Y^*$ est alors:
\begin{eqnarray}
CI_{1-\alpha}(Y^*) \quad &=& \quad \Bigg[\widehat{Y^*} \quad \pm \quad t_{n-2,1-\frac{\alpha}{2}} \times \widehat{\sigma}\sqrt{1+\frac1n + \frac{x^*-\overline{x}}{\sum_{i=1}^n(x_i-\overline{x})^2}}\Bigg]\\
&=& \quad \Bigg[\widehat{Y^*} \quad \pm \quad t_{n-2,1-\frac{\alpha}{2}} \times \sqrt{\widehat{var(\widehat{\varepsilon^*})}}\Bigg].
\end{eqnarray}
En calculant l'intervalle de prédiction sur tout le support de la droite de régression, nous obtenons deux hyperboles de prédiction.
```
```{remark}
* Nous souhaitons avoir des intervalles de prédiction aussi petits que possible. $\sum(x_i-\overline{x})^2$ doit donc être élevé, soit les $x_i$ doivent être largement répartis autour de $\overline{x}$.
* Plus $x^*$ est loin de $\overline{x}$, plus l'intervalle de prédiction est large. Il est donc dangereux de faire des prédictions pour un $x^*$ éloigné de $\overline{x}$.
```
<span style="color:Sienna">**Exemple:** </span>
* <span style="color:Sienna">La valeur prédite pour un diamètre orbital de 10 vaut $66.6307 - 10 \times 2.0982 = 45.6487$ caeca pyloriques. La valeur prédite pour un diamètre orbital de 20 vaut $66.6307 - 20 \times 2.0982 = 24.6663$ caeca pyloriques.</span>
* <span style="color:Sienna">L'intervalle de prédiction associé à un diamètre orbital de 10 est $[26.00,65.29]$, et celui associé à un diamètre orbital de 20 est $[3.75,45.58]$.</span>
```{r predictions2}
new <- data.frame(DO=c(10,20))
predictions <- predict.lm(mod1,newdata=new,interval="prediction")
predictions
```
* <span style="color:Sienna"> L'intervalle $[3.75,45.58]$ est plus large que $[26.00,65.29]$ car 20 est plus éloigné du DO moyen observé (9.19) que 10. </span>
* <span style="color:Sienna"> Pour visualiser les hyperboles de prédiction pour $Y^*$, voir Figure \@ref(fig:hyperconf).</span>
```{r, hyperconf,fig.width=7,fig.height=5,echo=TRUE,fig.cap="En bleu la droite de régression et en rouge les hyperboles de prédiction pour $Y^*$."}
predictions <- predict.lm(mod1,interval="prediction")
mydata <- cbind(lascaux2,predictions)
ggplot(mydata, aes(y = caec, x = DO)) +
xlim(4,20) +
geom_point() +
stat_smooth(method="lm", se=FALSE, color="blue",fullrange=TRUE) +
geom_line(aes(y = lwr), color = "red", linetype = "dashed")+
geom_line(aes(y = upr), color = "red", linetype = "dashed") +
labs(x ="Diamètre orbital", y = "Nombre de caeca pyloriques") +
ggtitle("Prediction hyperboles") +
theme(plot.title=element_text(size = 20, hjust=0.5), text = element_text(size=rel(4)))
```
### Vérifications des hypothèses sur les résidus observés {#verifhyp}
On a supposé les $\epsilon_i\overset{i.i.d.}{\sim} \mathcal{N}(0,\sigma^2)$, ce qui implique la normalité des résidus, l'homoscédasticité des résidus (ils ont tous même variance), ainsi que leur indépendance. Etant donné que l'inférence statistique est basée sur ces hypothèses, il est nécessaire de les vérifier, au moins par des graphiques.
En général ce ne sont pas les résidus observés qui sont représentés, mais les résidus studentisés, qui sont indépendants entre eux et suivent une $\mathcal{N}(0,1)$ (donc 95\% d'entre eux environ doivent être compris entre -1.96 et 1.96).
**Normalité** Faire un histogramme, un QQ-plot et y tracer la droite de Henry.
**Homoscédasticité** La variance des résidus ne doit pas augmenter en fonction de la valeur de l'observation, ou de la valeur d'une autre variable, comme c'est le cas par exemple sur le graphe (b) de la Figure \@ref(fig:graphhyps).
**Indépendance** Il y a non indépendance si la valeur de $Y$ en $x_i$ dépend des valeurs de $Y$ en d'autres $x_i$, ou bien dépend d'une autre variable. Pour observer cela, il faut tracer les résidus en fonction des variables explicatives. En cas de non indépendance, cela peut être dû:
* Soit à un mauvais choix de modèle: une variable a été oubliée ou la relation est non linéaire, voir les graphes (c) ou (d) de la Figure \@ref(fig:graphhyps).
* Soit à la nature (temporelle ou spatiale) des données, voir le graphe (d) de la Figure \@ref(fig:graphhyps). Par exemple:
+ Si on observe de nombreux oiseaux à un temps $t$, il y a de grandes chances d'en observer encore un grand nombre à $t+1$.
+ Si une espèce marine est beaucoup observée à une certaine profondeur, il est probable que cette espèce soit encore pas mal observée à des profondeurs voisines.
+ Il y a une non-indépendance spatiale et/ou temporelle concernant le nombre d'animaux infectés par une maladie dans une zone forestière.
```{r graphhyps,fig.width=10,fig.height=5,echo=FALSE,fig.cap="Les résidus studentisés sont représentés sur l'axe des ordonnées, et les indices des observations, les valeurs ajustées ou une autre variable sont représentés sur l'axe des abscisses. (a) les hypothèses semblent vérifiés. (b) les résidus sont hétéroscedastiques. (c) tendance des résidus, ils n'ont pas tous une moyenne nulle, une variable d'intérêt a dûe être oubliée dans le modèle. (d) les résidus ne semblent pas indépendants : une variable d'intérêt a dûe être oubliée dans le modèle ou les résidus suivent un modèle autorégressif."}
x <- seq(from=0,to=1000,length.out=1000)
y1 <- rnorm(1000,0,1)
tt <- seq(from=0.5,to=3,length.out=1000)
y2 <- rnorm(1000,rep(0,1000),tt)
y2 <- (y2-mean(y2))/sd(y2)
y3 <- -1 + 0.003 * x + rnorm(1000,mean=0,sd=0.5)
y3 <- y3/var(y3)
y4 <- -(x-500)^2
y4 <- (y4-mean(y4))/sd(y4) + rnorm(1000,mean=0,sd=1)
par(mfrow=c(1,4))
plot(y1 ~ x, xlab="(a)", ylab="studentized residuals",ylim=c(-4,4),cex.lab=1.5)
abline(h=2,lwd=2,col="red")
abline(h=-2,lwd=2,col="red")
plot(y2 ~ x, xlab="(b)", ylab="studentized residuals",ylim=c(-4,4),cex.lab=1.5)
abline(h=2,lwd=2,col="red")
abline(h=-2,lwd=2,col="red")
plot(y3 ~ x, xlab="(c)", ylab="studentized residuals",ylim=c(-4,4),cex.lab=1.5)
abline(h=2,lwd=2,col="red")
abline(h=-2,lwd=2,col="red")
plot(y4 ~ x, xlab="(d)", ylab="studentized residuals",ylim=c(-4,4),cex.lab=1.5)
abline(h=2,lwd=2,col="red")
abline(h=-2,lwd=2,col="red")
```
<span style="color:Sienna">**Exemple:** </span>
```{r graphhyps2,fig.width=6,fig.height=6,echo=TRUE,fig.cap="Vérification des hypothèses sur les résidus de l'exemple. Il faudrait tracer les résidus en fontion des valeurs ajustées et de toutes les variables explicatives possibles. Ici à titre d'exemple on ne les trace qu'en fonction des valeurs ajustées et de la longueur standard `LS`"}
library(tidyverse)
library(cowplot)
res <- mod1$res
fitted <- fitted(mod1)
resfit <- bind_cols(res = res,
fit = fitted,
LS = lascaux2$LS)
plot_grid(nrow = 2, ncol = 2,
ggplot(as_tibble(res), aes(x = res)) +
geom_histogram(bins = 30,
fill = "lightblue", color = "black")+
labs(title = "Histogramme des residus"),
ggplot(as_tibble(res), aes(sample = res)) +
stat_qq() +
stat_qq_line()+
labs(title = "Normal Q-Q Plot"),
ggplot(resfit, aes(y = res, x = fit)) +
geom_point()+
labs(title = "Valeurs ajustees"),
ggplot(resfit, aes(y = res, x = LS)) +
geom_point()+
labs(title = "En fonction de LS"))
```
## Régression linéaire multiple
### Modélisation et estimation
On généralise le modèle précédent \@ref(eq:reglinsimple). Pour chaque observation $i$, nous observons la variable à expliquer quantitative $Y$ ainsi qu'un ensemble de $p$ variables explicatives quantitatives. Le modèle linéaire s'écrit de la façon suivante:
\begin{equation}
Y_i=a + b_1 x_{i1} + b_2 x_{i2} + \ldots + b_p x_{ip} +\epsilon_i, \qquad i=1,\ldots,n, \qquad \textrm{avec} \quad \epsilon_i\overset{i.i.d.}{\sim} \mathcal{N}(0,\sigma^2). (\#eq:reglinmult)
\end{equation}
Avec
$x_i =
\begin{pmatrix}
1\\
x_{i1} \\
x_{i2}\\
\vdots\\
x_{ip}
\end{pmatrix}$ et
$\beta = \begin{pmatrix}
a \\
b_1\\
\vdots\\
b_p
\end{pmatrix}$, on peut écrire:
$$
Y_i=x_i'\beta +\epsilon_i, \qquad i=1,\ldots,n, \qquad \textrm{avec} \quad \epsilon_i\overset{i.i.d.}{\sim} \mathcal{N}(0,\sigma^2).
$$
Pour l'ensemble des observations on obtient:
\begin{eqnarray*}
Y_1 &=& a + b_1 x_{11} + b_2 x_{12} + \ldots + b_p x_{1p} +\epsilon_1\\
Y_2 &=& a + b_1 x_{21} + b_2 x_{22} + \ldots + b_p x_{2p} +\epsilon_2\\
&\vdots& \\
Y_n &=& a + b_1 x_{n1} + b_2 x_{n2} + \ldots + b_p x_{np} +\epsilon_n,
\end{eqnarray*}
ce qui nous amène à l'écriture matricielle suivante:
$$
Y = X \beta +\epsilon \qquad \textrm{avec} \quad \epsilon \sim \mathcal{N}(0,\sigma^2 I),
$$
avec $Y=\begin{pmatrix}
Y_1 \\
Y_2\\
\vdots\\
Y_n
\end{pmatrix}$,
$X = \begin{pmatrix}
1 & x_{11} & \ldots & x_{1p} \\
1 & x_{21} & \ldots & x_{2p} \\
\vdots & \vdots & \ldots & \vdots \\
1 & x_{n1} & \ldots & x_{np}
\end{pmatrix}$ et
$\epsilon=\begin{pmatrix}
\epsilon_1 \\
\epsilon_2\\
\vdots\\
\epsilon_n
\end{pmatrix}$.
```{remark}
Le modèle est linéaire en $\beta$. les variables explicatives peuvent quant à elles être des transformations non linéaires des variables d'origine. Par exemple:
$$
y_i=\beta_1x_{i1}+\beta_2\log(x_{i2})+\beta_3x_{i3}^2+\varepsilon_i, \qquad i=1,\ldots,n, \qquad \textrm{avec les } \varepsilon_i \underset{i.i.d.}{\sim}\mathcal{N}(0,\sigma^2).
$$
```
Pour ce modèle, les estimations sont toujours effectuées par les moindres carrés ordinaires (MCO) ou par maximum de vraisemblance (MV).
Cela donne les estimations suivantes par les moindres carrés ordinaires:
\begin{align*}
\widehat{\beta} &= (X'X)^{-1}X'Y\\
\widehat{\sigma}^2 &= \frac{\sum_{i=1}^n \widehat{\epsilon_i}^2}{n-(p+1)}=\frac{SCR}{n-(p+1)} \quad \textrm{non biaisé}.
\end{align*}
Cela donne les estimations suivantes par le maximum de vraisemblance:
\begin{align*}
\widehat{\beta}_{MV} &= (X'X)^{-1}X'Y\\
\widehat{\sigma}^2_{MV} &= \frac{\sum_{i=1}^n \widehat{\epsilon_i}^2}{n}=\frac{SCR}{n} \quad \textrm{biaisé}.
\end{align*}
Les interprétations des paramètres $\mu$ et $b_k$ sont similaires à l'interprétation de l'intercept et de la pente $b$ pour le modèle simple \@ref(eq:reglinsimple).
<span style="color:Sienna">**Exemple:** On veut expliquer le nombre de caeca pyloriques avec le diamètre orbital et la longueur standard. On a donc deux variables explicatives, soit $\mu$, $b_1$, $b_2$ et $\sigma^2$ à estimer.</span>
```{r reglinmult1}
mod2 <- lm(caec ~ DO + LS, data=lascaux2)
summary(mod2)
```
### Valeurs ajustées, résidus et prédictions
Les valeurs ajustées (pour $i=1,\ldots,n$) sont obtenues par:
\begin{eqnarray*}
\widehat{Y_i} &=& x_i' \widehat{\beta}\\
&=& \widehat{\mu} + \widehat{b_1}x_{1i} + \ldots \widehat{b_p}x_{pi}.
\end{eqnarray*}
Les résidus (pour $i=1,\ldots,n$) sont obtenus par:
\begin{eqnarray*}
\widehat{\epsilon_i} &=& Y_i - \widehat{Y_i} .
\end{eqnarray*}
Les prédictions pour de nouvelles valeurs non observées des variables explicatives $x*=(x_1^*,\ldots,x_p^*)$ sont obtenues par:
\begin{eqnarray*}
\widehat{Y^*} &=& x^{*'} \widehat{\beta}\\
&=& \widehat{\mu} + \widehat{b_1}x_1^* + \ldots \widehat{b_p}x_p^*.
\end{eqnarray*}
<span style="color:Sienna">**Exemple:**</span>
* <span style="color:Sienna">La valeur ajustée du nombre de caeca pyloriques pour un diamètre orbital de 8.20 mm et une longueur de 137.80 mm (truite 14) vaut 64.99129 - 3.91446*8.20 + 0.12660 *137.80 = 50.3382. </span>
* <span style="color:Sienna">La valeur observée pour cette truite 14 vaut 44, donc le résidu associé vaut 44-50.3382=-6.3382. </span>
* <span style="color:Sienna">La valeur prédite pour une truite qui aurait un diamètre orbital de 10 mm et une longueur de 150 mm est de 64.99129 - 3.91446*10 + 0.12660 *150 = 44.83669 caeca pyloriques.</span>
### Inférence statistique
Avec l'hypothèse $\epsilon \sim \mathcal{N}(0,\sigma^2 I)$ l'inférence statistique est possible.
En particulier, nous pouvons, comme dans le cas de la régression linéaire simple mener des tests sur les paramètres, donner des intervalles de confiance pour ces paramètres, et donner des intervalles de prédiction pour les prédictions effectuées.
On s'intéresse ici au test de l'hypothèse $H_0: b_j =0$, ce qui revient à tester l'absence de relation linéaire entre $Y$ et la $j^{\textrm{ème}}$ variable. On utilise pour cela la loi de l'estimateur $\widehat{b}_j$.
```{proposition}
Sous le modèle \@ref(eq:reglinmult), la loi de l'estimateur $\widehat{b}_j$ est la suivante:
\begin{displaymath}
\frac{\widehat{b}_j-b_j}{\widehat{\sigma}\sqrt{(X'X)^{-1}_{jj}}} \sim \mathcal{S}t_{n-p-1}.
\end{displaymath}
```
Les hypothèses du test sont les suivantes:
$$
H_0: b_j=0 \qquad vs \qquad H_1: b_j \neq 0.
$$
En utilisant la proposition ci-dessus, la statistique de test sous $H_0$ est donc:
$$
T_j = \frac{\widehat{b}_j}{\widehat{\sigma}\sqrt{(X'X)^{-1}_{jj}}} \sim \mathcal{S}t_{n-p-1}.
$$
La région critique se visualise sur la figure \@ref(fig:reglinmult1graph).
```{r reglinmult1graph,fig.width=8,fig.height=5,echo=FALSE,fig.cap="Région critique du test de Student pour le paramètre $b_j$."}
x <- seq(-2, 12, length = 5000)
plot(c(-2,12), c(-0.06,0.25), type="n",xlab="", ylab="", axes=F)
temp2 <- seq(qnorm(0.975,mean=5,sd=2),12,length=500)
x2 <- c(temp2,qnorm(0.975,mean=5,sd=2))
y2 <- c(dnorm(temp2,mean=5,sd=2),0)
polygon(x2, y2, col="khaki1",border=NA)
temp3 <- seq(-2,qnorm(0.025,mean=5,sd=2),length=500)
x3 <- c(temp3,qnorm(0.025,mean=5,sd=2))
y3 <- c(dnorm(temp3,mean=5,sd=2),0)
polygon(x3, y3, col="khaki1",border=NA)
abline(h=0, lwd=2)
lines(x,dnorm(x,mean=5,sd=2),col="tomato3",lwd=2)
segments(qnorm(0.5,mean=5,sd=2),0,qnorm(0.5,mean=5,sd=2),dnorm(qnorm(0.5,mean=5,sd=2),mean=5,sd=2),col="gray30",lty=2)
segments(qnorm(0.975,mean=5,sd=2),0,qnorm(0.975,mean=5,sd=2),dnorm(qnorm(0.975,mean=5,sd=2),mean=5,sd=2),col="gray30",lwd=2)
segments(qnorm(0.025,mean=5,sd=2),0,qnorm(0.025,mean=5,sd=2),dnorm(qnorm(0.025,mean=5,sd=2),mean=5,sd=2),col="gray30",lwd=2)
text(9.5, 0.03,expression(alpha/2), col = "darkblue",font=1, cex=1.5, pos=4)
text(-1, 0.03,expression(alpha/2), col = "darkblue",font=1, cex=1.5, pos=4)
text(qnorm(0.5,mean=5,sd=2), -0.01, 0, col = "darkblue",font=1)
segments(qnorm(0.025,mean=5,sd=2)+0.05,-0.04,qnorm(0.975,mean=5,sd=2)-0.05,-0.04,col="tomato3",lwd=2)
segments(qnorm(0.975,mean=5,sd=2)+0.05,-0.04,14,-0.04,col="seagreen",lwd=2)
segments(-12,-0.04,qnorm(0.025,mean=5,sd=2)-0.05,-0.04,col="seagreen",lwd=2)
text(5,-0.06, "H0 acceptée", col = "tomato3",font=1, cex=1.4)
text(11,-0.06, "H0 rejetée", col = "seagreen",font=1, cex=1.4)
text(-1,-0.06, "H0 rejetée", col = "seagreen",font=1, cex=1.4)
text(2,0.21, "sous H0", col = "tomato3", font=1, cex=1.4)
text(2,0.23, "distribution de Tj", col = "tomato3", font=1, cex=1.4)
text(qnorm(0.975,mean=5,sd=2)-1,-0.02,expression(t[paste("n-p-1,",1-alpha/2)]), col = "darkblue",font=1, cex=1.5, pos=4)
text(qnorm(0.025,mean=5,sd=2)-1.1,-0.02,expression(t[paste("n-p-1,",alpha/2)]), col = "darkblue",font=1, cex=1.5, pos=4)
```
<span style="color:Sienna">**Exemple:**
On souhaite tester la présence de relation linéaire entre le nombre de caeca pyloriques et le diamètre orbital. Nous obtenons $T_{DOobs}=-3.91/0.58=-6.71$.\\ La région critique de rejet est $]-\infty;-1.97] \cup [1.97,\infty[$.
Nous rejetons donc $H_0$: les données ne valident pas l'absence de relation linéaire entre le diamètre orbital et le nombre de caeca pyloriques.
Sur la sortie R précédente, les T-test sont effectués pour chaque composant du vecteur $\beta$. On voit que les paramètres associés au diamètre orbital et à la longueur sont tous les deux significativement non nuls (pvalues de 9.74e-11 et 0.000188): les données ne valident pas l'absence de relation linéaire entre le nombre de caeca pyloriques et le diamètre orbital d'une part, et entre le nombre de caeca pyloriques et la longueur d'autre part.
</span>
### Test entre modèles emboités {#test_modeles_emboites}
#### Hypothèses
Supposons un modèle $\mathcal{M}'$ avec $p$ variables explicatives $X_1,X_2,\ldots,X_p$, et un modèle $\mathcal{M}$ avec $q$ variables explicatives $X_1,X_2,\ldots,X_q$, $q<p$. On dit que $\mathcal{M}$ est emboité dans $\mathcal{M}'$.
On veut tester si l'ajout des variables $X_{q+1},\ldots,X_p$ dans $\mathcal{M}$ est pertinent.
Les hypothèses sont les suivantes:
\begin{eqnarray*}
H_0: \quad b_{q+1} = \ldots = b_{p} = 0 \qquad &vs& \qquad H_1: \quad \exists k \in \{q+1,\ldots,p\} \textrm{tel que } b_k \neq 0\\
\textrm{choix modèle $\mathcal{M}$} \qquad &\qquad& \qquad \textrm{choix modèle $\mathcal{M}'$}
\end{eqnarray*}
#### Principe du test et formule de l'anova
On effectue pour cela un test de Fisher. Le principe est de comparer la variabilité expliquée par ces variables supplémentaires à la variabilité résiduelle du modèle avec ces variables $\mathcal{M}'$ (variabilité non expliquée par $\mathcal{M}'$):
* Si les deux sont "proches", les variables supplémentaires expliquent peu la variabilité des données, et on considère qu'il ne faut pas les conserver dans notre modèle (on préfère $\mathcal{M}$).
* Si la variabilité expliquée par ces variables supplémentaires est bien supérieure à la variabilité résiduelle de $\mathcal{M}'$, on considère qu'il faut les conserver (on préfère $\mathcal{M}'$).
On va se baser sur les composants de la formule de l'analyse de la variance. Nous avons pour le modèle $\mathcal{M}$:
$$
y_{i}-\overline{y} = (\widehat{y_i}_{\mathcal{M}}-\overline{y}) + (y_{i}-\widehat{y_i}_{\mathcal{M}}).
$$
En élevant au carré les membres de cette relation et en sommant pour toutes les observations, on obtient:
$$
\sum_{i=1}^n (y_{i}-\overline{y})^2 = \sum_{i=1}^n (\widehat{y_i}_{\mathcal{M}}-\overline{y})^2 + \sum_{i=1}^n (y_{i}-\widehat{y_i}_{\mathcal{M}})^2
$$
Soit
$$ SCT = SCE_{\mathcal{M}} + SCR_{\mathcal{M}}$$
Avec:
* $\frac{SCT}{n-1}$ est un estimateur de la variance totale des données.
* $CME_{\mathcal{M}}=\frac{\textrm{SCE}_{\mathcal{M}}}{q+1}$ est un estimateur de la variance expliquée par $\mathcal{M}$.
* $CMR_{\mathcal{M}}=\frac{\textrm{SCR}_{\mathcal{M}}}{n-q-1}$ est un estimateur de la variance résiduelle du modèle $\mathcal{M}$.
On a le même type de formule pour le modèle $\mathcal{M}'$:
$$
SCT= SCE_{\mathcal{M}'} + SCR_{\mathcal{M}'},
$$
* $CME_{\mathcal{M}'}=\frac{\textrm{SCE}_{\mathcal{M}'}}{p+1}$ est un estimateur de la variance expliquée par $\mathcal{M}'$.
* $CMR_{\mathcal{M}'}=\frac{\textrm{SCR}_{\mathcal{M}'}}{n-p-1}$ est un estimateur de la variance résiduelle du modèle $\mathcal{M}'$.
#### Statistique de test de Fisher
Le principe est de comparer la part de variance expliquée par $X_{q+1},\ldots,X_p$ par rapport à celle non expliquée par $\mathcal{M}'$. On pose donc
\begin{align*}
F &=\frac{SCE_{X_{q+1},\ldots,X_p}/(p-q)}{SCR_{\mathcal{M}'}/(n-p-1)}\\
&=\frac{(SCE_{\mathcal{M}'}-SCE_{\mathcal{M}})/(p-q)}{SCR_{\mathcal{M}'}/(n-p-1)}\\
&=\frac{(SCE_{\mathcal{M}'}-SCE_{\mathcal{M}})/(p-q)}{CMR_{\mathcal{M}'}}\\
&=\frac{(SCR_{\mathcal{M}}-SCR_{\mathcal{M}'})/(p-q)}{CMR_{\mathcal{M}'}}\\.
\end{align*}
Sous $H_0$ $F \sim \mathcal{F}_{p-q,n-p-1}$, ce qui permet d'obtenir la forme de la région critique de rejet de $H_0$ pour un risque de première espèce $\alpha$: $RC=\{X,Y / F > \mathcal{F}_{p-q,n-p-1;1-\alpha}\}$.
<span style="color:Sienna">**Exemple:**
On a `mod1` qui ne contient que la variable explicative diamètre orbital. Le modèle `mod2 `contient en plus la variable longueur standard. Le modèle `mod3` contient en plus la variable nombre de points rouges sur la dorsale, et le modèle `mod4` contient en plus la variable nombre de rayons à la dorsale.
On effectue ci-dessous un test entre modèles emboités `mod2` et `mod4`, avec `anova(mod2,mod4)`. On peut retrouver la somme des carrés expliqués par les variables nombre de points rouges sur la dorsale et nombre de rayons à la dorsale avec `anova(mod4)` seulement si les variables testées ont été entrées en dernier dans la formule du modèle `mod4`.
</span>
```{r reglinmult2}
mod1 <- lm(caec ~ DO, data=lascaux2)
mod2 <- lm(caec ~ DO + LS, data=lascaux2)
mod3 <- lm(caec ~ DO + LS + PRD, data=lascaux2)
mod4 <- lm(caec ~ DO + LS + PRD + rd, data=lascaux2)
anova(mod2,mod4)
anova(mod2)
anova(mod4)
mod4 <- lm(caec ~ PRD + rd + DO + LS, data=lascaux2)
anova(mod4)
anova(mod2,mod4)
```
<span style="color:Sienna"> Dans le cas particulier où une seule variable diffère entre les deux modèles (et donc ici un seul paramètre diffère, ce ne sera pas forcément le cas en anova), on obtient le même résultat avec `anova(mod3,mod4)`, `summary(mod4)` et `drop1(mod4, test="F")` quelque soit l'ordre d'entrée de cette variable dans la formule du modèle `mod4`. On obtient par contre un résultat différent avec `anova(mod4)` si la variable supplémentaire n'a pas été entrée en dernier dans la formule du modèle `mod4`.</span>
```{r reglinmult3}
anova(mod3,mod4)
summary(mod4)
drop1(mod4, . ~ ., test="F")
```
### Vérification des hypothèses