A Bayesian Multilevel Ordinal Regression Model for Fish Maturity Data: Difference in Maturity Ogives of Skipjack Tuna (Katsuwonus pelamis) Between Schools in the Western and Central Pacific Ocean

The maturity ogive is vital to defining the fraction of a population capable of reproduction. In this study, we proposed a novel approach, a Bayesian multilevel ordinal regression (i.e., Bayesian continuation ratio model), to model the maturity ogive. The model assumes that the observed maturity stage originates from the categorization of latent continuous variables. We demonstrated this approach by testing whether there are differences in the maturity ogive of skipjack tuna (Katsuonus pelamis) in the western and central Pacific Ocean between two school types, i.e., free-swimming and floating-object-associated schools. The model results show that K. pelamis, given the same fork length, are more likely to have a higher maturity stage in a free-swimming school than those associated with floating objects. The gonadosomatic index revealed the same conclusion. Our results indicate that fish aggregation devices (FADs) could negatively affect the maturity of K. pelamis and consequently reduce the population reproductive potential. This study provides (1) an alternative approach to analyze fisheries ordinal data; (2) important quantitative evidence to evaluate the existing ecological hypotheses; and (3) implications for tuna fisheries management.


INTRODUCTION
The maturity ogive, which measures the proportion of mature fish at age or size, is instrumental to defining the fraction of the population that is capable of reproduction, commonly referred to as the spawning stock biomass (SSB) in a population dynamic context. It is highly influential upon the results of stock assessments and determining reference points for management (King and McFarlane, 2003). The adequate estimation of the maturity ogive is a critical step for assessment and management and typically requires (1) representative samples of ovaries from a well-designed sample adapted to a species-specific reproductive strategy; (2) precise determination of the stage of sexual maturation; (3) fitting a statistical model to the data from maturity staging. The maturity stage is commonly examined using either microscopic or macroscopic staging (Flores et al., 2019). Macroscopic staging is based on visual evaluation of the external gonad development and could be uncertain and inaccurate due to a lack of differentiating visible structures between stages; for instance, immature and regenerating/omitted spawning stages have very similar macroscopic aspects of the ovaries and are often misidentified. Microscopic staging is considered to be accurate and provides an unambiguous determination of maturity stages. However, histology is expensive, time consuming and impractical at sea during routine surveys (Flores et al., 2019). Other methods such as gonometric method that is based on the gonodasomatic index (GSI) has also been evaluated as an alternative to maturity staging (Flores et al., 2014(Flores et al., , 2019. Ordinal data like maturity stages arise when underlying continuous processes are censored in the observational process, where measurements are limited to discrete categories. This type of data is not uncommon in fisheries science. However, ordinal regression, which models the influence of latent effects on an ordinal outcome, is not commonly used in fisheries science. There are different types of ordinal regressions to model an ordinal response depending upon how their logits are formed and underlying assumptions (Bürkner and Vuorre, 2019): proportional odds model (POM), continuation ratio model (CRM), and adjacent category model (ACM). POMs, considered as cumulative higher categories verses cumulative remaining lower categories, are often used when assuming responses are from categorization of a latent continuous variable. However, CRMs (considered as cumulative higher category verses just lower category alone) and ACMs (between any of two consecutive categories) can assume that for every category there is a latent continuous variable that determines the transition of successive categories. The POM focuses on the cumulative odds of grouped categories, whereas the CRM estimates the conditional odds of being in a particular category, given that an individual has reached that category or above (Liu et al., 2011).
The maturity ogive is conventionally modeled as a function of age and/or size using logistic regression, often collapsing maturity stages into two levels (i.e., immature and mature). Such a dichotomization throws away important information in the data that are ordinal measurements, i.e., maturity stages, and the resulting estimates of the maturity ogive can be sensitive to the maturity scales and "cut point" used to determine whether an individual is matured or not (Núñez et al., 2015). Ordinal regression models making full use of the ordinality of the data is an approach that could fit maturity stages in a more natural way and potentially improve the accuracy of SSB calculation. However, to our knowledge, this method has not been used for modeling fish maturity data.
Skipjack tuna (Katsuwonus pelamis) are widely distributed in tropical and subtropical waters throughout the Western and Central Pacific Ocean (WCPO) and support important fisheries. In 2018, the catch of world's major commercial tunas was 5.2 million tons, of which the majority was skipjack tuna and caught via purse seining (ISSF, 2020). There are two types of tuna schools which purse seiners target: schools associated with floating objects including fish aggregation devices (FADs) (hereinafter referred to as an associated school) and free-swimming schools. About half of the total tuna catches in the WCPO were from associated schools, although the number of sets around FADs was much less than that on free-swimming schools (WCPFC, 2019). Capturing schools of tunas associated with floating objects was pioneered by tuna purse seiners operating in the eastern Pacific Ocean during the late 1950s (Hampton, 1993). Artificial FADs have been increasingly constructed and released to the global ocean since 1980s because floating objects significantly increase the catchability of tunas, i.e., purse seine sets around floating objects yield much higher success rate than those made on free-swimming schools (Suzuki et al., 2003).
The potential negative impacts of using FADs have gained a significant amount of scientific attention, including growth overfishing of tuna stocks (FAD catches are usually characterized by small individuals, whereas catches on free-swimming schools are typically dominated by large tuna), increased bycatch, and alteration of the life history traits of the species associated with FADs Bromhead, 2003;Hallier and Gaertner, 2008). The potential effects of floating objects on the behavior and biology of tunas have been widely studied, and hypotheses have been developed to explain ecological consequences, e.g., "ecological trap" (Marsac et al., 2000;Hallier and Gaertner, 2008). This hypothesis contends that releasing a large number of FADs in the ocean could prevent tuna from moving to places where they would go normally, lead them to low-quality habitats and consequently negatively impact their biology. For example, Jaquemet et al. (2011) found that K. pelamis and small yellowfin tuna (Thunnus albacores) caught around drifting FADs had a significantly higher frequency of empty stomachs and concluded that FADs could negatively impact the growth of K. pelamis and small T. albacores. Ashida et al. (2017) found that the relative condition factor and proportion of mature fish significantly differ between associated and free-swimming schools (lower in associated schools). However, K. pelamis and T. albacores are considered as income breeders (Zudaire et al., 2014;Grande et al., 2016), which means they rely on concurrent energy income from feeding to support the cost of reproduction. Therefore, the condition factor might not be directly related to maturation. Assuming the "ecological trap" hypothesis is true and considering tunas are income breeders, one would hypothesize that fish associated with floating objects for an adequate period of time would have a lower maturity stage compared with those in free-swimming schools, given the same age/size. However, there are few studies that have compared the size-at-maturity between tuna associated with FADs and free-swimming schools. This may be because most of the tuna catches (e.g., bigeye tuna (Thunnus obesus) and T. albacores) in floating object sets consist of juvenile individuals (Bromhead, 2003;Dagorn et al., 2013). Ashida et al. (2017) compared the difference in reproductive traits of K. pelamis females between the two schools, but they did not examine whether there was a difference in size-at-maturity between the two schools.
In this study, we proposed a novel approach: a Bayesian multilevel ordinal regression, to model the maturity ogive using maturity stage data. We demonstrated this approach by testing the hypothesis that there are differences in maturity characteristic between the associated and free-swimming schools. We also compared the results of ordinal regressions with those obtained from the conventional logistic regressions.

Data
A total of 2013 K. pelamis were sampled on Chinese commercial tuna fishing vessels in the WCPO during 2007-2014 from 154 fishing locations (Figure 1). The sampling date, location (latitude and longitude) and the school type (i.e., free swimming school or associated school) were recorded for all specimens. All specimens were measured for fork length (FL, mm to the nearest 0.1 mm) and gonads were removed and weighed. Body weight was measured when possible. The majority of the specimens were collected from associated schools (n = 1655). Gonadal stage of maturity was assigned for each specimen based on observation of the morphological features of the gonad within the body cavity. We followed the general maturity criteria for large pelagic gonads and established six developmental stages for both female and male ( Table 1): stage-1: Gonads are small ribbonlike, not possible to determine sex by gross examination; stage-2: Immature; stage-3: Early maturing; stage-4: Late maturing; stage-5: Ripe; and stage-6: Spent/Spawned. Sex was determined for each individual based on visual examination of the gonad. Of the total specimens, 929 specimens were female, 830 specimens were male. The sex of the rest 254 specimens (stage-1) could not be identified. We excluded the data of individuals identified as Spent/Spawned (i.e., stage-6) in this analysis due to the high sprobability of misidentification (Figure 2).

Bayesian Multilevel Ordinal Regression
The ordinal regression model is motivated by assuming underlying latent continuous variables that are related to the  ordinal response through the threshold concept (Hedeker and Gibbons, 1994). Because the maturity data are ordinal in nature, instead of a threshold value such as the cut-point to separate mature from immature fish as in the conventional logistic regression approach, there are several thresholds (depending on the number of maturity stages) that separate individuals into the various response categories, i.e., maturity stages. Gonadal maturation is a sequential biological process throughout the reproductive cycle, i.e., individuals must pass through lower stages in order to reach higher stages. Therefore, a CRM was considered in this study. Fisheries data, including maturity data, are often structured in a hierarchical manner. For example, observations of the maturity stage taken from an individual fish are nested within sampling sites. Maturity stages within a sampling site are likely to be more similar to each other compared to observations between sampling sites. Such multilevel structure would violate the assumption of independence and should be accounted for in the analyses.
Multilevel models have been extensively developed to model hierarchical data (Wagner et al., 2011). In this analysis, we proposed a multilevel CRM to estimate the maturity ogive. The CRM estimates the odds of being in a particular stage k relative to being above that stage and can be formulated as: where P Y = k|Xβ + Zµ is the conditional probability of being in stage k, conditional on being that stage or beyond, given the linear predictor (Xβ + Zµ), k = 1, 2, 3, . . . , K − 1, τ k are the thresholds or cut points, β are the regression coefficients at population-level (i.e., fixed effects), X are populationlevel covariates, Z are grouping variables, µ are group-level coefficients, which are also known as random effects that follow a normal distribution with a mean of zero. We assumed that β are constant across all response categories, i.e., the effect of each predictor is invariant across the ordinal responses. The covariates we considered in the ordinal regression were fork length (FL), school type (School), fishing location (location), year (Year) and month (Month). We excluded sex as a covariate in the model because sex could not be identified for stage-1 fish. The linear predictor of the full model was specified as where FL and School are population-level effects, and location, Year and Month are group-level effects, and Month is nested within Year. We included an interaction term to examine whether the effect of fork length on maturity is independent of school type. We conducted model comparisons of the full model to reduced models to better understand the relative importance of model terms (Supplementary Table 1).
All models were fitted with the "Bayesian Regression Models using Stan" (brms) R package (Bürkner, 2017). We chose a weakly informative prior, i.e., Normal(0,10) for β and τ k , and a Half-Cauchy prior, i.e., Cauchy(0,5) for the random effect standard deviations. We ran 4 chains for 1,500 warm-up iterations followed by an additional 3,000 iterations. We confirmed algorithm convergence with visual checks (e.g., traceplots) and the Rhat statistic. We chose the leave-one-out information criterion (LOOIC), a Bayesian information criterion based on out-of-sample predictive performance, to compare models fitted to the same data (Vehtari et al., 2017). The model with a lower LOOIC is expected to have better predictive accuracy (Hooten and Hobbs, 2015). We used the "LOO" R package to compute the LOOIC (v2.4.1; Vehtari et al., 2020). To test our hypothesis that fish in free-swimming schools on average have a higher maturity stage compared with those associated with floating objects given the same age/size, we calculated the evidence ratio of the posterior probability of School > 0 (i.e., school type has a positive effect on maturity) and the posterior probability of School < 0.

Sensitivity Analysis
We tested the sensitivity of the estimated maturity ogive to different types of models and data. First, we tested the sensitivity of results to the potential misidentification of maturity stages by running the analysis using the full maturity dataset including individuals identified as stage 6. Second, we fitted both ordinal and conventional logistic models with the same covariates to the same datasets. We chose a parsimonious model, M7, for the sensitivity analyses. The sensitivity analyses resulted in four model runs: ordinal and logistic regression fitted to full dataset and the dataset with stage-6 fish excluded, respectively. We considered fish in stages 4-6 to be matured fish. We also investigated that whether the logistic model with the same covariates as the full model would lead to the same conclusion regarding the effect of school type based on ordinal regression.

Gonadosomatic Index Difference Between Schools
The gonadosomatic index (GSI), measuring the ratio between gonad weight and body weight, can indicate maturity and stage of spawning and is commonly used to compare reproductive condition among individuals (Flores et al., 2014(Flores et al., , 2019. If there are differences in the maturity condition between the two schools, it is expected that the GSI also differs between the two schools. We therefore tested this hypothesis using a Bayesian generalized linear model where GSI is the ratio between gonad weight and gutted weight, β 0 , β 1 , β 2 and β 3 are population-level effects and are given a weakly informative prior, i.e., Normal(0,10); α location is the site-level effect, which follows a normal distribution with a mean of zero. The standard deviation is given a Half-Caughy prior, i.e., Caughy(0,5). The model was also fitted with the "Bayesian Regression Models using Stan" (brms) R package (Bürkner, 2017).

RESULTS
A total of 1,501 specimens (stages 1-5) were used in the analyses. The majority of them were from associated schools (1,202 specimens). The range of FL of skipjack tuna from associated and free-swimming schools were (27-73.3 cm) and (22.8-74.6 cm), respectively. The median FLs of skipjack tuna from associated and free-swimming schools were 44.8 and 53.2 cm, respectively (Figure 2). The median FL increased as tuna progressed to higher maturity stages, except stage 6 (Figure 2). This indicates that the likelihood of misidentifying spawned tuna is high. Model selection indicated that models with the random effect grouping variables year and sampling station tended to fit the data much better than models without them (Supplementary Table 2). This suggested that maturity of skipjack tuna varies across time and space. The full model had the lowest LOOIC, suggesting the most parsimonious fit. The mean and standard deviation of the posterior distribution for each parameter of the full model are summarized in Table 2. Rhat values suggested that the model had converged. The thresholds (Table 2) describe the four cut points of the model applied to the maturity stage data and indicate where the continuous latent variables were partitioned to produce the observed responses. The estimated mean of the coefficient of FL indicated that fork length had a positive effect on maturity stage, and its 95% credible interval did not contain zero. The effect of FL does not tend to be different between the two types of school as indicated by the posterior distribution of the coefficient of the interaction between FL and school type.
Because school type is coded as a factor with associated school as the reference category, the coefficient of school type indicates the extent to which fish of free-swimming schools differ from those of associated schools on the latent scale maturity stage. The estimated mean indicates that there was a strong effect of school type on the maturity of skipjack tuna, i.e., the maturity of fish from free-swimming schools appeared to be higher than that of fish from associated schools ( Table 2 and Figure 3). The 95% credible interval of this parameter was between 0.78 and 7.83, which did not contain zero. The evidence ratio suggested that the posterior probability of School>0 (i.e., free school has a positive effect on maturity) was 99%. We can therefore conclude with 99% probability that fish from free-swimming schools exist at a higher maturity stage than those associated with floating objects given the other covariates in the model. For any given fork length, K. pelamis from a free-swimming school are more likely to be in higher maturity stage than those associated with floating objects (Figure 3). For example, given a length of 50 cm, K. pelamis from free-swimming schools are more likely to be in mature stages (i.e., stages 4 and 5) than those from associated schools (Figures 4, 5).
Estimated maturity ogives were sensitive to the data used in the analyses. Excluding stage 6 fish in both ordinal and logistic regression analyses steepened the slope of the curve; proportion mature increased faster as fork length increased (Figure 6). This suggested that some of the fish in stage 3 (immature) were likely misidentified as stage 6 (mature). Therefore, the predicted proportion for small fish, e.g., 40 cm, was higher when stage 6 fish were included in the analysis. The ordinal and logistic regressions yielded very similar maturity ogives when using data without stage 6 fish. However, when stage 6 fish were included in the analyses, the logistic regression produced a higher proportion mature for large fish (e.g., >45 cm) compared with the ordinal regression. The logistic regression with full covariates (i.e., M0) did not detect a significant difference between the two school types (Supplementary Table 3).
The GSI model results indicated a significant difference in GSI between the two schools, i.e., the point estimate of the coefficient of school type is 1.04 and its 95% credible interval did not contain zero ( Table 3). The evidence ratio suggested that the posterior probability of School>0 (i.e., free school has a positive effect on GSI) is 99%. For a given fork length, fish from free-swimming schools have a larger GSI than those associated with floating objects (Figure 7). As expected, fork length had a significant effect on GSI. However, this effect is independent of school type, as the 95% credible interval of the coefficient of the interaction term between fork length and school type contained zero.

DISCUSSION
In this study, we proposed a novel Bayesian multilevel ordinal regression model to model maturity ogives and applied it to K. pelamis in the WCPO to examine whether maturity characteristics of K. pelamis differs between two school types: schools associated with floating objects and free-swimming Bulk_ESS is an estimate of the effective number of independent draws from the posterior distribution of the estimate of interest. Rhat measures the ratio of the average variance of samples within each chain to the variance of the pooled samples across chains. The Estimate column contains the posterior means of the parameters. Est.Error contains the parameters' posterior standard deviations. l-95% CI and u-95% CI are the 95% credible intervals.   schools. The model results showed that there is a significant difference in maturity conditions between the two schools. K. pelamis from free-swimming school were more likely to have a higher maturity stage than those associated with floating objects given the same fork length. This was supported by the results of the GSI analysis, which also showed a significant difference in GSI between the two schools, i.e., that the GSI of K. pelamis in free-swimming school was significantly higher given the same fork length. Together, our results provide evidence that floating objects, including FAD, affect the maturity ogive of K. pelamis, and individuals in free-swimming schools may be maturing earlier. Therefore, it is important to account for such FIGURE 6 | Fitted maturity ogives based on ordinal and logistic regressions using data with and without stage 6 fish.
differences when estimating maturity ogive for stock assessment and management. Our results also suggest that estimates of maturity ogive are sensitive to maturity data as well as models used in the analysis. The conventional logistic regression model did not detect a significant difference in maturity between the two school types and could underestimate the L 50 (length at which 50% of the fish are maturity), which leads to overestimation of maturity ogive, compared with the ordinal regression model. A previous study also found that reproductive biology of Red Snapper in the Gulf of Mexico differs significantly between natural and artificial habitats (Glenn et al., 2017). Red Snapper at natural habitats may be maturing earlier than those at artificial habitats because they can afford to trade somatic growth for reproductive potential. The impact of floating objects on K. pelamis's maturity could be a result of the negative impacts of FADs on the rates of food intake and the breeding strategy that K. pelamis utilize. Studies based on lipid and stable isotope analyses suggest that K. pelamis are income breeders (Grande et al., 2016). A income breeding strategy can adjust the reproductive schedule; for instance, the spawning season of K. pelamis in the western Indian Ocean matches with monsoon periods when food availability is high (Roger, 1994;Schott et al., 2002;Potier et al., 2008;Sardenne et al., 2016), and subsequent egg production of T. albacares increased in response to increasing daily rations (Margulies et al., 2007;Zudaire et al., 2015) relative to the current food availability. Studies have found that batch fecundity, egg size, and number of eggs spawned can change depending on the food availability (Ashida et al., 2017). Aggregation around FADs is believed to decrease the food intake for tuna species Hallier and Gaertner, 2008). In all the specimens that we have examined for stomach contents during 2007-2014, the percentage of empty stomachs for K. pelamis from associated schools was 94.7% (1,566 out of 1,653). This was far higher than that from free-swimming schools, which was 29.6% (106 out of 358). Such a difference between the two school types has also been observed in other studies Potier et al., 2001;Hallier and Gaertner, 2008;Jaquemet et al., 2011). However, such a large difference in percentages might be due in part to the time of day when the fishing operations were conducted. The commercial fishing vessels targeted associated schools at dawn when the food in stomach had already been digested because tuna have exceptionally high digestion rates (Olson and Boggs, 1986), whereas the fishing vessels targeted unassociated tuna schools primarily during the day and often during their feeding time. Therefore, the percentage of empty stomachs in this case might not accurately reflect the difference in food intake between these two types of schools. Our results suggested that K. pelamis associated with floating objects need more time for gonad development. This is likely due to poorer feeding conditions where FADs are typically deployed. There is a tendency for some fishermen deploy FADs in areas where productivity is low (Bromhead, 2003) or using materials other than natural logs (Leroy et al., 2013).
One caveat of this study is the absence of data on residence time around FADs. For tropical tunas, some tagging experiments revealed that the residence time around individual FADs are short (i.e., 2-8 days) (Hilborn, 1991;Schaefer and Fuller, 2013;Matsumoto et al., 2014). However, the residence time is subject to a degree of variability depending on species and the region (Baidai et al., 2020). Using echosounder buoys data, Baidai et al. (2020) found that time span of tuna aggregations at FAD followed a time-independent process with short-and long-term residence modes, and FADs remained occupied for a large proportion of time after being colonized. Therefore, it is entirely possible that skipjack reside around FADs continuously or discontinuously for an adequate period of time that could affect their maturity.
Scientific studies have raised the concerns about the potential negative impacts of using FADs on the biology of tuna species Potier et al., 2001;Hallier and Gaertner, 2008;Jaquemet et al., 2011;Wang et al., 2012Wang et al., , 2014Robert et al., 2014;Ashida et al., 2017). Our study also suggests that FADs could negatively affect the maturity of K. pelamis. From a management perspective, FADs may not only lead to growth overfishing but also reduce the population's reproductive potential. Our results are consistent with previous studies. Ashida et al. (2017) compared the difference in reproductive traits of female K. pelamis between schools in the same area as this study. They found that the proportion of mature fish was significantly different between associated and unassociated schools. Although difference in size at maturity was not examined, they proposed that poorer food availability near FADs may impact oogenesis in K. pelamis. The difference in size at maturity may further confirm this speculation. They also concluded that inherent spawning characteristics (sex ratio, observed minimum size at first maturity, and spawning season) were not affected between schools. This is not surprising given that a FAD is more likely to affect food availability and consequently reduce reproduction.
Marine fisheries management is often based on biological reference points, which are predicted on an assumed relationship between stock reproductive potential and subsequent recruitment. The estimate of spawning stock biomass, which is based on maturity ogive, is often used to approximate reproductive potential. Therefore, imprecise or biased estimates of maturity ogive can lead to fishery management decisions based on mis-specified biological reference points. The findings of this study have important implications for skipjack tuna assessment and management. The estimated maturity ogive based on samples collected from either free school or associated school in isolation would likely to be biased. Our results also suggested that it is important to account for spatiotemporal variability when estimating the maturity ogive. Therefore, a comprehensive survey design that allows for representative samples to be collected is essential to estimate the maturity ogive for assessment and management use.
The ordinal regression approach proposed in this study fits maturity data obtained from both microscopic and macroscopic approaches naturally and can estimate the probability of being in each of the maturity stages. Thus, the maturity ogive can be derived by summing the probability across maturity stages that are considered mature. This allows for sensitivity analyses on maturity ogive in the assessments without refitting the regression model. Since the ordinal regression approach makes full use of the ordinality of the data, we hypothesize that this approach may be more robust to the misspecifications of maturity stages. However, this needs to be further tested. Therefore, future studies should evaluate the performance of ordinal regressions compared with conventional logistic regression on estimating maturity ogive, especially based on the data from a macroscopic method where misspecification errors are more likely to present, using a simulation approach. Finally, we envision that this ordinal regression approach could be applied more broadly in fisheries science when observations are taken as ordinal measurements, e.g., feeding levels.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by Academic Committee of Shanghai Ocean University.