Linear mixed-effects models to describe length-weight relationships for yellow croaker (Larimichthys Polyactis) along the north coast of China

In this study, length-weight relationships and relative condition factors were analyzed for Yellow Croaker (Larimichthys polyactis) along the north coast of China. Data covered six regions from north to south: Yellow River Estuary, Coastal Waters of Northern Shandong, Jiaozhou Bay, Coastal Waters of Qingdao, Haizhou Bay, and South Yellow Sea. In total 3,275 individuals were collected during six years (2008, 2011–2015). One generalized linear model, two simply linear models and nine linear mixed effect models that applied the effects from regions and/or years to coefficient a and/or the exponent b were studied and compared. Among these twelve models, the linear mixed effect model with random effects from both regions and years fit the data best, with lowest Akaike information criterion value and mean absolute error. In this model, the estimated a was 0.0192, with 95% confidence interval 0.0178~0.0308, and the estimated exponent b was 2.917 with 95% confidence interval 2.731~2.945. Estimates for a and b with the random effects in intercept and coefficient from Region and Year, ranged from 0.013 to 0.023 and from 2.835 to 3.017, respectively. Both regions and years had effects on parameters a and b, while the effects from years were shown to be much larger than those from regions. Except for Coastal Waters of Northern Shandong, a decreased from north to south. Condition factors relative to reference years of 1960, 1986, 2005, 2007, 2008~2009 and 2010 revealed that the body shape of Yellow Croaker became thinner in recent years. Furthermore relative condition factors varied among months, years, regions and length. The values of a and relative condition factors decreased, when the environmental pollution became worse, therefore, length-weight relationships could be an indicator for the environment quality. Results from this study provided basic description of current condition of Yellow Croaker along the north coast of China.


Introduction
Length-weight relationships (LWR) are used for estimating the weight corresponding to a given length [1]. As most observations in fisheries surveys are typically obtained as the number of specimens and length of each sampled specimen, LWR are widely used to transform the survey data into estimates of the biomass which is important for modelling aquatic ecosystems. Therefore, estimation of the LWR is common practice and essential in fisheries science [2].
There are two parameters in the LWR model (W = a×L b ), the coefficient a and the exponent b. Parameter a is the condition factor, describing body shape which can be classified as four groups: eel-like, elongated, fusiform and short and deep [2]. Parameter b is the allometric growth parameter, which indicates isometric growth in body proportions if b!3 where fish have more girth as it grows longer; the species tends to be more streamlined if the exponent b<3 [2].
Numerous studies found that factors, such as geographic, seasonal, inter-annual, and environmental conditions, can affect the estimates of a and b in the LWR model [2][3][4][5]. Instead of constructing several models reflecting various situations (different regions and years), it is more reasonable to make use of generalized linear mixed model to cover all the spatial and temporal effects in a single model. Linear mixed-effects modelling, a mature model in the statistics community, has been used in multiple fields. Cnaan et al. (1997) provided two detailed case studies to sufficiently introduce the use of the general linear mixed effects model for the regression analysis of correlated data [6]. Xu et al. (2015Xu et al. ( , 2014 developed nonlinear mixedeffects model to study the individual-tree diameter growth and linear mixed effects model for individual-tree crown width of China-Fir trees in Southeast China [7,8]. Baayen et al. (2008) provided an introduction of mixed-effects models and illustrated its advantages, which would allow the researcher to simultaneously consider all factors that potentially contribute to the understanding of the structure of the data [9].
The linear mixed-effects model provides a powerful and versatile approach to analyze a wide variety of data structures, in which the linear predictor contains random effects in addi-specific models, rather than linear mixed-effects models, to detect the spatial and temporal variations of Yellow Croaker.
In this study, our objective was to analyze the biological characteristics of Yellow Croaker in recent years. Specifically, we estimated the LWR of Yellow Croaker used linear mixedeffects model, condition factors relative to reference years, and explore the possible relationship between parameters in LWR model and environmental factors, based on the observations along north coast of China among 2008 and 2011 to 2015.  Table 1, S1 Table, Fig 1), which are important spawning and feeding grounds of Yellow Croaker [19]. The maps of surveying area were plotted using package maps and mapdata of the R (version: R i386 3.3.1) [20,21].

Data collection
In Coastal waters of Northern Shandong, specimens were randomly collected from fishermen immediately after landing. For all other five regions, stratified random bottom trawling surveys were implemented to collect samples. In total, 3,275 individuals were collected, with sample size variations among years and regions (see Table 1 and S1 Table). For each individual, standard length was measured to the nearest mm and total weight was measured to the nearest gram.
The surveys, which held in the Seasonal Fishery Closure (June, July and August) were approved by the Department of Marine and Fishery of Shandong and Jiangsu Province. No specific permissions were required for all the other surveys, because these surveys are traditional surveys, which did not cover any marine protected area or private area and this study did not involve endangered or protected species.

LWR modelling
We fitted an exponential function to the weight at length data that took the form [1]: where W is the wet weight of an individual fish (g), L is the standard length (cm), a is the condition factor, and b is the allometric growth parameter. Because the variance of W increases when L increases, above equation was log-transformed and the equation becomes: We fit a generalized linear model (GLM), simple linear models for individuals in different regions and years (SLMR and SLMY for regions and years respectively), as well as nine linear mixed effect models (LMEM) that treated region and/or year as random effects to coefficient a and/or the exponent b to explain the relationship between length and weight (Table 2) [7]. Analysis of Variance performed by F-test between GLM and LMEM was computered to tests whether the temporal and spatial variations were significant. All these modelling processes were conducted, using package lme4 the R (version: R i386 3.3.1) [22]. Bootstrap the data with replacement for 1000 times was used to estimate the distributions and statistics of parameter estimates in these models.
We explored the influence of marine environmental status on variability of the LWR. Environmental data were drawn from Bulletin of Marine Environmental Status of China, which were reported by State Oceanic Administration People's Republic of China [23]. According to Sea Water Quality Standard (National Standard GB 3097-1997) [24], five levels (the first, second, third, fourth and worse than fourth levels) were used to evaluate the water quality, in which higher level indicated poorer quality. The polluted water area percentage (PWAP), which included the third, fourth and worse than fourth level water, was used to exemplify water quality. We explored the influence of PWAP on the LWR along north coast of China during 2008, 2011 to 2015.  Table 1.

Model comparison
Both the Akaike information criterion (AIC) and the mean absolute error (MAE) were calculated to compare the performance of the twelve candidate models. AIC was widely used to compare the quality of models [25,26]. Lower AIC value indicates a better model. MAE is a quantity used as a measurement on how close forecasts or predictions are to the eventual outcomes [27].
MAE was the average of the absolute errors: where f i was the estimation and y i the observed value.

Relative condition factor
Relative condition factor (K) was calculated for each individual [28]: where W and L were observed weight and standard length from our surveys, baseline model LWR with parameters a and b were drawn from historical references. Relative condition factor was recommended to measure the deviation of an individual from the baseline weight for length in the respective sample and to investigate changes in the stock over time [2,28]. Here the LWR models of Yellow Croaker in 1960,1986,2005,2007,2008.04~2009.03 (12 months) and 2010 were introduced to be the baselines (S2 Table) [3,11]. In 1960, Yellow Croaker were Table 2. Models for length-weight relationships of Yellow Croaker. abundant and the stock was not deeply disturbed by fishing [13]. However, Yellow Croaker experienced severe fishing pressure and reached its lowest biomass in 1986 [14]; after recovery, the catch of the stock was increasing in recent years [12,14].

Length and weight of Yellow Croaker
The

Recommended LWR
The LMEM (R&Y.I&S) with Regions and Years have random effects on both intercept and slope was recommended to be used as the best model after comparing with the other models based on AIC and MAE values (-4357 and 0.095, respectively) ( Table 2). According to this selected model, a and b were estimated as 0.0192 (95% CI = 0.0178, 0.0308) and 2.917 (95% CI = 2.731, 2.945), respectively ( Table 3). The estimates of a with the random effects from Region and Year ranged from 0.013-0.023, while the estimates of b with the random effects from Region and Year ranged from 2.835-3.017.

Spatial and temporal variations of LWR
The selected model with random effects from both regions and years applied to both a and b, indicated that there were substantial spatial and temporal variations of LWR for Yellow Croaker (Table 4,  In Yellow Sea, the PWAP occupied 3.2~7.1 percentage of the whole area during the survey time, while the PWAP covered 17.7~37.2% of the Bohai. In addition, the water quality of Chinese coast varied widely among years. When the water quality became poorer, condition factor a decreased, with the reverse condition for parameter b (Fig 4). Negative correlation (-0.47) and positive correlation (0.61) were found for PWAP with condition factor a and for PWAP parameter b respectively, while t-test indicated both the correlation coefficients were not significant (df = 4, P>0.05).  Table 1.    6). Without data in 2012, correlation analysis indicated negative correlations for PWAP with condtion factors relative to reference years except 1960 (from -0.68 to -0.61), while no significance were found in t-test (df = 3, P>0.05).

LWR of Yellow Croaker and its spatial-temporal variations
According to the results of LMEM (R&Y.I&S), the estimates of a (0.013~0.023), fell in the common range (0.001~0.05) of most fish species [2]. Study of Froese [29] shows that values of Fish spatial-temporal mixed effect growth model parameter a of fish species, with fusiform body shape, was 0.0112 (geometric mean) with 95% range of 0.00514 to 0.0245. The results are consistent with the body shape of this species, which is described as fusiform in Fishbase [30]. The point estimate of exponent b, was a little lower than values 3.03 (2.91-3.15) in Fishbase [29]. Values of parameter b for fishes were suggested by Carlander [31] and Froese [2] to normally fall within the range of 2.5~3.5 and 2.7~3.4, respectively, which cover the range of parameter b of Yellow Croaker in this study. The fixed value of allometric growth parameter b, indicated a very slight decrease in plumpness or elongation in form with increase in length. However, when random effects were used, no significant difference with isomeric growth was found for Yellow Croaker [32].
Random region effects for a demonstrated that values of a reduced gradually as the latitude decreased, from north to south, except in Coastal Waters of Northern Shandong. Samples in the Coastal Waters of Northern Shandong were collected from fishermen, and the origin of the samples were not clear, which could be the reason that parameters with the random effect of NS did not follow this trend. Yellow Croaker distributed on China's coast were considered to be divided into various populations, while individuals in northern Yellow Sea and Bohai Sea (NYBS) were believed to be different from those in southern Yellow Sea (SYS) [19,33,34]. The results of this study suggested that Yellow Croaker sampled from Yellow River Estuary and Coastal Waters of Northern Shandong should be included in NYBS, while those in other regions belonged to SYS. Different subpopulation maybe the main reason that led to their significant spatial variation in length-weight relationships.
Previous research has examined LWR variation for Yellow Croaker; including spatial, temporal, sexual, and length variations, which covered Bohai Sea, Yellow Sea and East China Sea from 1960 to 2010 [3,11,14,17,18,[35][36][37][38][39][40][41]. In these LWRs (51 models), values of parameter a Fish spatial-temporal mixed effect growth model ranged from 0.0061 to 0.1027, with mean of 0.028 ± 0.019; values of exponent b ranged from 2.32 to 3.35, with mean of 2.875 ± 0.217. This, accompanied with relative condition factors, could support the hypothesis, that individuals were becoming thinner and young-age in recent years.
Several influencing factors of LWR were evaluated in previous studies [3,42,43]. Pollutions such as inorganic nitrogen, active phosphate and oil pollutant, frequent pollutants in the north coast of China, are believed to make food availability and oxygen level deteriorate, which would negatively influence the growth of Yellow Croaker [44][45][46]. In this study, the data revealed that values of condition factor a decreased, when the environmental pollution increased. Therefore, LWR could be a rough indicator for the environment quality. Limited sample years maybe the reason that resulted in the non-significant correlation from statistics test, so further long term monitoring is suggested in the future. Fish spatial-temporal mixed effect growth model

Condition changes compared with historical records
Relative condition factor provided an effective approach to compare the observed weight of an individual with the baseline weight of that length [28,47]. Compared to individuals in previous decades [3,11], Yellow Croaker condition declined in recent years. K cur/1960 remained comparatively stable, without obvious variation. K cur/2005 , K cur/2007 , K cur/2008 9 and K 2010 had similar variation patterns, since the LWR of Yellow Croaker in these years were similar, while values of K cur/2005 were much lower. However, K cur/1986 have similar variation pattern but higher variation, since the stock of Yellow Croaker collapsed during the 1980s and had been in an unstable condition in 1986 [16,48].
Condition factors relative to all reference years, reached the highest value in 2012, as specimens in 2012 were collected in May and September (S1 Table), which was just before spawning and after feeding. A weak El Nino in 2012 might be another possible reason that led to elevated relative condition factors.
Monthly changes in condition factors are related to life history events of this fish, e.g., sexual maturation, spawning and feeding strategies [14,49,50]. Adult yellow croaker in these areas reach their full gonads maturity in spring, reproduce in June and then primarily feed in Fish spatial-temporal mixed effect growth model summer and autumn [49,51]. The condition factor (K) behaved in accordance to this series of events. It reached a high value in May, dropped precipitately in June, gradually recovered in August and reached the highest value in December.
Spatial variation of K discovered that Yellow Croaker was skinny in the south, with worse condition than in the north. It was widely reported the Yellow Croaker stock shifted in size and age structure to feature more small and young fish that mature earlier, in response to severe fishing pressure over the past half century [15,16,40,48]. This circumstance was also reflected in the length variation of condition factor relative to reference year of 1960.
Since there was only one sample of L = 21cm, the condition factor relative to all reference years became an exception in L = 21cm for length variation of Yellow Croaker. Moreover, just as values of parameter a, water pollution led to K worse. The exception value K in 2012 relative to all reference years, probably resulted from the more significant effects from monthly variation and nearly null El Nino in 2012 (Fig 5).
Applications of mixed effects model in studying spatial-temporal variations of LWR According to AIC and MAE, linear mixed effect model with effects from both regions and years on both parameters a and b was supported as the best one. The recommended model matches previous studies that the exact relationship between length and weight differs among regions and years, according to the habitat condition [4,5,52,53]. However, many studies on LWR for fishes used region specific or year specific models fitting to region or year specific dataset to detect the spatial and temporal variations of individuals [4,54,55], as seen from previous 51 LWR models for Yellow Croaker [3,11,14,17,18,[35][36][37][38][39][40][41]. The mixed effects model took account of regions and years as random effects in a single model, which was more convenient and reasonable to estimate the spatial and temporal variations. It also allow regions and years with limited samples to borrow strength from other regions and years [10].
Furthermore, mixed effects model provided an effective approach to estimate the effects of variables, observations of which could be unavailable [7,8]. In this study, our survey data were quite poor: there were 23 blanks in the 6×6 (Regions-Years) data matrix (Table 1), while mixed effects model effectively evaluated the effects of regions and years. In the future, more variations of LWR, such as sex, season, and growth stage should be considered in the mixed-effects model for Yellow Croaker.

Conclusions
According to this study, the condition of Yellow Croaker declined in recent years and was skinny in the south, with worse condition than in the north. The values of condition factor a in LWR decreased, when the environmental pollution increased. Therefore, LWR could be a rough indicator for the environment quality. The mixed effects model provided a more convenient and reasonable approach to estimate the spatial and temporal variations of LWR for Yellow Croaker.