Commit ddd47db1 authored by Benedicte FONTEZ's avatar Benedicte FONTEZ
Browse files

modification mise en forme convenu avec Léo

parent f7cdd2e5
Pipeline #41424 failed with stage
in 16 minutes and 5 seconds
This diff is collapsed.
......@@ -64,7 +64,7 @@ kable(summary(vignes),
### Définition
- Le modèle de régression linéaire simple de $Y/X=x$ s'écrit :
$$Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$$ où
$$Y_i = \mu + b_1 x_i + \varepsilon_i$$ où
+ $Y$ est la variable à expliquer, elle est dépendante de $X$ et aléatoire.
+ $X$ est la variable explicative, elle est indépendante et supposée fixe.
+ $(Y_1,\ldots,Y_n)$ désigne le n-uplet de variables aléatoires d'un échantillon de taille $n$ associé au design $(x_1,\ldots,x_n)$.
......@@ -88,17 +88,17 @@ res.lin <- lm(gs.rel~FTSW,data=vignes)
- Les estimateurs obtenus en minimisant le critère des moindres carrés
ordinaires sont les suivants :
- Pour l'intercept $\beta_0$ : $\hat{\beta_0} = \bar{Y}- \hat{\beta_1} \bar{x}$.
- Pour la pente $\beta_1$ : $\hat{\beta_1}= S_{xy}/S^2_x$.
- Pour l'intercept $\mu$ : $\hat{\mu} = \bar{Y}- \hat{b_1} \bar{x}$.
- Pour la pente $b_1$ : $\hat{b_1}= S_{xy}/S^2_x$.
Le critère des moindres carrés ordinaires s'écrit:
$$\min_{(\beta_0,\beta_1)} \sum_{i=1}^n (Y_i - \beta_0 - \beta_1 x_i)^2$$
$$\min_{(\mu,b_1)} \sum_{i=1}^n (Y_i - \mu - b_1 x_i)^2$$
- A partir de ces estimateurs, on définit une prédiction pour $Y_i$ : $\hat{Y}_i= \hat{\beta_0} + \hat{\beta_1} x_i$, qui correspond à la valeur donnée par la droite
- A partir de ces estimateurs, on définit une prédiction pour $Y_i$ : $\hat{Y}_i= \hat{\mu} + \hat{b}_1 x_i$, qui correspond à la valeur donnée par la droite
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{\beta_0} - \hat{\beta_1}x_i)^2$$
$$\hat{\sigma}^2 = \frac{1}{n-2} \sum_{i=1}^n (Y_i - \hat{\mu} - \hat{b}_1x_i)^2$$
```{r estimation}
# estimation et test des paramètres de la régression
......@@ -111,9 +111,9 @@ coef(res.lin)
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.
- Intercept: On connaît la distribution de $\hat{\beta_0} \sim N(\beta_0,\sigma^2/n + \bar{x}^2 (\sigma^2/ S^2_x)/(n-1))$. En remplaçant $\sigma$ par son estimation $\hat{\sigma}$, on peut calculer un intervalle de confiance $\hat{\beta_0} \pm t_{\alpha,(n-2)} \, \mid \bar{x}\mid (\hat{\sigma}/ S_x)/\sqrt{n-1}$.
- 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{\beta_1} \sim N(\beta_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{\beta_1} \pm t_{\alpha,(n-2)} \, (\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}
# Intervalles de confiance
......@@ -155,11 +155,11 @@ lines(FTSW,
Pour tester si la pente est significativement différente de 0, on utilise le test de Student décrit ci-après:
- Hypothèses de test :
+ $H_0 : \,\beta_1 = 0$, les variables $X$ et $Y$ ne sont pas liées selon le modèle linéaire.
+ $H_1 : \, \beta_1 \neq 0$, Les variables $X$ et $Y$ sont liées linéairement et dépendent l'une de l'autre.
+ $H_0 : \,b_1 = 0$, les variables $X$ et $Y$ ne sont pas liées selon le modèle linéaire.
+ $H_1 : \, b_1 \neq 0$, Les variables $X$ et $Y$ sont liées linéairement et dépendent l'une de l'autre.
- Statistique de test et sa loi sous l'hypothèse nulle ($H_0$):
$T = \frac{\hat{\beta_1} - 0}{\sqrt{(\hat{\sigma}^2 / S^2_x)/(n-1)}}\sim T_{n-2}$.
$T = \frac{\hat{b_1} - 0}{\sqrt{(\hat{\sigma}^2 / S^2_x)/(n-1)}}\sim T_{n-2}$.
- Risque de première espèce $\alpha = 5\%$
......@@ -201,9 +201,9 @@ Pour pouvoir interpréter les graphiques, certains termes doivent être définis
$$h_{ij} = (1,x_i)(\mathbf{X}'\mathbf{X})^{-1} \left( \begin{array}{c} 1 \\ x_j \end{array} \right)$$
- Distance de Cook: Cette distance mesure l'influence de l'observation $i$ sur l'estimation des paramètres $(\beta_0,\beta_1)$.
- Distance de Cook: Cette distance mesure l'influence de l'observation $i$ sur l'estimation des paramètres $(\mu,b_1)$.
$$DC_i = \frac{ \left(\hat{y}_i- (1,x_i) \left( \begin{array}{c} \hat{\beta_0}_{(-i)} \\ \hat{\beta_1}_{(-i)} \end{array} \right) \right)^2}{2 \hat{\sigma}^2}$$
$$DC_i = \frac{ \left(\hat{y}_i- (1,x_i) \left( \begin{array}{c} \hat{\mu}_{(-i)} \\ \hat{b_1}_{(-i)} \end{array} \right) \right)^2}{2 \hat{\sigma}^2}$$
On compare la distance de Cook au quantile 50 \% d'une loi de Fisher $F_{(2,n-2)}$.
......
......@@ -19,26 +19,26 @@ seed <- 1859; set.seed(seed)
### Le modèle de régression linéaire pondérée de Y/X=x
$$Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \; V(\varepsilon_i) = \sigma^2 \, \omega_i^2$$
$$Y_i = \mu + b_1 x_i + \varepsilon_i, \; V(\varepsilon_i) = \sigma^2 \, \omega_i^2$$
La condition d'homogénéité des variances n'est plus respectée et $\varepsilon_i$ iid $\sim N(0,\sigma^2 \, \omega_i^2)$. Dans ce contexte, on peut se ramener au cas de la régression linéaire ordinaire en appliquant la transformation suivante :
$$\frac{Y_i}{\omega_i} = \beta_0\frac{1}{\omega_i} + \beta_1 \frac{x_i}{\omega_i} + \frac{\varepsilon_i}{\omega_i}$$
$$\frac{Y_i}{\omega_i} = \mu\frac{1}{\omega_i} + b_1 \frac{x_i}{\omega_i} + \frac{\varepsilon_i}{\omega_i}$$
alors
$$Y_i^p = \beta_0 x_{i0}^p + \beta_1 x_{i1}^p + \varepsilon_i^p$$
$$Y_i^p = \mu x_{i0}^p + b_1 x_{i1}^p + \varepsilon_i^p$$
### Estimation
Les estimateurs sont obtenus en minimisant le critère des moindres carrés pondérés:
$$\min_{(\beta_0,\beta_1)} \sum_{i=1}^n (Y_i^p - \beta_0 x_{i0}^p - \beta_1 x_{i}^p)^2$$
$$\min_{(\mu,b_1)} \sum_{i=1}^n (Y_i^p - \mu x_{i0}^p - b_1 x_{i}^p)^2$$
soit
$$\min_{(\beta_0,\beta_1)} \sum_{i=1}^n (1/\omega_i^2)(Y_i - \beta_0 - \beta_1 x_i)^2$$
$$\min_{(\mu,b_1)} \sum_{i=1}^n (1/\omega_i^2)(Y_i - \mu - b_1 x_i)^2$$
On note $p_i = (1/\omega_i^2)$ le poids (weight) et il suffit de préciser le poids dans le programme de minimisation pour avoir les bons estimateurs et estimations pour $\beta_0$ et $\beta_1$. On utilisera l'estimateur sans biais habituel de la variance sur les données transformées $\hat{\sigma}^2 = \frac{1}{n-2}\sum_{i=1}^n(Y_i^P - \hat{Y}_i^P)^2$.
On note $p_i = (1/\omega_i^2)$ le poids (weight) et il suffit de préciser le poids dans le programme de minimisation pour avoir les bons estimateurs et estimations pour $\mu$ et $b_1$. On utilisera l'estimateur sans biais habituel de la variance sur les données transformées $\hat{\sigma}^2 = \frac{1}{n-2}\sum_{i=1}^n(Y_i^P - \hat{Y}_i^P)^2$.
_Quand les données ne sont disponibles que pour des valeurs moyennes (par exemple des valeurs moyennes par variété). On doit utiliser une régression linéaire pondérée où le poids $p_i = n_i$ dépend du nombre d'individus utilisés pour calculer les valeurs moyennes ($\bar{x}_i$ et $\bar{y}_i$) de l'observation $i$._
......@@ -144,16 +144,17 @@ http://stat.ethz.ch/R-manual/R-devel/library/nlme/html/corClasses.html
__definition__
L'analyse de la covariance permet de modéliser le lien entre une variable quantitative $Y$ et des variables explicatives de type quantitatif et de type qualitatif. On couple ensemble les notions d'anova et de regression linéaire. On se restreindra au cas de deux variables explicatives, l'une quantitative et l'autre qualitative. Le modèle s'écrit dans ce cas:
$$Y_{ij} = \beta_{0,i} + \beta_{1,i} x_{ij} +\varepsilon_{ij}, \; \varepsilon_{ij} \, \sim \, N(0,\sigma^2 \, \omega_{ij}^2)$$
$$Y_{ij} = (\mu +\alpha_{i}) + (b_{1} + \gamma_i) x_{ij} +\varepsilon_{ij}, \; \varepsilon_{ij} \, \sim \, N(0,\sigma^2 \, \omega_{ij}^2) \\
= \mu_{i} + b_{1,i} x_{ij} +\varepsilon_{ij}, \; \varepsilon_{ij} \, \sim \, N(0,\sigma^2 \, \omega_{ij}^2)$$
1.Indices: $i$ désigne la modalité de la variable qualitative et $j$
l'individu.
2. Intercept: sa valeur dépend de la modalité $i$. On peut poser $\beta_{0,i} = \mu + \alpha_i$, ce qui correspond à une autre écriture du même modèle pour laquelle on doit définir une contrainte d'identifiabilité (par exemple $\alpha_1 = 0$ ou $\sum \alpha_i=0$).
2. Intercept: sa valeur dépend de la modalité $i$. On peut poser $\mu_{i} = \mu + \alpha_i$, ce qui correspond à une autre écriture du même modèle pour laquelle on doit définir une contrainte d'identifiabilité (par exemple $\alpha_1 = 0$ ou $\sum \alpha_i=0$).
3. Pente: sa valeur dépend de la modalité $i$. On peut poser $\beta_{1,i} = \beta + \gamma_i$, ce qui correspond à une autre écriture du même modèle pour laquelle on doit définir une contrainte d'identifiabilité (par exemple $\gamma_1 = 0$ ou $\sum \gamma_i=0$).
3. Pente: sa valeur dépend de la modalité $i$. On peut poser $b_{1,i} = b_1 + \gamma_i$, ce qui correspond à une autre écriture du même modèle pour laquelle on doit définir une contrainte d'identifiabilité (par exemple $\gamma_1 = 0$ ou $\sum \gamma_i=0$).
L'hétérogénéité observé pour la régression linéaire précédente de l'ouverture des stomates (_gs.rel_) en fonction du stress hydrique (_FTSW_) peut s'expliquer par des relations linéaires différentes selon le couple
variété-année (_Label_).
......@@ -187,17 +188,17 @@ On a ajusté un modèle complexe, on va comparer ce modèle à des modèles plus
### Comparaison de modèles linéaires emboîtés : test F
__Définition modèles emboîtés__
le modèle $f(x,\theta)$ est emboîté dans le modèle
$g(x,(\theta,\delta))$, s'il existe une valeur de $\delta=\delta^*$ telle que
$f(x,\theta)=g(x,(\theta,\delta=\delta^*))$.
le modèle $f(x,\beta)$ est emboîté dans le modèle
$g(x,(\beta,\delta))$, s'il existe une valeur de $\delta=\delta^*$ telle que
$f(x,\beta)=g(x,(\beta,\delta=\delta^*))$.
Par exemple, le modèle linéaire $\beta_0 + \beta_1 x$ est emboîté dans le modèle de covariance $\beta_{0,i} + (\beta_{1,j}) x$ car il suffit de poser que $(\beta_{0,i},\beta_{1,i})=(\beta_0,\beta_1)$ pour tout $i$ pour retrouver le modèle linéaire.
Par exemple, le modèle linéaire $\mu + b_1 x$ est emboîté dans le modèle de covariance $\mu_{i} + (b_{1,j}) x$ car il suffit de poser que $(\mu_{i},b_{1,i})=(\mu,b_1)$ pour tout $i$ pour retrouver le modèle linéaire.
__Définition modèle linéaire__
Un modèle est dit linéaire s'il est linéaire en ses paramètres.
Par exemple, $f(x,(\beta_0,\beta_1))=\beta_0 + \beta_1 x^2$ est un modèle linéaire alors que
$f(x,(\beta_0,\beta_1))= \beta_0 + x ^{\beta_1}$ est un modèle non linéaire.
Par exemple, $f(x,(\mu,b_1))=\mu + b_1 x^2$ est un modèle linéaire alors que
$f(x,(\mu,b_1))= \mu + x ^{b_1}$ est un modèle non linéaire.
On peut constuire un test pour comparer des modèles linéaires emboîtés :
* Hypothèses de test
......@@ -225,18 +226,18 @@ la décomposition de la variance totale qui n'est pas la même si on travaille e
On emploi le terme de généralisation car ce test peut s'employer 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,\theta) + \varepsilon_i$ où $\varepsilon_i \; \text{indépendants} \; \sim N(0,\sqrt{\sigma^2 \, \omega_i^2})$
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})$
* Vraisemblance: Probabilité ou densité de l'échantillon vue comme une fonction des paramètres.Comme, $Y_i \sim N(f(x_i,\theta),\sigma \, \omega_i)$, et que les $Y_i$ sont supposés indépendants, on peut écrire la vraisemblance sous forme d'un produit:
* Vraisemblance: Probabilité ou densité de l'échantillon vue comme une fonction des paramètres.Comme, $Y_i \sim N(f(x_i,\beta),\sigma \, \omega_i)$, et que les $Y_i$ sont supposés indépendants, on peut écrire la vraisemblance sous forme d'un produit:
$$f(\theta,\sigma) = \frac{1}{\sqrt{2\pi \sigma^2 \omega_i^2}}e^{-\frac{(y_1-f(x_1,\theta))^2}{2\sigma^2 \omega_i^2}}\times \ldots \times
\frac{1}{\sqrt{2\pi \sigma^2 \omega_i^2}}e^{-\frac{(y_n-f(x_n,\theta))^2}{2\sigma^2 \omega_i^2}}$$
$$f(\beta,\sigma) = \frac{1}{\sqrt{2\pi \sigma^2 \omega_i^2}}e^{-\frac{(y_1-f(x_1,\beta))^2}{2\sigma^2 \omega_i^2}}\times \ldots \times
\frac{1}{\sqrt{2\pi \sigma^2 \omega_i^2}}e^{-\frac{(y_n-f(x_n,\beta))^2}{2\sigma^2 \omega_i^2}}$$
* Estimation par maximum de vraisemblance: On cherche la valeur des paramètres $(\theta,\sigma^2)$ qui maximise la vraisemblance. Maximiser la vraisemblance est équivalent à maximiser la log vraisemblance :
* Estimation par maximum de vraisemblance: On cherche la valeur des paramètres $(\beta,\sigma^2)$ qui maximise la vraisemblance. Maximiser la vraisemblance est équivalent à maximiser la log vraisemblance :
$$\ln(V) = -\frac{n}{2} \ln(2\pi \sigma^2 \omega_i^2) - \sum_{i=1}^n\frac{(y_i-f(x_i,\theta))^2}{2\sigma^2 \omega_i^2}$$
$$\ln(V) = -\frac{n}{2} \ln(2\pi \sigma^2 \omega_i^2) - \sum_{i=1}^n\frac{(y_i-f(x_i,\beta))^2}{2\sigma^2 \omega_i^2}$$
Quand $f(x_i,\theta) = \beta_0 + \beta_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,\theta))^2}{2\sigma^2 \omega_i^2} = \frac{1}{\sigma^2}\sum_{i=1}^n (1/\omega_i^2)(Y_i - \beta_0 - \beta_1 x_i)^2$. Ainsi, maximiser la vraisemblance est équivalent à minimiser les moindres carrés pour estimer le vecteur de paramètres $\theta=(\beta_0,\beta_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,\theta)$ 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_).
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}
......@@ -254,6 +255,11 @@ anova(res.lin,res.pon) # Test L.Ratio effectué car bonne prise en compte des df
anova(res.cov,res.cov.pon)
```
Le Schéma des comparaisons de modèles possibles est le suivant
![alt text](_bookdown_files/Cours_files/figure-html/SchemaTest.png)
_Pour mémoire_: basés sur la vraisemblance, des critères pour évaluer la qualité de l'ajustement
du modèle ont été proposés :
......@@ -263,6 +269,18 @@ du modèle ont été proposés :
Plus leur valeur est petite meilleur est l'ajustement.
Enfin, on peut maintenant vérifier les résidus pour valider le modèle retenu.
```{r}
# graphiques des résidus
ggplot(vignes, aes(fitted(res.cov.pon), residuals(res.cov.pon, type="normalized"))) +
geom_point(na.rm=TRUE) +
geom_hline(yintercept=0,col="red") +
xlab("Fitted Value") +
ylab("Standardized residuals") +
ggtitle("Validation de l'ancova pondérée")+theme_bw()
```
### Critères de comparaison de modèles pour la prédiction
* PRESS: Predictive residuals sum of squares. Cette valeur prédictive se calcule avec la formule suivante:
......
This diff is collapsed.
# Régression non linéaire
__Définition__
On pose le modèle suivant $Y_i = f(x_i,\theta) + \varepsilon_i$ où $\varepsilon_i \; iid \; N(0,\sigma)$ et où la fonction $f$ est définie par :
On pose le modèle suivant $Y_i = f(x_i,\beta) + \varepsilon_i$ où $\varepsilon_i \; iid \; N(0,\sigma)$ et où la fonction $f$ est définie par :
- Fonction non linéaire: _f_ est une fonction non linéaire du vecteur de paramètres $\theta$. exemple : $f(x,\theta) = \frac{\beta_1 x}{\beta_2 + x}$
- Fonction non linéaire: _f_ est une fonction non linéaire du vecteur de paramètres $\beta$. exemple : $f(x,\beta) = \frac{b_1 x}{b_2 + x}$
- Fonction intrinsèquement non linéaire: $f$ est une fonction du vecteur de paramètres qui ne peut être linéarisée (Seber et Wild). Par exemple, le modèle précédent peut être linéarisé car $1/f(x,\theta) = 1/\beta_1 + \beta_3 \times \frac{1}{x}$.
- Fonction intrinsèquement non linéaire: $f$ est une fonction du vecteur de paramètres qui ne peut être linéarisée (Seber et Wild). Par exemple, le modèle précédent peut être linéarisé car $1/f(x,\beta) = 1/b_1 + b_3 \times \frac{1}{x}$.
Une autre définition existe (Bates et Watts), indiquant que $f$ est intrinsèquement non linéaire si $f(x_i,\theta)$ n'est pas linéaire dans le voisinage de $\hat{\theta}$. Cette définition est utilisée pour définir les propriétés de la méthode d'estimation.
Une autre définition existe (Bates et Watts), indiquant que $f$ est intrinsèquement non linéaire si $f(x_i,\beta)$ n'est pas linéaire dans le voisinage de $\hat{\beta}$. Cette définition est utilisée pour définir les propriétés de la méthode d'estimation.
## Linéarisation d'un modèle
Les motivations pour linéariser un modèle sont les suivantes:
Une première approche peut consister en la transformation des données afin de **linéariser le modèle**. Les motivations pour linéariser un modèle sont les suivantes:
- Linéariser le modèle et utiliser les outils de régression linéaire.
- Normaliser les erreurs $\varepsilon_i$.
......@@ -20,25 +20,25 @@ Quand on applique une transformation, il faut toujours vérifier l'impact de cet
transformation $\ln$ :
Si on observe une relation exponentielle entre la variable à expliquer $Y$ et la variable explicative $X$ et qu'en plus la variance des $Y_i$ est hétérogène, on peut supposer que l'erreur dans le modèle est multiplicative. Dans ce contexte la transformation $\ln$ permet de stabiliser la variance des $Y_i$ tout en linéarisant le modèle:
$$Y_i = exp(\beta_0 + \beta_1 x_i) \varepsilon_i$$
$$Y_i = exp( + b_1 x_i) \varepsilon_i$$
soit en échelle log
$$\ln(Y_i) = \beta_0 + \beta_1 x_i + \ln(\varepsilon_i)$$
$$\ln(Y_i) = \mu + b_1 x_i + \ln(\varepsilon_i)$$
Si la variance n'est pas hétérogèneet que l'erreur est de type additive, alors la transformation $\ln$ rendra la variance de $Y_i$ dépendante de $x_i$, car:
$$Y_i = exp(\beta_0 + \beta_1 x_i) + \varepsilon_i$$
$$Y_i = exp(\mu + b_1 x_i) + \varepsilon_i$$
devient avec la transformation logarithme:
$$\ln(Y_i) = \beta_0 + \beta_1 x_i + \ln(1+\varepsilon_i /exp(\beta_0 + \beta_1 x_i))$$
$$\ln(Y_i) = \mu + b_1 x_i + \ln(1+\varepsilon_i /exp(\mu + b_1 x_i))$$
et le terme d'erreur est lié à la variable explicative $X$.
On peut observer l'impact de la transformation sur la variance avec le graphique _scale-location_ des résidus sous $R$ et l'impact de la transformation sur la normalité avec le graphique _qqnorm_ des résidus sous $R$.
On peut observer **l'impact de la transformation sur la variance avec le graphique _scale-location_ des résidus** sous $R$ et **l'impact de la transformation sur la normalité avec le graphique _qqnorm_ des résidus** sous $R$.
Nous allons maintenant étudier comment appliquer la transformation inverse quand après estimation, on souhaite revenir aux données initiales.
On considère une réaction enzymatique (Bates et Watts) impliquant la Puromycin. On observe la vitesse de réaction ($nb.min^{-1}$) en fonction de la quantité de substrat ($ppm$). Le modèle habituellement utilisé est celui de
Michaelis Menten :
$$f(x,\theta) = \frac{\beta_1 x}{\beta_2 + x}$$
$$f(x,\beta) = \frac{b_1 x}{b_2 + x}$$
Linéarisable en
$$1/ f(x,\theta) = \beta_0^l + \beta_1^l \frac{1}{x}$$
$$1/ f(x,\beta) = \mu^l + b_1^l \frac{1}{x}$$
1. Représentons graphiquement les données transformées.
```{r, fig.height=3}
......@@ -46,8 +46,16 @@ $$1/ f(x,\theta) = \beta_0^l + \beta_1^l \frac{1}{x}$$
X <- c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.1,1.1)
Y <- c(76,47,97,107,123,139,159,152,191,201,207,200)
data <- bind_cols(
tibble(X),
tibble(Y)
)
#graphique
plot(1/Y~I(1/X),xlab="substrat",ylab="vitesse de réaction",main="Données transformées")
ggplot(data, aes(y = 1/Y, x = 1/X)) +
geom_point(size = 2) +
labs(x = "Substrat", y = "Vitesse de réaction", title = "Données transformées")+
theme_bw()
```
......@@ -57,8 +65,12 @@ plot(1/Y~I(1/X),xlab="substrat",ylab="vitesse de réaction",main="Données trans
res.l <- lm(1/Y~I(1/X)) # modèle linéarisé
#graphique
plot(1/Y~I(1/X),xlab="substrat",ylab="vitesse de réaction",main="Données transformées")
abline(coef(res.l),col="red",lwd=2)
ggplot(data, aes(y = 1/Y, x = 1/X)) +
geom_point(size = 2) +
geom_abline(intercept = coef(res.l)[1], slope = coef(res.l)[2],
color = "darkred", size = 1) +
labs(x = "Substrat", y = "Vitesse de réaction", title = "Données transformées")+
theme_bw()
```
......@@ -70,51 +82,63 @@ summary(res.l)
4. Représentons graphiquement les données non transformées. Rajoutons la courbe obtenue avec la transformation inverse. On constate que la courbe détransformée est au niveau du plateau très en dessous des valeurs observées (mauvais ajustement).
```{r, fig.height=3}
# transformation inverse sans correction
beta1 <- 1/coef(res.l)[1]
beta2 <- coef(res.l)[2]/coef(res.l)[1]
b1 <- 1/coef(res.l)[1]
b2 <- coef(res.l)[2]/coef(res.l)[1]
data2 <- data %>% mutate(courbe = (X*b1)/(X+b2))
plot(Y~X,xlab="substrat",ylab="vitesse de réaction",main="Puromyin experiment")
lines(X,(X*beta1)/(X+beta2),col="red",lwd=2)
ggplot(data2, aes(y = Y, x = X)) +
geom_point(size = 2) +
geom_line(aes(y = courbe),
color = "red", size = 1) +
labs(x = "Substrat", y = "Vitesse de réaction", title = "Données transformées")+
theme_bw()
```
En effet la fonction $f(x_i,\theta)$ est la valeur observée en moyenne, c'est l'espérance de $E(Y_i \mid x_i)$. Or, on ne peut inverser l'ordre entre espérance et transformation non linéaire: $E(1/Y_i) \neq 1/(E(Y_i))$. Donc, pour calculer $E(Y_i)$ à partir de $E(1/Y_i)$, il faut appliquer une autre méthode. Si, comme dans notre cas, la transformation utilisée est monotone alors, on peut utiliser
En effet la fonction $f(x_i,\beta)$ est la valeur observée en moyenne, c'est l'espérance de $E(Y_i \mid x_i)$. Or, on ne peut inverser l'ordre entre espérance et transformation non linéaire: $E(1/Y_i) \neq 1/(E(Y_i))$. Donc, pour calculer $E(Y_i)$ à partir de $E(1/Y_i)$, il faut appliquer une autre méthode. Si, comme dans notre cas, la transformation utilisée est monotone alors, on peut utiliser
la définition des quantiles car (par exemple) :$P(Y_i <q_{0.975}) = 0.975$ et $P(1/Y_i<1/q_{0.975}) = 0.975$. Donc, si pour le modèle linéarisé (en $1/Y$), on calcule l'intervalle de confiance $[lower_i,upper_i]$ de la prédiction moyenne pour tous les $x_i$ observés, on sait d'après la constatation précédente que $[1/lower_i,1/upper_i]$ est l'intervalle de confiance pour $\hat{Y}_i$ et on pose $\hat{y}_i = (1/lower_i + 1/upper_i)/2$.
```{r}
# correction
est_cor <- (1/predict(res.l,interval="confidence")[,"lwr"] +
1/predict(res.l,interval="confidence")[,"upr"])/2
data3 <- data2 %>%
mutate(est_cor = (1/predict(res.l,interval="confidence")[,"lwr"] +
1/predict(res.l,interval="confidence")[,"upr"])/2) %>%
pivot_longer(-c(X,Y), names_to = "variable", values_to = "value") %>%
mutate(legende = ifelse(variable == "courbe", "1/Moyenne", "Milieu de 1/IC"))
# Graphique
plot(Y~X,xlab="substrat",ylab="vitesse de réaction",main="Puromyin experiment")
lines(X,(X*beta1)/(X+beta2),col="red",lwd=2)
lines(X,est_cor,col="green",lwd=2,lty=2)
legend("bottomright",legend=c("1/Moyenne","Milieu de 1/IC"),col=c(2,"green"),lty=1:2)
ggplot(data3, aes(y = Y, x = X)) +
geom_point(size = 2) +
geom_line(aes(y = value, color = legende, linetype = legende),
size = 1) +
labs(x = "Substrat", y = "Vitesse de réaction", title = "Données transformées")+
theme_bw()+
theme(legend.justification = c(1, 0),
legend.position = c(1, 0),
legend.background = element_rect(colour = "black"))
```
## Estimation d'un modèle de régression non linéaire par la méthode des moindres carrés
On minimise en fonction de $\theta$, la fonction $S(\theta)$:
On minimise en fonction de $\beta$, la fonction $S(\beta)$:
$$S(\theta) = \sum_{i=1}^n (y_i - f(x_i,\theta))^2$$
$$S(\beta) = \sum_{i=1}^n (y_i - f(x_i,\beta))^2$$
Malheureusement, il existe parfois des minima locaux en plus du minimum global en régression non linéaire. Malgré tout, comme dans le cas linéaire, on dérive par rapport à $\theta$ et on cherche les valeurs du vecteur $\theta$ qui annule la dérivée:
Malheureusement, il existe parfois des minima locaux en plus du minimum global en régression non linéaire. Malgré tout, comme dans le cas linéaire, on dérive par rapport à $\beta$ et on cherche les valeurs du vecteur $\beta$ qui annule la dérivée:
$$\sum_{i=1}^n (y_i - f(x_i,\theta)) \frac{\delta f(x_i,\theta)}{\delta \theta_r} \mid_{\hat{\theta}} \; =0 \; r=1,\ldots,p$$
$$\sum_{i=1}^n (y_i - f(x_i,\beta)) \frac{\delta f(x_i,\beta)}{\delta \beta_r} \mid_{\hat{\beta}} \; =0 \; r=1,\ldots,p$$
où $p$ est le nombre de paramètres (longueur du vecteur $\theta$).
où $p$ est le nombre de paramètres (longueur du vecteur $\beta$).
On obtient un système de $p$ équations qui n'a pas de solution analytique. Un algorithme iteratif est utilisé généralement.
### Principe de l'algorithme du gradient
Le principe de cet algorithme est de partir d'un point initial $\theta^0$ et de se laisser glisser vers le point le plus bas. On calcule la pente du point initial $S'(\theta^0)$. Si cette pente est positive (ça monte), alors on recule d’un pas. Si cette pente est négative, alors on avance d’un pas.
$$\theta^1 = \theta^0 - \text{pas} \times S'(\theta_0)$$
Le principe de cet algorithme est de partir d'un point initial $\beta^0$ et de se laisser glisser vers le point le plus bas. On calcule la pente du point initial $S'(\beta^0)$. Si cette pente est positive (ça monte), alors on recule d’un pas. Si cette pente est négative, alors on avance d’un pas.
$$\beta^1 = \beta^0 - \text{pas} \times S'(\beta_0)$$
Illustrons cet algorithme avec la fonction simplifiée suivante $S(\theta) = \theta^2/2$. Dans ce cas, la dérivée de $S(\theta)$ par rapport à $\theta$ vaut $S'(\theta)=\theta$ et $\theta^1 = \theta^0 - \text{pas} \times \theta_0$
Illustrons cet algorithme avec la fonction simplifiée suivante $S(\beta) = \beta^2/2$. Dans ce cas, la dérivée de $S(\beta)$ par rapport à $\beta$ vaut $S'(\beta)=\beta$ et $\beta^1 = \beta^0 - \text{pas} \times \beta_0$
```{r, fig.height=3}
theta <- seq(-2,2,by=0.1)
......@@ -147,23 +171,23 @@ theta_1
Il peut se schématiser comme suit :
1. Donner une valeur initiale $\theta^0$
2. Calculer une nouvelle valeur pour $\theta$ qui minimise approximativement $S(\theta)$, soit
$$\theta^1 = \theta^0 + (F(\theta^0)'F(\theta^0))^{-1}F(\theta^0)'(\mathbf{y}-f(\mathbf{x},\theta^0))$$
1. Donner une valeur initiale $\beta^0$
2. Calculer une nouvelle valeur pour $\beta$ qui minimise approximativement $S(\beta)$, soit
$$\beta^1 = \beta^0 + (F(\beta^0)'F(\beta^0))^{-1}F(\beta^0)'(\mathbf{y}-f(\mathbf{x},\beta^0))$$
où $F(\theta^0)$ est la matrice $n \times p$ des dérivées de $f(x_i,\theta)$ par rapport à $\theta$ prises au point $\theta^0$.
où $F(\beta^0)$ est la matrice $n \times p$ des dérivées de $f(x_i,\beta)$ par rapport à $\beta$ prises au point $\beta^0$.
3. Donner la valeur de $\theta_1$ à $\theta_0$ et recommencer jusqu'à convergence.
3. Donner la valeur de $\beta^1$ à $\beta^0$ et recommencer jusqu'à convergence.
D'autres algorithmes existent dont celui de Newton qui est basé sur un développement de Taylor à l'ordre 2 (quadratique) de $S(\theta)$. Ces algorithmes sont convergents pour des fonctions convexes, sous réserve que le point initial ne se trouve pas trop loin du minimum global (sinon, l'algorithme peut rester bloqué dans un minimum local).
D'autres algorithmes existent dont celui de Newton qui est basé sur un développement de Taylor à l'ordre 2 (quadratique) de $S(\beta)$. Ces algorithmes sont convergents pour des fonctions convexes, sous réserve que le point initial ne se trouve pas trop loin du minimum global (sinon, l'algorithme peut rester bloqué dans un minimum local).
Sur les logiciels comme $R$, on doit préciser une valeur initiale pas trop mauvaise pour avoir convergence de l'algorithme vers le vrai minimum global et on peut aussi donner la précision minimum pour le pas de l'itération. Enfin, cet algorithme itératif ne sera performant que si le modèle est intrinsèquement linéaire (au sens de Bates et Watts).
### Propriété des estimateurs
- Estimateur $\hat{\theta}$: Si $\varepsilon_i \sim N(0,\sigma)$ et les $\varepsilon_i$ sont indépendants, alors pour $n$ suffisamment grand on a approximativement $\hat{\theta} - \theta^* \sim N(0,\sigma^2 (F(\theta^*)'F(\theta^*))^{-1} )$ où $\theta^*$ est la vraie valeur inconnue de $\theta$.
- Estimateur $\hat{\beta}$: Si $\varepsilon_i \sim N(0,\sigma)$ et les $\varepsilon_i$ sont indépendants, alors pour $n$ suffisamment grand on a approximativement $\hat{\beta} - \beta^* \sim N(0,\sigma^2 (F(\beta^*)'F(\beta^*))^{-1} )$ où $\beta^*$ est la vraie valeur inconnue de $\beta$.
- Estimateur $\hat{\sigma}^2$: Sous les mêmes conditions, $\hat{\sigma}^2 = S(\hat{\theta})/(n-p)$ est tel que $(n-p)\hat{\sigma}^2/\sigma^2 \sim \chi^2_{n-p}$.
- Estimateur $\hat{\sigma}^2$: Sous les mêmes conditions, $\hat{\sigma}^2 = S(\hat{\beta})/(n-p)$ est tel que $(n-p)\hat{\sigma}^2/\sigma^2 \sim \chi^2_{n-p}$.
## Mise en oeuvre avec R
......@@ -179,7 +203,7 @@ summary(don2)
Le modèle proposé par Pellegrino, Gozé, Lebon et Wery (2006) pour relier l'ouverture des stomates (gs.rel) au stress hydrique (FTSW) est en fait non linéaire :
$$E(gs.rel) = f(FTSW,\theta) = \frac{a}{1+exp[-b_0(FTSW - b_1)]}$$.
$$E(gs.rel) = f(FTSW,\beta) = \frac{a}{1+exp[-b_0(FTSW - b_1)]}$$.
Mais, cette relation peut se linéariser comme suit :
$$\ln \left( \frac{E(gs.rel)/a}{1-E(gs.rel)/a} \right) = -b_0b_1 + b_0 FTSW$$
......@@ -248,7 +272,7 @@ legend("topleft", legend=c("Initiale","Finale"), col=c(2,"green"), lty=2:1)
### Estimation en utilisant les fonctions prédéfinies de R
Les auteurs du package _nlme_ ont créé des fonctions prédéfinies pour faciliter l'utilisation et l'initialisation de la régression non linéaire avec la fonction _nls_. Notamment, pour notre exemple on peut utiliser la fonction _SSlogis_ qui correspond à:
$$f(x_i,\theta) = \frac{\text{Asym}}{1+e^{(\text{xmid} - x_i)/ \text{scal}}}$$
$$f(x_i,\beta) = \frac{\text{Asym}}{1+e^{(\text{xmid} - x_i)/ \text{scal}}}$$
associée à la commande _getInitial_ pour obtenir des valeurs initiales pour les paramètres. Attention, dans notre exemple, le modèle utilisé par Pellegrino _et al_ n'a pas la même paramétrisation. Il faut ensuite convertir, avant d'utiliser les valeurs dans le start de _nls_.
```{r}
param_SS <- getInitial(gs.rel~SSlogis(don2$FTSW, Asym, xmid, scal), data = don2)
......
......@@ -39,7 +39,7 @@ On constate, d'après le graphique, que plus la valeur du PARint est élevée, p
On veut modéliser
$$y_i=\alpha+\beta x_i+\varepsilon_i$$
$$y_i=\mu+b_1 x_i+\varepsilon_i$$
où $\varepsilon_i \sim N(0,\sigma^2 \, \omega_i^2)$,
$y_i$ est la biomasse en $g.m^{-2}$, $x_i$ est le PARint en $MJ.m^{-2}$.
......
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