Exploring the potential of ETA on different multivariate spaces

Anthony Sturbois, Miquel De Cáceres, Nicolas Djeghri

VivArmor Nature, EMF-CREAF, LEMAR

Outline

Introduction
Taxonomic trajectory analysis
Functional trajectory analysis
Stable isotope trajectory analysis

1. Introduction

About this lesson

Goal

This lesson focuses on the diversity of multivariate \(\Omega_0\) spaces that can be used to answer different ecological questions using ETA concepts and metrics.

Data sets

We focus the presentation on three applications based on taxonomic, functional and stable isotope data sets looking for dynamics though trajectory patterns:

Taxonomic and functional trajectory analysis on the intertidal soft-bottom benthic community in the Bay of Saint-Brieuc (Britanny, France), published in:

Sturbois A., Cormy G., Schaal G., Gauthier O, Ponsero A., Le Mao P., Riera P., Desroy N., 2021, Characterizing spatiotemporal changes in benthic communities: Taxonomic and functional trajectories of intertidal assemblages in the bay of Saint-Brieuc (English Channel). Estuarine, Coastal and Shelf Science, Issue 262, 107603. https://doi.org/10.1016/j.ecss.2021.107603

Stable isotope trajectory analysis for the spatial and temporal resource partitioning in fur seals, published in:

Sturbois, A., Cucherousset, J., De Cáceres, M., Desroy, N., Riera, P., Carpentier, A., Quillien, N., Grall, J., Espinasse, B., Cherel, Y., Schaal, G. (2021). Stable Isotope Trajectory Analysis (SITA) : A new approach to quantify and visualize dynamics in stable isotope studies. Ecological Monographs, 92, e1501. https://doi.org/10.1002/ecm.1501.

Loading libraries

First of all, we load the required libraries, including ecotraj:

library(ecotraj)
library (ape)
library(mapplots)
library(factoextra)
library(FactoMineR)
library(reshape2)
library(sp)

2. Taxonomic trajectory analysis

Taxonomic data set

Description

The taxonomic data set consists in a Hellinger transformed abundance matrix composed of 105 species, 42 entities (sampling sites located in the bay of Saint-Brieuc) and 3 times over a 30 years period (sampling years: 1987-2001-2019).

Loading

We first load the 3 data sets necessary for this application. The taxonomic dataset ‘dataI’:

load(file="../exercises/StudentRdata/dataI.Rdata")

The GPS location of entities ‘loc_entities’:

load(file="../exercises/StudentRdata/loc_entities.Rdata")

And the shapefile for the map representation ‘map’:

load(file="../exercises/StudentRdata/map.Rdata")

Taxonomic space

Here we choose to derive a distance matrix from a PCA points cloud as input for trajectory analysis.

We first perform the PCA, look for the influence of variables (i.e. species) in the two first dimensions and collect the coordinates of the PCA points cloud to derive the distance matrix \(\mathbf{D}\) used for trajectory analysis:

pca <- PCA(as.data.frame(dataI[,-c(1:2)]), ncp=105, graph=FALSE)
var <- get_pca_var(pca)

We then make the factor map:

fviz_pca_var(pca, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

Trajectory definition

We use PCA axes and Euclidean distances to define the \(\Omega_0\) space:

ind <- get_pca_ind(pca)
xy <- matrix(ind$coord, nrow=126,ncol=105)
d <- dist(xy)

Next, we use vectors in dataI to define the trajectory metadata:

entities <-dataI[,1]
times <-dataI[,2]

We can now define the trajectories with the function defineTrajectories() as an object of class trajectory containing the taxonomic distance matrix, and the vectors of entities and times:

x <- defineTrajectories(d, entities, times = times)

Trajectory plot

Let’s have a look at a basic representation of trajectory though the taxonomic \(\Omega_0\) space using the function trajectory PCoA():

trajectoryPCoA(x)

Trajectory metrics

Let’s calculate some trajectory metrics. The idea is to measure taxonomic changes over time with the two consecutive segment length (SL) values (i.e. 1987-2001 and 2001-2019) and to use the net change (NC) values (i.e 1987-2019) to separate recovering from departing patterns. We use the function trajectoryLengths() and fill the argument relativeToinitial to calculate SL (FALSE) or NC (TRUE):

SL <- trajectoryLengths(x)
NC <- trajectoryLengths(x, relativeToInitial = TRUE)

Trajectory clusters

We also aim to define trajectory clusters to look for similar taxonomic dynamics among entities. For that we compute a distance matrix among trajectories with the function trajectoryDistances():

D <- trajectoryDistances(x, distance.type = "DSPD")

Here we performed a dendrogram derived from D, the distance matrix among trajectories using the function hclust(). We then collect entities membership in the vector grp:

Hst <- 25
colstd<-c("lightskyblue4","lightskyblue1","lightskyblue3","navyblue")
cols<-c("lightskyblue1","lightskyblue3","navyblue","lightskyblue4")

hsxy <- hclust(D, "ward.D2")
plot(hsxy,hang = -1, main="distancehellinger 87_19", cex=.6)
y<-rect.hclust (hsxy, h=Hst,
                border = colstd)

Plots with clusters

The trajectory diagram can be now customized using a vector of color depending on the four different trajectory clusters using the vector grp:

grp <- cutree(hsxy, h=Hst)
trajectoryPCoA(x, traj.colors=cols[grp])

Plots with clusters

We can also map the entities and their trajectory clusters memberships:

plot(loc_entities$X,loc_entities$Y,col = cols[grp], bg = cols[grp],pch = 16, cex=2,xlim=c(282570-4400,282570+4400),
     ylim=c(6840592-4400,6840592+4400),xlab="Longitude",ylab="Latitude", main="Taxonomic trajectory clusters (1987,2001,2019)")
plot(map,add=T, col="transparent")

Synthetic maps

The objective now is to represent trajectory metrics though a synthetic maps including the two segment length (i.e. 1987-2001 and 2001-2019) and the net change (i.e. 1987-2019) values with a symbol composed of one circle (NC) and two peripheric triangles (SLs), the top one (SL2001_2019) being directed depending on recovering vs departing pattern.

We start assembling the NC and SL metrics:

SL_NC <- data.frame(Entities = entities[1:42],
                  SL87_01 = SL[,1],
                  SL01_19 = SL[,2],
                  NC87_19 = NC[,2])

Next we calculate the RDT for the distinction between recovering and departing trajectories (RDT):

SL_NC$RDT <- ifelse(SL_NC$SL87_01 - SL_NC$NC87_19 > 0, 1, 2)

We correct the geographical positions for top and down triangles to avoid overlaps. Here the value of 200 was adapted but users should adjust this value depending of their map properties.

loc_entities_upY<-c(loc_entities$Y+200)
loc_entities_dwnY<-c(loc_entities$Y-200)

Synthetic maps

We finally create the map using:

col<-c("black","grey")
cex<-c(0.8,0.8)
pch<-c(24,24)
pch2<-c(25,24)

plot(loc_entities$X,loc_entities$Y,col = "blue",pch = 1 ,cex=SL_NC$NC87_19/7,xlim=c(282570-4400,282570+4400),
     ylim=c(6840592-4400,6840592+4400),xlab="Longitude",ylab="Latitude", main="Taxonomic changes (1987,2001,2019)")
plot(map,add=T, col="transparent")
points(loc_entities$X,loc_entities_dwnY,bg=col[SL_NC$RDT],col =col[SL_NC$RDT] ,
       pch = pch ,cex=SL_NC$SL87_01/7,xlim=c(282570-4400,282570+4400),
       ylim=c(6840592-4400,6840592+4400))
points(loc_entities$X,loc_entities_upY,bg=col[SL_NC$RDT],col = col[SL_NC$RDT],
       pch = pch2[SL_NC$RDT] ,cex=SL_NC$SL01_19/7,xlim=c(282570-4400,282570+4400),
       ylim=c(6840592-4400,6840592+4400))
points(loc_entities$X,loc_entities$Y,col = "blue",pch = 1 ,cex=SL_NC$NC87_19/7,xlim=c(282570-4400,282570+4400),
       ylim=c(6840592-4400,6840592+4400))

Net change ratio maps

We can also represent Pie using Path vs Net changes values to represent the intensity of recovering or departing patterns

X<-rep(loc_entities$X,2)
Y<-rep(loc_entities$Y,2)
Length<-c(rep("Path",42),rep("Net",42))
TL<-as.data.frame(cbind(X,Y,Length))
TL$Dist<-c(SL[,3]-NC[,2],NC[,2])
TL$Length<-Length
xyz <- make.xyz(TL$X,TL$Y,TL$Dist,TL$Length)

plot(loc_entities$X,loc_entities$Y,type="n",xlab="",ylab="", main="Taxonomic Path vs Netchange",xlim=c(282570-4400,282570+4400),
     ylim=c(6840592-4400,6840592+4400))
plot(map,add=T, col="transparent")
draw.pie(xyz$x, xyz$y, xyz$z, radius = 300, col=c("blue","white"))
legend('bottomright',legend="Net changes between 1987 and 2019", pch=15, col="blue", cex=0.75)
legend.z <- round(max(xyz$z,na.rm=TRUE))
legend.bubble("bottomleft", z=legend.z,round=0,maxradius=300,bty="n",txt.cex=0.6, col="darkgrey")
text(278100,6836700,"Path",cex=0.7)

3. Functional trajectory analysis

Functional approach

Taxonomic dynamics do not necessary lead to functional community dynamics. We were interested by this aspect for the intertidal part of the bay of Saint-Brieuc.

Here, we calculate Community Weighted Means (CWM, community trait values weighted by abundance of species) and use this new CWM data set derived from the taxonomic data as input for ETA.

We run functional trajectory analysis with the same method used for taxonomic analysis.

Functional data set

The functional data set was composed of six traits divided in 23 categories. As for the taxonomic data sets, the functional one consists in 42 entities (sampling sites located in the bay of Saint-Brieuc) and three 3 times over a 30 years period (sampling years: 1987-2001-2019)).

Lets load the data set:

load(file="../exercises/StudentRdata/data_traits.Rdata")

These traits characterize the morphology (body size, flexibility, fragility) and behavioral traits (feeding behavior, living habit, tolerance). This set of traits is related to the vulnerability of species to mechanical disturbances (associated to recreational and professional fishing activity and the circulation of vehicles) and organic enrichment.

Functional space

Performing the PCA:

pcaf<-PCA(as.data.frame(data_traits), ncp=23, graph = F)

var <- get_pca_var(pcaf)
fviz_pca_var(pcaf, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

Trajectory definition

Collecting the coordinates of point within the PCA points cloud to define \(\Omega_0\):

ind <- get_pca_ind(pcaf)
xy <- matrix(ind$coord, nrow=126,ncol=23)

Defining the trajectories:

xf <- defineTrajectories(dist(xy), entities, times = times)

Trajectory metrics and trajectory clustering

Calculation of trajectory metrics and dissimilarities among trajectories:

SLf <- trajectoryLengths(xf)
NCf <- trajectoryLengths(xf, relativeToInitial = TRUE)
Df <- trajectoryDistances(xf, distance.type = "DSPD")

We can now define trajectory clusters from Df:

Hst<-15
colstd<-c("chartreuse", "chartreuse3", "chartreuse4","darkgreen")
cols<-c("chartreuse", "chartreuse3", "chartreuse4","darkgreen")

hsxyf <- hclust(Df, "ward.D2")
plot(hsxyf,hang = -1, main="functionnal trajectory clusters 87_19", cex=.6)
yf<-rect.hclust (hsxyf, h=Hst, border =colstd)

Trajectory plot

We use a dendrogram cut to define groups:

grpf <- cutree(hsxyf, h=Hst)

which can be used to distinguish trajectories in the trajectory plot:

trajectoryPCoA(xf, traj.colors=cols[grpf])

Trajectory map

Map of trajectory clusters memberships

plot(loc_entities$X,loc_entities$Y,col = cols[grpf], bg = cols[grpf],pch = 16, cex=2,xlim=c(282570-4400,282570+4400),
     ylim=c(6840592-4400,6840592+4400),xlab="Longitude",ylab="Latitude", main="Functional trajectory clusters (1987,2001,2019)")
plot(map,add=T, col="transparent")

Synthetic maps

Metrics calculations and distinction between recovering and departing trajectories for the synthetic map:

SL_NCf <- data.frame(Entities = entities[1:42],
                     SLf87_01 = SLf[,1],
                     SLf01_19 = SLf[,2],
                     NCf87_19 = NCf[,2])
SL_NCf$RDTf <- ifelse((SL_NCf$SLf87_01 - SL_NCf$NCf87_19) > 0, 1, 2)

Synthetic maps

Synthetic maps including the two segment length and the net change values with the symbol composed of one circle and two triangles:

#coordinate of triangles
loc_entities_upY<-c(loc_entities$Y+200)
loc_entities_dwnY<-c(loc_entities$Y-200)

col<-c("black","grey")
cex<-c(0.8,0.8)
pch<-c(24,24)
pch2<-c(25,24)

plot(loc_entities$X,loc_entities$Y,col = "green",pch = 1 ,cex=SL_NCf$NCf87_19/4,xlim=c(282570-4400,282570+4400),
     ylim=c(6840592-4400,6840592+4400),xlab="Longitude",ylab="Latitude", main="Functional changes (1987,2001,2019)")
plot(map,add=T, col="transparent")
points(loc_entities$X,loc_entities_dwnY,bg=col[SL_NCf$RDTf],col =col[SL_NCf$RDTf] ,
       pch = pch ,cex=SL_NCf$SLf87_01/4,xlim=c(282570-4400,282570+4400),
       ylim=c(6840592-4400,6840592+4400))
points(loc_entities$X,loc_entities_upY,bg=col[SL_NCf$RDTf],col = col[SL_NCf$RDTf],
       pch = pch2[SL_NCf$RDTf] ,cex=SL_NCf$SLf01_19/4,xlim=c(282570-4400,282570+4400),
       ylim=c(6840592-4400,6840592+4400))
points(loc_entities$X,loc_entities$Y,col = "green",pch = 1 ,cex=SL_NCf$NCf87_19/4,xlim=c(282570-4400,282570+4400),
       ylim=c(6840592-4400,6840592+4400))

Net change ratio maps

The second map of pie Path vs Net changes (intensity of recovering or departing functional patterns)

X<-rep(loc_entities$X,2)
Y<-rep(loc_entities$Y,2)
Length<-c(rep("Path",42),rep("Net",42))
TLf<-as.data.frame(cbind(X,Y,Length))
TLf$Dist<-c(SLf[,3]-NCf[,2],NCf[,2])
TLf$Length<-Length
xyz <- make.xyz(TLf$X,TLf$Y,TLf$Dist,TLf$Length)

plot(loc_entities$X,loc_entities$Y,type="n",xlab="",ylab="", main="Functional Path vs Netchange",xlim=c(282570-4400,282570+4400),
     ylim=c(6840592-4400,6840592+4400))
plot(map,add=T, col="transparent")
draw.pie(xyz$x, xyz$y, xyz$z, radius = 300, col=c("green","white"))
legend('bottomright',legend="Net changes between 1987 and 2019", pch=15, col="green", cex=0.75)
legend.z <- round(max(xyz$z,na.rm=TRUE))
legend.bubble("bottomleft", z=legend.z,round=0,maxradius=300,bty="n",txt.cex=0.6, col="darkgrey")
text(278100,6836700,"Path",cex=0.7)

4. Stable isotope trajectory analysis

Fur seal stable isotope dataset

This is a subset of the dataset provided in:

Kernaléguen, L., Arnould, J.P.Y., Guinet, C., Cherel, Y., 2015. Determinants of individual foraging specialization in large marine vertebrates, the Antarctic and subantarctic fur seals. Journal of Animal Ecology 1081–1091.

Briefly, fur seals [the Antarctic fur seal Arctocephalus gazella (AFS) and subantarctic fur seal A. tropicalis (SAFS)] whisker SI values yield unique long-term information on individual behaviour which integrates the spatial, trophic and temporal dimensions of the ecological niche. The foraging strategies of this two species of sympatric fur seals were examined in the winter 2001/2002 at Crozet, Amsterdam and Kerguelen Islands (Southern Ocean) using the stable isotope values of serially sampled whiskers. The method consists in the analysis of consecutive whisker sections (3 mm long) starting from the proximal (facial) end, with the most recently synthesized tissue remaining under the skin. Only individuals (n = 47) with whiskers totalizing at least 30 sections were selected in the initial data, and only those 30 sections were considered herein, from t1 (more recent values) to t30 (oldest values).

We begin by loading the package dataset furseals:

data("furseals")

Trajectory metrics

In this section, we illustrate how to calculate trajectory metrics to characterize the foraging strategy of each fur seal. In the following sections, we show how to use these metrics as data to create plots.

First, we calculate net changes relative to the initial state (i.e. the distance between stable isotope compositions (i.e state) of each whisker section and the initial stable isotope composition) Note that we use here the 2D funtion of ETA for the calculation of trajectory lengths:

Net_changes<-trajectoryLengths2D(furseals[,c("d13C","d15N")],
                                 furseals$ID_SITA,
                                 furseals$Time, relativeToInitial=TRUE) 
head(Net_changes, 3)

      Lt1_t2    Lt1_t3    Lt1_t4    Lt1_t5     Lt1_t6    Lt1_t7    Lt1_t8
1 0.08246211 1.0892314 2.0775851 2.3714985 2.20626585 1.7839633 1.6819920
2 0.59304131 0.9583074 0.8218643 0.2677928 0.08105554 0.7460168 0.4970191
3 0.23256182 1.4622975 3.0029069 2.8817373 2.55906428 1.6445066 1.6734124
    Lt1_t9   Lt1_t10   Lt1_t11   Lt1_t12   Lt1_t13  Lt1_t14   Lt1_t15   Lt1_t16
1 1.344206 1.2735702 1.0631895 0.7915611 0.8486489 1.050064 1.1961910 1.5432880
2 1.092157 0.6772156 0.7744366 0.2702240 0.6136685 1.026186 0.7000864 0.5592397
3 2.518128 2.0782324 1.7294097 1.1594387 1.7649975 2.776892 2.4641059 0.2755014
    Lt1_t17   Lt1_t18   Lt1_t19  Lt1_t20  Lt1_t21   Lt1_t22  Lt1_t23   Lt1_t24
1 1.6484832 1.6042908 1.6566572 1.108358 0.381577 0.3006742 1.477214 1.5890755
2 0.4460235 0.6495922 0.7966461 1.787011 1.591433 1.5776276 1.630589 0.2046069
3 0.4497755 0.6170292 0.2879392 1.103359 2.145421 2.1758332 2.113986 2.2036883
    Lt1_t25  Lt1_t26   Lt1_t27   Lt1_t28   Lt1_t29   Lt1_t30     Path
1 1.0783103 1.117238 2.0638181 2.5204650 2.2934779 1.6320787 40.87543
2 0.3937309 0.384948 0.7032112 0.8949106 0.3901513 0.7480241 21.87682
3 1.8413856 1.635446 1.9473328 1.5905260 1.3147125 1.6315076 49.28114

Trajectory metrics

We then calculate trajectory segment lengths, i.e. the distance between the stable isotope composition of consecutive whisker sections in the stable isotope space:

Segment_lengths<-trajectoryLengths2D(furseals[,c("d13C","d15N")],
                                     furseals$ID_SITA,
                                     furseals$Time, relativeToInitial=FALSE) 
head(Segment_lengths, 3)

          S1        S2       S3        S4        S5        S6        S7
1 0.08246211 1.0799190 1.109759 0.3970806 0.1660482 0.6954279 0.1332817
2 0.59304131 0.4019664 0.165463 0.6287830 0.2497779 0.8109359 0.4377728
3 0.23256182 1.2719328 1.561554 0.1320606 0.8104443 1.2189323 0.3748826
         S8        S9       S10       S11        S12       S13       S14
1 0.3520696 0.1380036 0.2119929 0.3940774 0.08004998 0.2026154 0.2040221
2 0.9468030 0.4481886 0.6290191 0.5362695 0.52697249 0.4552637 0.3411231
3 1.0407238 0.4410726 0.3592548 0.7240691 0.66598649 1.0131856 0.4798427
        S15       S16       S17       S18       S19       S20        S21
1 0.3703039 0.1632330 1.7046178 1.0355916 1.0132186 0.7423544 0.11970798
2 0.1554670 0.2837781 0.2257632 0.6603492 0.9904383 0.5893293 0.03956008
3 2.3198028 0.5979908 0.3078262 0.3993219 1.3487198 1.1301832 0.04002499
         S22       S23       S24       S25       S26       S27       S28
1 1.76281735 0.1147911 0.6148642 0.3162910 0.9585927 0.5338661 0.2521686
2 0.27828942 1.4322639 0.5912495 0.7604426 0.3286716 0.2085977 0.5259743
3 0.09244458 0.1509503 0.6998436 0.3191520 0.6441894 0.5047336 0.3805181
        S29     Path
1 0.6836988 15.63293
2 0.3593800 14.60093
3 0.3359345 19.59814

Identification and characterization of trajectory clusters

Here we aim to define groups of fur seals depending on the similarity of their foraging strategy. We need first to calculate distances between pairs of complete trajectories in the stable isotope space:

dsi <- dist(furseals[,c("d13C","d15N")])
entities <- furseals$ID_SITA
times <- furseals$Time

xsi <- defineTrajectories(dsi, entities, times = times)
  
Ds <- trajectoryDistances(xsi, distance.type = "DSPD")

Then, we can use function hclust() to conduct a hierarchiacal clustering on the symmetric matrix Ds:

colstd<-c("black","yellow","green","blue","grey","red")
pt<-c(16,16,16,16)
hsxy <- hclust(Ds, "ward.D2")
plot(hsxy,hang = -1, main="distance Fur Seals", cex=.6)
Hst=2 # Cutting height
x<-rect.hclust(hsxy, h=Hst,
               border = colstd)

Identification and characterization of trajectory clusters

We cut the dendrogram at height Hst to obtain a vector of cluster membership and copy it in furseals as a factor:

groups <- cutree(hsxy, h=Hst)
furseals$cluster <- as.factor(groups)

Individual trophic trajectories for males and females of A. gazella and A. tropicalis

Here we display trophic trajectories of all individuals while identifying species and gender. Specifically, we create of a vector with the combination of species and gender and create a diagram to display trophic trajectories in the stable isotope space:

furseals$sp_gender<-paste(furseals$Sexe, furseals$Species, sep=" ")
ggplot(data=furseals,aes(x=d13C,y=d15N,color=cluster,shape=Place))+
  geom_point()+
  geom_path(aes(x=d13C,y=d15N,group=ID_SITA,color=cluster),arrow = arrow(length = unit(0.10, "cm")))+
  xlab(expression(delta^13*"C"))+ ylab(expression(delta^15*"N"))+ facet_wrap(~sp_gender) +
  theme_classic()

In each panel, X-Y axes are defined by d13C and d15N stable isotope values. Arrows connects all whiskers section SI states from t1 to t30 (i.e. most recent to oldest SI state). Colors corresponds to trajectory clusters and shape to breeding sites.

Net changes time series for males and females of both fur seal species

In this sub-section we display net changes time series for all individuals, in plots corresponding to combinations of species and gender We prepare a subset of the data called NC:

NC<-Net_changes[,-30]
NC$cluster<-furseals$cluster[1:47]
NC$ID<-as.numeric(rownames(NC))
colnames(NC)<-c(2:30,"cluster","ID")

We then prepare the subset. We notably transform NC in a line structure to compute NCline using the function “melt”, order the data set and add the vector sp_gender:

NCline<-as.data.frame(melt(NC,id=c("ID","cluster")))
colnames(NCline)<-c("ID","Clusters","Time_from_present","Net_changes")
NCline[,3]<-as.numeric(NCline[,3])
NCline[,2]<-as.factor(NCline[,2])
NCline<-NCline[order(NCline[,3],decreasing=F), ]
NCline$sp_gender<-c(furseals$sp_gender[1:47])

Net changes time series for males and females of both fur seal species

We fineally create the plot to display net changes time series for all individuals in panel corresponding Arrows connects all whiskers section stable isotope values from t1 to t30 (i.e. most recent to oldest stable isotope values). Colours corresponds to trajectory clusters:

ggplot(data=NCline,aes(x=Time_from_present,y=Net_changes,color=Clusters))+
  geom_path(aes(x=Time_from_present,y=Net_changes,group=ID,color=Clusters),arrow = arrow(length = unit(0.10, "cm")))+
  facet_wrap(~sp_gender)+
  theme_classic()

M.C. Escher - Sky and Water I, 1938