Commit 308f52e6 authored by sanchezi's avatar sanchezi
Browse files

correc 04

parent 5e8edc93
Pipeline #42105 failed with stage
in 14 minutes and 38 seconds
......@@ -17,7 +17,6 @@ image: rocker/verse:4.0.4
- Rscript -e "install.packages(c('pscl'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('DHARMa'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('ade4'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('car'), repos='https://cloud.r-project.org')"
- Rscript -e "bookdown::render_book('index.Rmd', 'all', output_dir = 'public')"
artifacts:
paths:
......
......@@ -41,6 +41,7 @@ Linéarisable en
$$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}
# Données de puromycin
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)
......@@ -58,8 +59,8 @@ ggplot(data, aes(y = 1/Y, x = 1/X)) +
theme_bw()
```
2. Effectuons la régression linéaire sur données transformées et ajoutons la droite de régression au graphique précédent.
```{r, fig.height=3}
# régression sur données transformées
res.l <- lm(1/Y~I(1/X)) # modèle linéarisé
......@@ -73,13 +74,14 @@ ggplot(data, aes(y = 1/Y, x = 1/X)) +
theme_bw()
```
3. Observons la sortie du _summary_.
```{r}
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
b1 <- 1/coef(res.l)[1]
......@@ -136,6 +138,7 @@ On obtient un système de $p$ équations qui n'a pas de solution analytique. Un
### Principe de l'algorithme du gradient
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(\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$
......@@ -189,7 +192,6 @@ Sur les logiciels comme $R$, on doit préciser une valeur initiale pas trop mauv
- 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
```{r}
......@@ -213,28 +215,30 @@ $$\ln \left( \frac{E(gs.rel)/a}{1-E(gs.rel)/a} \right) = -b_0b_1 + b_0 FTSW$$
Pour l'estimation du modèle non linéaire, il faut fournir des valeurs initiales les plus pertinentes possibles. Quand un modèle est linéarisable, on utilise cette propriété pour estimer des valeurs initiales comme suit :
1. Estimons une valeur pour l'asymptote, notée $aE$, visuellement à partir du graphique ou numériquement à partir des données.
```{r}
# Estimation de valeurs initiales : linéarisation du modèle
aE <- max(don2$gs.rel)+0.0001
print(aE)
```
2. Estimons $b_0$ comme le coefficient de pente de la régression de $\ln(gs.rel/(aE-gs.rel))$ en fonction de FTSW.
```{r}
res.prepa <- lm(log(gs.rel/(aE-gs.rel))~FTSW,data=don2)
b0E <- as.numeric(coef(res.prepa)[2])
print(b0E)
```
3. Estimons $b_1$ à partir du rapport de l'intercept sur le coefficient de pente estimés par la régression de $\ln(gs.rel/(aE-gs.rel))$ en fonction de FTSW.
```{r}
b1E <- -as.numeric(coef(res.prepa)[1]/coef(res.prepa)[2])
print(b1E)
```
4. Représentons les données et la courbe avec les valeurs initiales des paramètres. L'ajustement n'est pas parfait mais pas trop éloigné de la réalité, les estimations précédentes pourront servir comme valeurs initiales.
```{r}
# graphique des données et du modèle ajusté par les estimations grossières
modele1 <- function(x,a,b0,b1){a/(1+exp(-b0*(x-b1)))}
......@@ -246,8 +250,8 @@ lines(don2$FTSW[order(don2$FTSW)],
modele1(don2$FTSW[order(don2$FTSW)], aE, b0E, b1E), col="red", lty=2)
```
5. Ajuster le modèle non linéaire avec la fonction _nls_. Commenter les résultats et le graphique des résidus.
```{r}
# function logistique et dérivée
formula_perso <- as.formula(gs.rel ~ a/(1+exp(-b0*(FTSW-b1))))
......@@ -257,6 +261,7 @@ summary(res.nls)
```
6. Rajouter la courbe estimée par la fonction _nls_ au graphique précédent.
```{r}
plot(gs.rel~FTSW, data=don2, asp=1, ylim = c(0, max(don2$gs.rel)),
ylab="Ouverture des stomates", pch=15)
......@@ -271,9 +276,11 @@ 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,\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)
param_SS
......@@ -294,24 +301,31 @@ abline(h=0,col="green",lwd=2)
```
Pour compléter cette validation, on peut utiliser le package _nlstools_ qui facilite la production de graphiques de validation.
```{r,fig.cap="Nuage de points et valeurs prédites, graphique réalisé avec nlstools", fig.height=3}
library(nlstools)
plotfit(res.nls,variable=6)
```
On peut produire des graphiques similaires à ceux de la régression linéaire. Dans notre exemple, l'inspection visuelle révèle une possible hétérogénéité de la variance (une variablité croissante des résidus en fonction de la valeur prédite).
```{r}
plot(nlsResiduals(res.nls))
```
Ce package permet aussi d'identifier les observations influentes pas Jacknife (méthode de validation croisée par leave one out)
Ce package permet aussi d'identifier les observations influentes pas Jacknife (méthode de validation croisée par leave one out
```{r}
plot(nlsJack(res.nls))
```
Le package permet aussi de calculer des intervalles de confiance par une approche bootstrap (tirage d'échantillon avec remise) en complément de la commande habituelle _confint_.
```{r}
confint(res.nls)
```
```{r, eval=FALSE}
summary(nlsBoot(res.nls))
```
......@@ -319,15 +333,19 @@ summary(nlsBoot(res.nls))
## Hétérogénéité de la variance
Au vu du graphique des résidus Fig. \@ref{fig:resnls}, on peut penser que la variance est hétérogène. On va modéliser la variance comme proportionnelle à une fonction puissance de la valeur prédite.
```{r}
res.gnls <- gnls(formula_perso, start = list(a = aE,b1 = b1E, b0 = b0E), data=don2,
weights=varPower())
summary(res.gnls)
```
Le graphique des résidus normalisés montre que l'impression d'hétérogénéité de la variance a été corrigée.
```{r,fig.cap="Régression non linéaire généralisée de l'ouverture des stomates (gs.rel) en fonction du stress hydrique (FTSW): graphique des résidus.", fig.height=2.5}
plot(predict(res.gnls),residuals(res.gnls,type="normalized"), ylab="Prédictions",
plot(predict(res.gnls),residuals(res.gnls,type="normalized"),
ylab="Prédictions",
xlab="Résidus Normalisés")
abline(h=0,col="green",lwd=2)
```
......@@ -335,6 +353,7 @@ abline(h=0,col="green",lwd=2)
## Modèle à coefficients variables (les paramètres peuvent dépendre linéairement d'une ou plusieurs covariables)
Dans notre exemple, l'impression d'hétérogénéité de la variance est peut être due au fait que la relation non linéaire est différente selon la localité. La forme globale reste la même, mais les valeurs des coefficients/paramètres changent selon la localité. Pour le vérifier,on commence par agglomérer les parcelles de Sauteyrargues sur les 3 années.
```{r}
don2$Label <- factor(substring(as.character(don2$Label),1,3))
```
......@@ -351,6 +370,7 @@ summary(res.gnls2)
```
On a ajusté une courbe pour chaque Label (localité).
```{r}
plot(gs.rel~FTSW, data=don2, asp=1, ylim = c(0, max(don2$gs.rel)),
ylab="Ouverture des stomates", pch=15,col=don2$Label)
......@@ -373,15 +393,19 @@ D'après le critère d'AIC, le modèle retenu est celui avec le plus petit AIC.
```{r}
AIC(res.nls,res.gnls,res.gnls2)
```
Les résidus du modèle
```{r,fig.cap="Régression non linéaire de l'ouverture des stomates (gs.rel) en fonction du stress hydrique (FTSW) et de la localité: graphique des résidus.", fig.height=3}
plot(predict(res.gnls2),residuals(res.gnls2,type="pearson"),ylab="Prédictions",
plot(predict(res.gnls2),residuals(res.gnls2,type="pearson"),
ylab="Prédictions",
xlab="Résidus Pearson")
abline(h=0,col="green",lwd=2)
```
- On combine les deux approches: coefficients variables et modélisation de la variance
```{r}
res.gnls3 <- gnls(gs.rel ~ a/(1+exp(-b0 * (FTSW-b1))), params=list(a~1,b0~1, b1 ~ Label),
start = list(a = aE, b0 = b0E, b1 = rep(b1E,3)),
......
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