Creating functional groups of marine fish from categorical traits

Select the functional group to be defined

The type of functional group defined will be dependent on the ecosystem that is being modelled. Different ecosystems require different functions in order for their production to be exploited by its inhabitants (Fonseca & Ganade, 2001). For example, coral reef fishes need strong, sharp teeth in order to exploit polyps, while large pelagic species need to be fast moving in order to capture prey. Functional groups of species should be defined by how the species use their environment and its resources as ecosystem models attempt to model the entire process of an ecosystem spatially and temporally (Fulton et al., 2004). As we are modelling an open ecosystem, where species can enter and leave, it is important to try and capture some of the diversity of how species use an ecosystem daily, seasonally and yearly. The final groupings of species should exhibit similar responses to environmental conditions and have similar effects on the ecosystem processes (Fonseca & Ganade, 2001), though a good way to test these characteristics is yet to be found.

Select species to include

The species selected to include should represent the taxonomy, time and space that the functional groups are trying to capture (Fonseca & Ganade, 2001). That is, species that rarely occupy the area of interest, or species with greatly differeing biomasses should be included in the analysis. This is because including many species in functional groups better explains changes in the biodiversity of a given system (Naeem & Li, 1997). For this study, a comprehensive list of species of fishes from TBGB was made from the latest published account of trawl data (Stevenson & MacGibbon, 2015) and from published accounts of species known to inhabit TBGB (Roberts, Stewart & Struthers, 2015). While there are obvious functional differences between adults and juveniles of many species, that should be addressed and incorporated, such a delineation was beyond the scope of this project.

Select functions of interest

To avoid functional redundancy more functions can be selected to increase the chances of species having unique roles within the ecosystem, while ensuring that species who display the same traits across a number of functions truly belong to the same functional group. We selected four different functions to represent how the species of interest utilise their environment: diet, morphology, habitat use and life history traits (Villéger, Maire & Leprieur, 2017; Gravel, Albouy & Thuiller, 2016; Costello et al., 2015). Diet determines a species influence on other organisms in the environment and its position in the food web (Costello et al., 2015). Habitat preference allow us to understand how the different species might aggregate in the environment and can provide information about the likely lifestyle of the species (Chan, 2001; Vadas Jr & Orth, 1997). Morphology traits are important in defining the range of food sources, behaviour, adaptation and habitat use available to a certain species (Sibbing & Nagelkerke, 2000). Life history primarily reflected the reproductive strategies of the species which may be indicative of their abundance and resilience in the environment (Villéger, Maire & Leprieur, 2017).

Trait selection

Traits should reflect the functions of interest. A literature review was conducted that identified 94 potential traits that could be recorded from fish species. As cost and time are often significant motivators for conducting research, it was a goal of this study to record functional trait information only from published resources or from photographs, rather than collecting and measuring specimens. We identified 40 traits that could be recorded without measuring species directly (Table S1). For some cases, variables that previously required a specimen to be measured were able to be categorised into nominal variables. For example, caudal peduncle aspect ratio was recorded as caudal fin shape. Where information differed ontogenetically within species, the information for adult females was recorded. The final list of recorded traits is provided in Table 1.

Morphology traits describe how species move around their environment and can potentially be used as an indicator of prey preferences (Albouy et al., 2011). Most of the traits recorded for morphology were determined from pictures of the species. Descriptions of the species fins were recorded either as their position on the body (pelvic), the shape of the fin (caudal) or the fin composition (soft ray or spines—dorsal). The shape of the caudal fin is important in determining the ability of a species to transition between vertical habitats (Bridge et al., 2016). The swimming mode of the species was recorded as either body caudal fin (BCF) locomotion or median paired fin (MPF) locomotion that is an indicator of the evasiveness of the food types targeted (Sfakiotakis, Lane & Davies, 1999; Webb, 1984b). The body form of the species was recorded as either fusiform, flat, cylindrical or compressed which is an indicator of how species acquire their food (Webb, 1984a). Eye position indicates the likely location of the species in the water column (Mindel et al., 2016). The spiny dorsal fin type may be an indicator of protection (i.e., from the number of spines—another recorded variable) but can also indicate the manoeuvrability of the species. The soft dorsal fin can help a fish to remain stable while swimming but is also able to generate thrust (Lauder & Drucker, 2004). Oral gape position can indicate feeding position in the water column (Albouy et al., 2011) and prey types that may be acquired (Zhao et al., 2014). Teeth shape indicate the type of prey consumed and the substrate on which a species may be feeding (Bellwood et al., 2014). Body length is an indicator of potential prey available and it correlates with size at maturity, fecundity, growth rate and longevity (Sibbing & Nagelkerke, 2000; Mindel et al., 2016). Physical protection was recorded as present or absent as an indicator of how difficult the species would be to use as prey (Reecht et al., 2013).

The life history traits selected primarily reflect the reproductive strategies for each of the species. Parental care (care, no care) was included as it can indicate where a species chooses to breed as well as the size and amount of the offspring (Franco et al., 2008). The spawning season and location (river, bay, ocean) were also recorded as it indicates when species would be expected to be found together and their potential seasonal movements. Gregariousness or schooling type was defined as solitary, faculative or obligative which help to explain how species aggregate and how often. Fish that are obligative schoolers (highly gregarious) tend to be preferred prey of large and fast predators (Spitz, Ridoux & Brind’Amour, 2014). Mortality and maximum age are indicators of population turnover rates and longevity and may also be an indicator of population size (Palomares & Pauly, 1998). Age or length at maturity affects the resilience of a population, as species that mature younger are more resilient (Froese & Binohlan, 2000). Number of eggs or brood size is an indicator of fecundity (Clavel et al., 2013). Spawning frequency was recorded as singular (semalparous), batch or serial spawning and annual which can indicate stability of stocks between years, where species that spawn more often tend to have more stable populations (Longhurst, 2002). Fish that provide parental care or give birth to live young (viviparous) tend to give birth to fewer, larger offspring, often in more sheltered habitats such as estuaries.

Habitat traits are important in defining how a species uses their environment. As we focused on a small ecosystem the habitat variables of a given species must match the available habitat of that ecosystem. We included the minimum and maximum known depth of the species as TBGB is a relatively shallow bay (max depth 200 m). Knowing the vertical space that the species occupy informs us of potential intraspecific competition (Munday, Jones & Caley, 2001). We included the preferred temperature gradient (tropical, subtropical, temperate or deep) as temperature is an important indicator of how species use the ecosystem (Malavasi et al., 2004). Horizontal habitat (coastal, neritic or ocean) was used as another indicator of how species may group together in similar habitats.

Diet traits allow us to understand a species position within a food web. Diet can be recorded in a number of ways, but for our purposes we sought a simple classification of diet. Therefore we have two diet variables only; diet category (omnivore; invertivore, piscivore, herbivore and gelatinous invertebrate feeders) and trophic level (obtained from FishBase for consistency).

Data sources

Functional traits of species were sourced primarily from FishBase—a global information system on fishes (Froese & Pauly, 2017) and from ‘The fishes of New Zealand’—a comprehensive text with citations of all known fish species in New Zealand (Roberts, Stewart & Struthers, 2015). Additional trait data were obtained from a combination of published research and reports. When data was obtained from sources other than FishBase or Roberts, Stewart & Struthers (2015) the source is referenced. To obtain traits from FishBase we utilised the R package rfishbase (Boettiger, Lang & Wainwright, 2012).

Data preparation

Only 22 of the 40 recorded traits had less than 25% missing data and were retained for analysis. 25% was selected as the cutoff as the accuracy of imputed datasets is seriously degraded above 20–25% for small datasets (Clavel, Merceron & Escarguel, 2014). Distance matrix calculations require complete information, therefore we choose to impute the missing data in these 22 variables. Numerous methods exist for imputing data, and many of these have been examined for their precision in imputing continuous variables (Penone et al., 2014; Clavel, Merceron & Escarguel, 2014). What is unknown is how well these packages perform for nominal variables. To find the most accurate imputation method for nominal data we used three different approaches (all implemented in R packages): random forests implemented in missForest (Stekhoven & Bühlmann, 2012), multiple correspondence analysis (MCA) implemented in missMDA (Josse & Husson, 2016) and polytomous logistic regression implemented in MICE (Van Buuren & Groothuis-Oudshoorn, 2011) (described in Data S1). We also selected a simple imputation method using the mode value for each variable to serve as a baseline. In the mode replacement method, all missing values are replaced with the same value that is most frequently occurring. This method was used to compare against other imputation methods that use more information to inform the imputation (Taugourdeau et al., 2014). To test the accuracy of the different imputation methods we first selected all 13 variables from the database with complete information (Table 1). For each method, we ran a simulation in which data were randomly deleted and imputed 100 times. The probability of the method correctly imputing values were tested over a range of proportions of missing data ranging from 0.05 to 0.45, increasing in steps of 0.05. The final accuracy was calculated as the number of incorrect imputations divided by the number of possible imputations.

Four of the 22 trait variables were continuous and were discretized to turn them into categorical variables. It was a goal of the discretization process to maintain the underlying distribution of the data while creating similar number of categories in each variable (Teletchea et al., 2009). Each continuous variable was plotted on a histogram and bins were selected such that the distribution of the variable was maintained using four or five bins (Fig. S1). The final categories for each continuous variable and their values are reported in Table 1. The final trait matrix consisted of m = 22 traits and n = 116 fish species.

Distance matrices

Hierarchical clustering methods utilise distance matrices to make groups. A distance matrix in this context is a measure of pairwise similarities or dissimilarities between species (rows) based on their trait values (columns). There are a wide range of distance matrices and clustering methods available to cluster nominal data, and the combinations selected will influence the resulting groups. Having nominal data prevents us from using some measures, such as Euclidean distances, as they assume an inherent ordering within variables. For binary data, treating data (0 or 1) as continuous is a valid measure of difference, but for variables with more than two categories the various distances between values do not represent meaningful differences. Boriah, Chandola & Kumar (2008) evaluate 14 alternative measures of calculating distance matrices for nominal data and here we evaluate five: simple matching (SM—as in Gower’s distance), Eskin, Lin, inverse frequency of occurrence (IOF), and Goodall’s, available in the R package nomclust which are described in Data S2. The other measures available are derivatives of these measures and were not shown to improve performance in preliminary analyses. Briefly, the five distance matrices are described. The SM distance, which is the simplest approach to creating a distance matrix, awarding 1 to observations that are the same and 0 if not. This is the approach used for Gower’s similarity measure of nominal data (Gower, 1967). Eskin’s distance, which uses a SM criteria that gives more weight to mismatches on variables that have more categories (Eskin et al., 2002). The inverse occurrence frequency (IOF) distance has the same approach as Eskin but gives less weight to mismatches on variables that have more categories (Sparck Jones, 1972). This uses the absolute frequencies of observed categories. Goodall’s distance, which when comparing two observations of a given variable, takes into account relative frequencies of categories (Goodall, 1966). A similarity value is assigned based on the normalised similarity between the two observations, where the similarity value is higher if a category occurs infrequently. This method takes into account that individuals attributes occur stochastically and independently in a population. Lin’s distance is an information theoretic definition of similarity based on relative frequencies (Lin, 1998). Matches are given higher weightings when they occur infrequently, and conversely mismatches are given higher weightings when they occur frequently.

Clustering methods

As we do not know the number of functional groups in the ecosystem a priori, we used hierarchical clustering to visualise group association given our chosen distance metric. Hierarchical clustering first places all n objects in n separate single member clusters, and larger clusters are formed by sequentially joining first individual observations and then groups of observations until at last all observations are in a single group. The closeness of pairs of observations or groups of observations to another are determined by a measure of distance calculated in the preceding step. In linkage, all pairwise inter-cluster dissimilarities are calculated. The pair of clusters that are least dissimilar (that is, most similar) is identified and these two clusters are fused. Once observations or clusters are joined to a group they remain as a part of that cluster for the remainder of the analysis. There are a number of linkage methods that can be used for this type of data (Blashfield, 1976) and here we explore three methods available in the R package nomclust. To describe the linkage methods we use the following notation: D(A, B) is the distance between clusters A and B, which have sizes n_A and n_B respectively. In single linkage (minimising inter-cluster dissimilarity), the dissimilarity between two clusters is the smallest of all pairwise distances between the observations in the two clusters: (1)${\displaystyle D\left(A,B\right)=\text{min}\left[d\left(x,y\right):x\in A,y\in B\right].}$ In complete linkage (maximises inter-cluster dissimilarity), the dissimilarity between two clusters is the largest of all pairwise distances between the observations in the two clusters: (2)${\displaystyle D\left(A,B\right)=\text{max}\left[d\left(x,y\right):x\in A,y\in B\right].}$ In average linkage, the dissimilarity between two clusters is the average of all pairwise distances between observations in the two clusters: (3)${\displaystyle D\left(A,B\right)=\frac{1}{n_A{n}_B}\sum _{x\in A}\sum _{y\in B}d\left(x,y\right).}$

Selection of distance matrices, clustering methods and number of clusters

Evaluating clustering outputs can occur in two ways; external, where the resulting clusters are compared against known groupings (as in supervised learning), or internal evaluation, where some metric (there are many) is used to evaluate cluster separation and compactness. Since in our case the true groupings are unknown only internal evaluation is considered. To select the best distance matrix and clustering method for our data we utilised internal evaluation measures available from nomclust (Šulc & Řezankovà, 2015). The within-cluster entropy coefficient (WCE) is a measure of compactness which evaluates the variability of each cluster by calculating a measure of normalised entropy (the number of variables that have the same categories from each of the variables evaluated) (Šulc, 2016). WCE is measured from 0 to 1, where a lower value indicates intra-cluster homogeneity. Due to the way that these values are calculated they will generally always improve by adding clusters to the solution because the within cluster variability decreases: (4)${\displaystyle WCE\left(k\right)=\sum _{g=1}^k\frac{n_g}{n\times m}\sum _{c=1}^m\left(-\sum _{u=1}^{K_c}\left(\frac{n_{gcu}}{n_g}ln\frac{n_{gcu}}{n_g}\right)\right)}$ where n is the total number of objects (species), m is the number of variables (traits), n_g is the number of objects in the g^th cluster (g = 1, …, k) and n_gcu is the number of objects in the g^th cluster by the c^th variable with the u^th category (u = 1, …, K_c).

To select the number of groups we use the pseudo F coefficient based on the entropy (PSFE), a measure of separation (Šulc, 2016). The PSFE is a measure of entropy of the between- and within-cluster variability adjusted for the number of clusters and number of objects in the cluster where a higher value indicates a better grouping: (5)${\displaystyle PSFE\left(k\right)=\frac{\left(n-k\right)\left[nWCE\left(1\right)-nWCE\left(k\right)\right]}{\left(k-1\right)nWCE\left(k\right)}}$ where n is the number of observations and k is the number of clusters, nWCE(1) is the variability in the whole dataset, and nWCE(k) the within-cluster variability in the k-cluster solution.

Therefore, a more informative measure of performance is the degree of improvement with increasing number of clusters. Results from these measures are therefore presented as the difference between the k_th cluster and the k_th−1. Equivalent measures of all the aforementioned evaluation techniques are available in nomclust using the Gini coefficient instead of entropy and are provided in Fig. S2 as a reliability measure of our results.

We use t-Distributed Stochastic Neighbour Embedding (t-SNE) (Van Der Maaten, 2014) to construct a two-dimensional scatter plot in which each point represents a species. t-SNE minimises the distance between two distributions, one that was derived from a similarity matrix, and one that is derived from embedding the same matrix. To do so, a principal components analysis (PCA) is constructed from a dissimilarity matrix which allows species with similar trait profiles to be mapped in two-dimensions. These graphs provide a visual demonstration of similar species by the closeness of their points, and we use these graphs to evaluate our final group clustering. We evaluated a cluster-wise measure of cluster stability through a bootstrapping procedure available in the R package fpc (Hennig, 2013). The clusterboot function draws a sample of size N from the original data set, computes the clustering using partitioning around mediods (PAM), then calculates the maximum Jaccard coefficient between the most similar cluster in the bootstrapped data sets (Hennig, 2007). PAM is an agglomerative clustering approach that moves a pre-defined number of centres, here three, five, seven and nine, around a group of data to find the total minimum distance between the centres and the observations (Brock et al., 2008). This is repeated 100 times and an average Jaccard coefficient is found for each of the clusters which is representative of cluster stability. The more stability in clusters the less deviation evident in the Jaccard coefficient and as such the results are plotted as error bars. Finally, the adjusted Rand index (Hubert & Arabie, 1985) is used to compare the partitioning of groups for the different combinations of distance matrix and cluster algorithm for three, five, seven and nine cluster groups. The index can range from 0 (groups are completely random, may be negative if the index is less than the expected index) to 1 (each group contains the same observations).

The results of these analyses are discussed in terms of connectedness, compactness, separation and stability. Compact groups are those which minimise the spread of observations within a cluster and are assessed with the WCE and visually with the t-SNE graphs. Separation refers to the between cluster distances, which ideally should be maximised and are assessed by the PSFE, and visually with the t-SNE graphs. A well connected cluster is one where an observations nearest neighbour is from its own cluster and is assessed solely by visualisation. Stability is assessed via the results of the bootstrapping procedure.