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, union
library(ggplot2) # Pour les graphiques
library(emmeans) # Pour les comparaisons de moyennes

On 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: 4

  • Sur 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 ' ' 1

  • Sur 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:Eau n’est pas significative (P-value > 0.05). Cependant l’information apportée par la variable Sol seule 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 variable Eau seule dans le modèle est à présent considérée comme non-significative. Cela s’explique par l’introduction de la variable Sol, qui rend la variable Eau redondante (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\).

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.2708
AIC(modele_logistique_sol) #modèle contenant seulement les effets du Sol
## [1] 218.8544

On 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 ' ' 1

  • Sur 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:Biomass est 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\).

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.557095
exp(valeur_1)
## [1] 35.0612

Or, 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.137235
exp(valeur_2)
## [1] 23.04007

Pour 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.29541
exp(valeur_3)
## [1] 9.928506

Pour 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 tests

On 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.