Rendu du devoir sur le GLM
BRUSA PASQUE Cecile & CABANNE Faustine
Janvier 2021
1 Librairies utilisées et importation des données
On importe les librairies dont on a besoin :
library(dplyr) # Pour la manipulation de données##
## Attachement du package : 'dplyr'## Les objets suivants sont masqués depuis 'package:stats':
##
## filter, lag## Les objets suivants sont masqués depuis 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2) # Pour les graphiques
library(emmeans) # Pour les comparaisons de moyennesOn importe les bases de données qui seront utilisées dans ce TP :
#Exercice 1
presence_caesia <- read.table("presence_caesia.txt",
header = TRUE, # 1ere ligne donne le nom des colonnes
sep = ";", # Separateur de champ
encoding = "UTF-8") # Encodage du fichier initial
#Exercice 2
especes <- read.table("species.txt", sep = " ", header = TRUE)2 Présence de dehaasia caesia
2.1 Inclusion de la seule variable eau
2.1.1 Approche exploratoire
ggplot(presence_caesia) + # Données à représenter
aes(x = Presence, y = Eau) + # Ce qu'on représente (Présence en abscisse, Eau en ordonnee)
geom_boxplot() # Comment on le représente (boxplot)
Dâaprès la lecture de ce graphique, malgré lâobservation dâune tendance (moins de présence dans les milieux plus humides), nous ne pouvons pas conclure quant à la significativité de lâinfluence de lâeau sur la présence ou non de dehaasia caesia car les boxplot se chevauchent.
Nous avons également fait un graphique cherchant à voir lâinfluence des différents types de sols :
ggplot(presence_caesia) + # Données à représenter
aes(x = Sol, y = Eau, fill=Presence) +
geom_boxplot() # Comment on le représente (boxplot)
Il nây a toujours pas de tendance significative observable entre la présence ou lâabsence de la plante selon la quantité dâeau contenue dans le sol. Il peut être intéressant cependant de continuer à explorer lâhypothèse de lâinfluence de lâeau sur la présence de Dehaasia caesia, avec un test pour mesurer la significativité de lâimpact de cette variable dâaprès un modèle de régression logistique.
2.1.2 Modèle de régression logistique
On suppose que \(y_k\) est la réalisation dâune variable aléatoire \(Y_k\) telle que: \[Y_k\sim Bernoulli(p_k),\] où \[\mathrm{logit}(p_k)=\log\left(\frac{p_k}{1-p_k}\right)=\beta_0+\beta_1 x_k,\] et où les paramètres \(\beta_0\) et \(\beta_1\) sont des inconnues à estimer.
2.1.3 Ajustement avec R
# La fonction mutate permet de créer ou modifier une colonne
# Ce code est strictement equivalent Ã
# presence_caesia$Y <- ifelse(presence_caesia$Presence == "oui", 1, 0)
presence_caesia <- mutate(presence_caesia,
# La fonction ifelse permet de binariser
Y = ifelse(Presence == "oui", # Condition de test
1, # Si TRUE, on renvoie 1
0) # Si FALSE, on renvoie 0
)
#On ajuste le modèle :
modele_eau <- glm(Y ~ Eau, # Modèle (avec la variable Y créée plus haut)
data = presence_caesia,
family = binomial(link = "logit") # Famille de GLM avec
# fonction de lien
)
summary(modele_eau)##
## Call:
## glm(formula = Y ~ Eau, family = binomial(link = "logit"), data = presence_caesia)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6115 -1.2758 0.8899 1.0218 1.5643
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1940 0.4253 2.808 0.00499 **
## Eau -0.2551 0.1223 -2.086 0.03702 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 242.03 on 178 degrees of freedom
## Residual deviance: 237.52 on 177 degrees of freedom
## AIC: 241.52
##
## Number of Fisher Scoring iterations: 4Sur la ligne
Eau, on lit le résultat du test dont les deux hypothèses sont les suivantes :
\(H_0 : \mathrm{logit}(p_k)=\beta_0\) vs \(H_1 : \mathrm{logit}(p_k)=\beta_0+\beta_1 x_k\)
ou plus simplement :
\(H_0 : \beta_1=0\) vs \(H_1 : \beta_1\neq 0\)On utilise la statistique de test du Z-score : \[Z=\frac{\hat\beta_1}{\sqrt{\mathbb{V}(\hat\beta_1)n}}\]
On a vu que sous \(H_0\), \(Z\xrightarrow[n\to+\infty]{\mathcal{L}}\mathcal{N}(0,1)\) Ce nâest pas la loi exacte, câest une loi asymptotique. Le nombre dâobservation est néanmoins assez élevé ici \(n = 146\) pour que cette approximation soit bonne.
Au niveau \(5\%\), on conclue donc à la significativité de lâeffet de la teneur en eau sur la présence de Dehaasia caesia, la P-value étant égale à \(0.03702\) \((<0,05)\).
2.2 Inclusion du type de sol
2.2.1 Ecriture du modèle
On dispose de \(n=179\) sites répartis sur trois types de sols quâon codera de \(1\) à \(3\) (\(1\) pour Alluvial, \(2\) pour Dunaire, et \(3\) pour Grès). Sur chaque sol, on dispose respectivement de \(n_1=60\), \(n_2=60\) et \(n_3=59\) observations.
Pour tout \(1\leq i\leq 3\), \(1\leq k\leq n_i\), on note \(y_{i,k}\) la présence ou lâabsence de lâespèce sur le \(k\)ième site de type de sol \(i\), avec la convention \(y_{i,k}=1\) si lâespèce est présente et \(y_{i,k}=0\) sinon. De plus, on note \(x_{i,k}\) la teneur en eau sur ce site.
Pour tout \(1\leq i\leq 3\), \(1\leq k\leq n_i\), on suppose que \(y_{i,k}\) est la réalisation dâune variable aléatoire \(Y_{i,k}\) telle que : \[Y_{i,k}\sim Bernouilli(p_{i,k}),\] où \[\mathrm{logit}(p_{i,k})=\beta_0+\alpha_i+(\beta_1+\gamma_i)x_{i,k},\]
et où les paramètres \(\beta_0, \beta_1, \alpha_1, \alpha_2, \alpha_3, \gamma_1, \gamma_2, \gamma_3\) sont des inconnues à estimer.
2.2.2 Ajustement avec R
modele_logistique_complet <- glm("Y ~ Sol + Eau + Sol:Eau",
data = presence_caesia,
# On précise la famille
family = binomial(link = "logit"))On effectue le test de type II des effets :
car::Anova(modele_logistique_complet)## Analysis of Deviance Table (Type II tests)
##
## Response: Y
## LR Chisq Df Pr(>Chisq)
## Sol 24.7553 2 4.212e-06 ***
## Eau 0.0896 1 0.76474
## Sol:Eau 5.4940 2 0.06412 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sur la ligne
Sol:Eau, on teste :- \(H_0 : \mathrm{logit}(p_{i,k})=\beta_0+\alpha_i+\beta_1 x_{i,k}~(M_0)\)
- \(H_1 : \mathrm{logit}(p_{i,k})=\beta_0+\alpha_i+(\beta_1+\gamma_i) x_{i,k}~(M_1)\)
Il sâagit ici dâun test de modèles emboîtés, la statistique de test est la statistique du rapport de vraisemblance \(LRT\) :
\[LRT(M_1,M_0)=2(\ell(Y,\hat\theta_1)-\ell(Y,\hat\theta_0)),\]
où \[\ell(Y,\hat\theta_1)=\sum_{i=1}^3\sum_{k=1}^{n_i}\left[Y_{i,k}x_{i,k}\hat\theta_1+\ln(1+e^{x_{i,k}\hat\theta_1})\right],\] avec \(x_{i,k}\hat\theta_1=\hat\beta_0+\hat\alpha_i+(\hat\beta_1+\hat\gamma_i)x_{i,k}\), et \[\ell(Y,\hat\theta_0)=\sum_{i=1}^3\sum_{k=1}^{n_i}\left[Y_{i,k}x_{i,k}\hat\theta_0+\ln(1+e^{x_{i,k}\hat\theta_0})\right],\] avec \(x_{i,k}\hat\theta_0=\hat\beta_0+\hat\alpha_i+\hat\beta_1 x_{i,k}\)On a vu que sous \(H_0\) (lorsque \(H_0\) est vérifiée), \[LRT(M_1,M_0)\xrightarrow[n\to+\infty]{\mathcal{L}}\chi^2_q,\] où \(q=(8-2)-(5-1)=2\) est le nombre de régresseurs présents dans \(M_1\) et pas dans \(M_0\)..
Au risque \(\alpha = 5\%\), on conclue que lâinformation apportée par lâinteraction
Sol:Eaunâest pas significative (P-value > 0.05). Cependant lâinformation apportée par la variableSolseule est très significative : sa P-value est très inférieure à 0.001. On remarque par ailleurs que dans ce modèle, il sâagit de la seule variable ayant un apport significatif dans lâexplication de la présence de Dehaasia caesia.Sur la ligne
Eau, on teste : \[H_0 : \mathrm{logit}(p_{i,k})=\beta_0+\alpha_i~(M_0)\] \[H_1 : \mathrm{logit}(p_{i,k})=\beta_0+\alpha_i+\beta_1 x_{i,k}~(M_1)\]
ou encore \[H_0 : \beta_1=0~(M_0)\] \[H_1 : \beta_1\neq0~(M_1)\]
On obtient un résultat en contradiction avec celui obtenu précédemment en 2.1.3. En effet, lâapport dâinformation par la variableEauseule dans le modèle est à présent considérée comme non-significative. Cela sâexplique par lâintroduction de la variableSol, qui rend la variableEauredondante (la teneur en eau dépendant déjà du type de sol).
summary(modele_logistique_complet)##
## Call:
## glm(formula = "Y ~ Sol + Eau + Sol:Eau", family = binomial(link = "logit"),
## data = presence_caesia)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0782 -1.0538 0.5395 1.0252 1.6634
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5614 1.1977 1.304 0.1923
## SolDunaire -2.2001 1.3734 -1.602 0.1092
## SolGrès 0.7256 1.6628 0.436 0.6625
## Eau -0.4965 0.2888 -1.719 0.0856 .
## SolDunaire:Eau 0.8319 0.3834 2.170 0.0300 *
## SolGrès:Eau 0.2994 0.4687 0.639 0.5229
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 242.03 on 178 degrees of freedom
## Residual deviance: 207.27 on 173 degrees of freedom
## AIC: 219.27
##
## Number of Fisher Scoring iterations: 4- On a :
- \(\hat\beta_0 = 1.5614\) ;
- \(\hat\beta_1 = -0.4965\) ;
- \(\hat\alpha_1 = 0\) ;
- \(\hat\alpha_2 = -2.2001\) ;
- \(\hat\alpha_3 = 0.7256\) ;
- \(\hat\gamma_1 = 0\) ;
- \(\hat\gamma_2 = 0.8319\) ;
- \(\hat\gamma_3 = 0.2994\).
- \(\hat\beta_0 = 1.5614\) ;
Si les résultats peuvent en premier lieu avoir lâair contradictoires (les P-values de SolDunaire et SolGrèssont supérieures à 0.05), il nâen est rien puisque nous réalisons ici un test de type II. Cela signifie que sur la ligne SolGrès, le modèle considère déjà la modalité précédente SolDunaire. La variable Sol peut être significative (comme vu précédemment) sans que lâapport des effets des modalités SolDunaire, puis SolDunaire et SolGrès au modèle, soit significatif.
modele_logistique_sol<-glm("Y~Sol",
data=presence_caesia,
family=binomial(link="logit"))
#summary(modele_logistique_sol)
#car::Anova(modele_logistique_sol)
#Affichage des critères d'AIC des modèles
AIC(modele_logistique_complet) #modèle contenant les effets du Sol, de l'Eau et des interactions## [1] 219.2708AIC(modele_logistique_sol) #modèle contenant seulement les effets du Sol## [1] 218.8544On voit que le critère dâAkaike est meilleur (le minimum entre les 2 valeurs obtenues) pour le modele_logistique_sol. Il sâagit donc du meilleur modèle logistique.
Suite à ces analyses, nous pouvons conclure que la présence de Dehaasia caesia est corrélée à la teneur en eau du sol, car cette teneur en eau dépend du type de sol. Cela se vérifie par la significativité du paramètre Sol, qui, lorsquâon lâintroduit dans le modèle, entraîne une perte de significativité pour le paramètre Eau (cf. summary du modèle logistique complet).
2.3 Prédictions
probabilites_predites <- predict(modele_logistique_complet,
type = "response") # On prédit l'espérance de Y
presence_predite <- as.numeric(probabilites_predites > 0.5) # On predit 1
# si la proba predite est supérieure à 0.5
# table de contingence
table(presence_predite, presence_caesia$Y,
dnn = c("Prédiction", "Vérité"))## Vérité
## Prédiction 0 1
## 0 43 27
## 1 30 79- Le taux dâerreur global (présence prédite alors quâil y a absence et inversement/nombre total de prédictions) correspond à :
(30+27)/(43+27+30+79)*100## [1] 31.84358\(32\%\).
- Quant à la probabilité quâon se trompe si on prédit que lâespèce est présente, elle est de :
30/(30+79)*100## [1] 27.52294\(28\%\).
3 Biodiversité et biomasse
3.1 Approche descriptive
ggplot(especes) + # Données à représenter
aes(x = Biomass, y = Species,# Ce qu'on représente
color = pH, size = surface) + # La taille donne la surface couverte
geom_point() +
labs(x = "Biomasse", y = "Nombre d'espèces")
Le graphe montre que, plus la biomasse du sol est élevée, moins on observe dâespèces de plantes différentes, il semble donc y avoir un lien décroissant entre biomasse et nombre dâespèces. Parallèlement, le nombre total dâespèces semble être lié au type de sol : plus le pH est élevé, plus le nombre dâespèces végétales est élevé.
On distingue 3 groupes distincts qui semblent se répartir selon 3 droites parallèles. Ainsi, le lien entre biomasse et nombre dâespèces ne dépendrait pas du pH.
3.2 Modèle
On dispose de \(n=90\) observations effectuées sur des de sols de gammes de pH différentes et indicées de \(1\) à \(3\) (\(1\) pour basique, \(2\) pour neutre, et \(3\) pour acide). Pour chaque gamme de pH, on dispose de \(n_1=n_2=n_3=30\) observations.
Pour tout \(1\leq i\leq 3\), \(1\leq k\leq n_i\), on note \(y_{i,k}\) le nombre dâespèces de plantes répertoriées sur le \(k\)ième site dans la gamme de pH \(i\). De plus, on note \(x_{i,k}\) la biomasse mesurée sur ce site.
Pour tout \(1\leq i\leq 3\) et \(1\leq k\leq n_i\), on suppose que \(y_{i,k}\) est la réalisation dâune variable aléatoire \(Y_{i,k}\) telle que : \[Y_{i,k}\sim Poisson(\lambda_{i,k}\times w_{i,k}),\]
où
\[\mathrm{ln}(\lambda_{i,k})=\beta_0+\alpha_i+(\beta_1+\gamma_i)x_{i,k},\]
avec \(\lambda_{i,k}\) le paramètre de la loi de Poisson, \(w_{i,k}\) la surface sur laquelle lâétude a été réalisée, qui ici fait office de variable dâeffort dâéchantillonnage (offset) et où les paramètres \(\beta_0, \beta_1, \alpha_1, \alpha_2, \alpha_3, \gamma_1, \gamma_2, \gamma_3\) sont des inconnues à estimer.
3.3 Ajustement avec R
modele_poisson <- glm(Species ~ pH * Biomass , data = especes,
family = poisson(link = "log"),
offset = log(surface)) # En R, on met l'offset dans
# l'espace de la combinaison linéaire (donc en log, ici)On effectue les tests en type II :
car::Anova(modele_poisson)## Analysis of Deviance Table (Type II tests)
##
## Response: Species
## LR Chisq Df Pr(>Chisq)
## pH 344.09 2 < 2.2e-16 ***
## Biomass 176.28 1 < 2.2e-16 ***
## pH:Biomass 24.72 2 4.288e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sur la ligne
pH:Biomass, on teste :- \(H_0 : \mathrm{ln}(\lambda_{i,k})=\beta_0+\alpha_i+\beta_1 x_{i,k}~(M_0)\)
- \(H_1 : \mathrm{ln}(\lambda_{i,k})=\beta_0+\alpha_i+(\beta_1+\gamma_i) x_{i,k}~(M_1)\)
Il sâagit ici dâun test de modèles emboîtés, la statistique de test est la statistique du rapport de vraisemblance \(LRT\) :
\[LRT(M_1,M_0)=2(\ell(Y,\hat\theta_1)-\ell(Y,\hat\theta_0)),\]
où \[\ell(Y,\hat\theta_1)=\sum_{i=1}^3\sum_{k=1}^{n_i}\left[Y_{i,k}x_{i,k}\hat\theta_1+\ln(1+e^{x_{i,k}\hat\theta_1})\right],\] avec \(x_{i,k}\hat\theta_1=\hat\beta_0+\hat\alpha_i+(\hat\beta_1+\hat\gamma_i)x_{i,k}\), et \[\ell(Y,\hat\theta_0)=\sum_{i=1}^3\sum_{k=1}^{n_i}\left[Y_{i,k}x_{i,k}\hat\theta_0+\ln(1+e^{x_{i,k}\hat\theta_0})\right],\] avec \(x_{i,k}\hat\theta_0=\hat\beta_0+\hat\alpha_i+\hat\beta_1 x_{i,k}\)On a vu que sous \(H_0\) (lorsque \(H_0\) est vérifiée), \[LRT(M_1,M_0)\xrightarrow[n\to+\infty]{\mathcal{L}}\chi^2_q,\] où \(q=(8-2)-(5-1)=2\) est le nombre de régresseurs présents dans \(M_1\) et pas dans \(M_0\). Au risque \(\alpha = 5\%\), on conclue que lâinformation apportée par lâinteraction
pH:Biomassest très significative : sa P-value est très inférieure à 0.001. On a donc tout intérêt à conserver cette intéraction dans le modèle. De plus, tous les autres paramètres sont également significatifs, on les conserve aussi dans le modèle.
Tableau_valeurs <- summary(modele_poisson)
Tableau_valeurs##
## Call:
## glm(formula = Species ~ pH * Biomass, family = poisson(link = "log"),
## data = especes, offset = log(surface))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1584 -0.9719 0.0591 0.8544 3.2951
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.81462 0.06090 62.632 < 2e-16 ***
## pHlow -0.79766 0.10282 -7.758 8.62e-15 ***
## pHmid -0.27966 0.09306 -3.005 0.00266 **
## Biomass -0.10301 0.01230 -8.375 < 2e-16 ***
## pHlow:Biomass -0.18561 0.04031 -4.605 4.13e-06 ***
## pHmid:Biomass -0.05608 0.02362 -2.375 0.01757 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 559.61 on 89 degrees of freedom
## Residual deviance: 143.70 on 84 degrees of freedom
## AIC: 574.89
##
## Number of Fisher Scoring iterations: 4- On a :
- \(\hat\beta_0 = 3.81462\) ;
- \(\hat\beta_1 = -0.10301\) ;
- \(\hat\alpha_1 = 0\) ;
- \(\hat\alpha_2 = -0.27966\) ;
- \(\hat\alpha_3 = -0.79766\) ;
- \(\hat\gamma_1 = 0\) ;
- \(\hat\gamma_2 = -0.05608\) ;
- \(\hat\gamma_3 = -0.18561\).
- \(\hat\beta_0 = 3.81462\) ;
On a \(w_i=1\), ainsi lâintensité est simplement égale à \(\lambda_i\) avec \(\mathrm{ln}(\lambda_{i})=\hat\beta_0+\hat\alpha_i+(\hat\beta_1+\hat\gamma_i)x_{obs}\)
valeur_1 = 3.81462+0+(-0.10301+0)*2.5
valeur_1## [1] 3.557095exp(valeur_1)## [1] 35.0612Or, pour un sol au pH basique, on a : \(\hat\beta_0+\hat\alpha_1+(\hat\beta_1+\hat\gamma_1)x\approx3.56\) donc \(\lambda_1\approx35\)
valeur_2 = 3.81462-0.27966+(-0.10301-0.05608)*2.5
valeur_2## [1] 3.137235exp(valeur_2)## [1] 23.04007Pour un sol au pH neutre, on a : \(\hat\beta_0+\hat\alpha_2+(\hat\beta_1+\hat\gamma_2)x\approx3.14\) donc \(\lambda_2\approx23\)
valeur_3 = 3.81462-0.79766+(-0.10301-0.18561)*2.5
valeur_3## [1] 2.29541exp(valeur_3)## [1] 9.928506Pour un sol au pH acide, on a : \(\hat\beta_0+\hat\alpha_3+(\hat\beta_1+\hat\gamma_3)x\approx2.3\) donc \(\lambda_3\approx10\)
On réalise ensuite un test de comparaison deux à deux des log intensités selon les pH.
emmeans(modele_poisson, specs = "pH") %>%
contrast(method = "pairwise", adjust = "bonferroni")## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df z.ratio p.value
## high - low 1.458 0.0977 Inf 14.915 <.0001
## high - mid 0.479 0.0561 Inf 8.533 <.0001
## low - mid -0.979 0.1007 Inf -9.715 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: bonferroni method for 3 testsOn teste les hypothèses suivantes :
\(H_0 : \mathrm{ln}(\lambda_i)=\mathrm{ln}(\lambda_{i'})\) vs \(H_1 : \mathrm{ln}(\lambda_i)\neq \mathrm{ln}(\lambda_{i'})\)
où
\(1\leq i,i'\leq 3, i\neq i'\).
Câest-à -dire, a-t-on une différence significative du lien entre biomasse et nombre dâespèces végétales selon les différentes gammes de pH ?
Ici toutes les P-values sont inférieures à 0.001, donc au risque \(\alpha=5\%\), on peut dire que pour toutes les lignes on rejette lâhypothèse \(H_0\). Les différences sont significatives entre toutes les gammes de pH.
Pour conclure, le nombre dâespèces de plantes sur un site est très bien expliqué par la biomasse présente et le pH du sol de ce dernier. Par ailleurs, même si sur le graphique obtenu initialement cela ne semblait pas le cas, le lien entre le nombre dâespèces végétales et la biomasse sur un site varie selon le pH.