Commit 905453e2 authored by Isabelle Sanchez's avatar Isabelle Sanchez
Browse files

test sans fichier 05_applications

parent 5a1d6a49
Pipeline #41554 failed with stage
in 16 minutes and 46 seconds
# Applications
```{r, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE, cache=FALSE,eval=TRUE)
suppressPackageStartupMessages(library(nlme))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(sandwich))
suppressPackageStartupMessages(library(msm))
suppressPackageStartupMessages(library(GGally))
```
## GLM - variance hétérogène
Le but de cet exercice est de déterminer la relation qui existe entre la **biomasse** d'une plante et son **efficience biologique** (RUE ou capacité à convertir une lumière interceptée en biomasse). Nous allons tout d'abord observer nos données.
### Chargement et description des données:
On charge les données de _RUE.csv_.
```{r}
RUE <- read.csv("data/RUE.csv", header=TRUE, sep=";", dec =".")
summary(RUE)
```
On va alors **représenter les données**
```{r}
library(ggplot2)
ggplot(data=RUE, aes(y=Biomass_g_m2, x=PARint_MJ_m2)) +
geom_point() +
geom_smooth(method = "lm")
```
On constate, d'après le graphique, que plus la valeur du PARint est élevée, plus la variance est grande. Cela signifie que la variance des résidus est hétérogène. Pour représenter la relation entre les deux variables Biomasse et Parint, on devra non pas utiliser la _regression linéaire simple_, mais la **régression linéaire pondérée**.
### Modèle
On veut modéliser
$$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}$.
On pose $\omega_i = x_i^a$, on veut connaitre la puissance $a$ du poids $\omega_i$.
```{r}
library(nlme)
RUE <- dplyr::rename(RUE,PARint= PARint_MJ_m2 ,Biomass = Biomass_g_m2 )
reg3 <-gls(Biomass~PARint,weights=varPower(form=~PARint),data=RUE)
summary(reg3)
```
La puissance du poids est $\hat{a} = 0.74$ et notre équation de régression est:
$$\hat{y}_i= -2.53 + 1.46 x_i$$
### Validation du modèle
Le graphique des résidus standardisés montrent une stabilisation de la variance, mais une surestimation du modèle pour des faibles valeurs de biomasse. Le modèle linéaire n'est peut être pas le plus adapté ou il manque une variable explicative en plus de _Parint_.
```{r}
# graphiques des résidus
ggplot(RUE,aes(x=predict(reg3),y=residuals(reg3,type="normalized"))) +
geom_point()+
geom_abline(intercept=0,slope=0,col="red") +
xlab("Prédictions") + ylab("Normalized residuals")
```
### Conclusion
Pour conclure, comparons les prédictions d'une régression linéaire simple avec celles d'une régression pondérée.
```{r}
reg <- lm(Biomass~PARint,data=RUE)
ggplot(data=RUE, aes(y=Biomass, x=PARint)) +
geom_point() +
geom_abline(intercept=coef(reg3)[1],slope=coef(reg3)[2],col="magenta",lwd=1.2) +
geom_abline(intercept=coef(reg)[1],slope=coef(reg)[2],col="green",lwd=1.2)
```
Dans les 2 cas, on sur-estime la Biomasse pour de faibles valeurs de _Biomass_ (ou de _Parint_). Le modèle linéaire n'est pas pertinent.
## GLM - loi de Poisson
### Importation et mise en forme des données
Des lignées de sojas sont mises en présence d’une source de pucerons (toutes les lignées à égale distance de la source). Après 14 jours, une notation de l’infestation est réalisée sur chaque plante, selon une échelle de 0 à 6.
La variable d'intérêt est la variable _rating_ qui désigne le score (entre 0 et 6) d'infestation par les pucerons.
- Importation des données (fonction Import Dataset)
```{r}
mydata <- read.csv2("data/data_aphids.csv",stringsAsFactors = TRUE)
```
- Mise en forme d'une variable: _test_ et _rep_ sont des variables qualitatives, mais considérées comme des variables quantitatives par R. Comme ces variables n'ont pas été importées avec le bon type, il faut les transformer en variable qualitative:
```{r}
mydata$test <- factor(mydata$test,levels=1:25,labels=paste0("t", 1:25))
mydata$rep <- factor(mydata$rep,levels=1:8,labels=paste0("r", 1:8))
```
Vérification que la transformation a bien fonctionné:
```{r}
summary(mydata)
#aphid2 <- na.omit(mydata)
#summary(aphid2)
```
On constate avec le _summary_ qu'il y 2 données manquantes dans nos données pour la variable _rating_.
### Description de la covariable (ou variable explicative)
La covariable qualitative que nous allons utiliser est _acno_ le numéro de l'accession (ou variété de soja).
```{r}
#selection <- aphid2$acno %in% c("507776","507778","507771")
aphid2 <- na.omit(mydata)
selection <- aphid2$acno %in% levels(aphid2$acno)[1:20]
aphid3 <- aphid2[selection,]
aphid3$acno <- factor(aphid3$acno)
plot(rating~acno,data=aphid3)
```
```{r}
ggplot(mydata, aes(x=acno, y=rating)) +
stat_summary(fun=mean,na.rm=TRUE)
```
Lorsqu'on observe nos données, on se rend compte que le score de pucerons semble lié à la variété de soja.
Certaines variétés de soja ont des effectifs faibles (8 points de mesure).
```{r}
min(table(mydata$acno)) # ne pas utiliser le summary
```
### Analyse descriptive de la variable d'intérêt
Le modèle de régression de Poisson suppose que la valeur moyenne et la variance du score de pucerons est la même (approximativement), conditionnellement à la variété.
Affichons un résumé statistique par type de variété:
```{r}
mydata %>% group_by(acno) %>% summarise(avg=mean(rating,na.rm=TRUE), variance = var(rating,na.rm=TRUE))
```
Pour certaines variétés, les faibles effectifs ont conduit à une variance estimée nulle !
Pour mieux visualiser le lien entre la distribution du nombre de pucerons prédatés et la session, réalisons un histogramme conditionnel.
```{r, fig.cap="Histogramme du score de puceron observé par variété"}
monjeux <- levels(mydata$acno)[table(mydata$acno)>30]
ggplot(mydata[mydata$acno %in% monjeux,], aes(rating, fill = acno)) +
geom_histogram(breaks=c(0,1,2,3,4,5,6), position="dodge")
```
Certains des histogrammes obtenus ressemblent à des histogrammes que pourraient produire une loi de Poisson (décroissant comme _549046_ ou unimodaux asymétriques). Mais certains posent question (comme le _464889_ ou _549035_) avec une forme bimodale ou atypique (comme le _522212_).
### Mise en oeuvre d'une régression de type Poisson
Maintenant, nous allons réaliser une régression de Poisson sous _R_ à l'aide de la fonction _glm_.
```{r}
m1 <- glm(rating ~ acno , family="poisson", data=mydata)
summary(residuals(m1,type="deviance"))
ggplot(na.omit(mydata), aes(y=residuals(m1,type="deviance"), x=predict(m1, type="response"),col=rep)) +
geom_point()
```
L'information sur les résidus de la déviance montrent qu'il y a plutôt symétrie : la valeur médiane est nulle et les quartiles et les extrèmes plutôt symmétriques.
Nous effectuons un test statistique sur les déviances résiduelles afin de savoir si elles sont faibles et ainsi vérifier que le modèle ajuste correctement les données.
```{r}
#test statistique associé à la déviance résiduelle
with(m1, cbind(res.deviance = deviance, df = df.residual,
p = pchisq(deviance, df.residual, lower.tail=FALSE)))
```
De plus,
```{r}
## test model differences with chi square test
anova(m1, test="Chisq")
```
On remarque que la Residual deviance (comparaison de notre modèle avec le modèle saturé) est inférieure à la Null deviance (comparaison de notre modèle avec un modèle où il n'y aurait que l'intercept sans les variables explicatives). Notre variable explicative variété (_acno_) a donc bien un rôle à jouer.
### Conclusion sur la régression de Poisson
Si on regarde la prédiction (et que l'on compare aux observations). On constate que les variétés dont l'histogramme était biomodal donnent des prédictions discutables.
```{r, fig.cap="Histogramme du score de puceron prédit par variété"}
newdata <- cbind(na.omit(mydata),predict=predict(m1,type="response"))
monjeux <- levels(newdata$acno)[table(newdata$acno)>30]
ggplot(newdata[newdata$acno %in% monjeux,], aes(predict, fill = acno)) +
geom_histogram(breaks=c(0,1,2,3,4,5,6), position="dodge")
```
```{r}
ggplot(aphid2,aes(y=residuals(m1,type="deviance"),x=predict(m1,type="response"),col=design)) +
geom_point()
```
## NLM
### Importation et description des données
```{r}
# importation des données
don <- read.csv2(file="data/reg_non_lin_FintPois_chiche.csv", dec=".")
```
- Résumé des données
```{r}
summary(don)
```
Notre variable d'intérêt est le _DAS_ que l'on va étudier en fonction du temps thermique (ou _Thermal.time_) ou de _f_. La relation sera surement différente selon les génotypes. Nous allons supposer que la forme de la relation reste la même (fonction logistique), mais la valeur des coefficients/paramètres varient en fonction du génotype.
- Analyse graphique
```{r}
GGally::ggpairs(na.omit(don[,5:7]))
```
Nous reprenons le schéma d'analyse proposé en cours:
1. Estimons une valeur pour l'asymptote, notée $aE$, visuellement à partir du graphique ou numériquement à partir des données.
```{r}
geno1 <- don[don$Geno==1,]
# Estimation de valeurs initiales : linéarisation du modèle
aE <- max(geno1$f,na.rm=TRUE)+0.0001
print(aE)
```
2. Estimons $b_0$ comme le coefficient de pente de la régression de $\ln(DAS/(aE-DAS))$ en fonction de _f_.
```{r}
res.prepa <- lm(log(f/(aE-f))~DAS,data=geno1)
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(DAS/(aE-DAS))$ en fonction de _f_.
```{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)))}
plot(f~DAS,data=geno1,ylim=c(0,max(geno1$f,na.rm=TRUE)),
ylab="f",pch=15)
lines(geno1$DAS[order(geno1$DAS)],
modele1(geno1$DAS[order(geno1$DAS)], aE, b0E, b1E), col="red", lty=2)
```
```{r}
init <- getInitial(f~SSlogis(DAS,Asym,xmid,scal),data=don)
plot(f~DAS,data=geno1,ylim=c(0,max(geno1$f,na.rm=TRUE)),
ylab="f",pch=15)
lines(geno1$DAS[order(geno1$DAS)],
SSlogis(geno1$DAS[order(geno1$DAS)], init[1], init[2], init[3]), 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
res.nls <- nls(f ~ a/(1+exp(-b0*(DAS-b1))), start = list(a = aE,b1 = b1E, b0 =b0E),
data=geno1)
summary(res.nls)
```
6. Rajouter la courbe estimée par la fonction _nls_ au graphique précédent.
```{r}
plot(f~DAS, data=geno1, ylim = c(0, max(geno1$f)),
ylab="f", pch=15)
lines(geno1$DAS[order(geno1$DAS)], modele1(geno1$DAS[order(geno1$DAS)],
aE, b0E, b1E), col="red", lty=2)
lines(geno1$DAS[order(geno1$DAS)], modele1(geno1$DAS[order(geno1$DAS)],
coef(res.nls)[1], coef(res.nls)[3],coef(res.nls)[2]), col="green", lwd=2)
legend("topleft", legend=c("Initiale","Finale"), col=c(2,"green"), lty=2:1)
```
### Validation du modèle : inspection visuelle des résidus
```{r,fig.cap="Régression non linéaire de f en fonction du DAS: graphique des résidus."}
plot(predict(res.nls),residuals(res.nls,type="pearson"), ylab="Prédictions", xlab="Résidus Pearson")
abline(h=0,col="green",lwd=2)
```
### Modèle à coefficients variables
```{r}
don$Geno <- factor(don$Geno)
plot(f~DAS, data=don, ylim = c(0, max(don$f,na.rm=TRUE)),
ylab="f", pch=15,col=don$Geno)
legend("topleft",legend=levels(don$Geno),lty=2:6,col=1:nlevels(don$Geno))
```
Supposons que les paramètre $a$, $b_0$ et $b_1$ dépendent du génotype, on change notre modèle:
```{r}
for (i in 1:nlevels(don$Geno)) {
res.nls <- nls(f ~ a/(1+exp(-b0*(DAS-b1))), start = list(a = aE,b1 = b1E, b0 =b0E),
data=na.omit(don[don$Geno==i,]))
print(coef(res.nls))
}
```
$$f_{ij} = \frac{a_{i}}{1+exp[-b_{0, \, i}(DAS_{ij} - b_{1, \, i})]} + \varepsilon_{ij}$$
```{r}
don.reduit <- na.omit(don[don$Geno %in% c("1","2","3","4","5","6","7","8"),])
don.reduit$Geno <- factor(don.reduit$Geno)
res.prepa <- lm(log(f/(aE-f))~-1 + DAS*Geno,data=don.reduit)
b0E <- as.numeric(coef(res.prepa)[1] +
c(0,coef(res.prepa)[(2+nlevels(don.reduit$Geno)):(2*nlevels(don.reduit$Geno))]))
print(b0E)
b1E <- -as.numeric(coef(res.prepa)[2:(1+nlevels(don.reduit$Geno))]/b0E)
print(b1E)
res.gnls2 <- gnls(f ~ a/(1+exp(-b0 * (DAS-b1))), params=list(a~Geno,b0~ Geno, b1 ~ Geno),
start = list(a = rep(1,nlevels(don.reduit$Geno)), b0 = b0E , b1 = b1E), data=don.reduit)
coef(res.gnls2)
```
D'après le critère d'AIC, le modèle retenu est celui avec le plus petit AIC.
```{r}
AIC(res.nls,res.gnls2)
```
Les résidus du modèle
```{r,fig.cap="Régression non linéaire de f en fonction de DAS et du génotype: graphique des résidus."}
plot(predict(res.gnls2),residuals(res.gnls2,type="pearson"),ylab="Prédictions",
xlab="Résidus Pearson")
abline(h=0,col="green",lwd=2)
```
......@@ -10,7 +10,6 @@ rmd_files: [
"02-GLM-variance-reduit.Rmd",
"03-GLM-loi.Rmd",
"04-NLR.Rmd",
"05-Applications.Rmd",
"06-summary.Rmd",
"07-references.Rmd",
]
# Exercice - Régression non linéaire
```{r, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
suppressPackageStartupMessages(library(nlme))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))
```
## Bibliographie
tiré de Amani Louis Kouadio, Bakary Djaby, Grégory Duveiller, Moussa El Jarroudi & Bernard Tychon
_Cinétique de décroissance de la surface verte et estimation du rendement du blé d'hiver_
Pour caractériser l'évolution temporelle du LAI (calculé à partir d'indices de végétation) d'un couvert de blé, Baret et al. (1986) ont proposé un modèle paramétrique combinant un modèle logistique simple modélisant la phase de croissance de la culture et un modèle exponentiel négatif caractérisant la phase de sénescence de la plante. Gooding et al. (2000) ont montré que la sénescence des dernières feuilles émergées du blé est bien caractérisée par une fonction Gompertz modifiée. Étant donné que l'objectif de notre étude est de valoriser une des conclusions de l'étude précédente, et pour apprécier leur capacité à décrire l'évolution décroissante du GAI chez le blé, les fonctions Gompertz modifiée (Équation 1) et logistique modifiée (Équation 2) proposées par Gooding et al. (2000) seront utilisées pour décrire la phase de décroissance du GAI.
![](img-5.png)
où A désigne le pourcentage initial du GAI, m (exprimé en degrés par jour) est le moment où le point d'inflexion de la courbe apparait, k exprime le taux de dégénérescence du GAI et tt, le temps thermique exprimé en degrés-jours.
Lorsque A est fixé à 100 % (GAI est à sa valeur maximale), l'inflexion de la courbe lors de la phase décroissante apparait à GAI = 37 % pour la formule Gompertz modifiée (50 % pour la fonction logistique modifiée). Le temps thermique, tt, sur n jours se calcule comme suit :
![](img-6.png)
où Tai (°C·j-1) est égale à la moyenne journalière de la température de l'air. [Tai = 0,5*(Tmax + Tmin)] ; Tb (°C·j-1) étant égale à la température de base. Dans le calcul des cumuls, (Tai - Tb) est égal à 0 lorsque Tai < Tb. Pour le blé, Tb est fixée à 0 °C (Baker et al., 1980 ; Baker et al., 1983 ; Davidson et al., 1983 ; Masle et al., 1989). La période de calcul du temps thermique commence à la valeur maximale de GAI (GAImax). Il est alors fixé à 0 à cette date et le cumul est effectué pour les jours suivants. Ainsi, le metric m représente le temps thermique nécessaire pour atteindre soit 37 % (Équation 1), soit 50 % (Équation 2) de la valeur de GAI restant depuis sa valeur maximale. Les variables climatiques, températures minimale (Tmin) et maximale (Tmax), sont issues de stations météorologiques proches des parcelles suivies (entre 1 et 5 km).
## Importation des données
```{r}
mydata <- read.csv2("data/reg_non_lin_GPAI.csv", dec=".",stringsAsFactors = TRUE)
summary(mydata)
```
## Exploration du jeu de données
```{r}
ggplot(mydata,aes(x=GPAI,y=DAP,col=densite)) +
geom_point()
```
## Modélisation
### Valeurs initiales
```{r}
aE <- max(mydata$DAP,nr.rm=TRUE)+0.0001
res.prepa <- lm(log(DAP/(aE-DAP))~GPAI,data=mydata)
b0E <- as.numeric(coef(res.prepa)[2])
b1E <- -as.numeric(coef(res.prepa)[1]/coef(res.prepa)[2])
```
### Modèle logistique
```{r, fig.height=4}
res.nls <- nls(DAP ~ a/(1+exp(-b0*(GPAI-b1))), start = list(a = aE,b1 = b1E, b0 =b0E),
data=mydata)
summary(res.nls)
library(nlstools)
plotfit(res.nls,variable=4)
```
Nous allons représenter les résidus du modèle:
```{r}
plot(nlsResiduals(res.nls))
```
On constate une dépendance dans les résidus (graphique d'autocorrélation). Ce phénomène est fréquent pour des données longitudinales. Nous allons rajouter au modèle une structure de corrélation de type Auto Régressif d'ordre 1. On suppose que pour une parcelle donnée (_id_parcelle_) les données sont corrélées dans le temps selon le modèle suivant:
$$cor(DAP_{t_i},DAP_{t_j}) = \rho^{|t_j-t_i|}$$
Normalement, ce type de modèle dans son implémentation dans le package _nlme_ suppose que les mesures sont espacées du même laps de temps, ce qui n'est pas tout à fait le cas:
```{r}
unique(mydata$date)
```
Malgré tout, nous testons ce modèle sur nos données pour tenir compte de l'autocorrélation détectée dans les résidus:
```{r}
res.gnls <- gnls(DAP ~ a/(1+exp(-b0 * (GPAI-b1))), start = list(a = aE,b1 = b1E, b0 = b0E),
data=mydata, correlation=corAR1(form=~1|id_parcelle),
control = list(tolerance=10)) #Augmentation de la tol?rance pour que cela fonctionne
summary(res.gnls)
AIC(res.nls,res.gnls)
plot(predict(res.gnls),residuals(res.gnls,type="normalized"),ylab="Prédictions",
xlab="Résidus Normalisés",col=mydata$date)
abline(h=0,col="green",lwd=2)
```
Le modèle ajusté et les données ne montrent pas une meilleure adéquation, même si le critère AIC nous incité à choisir le modèle avec la modélisation de l'autocorrélation.
```{r}
modele1 <- function(x,a,b0,b1){a/(1+exp(-b0*(x-b1)))}
plot(DAP~GPAI, data=mydata, ylim = c(0, max(mydata$DAP)),
ylab="DAP", pch=15,col=mydata$densite)
# densite1
lines(mydata$GPAI[order(mydata$GPAI)], modele1(mydata$GPAI[order(mydata$GPAI)],
coef(res.gnls)[1], coef(res.gnls)[3],coef(res.gnls)[2]), col=1, lwd=2)
```
### Modèle logistique généralisé
```{r}
res.prepa <- lm(log(DAP/(aE-DAP))~-1 + GPAI*densite,data=mydata)
b0E <- as.numeric(coef(res.prepa)[1] +
c(0,coef(res.prepa)[(2+nlevels(mydata$densite)):(2*nlevels(mydata$densite))]))
print(b0E)
b1E <- -as.numeric(coef(res.prepa)[2:(1+nlevels(mydata$densite))]/b0E)
print(b1E)
res.gnls2 <- gnls(DAP ~ a/(1+exp(-b0 * (GPAI-b1))), params=list(a~1,b0~ 1, b1 ~ densite),
start = list(a = aE, b0 = b0E[1] , b1 = b1E), data=mydata)
summary(res.gnls2)
```
```{r}
modele1 <- function(x,a,b0,b1){a/(1+exp(-b0*(x-b1)))}
plot(DAP~GPAI, data=mydata, ylim = c(0, max(mydata$DAP)),
ylab="DAP", pch=15,col=mydata$densite)
# densite1
lines(mydata$GPAI[order(mydata$GPAI)], modele1(mydata$GPAI[order(mydata$GPAI)],
coef(res.gnls2)[1], coef(res.gnls2)[2],coef(res.gnls2)[3]), col=1, lwd=2)
# densite2
lines(mydata$GPAI[order(mydata$GPAI)], modele1(mydata$GPAI[order(mydata$GPAI)],
coef(res.gnls2)[1],coef(res.gnls2)[2],
(coef(res.gnls2)[3] + coef(res.gnls2)[4])), col=2, lwd=2)
# densite3
lines(mydata$GPAI[order(mydata$GPAI)], modele1(mydata$GPAI[order(mydata$GPAI)],
coef(res.gnls2)[1], coef(res.gnls2)[2],
(coef(res.gnls2)[3] + coef(res.gnls2)[5])), col=3, lwd=2)
legend("bottomright",legend=c("densite1","densite2","densite3"),lty=2:6,col=1:3)
```
Le graphique précédent montre un meilleur ajustement avec le modèle à coefficients variables. On a l'impression visuelle qu'en plus du paramètre $b_1$, le paramètre $b_0$ lié à la pente-tangente au point d'inflexion $(b_1/2,a)$ est variable selon la densité.
```{r}
res.gnls3 <- gnls(DAP ~ a/(1+exp(-b0 * (GPAI-b1))),
start = list(a = aE, b0 = b0E, b1 = b1E),
params=list(a~1,b0~densite, b1~densite),
data=mydata, control = list(tolerance=10))
summary(res.gnls3)
```
## Validation
D'après le critère d'AIC, le modèle retenu est celui avec le plus petit AIC.
```{r}
AIC(res.nls,res.gnls2, res.gnls3)
```
Au sens du critère de l'AIC, le meilleur modèle est le modèle avec les coefficients $b_0$ et $b_1$ variables. Les résidus de ce modèle montrent toujours une structure en fonction de la date. Il faudrait tenir compte de l'autocorrélation dans le jeu de données mais avec un modèle plus pertinent que l'AR1 (car les mesures ne sont pas espacées régulièrement dans le temps).
```{r,fig.cap="Régression non linéaire de f en fonction de DAS et du génotype: graphique des résidus."}
plot(predict(res.gnls3),residuals(res.gnls3,type="pearson"),ylab="Prédictions",
xlab="Résidus Pearson",col=mydata$date)
abline(h=0,col="green",lwd=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