Introduction to Ecological Quality Assessment (EQA)
Anthony Sturbois/Miquel De Cáceres
2023-10-30
IntroductionEQA.Rmd1. About this vignette
This vignette present the ecological applications used by Sturbois et al. (under review) for the presentation of the Ecological Quality Assessment (EQA) framework. The vignette illustrates how to conduct an Ecological Quality Assessement through state-based and trajectory-based variants, and how to presentation the results graphically. The two data sets used in this vignette have been included in the ecotraj package for reproducibility.
2. Loading libraries
First of all, we load the libraries required for the examples, including ecotraj:
## Loading required package: Rcpp## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-43. Heathland conservation status in Landes et Marais de Glomel
3.1 Context
The nature reserve of Landes et Marais de Glomel (Brittany, France) is composed of temperate Atlantic wet heaths whose reference state is commonly considered dominated by plant communities associated to acid, nutrient poor soils that are at least seasonally water logged and dominated by Erica tetralix and E. ciliaris. This habitat is considered of community interest as part of the European directive on the conservation of natural habitats and of wild fauna and flora (97/62/CEE). For the definition of the reference envelope, surveys integrated the variability of reference states from youngest to early senescent stages, based on expert assessment with respect to the European Directive (97/62/CEE) requirements.
3.2 Loading data
The data set consists of 23 rows (stations) and 46 columns (species). The percent cover values derived from Braun-Blanquet ordinal scale were estimated for 43 species of vascular plants. The first five stations (rows) were used to define the reference envelope, and the next 18 stations (rows) where those for which the conservation status was to be assessed.
## ID Ref Calluna vulgaris Erica ciliaris Erica tetralix Ulex gallii
## 2 1A TRUE 17.5 1 37.5 17.5
## 4 1B TRUE 5.0 0 37.5 5.0
## 6 1C TRUE 0.0 0 17.5 37.5
## 8 5A TRUE 37.5 0 17.5 5.0
## 10 13 TRUE 37.5 0 17.5 1.0
## 12 2 FALSE 5.0 0 37.5 5.0
## Molinia caerulea Betula pubescens
## 2 37.5 1.0
## 4 37.5 0.0
## 6 37.5 0.0
## 8 37.5 0.5
## 10 37.5 5.0
## 12 37.5 1.03.3 Ecological Quality Assessment
We first create a compositional data matrix, by extracting species data from the initial data table:
glomel_comp <- as.matrix(glomel[,!(names(glomel) %in% c("ID", "Ref"))])
rownames(glomel_comp) <- glomel$ID
dim(glomel_comp)## [1] 23 43Using the compositional data, we can use function
vegdist() from package vegan to calculate
of Bray Curtis distances between ecosystem states:
glomel_bc <- vegan::vegdist(glomel_comp, method = "bray")The distance matrix glomel_bc defines the multivariate
space that represents the ressemblance between ecological states.
To conduct EQA, we also need to define the conservation targets.
Column Ref in glomel identifies the stations
considered as reference. We define the state reference
envelope using reference station IDs as follows:
glomel_env <- glomel$ID[glomel$Ref]We are now ready to compare the assessed stations with respect to the
state reference envelope using function
compareToStateEnvelope():
glomel_assess <- compareToStateEnvelope(glomel_bc, glomel_env, m=1.7, distances_to_envelope = TRUE)
head(glomel_assess)## Observation Envelope SquaredDist Q
## 1 1A TRUE 0.03790333 1.0000000
## 2 1B TRUE 0.07132734 0.8314306
## 3 1C TRUE 0.08877350 0.6846406
## 4 5A TRUE 0.04364347 1.0000000
## 5 13 TRUE 0.06748066 0.8701426
## 6 2 FALSE 0.08788512 0.6911247Function compareToStateEnvelope returns
SquaredDist, the squared distance to the envelope centroid,
and the Q statistic, which defines the quality of the
tested stations and is comprised between 0 and 1. We can transform those
values into a qualitative assessment using:
And the assessment for the 18 tested stations is:
glomel_assess[!glomel_assess$Envelope,]## Observation Envelope SquaredDist Q Status
## 6 2 FALSE 0.08788512 0.69112475 Inside
## 7 3A FALSE 0.10690364 0.57054500 Inside
## 8 3B FALSE 0.03193706 1.00000000 Inside
## 9 3C FALSE 0.10365540 0.58868886 Inside
## 10 5B FALSE 0.12638739 0.47819512 Oustside
## 11 6A FALSE 0.32161421 0.15285304 Oustside
## 12 6C FALSE 0.11525495 0.52775487 Inside
## 13 6D FALSE 0.39062446 0.11797552 Oustside
## 14 7A FALSE 0.30326686 0.16512832 Oustside
## 15 8A FALSE 0.17386576 0.33227568 Oustside
## 16 8B FALSE 0.20163739 0.27768273 Oustside
## 17 8C FALSE 0.45790315 0.09515562 Oustside
## 18 9 FALSE 0.21340883 0.25886134 Oustside
## 19 10 FALSE 0.07104764 0.83415862 Inside
## 20 11 FALSE 0.08829721 0.68810525 Inside
## 21 12 FALSE 0.03390049 1.00000000 Inside
## 22 14A FALSE 0.08051443 0.74876145 Inside
## 23 14B FALSE 0.21297195 0.25952180 Oustside3.4 Representing EQA results graphically
To represent the results in a graphical way, we start by conducting a Principal Coordinates Analysis (PCoA) using package ape:
pcoa_glom<-ape::pcoa(glomel_bc)We can use function biplot() from package
ape to show the ordination of stations with the species
projected as arrows:
biplot(pcoa_glom, glomel_comp)
To generate a different representation, we first get coordinates of stations in the two first dimensions of the PCoA and copy them into the assessment data frame:
PCOA_DIM1_2 <- pcoa_glom[["vectors"]][,1:2]
glomel_assess$Dim1<-PCOA_DIM1_2[,1]
glomel_assess$Dim2<-PCOA_DIM1_2[,2]
head(glomel_assess)## Observation Envelope SquaredDist Q Status Dim1 Dim2
## 1 1A TRUE 0.03790333 1.0000000 Inside -0.2185262 0.05500827
## 2 1B TRUE 0.07132734 0.8314306 Inside -0.1897143 0.19276850
## 3 1C TRUE 0.08877350 0.6846406 Inside -0.1123913 0.02290741
## 4 5A TRUE 0.04364347 1.0000000 Inside -0.1842064 -0.14315651
## 5 13 TRUE 0.06748066 0.8701426 Inside -0.1594965 -0.17397305
## 6 2 FALSE 0.08788512 0.6911247 Inside -0.1840068 0.20520893We then use package ggplot to achieve a nicer representation of the assessment. Symbol colors are used to illustrate stations that do achieve conservation objectives (in green) from those that do not achieve them (in red). The size of dots represents the squared distance to the centroid of the state reference envelope:
p<-ggplot(glomel_assess,
mapping=aes(x=Dim1,y=Dim2,size=SquaredDist, color=Status, shape=Status))+
geom_point()+
scale_colour_manual(values=c("#00BFC4", "#F8766D"))+
geom_text(glomel_assess,
mapping=aes(x=Dim1,y=Dim2,label=Observation),
hjust=-0.3, vjust=-0.3,size=2.5, color="Black")+
xlab(expression("PCoA1 (27.7%)"))+
ylab(expression("PCoA2 (13.7%)"))
p<-p + theme_minimal()+
theme(
# Hide some graphical elements
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))
p
4. Impact of fishing activities on marine habitats
4.1 Context
We used an experimental data set built by Tauran et al. (2020) to study the impact of fishing dredges and varying fishing pressures on maerl beds, in the bay of Brest (Brittany, France). The study follows a Before-After-Control-Impact design (Stewart-Oaten et al., 1986). Briefly, three control stations were surveyed and compared to five treatment stations characterized by different fishing dredges and pressure levels (Tauran et al., 2020). Three fishing dredges are compared: (1) a clam dredge (CD), 70 to 90 kg, 1.5 m wide, 40 teeth of 11 cm each; (2) a queen scallop dredge (QSD), 120 kg,1.8 m wide, with a blade; and (3) a king scallop dredge (KSD), 190 kg, 1.8 m wide, 18 teeth of 10 cm each every 9 cm. Pressure levels were measured as the number of dredge tows performed on the zone during the experimental dredging session: 0 (i.e. control), 10, or 30 dredge tows. Samples were collected from April 2016 to April2017. Session 1 was sampled just before the experimental dredging (t0); session 2 at t0 + 1 week; session 3 at t0 + 1 month and session 4 at t0 + 12 months. Nine replicates were sampled for all treatments and sessions with a Smith-McIntyre grab (0.1m2). Abundance data were pooled at the treatment levels (i.e. CTRL1, CTRL2, CTRL3, KSD_10, CD_10, CD_30, QSD_10 and QSD_30; Appendix_D).
4.2 Loading data
The data set has 32 observations (i.e. treatments*surveys) and 250 species. We begin by loading the data:
## Abundance.Abludomelita gladiosa Abundance.Abra alba
## CD_10-S1 60 15
## CD_10-S2 54 14
## CD_10-S3 42 15
## CD_10-S4 45 8
## CD_30-S1 58 18
## CD_30-S2 39 19
## Abundance.Acanthochitona fascicularis Abundance.Achelia hispida
## CD_10-S1 2 0
## CD_10-S2 1 0
## CD_10-S3 0 0
## CD_10-S4 0 0
## CD_30-S1 4 0
## CD_30-S2 3 0
## Abundance.Alcyonidium diaphanum Abundance.Alentia gelatinosa
## CD_10-S1 0 1
## CD_10-S2 0 0
## CD_10-S3 0 0
## CD_10-S4 0 0
## CD_30-S1 0 0
## CD_30-S2 0 04.3 Definition of the ressemblance space
We separate the columns that identify the treatment and surveys.
#vector Treatment
Treatment<-glenan$Treatment
#vector Surveys
Surveys<-glenan$SurveysWe also extract the compositional data and log-transform it as initially performed in Tauran et al. (2020).
Using the compositional data, we can use function
vegdist() from package vegan to calculate
of Bray Curtis distances between ecosystem states:
grab_bc <- vegan::vegdist(grab_comp, method = "bray")Like in the first example, the distance matrix grab_bc
defines the multivariate space that represents the ressemblance between
ecological states of maerl beds.
4.4 State-based EQA
In addition to the resemblance space, to conduct the EQA we need to define the conservation targets. In this case, we define the states conforming our state reference envelope manually among the observation IDs:
#Definition of the state reference envelope
ID<-rownames(grab_comp)
grab_env <- ID[c(1,9:20,21,25,29)]We are now ready to compare the assessed stations with respect to the
state reference envelope using function
compareToStateEnvelope():
#Comparing assessed stations with respect to the state reference envelope
grab_assess <- compareToStateEnvelope(grab_bc, grab_env, m=1.5,
distances_to_envelope = TRUE)As before we derive a qualitative assessment:
grab_assess$Status<-c(ifelse(grab_assess$Q>=0.5,"Inside", "Outside"))
grab_assess$Status[c(1,9:20,21,25,29)]<-"Reference"
grab_assess$Status<- factor(grab_assess$Status, levels = c("Reference","Inside","Outside"))
head(grab_assess)## Observation Envelope SquaredDist Q Status
## 1 CD_10-S1 TRUE 0.03374086 0.5587631 Reference
## 2 CD_10-S2 FALSE 0.02316547 0.9025909 Inside
## 3 CD_10-S3 FALSE 0.04174879 0.4041249 Outside
## 4 CD_10-S4 FALSE 0.03802341 0.4677650 Outside
## 5 CD_30-S1 FALSE 0.02940983 0.6757536 Inside
## 6 CD_30-S2 FALSE 0.04151217 0.4078031 OutsideAs before, we can represent the results of the state-based EQA on in the two first dimensions of a Principal Coordinates Analysis:
#Get coordinates of stations in the two first dimensions of the PCoA
pcoa<-pcoa(grab_bc)
PCOA_DIM1_2<-pcoa[["vectors"]][,1:2]
grab_assess$Dim1<-PCOA_DIM1_2[,1]
grab_assess$Dim2<-PCOA_DIM1_2[,2]
#Information about Treatment and Surveys
grab_assess$Treatment<-Treatment
grab_assess$Surveys<-Surveys
p<-ggplot(grab_assess, aes(x=Dim1, y=Dim2, group=Status)) +
geom_point(aes(shape=Status, color=Status, size=SquaredDist))+
scale_shape_manual(values=c(1, 16, 17))+
scale_color_manual(values=c("#00BFC4", "#00BFC4", "#F8766D"))+
geom_path(aes(x=Dim1,y=Dim2,group=Treatment),color="grey", arrow = arrow(length = unit(0.03, "cm")))+
geom_text(grab_assess,mapping=aes(x=Dim1,y=Dim2,label=ID),hjust=-0.3, vjust=-0.3,size=2.5, color="Black")+
xlab(expression("PCoA1 (13.3%)"))+
ylab(expression("PCoA2 (9.9%)"))
p<-p+ theme_minimal()+
theme(
# Hide some graphical elements
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"))
p
Since we are dealing with temporal data, we can create a bar plot showing the evolution of the squared distance to the centroid of the state reference envelope during the study period for each treatment:
grab_assess$Treatment<-factor(grab_assess$Treatment,c("CTRL1","CTRL2","CTRL3",
"CD_10","CD_30",
"QSD_10","QSD_30",
"KSD_10"))
grab_assess$Status<-c(ifelse(grab_assess$Q>=0.5,"Inside", "Outside"))
grab_assess$Status<- factor(grab_assess$Status, levels = c("Inside","Outside"))
ggplot(grab_assess, aes(x=Surveys, y=SquaredDist, color=Status, fill=Status)) +
geom_bar(stat = "identity",width=0.5)+
scale_colour_manual(values=c("#00BFC4", "#F8766D"))+
scale_fill_manual(values=c("#00BFC4", "#F8766D"))+
scale_y_continuous("Squared distance to the reference envelope")+
theme_minimal()+
facet_wrap(~Treatment)+
theme(axis.text.x = element_text(angle=90))
In addition, we can create the heatmap showing the evolution of the distance between control and dredged stations during the experimentation:
grabbc_sel <- melt(as.matrix(grab_bc),rnames = c("row", "col"))
#selection of the observation of interest (i.e. control stations for heat map row, and dredged stations for heat map column)
grabbc_sel<-grabbc_sel[!grabbc_sel$Var2 %in% c("CTRL1-S1","CTRL2-S1","CTRL3-S1","CTRL4-S1",
"CTRL1-S2","CTRL2-S2","CTRL3-S2","CTRL4-S2",
"CTRL1-S3","CTRL2-S3","CTRL3-S3","CTRL4-S3",
"CTRL1-S4","CTRL2-S4","CTRL3-S4","CTRL4-S4"),]
grabbc_sel<-grabbc_sel[grabbc_sel$Var1 %in% c("CTRL1-S1","CTRL2-S1","CTRL3-S1","CTRL4-S1",
"CTRL1-S2","CTRL2-S2","CTRL3-S2","CTRL4-S2",
"CTRL1-S3","CTRL2-S3","CTRL3-S3","CTRL4-S3",
"CTRL1-S4","CTRL2-S4","CTRL3-S4","CTRL4-S4"),]
grabbc_sel$Var1<-factor(grabbc_sel$Var1,
c("CTRL1-S1","CTRL2-S1","CTRL3-S1","CTRL4-S1",
"CTRL1-S2","CTRL2-S2","CTRL3-S2","CTRL4-S2",
"CTRL1-S3","CTRL2-S3","CTRL3-S3","CTRL4-S3",
"CTRL1-S4","CTRL2-S4","CTRL3-S4","CTRL4-S4"))
grabbc_sel$Var2<-factor(grabbc_sel$Var2,
c("KSD_10-S1","CD_10-S1","CD_30-S1","QSD_10-S1","QSD_30-S1",
"KSD_10-S2","CD_10-S2","CD_30-S2","QSD_10-S2","QSD_30-S2",
"KSD_10-S3","CD_10-S3","CD_30-S3","QSD_10-S3","QSD_30-S3",
"KSD_10-S4","CD_10-S4","CD_30-S4","QSD_10-S4","QSD_30-S4"))
ggplot(grabbc_sel, aes(x=Var2, y=Var1, fill=value)) +
geom_raster() +
scale_fill_viridis_c()+
theme(axis.text.x=element_text(angle = 90, hjust = 0))
4.5 Trajectory-based EQA
One of the nice features of EQA is the possibility to conduct the quality assessment taking into account ecosystem dynamics. In this case, the conservation target is defined by a set of trajectories conforming a trajectory reference envelope, which in our case are the three control stations:
#Definition of the trajectory reference envelope
grab_dynenv<-c("CTRL1","CTRL2","CTRL3")We are now ready to compare the assessed trajectories (those subject
to treatment) with respect to the trajectory reference envelope using
function compareToTrajectoryEnvelope():
##Comparing assessed trajectories with respect to the trajectory reference envelope
grab_dynassess <- compareToTrajectoryEnvelope(grab_bc,Treatment, grab_dynenv,
m=1.5, distances_to_envelope = TRUE)As before, we can create a bar plot showing the squared distance to the centroid of the trajectory reference envelope during the study period for each treatment:
grab_dynassess$Treatment<-factor(grab_dynassess$Site,c("CTRL1","CTRL2","CTRL3",
"CD_10","CD_30",
"QSD_10","QSD_30",
"KSD_10"))
grab_dynassess$Status<-c(ifelse(grab_dynassess$Q>=0.5,"Inside", "Outside"))
grab_dynassess$Status<- factor(grab_dynassess$Status, levels = c("Inside","Outside"))
ggplot(grab_dynassess, aes(x=reorder(Treatment, SquaredDist),
y=SquaredDist, color=Status, fill=Status)) +
geom_bar(stat = "identity",width=0.5)+
scale_colour_manual(values=c("#00BFC4", "#F8766D"))+
scale_fill_manual(values=c("#00BFC4", "#F8766D"))+
scale_y_continuous("Squared distance to the reference envelope")+
theme_minimal()+
theme(axis.text.x = element_text(angle=90))
5. References
Sturbois, A., De Cáceres, M., Bifolchi, A., Bioret, F., Boyé, A., Gauthier, O., Grall, J., Grémare, A., Labrune, C., Robert, A., Schaal, G., Desroy, N. (2023). Ecological Quality Assessment: a general multivariate framework to report the quality of ecosystems and their dynamics with respect to reference conditions. Ecosphere.
Tauran, A., Dubreuil, J., Guyonnet, B., Grall, J., 2020. Impact of fishing gears and fishing intensities on maerl beds: An experimental approach. Journal of Experimental Marine Biology and Ecology 533, 151472. https://doi.org/10.1016/j.jembe.2020.151472