Computation of alpha functional indices: A temporal approach using the mFD package

Carlos E. Paz-Ríos
Postdoctoral fellow at Universidad Autonoma de Campeche
Instituto de Ecologia, Pesquerias y Oceanografia del Golfo de Mexico

September 29, 2023

Abstract

The present analysis sequence follows the ‘General Workflow’ by Magneville et al. (2022). For more details on using the mFD package, see this website. This guide is intended to shows the basic procedures for determining alpha functional indices, so interpretation of results shall be at minimum.

1. Load data

A simple but important step for data processing in R, set working directory.

Traits

The data set used here for fish traits was taken from Paz-Ríos et al. (2023); these are described and freely available in a zenodo open data repository, click here. This data set should be in the form of a data frame or matrix, where species are in rows and traits are in columns.

traits <- read.csv('fish_traits.csv')
rownames(traits) <- traits$'species'
traits <- traits[ , -1]

Define nature of traits (i.e. factor, ordinal, or fuzzy)

Samples

The study design and fish community descriptors used here are described in Paz-Ríos et al. (2022); to redirect to this paper, follow this link. This data set must be in the form of a matrix, where samples are in rows and species are in columns.

samples <- read.csv('fish_sample.csv') 
rownames(samples) <- samples$samlple
samples <- samples[ , -c(1)]
samples <- as.matrix(samples)

Trait type

Define trait type (i.e. numeric, factor, ordinal, or fuzzy)

type <- read.csv('fish_types.csv')

2. Processing traits data

Only the datframe of traits is used in this section.

Species traits summary

fish_traits_summ <- mFD::sp.tr.summary(
  tr_cat     = type,   
  sp_tr      = traits, 
  stop_if_NA = TRUE)
fish_traits_summ

Compute trait-based distances

sp_dist_fish <- mFD::funct.dist(
  sp_tr         = traits,
  tr_cat        = type,
  metric        = 'gower',
  scale_euclid  = 'scale_center',
  ordinal_var   = 'metric',
  weight_type   = 'equal',
  stop_if_NA    = TRUE)

Computing functional spaces and their quality

fspaces_quality_fish <- mFD::quality.fspaces(
  sp_dist             = sp_dist_fish,
  maxdim_pcoa         = 10,
  deviation_weighting = 'absolute',
  fdist_scaling       = TRUE,
  fdendro             = 'average')

## Registered S3 method overwritten by 'dendextend':
##   method     from 
##   rev.hclust vegan

Check the ‘mad’ value to inspect the quality of spaces

round(fspaces_quality_fish$'quality_fspaces', 3)

Plot functional spaces and their quality

mFD::quality.fspaces.plot(
  fspaces_quality            = fspaces_quality_fish,
  quality_metric             = 'mad_scaled',
  fspaces_plot               = c('tree_average', 'pcoa_1d', 'pcoa_2d', 
                                 'pcoa_3d', 'pcoa_6d'),
  name_file                  = NULL,
  range_dist                 = NULL,
  range_dev                  = NULL,
  range_qdev                 = NULL,
  gradient_deviation         = c(neg = 'darkblue',
                                 nul = 'grey80',
                                 pos = 'darkred'),
  gradient_deviation_quality = c(low = 'yellow',
                                 high = 'red'),
  x_lab                      = 'Trait-based distance')

Extract the values of the PCOs axes

sp_faxes_coord_fish <- fspaces_quality_fish$'details_fspaces'$'sp_pc_coord'

Plot functional space

big_plot <- mFD::funct.space.plot(
  sp_faxes_coord  = sp_faxes_coord_fish[ , c('PC1', 'PC2',
                                             'PC3', 'PC6')],
  faxes           = c('PC1', 'PC2', 'PC3', 'PC6'),
  name_file       = NULL,
  faxes_nm        = NULL,
  range_faxes     = c(NA, NA),
  color_bg        = 'grey95',
  color_pool      = 'darkgreen',
  fill_pool       = 'white',
  shape_pool      = 21,
  size_pool       = 1,
  plot_ch         = TRUE,
  color_ch        = 'black',
  fill_ch         = 'white',
  alpha_ch        = 0.5,
  plot_vertices   = TRUE,
  color_vert      = 'blueviolet',
  fill_vert       = 'blueviolet',
  shape_vert      = 23,
  size_vert       = 1,
  plot_sp_nm      = NULL,
  nm_size         = 3,
  nm_color        = 'black',
  nm_fontface     = 'plain',
  check_input     = TRUE)

Plotting

big_plot$patchwork

3. Compute functional alpha diversity indices

Here the extracted values of the PCOs axes (the multidimensional space) and structure of samples are used.

alpha_fd_indices_fish <- mFD::alpha.fd.multidim(
  sp_faxes_coord   = sp_faxes_coord_fish[ , c('PC1', 'PC2',
                                              'PC3', 'PC6')],
  asb_sp_w         = samples,
  ind_vect         = c('fdis', 'fmpd', 'fnnd',
                       'feve', 'fric', 'fdiv',
                       'fori', 'fspe', 'fide'),
  scaling          = TRUE,
  check_input      = TRUE,
  details_returned = TRUE)

## 1980_1981feb done 2.2%

## 1980_1981mar done 4.3%

## 1980_1981apr done 6.5%

## 1980_1981may done 8.7%

## 1980_1981jun done 10.9%

## 1980_1981jul done 13%

## 1980_1981aug done 15.2%

## 1980_1981sep done 17.4%

## 1980_1981oct done 19.6%

## 1980_1981nov done 21.7%

## 1980_1981dec done 23.9%

## 1980_1981jan done 26.1%

## 1998_1999feb done 28.3%

## 1998_1999mar done 30.4%

## 1998_1999apr done 32.6%

## 1998_1999may done 34.8%

## 1998_1999jun done 37%

## 1998_1999jul done 39.1%

## 1998_1999aug done 41.3%

## 1998_1999sep done 43.5%

## 1998_1999oct done 45.7%

## 1998_1999nov done 47.8%

## 1998_1999dec done 50%

## 1998_1999jan done 52.2%

## 2010_2011nov done 54.3%

## 2010_2011dec done 56.5%

## 2010_2011jan done 58.7%

## 2010_2011feb done 60.9%

## 2010_2011mar done 63%

## 2010_2011apr done 65.2%

## 2010_2011may done 67.4%

## 2010_2011jun done 69.6%

## 2010_2011jul done 71.7%

## 2010_2011aug done 73.9%

## 2010_2011sep done 76.1%

## 2010_2011oct done 78.3%

## 2016_2017sep done 80.4%

## 2016_2017oct done 82.6%

## 2016_2017dec done 84.8%

## 2016_2017feb done 87%

## 2016_2017mar done 89.1%

## 2016_2017apr done 91.3%

## 2016_2017may done 93.5%

## 2016_2017jun done 95.7%

## 2016_2017jul done 97.8%

## 2016_2017aug done 100%

Extract values as data frame

fd_ind_values_fish <- alpha_fd_indices_fish$'functional_diversity_indices'
head(fd_ind_values_fish)

Save data frame

write.csv(fd_ind_values_fish, file = 'res_alpha_mFD.csv')

4. Exploring trends in alpha diversity indices

Hereafter, the analysis sequence was adapted to examine the distribution of functional indices regarding a specific sample structure. To do it, a new data set was called to include classification criteria (factors) and environmental variables. It is also included the Taxonomic distinctness analyzed in Paz-Ríos et al. (2022).

factoenvi <- read.csv('fish_envi.csv')
head(factoenvi)

Merge data frames (functional indices + environment)

dive_envi <- cbind(fd_ind_values_fish, factoenvi)
names(dive_envi)

Assign factors

dive_envi$campaing <- as.factor(dive_envi$campaing)
dive_envi$season <- as.factor(dive_envi$season)

Then, call packages to be used, and set a general kind of plot features

library(plyr)
library(ggplot2)
library(ggpubr)

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:plyr':
## 
##     mutate

theme_set(
  theme_bw() +
    theme(legend.position = 'top')
)

Plotting functional indices versus different sampling years

p1 <- ggplot(data = dive_envi, aes(x=campaing, y=sp_richn))+
  rremove('x.text') +
  geom_boxplot(color="blue", fill="lightblue", alpha=0.2)

p2 <- ggplot(data = dive_envi, aes(x=campaing, y=Taxo_dis))+
  rremove('x.text') +
  geom_boxplot(color="blue", fill="lightblue", alpha=0.2)

p3 <- ggplot(data = dive_envi, aes(x=campaing, y=fric))+
  rremove('x.text') +
  geom_boxplot(color="blue", fill="lightblue", alpha=0.2)

p4 <- ggplot(data = dive_envi, aes(x=campaing, y=fdis))+
  geom_boxplot(color="blue", fill="lightblue", alpha=0.2)

p5 <- ggplot(data = dive_envi, aes(x=campaing, y=feve))+
  geom_boxplot(color="blue", fill="lightblue", alpha=0.2)

p6 <- ggplot(data = dive_envi, aes(x=campaing, y=fdiv))+
  geom_boxplot(color="blue", fill="lightblue", alpha=0.2)

plot1 <- ggarrange(p1, p2, p3, p4, p5, p6, 
                   labels = c('A', 'B', 'C', 'D', 'E', 'F'),
                   ncol = 3, nrow = 2)
plot1

Biplots

A correlation between the functional richness and the species richness and taxonomic distinctness was performed.

Relationship between Species richnnes vs Functional richnnes

cor.test(dive_envi$sp_richn, dive_envi$fric)

## 
##  Pearson's product-moment correlation
## 
## data:  dive_envi$sp_richn and dive_envi$fric
## t = 8.046, df = 44, p-value = 3.466e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6202185 0.8675636
## sample estimates:
##       cor 
## 0.7715927

p7 <- ggplot(dive_envi, aes(sp_richn, fric)) +
  geom_point()
p7 <- p7 + geom_smooth(method = lm)

data_sum <- ddply(dive_envi, c('campaing', 'season'), summarise,
                  sp_richnmn=mean(sp_richn), 
                  sp_richnsd=sd(sp_richn)/sqrt((length(sp_richn))),
                  fricmn=mean(fric),
                  fricsd=sd(fric)/sqrt((length(fric))))

Ylims<-aes(ymax= fricmn+fricsd, ymin= fricmn-fricsd)
Xlims<-aes(xmax= sp_richnmn+sp_richnsd, xmin= sp_richnmn-sp_richnsd)

p8 <-ggplot(data_sum, aes(x=sp_richnmn, y=fricmn,
                          colour=season, shape =season)) + 
  geom_point(size=3) + 
  geom_errorbar(Ylims, width=0.01) + 
  geom_errorbarh(Xlims, height=0.01) +
  ylab(expression('Functional richness')) +
  xlab(expression('Species richness')) +
  facet_grid(. ~ campaing)+
  theme_bw()

plot2 <- ggarrange(p7, p8, 
                   labels = c('A', 'B'),
                   ncol = 1, nrow = 2)

## `geom_smooth()` using formula = 'y ~ x'

plot2

Relationship between Taxonomic distinctness vs Functional richness

cor.test(dive_envi$Taxo_dis, dive_envi$fric)

## 
##  Pearson's product-moment correlation
## 
## data:  dive_envi$Taxo_dis and dive_envi$fric
## t = -0.90918, df = 44, p-value = 0.3682
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4099325  0.1608444
## sample estimates:
##        cor 
## -0.1357943

p9 <- ggplot(dive_envi, aes(Taxo_dis, fric)) +
  geom_point()
p9 <- p9 + geom_smooth(method = lm)

data_sum.tax<-ddply(dive_envi, c('campaing', 'season'), summarise,
                    tax_dismn=mean(Taxo_dis), 
                    tax_dissd=sd(Taxo_dis)/sqrt((length(Taxo_dis))),
                    fricmn=mean(fric),
                    fricsd=sd(fric)/sqrt((length(fric))))

Ylims<-aes(ymax= fricmn+fricsd, ymin= fricmn-fricsd)
Xlims<-aes(xmax= tax_dismn+tax_dissd, xmin= tax_dismn-tax_dissd)

p10 <-ggplot(data_sum.tax, aes(x=tax_dismn, y=fricmn, colour=season, shape =season)) + 
  geom_point(size=3) + 
  geom_errorbar(Ylims, width=0.01) + 
  geom_errorbarh(Xlims, height=0.01) +
  ylab(expression('Functional richness')) +
  xlab(expression('Taxonomic distinctness')) +
  facet_grid(. ~ campaing)+
  theme_bw()

plot3 <- ggarrange(p9, p10, 
                   labels = c('A', 'B'),
                   ncol = 1, nrow = 2)

## `geom_smooth()` using formula = 'y ~ x'

plot3

5. Trait-environment relationships

The next section is based on the tutorial from Frelat et al. (2022), which is freely available in a zenodo open data repository, click here. Analysis rely on the mgcv package that perform Generalized Additive Models (GAMs) by building here a function to display GAMs partial smooth plots and return summary statistics for the fitted model.

library(mgcv)

## Loading required package: nlme

## This is mgcv 1.8-41. For overview type 'help("mgcv-package")'.

funGam<-function(Var){
  colnames(dive_envi)[which(colnames(dive_envi)==Var)]<-"Var"
  fit<-gam(Var~s(Dep,k=3)+
             s(Secchi,k=3)+
             s(Tem,k=3)+
             s(Sal,k=3)+
             s(pH,k=3)+
             s(Oxy, k=3),
           na.action=na.exclude, data=dive_envi)
  par(mfrow=c(4,2),mar=c(4,4,0.2,0.2))
  plot(fit,shade=T,shade.col="grey",res=T,rug=F,pch=20,ylab="")
  par(mfrow=c(1,1),mar=c(5,4,2,2))
  return(summary(fit))
}

Make a list of functional indices from mFD calculation

indices <- colnames(fd_ind_values_fish)
indices

##  [1] "sp_richn" "fdis"     "fmpd"     "fnnd"     "feve"     "fric"    
##  [7] "fdiv"     "fori"     "fspe"     "fide_PC1" "fide_PC2" "fide_PC3"
## [13] "fide_PC6"

Modelling functional richness (fric) for instance

funGam(indices[6])

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Var ~ s(Dep, k = 3) + s(Secchi, k = 3) + s(Tem, k = 3) + s(Sal, 
##     k = 3) + s(pH, k = 3) + s(Oxy, k = 3)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.45793    0.02016   22.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df     F p-value   
## s(Dep)    1.040  1.077 0.866 0.33507   
## s(Secchi) 1.000  1.000 4.015 0.05207 . 
## s(Tem)    1.000  1.000 2.032 0.16201   
## s(Sal)    1.405  1.645 6.608 0.00364 **
## s(pH)     1.000  1.000 8.429 0.00605 **
## s(Oxy)    1.000  1.000 0.029 0.86617   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.355   Deviance explained = 44.7%
## GCV = 0.022307  Scale est. = 0.018697  n = 46

6. Acknowledgments

The present analysis sequence forms part of a postdoctoral research funded by the National Council of Humanities, Sciences and Technologies, Mexico (CONAHCYT).

7. References

Frelat R, Beukhof E, Weigel B, Lindegren M. 2022, June 23. Tutorial for analysing trait-environment relationships. Zenodo. https://doi.org/10.5281/zenodo.6712534

Magneville C, Loiseau N, Albouy C, Casajus N, Claverie T, Escalas A, Leprieur F, Maire E, Mouillot D, Villéger S. 2022. mFD: an R package to compute and illustrate the multiple facets of functional diversity. Ecography. https://onlinelibrary.wiley.com/doi/10.1111/ecog.05904

Paz-Ríos CE, Sosa-López A, Torres-Rojas YE, del Río-Rodríguez RE. 2022. Long-term multiscale analysis on temporal variability of the fish community in Terminos Lagoon. Estuarine Coastal and Shelf Science, 277: 108066. https://doi.org/10.1016/j.ecss.2022.108066

Paz-Ríos CE, Aguilar-Medrano R, Sosa-López A. 2023. Functional traits of the ichthyofauna from Terminos Lagoon, Campeche, Mexico [Data set]. Zenodo. https://doi.org/10.5281/zenodo.8367849