Commit 0609ce3d authored by sanchezi's avatar sanchezi
Browse files

new deploy - bug free

parent 308f52e6
Pipeline #42338 failed with stage
in 17 minutes and 29 seconds
......@@ -3,21 +3,27 @@ image: rocker/verse:4.0.4
.bookdown:
stage: deploy
script:
- Rscript -e "install.packages(c('dplyr'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('gstat'), repos='https://cloud.r-project.org')"
- apt-get update
- apt-get install -y libblas-dev
- apt-get install -y liblapack-dev
- Rscript -e "install.packages(c('ggplot2'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('sp'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('nlme'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('tidyverse'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('dplyr'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('ade4'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('sandwich'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('msm'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('nlstools'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('kableExtra'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('tidyverse'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('DHARMa'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('AER'), repos='https://cloud.r-project.org')"
- 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 "bookdown::render_book('index.Rmd', 'all', output_dir = 'public')"
- Rscript -e "install.packages(c('nlstools'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('gstat'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('cowplot'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('kableExtra'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('car'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('sp'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('nlme'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('scales'), repos='https://cloud.r-project.org')"
- Rscript -e "install.packages(c('ipred'), repos='https://cloud.r-project.org')"
artifacts:
paths:
- public
......
This diff is collapsed.
......@@ -29,6 +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 f01-exocor}
x <- sort(round(runif(15,max=1.5),2))
y1 <- exp(x)
......@@ -52,12 +53,13 @@ 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ées dans la Table ref{tab:f01-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 le tableau suivant:
```{r f01-importation}
vignes <- read.csv2("data/vignes.csv",stringsAsFactors = TRUE)
kable(summary(vignes),
booktabs=TRUE)
booktabs=TRUE,
caption="Résumé des données")
```
......@@ -68,7 +70,7 @@ $$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)$.
+$\varepsilon_i$ représente la variabilité des individus associés à $x_i$, souvent noté aléas individuel ou erreur de modélisation.
+ $\varepsilon_i$ représente la variabilité des individus associés à $x_i$, souvent noté aléas individuel ou erreur de modélisation.
- On suppose que les conditions suivantes sont vérifiées :
+ Indépendance, échantillonnage: Les paires d'observations $(Y_i,x_i)$ résultent de $n$ tirages indépendants et aléatoires (échantillonnage aléatoire simple) dans une population déterminée.
......@@ -124,7 +126,7 @@ confint(res.lin)
- Prédiction individuelle: D'après le modèle, une nouvelle observation
$Y_*$ pour $x_*$ s'écrit comme $Y_* = \hat{Y}_* + \hat{\varepsilon}_*$ et suit
une loi normale. On peut calculer un intervalle de confiance
$$ \hat{Y}_* \pm t_{\alpha,(n-2)} \sqrt{\hat{\sigma}^2 + \hat{\sigma}^2/n + (x_* - \bar{x})^2 (\hat{\sigma}^2/S^2_x)/(n-1)}$$
$$\hat{Y}_* \pm t_{\alpha,(n-2)} \sqrt{\hat{\sigma}^2 + \hat{\sigma}^2/n + (x_* - \bar{x})^2 (\hat{\sigma}^2/S^2_x)/(n-1)}$$
- Prédiction moyenne: D'après le modèle, une nouvelle observation moyenne
$\bar{Y}_*$ pour $x_*$ s'écrit comme $\bar{Y}_* = \hat{Y}_*$ et suit
......@@ -181,7 +183,7 @@ summary(res.lin)
Le coefficient de détermination est un critère pour juger de l'adéquation d'un modèle aux données. Il représente le pourcentage de variance expliquée par le modèle linéaire. Dans le cas d'un modèle avec intercept, il se calcule comme suit :
$$ R^2 = 1 - \frac{\sum_{i=1}^n (y_i -\hat{y}_i)^2}{\sum_{i=1}^n (y_i -\bar{y})^2} = 1 - \frac{\hat{\sigma}^2 \times (n-2)}{S^2_y \times (n-1)}$$
$$R^2 = 1 - \frac{\sum_{i=1}^n (y_i -\hat{y}_i)^2}{\sum_{i=1}^n (y_i -\bar{y})^2} = 1 - \frac{\hat{\sigma}^2 \times (n-2)}{S^2_y \times (n-1)}$$
Dans le cas de l'exemple, le $R^2$ est indiqué dans le summary sous
l'appellation Multiple R-squared et vaut 64.62\%. Ainsi, le modèle proposé n'explique qu'environ 65 \% de la variabilité observée sur l'ouverture des stomates et peut encore être amélioré.
......
......@@ -264,7 +264,7 @@ anova(res.cov,res.cov.pon)
Le Schéma des comparaisons de modèles possibles est le suivant
![alt text](SchemaTest.png)
![Schema de tests](SchemaTest.png)
......@@ -286,7 +286,8 @@ ggplot(vignes, aes(fitted(res.cov.pon), residuals(res.cov.pon, type="normalized"
geom_hline(yintercept=0,col="red") +
xlab("Fitted Value") +
ylab("Standardized residuals") +
ggtitle("Validation de l'ancova pondérée")+theme_bw()
ggtitle("Validation de l'ancova pondérée") +
theme_bw()
```
......@@ -300,22 +301,29 @@ où $\hat{y}_{i,(-i)}$ est l'estimation pour l'observation $i$ à partir d'un mo
**Remarque** : la validation croisée permet aussi d'enlever les individus par paquets.
```{r}
require(qpcR) # fonctionne avec les sorties de lm, glm, nls
res.lin <- lm(gs.rel~FTSW,data=vignes)
PRESS(res.lin,verbose=F)$stat
res.cov <- lm(gs.rel~FTSW*Label,data=vignes)
PRESS(res.cov,verbose=F)$stat
```
```{r,eval=FALSE}
library(qpcR) # fonctionne avec les sorties de lm, glm, nls
PRESS(res.lin,verbose=F)$stat
PRESS(res.cov,verbose=F)$stat
```
* RMSE et RMSEP: Root mean squared error ou racine carrée de l'erreur quadratique moyenne. Elle se calcule à partir de la formule suivante :
$$RMSEP = \sqrt{\sum_{i=1}^{n} \frac{1}{n}(y_i -\hat{y}_i)^2}$$
Afin de l'utiliser dans une perspective de prédiction (RMSEP), la formule s'applique sur des données nouvelles (dites de test) ou simulées et non sur le même jeu de données (dit d'entraînement) qui a servi à ajuster le modèle. On peut aussi utiliser la validation croisée pour calculer le RMSEP comme pour le PRESS. En général, le RMSEP est plus grand que le RMSE.
```{r}
```{r rmse,eval=FALSE}
# fonction de qpcR
RMSE(res.lin)
RMSE(res.cov)
```
```{r errortest}
library(ipred)
errorest(gs.rel ~ FTSW, data=vignes, model=lm,
est.para=control.errorest(k=20,random=FALSE))
......@@ -323,7 +331,7 @@ errorest(gs.rel ~ FTSW*Label, data=vignes, model=lm,
est.para=control.errorest(k=20,random=FALSE))$error
```
Il faut choisir le nombre de folders $k$ suffisament grand pour avoir dans chaque folder des individus des différents _Label_ pour l'analyse de la covariance (ancova).
Il faut choisir le nombre de folders $k$ suffisamment grand pour avoir dans chaque folder des individus des différents _Label_ pour l'analyse de la covariance (ancova).
### Comparaison de modèles emboîtés : critère du $R^2$
......
Cet enseignement utilise les packages suivants
```{r, include=TRUE}
# Introduction et Visualisation de données simulées avec une loi de Poisson
Cet enseignement utilise les packages suivants:
```{r f03-setting, include=TRUE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE,
message=FALSE, cache=FALSE,eval=TRUE)
suppressPackageStartupMessages(library(ggplot2))
......@@ -8,14 +12,12 @@ suppressPackageStartupMessages(library(msm))
suppressPackageStartupMessages(library(dplyr))
```
# Introduction et Visualisation de données simulées avec une loi de Poisson
## Histogramme d'une loi de Poisson
La loi de Poisson permet de modéliser la distribution d'une variable de comptage $Y$. Elle a comme propriété principale que son espérance est égale à sa variance $E(Y) = V(Y) =\lambda$. Sa loi ou distibution s'écrit comme:
La loi de Poisson permet de modéliser la distribution d'une variable de comptage $Y$. Elle a comme propriété principale que son espérance est égale à sa variance $E(Y) = V(Y) =\lambda$. Sa loi ou distribution s'écrit comme:
$$P(Y=y) = \frac{\lambda^y}{y!} \, e^{-y}$$
```{r,fig.height=3}
```{r f03-simu,fig.height=3}
y <- rep(c(0:10),3)
distribution <- c(dpois(y[1:11],1),dpois(y[1:11],4),dpois(y[1:11],10))
mu <- rep(c(1,4,10),each=11)
......@@ -25,7 +27,7 @@ mydata <- transform(data.frame(x=y,distribution=distribution,mu=factor(mu)))
dodge <- position_dodge(width = 0.1)
ggplot(mydata, aes(y, distribution,fill = mu)) +
geom_bar(stat="identity",position = position_dodge())+
geom_bar(stat="identity",position = position_dodge()) +
guides(fill = guide_legend(title =expression(mu))) +
theme_bw()
```
......@@ -34,7 +36,7 @@ ggplot(mydata, aes(y, distribution,fill = mu)) +
1. Relation décroissante où $Y \sim P(\lambda(x))$ et $\lambda(x)=e^{1-0,2x}$
```{r,fig.height=3}
```{r f03-simu2,fig.height=3}
Var_exp <- seq(0,50,by=1)
mu <- exp(1-0.2*Var_exp)
n=11
......@@ -50,28 +52,34 @@ ggplot(mydata,aes(x,y)) +
theme_classic()
```
En première approche, comme $E(y)=\lambda(x)=e^{1-0,2x}$, on peut décider d'utiliser la transformation logarithme ($\ln(y+1)$) pour estimer les paramètres de la relation entre $Y$ et $X$. On rajoute $+1$ pour éviter $\ln(0) = -\infty$ et de plus en décalant de $+1$ on conserve les observations nulles après transformation (car $\ln(0+1)=0$).Ensuite, on espère retrouver les valeurs qui ont servi à la simulation (intercept=1 et pente = -0,2) à partir d'une régression linéaire simple (classique):
```{r}
En première approche, comme $E(y)=\lambda(x)=e^{1-0,2x}$, on peut décider d'utiliser la transformation logarithme ($\ln(y+1)$) pour estimer les paramètres de la relation entre $Y$ et $X$. On rajoute $+1$ pour éviter $\ln(0) = -\infty$ et de plus en décalant de $+1$ on conserve les observations nulles après transformation (car $\ln(0+1)=0$). Ensuite, on espère retrouver les valeurs qui ont servi à la simulation (intercept=1 et pente = -0,2) à partir d'une régression linéaire simple (classique):
```{r f03-lm}
mydatalog <- data.frame(y_log=log(y.sim+1), Var_exp=Var_exp)
res <- lm(y_log~ Var_exp,data=mydatalog)
print(coef(res))
```
On voit que l'estimation des paramètres est mauvaise et la relation est mal modélisée.
```{r}
```{r f03-ggsim}
ggplot(mydata, aes(y=y.sim,x= Var_exp)) +
geom_function(fun = function(x) exp(1-0.2*x), col="green") +
geom_function(fun = function(x) exp(0.5122 - 0.0142*x)-1, col="red")+
geom_point(col="darkred") +
ggtitle(expression(paste(lambda(x), " et son estimation par transformation logarithme"))) +
geom_function(fun = function(x) exp(1-0.2*x), col="green") +
geom_function(fun = function(x) exp(0.5122 - 0.0142*x)-1, col="red")+
geom_point(col="darkred") +
ggtitle(expression(paste(lambda(x), " et son estimation par transformation logarithme"))) +
theme_classic()
```
Alors qu'avec la régression de Poisson que nous allons étudier dans ce chapitre:
```{r}
```{r f03-glm}
mydata <- data.frame(y=y.sim, Var_exp=Var_exp)
res <- glm(y~ Var_exp,family="poisson",data=mydata)
print(coef(res))
```
```{r}
```{r f03-glm2}
ggplot(mydata, aes(y=y.sim,x= Var_exp)) +
geom_function(fun = function(x) exp(1-0.2*x), col="green") +
geom_function(fun = function(x) exp(1.107 - 0.215*x), col="red")+
......@@ -83,7 +91,7 @@ ggplot(mydata, aes(y=y.sim,x= Var_exp)) +
2. Relation linéaire croissante où $Y \sim P(\lambda(x))$ et $\lambda(x)=e^{0.1+ 0.1 x}$
```{r,fig.height=3}
```{r f03-simu3,fig.height=3}
Var_exp <- seq(0,50,by=1)
lambda <- exp(0.1+0.1*Var_exp)
n=11
......@@ -95,8 +103,10 @@ ggplot(data.frame(y=y.sim,x=Var_exp),aes(x,y)) +
geom_point(col="darkred")+
theme_classic()
```
L'ajustement avec $glm$ donne aussi de très bon résultat sur ces données simulées:
```{r}
```{r f03-glm3}
mydata <- data.frame(y=y.sim, Var_exp=Var_exp)
res <- glm(y~ Var_exp,family="poisson",data=mydata)
print(coef(res))
......@@ -120,7 +130,7 @@ _La loi de la variable d'intérêt n'est plus une loi Normale mais une loi de Po
_Exemple: importation des données_
- Importation des données (fonction Import Dataset)
```{r}
```{r f03-import}
Weeds <- read.csv("data/weeds_simp.csv", sep=";",stringsAsFactors = TRUE)
```
......@@ -130,7 +140,7 @@ Weeds <- read.csv("data/weeds_simp.csv", sep=";",stringsAsFactors = TRUE)
+ id_field: identifiant champ
+ zone zone: d'échantillonga dans la champ (2 zone par champ)
```{r}
```{r f03-dm}
don <- Weeds %>%
mutate(id = paste(session,parcelle,systeme,transect,quadrat)) %>%
column_to_rownames("id")
......@@ -157,19 +167,19 @@ _Exemples_
- Abondance _(ab_tot)_
```{r}
ggplot(data=don,aes(x=don$ab_tot)) +
geom_histogram(bins=20,col="magenta",
fill="blue")
```{r f03-ggdon}
ggplot(data=don,aes(x=ab_tot)) +
geom_histogram(bins=20,col="magenta",fill="blue")
```
- Diversité _(rich_sp)_
```{r}
gg <- ggplot(don, aes(x=don$rich_sp))
```{r f03-ggdon2}
gg <- ggplot(don, aes(x=rich_sp))
gg <- gg + geom_histogram(binwidth=1, colour="black",
aes(y=..density.., fill=..count..))
gg <- gg + scale_fill_gradient("Count", low="#DCDCDC", high="#7C7C7C") + labs(x="Diversité",title="Distribution variable d'intérêt")
gg <- gg + scale_fill_gradient("Count", low="#DCDCDC", high="#7C7C7C") +
labs(x="Diversité",title="Distribution variable d'intérêt")
gg
```
......@@ -193,7 +203,7 @@ _Exemple : Mise en forme des variables explicatives_
+ _session_: on définit deux modalités s1 et s2
Comme cette variable n'a pas été importée avec le bon type, il faut la transformer en variable qualitative:
```{r}
```{r f03-dm2}
don$session <- factor(don$session,levels=1:2,labels=c("s1","s2"))
```
......@@ -202,12 +212,12 @@ don$session <- factor(don$session,levels=1:2,labels=c("s1","s2"))
Parmi les variables quantitatives du jeu de données, nous allons brièvement et visuellement étudier leurs relations pour identifier au moins une variable prédictive.
```{r}
```{r f03-summary}
summary(don[,c(5,7,8)])
```
```{r}
```{r f03-plot}
plot(don[,c(5,7,8)])
```
......@@ -244,7 +254,7 @@ Le modèle de régression de Poisson suppose que la valeur moyenne et la varianc
Affichons un résumé statistique par système:
```{r}
```{r f03-resum}
don %>% group_by(systeme) %>%
summarise(avg=mean(rich_sp,na.rm=TRUE),
variance = var(rich_sp,na.rm=TRUE))
......@@ -254,7 +264,7 @@ On observe un effet du système avec la _AF_ qui présente une densité plus él
Pour mieux visualiser le lien entre la distribution de la diversité et le système, réalisons un histogramme conditionnel.
```{r, fig.cap="Histogramme de la diversité par système"}
```{r f03-hist, fig.cap="Histogramme de la diversité par système"}
ggplot(don, aes(rich_sp, fill = systeme)) +
geom_histogram(binwidth=5, position="dodge")
```
......@@ -325,7 +335,7 @@ Les outils habituels d'inférence (intervalle de confiance et tests) restent val
Maintenant, nous allons réaliser une regression de Poisson sous _R_ à l'aide de la fonction _glm_.
```{r}
```{r f03-glm4}
m1 <- glm(rich_sp ~ systeme * dist_LSA, family="poisson", data=don)
summary(m1)
```
......@@ -333,13 +343,14 @@ summary(m1)
Cameron and Trivedi (2009) recommande l'utilisation d'erreurs types robustes des estimateurs lors de l'estimation des paramètres dans le cas où la condition d'égalité entre moyennes et variances conditionnelles ne seraient pas tout à fait réalisée. L'utilisation du package _sandwich_ de _R_ permet le calcul des erreurs types robustes et ensuite on peut en déduire les p-valeurs et intervalles de confiance à 95% associés.
```{r}
```{r f03-conf}
cov.m1 <- vcovHC(m1, type="HC0")
std.err <- sqrt(diag(cov.m1))
r.est <- cbind(Estimate= coef(m1), "Robust SE" = std.err,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err,
UL = coef(m1) + 1.96 * std.err)
r.est <- cbind(Estimate= coef(m1),
"Robust SE" = std.err,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err,
UL = coef(m1) + 1.96 * std.err)
r.est
```
......@@ -350,12 +361,12 @@ Etudions la sortie de la fonction _glm_ plus attentivement.
- On peut vouloir exprimer les paramètres (coefficients) dans leur unité d'origine et non plus en échelle log. Pour cela on utilise la Delta métode qui est implémenté dans la fonction _deltamethod_ du package _msm_. Principe : pour effectuer la transformation inverse rappelons que l'espérance d'un log n'est pas égale au log de l'espérance, on utilise alors un développement limité à l'ordre 1 pour approximer la loi du paramètre estimé détransformé, notamment sa moyenne et sa variance.
```{r}
```{r f03-deltam}
s <- deltamethod(list(~ exp(x1), ~ exp(x2), ~ exp(x3), ~exp(x4)), coef(m1), cov.m1)
```
```{r}
```{r f03-deltam2}
## exponentiate old estimates dropping the p values
rexp.est <- exp(r.est[, -3])
## replace SEs with estimates for exponentiated coefficients
......@@ -369,18 +380,18 @@ Les résultats précédents indiquent que la diversité moyenne pour le système
Pour comparer les valeurs moyennes de diversité attendues pour la distance moyenne observée, on peut utiliser la fonction _predict_ comme suit:
```{r}
```{r f03-resu4}
(s1 <- data.frame(dist_LSA = mean(don$dist_LSA),
systeme = factor(1:2, levels = 1:2, labels = levels(don$systeme))))
```
```{r}
```{r f03-predict4}
predict(m1, s1, type="response", se.fit=TRUE)
```
On retrouve les valeurs calculées avec la delta méthode, pour _systemeTA_ = 6.872467/5.881894 = `r 6.872467/5.881894`.
```{r, fig.cap="Prédiction de la diversité moyenne par système en fonction de la distance à l'arbre"}
```{r f03-gg4, fig.cap="Prédiction de la diversité moyenne par système en fonction de la distance à l'arbre"}
## calculate and store predicted values
don$rich_hat <- predict(m1, type="response")
......@@ -416,37 +427,37 @@ Si la déviance est élevée, cela signifie que le modèle utilisé ne décrit p
### Comparaison de modèles emboités
On va comparer deux modèles emboîtés $H_0 : \mathcal{M}_0$ et $H_1 : \mathcal{M}_1$, ayant mêmes distributions et mêmes fonctions de lien. Pour cela, on va comparer les déviances des deux modèles:
\begin{eqnarray*}
\begin{eqnarray}
\Delta D &=& D_0-D_1=2[\mathcal{L}(\widehat{\beta}_{saturé};y)-\mathcal{L}(\widehat{\beta}_{\mathcal{M}_0};y)]-2[\mathcal{L}(\widehat{\beta}_{saturé};y)-\mathcal{L}(\widehat{\beta}_{\mathcal{M}_1};y)]\\
&=& 2[\mathcal{L}(\widehat{\beta}_{\mathcal{M}_1}\mathcal{M}_;y)-\mathcal{L}(\widehat{\beta}_{\mathcal{M}_0};y)]
\end{eqnarray*}
\end{eqnarray}
- Si $\mathcal{M}_0$ et $\mathcal{M}_1$ décrivent bien les données, alors $D_0 \sim \chi^2_{n-q}$, $D_1 \sim \chi^2_{n-p}$, et $\Delta D \sim \chi^2_{p-q}$. On doit donc avoir pour $\Delta D$ une valeur plausible pour une $\chi^2_{p-q}$. Si c'est le cas, par principe de parcimonie on choisira $\mathcal{M}_0$.
- Si $\mathcal{M}_0$ décrit mal les données et $\mathcal{M}_1$ les décrit bien, $D_0$ sera plus grande qu'attendu pour une $\chi^2_{n-q}$, et $\Delta D$ sera plus grand qu'attendu pour une $\chi^2_{p-q}$.
- De manière générale, on effectuera le test de la façon suivante: si la valeur de $\Delta D$ est dans la région critique $[\chi^2_{p-q;1-\alpha};+\infty[$, on rejette $\mathcal{M}_0$ en faveur de $\mathcal{M}_1$. En effet, on considèrera que $\mathcal{M}_1$ décrit mieux les données que $\mathcal{M}_0$, de manière significative (même s'il ne les décrit pas parfaitement).
_Exemple: Déviance_
```{r}
```{r f03-dev}
summary(m1)
```
La déviance résiduelle est aussi calculée dans la sortie _R_. Elle peut être utilisée pour effectuer le test d'adéquation du modèle. La déviance résiduelle fournie par _R_ est la différence entre la déviance du modèle utilisé et la déviance du modèle saturé ($D_1$). En la soustrayant à la déviance nulle ($D_0$ où $\mathcal{M}_0$ ne contient qu'un intercept et aucune variable explicative), on peut construire un test statistique $\Delta D$. Si le test statistique associée à la déviance résiduelle est significatif, alors les différences résiduelles sont importantes et le modèle n'ajuste pas correctement les données. Il faudra ajouter d'autres variables explicatives dans le modèle ou tenir compte de la surdispersion des données.
On peut utiliser un test séquentiel qui va tester la significativité de l'ajout d'un terme après l'autre:
```{r}
```{r f03-anova}
## test model differences with chi square test
anova(m1, test="Chisq")
```
On peut aussi comparer directement deux modèles et par exemple étudier l'effet global du systeme et son interaction avec la distance à l'arbre en comparant le modèle avec et le modèle sans cette variable d'interaction (cas de modèles emboîtés).
```{r}
```{r f03-update}
## update m1 model dropping prog
m2 <- update(m1, . ~ . - systeme:dist_LSA)
## test model differences with chi square test
......@@ -455,7 +466,6 @@ anova(m2, m1, test="Chisq")
Quelle que soit l'approche (séquentielle croissante ou décroissante), tous les effets et facteurs sont significatifs.
### Comparaison de modèles non emboités
Pour des modèles non emboités ou avec des distributions ou des fonctions de liens différentes, on peut utiliser les critères AIC et BIC.
......@@ -472,7 +482,7 @@ On préfère les modèles ayant les plus petites valeurs pour ces critères.
On peut utiliser les résidus habituels (standardisés, studentisés), ceux de Pearson etc., comme en régression linéaire.
On définit aussi des résidus issu du calcul de la déviance. Cette dernière peut en effet être vue comme la somme de petites quantités, chacune amenée par une observation: $\sum_{i=1}^n d_i=D$. On définit alors le résidu associé à la $i^{\textrm{ème}}$ observation:
$$ e_{Di} = sign(y_i-\widehat{\lambda_i})\sqrt{d_i} $$
$$e_{Di} = sign(y_i-\widehat{\lambda_i})\sqrt{d_i} $$
On a $\sum_{i=1}^n e_{Di}^2=D$ la déviance. Ce sont les résidus utilisés par défaut par R. McCullagh et Nelder recommandent leur utilisation, car leur distribution se rapproche plus de la normalité.
......@@ -480,7 +490,7 @@ _Exemple: analyse des résidus_
L'information sur les résidus de la déviance dans le _summary_ montrent une légère asymétrie (valeur médiane non nulle et différence entre le 1er et le 3eme quartile).
```{r}
```{r f03-resid}
m1Diag <- data.frame(don,
link=predict(m1, type="link"),
fit=predict(m1, type="response"),
......@@ -493,7 +503,7 @@ m1Diag <- data.frame(don,
- Analyse du biais
```{r, eval=FALSE}
```{r f03-biais, eval=FALSE}
ggplot(data=m1Diag, aes(x=fit, y=resid)) +
geom_point() +
geom_abline(intercept = 0, slope = 0,col="red") +
......@@ -502,7 +512,7 @@ ggplot(data=m1Diag, aes(x=fit, y=resid)) +
Les résidus de la déviance de même que ceux de Pearson ne conviennent pas pour les cas de faibles effectifs. Ils sont visuellement non homogènes en général. Dans ce cadre, des auteurs proposent de simuler les quantiles des résidus (Dunn and Smyth, 1996). La méthode est implémenté dans le package _DHARMa_.
```{r, eval=TRUE}
```{r f03-dharma, eval=TRUE}
library(DHARMa)
res = simulateResiduals(m1)
plot(res)
......@@ -512,7 +522,7 @@ plot(res)
Pour cela, on essaie de représenter le lien entre la variance et la moyenne à travers la quantité $(y_i - \hat{y}_i)^2$ en fonction de la moyenne $\hat{y}_i$. On ajoute une droite verte sur ce graphique pour représenter la relation attendue dans le cas d'une loi de Poisson. On ajoute aussi une moyenne mobile en bleu des carrés résiduels ($(y_i - \hat{y}_i)^2$).
```{r}
```{r f03-diag}
ggplot(data=m1Diag, aes(x=fit, y=residSqr)) +
geom_point() +
geom_abline(intercept = 0, slope = 1,color="green") +
......@@ -531,7 +541,7 @@ Deux approches sont alors possibles:
- Tenir compte de la surdispersion en estimant le paramètre $\phi$ et en corrigeant les tests et le calcul des résidus en conséquence
```{r, eval=TRUE}
```{r f03-surdisp, eval=TRUE}
m2 <- glm(rich_sp~systeme + dist_LSA,
family="quasipoisson",
data=don)
......@@ -545,10 +555,9 @@ summary(m2)
- Quand il semble y avoir sur-dispersion, il faut s'assurer que le modèle est correctement spécifié, et notamment que des variables explicatives n'ont pas été oubliées. En effet, un modèle mal spécifié peut donner des résultats similaires à un problème de surdispersion.
- Sous la condition que le modèle soit bien spécifié, il faut tester si la condition d'égalité de la moyenne et de la variance est vérifiée. En cas de surdispersion, on peut ajuster un modèle de loi binomiale négative au lieu de Poisson. Pour tester la surdispersion, on peut utiliser le package _AER_ ou le package _pscl_ (Political Science Computational Laboratory, Stanford University) qui contiennent respectivement un dispersiontest et un odTest.
```{r, eval=TRUE}
```{r f03-aer, eval=TRUE}
library(AER)
deviance(m1)/m1$df.residual
dispersiontest(m1)
......@@ -557,15 +566,15 @@ m3 <- MASS::glm.nb(rich_sp~systeme + dist_LSA, data=don)
odTest(m3)
```
- Si la surdispersion est due à un excès de zéros, il faut utiliser un modèle adapté (zéro inflated modèle).
- Si la surdispersion est due à un excès de zéros, il faut utiliser un modèle adapté ("zéro inflated model").
```{r, eval=FALSE}
```{r f03-pscl, eval=FALSE}
library(pscl)
m4 <- zeroinfl(var.int~ var.exp, data=don, dist="poisson")
AIC(m1, m4)
```
- Si les données de comptage sont associées à une variable dite d'exposition, qui précise le nombre de fois où la donnée est exposé, il faut en tenir compte en l'ajoutant comme offset au modèle.
- Si les données de comptage sont associées à une variable dite d'exposition, qui précise le nombre de fois où la donnée est exposée, il faut en tenir compte en l'ajoutant comme offset au modèle.
- Deux règles à retenir sur le modèle de régression de Poisson : la variable d'intérêt ne peut pas prendre de valeurs négatives. La variable d'exposition ne peut prendre de valeurs nulles (pour l'offset).
......
# Régression non linéaire
```{r, include=TRUE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE,
message=FALSE, cache=FALSE,eval=TRUE)
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(nlme))
```
__Définition__
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 :
......
# AgroDataScience
Science des données pour l’agronomie et l’agroécologie : notes de cours
Livre bookdown de Data Science appliquée à l'agronomie sens large.
http://benedicte.fontez.pages.mia.inra.fr/agrodatascience/
Licence GPL-3 - l'institut Agro - Montpellier SupAgro - 2021
UMR MISTEA & XXX
UMR MISTEA & UMR ABsys
* Meïli Baragatti [aut]
* Bénédicte Fontez [cre, aut]
......@@ -14,7 +16,10 @@ UMR MISTEA & XXX
Livré créé avec le package **bookdown** dans R: https://bookdown.org/yihui/bookdown/ et compilé dans la forgemia Gitlab: https://about.gitlab.com/
IS 24/09/2021:
Temp.R contient le code du fichier 00- à corriger des erreurs latex de syntaxe: je fais ça progressivement.
IS 05/10/2021:
* correction de syntaxe dans 00: bcp de problèmes latex principalement
* suppression de l'appel de pcR dans 02-: l'image docker n'arrive pas à l'installer... j'ai mis les chunk en eval=FALSE afin que le code aparaisse mais qu'il ne soit pas executé
* rajout de la figure SchemaTest.png
* suppression d'une reference à une table kable qui ne fonctionnait pas
revoir l'appel de la figure SchemaTest.png dans 02-
book_filename: "agrodatascience"
book_filename: "testbd"
delete_merged_file: true
language:
ui:
chapter_name: "Chapter "
rmd_files: [
"index.Rmd",
"01-intro.Rmd",
"00-modele_lineaire.Rmd",
"01-intro.Rmd",
"03-GLM-loi.Rmd",
"02-GLM-variance-reduit.Rmd",
"04-NLR.Rmd",
"06-summary.Rmd",
"07-references.Rmd",
......
---
title: "PVD, promotion 2019-2020: Notes de cours "
title: "Science des données pour l’agronomie et l’agroécologie : notes de cours"
author: "Bénédicte Fontez, Meïli Baragatti & Léo Garcia - Montpellier SupAgro"
date: "`r Sys.Date()`"
......@@ -10,7 +10,7 @@ always_allow_html: true
url: http\://benedicte.fontez.pages.mia.inra.fr/agrodatascience
description: "Ce document comprend la compilation des notes de cours pour la régression linéaire généralisée et la régression non linéaire. Il a été réalisé grace à Rmarkdown et le package bookdown::gitbook. Les données et interprétations ont été fournies par Hélène Marroux et Léo Garcia."
description: "Ce document comprend la compilation des notes de cours pour la régression linéaire généralisée et la régression non linéaire. Il a été réalisé grace à **Rmarkdown** et le package **bookdown::gitbook**. Les données et interprétations ont été fournies par Hélène Marroux et Léo Garcia."
site: bookdown::bookdown_site
documentclass: book
......
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