Commit 072ad737 authored by Isabelle Sanchez's avatar Isabelle Sanchez
Browse files

correction de qqs erreurs de syntaxe

parent 633db2a5
Pipeline #41813 failed with stage
in 18 seconds
# Modèle linéaire {#ModLin}
```{r, include=FALSE}
```{r f00-setting, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE, cache=FALSE,eval=TRUE)
library(knitr)
options(width=80)
......@@ -9,8 +9,11 @@ library(ggplot2)
**Exemple de la truite commune dans les cours d'eau du bassin pyrénéen méditerranéen:**
Dans ce chapitre nous utiliserons l'exemple suivant, constitué de données issues de la thèse de doctorat en sciences agronomiques de J.M. Lascaux, soutenue en 1996 et intitulée \textit{Analyse de la variabilité morphologique de la truite commune (Salmo trutta L.) dans les cours d'eau du bassin pyrénéen méditerranéen}. Une description détaillée des données à disposition se trouve dans le document URL http://pbil.univ-lyon1.fr/R/pdf/pps022.pdf.Nous disposons, pour chacune des 306 truites communes de la base de données:* Le **nom de la rivière** dans laquelle elle a été capturée (12 rivières en tout).
Dans ce chapitre nous utiliserons l'exemple suivant, constitué de données issues de la thèse de doctorat en sciences agronomiques de J.M. Lascaux, soutenue en 1996 et intitulée *Analyse de la variabilité morphologique de la truite commune (Salmo trutta L.) dans les cours d'eau du bassin pyrénéen méditerranéen*. Une description détaillée des données à disposition se trouve dans le document URL http://pbil.univ-lyon1.fr/R/pdf/pps022.pdf.
Nous disposons, pour chacune des 306 truites communes de la base de données:
* Le **nom de la rivière** dans laquelle elle a été capturée (12 rivières en tout).
* Son **sexe**.
* De **5 variables méristiques**: nombre de rayons à la dorsale, nombre de rayons de la nageaoire anale, nombre de rayons de la nageoire pelvienne gauche, nombre de rayons de la nageaoire pectorale gauche, nombre de caeca pyloriques.
* De **15 variables de coloration de la robe**. Par exemple, nombre de points rouges avant l'aplomb de la dorsale, nombre de points noirs avant l'aplomb de la dorsale, nombre de points rouges après l'aplomb de l'anale,...
......@@ -22,7 +25,7 @@ En particulier, nous nous intéresserons dans ce cours à la variable *nombre de
Les données se trouvent dans le package `ade4`. Les variables qui nous intéressent sont réparties dans différentes `data.frame`, on les met toutes dans une même `data.frame` appelée `lascaux2`.
```{r lascaux, message=FALSE, warning=FALSE}
```{r f00-lascaux, message=FALSE, warning=FALSE}
library(ade4)
data(lascaux)
library(tidyverse)
......
# Fiche résumée du modèle linéaire simple {#intro}
```{r, include=FALSE}
```{r f01-setting, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE, cache=FALSE,eval=TRUE)
library(knitr)
```
......@@ -29,7 +29,7 @@ Le coefficient de corrélation linéaire est toujours compris entre -1 et 1, ave
- Illustration
On simule trois jeux de données de 15 observations $(x_i,y_i)_{i=1,\ldots,15}$, le premier avec une relation exponentielle, le second avec une relation racine carré et le dernier avec indépendance entre $x$ et $y$.
```{r exocor}
```{r f01-exocor}
x <- sort(round(runif(15,max=1.5),2))
y1 <- exp(x)
y2 <- sqrt(x)
......@@ -37,7 +37,7 @@ y3 <- c(round(runif(15,max=1.5),2))
```
Le résultat est le suivant:
```{r illustration,fig.height=3.5}
```{r f01-illustration,fig.height=3.5}
plot(x,xlim=c(0,1.5),ylim=c(0,exp(1.5)),type="n",xlab="X",ylab="Y",
main="Coefficient de corrélation linéaire pour trois jeux de données")
points(x,y1,col="red", pch=16)
......@@ -52,9 +52,9 @@ legend("topleft",legend=c(paste("y=exp(x), R^2 = ",round(cor(x,y1),2)),
## Régression linéaire simple : des conditions sur le processus d'erreurs
### Exemple
Pour illustrer la notion de régression linéaire, la modélisation de l'ouverture des stomates (variable: gs/max(gs)) en fonction du stress hydrique (variable: FTSW, fraction of the transpirable soil water) servira d'exemple. Les données sont résumés dansla Table \@ref{tab:importation}
Pour illustrer la notion de régression linéaire, la modélisation de l'ouverture des stomates (variable: gs/max(gs)) en fonction du stress hydrique (variable: FTSW, fraction of the transpirable soil water) servira d'exemple. Les données sont résumées dans la Table \@ref{tab:f01-importation}
```{r importation}
```{r f01-importation}
vignes <- read.csv2("data/vignes.csv",stringsAsFactors = TRUE)
kable(summary(vignes),
booktabs=TRUE)
......@@ -77,7 +77,7 @@ $$Y_i = \mu + b_1 x_i + \varepsilon_i$$ où
Conditions qui se résument ainsi : $\varepsilon_i \, iid \, \sim N(0,\sigma)$.
```{r lm, fig.height=3}
```{r f01-lm, fig.height=3}
# graphique des données
plot(gs.rel~FTSW, data=vignes, ylim=c(0,max(vignes$gs.rel)), ylab="Ouverture des stomates")
# régression linéaire
......@@ -100,7 +100,7 @@ de régression.
- L'estimateur sans biais de la variance $\sigma^2$ s'écrit :
$$\hat{\sigma}^2 = \frac{1}{n-2} \sum_{i=1}^n (Y_i - \hat{\mu} - \hat{b}_1x_i)^2$$
```{r estimation}
```{r f01-estimation}
# estimation et test des paramètres de la régression
coef(res.lin)
```
......@@ -109,13 +109,13 @@ coef(res.lin)
### Intervalles de confiance
Par la suite, pour simplifier les notations, l'estimation et l'estimateur d'un
paramètre du modèle s'écriront de la même façon, avec un chapeau. Sous les conditions précédentes et pour calculer les intervalles de confiance.
paramètre du modèle s'écriront de la même façon, avec un chapeau. Sous les conditions précédentes et pour calculer les intervalles de confiance.
- Intercept: On connaît la distribution de $\hat{\mu} \sim N(\mu,\sigma^2/n + \bar{x}^2 (\sigma^2/ S^2_x)/(n-1))$. En remplaçant $\sigma^2$ par son estimation $\hat{\sigma}^2$, on peut calculer un intervalle de confiance $\hat{\mu} \pm t_{\alpha,(n-2)} \, \mid \bar{x}\mid (\hat{\sigma}/ S_x)/\sqrt{n-1}$.
- Pente: On connaît la distribution de $\hat{b}_1 \sim N(b_1, (\sigma^2 / S^2_x)/(n-1))$. En remplaçant $\sigma$ par son estimation $\hat{\sigma}$, on peut calculer un intervalle de confiance $\hat{b_1} \pm t_{\alpha,(n-2)} \, (\hat{\sigma} / S_x)/\sqrt{n-1}$.
```{r estimationIC}
```{r f01-estimationIC}
# Intervalles de confiance
confint(res.lin)
```
......@@ -131,7 +131,7 @@ $$ \hat{Y}_* \pm t_{\alpha,(n-2)} \sqrt{\hat{\sigma}^2 + \hat{\sigma}^2/n + (x_*
une loi normale. On peut calculer un intervalle de confiance
$$\hat{Y}_* \pm t_{\alpha,(n-2)} \sqrt{\hat{\sigma}^2/n + (x_* - \bar{x})^2 (\hat{\sigma}^2/S^2_x)/(n-1)}$$
```{r predictions, fig.height=3}
```{r f01-predictions, fig.height=3}
# graphique des données
plot(gs.rel~FTSW, data=vignes, ylim=c(0,max(vignes$gs.rel)), ylab="Ouverture des stomates")
......@@ -169,7 +169,7 @@ Pour tester si la pente est significativement différente de 0, on utilise le te
+ si p-value < $\alpha$, on rejette $H_0$ et on conclut que la pente est significativement différente de 0 avec un risque $\alpha$ de se tromper.
+ si p-value > $\alpha$, on ne rejette pas $H_0$.
```{r test}
```{r f01-test}
# estimation et test des paramètres de la régression
summary(res.lin)
```
......@@ -207,11 +207,13 @@ $$DC_i = \frac{ \left(\hat{y}_i- (1,x_i) \left( \begin{array}{c} \hat{\mu}_{(-i)
On compare la distance de Cook au quantile 50 \% d'une loi de Fisher $F_{(2,n-2)}$.
```{r grapheresidus}
par(mfrow=c(2,2)); plot(res.lin)
```{r f01-gpresidus, fig.cap="Graphiques diagnostics"}
par(mfrow=c(2,2))
plot(res.lin)
par(mfrow=c(1,1))
```
Objectifs des graphiques des résidus fournis par $R$, Fig. \@ref(fig:grapheresidus):
Objectifs des graphiques des résidus fournis par $R$, Fig. \@ref(fig:f01-gpresidus):
- Le graphique du quadrant supérieur gauche sert à vérifier si les résidus sont centrés en 0 et sans structure.
......
......@@ -78,7 +78,7 @@ ggplot(res.pon, aes(.fitted, sqrt(abs(.stdresid)))) +
ggtitle("Scale-Location")+theme_bw()
```
Le problème d'hétérogénéité de la variance semble gommer par l'utilisation du poids.
Le problème d'hétérogénéité de la variance semble gommé par l'utilisation du poids.
- Comparons les résultats d'estimation des paramètres en régressions linéaires simple et pondérée.
......@@ -109,6 +109,7 @@ On constate que ce poids est trop fort et la variance au lieu d'être stabilisé
- Option poids fixe
Le paramétrage de l'argument _weights_ n'est pas le même avec la function _gls_ du package _nlme_. Dans le modèle de régression linéaire pondérée, la variance $V(Y_i) = V(\varepsilon_i) = \sigma^2 \omega_i^2$ dépend linéairement d'une variable $\omega_i^2$. Cette dépendance est modélisée comme une équation linéaire sans intercept où $\sigma^2$ joue le rôle du coefficient de pente associé à $\omega_i^2$. Dans cette vision, l'argument $weights = ~ \omega_i^2$ définit bien la relation linéaire (avec le symbol $~$, habituel sous _R_). Il s'agit d'un jeu d'écriture et on obtient le même résultat:
```{r}
# Donne les mêmes résultats que la régression pondérée de lm
res.gls <- nlme::gls(gs.rel~FTSW,weights=~FTSW,data=vignes)
......@@ -123,12 +124,14 @@ Une autre façon de paramétriser les poids est d'utiliser des fonctions prédé
$V(\varepsilon_i) = \sigma^2 (variable^{puissance})^2$
On peut choisir de fixer la valeur de la puissance ou estimer ce terme à partir des données. Dans un premier temps, fixons le à 0.5 pour avoir une définition identique aux essais de régression pondérée précédents:
```{r}
res.fix <- gls(gs.rel~FTSW,weights=varPower(form=~FTSW,fixed=0.5),data=vignes)
summary(res.fix)
```
Maintenant, laissons $R$ déterminer la puissance optimum de FTSW dans la définition du poids. Dans ce contexte, l'estimation des paramètres (y compris ceux qui modélisent la variance) est faite par maximum de vraisemblance (_ML_) ou par maximum de vraisemblance restreint (_REML_).
```{r}
res.var <-gls(gs.rel~FTSW,weights=varPower(form=~FTSW),data=vignes, method="ML")
summary(res.var)
......@@ -226,7 +229,7 @@ la décomposition de la variance totale qui n'est pas la même si on travaille e
### Généralisation: comparaison des vraisemblances des modèles, test du maximum de vraisemblance
On emploi le terme de généralisation car ce test peut s'employer dans beaucoup plus de cas, il est générique.
On emploie le terme de généralisation car ce test peut s'appliquer dans beaucoup plus de cas, il est générique.
Définition des termes dans le cas du modèle suivant: $Y_i= f(x_i,\beta) + \varepsilon_i$ où $\varepsilon_i \; \text{indépendants} \; \sim N(0,\sqrt{\sigma^2 \, \omega_i^2})$
......@@ -242,6 +245,7 @@ $$\ln(V) = -\frac{n}{2} \ln(2\pi \sigma^2 \omega_i^2) - \sum_{i=1}^n\frac{(y_i-f
Quand $f(x_i,\beta) = \mu + b_1 x_i$, on voit apparaître le critère des moindres carrés pondérés dans la log vraisemblance car $\sum_{i=1}^n\frac{(y_i-f(x_i,\beta))^2}{2\sigma^2 \omega_i^2} = \frac{1}{\sigma^2}\sum_{i=1}^n (1/\omega_i^2)(Y_i - \mu - b_1 x_i)^2$. Ainsi, maximiser la vraisemblance est équivalent à minimiser les moindres carrés pour estimer le vecteur de paramètres $\beta=(\mu,b_1)$. Pour une même structure de variance, le _LR_ (liklihood ratio test) ou test du rapport des vraisemblances ($\ln V_{H_0}/V_{H_1}$) pour comparer deux modèles linéaires emboîtés est équivalent au test $F$. Il permet de généraliser la comparaison de deux modèles à des cas plus généraux, comme des modèles non linéaires emboîtés, ou des modèles ayant des structure de variance-covariance emboîtées. Ainsi, pour une même fonction $f(x,\beta)$ dans les deux modèles, on peut comparer l'intérêt d'ajouter une structure de variance (par exemple la pondération avec la fonction _varPower_).
Avant il faut s'assurer d'avor estimer les modèles par la même méthode _REML_ (Maximum de vraisemblance restreint), variante de _ML_ qui assure une estimation sans biais de la variance. On constate tout d’abord qu’en comparant des modèles construits avec des méthodes d’estimation différentes, le test L.Ratio n’est pas effectué :
```{r}
anova(res.lin,res.cov) # Test L.Ratio non effectué!
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment