ENV2015 -Gestion de l’environnement, Hiver 2022, Valérie BÉRUBÉ
QUESTION 1
a) Hypothèses et choix d’analyse
Hypothèse biologique: Les terreaux à base de tourbe de sphaigne favoriseront les éricacées.
Hypothèses statistiques: H0: µsoil = µpeat = µmix =µmine et H1: au moins 1 µ diffère des autres (α = 0.05)
Afin de tester ces hypothèss, le dispositif en bloc et l’analyse de variance ANOVA devraient être utilisés. En effet, l’ANOVA permettra de comparer les moyennes de plus de deux groupes afin de déterminer si au moins une moyenne diffère. Étant donné que les dispositifs expérimentaux sont répartis sur plusieurs sites, une ANOVA avec blocs complets aléatoires permettra d’exercer un contrôle sur l’hétérogénéité du milieu naturel dans lequel les dispositifs expérimentaux se trouvent. Cette démarche permettra de diminuer le risque de variabilité entre les sites.
b) Structure des données brutes
#---Specifier le repertoire de travail
setwd("C:\\Users\\Pounie\\Documents\\Universite_Laval_21-22\\Gestion\\Berube_ex2")
#---Importer les données
mine<- read.table("mine.txt", header= TRUE)Sommaire des données brutes
#---Structure de l'objet "mine"
head(mine)## Site Substrate Biomass
## 1 site1 mine 0.62
## 2 site1 mix 2.03
## 3 site1 peat 0.83
## 4 site1 soil 0.55
## 5 site2 mine 2.39
## 6 site2 mix 3.14summary(mine)## Site Substrate Biomass
## Length:40 Length:40 Min. :0.0100
## Class :character Class :character 1st Qu.:0.8975
## Mode :character Mode :character Median :1.7450
## Mean :1.7432
## 3rd Qu.:2.4175
## Max. :4.5500Distribution des données à l’aide de graphiques
Deux graphiques ont été produits afin de présenter les données brutes. L’utilisation de boxplots permet d’obtenir un apperçu rapide de l’étendue des variances pour chacun des groupes. Le premier graphique met en relation les valeurs de biomasses relevées pour chacun des traitments (substrats).
#---Boxplot de la biomasse produite selon le substrat
boxplot(Biomass~Substrate, data=mine, ylab = "Biomasse (kg/m2)", xlab="Substrats", col="chartreuse4", main="Distribution des valeurs de biomasse selon les traitements")
L’analyse d’un tel graphique en vient à rechercher des défaillances quant à la symétrie. Quelques observations sont à souligner pour cette figure. D’abord, il est posssible de voir que le substrat “peat” génère une plus grande variablité de par l’étendue de l’écart interquartile. Ensuite, les quartiles de “mine” sont peu proportionnels et peuvent indiquer une différence entre la concentration des données à l’intérieur du groupe. Enfin, les médianes de “soil” et “mine” sont ni centrées dans le box plot, ni relativement au même niveau que les autres traitements. Ces caractéristiques peuvent induire que certains traitements fournissent des résultats différents. Enfin, il y a peu de valeurs aberrantes, mis à part pour “soil”, ce qui contribue à l’homogénéité des données.
Le second graphique met en relation la biomasse produite par chacun des sites. Il convient de présenter les données brutes selon le site afin de prendre en compte une possible variabilité induite par la présence d’échantillons répartis sur plusieurs sites.
#---Boxplot de la biomasse produite selon le site
boxplot(Biomass~Site, data=mine, ylab = "Biomasse (kg/m2)", xlab="Substrats", col="orange", main="Distribution des valeurs de biomasse selon les sites",las=2)
L’envergure irrégulière des box plots (ex: site 4 et 5),la différence dans l’étendue des données (moustaches) (ex: site 10 et site 4) et des médianes différentes entre chacun des groupe peut supposer une variabilité inter et intra-groupes. .
c) Analyse
#---Convertir la variable "Site" en facteur afin d'inclure l'effet du block
setwd("C:\\Users\\Pounie\\Documents\\Universite_Laval_21-22\\Gestion\\Berube_ex2")
mine.facteur<- read.table("mine.txt", header= TRUE)Vérification des conditions d’application
#---Réaliser l'ANOVA
mine.aov<- aov(Biomass~Substrate+Site,data=mine.facteur)
#---Organisation des graphiques
par(mfrow = c(2, 2))
#---Analyse graphique
plot(mine.aov)## hat values (leverages) are all = 0.325
## and there are no factor predictors; no plot no. 5
D’après les informations disponibles, il est possible de conclure que l’ANOVA respecte les conditions d’application. D’abord,l’échantillonnage est indépendant et aléatoire (confirmé par les consignes de l’exercice) . Ensuite,chaque variable échantillonnée semble être distribuée normalement dans sa population. Notamment, le graphique de normalité des résidus (QQ plot) illustre que la majeure partie des points sont distribués en suivant la droite de normalité. Aussi, il est possible d’y voir qu’il y a peu de valeurs extrêmes. Par la suite, la condition d’homoscédasticité est confirmée par le graphique Residuals vs Fitted qui illustre que les résidus sont répartis aléatoirement et de manière homogène autour de l’horizon 0.
d) Interprétation des résultats de l’analyse
#---Visualiser les résultats du test
summary(mine.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## Substrate 3 8.983 2.9942 5.427 0.004715 **
## Site 9 24.486 2.7207 4.931 0.000582 ***
## Residuals 27 14.896 0.5517
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pour la varibale Substrate, le résultat d’ANOVA indique une p-value de 0.004715, soit 0.004715 < α=0.05. Pour Site, la p-value est de 0.000582, soit 0.004715 < α=0.05. Ces résultats de p-value plus petites que le seuil de 0.05 indique qu’il est possible de rejeter H0. Ainsi,le test d’ANOVA indique qu’il existe au moins une différence significative entre les moyennes des groupes.
Analyse des différences-Comparaisons multiples
Le test de Tukey est utilisé sur la variable Substrat. Elle permet de vérifier quelles moyenness sont différentes et quelle est l’amplitude de la différence. L’effet du bloc (Site) ne nécessitera pas une comparaison multiple.
#---Test de Tukey
TukeyHSD(mine.aov, which = "Substrate", order = TRUE)## Tukey multiple comparisons of means
## 95% family-wise confidence level
## factor levels have been ordered
##
## Fit: aov(formula = Biomass ~ Substrate + Site, data = mine.facteur)
##
## $Substrate
## diff lwr upr p adj
## soil-mine 0.713 -0.196024 1.622024 0.1642967
## peat-mine 1.077 0.167976 1.986024 0.0156540
## mix-mine 1.227 0.317976 2.136024 0.0051488
## peat-soil 0.364 -0.545024 1.273024 0.6950602
## mix-soil 0.514 -0.395024 1.423024 0.4245062
## mix-peat 0.150 -0.759024 1.059024 0.9687533D’après les résultats, “peat” et “mix” ont une différence significative en raison de leur p-value < 0.05. Ces résultats justifient le rejet de l’hypothèse nulle, voulant que les moyennes soient égales dans tous les groupes. L’hypothèse biologique est d’autant plus infirmée puisque les résultats indiquent que ce sont les terreaux à base de tourbe de sphaigne et de terre arable “mix” qui favorisent les éricacées.
Un graphique de Tukey confirme cette conclusions en illustrant la relation entre chacune des variables.
#---Disposition des marges
par (mar= c(6,6,6,6))
#---Visualisation graphique du test de Tukey
Tukey<-TukeyHSD(mine.aov, which = "Substrate", order = TRUE)
plot(Tukey,las=1)
Le constat est que “peat-mine” et “mix-mine” ne comprennent pas “0” dans leur intervalle de confiance. Cela confirme qu’ils sont significativement différents.
e) Résultats de l’analyse d’ANOVA
Le graphique suivant présente les résulats de l’analyse d’ANOVA dans un graphique avec des barres d’erreurs.
#---Moyenne par traitement
moy.groupe<-tapply(mine.facteur$Biomass,INDEX=mine.facteur$Substrate, FUN=mean)
#---Barres d'erreurs
groupe.res<-sqrt(0.5517/10)
#---Calcul de l'intervalle de confiance
low95<-moy.groupe+qt(p=0.025,df=mine.aov$df.residual)* groupe.res
up95<-moy.groupe+qt(p=0.975,df=mine.aov$df.residual)* groupe.res
#---Assignation des valeurs de x
val.x<-1:4
#---Graphique avec intervalles de confiance
plot(moy.groupe ~ val.x, type = "p", ylab = "Biomasse (kg/m2)", xlab = "Substrats", cex.lab = 1.2, cex.axis = 1.2, ylim = range(c(low95, up95)), xaxt = "n", main= "Moyennes± IC 95%")
axis(side = 1, at = val.x,
labels = names(moy.groupe))
segments(x0 = val.x,
y0 = low95,
x1 = val.x,
y1 = up95)
Étant donné que la barre d’erreur de mine ne se chevauche pas avec celle de mix ni peat, les groupes mix et peat seraient significativement différents, sachant que mine est une variable de contrôle.
QUESTION 2
a) Représentation graphique des données brutes
#---Specifier le repertoire de travail
setwd("C:\\Users\\Pounie\\Documents\\Universite_Laval_21-22\\Gestion\\Berube_ex2")
#---Importer les données
richesse<- read.table("richesse.txt", header= TRUE)
#---Structure de l'objet "richesse"
head(richesse)## Site Time.abandon Mean.richness
## 1 site1 43 100.2
## 2 site2 24 53.4
## 3 site3 12 8.0
## 4 site4 32 65.6
## 5 site5 11 30.7
## 6 site6 2 11.4#---Graphique des données brutes
plot(Mean.richness~Time.abandon, data=richesse, xlab="Temps depuis l'abandon du site(années)",ylab="Nbre moyen d'espèces par point d'écoute",cex.axis=1.2,cex.lab=1.2, main="Nombre moyen d'espèces par point d'écoute selon le temps depuis l'abandon du site", cex.main=.75)
La disposition des données révèle une relation linéaire entre celles-ci puisque les points se concentrent autour d’une même ligne directrice imaginaire. Il ne semble pas y avoir de données extrêmes.
Hypothèses
Les hypothèses statistiques sont les suivantes:
HO: βTime.abandon = 0. La pente est nulle. Il n’y a pas de relation linéaire entre la variable X et Y.
H1: βTime.abandon ≠ 0. Il y a une relation linéaire entre la variable X et Y. (α = 0.05)
b) Régression linéaire
Le graphique suivant présente la régression linéaire produite par une mise en relation des deux variables.
#--- Graphique
regression<-plot(Mean.richness~Time.abandon, data=richesse, xlab="Temps depuis l'abandon du site(années)", ylab="Nombre moyen d'espèces", main="Nombre moyen d'espèces par point d'écoute selon le temps depuis l'abandon du site", cex.main=.75)
#--- Droite de régression
droite.regression<- lm(Mean.richness~Time.abandon, data=richesse)
abline(droite.regression)
Vérification des conditions d’application
setwd("C:\\Users\\Pounie\\Documents\\Universite_Laval_21-22\\Gestion\\Berube_ex2")
#--- Disposition des graphiques
par(mfrow = c(2, 2))
#---Vérification graphique des conditions d'application
plot(droite.regression)
Voici le constat quant aux conditions d’application de la régression:
La condition d’existence est validée par la présence d’une valeur de Y pour chaque valeur de X (possibilité d’en faire le constat en observant les données brutes).
La normalité des données peut se vérifier par le graphique Normal Q-Q. Le graphique représente la normalité de la distribution des résidus. Afin que la condition de distribution normale des résidus soit respectée, les points devraient se trouvés alignés le long de la droite. L’analyse des résidus permet de constater qu’une donnée est plutôt extrême et que, de manière générale, les données ne sont pas toutes très alignées le long de la droite. Cela laisse supposer que la condition d’application pourrait ne pas être respectée. Cependant, aucune conclusion n’est encore tirée et l’analyse des données avec la distance de Cook permettra de confirmer si les données influencent significativement la droite de régression.
L’homoscédasticité des résidus veut que les résidus soient normalement distribués autour des valeurs prédites par la régression. Cela peut être confirmé par les graphiques ci-haut. Le graphique Residual vs Fitted présente des données répartis aléatoirement de façon homogène autour de l’horizon 0, ce qui suppose que cette condition est respectée.
La condition d’indépendance exige l’absence de liens entre les Yi. cela est confirmé du fait qu’il y a trois points d’écoutes qui sont distincts pour chacun des sites.
La condition de linéarité exige une relation linéaire entre Y et X (condition confirmée par l’observation du graphique de régression)
Les valeures de X doivent être mesurées avec exactitude.
Avec ce graphiqueResidual vs Leverage , la distance de Cook est calculée pour chaque point afin de mesurer son influence relative sur la régression. Aucune donnée ne se trouve au-dessus de la limite de 0.5. Cela confirme qu’aucun point n’exerce une influence disproportionnée sur la droite de régression.
#---Régression
lm<-lm(Mean.richness~Time.abandon, data=richesse)
#---Distance de cook
cookd<-cooks.distance(lm)
#---Représentation graphique
plot(cookd)
Cette figure confirme qu’aucune mesure possède une influence démesurée sur la régression puisque les valeurs sont inférieures à 1.
Afin de vérifier la validité de cette interprétation des conditions d’application, le test gvlma permet de conclure formellement sur la satisfaction des conditions d’application
#---Test diagnistique gvlma
library(gvlma)
lm<-lm(Mean.richness~Time.abandon, data=richesse)
gvlma(lm)##
## Call:
## lm(formula = Mean.richness ~ Time.abandon, data = richesse)
##
## Coefficients:
## (Intercept) Time.abandon
## 0.5789 2.1080
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = lm)
##
## Value p-value Decision
## Global Stat 1.06109 0.9004 Assumptions acceptable.
## Skewness 0.11025 0.7399 Assumptions acceptable.
## Kurtosis 0.57635 0.4477 Assumptions acceptable.
## Link Function 0.35624 0.5506 Assumptions acceptable.
## Heteroscedasticity 0.01825 0.8926 Assumptions acceptable.Les conditions d’applications sont respectées en raison du verdict “Assumptions acceptable”
c) Analyse des résultats de la régression
#--- Régression
lm<-lm(Mean.richness~Time.abandon, data=richesse)
#---Sommaire des résultats
summary (lm)##
## Call:
## lm(formula = Mean.richness ~ Time.abandon, data = richesse)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.347 -9.195 1.159 8.105 21.721
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5789 4.3525 0.133 0.895
## Time.abandon 2.1080 0.1471 14.329 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.43 on 38 degrees of freedom
## Multiple R-squared: 0.8438, Adjusted R-squared: 0.8397
## F-statistic: 205.3 on 1 and 38 DF, p-value: < 2.2e-16Les éléments de la sortie summary() permettent de conclure les informations suivantes:
β0= 0.58. Lorsqu’aucune année n’est écoulée depuis l’abandon du site, le nombre moyen d’espèces d’oiseaux par point d’écoute est de 0,58 .
β1= 2.10. Lorsque le temps d’abandon augmente d’une année, la richesse en espèces augmente de 2.1080.
C’est 83.97% de la variance de le variable « richesse en espèces » qui est expliquée par la régression, car le coefficient de régression est de 0.8397.
Après 10 ans, la richesse prédite serait de 21.66. Cela se justife par la formule de prédiction suivante:
#---Régression
lm<-lm(Mean.richness~Time.abandon, data=richesse)
#---Prédiction
pred10 <- predict( lm,
newdata = data.frame( Time.abandon = c(10) ),
interval = "prediction" )
show(pred10)## fit lwr upr
## 1 21.65886 -4.274417 47.59213- La force de la corrélation est de 0.8438 et peut être observée dans le summary() ou selon la fonction suivante:
#---Corrélation
summary(lm)$r.squared## [1] 0.8438338d) Résultat de la régression et barres d’erreurs
#---Ajouter le graphique et la droite de régression
plot(Mean.richness~Time.abandon, data=richesse, xlab="Temps depuis l'abandon du site (années)",ylab="Nombre d'espèces",cex.axis=1.2,cex.lab=1.2, main="Nombre moyen d'espèces par point d'écoute selon le temps depuis l'abandon du site", cex.main=.75)
lm<-lm(Mean.richness~Time.abandon, data=richesse)
abline (lm)
#---Niveau de confiance
confint( droite.regression, level = 0.95 )## 2.5 % 97.5 %
## (Intercept) -8.232188 9.389971
## Time.abandon 1.810187 2.405806summary(richesse$Time.abandon)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.0 14.0 26.5 26.4 37.5 49.0#---Modification de l'axe des x
newx<- seq(2.0,49.0, by=0.05)
#---Intervalle de 95% de confiance
conf_interval<- predict(lm,newdata=data.frame(Time.abandon=newx),interval="confidence", level= 0.95)
#---Ajout des barres d'erreur
lines(newx, conf_interval[,2], col="blue")
lines(newx, conf_interval[,3],col= "blue")