Performance of maximum likelihood mixture models to estimate nursery habitat contributions to fish stocks: a case study on sea bream Sparus aurata

View article
Note that a Preprint of this article also exists, first published May 3, 2016.


Evaluating the contribution of different sources to a mixture is a common problem in ecology, biology and natural resource management (Kimura & Chikuni, 1987; Smouse, Waples & Tworek, 1990; Van Dongen, Lens & Molemberghs, 1999; Fleischman & Burwen, 2003; Manel, Gaggiotti & Waples, 2005; Phillips, Newsome & Gregg, 2005; Newman & Leicht, 2007). Fish ecologists and fisheries scientists, for example, are frequently interested in estimating the contribution from different nursery habitats (sources) to adult aggregations, demographic units or stocks (mixtures). Beyond its inherent scientific interest, this kind of assessment has practical relevance for both management and conservation purposes (Kerr, Cadrin & Secor, 2010). Assessing the accuracy and precision of parameters resulting from mixture analysis is a fundamental, still often neglected step, required to facilitate the incorporation of these results into modern management models (Kritzer & Liu, 2014).

Mixture analysis in fish ecology and other disciplines relies heavily on the use of artificial and natural tags suitable for tracking or identifying the different sources (origins) contributing to a mixture (Gillanders, 2009). Within natural tags, the elemental and isotopic composition of teleost fish otoliths has been an increasingly common choice for this type of studies during the last decades (Kerr & Campana, 2014). Otoliths grow throughout lifetime by a regular deposition of calcium carbonate and protein layers, which, unlike bones, are not reabsorbed (Panfili et al., 2002). While calcium can be partially replaced by other metals (including Sr, Mn and Ba), dominant carbon and oxygen isotopes (12C and 16O) can be replaced by their less frequent alternatives 13C and 18O. When these substitutions are under weak internal control, they may reflect environmental and/or physiological variability (Panfili et al., 2002), and the elemental/isotopic otolith signatures can be considered “fingerprints” for the water masses inhabited by fish at carbonate deposition time (Elsdon et al., 2008). As deposition time can be often inferred from the same otolith through ageing techniques, a retrospective identification of nursery or feeding habitats, demographic units (∼stocks) and/or migration patterns becomes possible (Campana & Thorrold, 2001; Rooker & Secor, 2004; Elsdon et al., 2008; Rachel et al., 2008; Arkhipkin, Schuchert & Danyushevsky, 2009; Darnaude et al., 2014; Niklitschek et al., 2014).

Two main statistical approaches have been used to estimate the contribution of different sources to a mixture: Discriminant Functions (DF) and Mixture Models (MM) (Millar, 1990a; Koljonen, Pella & Masuda, 2005). DF approaches include linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), multinomial regression (MNR) and random forest analysis (RM), among several other related methods (Edmonds, Caputi & Morita, 1991; Elsdon & Gillanders, 2003; Pella & Masuda, 2005; Mercier et al., 2011; Jones, Palmer & Schaffler, 2016). Although some parametric DF can be seen as special cases of MM, they have some important differences in focus and estimation procedures. DF focus on developing discriminant algorithms, which are fit (“trained”) using samples from known origins (i.e. pre-migratory juveniles sampled at their nursery-sources), and applied, at a subsequent step, to assign putative origins to older (adult) individuals sampled from the mixed-stock. Therefore, mixing proportions are not estimated directly as model parameters, but quantities computed afterwards from the putative origins assigned by the model to the individuals present in the mixed-stock dataset. MM approaches focus, instead, on estimating these mixing proportions, which are explicit and fundamental model parameters, estimated directly from the mixed-stock dataset. Baseline nursery-signatures are also explicit MM parameters, which are commonly of high scientific interest on their own.

Described in detail by Everitt & Hand (1981), MM were probably introduced into fisheries science by Cassie (1954). Applications to mixed stock analysis were first presented by Fournier et al. (1984) and increased largely after the HISEA software was made available by Millar (1990b). Under their finite mixture distribution form (Everitt & Hand, 1981), MM are defined as, f(x)=k=1Kpkg(xi;θk) where, the density function f(x) is defined by three groups of parameters: the number of components or sources (K), the mixing proportions (pk) and the set of baseline parameters θk that characterize each source k, given the probability distribution function g(). This function is frequently, although not necessarily, assumed multivariate normal. Thus, θk can be decomposed in a vector of means (μk) and a covariance matrix (Σk) for all response variables used to characterize each source k. Translating this terms into the lexicons of otolith chemistry and mixed-stock analysis, K corresponds to the number of nursery or spawning sources, pk to the proportional contribution made by each of these sources to the mixed stock, and θk to the baseline chemical signatures (i.e. means and covariances of elemental or isotopic ratios) that characterize otolith material formed at each nursery source k.

Most MM applications to fisheries during the last four decades have used maximum likelihood methods (Millar, 1987; Reynolds & Templin, 2004). Within this framework, the expectation-maximization algorithm (EM) (Dempster, Laird & Rubin, 1977) has been used as the dominant likelihood maximization procedure. In more recent years, and following an evident worldwide trend in statistical methods, an important development of Bayesian approaches has been reflected in an increasing number of mixed stock applications (Pella & Masuda, 2001; Koljonen, Pella & Masuda, 2005; Munch & Clarke, 2008; White et al., 2008; Smith & Campana, 2010; Standish, White & Warner, 2011), including parametric and non-parametric approaches and important software development efforts (Neubauer, Shima & Swearer, 2013). Despite of these promising developments, MM methods probably remain as the most common approach being used for stock mixture analysis at scientific and management organizations.

Most MM applications to mixed stock analysis tend to focus on estimating pk, given all potential nursery sources have been previously identified (i.e. K is known) and sampled for pre-migratory juveniles to produce ex-ante θk estimates, which are then supplied to the MM as fixed quantities. This conditional MM approach follows Millar (1987) and tends to be dominant in the MM literature (Hamer, Jenkins & Gillanders, 2005; Crook & Gillanders, 2006; Schloesser et al., 2010; Secor, Gahagan & Rooker, 2012). Some drawbacks of this approach, shared by DF methods, are that it fails if not all nursery-sources are known or sampled, and that it neglects all information about θk being contained in the mixed-stock data. Under an unconditional MM approach, θ^k are produced or updated using the information contained in the mixed-stock data. Thus, unconditional models benefit (greatly) from nursery-source sampling, but can still be fit if such sampling fails for some or all nursery-sources. Moreover, these models are expected to be less sensitive to small sampling sizes (Koljonen, Pella & Masuda, 2005).

Failing to sample some or all nursery-sources is a common problem in fish ecology, particularly for open marine populations (Campana et al., 2000), where the location of nursery habitats may be unknown, remote or inaccessible, or where juvenile fish may be to cryptic or invulnerable to the sampling gear. Under these scenarios, there may not be other option than the simultaneous (unconditional) estimation of both pk and θk parameters from the mixed-data (Smouse, Waples & Tworek, 1990; Niklitschek et al., 2010; Smith & Campana, 2010). Furthermore, if not even K (the total number of nursery-sources) is known, all three sets of parameters (pk, θk and K) need to be estimated from the mixed-stock data. Estimating all three sets of parameters within the same likelihood maximization fit may lead however to identifiability issues. As an alternative, a model comparison approach can be used (Everitt & Hand, 1981), to select the most informative value of K according to some criterion such as Akaike (1973) and Schwarz (1978), entropy (Celeux & Soromenho, 1996), deviance (White et al., 2008) or some other information criterion, depending on the modelling framework being used.

The simultaneous estimations of pk, θk and/or K from the mixed-stock data might introduce bias related to the existence of multiple solutions and/or to multiple local maxima (McLachlan & Peel, 2004; Reynolds & Templin, 2004), as well as to the presence of confounding covariates affecting the mixed-stock data, which may either blur or spuriously enhance variability in nursery signatures. Among the three sets of MM parameters, K seems to be the most prone to bias, as shown by several theoretical and practical studies (Titterington, 1990; Celeux & Soromenho, 1996; White et al., 2008; Neubauer, Shima & Swearer, 2013). Assessing the magnitude of such bias under incomplete sampling scenarios is not frequently reported (Wood et al., 1987) as no reference data exists to contrast the parameters estimated by the model. Indirect assessment approaches can be conducted, however, using either simulated or empirical datasets whose true pk, θk and/or K parameters were actually known.

In this article, we evaluate the performance of maximum likelihood mixture models, from now on “ML-MM,” to estimate pk, θk and K parameters under several scenarios that simulated incomplete sampling and different degrees of separation among nursery signatures. Departing from the mainstream of ML-MM, we adopted an unconditional approach to estimate and/or update θ^k using all available nursery-source and mixed-stock data. To conduct this evaluation, we follow a case study approach focused on a comprehensive spatio-temporal dataset containing baseline chemical signatures from young-of-the-year Sparus aurata collected in four separate nursery habitats (Mediterranean lagoons), in three highly contrasting years (Tournois et al., 2013). By sub-setting, resampling and manipulating this dataset we evaluated bias and uncertainty in pk, θk and K as a function of (i) the number of nursery sources identified and sampled for pre-migratory juveniles to estimate nursery-signature parameters, (ii) the observed variability in nursery-signatures among cohorts, and (iii) the degree of separation among nursery-signature centroids.

Materials and Methods

Data set description

The dataset used in the present work, obtained from Tournois et al. (2013), included otolith fingerprints from 301 young-of-the-year Sparus aurata, collected in the Gulf of Lions (NE Mediterranean Sea) from four Mediterranean lagoons: Bages-Sigean, Mauguio, Salses-Leucate and Thau, in three different years (= cohorts): 2008, 2010, and 2011. Collection always occurred in the late summer, just before they returned to mix with individuals from nearby lagoons in the open sea.

The chemical composition analysis of otolith samples was performed using Solution Based Inductively Coupled Plasma Mass Spectrometry, and included 43Ca and another 11 elements (Tournois et al., 2013). For this study, we only kept the seven most discriminant ones: 7Li, 11B, 25Mg, 85Rb, 86Sr, 89Y and 138Ba. Their concentrations in the otoliths were expressed as elemental ratios to Ca, and standardized to mean = 0, and SD = 1 to scale all elements equally and facilitate bias analysis. Three obvious outliers were discarded, working with a depurated sample size of 298 otoliths. Data was normalized using a multivariate Box & Cox (1964)’s transformation although it failed to fully normalize three of the seven elemental ratios (Mg, Rb and Ba).

The four lagoons sampled by Tournois et al. (2013) differ greatly in size, depth, freshwater input and degree of connection with the sea, leading to physical and chemical differences in the water and, therefore, in otolith signatures of juvenile S. aurata (Tournois et al., 2013). Nonetheless, these lagoons are strongly influenced by rainfall, wind and other environmental forces (Sara, Leonardi & Mazzola, 1999; Martins et al., 2001), leading to high interanual variability in otolith signatures (descriptive statistics provided in Appendix 2). For instance, the average squared Mahalanobis distance among nursery sources decreased from 3.29 to 1.18 SD, between 2008 and 2011, while its multivariate spread (effective general variance) ranged between 0.54–0.86 within the same interval (Table 1). A multivariate analysis of variance showed significant effects of cohort, nursery-source and their interaction upon these nursery signatures (p < 0.001).

Table 1:
Average Mahalanobis distance and effective standard deviation for elemental compositions of selected metals in otoliths of juvenile Sparus aurata.
Average Mahalanobis distance and effective standard deviation for elemental compositions of selected metals in otoliths of juvenile Sparus aurata. Average distances within nursery source and cohort computed from vectors of observations. Average distances within years and within sources computed from vectors of means corresponding to each source or year, respectively. All values computed after standardizing all data to = 0 and = 1.
Cohort Bages-Sigean Mauguio Salses-Leucate Thau Within cohorts
Average squared Mahalanobis distance (Δ2)
2008 2.52 2.53 2.53 2.55 3.29
2010 2.47 2.52 2.47 2.49 2.78
2011 2.48 2.54 2.51 2.50 1.18
Within sources 3.49 3.21 2.20 1.88 3.50
Effective standard deviation |S|1/14
2008 0.40 0.39 0.35 0.40 0.65
2010 0.41 0.51 0.63 0.53 0.86
2011 0.38 0.39 0.46 0.49 0.54
Within sources 0.63 0.66 0.73 0.71 0.85
DOI: 10.7717/peerj.2415/table-1

Simulation approach and scenarios

All simulations and analyses were based upon the construction of two datasets: (i) a “nursery-source dataset” that represented otolith composition data from pre-migratory juveniles, sampled at their corresponding nursery-origins, and (ii) a “mixed-stock dataset” that represented otolith composition data from older fish sampled after they had mixed with fish from the other four origins (i.e. at the open sea). Besides the observed variability among the three sampled cohorts, we simulated additional variability in two dimensions: (i) the number of sources being sampled and included in the “nursery-source dataset” and/or (ii) the degree of separation among nursery-signature centroids (Table 2).

Table 2:
Main configuration of the simulation and resampling procedures used for assessing the performance of maximum likelihood mixed models.
Main configuration of the simulation and resampling procedures used for assessing the performance of maximum likelihood mixed models. Observed cohorts corresponded to juvenile Sparus aurata collected from four Mediterranean lagoons in three highly contrasting years. Virtual cohorts corresponded to artificial data, aimed to expand the observed range of separation among nursery signatures k and y subscripts represent nursey-sources and cohorts, respectively. Observed mean vectors and covariance matrices available in Appendix 2.
Observed cohorts Virtual cohorts
Number of nursery sources included in nursery-source datasets (KS) 0–4 0–4
Number of cohorts 3 6
Nursery-signature mean vectors (μk,y) Observed μk,y μk,2010 scaled to match target separation
Nursery-signature covariance matrices Σk,y Observed Σk,y Observed Σk,2010
Average Mahalanobis distance among nursery-signature centroids (Σy2) Observed Δy2=1.183.29 Simulated Δy2 = {0.5, 1.5, 2.5, 3.5, 4.5, 5.5}
Mixing proportion of nursery-sources in the mixed-stock dataset pk,y pk,y = {0.1, 0,2, 0.3, 0.4}, randomly assigned to sources within runs pk,y = {0.1, 0,2, 0.3, 0.4}, randomly assigned to sources within runs
DOI: 10.7717/peerj.2415/table-2

Nursery-source sampling scenarios

Five scenarios were defined by the number of nursery-sources included in the nursery-source dataset. At each run, KS = 0–4 sources were randomly selected to be included in the nursery-source dataset as “known” nursery habitats, which had been sampled for pre-migratory juveniles. These data had been used as baseline data to produce initial estimates for θ^k and then mixed with the mixed-stock data to produce final θ^k parameters, following an unconditional likelihood maximization approach. All remaining “unknown” nursery sources (KU = K−KS) were excluded from the nursery-source dataset and lacked of initial θ^k values.

Separation among nursery signatures

To improve our empirical understanding about the effects the separation among nursery-signature centroids may have on bias and uncertainty in ML-MM parameters, we applied the five sampling scenarios to (i) the three observed cohorts (2008, 2010 and 2011); and (ii) six virtual cohorts created to expand the range of Mahalanobis distances Δ2 among nursery-signature centroids observed in the three sampled cohorts (Δ2 = 1.18–3.29), Table 1. To build these six virtual cohorts, covariance matrices were set equal to those estimated in 2010 for each nursery-source (Appendix 1), while the vectors of means correspondng to this same year (Appendix 1) were re-scaled to get Δ2 values of 0.5, 1.5, 2.5, 3.5, 4.5 and 5.5 (Table 2).

Resampling procedures

Datasets for each tested scenarios and independent run were produced by parametric bootstrapping. Nursery-source datasets included 25 observations drawn from each of the KS known nursery-sources and from each of the three cohorts. Mixed-stock datasets included a total of 100 observations per cohort, drawn from all four nursery-sources, mixed using uneven mixing proportions (m) of 0.1, 0.2, 0.3 or 0.4. These four proportions represented the true value of pk and were randomly allocated to the four nursery-sources, within each resampling run. Resampling was followed by a standard modelling and fitting procedure, detailed in Appendix 1. This was a 12-steps sequence, which was repeated 1,000 times for each cohort and scenario. Each repetition was labelled as a “resampling run.”

Mixing proportions

To evaluate bias and uncertainty of mixing proportion estimates (p^k), we assumed the true number of nursery sources was known and fixed (K = 4) across all scenarios. Bias in p^ was computed as the sum of the average differences between the estimated (p^mr) and the true mixing proportion (pmr) assigned to each nursery-source within each resampling run. The subscript m = {0.1, 0.2, 0.3, 0.4} represents here the vector of mixing proportions randomly allocated to the four nursery-sources, within each of the R = 1,000 resampling runs. Therefore, BIp^ = m=1M1Rr=1R| p^mrpmr |

Uncertainty in p^k was indexed as the average of the four empirical standard errors of p^ computed for each of the four possible values of m, across the R = 1,000 resampling runs, SEp^ = 1Mm=1M(p^mrpmr)2R

Nursery-signature parameters

Under the assumption of multivariate normal distribution, each set of estimated nursery-signature parameters θ^k was composed by a vector of means μ^k and a covariance matrix Σ^k, which described the multivariate distribution of the seven chemical elements measured in the otoliths included in the dataset. Assessing bias in Σ^k is a complex task, which we considered that exceeded the scope of this paper. Therefore, all bias measures provided hereafter for θ^k are strictly referred to μ^k.

As done for p^k, the assessment of bias and uncertainty in θ^ was conducted assuming a known and fixed value of K = 4. Overall bias in θ^k was indexed by averaging the absolute mean differences between estimated (μ^kr) and true (μkr) vectors of means, across all elemental ratios (J = 7) and nursery-sources (K = 4). As all elemental ratios were previously standardized, bias units were equivalent to standard deviations and computed as follows, BIθ^ = 1JKj=1Jk=1K| 1Rr=1Rμ^krμk |

Overall uncertainty in θ^i, was indexed by its effective standard deviation (Peña & Rodríguez, 2003), defined as, SEθ^=1Kk=1K| ^k |1/2J where, Σ^k is the covariance matrix computed from all μ^kr estimated for each nursery-source k across the R = 1,000 resampling runs.

Number of contributing nursery sources

For assessing bias and uncertainty in K^, this parameter was not set to a constant as before, but estimated by a model selection approach as described in the next section. Bias in K^ was computed as BIK^ = K^4, and uncertainty (SEK^) as the standard error of K^ computed from all resampling runs, within each cohort and scenario. The strength of the selection was indexed by ΔBICK, defined as the average of the differences between the lowest and the median Bayesian Information Criterion (BIC) values computed for all competing models, within each resampling run.

ML-MM parameter estimation

Parameters p^k and θ^k (μ^k and Σ^k) were estimated by maximum likelihood, using the Expectation-Maximization (EM) algorithm (Dempster, Laird & Rubin, 1977), modified to follow an unconditional approach where initial θ^k estimates were updated to maximize the joint likelihood of both nursery-source and mixed-stock datasets (Appendix 1). The M-step was constrained to produce definite positive covariance matrices, with det(Σ) > 109. Starting μ^k and Σ^k parameters for the KS known nursery-sources were computed directly from the nursery-source dataset. Starting μ^k for the KU unknown nursery-sources were obtained from the mixed-stock dataset trough a semi-supervised K-means clustering procedure, implemented using the R-package “vegclust” (De Cáceres, Font & Oliva, 2010), which allowed for combining “fixed” and “mobile” centroids. Fixed centroids corresponded to μ^k estimated at the previous step for the KS known nursery-sources. KU additional mobile centroids, which represented the KU unknown nursery sources, were estimated as the values that minimized the mean square distance between the mixed-stock data and all (fixed and mobile) model centroids. Starting Σ^k for the KU unknown nursery-sources were computed, at a subsequent step, from the mixed-stock data clustered into each of these KU additional clusters (Appendix 1). Starting p^k were calculated as the empirical proportion of individuals represented in the mixed-stock dataset assigned to each putative nursery-source k in order to maximize de probability density of each observation p(xi|θ^k), assuming xiMVN(θ^k).

Parameter K^ was estimated following a model selection procedure, where multiple competing models were fit to each pair of nursery-source and mixed-stock datasets. These competing models considered a range of K values, between a minimum of Kmin = KS and a maximum defined arbitrarily as Kmax = 8. As result, within each resampling run, and depending upon the value of KS, a total of 4–9 competing models were fit and compared. Model comparisons were performed using Schwarz (1978)’s BIC, where the most informative value of K was addressed as the “estimated” number or nursery-sources (K^).


Mixing proportions

Bias in p^(BIp^) ranged between 0 and 0.24 across all data availability scenarios and observed cohorts. Relatively unbiased p^ estimates (BIp^ < 0.1) were obtained under most data availability scenarios (KS = 1–4) for cohorts 2008 and 2010 (Table 3; Fig. 1), but exceeded 0.11 across all scenarios for cohort 2011. The highest values of BIp^ (0.12–0.24) were found at KS = 0, when none of the nursery-sources were known (Table 3; Fig. 1). When all nursery-sources were known (KS = 4), BIp^ approached zero for cohorts 2008 and 2010, but remained relatively high (BIp^∼0.11) for cohort 2011. Such a decrease in bias was near one order of magnitude for cohorts 2008 and 2010, and greater than 50% for cohort 2011 (Fig. 1). Uncertainty in p^(SEp^) ranged between 0.06 and 0.28, with much higher values observed for cohort 2011 (SEp^ = 0.03–0.11) than for both cohort 2010 (SEp^ = 0.007–0.05) and for cohort 2008 (SEp^ = 0.01–0.07). Following a pattern somewhat similar to BIp^, we found that SEp^ decreased rapidly as KS increased (Fig. 1).

Table 3:
True and estimated mixing proportions of nursery-sources in the mixed-stock dataset (pk).
True and estimated mixing proportions of nursery-sources in the mixed-stock dataset (pk). Data from all nursery-sources combined for simplicity. True values corresponded to the proportion of bootstrap samples drawn from each cohort and nursery-source at each resampling runs (R = 1,000). These nursery-source proportions were assigned randomly from the vector m = {0.1, 0.2, 0.3, 0.4}.
Cohort True proportion in mixed-stock dataset (pk) Sampling scenario (number of habitats represented in nursery-source datasets)
KS = 4 KS = 3 KS = 2 KS = 1 KS = 0
2008 0.1 0.10 0.10 0.11 0.13 0.15
0.2 0.20 0.20 0.20 0.21 0.21
0.3 0.30 0.30 0.30 0.29 0.29
0.4 0.40 0.40 0.39 0.37 0.35
2010 0.1 0.11 0.11 0.12 0.14 0.16
0.2 0.20 0.20 0.21 0.21 0.22
0.3 0.30 0.30 0.30 0.29 0.28
0.4 0.40 0.39 0.38 0.36 0.35
2011 0.1 0.14 0.15 0.16 0.17 0.19
0.2 0.21 0.22 0.22 0.23 0.23
0.3 0.28 0.29 0.28 0.28 0.27
0.4 0.36 0.34 0.34 0.33 0.31
DOI: 10.7717/peerj.2415/table-3
Bias and uncertainty in mixing proportions for observed cohorts.

Figure 1: Bias and uncertainty in mixing proportions for observed cohorts.

Bias (A) and uncertainty (B) in mixing proportions (p^) of four nursery-sources to artificial mixed-stocks of Sparus aurata built using data from cohorts 2008, 2010 and 2011. Different sampling scenarios are represented by the number of nursery-sources (KS) simulated to be known and sampled for pre-migratory juveniles.

The rank order of BIp^ and SEp^ values among the three observed cohorts was inverse to the rank order of their average distance among nursery-signature centroids (Table 1). This inverse relationship was also observed in the six nursery-signature separation scenarios, where BIp^ and SEp^ decreased rapidly as the distance among nursery signatures increased (Fig. 2). Such a decrease tended to evolve from a linear pattern at the worse 1–2 scenarios (KS = 0–1) to a more exponential decay pattern as KS approached its maximum (Fig. 2) There was also an evident trend to observe positive bias at lower p^ values, and negative bias at higher p^ values, which was more pronounced as KS decreased (Table 3).

Nursery-signature parameters

Estimated nursery-signature parameters provided relatively unbiased and consistent (similar shape and orientation) fits to the “true” distribution of means, even at the KS = 0 scenario (Fig. 3). Considering the observed cohorts, BIθ^ ranged between 0.005 and 0.13 at KS = 4 and KS = 0, respectively (Fig. 4). As observed before for mixing proportions, BIθ^ tended to be much higher for cohort 2011 than for cohort 2008, across all scenarios. However, the values of BIθ^ for cohort 2010 tended to be much closer to those computed for cohort 2011, than to the ones computed for cohort 2008 (Fig. 4). As for the six nursery-signature separation scenarios (Fig. 2), BIθ^ tended to increase with distance from minimum values at Σ2 = 0.5 towards maximum values and highest sensitivity to incomplete sampling Δ2 values of 3.5 or 4.5 depending on the number of known nursery sources (Fig. 2).

Bias and uncertainty in mixture model parameters for virtual cohorts of Sparus aurata.

Figure 2: Bias and uncertainty in mixture model parameters for virtual cohorts of Sparus aurata.

Average bias and uncertainty in mixing proportions p^ (A and B), nursery-signature estimates θ^ (C and D) and number of contributing sources K^ (E and F) as a function of the simulated distance among nursery-signature signature centroids and the number of nursery-sources (KS) simulated to be known and sampled for pre-migratory juveniles.
Empirical and estimated distributions of elemental ratios in juvenile Sparus aurata otoliths.

Figure 3: Empirical and estimated distributions of elemental ratios in juvenile Sparus aurata otoliths.

Principal component diagrams representing the distribution of otolith elemental ratios in juvenile Sparus aurata. Coloured points represent empirical means per nursery source and cohort (2008, 2010 and 2011), corresponding to 1,000 bootstrap samples (size = 25). Ellipses represent 67% confidence intervals for estimated means under two selected sampling scenarios: KS = 0, where no nursery-sources were sampled for pre-migratory juveniles (2008-A, 2010-A and 2011-A) and KS = 3, where three out of the four nursery-sources were simulated as known and sampled for pre-migratory juveniles (2008-A, 2010-A and 2011-A).
Bias and uncertainty in nursery-signature estimates.

Figure 4: Bias and uncertainty in nursery-signature estimates.

Average bias (A) and uncertainty (B) in estimated nursery signatures corresponding to juvenile Sparus aurata sampled from four nursery-sources in three climatically contrasting years, bootstrapped and combined into artificial mixed-stocks. Different sampling scenarios were represented by the number of nursery-sources (KS) simulated to be known and sampled for pre-migratory juveniles. Data obtained after 1,000 resampling runs per scenario.

Uncertainty in θ^k(SEθ^) was less sensitive to data availability (KS) than BIθ^, showing a moderate although nearly constant decrease with KS, in cohorts 2008 and 2010. For cohort 2011, however, SEθ^ was similarly high between KS = 0 and KS = 2, decreasing afterwards. Uncertainty in θ^k tended to be higher for cohorts 2010 and 2011 than for cohort 2008, across all scenarios, but particularly between KS = 1 and KS = 3. Results from simulated nursery-signature separation scenarios (Fig. 2) showed SEθ^ was not affected by distance among nursery signatures when all nursery-sources were previously known and sampled (KS = 4). Otherwise, SEθ^ tended to decrease with distance among nursery signatures, although for KS = 0 this decrease was only evident at the two highest distances among nursery signatures (Fig. 2).

Number of contributing nursery sources

ML-MM tended to underestimate the true value of K for all observed cohorts, under most data availability scenarios (Fig. 5), with negative bias (BIk^), between −0.56 and −2, across all incomplete sampling scenarios. Only under the ideal scenario (KS = 4) BIk^ became zero and the true value of K was correctly estimated in 100% of all simulations. Nonetheless, it must be recalled that under KS = 4, K^ was constrained to values greater or equal to Kmin = 4. As for the remaining scenarios (KS < 4), the absolute value of BIk^ tended to decrease as KS increased, at least for cohorts 2008 and 2010 (Fig. 5).

Estimated number of contributing sources.

Figure 5: Estimated number of contributing sources.

Estimated number of nursery-sources (K^) contributing to simulated mixed-stock of Sparus aurata obtained following a model selection approach based on Bayesian Selection Criterion. Different sampling scenarios were represented by the number of nursery-sources (KS) simulated to be known and sampled for pre-migratory juveniles. Relative frequencies computed after 1,000 resampling runs per tested scenario. BIk^ = mean bias for K^; SEK^ = mean standard error of K^; ΔBICK = difference between minimum and median values of the Bayesian Information Criterion.

The under-estimation of K^ was a highly consistent pattern, observed across all nursery-signature separation scenarios, were no estimated values of K^ ≥ 5 were ever obtained, regardless of the data availability scenario or cohort (Fig. 2). Nonetheless, the magnitude of this bias decreased by 50–70% as the distance among nursery signatures increased (Fig. 2). Uncertainty in K^(SEK^) was zero at KS = 4, but relatively high (SEK^ > 0.4) in cohorts 2008 and 2010 for all KS < 4, reaching maximum values at KS = 2 (Fig. 5). The strength of the model selection measured by ΔBICK tended to increase towards KS = 4, where it was maximum for cohorts 2008 and 2010 (Fig. 5). The consistent, although biased, values of K^(BIK^ = 1–2) lead to particularly low values of SEK^ for cohort 2011, across all scenarios. Consistent values of SEK^ ≈ 0 at KS = 4, and a weak relationship between KS and SEK^ for all KS < 4 were observed across all nursery-signature separation scenarios. Nonetheless, it was evident that SEK^ tended to increase with distance among nursery signatures for all incomplete sampling scenarios.


In this article, we combined real-world and virtual datasets to evaluate the performance of unconditional ML-MM when used to estimate the three fundamental sets of mixed stock parameters (pk, θk and K), under a range of data availability and distance among nursery-signature scenarios. Although using a single real-world dataset might limit the generalization of our results, which may not be transferable to other stocks, we believe the large variability observed among cohorts in the real-world dataset, along with the additional variability included in the simulated datasets might represent a relevant part of the variability that could be found in other populations and geographical areas.

Mixing proportions estimated by ML-MM (p^k) showed low bias and variability when at least one nursery source was included in the nursery-source dataset. Large variability in bias and uncertainty of p^k was found, however, among cohorts, with cohort 2011 exhibiting the highest bias and uncertainty values in p^k. This variability in bias and uncertainty among cohorts was consistent with inter-annual differences in sensitivity to incomplete sampling resulting from variable degrees of separation among nursery signatures. Results from simulations showed, for instance, that a minimum of three sampled nursery sources had been required to reduce BIp^ below 0.1, given a Mahalanobis distance among nursery-signature centroids ≤ 1.5.

Overall our results confirmed the suitability of using ML-MM for estimating unbiased mixing proportions, given the number of nursery-sources (K) was known and there was a proper balance between the number of nursery-sources sampled for prior nursery-signature parameters and the actual degree of separation among nursery-signature centroids. Although results from our simulated separation scenarios cannot be turn into a prescriptive guideline, they explore such trade-offs and provide a general idea about which combinations might work, at least for Mediterranean stocks of S. aurata.

When incomplete sampling and/or reduced distance among nursery signatures biased mixing proportions estimates, bias tended to be greater for the most extreme p^k values, which were shifted towards intermediate values. Therefore, the smallest nursery contributions tended to be overestimated, while the largest ones to be underestimated. This behaviour is probably related to the model constrain that all proportions must sum one, which limits the parameter space and tend to increase negative correlation between the two most extreme values. It might be also somewhat related to the EM algorithm, which may converge to unsatisfactory local maxima (Marin, Mengersen & Robert, 2005).

Nursery-signature parameter estimates were relatively unbiased for all evaluated cohorts, under most data availability scenarios. When two or more nursery-sources were known and previously sampled, BIθ^ dropped below 0.10. Uncertainty in θ^k, however, was relatively high under incomplete sampling scenarios, particularly in cohorts 2010 and 2011, where SEθ^ exceeded 0.20 at all KS < 3 scenarios. The important differences in bias and uncertainty we found among cohorts, were likely related to large environmental variability in the study area (Tournois et al., 2013), which were reflected in highly variable distribution patterns of nursery signatures among cohorts (Fig. 3). Simulated nursery-signature separation scenarios showed a dome-shaped relationship between BIθ^ and distance among nursery-signature centroids, with maximum bias at intermediate distances and minimum bias at the shortest distance we tested (Δ2 = 0.5). This counter-intuitive pattern may have emerged from the antagonistic effects of less separable nursery signatures but more constrained parameter spaces at shorter distances among their centroids.

Unlike what we found for p^k and θ^k, the performance of our model selection approach for estimating the true number of nursery-sources (K) was poor and exhibited an evident trend to underestimate the true value of K by 1 or 2 nursery-sources in all observed and most of the simulated cohorts. As the magnitude of this underestimation bias was constrained by the number of known nursery-sources, it resulted obvious that enhancing our knowledge about the minimum number of contributing sources had reduced the risk of underestimating K. Given this biasing trend, the model selection approach we followed provided a lower bound rather than an accurate estimate of K.

Exploratory comparisons (not shown in results) between BIC and Akaike (1973)’s Information Criteria (AIC) yielded results that agreed with the idea that BIC produces underestimated but less variable estimates of K, while the opposite would be true for AIC (Koehler & Murphree, 1988). Alternative model selection criteria have been proposed by several authors and deserve further testing under a greater range of scenarios (Celeux & Soromenho, 1996). Departing from model selection approaches, bootstrapping has been proposed and used to evaluate consistency in K^ (McLachlan & Peel, 2004). Neubauer, Shima & Swearer (2013) followed, instead, a Bayesian approach to internalize the estimation of K into a Dirichlet process mixture model, which had the advantage of producing marginal distributions over a range of plausible K^ values and allowed for a direct probabilistic interpretation of model results. There is, however, a generalized view that estimating K is one of the most difficult tasks within MM, for which more satisfactory solutions are still needed (Celeux & Soromenho, 1996; McLachlan & Peel, 2004; White et al., 2008; Neubauer, Shima & Swearer, 2013).

Although we focused on testing the performance of unconditional ML-MM, we must acknowledge the growing importance of Bayesian approaches, observed in relatively recent years (Pella & Masuda, 2001; Marin, Mengersen & Robert, 2005; Munch & Clarke, 2008; White et al., 2008; Smith & Campana, 2010; Standish, White & Warner, 2011; Neubauer, Shima & Swearer, 2013). They represent obvious alternatives to conventional ML-MM which may be particularly advantageous for considering missing or incomplete data scenarios. Other approaches used for estimating mixing proportions under partial or complete lack of nursery-source data have been based upon unsupervised clustering followed by discriminant analysis (Arkhipkin, Schuchert & Danyushevsky, 2009; Shima & Swearer, 2009; Schuchert, Arkhipkin & Koenig, 2010). Nonetheless, no independent assessments of bias or uncertainty seem to be available for this clustering approaches when applied to mixed stock analysis under incomplete sampling scenarios.

The large differences in chemical nursery signatures we found among cohorts in this dataset reflected large interanual variability in these nursery habitats (Tournois et al., 2013), which may be common to most shallow water and estuarine nursery areas (Secor, 2015). The high sensibility of ML-MM reliability to this inter-annual variability highlights the need to assure true independence among individual samples being used to build nursery-signature parameter estimates. Otherwise, variability among nursery-sources might be easily confounded with variability among years, schools, sampling events or other sources of correlation, commonly neglected in fisheries and ecological studies (Zuur, Ieno & Smith, 2007). Adding the effects of these random sources of correlation through a mixed or hierarchical approach (Bolker et al., 2009) has been already explored as extension of the EM algorithm by Heinzl & Tutz (2013). A very intuitive step here would be to combine data from multiple cohorts to improve the estimation of K, although this task could be also achieved using the joint likelihood from individual ML-MM models fit to each yearly dataset.

It is important to acknowledge that we did not explore sample size effects, but used constant sample sizes of 25 individuals per known nursery-source and 100 individuals per mixed-stock and cohort. While no major changes in bias would be expected at different sample sizes, an obvious reduction in uncertainty would be expected if the number of fish included in the nursery-source and/or mixed-stock datasets were larger. It would be of particular interest to consider intermediate scenarios where some small level of sampling existed for some/all nursery-sources, which maybe a common real-life situation when juvenile sampling is rather opportunistic. Under the unconditional approach we have followed, these pieces of information could be used without the risk of giving them full weight, as it would occur if they were used for producing fixed nursery-signature parameters, under the conditional approach. Moreover, these data, although limited, could be used to sustain some resampling procedures aimed to explore the risk of biased conclusions under plausible scenarios, as done in the present work.

In conclusion, unconditional ML-MM showed to be a suitable tool for estimating mixing proportions and nursery signatures of S. aurata under multiple scenarios that included incomplete sampling and a range of chemical signature separations among nursery-sources. In contrast, our approach yielded rather discouraging results regarding the estimation of the true number of nursery-sources (K), under incomplete sampling and/or identification of nursery-sources. Therefore, new efforts aimed to develop new mixed-stock analysis tools and/or to evaluate the performance of the existing ones are required. Such evaluations should be conducted over the widest possible range of species, habitats and biological scenarios, both to improve the reliability of these tools and to enhance our understanding about the structure and connectivity of exploited and threatened fish populations.

Supplemental Information

Methodological details.

DOI: 10.7717/peerj.2415/supp-1

Descriptive statistics.

DOI: 10.7717/peerj.2415/supp-2

Appendix 3.

Main analyses (R script in pdf format).

DOI: 10.7717/peerj.2415/supp-3

Elemental fingerprints data.

DOI: 10.7717/peerj.2415/supp-4
5 Citations   Views   Downloads