A Revision of the Traditional Analysis Method of Allometry to Allow Extension of the Normality-Borne Complexity of Error Structure: Examining the Adequacy of a Normal-Mixture Distribution-Driven Error Term

Huxley's model of simple allometry provides a parsimonious scheme for examining scaling relationships in scientific research, resource management, and species conservation endeavors. Factors including biological error, analysis method, sample size, and overall data quality can undermine the reliability of a fit of Huxley's model. Customary amendments enhance the complexity of the power function-conveyed systematic term while keeping the usual normality-borne error structure. The resulting protocols bear multiple-parameter complex allometry forms that could pose interpretative shortcomings and parameter estimation difficulties, and even being empirically pertinent, they could potentially bear overfitting. A subsequent heavy-tailed Q-Q normal spread often remains undetected since the adequacy of a normally distributed error term remains unexplored. Previously, we promoted the advantages of keeping Huxley's model-driven systematic part while switching to a logistically distributed error term to improve fit quality. Here, we analyzed eelgrass leaf biomass and area data exhibiting a marked size-related heterogeneity, perhaps explaining a lack of systematization at data gathering. Overdispersion precluded adequacy of the logistically adapted protocol, thereby suggesting processing data through a median absolute deviation scheme aimed to remove unduly replicates. Nevertheless, achieving regularity to Huxley's power function-like trend required the removal of many replicates, thereby questioning the integrity of a data cleaning approach. But, we managed to adapt the complexity of the error term to reliably identify Huxley's model-like systematic part masked by variability in data. Achieving this relied on an error term conforming to a normal mixture distribution which successfully managed overdispersion in data. Compared to normal-complex allometry and data cleaning composites present arrangement delivered a coherent Q-Q normal mixture spread and a remarkable reproducibility strength of derived proxies. By keeping the analysis within Huxley's original theory, the present approach enables substantiating nondestructive allometric proxies aimed at eelgrass conservation. The viewpoint endorsed here could also make data cleaning unnecessary.


Introduction
Julian Huxley envisioned the notion of constant relative growth between the size of a trait y and overall body size x [1][2][3]. Concurring formulation ordinarily referred to as Huxley's model of simple allometry is expressed through the power-function law: The power function form bearing Huxley's model sustains both theoretical and empirical approaches in many research endeavors, e.g., biology [4][5][6], physics [7], economy [8], earth and atmospheric sciences [9], ecology [10], and resource management [11]. Allometric methods are particularly relevant in seagrass research. Seagrass species provide valuable ecological services in estuaries and nearshore environments, for instance, by offering food and shelter for a myriad of ecologically and economically valued marine organisms [12][13][14], contributing to nutrient cycling [15,16], favoring the stabilization of the shoreline as roots and rhizomes compact the substrate, preventing erosion [17,18], participating in the foundation of the detrital food web [19], and also playing a fundamental role in carbon sequestration [20]. In seagrass research, allometric methods mainly aim at predicting response to changing environmental conditions or analyzing growth patterns, for example, the relationship between the width of the leaves and their dry weight [21], the relationship between the length of stems and their density [22], and the relationship between the size of the leaves and their dry weight [23].
Zostera marina L. also known as eelgrass is an essential seagrass species providing vital ecological services in estuaries and nearshore environments. In addition to the aforelisted seagrass benefits, eelgrass offers a nursery for waterfowl and fish species and nutrient recycling. But despite the ecological relevance of eelgrass meadows, deleterious anthropogenic influences currently threaten their permanence [24]. Eelgrass remediation efforts mainly rely on transplanting endeavors [25]. Assessing the success of concurring plots depends on nondestructive estimations of standing stock from which total leaf biomass is an important constituent. When Huxley's model produces a reliable fit to an eelgrass leaf biomass and area data set, it could provide reliable surrogates of eelgrass leaf biomass based on direct nondestructive measurements of leaf area. Conceiving present allometric examination methods aim to enhance the efficiency of Huxley's model-based constructs for eelgrass conservation.
Despite the pertinence of Huxley's model, some factors limit the accuracy of deriving projections. Firstly, a response variable y expressed as a function of its covariate x through Equation (1) is extremely sensitive to the variation of estimates of the parameters α and β [26]. Then, error propagation could undermine the precision of Huxley's model-based projections of response values. Other prime influencers on the accuracy of estimates of parameters in Huxley's model are the analysis method, sample size, and data quality [27][28][29]. Mainly, envisioning a suitable analysis method to get estimates of parameters in Huxley's model relies primarily on detecting the implicit variation pattern in the original data scales. In some settings, the data spread displays a heteroscedastic pattern characterized as an increasing variation in the response's replicates concerning or relative to the predictive variable. The traditional approach to assembling a regression scheme involves Huxley's power function form as the systematic part and a multiplicative error term specified as a lognormally distributed random variable. The resultant method entails parameter identification through nonlinear regression in the direct scales of data. Concurrent to this approach concerns contemplating a logarithmic transformation that allows the analysis transference from natural scales into geometrical space. Their examination involves a regression model including a systematic linear part and an additive and normally distributed error term. The method completes by performing a back transformation step determining the identified form of Huxley's model of simple allometry in the direct scales. The last phase requires using a factor seeking to correct the bias of retransformation [30][31][32]. We will refer to this log transformation approach as the traditional analysis method of allometry (TAMA). Oppositely to a heteroscedastic spread pattern, data in the original scales could adapt to a homoscedastic dispersion outline. Accordingly, it may be pertinent to assume the appropriateness of a regression model composing Huxley's power function form at the systematic part along with an additive and normally distributed error term. Such a scheme precludes a log transformation step and sets nonlinear regression in the direct data scales as a necessary protocol for parameter identification tasks. So built regression model refers here as a direct nonlinear regression scheme (DNLR).
In allometric examination, when addressing the traditional multiplicative, log transformation, or direct nonlinear regression schemes, the essential task should simultaneously identify the systematic part and the error term in the associating regression protocol. The systematic part and the error term are of chief importance since the first determines the trend and the second, the dispersion pattern of the data. Nevertheless, usually, the virtual run concerns the identification of the systematic component of the model. Such a drive embodies attempts to improve goodness-of-fit by modifying Huxley's systematic part while keeping a normality-borne shape of the error term. Such a move conceives constructs referred to as complex allometry forms [33][34][35]. Moreover, commonly, concurrent examination attempts to carry out the normality of errors without verifying the assumptions of this model. Such a tactic somehow averts due attention to questioning the adequacy of the complexity of the assumed distribution of the error term. But Montesinos-López et al. [36] stirred away from this practice. They dealt with a data set including pairs of measurements of eelgrass leaf biomass and area. A TAMA fit resulted inconsistent. Montesinos-López et al.'s [36] amendment conformed to a systematic linear part but switched to a nonnormal distribution-brought error term. Moreover, Montesinos-López et al. [36] found that clinging to a logistic distribution-based try, not only lead to a consistent residual spread but also, remarkably improved the reproducibility strength of proxies for average leaf biomass derived from allometric projections established by a TAMA approach.
Here, we analyzed eelgrass leaf biomass and area data adapted from a sample originally reported in Echavarria-Heras et al. [37]. Compared to [36], present data exhibited overdispersion and marked size-related heterogeneity, perhaps explained by a lack of systematization in data gathering. Linked spread pattern precluded consistency of a TAMA fit. Moreover, the modification based on a logistic error term addressed in [36], fitted hereafter data cleaning procedures, also failed to produce consistent results. A normal mixture distribution provides a convenient model for data sets displaying high variability and heterogeneity. Accordingly, we conjectured that while analyzing present data, keeping the systematic linear term inherent to the TAMA approach but bearing a mixture of two normal distributions as a candidate for error term distribution could be suitable. Compared to composites including multiple parameter-complex allometry forms, normality of errors, and data cleaning procedures, the consistency of the present arrangement fitted in raw data delivered a notable reproducibility strength of proxies of eelgrass leaf biomass. Identifying the referred allometric relationship upon Huxley's power function model framework avoids complications from complex allometry forms during verifying parameter invariance. This feature is crucial for genuinely nondestructive assessments since previously fitted parameters could be used to get allometric projections of eelgrass leaf biomass values.
Because of difficulties tied to the implementation of complex allometry constructs, the present findings certainly enhance the perception of looking for the appropriateness of the error distribution as a mechanism to achieve a better fit of Huxley's model as suggested by Montesinos-López et al. [36]. Present findings exhibit the strength of a normal mixture distribution-borne error term as a device to produce a consistent fit of Huxley's model in a scenario of marked variability in data. It is also worth emphasizing that by keeping the analysis within the confines of Huxley's original theory, the current approach enables substantiating nondestructive allometric proxies aimed at eelgrass conservation. Besides, the scheme endorsed here could also make data cleaning unnecessary. And, since we provide a detailed explanation of the implementation, including mathematical, statistical, and computational aspects, the offered scheme can be straightforwardly adapted to other allometric examination endeavors. Also, our mixture distribution assumption on driving the error term bears a path yet not undertaken within the traditional analytical assortment of allometry. Therefore, we considered it worth reporting its suitability, and this manuscript devotes itself to that aim.
This work is structured as follows: In Section 2, we explain the formalities of the basic regression schemes deriving from Equation (1) and that circumscribe to a normal distribution-borne error term. For comparison aims in Section 3, we include the spread plots of the present data and that analyzed by Montesinos-López et al. [36]. Section 3 also elaborates on the necessary modifications of the basic regression schemes introduced in Section 2 that allow consid-eration of error terms driven by either a logistic or normal mixture distributions. Section 3 presents the results of conceived regression protocols fitted to the present data and compares their reproducibility power. The results section also incorporates a simulation study aimed to establish the strength of the approach under a known scenario. Section 4 pertains to the discussion that stresses the strengths and weaknesses of the current approach. We strain on that by keeping allometric examination within the confines of Huxley's original theory; the present process bears advantages while substantiating nondestructive assessment proxies aimed at eelgrass conservation. Section 5 presents the conclusions of this study and suggests future work. Appendix A elaborates on reproducibility measures to assess the suitability of the allometric projection methods offered here. Appendix B presents the formalities of the AIC index-based comparison of models fitted on different scales.

Huxley's Multiplicative Error
Model. There are settings in allometric analysis with a spread in the original scales of data displaying a pattern of increasing variation in the response concerning itself or relative to the predictive variable. Assembling a candidate regression scheme usually undertakes a Huxley's multiplicative error model. It involves a systematic part acquiring a power function form and an error term specified as a function δðϵÞ of a random variable ϵ that acts in a multiplicative way, namely, where y stands for the response variable, x for the covariate, α and β are parameters, δðϵÞ = exp ðϵÞ and ϵ taken as a normally distributed random variable having zero mean and deviation σ; that is, ϵ~Nð0, σÞ [38][39][40]. Proposed form sets δðϵÞ as a lognormally distributed random variable with zero log-mean and log-deviation σ, that is, δðϵÞ~lognormð0, σÞ: The likelihood function takes the form where μ i = log ðβÞ + α log ða i Þ. Besides, as it is set by Equation (2), the variability of the response at a given value of the covariate regulates by the contribution of the random error δðϵÞ and the value of the systematic part βx α . So, for large covariate values, the named range of variation grows, thus resulting in a heteroscedastic statistical model. For present aims, the model set by Equation (2) refers as Huxley's multiplicative error model (HMEM) or simply as a MEM protocol.
Assuming suitability of the model of Equation (2), we usually address the problem of acquiring the mean of the response y conditioned on a covariate value x. Associating form denoting here through EðyjxÞ is gotten by taking the 3 BioMed Research International expected value on both sides of Equation (2) conditioned by the explanatory variable x, namely, Then, to obtain the mean response EðyjxÞ, in addition to the power function-like systematic term βx α , we must consider a factor EðδÞ, which interprets as a correction factor (CF) for bias of allometric projection of the mean response through the estimated form of the power function βx α . Only in the case of Huxley's lognormal multiplicative error model the correction factor EðδÞ takes on the form given by Equation (4). In the general settings given the distribution of the random variable δðϵÞ we could attempt to obtain a closed-form for the correction factor EðδÞ by evaluating the expectation of the response variable y. It is worth emphasizing that in getting a closed form for EðδÞ, it is essential to identify the form of the distribution acquired by the error term δðϵÞ.
2.2. The Traditional Analysis Method of Allometry. Concurrent to Huxley's lognormal multiplicative error model of Equation (2), there is an approach relying on a log transformation: ðx, yÞ ⟶ ðu, vÞ = ðlnx, lnyÞ, that allows the contemplation of a linear regression model in the geometrical scales, namely, where β 0 = lnβ and with an additive error term ϵ = ln ðδÞ expressing as a normally distributed random variable having zero mean and deviation σ, that is, ϵ~Nð0, σÞ. The likelihood function is as follows: and with μ i = β 0 + αu i . Based on the identified form of Equation (6), we perform a back transformation step to get estimated form of the mean response function of Equation (4) and that of the correction factor (5) [30][31][32]. Afterwards, we refer to the protocol of Equation (4) as the traditional analysis method of allometry (TAMA) [41][42][43].
2.3. The Direct Nonlinear Regression Protocol. Oppositely to a circumstance described by the Huxley's lognormal multiplicative error model of Equation (2), it may be pertinent to conceive a regression model where the systematic power function-like term maintains, but that the random error contributes additively to the variability of the response y, that is, with ϵ usually assumed as a normally distributed random variable having a zero mean and a deviation σ, that is, ϵÑ ð0, σÞ. Therefore, oppositely to the heteroscedastic spread entailed by the multiplicative error model of Equation (2) for the additive error model of Equation (6), the contribution of ϵ to y variability is its value itself, being this null when ϵ vanishes. The likelihood function turns out to be with μ i = βx α i . We further on refer to the scheme of Equation (6) as direct nonlinear regression (DNLR) [40,[44][45][46].
In what follows, we will refer generically to the error structure of a given regression scheme as a merge of the way the error term enters into the model and the distribution that drives its stochasticity. Accordingly, we say that the model of Equation (2) bears a multiplicative-lognormal error structure and that the models of Equation (6) and Equation (8) both share an additive-normal error structure. For the aim of exploring the extent of modifying the error structure of the MEM, TAMA, or DNLR schemes in what follows, we conceive composite regression schemes that maintain the involved systematic terms but modify the assumption on the error-shaping random variable ϵ from normality to being logistically distributed or else, according to a mixture of two normal distributions of common zero mean but different deviations. Particularly, for the DNLR scheme we adapt a Breusch-Pagan [47] type variance function form, aimed to take over heteroscedasticity. For comparison, we include the polynomial modification to TAMA's scheme undertaken by Echavarría-Heras et al. [37]. Formal expressions of the composite regression schemes addressed here appear in the results section.

Results
3.1. Data. Present examination relies in a data set comprising pairs of measurements of leaf weight y ½g and relating area x ½mm 2 , adapted from a sample reported in Echavarría Heras et al. [37] obtained by a 13-month sampling performed on an eelgrass meadow in Ensenada, B.C., Mexico. Figure 1(a) displays data spread in the original arithmetical scales. We can be aware of noticeable variability of replicates as well as of marked heterogeneity of patterns among the pools of smaller and larger leaf area values in the sample. Figure 1(b) pertains to spread corresponding to log scales. To compose present data set we removed two of the 10412 pairs reported. We further refer to the resulting 10,410 pairs as the present data set that could be also indistinctly referred by means of the symbol EHDS, for Echavarría-Heras et al. [37] data set.
For comparison aims, we depend on a second eelgrass leaf biomass to area data set examined by Montesinos-López et al. [36] and collected at the same meadow as the EHDS. It composes a total of 537 pairs of measurements of leaf weight ðyÞ and relating area ðxÞ, also gathered in the Echavarria-Heras et al. [37] study site but only over a onemonth sampling. Ahead this data set identifies by the label MLDS, for Montesinos-López et al. [36] data set. Figure 2(a) displays associating spread in the original arithmetical scales. Figure 2(b) presents that one corresponding to geometrical scales.

3.2.
Performance of the MEM and TAMA Schemes. The heteroscedastic pattern shown in Figure 1(a) suggests exploring the suitability of a MEM assumption. Therefore, we could firstly analyze the EHDS according to the regression model of Equations (2) and (3) or equivalently by calling on to the concurring TAMA scheme in geometrical scales appointed by Equations (6) and (7). Including this last fit provides insight at envisioning the actual distribution of the error term contemplated in the regular MEM scheme. Also incorporating the TAMA-based fits allows comparing present findings to results reported by Echavarría-Heras et al. [37], as well as to those by Montesinos-López et al. [36] when analyzing the MLDS. Table 1 presents fitting statistics of a MEM try, and those relating to a TAMA fit do in Table 2. Spread plots on both the MEM and the TAMA fits appear in Figures 3 and 4 one to one. Figure 3(a) displays the spread about the fitted MEM's power function-like systematic part. Figure 3(b) suggests that the error term does not match the expected lognormal distribution pattern. Besides, Figure 3(c) shows the associating Q-Q lognormal plot of the residuals. Vertical lines sketched in Figure 3(c) delimit the linear part of the Q-Q lognormal plot. Such a sector places in the interval (0.240, 3.883), between the ordered observations, num. 700 and the 9500 of the 10410 original data pairs. Lowermost panels in Figure 3 display a close-up split of the Q-Q lognormal diagram in Figure 3 (6) and (7) fitted to the EHDS in geometrical scales. Moreover, Figure 4(b) already suggest that the error term does not match the expected normal distribution pattern. Besides, opposing to a masked heavy tail at the left extreme of Figure 3(c) the spread in Figure 4(c) clearly reveals that while the normal distribution fits adequately in the central part, the pattern at the extremes departs from that corresponding to a normal distribution; namely, we have heavy tails in TAMA's set of residuals. Concerning this fit, calculating the value of the kurtosis coefficient provides additional evidence of the existence of heavy tails in the distribution of the error term. Certainly, the associating kurtosis coefficient attained a value of kurt = 14,176, which is a much larger value than the one corresponding to a normal distribution ðkurt = 3:0Þ:According to Wheeler [48], such a big value indicates that the distribution of the residues has heavy tails since kurtosis provides information on the extremes rather than the central part of the distribution. Therefore, a TAMA scheme turns out to be very simple to produce a reliable fit on the EHDS. Indeed, the analyzed data exhibit a variability pattern suggesting that the complexity of an error structure beyond that one bearing to normality turns out to be necessary to grant a coherent fit. Thus, the resulting residuals should instead model through a distribution with a considerably greater overdispersion than the normal one in the geometric scales could explain. In summary, results suggest that a normal distribution assumption for the random variable ϵ inherent to the MEM and TAMA fits does not support a suitable model for the variability pattern inherent to the EHDS.
3.3. Implementation of the Breusch-Pagan Modification on the Regular DNLR Protocol. As it conceives here, a DLNR-BP protocol stands for a Breusch-Pagan [47] modification of the basic DLNR scheme of Equation (6) envisioned to account for the heteroscedastic pattern shown in Figure 1(a). Formally, a DNLR-BP adaptation acquires a form: with ϵ taken as normally distributed random variable, having a zero mean, and a covariate dependent deviation σðxÞ, that is, ϵ~Nð0, σðxÞÞ. To offer a suitable candidate form for σðxÞ, we recall the procedure yielding the Breusch-Pagan [47] test, so we set the following: Additionally, the likelihood function becomes the following: with μ i = βx α i and σðy i jx i Þ = σð1 + kx i Þ.  Table 2: Estimated parameters and related statistics associating with a TAMA scheme as given by Equations (6) and (7) and fitted on the EHDS. We include parameter estimates with their projected uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one. LogLikMxA and AICA correspond one to one to values of the LogLikMx and AIC statistics expressed in direct arithmetical scales.  Table 3 presents parameter estimates and associating fitting statistics of the DNLR-BP scheme of Equations (10) and (12) performing on the EHDS. The resulting spread plots appear in Figure 5. In Figure 5(a), we exhibit the dispersion about fitted Huxley's power function form. Figure 5(b) shows the scatter diagram of residuals against leaf area. In turn, Figure 5(c) shows a Q-Q normal plot of the residuals on direct arithmetical scales. Although spreads in (a) and (b) suggest the presence of two correspondence rules conforming to the global area-weight relationship, the fitted power function deems roughly consistent. But still, Figure 5(c) displays a heavy tails pattern of the error term, although now being asymmetrical. Such a Q-Q normal diagram perhaps explains by the fact that even by appointing a Breusch-Pagan [47] variance function form, that move failed to provide a sound model for the heterogeneity of the heteroscedastic spread. But said asymmetrical spread could additionally explain by the Lai et al. [40] observation that direct nonlinear regression can produce bias for large covariate values. Moreover, comparison of AIC values among the MEM (AIC = −94745:76) and the DNLR-BP ð AIC = −84528:9Þ fit produces ΔAIC = −10217, favoring the MEM fit. Therefore, a multiplicative error structure seems more fit at analyzing the EHDS. The relatively fair spread about the fitted systematic term hints on adequacy of Huxley's model at describing the trend in the EHDS but the spreads in the Q-Q normal diagrams accompanying the MEM, TAMA, and DNLR-BP fits suggest that assuming the error shaping random variable ϵ as being normally distributed is statistically unsupported.   [36], while analyzing the EHDS, we attempted to model overdispersion or residuals produced by the regular MEM or TAMA fits by agreeing to a complexity of error structure adaptation approach. We assumed first that the basic random variable ϵ better conforms to a logistic distribution, since this last one works fine when the overdispersion is low. Resulting analytical arrangements describe by the acronyms MEM-loglogistic and TAMA-logistic. Formally, the MEM-loglogistic construct stands for a modification of the regression model of Equation (2), which establishes through with δðϵÞ = exp ðϵÞ and ϵ taken as a logistically distributed random variable having zero location and scale σ, that is,  Figure 4: Dispersions about the fitted systematic part, residuals, and Q-Q diagram of a TAMA scheme as given by Equations (6) and (7) and fitted on the Echavarría-Heras et al. [37] data.   (10) and (12)) fitted on the EHDS. We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one.
being μ i = βa α i . The mean response function becomes the following: where Meanwhile, according to Equation (5), a TAMA-logistic scheme sets by the model with ϵ~logisticð0, σÞ: The likelihood function takes on a form: where μ i = β 0 + αu i . The back transformation step requires using the correction factor Estimated parameter values and relating statistics for a MEM-loglogistic protocol of Equations (13) and (14) fitted on the EHDS appear in Table 4. Correspondingly, those pertaining to a TAMA-logistic fit (cf. Equations (17) and (18)) display in Table 5. In Figure 6, we include spread plots on both the MEM-loglogistic and the TAMA-logistic fits. Besides, upper panels of Figure 6 present spreads associating with MEM-loglogistic fit, that is, Figure 6(a) includes spread about fitted Huxley's power function, Figure 6(b) corresponding to residual's dispersion and Figure 6(c) devoting to the Q-Q loglogistic diagram spread of residuals in direct arithmetical scales. As we arranged for Figure 3(c), by splitting the diagram in Figure 6(c) and then amplifying we can also be aware of that the pattern at the extremes departs from that corresponding to a loglogistic distribution; namely, we have heavy tails in the set of residuals corresponding to the MEM-loglogistic fit. Correspondingly, lower panels in Figure 6 present dispersion patterns resulting from a TAMA-logistic fit. Figure 6(d) portraits dispersion about fitted TAMA-logistic line. Figure 6(e) exhibits residual dispersion led by the TAMA-logistic fit.
And Figure 6(f) portraits the corresponding Q-Q logistic plot of residuals in geometrical space. Again, a heavy tails pattern shows. For the sake of conciseness when comparing to the Montesinos-López et al. [36] fit, we only discuss the implications of the TAMA-logistic spread plots. Compared to a regular TAMA fit present residual spread still shows an uneven pattern. Moreover, the improvement in consistency Table 4: Estimated parameters, related statistics, and AIC values associating with a MEM-loglogistic scheme (cf. Equations (13) and (14)) fitted on the EHDS. We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one.   (17) and (18)) fitted on the EHDS. We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one. LogLikMxA and AICA correspond one to one to values of the LogLikMx and AIC statistics expressed in direct arithmetical scales.   (13) and (14)) and a TAMA-logistic (cf. Equations (17) and (18)  The logistic distribution can reasonably model a higher dispersion than the normal one [36]. But the heavy tails arrangement displaying in the Q-Q loglogistic diagram of Figure 6(c) refers to an unfair MEM-loglogistic fit (Equation (9)), and equally, the Q-Q logistic plot of the residuals of a TAMA-logistic fit (Figure 6(f)) on geometrical space showing heavy tails also point to an inconsistent TAMA-logistic (Equation (10)) fit.

Performing the Breusch-Pagan and Logistic Distributed
Error Modification on the DNLR Protocol. As it conceives here, a DLNR-BP-logistic protocol stands for a Breusch-Pagan [47] along with a logistically distributed error modification of the basic DLNR scheme of Equation (6) envisioned to account for the heteroscedastic pattern shown in Figure 1(a). Formally, a DNLR-BP-logistic adaptation admits a form: with ϵ taken as a logistically distributed random variable, having a zero location and a covariate dependent scale σðxÞ, that is, ϵ~logisticð0, σðxÞÞ, where we set the following: The likelihood function becomes with μ i = βx α i , and σ i = σð1 + kx i Þ. Table 6 presents parameter estimates and associating fitting statistics of the DNLR-BP-logistic scheme as given by Equations (20) and (22) Figure 7(c) displays a heavy tails pattern of the error term, once again being asymmetrical as in the case of a DNLR-BP fit. Again, such an asymmetrical spread could explain by the embedding of a real multiplicative error structure into an additive error counterpart. An on top of that, the Lai et al. [40] observation of biased direct nonlinear regression output for large covariate values could be also pertinent. In any event, the spread in the Q-Q logistic in Figure 7(c) bears that assuming logistically distributed residuals turns out to be unfeasible. At this point of the matter, by looking at a relatively fair spread about the systematic parts of fitted MEM, TAMA, DNLR-BP schemes and their subsequent modifications to contemplate a logistic distribution-borne error term, it is not idle to say that the detected inconsistencies in Q-Q diagram spreads, could already hint at a lack of suitability of either a normal or a logistic distribution as compatible models for the residual dispersion in present data.
3.6. Implementation of the TAMA Scheme Fitted on Processed Data. Heterogeneity of spread in Figure 1 could perhaps explain by the participation of multiple agents that contributed to data gathering without standardization of routines [37]. Then, at first glance, proliferation of unduly replicates could generate issues at data quality that elucidate the lack of fit of schemes based on Huxley's model and accompanying normal or logistic distributions-borne error terms. Echavarria-Heras et al. [28] adapted a median absolute deviation procedure to remove anomalous replicates in a sample of similar allometric eelgrass data. We engaged these procedures on present data to explore the extent of data quality influences in determining the inconsistencies of the fits above. Table 7 presents parameter estimates and Table 6: Estimated parameters, related statistics, and AIC values associating with a DNLR-BP-logistic scheme of Equations (20) and (22) fitted on the EHDS. We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one.  (6) and (7) and fitted on the processed EHDS. These statistics on the processed data are only provided for completion of the presentation and do not intend to support comparison to the remaining fits as those relied upon the crude EHDS. Figure 8 displays the spread plots associated with a TAMA fit on processed data. Figure 8(a) displays dispersion about fitted Huxley's linear systematic part. Figure 8(b) shows the scatter diagram of residuals against leaf area. Figure 8(c) presents the Q-Q normal plot of the residuals that shows heavier tails than those expected for a normal distribution. Even though the addressed median absolute deviation procedure removed a large share (25%) of the original data, the TAMA fit could not normalize the residual dispersion. Moreover, differentiation of spread patterns among the pools' smaller and larger leaf sizes barely hinted by plots in Figure 1 seems to portray undoubtedly once data processing completes, as shown by (a) and (b) in Figure 8. Therefore, considering this and for improving the quality of fits of schemes involving Huxley's model-driven systematic term, it seems reasonably calling in distributions that allow a greater dispersion than the normal or the logistic ones at assembling the accompanying error term.

3.7.
Assessments of the Normal Mixture Amendments to the Regular MEM and TAMA Schemes. Constituents in the family of finite mixtures of distributions are highly flexible due to the diversity of forms that they can acquire. Particularly, the distributions that we can construct through finite mixtures of normal distributions are very varied since we can obtain multimodal, skewed, and distributions with excess kurtosis. Particularly, at adapting the symmetric distribution with zero mean and heavy tails suggested by the spread in Figure 4(c), we could explore the appropriateness of a mixture of two individual normal distributions; N 1 ð0, σ 1 Þ and N 2 ð0, σ 2 Þ having a common zero mean but different deviations σ 1 and σ 2 .We assume also, that the weight through which N 1 ð0, σ 1 Þ and N 2 ð0, σ 2 Þ participate in the mixture designates by means of p: For simplicity further on we employ the symbol M2N to designate so conceived mixture distribution. Assume that the error shaping random variable ϵ in Equation (9)   Heavy tails appear through a markedly asymmetrical pattern which refers to an overall inconsistent DNLR-BP-logistic fit.

BioMed Research International
Then, denoting associating probability density function through the symbol f ðϵ ; p, σ 1 , σ 2 Þ, we have the following: where f 1 ðϵ ; 0, σ 1 Þ and f 2 ðϵ ; 0, σ 2 Þ stand for one to one the density functions of the components N 1 ð0, σ 1 Þ and N 2 ð0, σ 2 Þ. It follows, directly from Equation (14) that EðϵÞ = 0 and that joining variance σ 2 ðϵÞ takes on the form: Therefore, the assumption of ϵ being distributed according to presently conceived normal mixture M2N symbolizes through ϵ~M2Nðp, 0, σÞ where σ derives from Equation (16). Whenever ϵ~M2Nðp, 0, σÞ, the random variable δðϵÞ = exp ðϵÞ acquires a log-mixture of two normal distribution of zero log-mean, log-deviation σ and weight p, that is, δðϵÞ LM2Nðp, 0, σÞ. Moreover, the associating density function denoted through f δ ðδ ; 0, p, σ 1, σ 2 Þ is given by the following: where f 1 ðδ ; 0, σ 1 Þ and f 2 ðδ ; 0, σ 2 Þ are the two lognormal density functions, with common log-mean zero and logdeviations σ 1 and σ 2 , respectively. Correspondingly, the adaptation of the MEM protocol that replaces a lognormally distributed error term by a LM2N alternate denotes by MEM-LogM2N. Formally, such a scheme takes on a form: with δðϵÞ = exp ðϵÞ and ϵ taken as a M2N distributed random variable, that is, ϵ~M2Nðp, 0, σÞ. Therefore, δðϵÞL ogM2Nð0, σÞ. The likelihood function becomes the following: where μ i = log ðβÞ + α log ðx i Þ. And the mean response function takes on a form: Then, Equation (18) implies the following: Besides, we use a composite TAMA-M2N to distinguish the regression arrangement that modifies the regular TAMA scheme to consider a M2N distributed error term. Formally, the scheme stands for a modification of the regression model of Equation (5) that establishes through the following: with ϵ being a M2N distributed random variable with zero mean and deviation σ, that is, ϵ~M2Nð 0, σÞ. The corresponding likelihood function is given by the following: where μ i = β 0 + αu i . Table 7: Estimated parameters and related statistics associating with a TAMA scheme as given by Equations (6) and (7) and fitted on the processed EHDS. We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one. LogLikMxA and AICA correspond one to one to values of the LogLikMx and AIC statistics expressed in direct arithmetical scales. Statistics on the processed data are only provided for completion of the presentation and do not intend to support comparison to the remaining fits as those relied upon the crude EHDS. Correspondingly, according to Equation (19), the correction factor for retransformation aims becomes the following: Table 8 presents fitting statistics relating to a MEM-LogM2N protocol (cf. Equations (26) and (27)) performing on the EHDS. Table 9 provides fitting statistics associating with a TAMA-M2N fit (cf. Equations (30) and (31)). Figure 9(a) presents the spread about the MEM-LogM2N fitted Huxley's form of a systematic part. Figure 9(b) exhibits the corresponding residual spread. Figure 9 .We may be also aware of improved residual spread relative to the severe heavy tail patterns underlying both the MEM and the MEM-loglogistic schemes. Figure 9(d) shows the spread about the linear mean response in geometrical space produced by a TAMA-M2N scheme fitted on the EHDS. Figure 9(e) displays TAMA-M2N residual spread and Figure 9(f) presents the link up with Q-Q M2N plot. The M2N modified TAMA scheme shows a better fit than either the regular TAMA or its logistic-borne counterpart since as we can be ascertained by comparing Q-Q diagrams in Figure 3 through Figure 7, the consistencies of the Q-Q M2N-borne spreads are remarkable, compared to the severe heavy tails patterns displayed by either the Q-Q normal and Q-Q logistic plots.
Besides when comparing the fits with error term produced by the M2N to those driven by the normal or logistic distributions on the arithmetic scale, using the Akaike information criterion (AIC), where the lower value of the AIC index represents a better fit, we can be aware from with ϵ~M2Nðp, 0, σÞ. For a M2N distribution, the standard deviation becomes the following: In order to adapt the model of Equation (19) for the heteroskedastic case, similarly to the Breusch-Pagan [47] test, we modify the standard deviations σ 1 and σ 2 , such that one to one take forms σ 1 ð1 + kxÞ and σ 2 ð1 + kxÞ, to obtain the following: or equivalently The likelihood function of this model can be written as follows: where μ i = βx α i and σ ji = σ j ð1 + kx i Þ, j = 1, 2.  (26) and (27)) fitted on the EHDS We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one.   Table 10 presents the parameter estimates and fitting statistics relating to a DNLR-BP-M2N scheme (cf. Equations (33) and (37)) fitted on the EHDS. Figure 10 exhibits corresponding spread plots. Figure 10   function and that of residuals around the zero lines seem fair. Compared to the fit of a regular DNLR scheme of Equations (8) and (9), the Q-Q M2N diagram appearing in Figure 10(c) displays improvement regarding the expected pattern for a heavy tailed distribution. But anyhow, an asymmetrical heavy tailed spread persists. Table 5 explains the  (33) and (37)). We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMx, and AIC stand for maximum loglikelihood and Akaike information index values one to one.   improvement of the quality of the fit in terms of the AIC value achieved by the DNLR-M2N fit. In this table, we note that there is a greater reduction when we change from the normal distribution borne DNLR to the DNLR-M2N composite; that is, we find an important change when going from the normal to the M2N distribution, this of course due to the heavy-tailed pattern that exists in the data.

Implementation of the TAMA-Poly(K): The Kth-Degree
Polynomial Modification to the Linear Systematic Term of a Regular TAMA Scheme. With the aim of comparing the reproducibility strength of a M2N-driven modification of a TAMA scheme against that entailed by the Kth-degree polynomial variation to a regular TAMA undertaken by Echavarría-Heras et al. [37]. So, conceived scheme refers through the symbol, TAMA-PolyðKÞ, and entails a form of a complex allometry and normally distributed error composite scheme in geometrical space. Formally, a TAMA-PolyðKÞ scheme stands for a modification of the regression model of Equation (2) that establishes through the following: with ϵ~Nð0, σÞ, and so, the likelihood function of this model is given as follows: where This time, the correction factor for retransformation aims becomes the following: Table 11 presents fitting statistics of a TAMA-Poly(6) scheme (cf. Equations (38) and (39)) performing on the EHDS. Figure 11(a) shows the spread about the systematic term, in geometrical space of a TAMA-Poly(6) scheme fitted on the EHDS. Figure 11(b) shows the scatter diagram of the TAMA-Poly(6) residuals against the logarithm of the area. Figure 11(c) displays the TAMA-Poly(6) associating Q-Q normal plot of the residuals on a geometric scale. Even though, fitted TAMA-Poly(6) systematic part seems to better describe a phase-like heterogeneity shown in Figure 1(b), still the Q-Q normal plot of Figure 11(c) displays a heavy tails spread pattern. AIC values in Table 8 and Table 10 allow comparison of the TAMA-M2M and the TAMA-Poly(6) models fitted on the EHDS. We can be aware of a difference ΔAIC = 3139 favoring selection of the TAMA-M2M over the TAMA-Poly(6) counterpart. This try conveniently explains the failure of a polynomial-like complex allometry form to normalize the distribution of residuals.

Assessment of Reproducibility Strength of Addressed
Methods. We present statistics that allow the assessment of the reproducibility strengths of adapted allometric proxies while projecting values for the assessment of reproducibility strength of allometric proxies for observed monthly average leaf biomass in shots reported in the EHDS (Table 12). We provide AIC index and other model performance metrics, such as the standard error of estimate (SEE) and mean prediction error (MPE) [49][50][51][52]. Likewise, we include Lin's concordance correlation coefficient, denoted here through the CCC symbol [53]. Agreement between allometrically Table 11: Fitting results of a TAMA-Poly(6) regression protocol based on a 6th degree polynomial systematic term modification to the regular TAMA scheme (cf. Equations (38) and (39)) performing on the Echavarría-Heras et al. [37] data set. We include parameter estimates with their estimated uncertainties, t value, p value, LogLikMax, and AIC stand for maximum loglikelihood and Akaike information index values one to one. LogLikMxA and AICA correspond one to one to values of the LogLikMx and AIC statistics expressed in direct arithmetical scales.   Figure 11: Spreads in geometrical space produced by the TAMA-Poly(6) and the polynomial modified TAMA scheme proposed by Echavarría-Heras et al. [37] when fitted on the EHDS (cf. Equations (38) and (39)). (a) Spread about fitted 6th degree polynomial systematic term, in geometrical space. (b) The scatter diagram of residuals against the logarithm of the area. (c) The associating Q-Q normal plot of the residuals on a geometric scale. Even though fitted systematic part seems to better describe a phase-like heterogeneity, still a heavy tails spread pattern in the Q-Q normal plot shows up. Table 12: Statistics for the assessment of the reproducibility strengths of the presently acquired allometric proxies while projecting monthly averages of observed eelgrass leaf biomass in shoots as reported in the EHDS. We include Akaike information criterion (AIC) index, Lin's concordance correlation coefficient ðCCCÞ, standard error of estimate (SEE), mean prediction error (MPE), and relative deviation (RD). Comparison restricts to better performing proxies identified in arithmetical scales as well as the retransformed TAMA-Poly(6) scheme. between the overall mean of monthly averages of observed leaf biomass values and the corresponding one derived allometrically, and then, we divide by the overall mean of monthly averages of observed individual leaf biomass values.

Model
According to Tables 1 and 3 for the normal distribution produced models, it was the MEM one which performed the better. Nevertheless, a MEM row in Table 12 allows assessing the feasibility of assuming that the basic error random variable ϵ distributes normally. As we can be aware of, the MEM row implies the less accurate proxies for monthly average leaf biomass in shoots .Certainly among the compared methods, the MEM one arrangement returned the smaller CCC entry and the greater SSE, MPE, MPSE, and RD values. We can also ascertain that enhancing complexity of error structure to contemplate a M2N distribution led to a remarkable reproducibility strength of MEM-LogM2N. Moreover, a difference between a CCC value of the TAMA-Poly (6)  3.11. Simulation Study. In this section, we carry out a study aimed to evaluate the performance of the proposed protocol under a known scenario. For this aim, we adapt a simulation procedure based on the multiplicative allometric model with a mean response of the Huxley power function type and error term according to a logM2Nðp, 0, σÞ distribution. The joining formulation is given by Equation (26). The procedure yields pairs including the foliar area values x i in the EHDS and simulated values of their corresponding leaf biomass replicates y ij , j = 1, 2, ⋯, rðiÞ. To produce the leaf biomass replicates, we use the formula: where β and α with are the values of the allometric parameters given in Table 8 and for a fixed x i i = 1, 2, ⋯, m , and for j = 1, 2, ⋯, rðiÞ,ϵ ij drawn from a M2Nðp, 0, σÞ distribution with σ = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pσ 2 1 + ð1 − pÞσ 2 2 p for p , σ 1 and σ 2 as listed in Table 8. The EHDS compose n = 10450 pairs of measurements of leaf area and related weight. The number of different leaf area values recorded in the sample is m = 2295 . Therefore, the vector ðrð1Þ, rð2Þ, ⋯, rðmÞÞ satisfies ∑ m 1 rðiÞ = n.
Once we produced the n simulated data pairs, we generated samples obtained by randomly removing ten pairs out of the simulated pool. We returned the extracted data so we could deal with the complete collection of simulated data before getting the following sample. The process iterated until the completion of a number of 1000 random samples. For the ith sample, we achieved a fit of the model of Equations (26) and (27) to obtain a vector ðβ si , α si , p si , σ 1si , σ 1si Þ of estimated parameters. Then, we acquired the root mean squared errors (RMSE) of so estimated parameters and their reference values ðβ, α, p, σ 1 , σ 2 Þ listed in Table 8. Calculation of RMSE values achieved through the formula: where θ stands for the value of a reference parameter and θ s for its estimated value through the simulation procedure. Correspondingly θ si stands for the proxy of θ s deriving from the associating ith sample of data pairs drawn from the simulated pool. We present the statistics resulting from simulation runs in Table 13. Figure 12 allows visual assessment of reproducibility features of the mean response function fitted from simulated data. Table 14 allows a comparison of referred reproducibility strength in quantitative grounds. Entries in Table 13 demonstrate that reference values of the allometric parameters in the model of Equations (25) and ((26)) can be efficiently approximated through the addressed simulation procedure. Histograms in Figure 2(a) and Figure 2(b) suggest consistency of an M2N distribution for both the EHDS and the simulated data set. We can also be aware from Figure 12(c) and Figure 12(d) that the simulation procedure failed to reproduce the domed-like spread at the beginning of the distribution in the EHDS. This fact readily explains the drop in AIC deriving from simulated data. It also strengthens a perception of a lack of standardization in data gathering while conforming to the EHDS sample. Spreads around the mean response function fitted from simulated data (blue lines) show fairness for both the EHDS and the simulated data. Whatever bias among projections produced by mean response curves explains by the lack of the referred domed-like spread in simulated data. The closeness of reproducibility index values in Table 14 corroborates this judgment. Altogether, the result of addressed simulation runs allows confirming the fairness of the offered M2N-based allometric analysis protocol.
3.12. Summary of Results. We extended the approach by Montesinos-López et al. [36] to consider a setup where a mixture of normal distributions drives stochasticity of the error term of the regular TAMA fit. Present data composes (n = 10,410) pairs of measurements of eelgrass leaf biomass and corresponding area adapted from a sample originally reported in Echavarria-Heras et al. [37] by removing two data pairs that we considered as unduly outliers (Figure 1). Moreover, at a first glance, we can be aware that present data stands for bigger an extent and variability than contemplated by Montesinos-López et al. [36] (Figure 2).
By looking at a relatively fair spread about the systematic parts of fitted MEM, TAMA, and DNLR-BP schemes (Figure 3(a)), (Figure 4(a)), and ( Figure 5(a)), it is reasonably assuming the reliability of Huxley's model in detecting 20 BioMed Research International the actual trend in present data. But in spite of this fact, the Q-Q lognormal and Q-Q normal diagrams, respectively, shown in Figure 3(c), Figure 4(c), and Figure 5(c), reveal that a normal distribution assumption for the random variable ϵ that determines the error terms to the MEM, TAMA, and DNLR-BP fits (Equations (2), (6), and (10) one to one) does not support a suitable model for the variability pattern inherent to the EHDS. To highlight the relevance of conveying a proper error structure, it is worth stating that assuming an additive error model in direct scales (Equation (10)) did not provide suitability of fit. Indeed, the ensuing DNLR-BP model resulted in a noticeably larger AIC value ðAIC = − 84528:9Þ than its MEM equivalent ðAIC = −94745:76Þ: Supposedly, the DNLR-BP should entail the same reproducibility features as a MEM counterpart. Nevertheless, when fitted to present data, an inherent multiplicative error structure of the MEM entailed a better performance, independently of the BP adaptation for heteroscedasticity that strengthened the additive DNLR scheme. Moreover, agreeing to Table 15, the DNLR-BP composites registered the most Table 13: Values of the reference parameters β, α, p , σ 1 , and σ 2 obtained by fitting the model of Equations (25) and (26) to the EHDS (EHDS row) and the corresponding ones; β s , α s , p s , σ 1s , and σ 1s deriving from simulated data (simulated row). We include RMSE values calculated through Equation (42)   The logistic distribution can reasonably model a higher dispersion than the normal one [36]. But, heavy tails displayed in the Q-Q loglogistic diagram of Figure 6(c) refer to an unfair MEM-loglogistic fit (Equation (9)). Equally, the Q-Q logistic plot of the residuals of a TAMA-logistic fit (Figure 6(f)) on geometrical space showing heavy tails also points to an inconsistent TAMA-logistic (Equation (13) (Table 15). Therefore, even switching to a logistic distribution, a perception of the unsuitability of an additive error structure maintains. Certainly, Figure 7(c) displays a heavy-tailed pattern, asymmetrical as in a DNLR-BP fit (cf. Figure 5(c)). Again, such an asymmetrical spread Figure 7(c) for the DNLR-BP-logistic scheme could explain by the attempt to model the actual heteroscedastic spread of Figure 1 through a BP-logistic adapted DNLR additive error model. And on top of that, Lai et al. [40] observation of biased direct nonlinear regression for large covariate values could also be pertinent. In any event, the spreads in the Q-Q loglogistic in Figure 6(c) and QQ-logistic plots in Figure 6(f) and Figure 7(c) bear that assuming a logistical distribution-driven error term also turns out to be unfeasible.
The MAD technique adapted by Echavarria-Heras et al. [28] delivered a quality-controlled data set missing about 25% of the total replicates in the crude EHDS. But, even though the MAD procedure removed a large share of the original data, the spread in Figure 8(a) and Figure 8(b) suggests the unreliability of a TAMA protocol when fitted to processed data. Still, regardless of an apparent curvature indicated by these spread plots, Figure 8(c) depicts a symmetric heavy tails arrangement. Therefore, that unravels that a TAMA fit failed to normalize the residual dispersion on processed data. Consequently, due to this fact and being aware of the fitting results of the protocols addressed so far, a conclusion seems indubitable. If the aim is to explore the possibility of maintaining Huxley's model-like systematic term in the analysis of the EHDS, we definitely ought to adapt the error driving random variable ϵ to conform to a symmetric distribution with zero mean and heavy tails. This perspective endorses consideration of a mixture of two normal distributions with zero mean and different deviations as a reasonable model for the randomness of ϵ: Table 15 allows comparison of the fits based on an M2N distribution to those related to the normal or logistic ones on the arithmetic scales, using the Akaike information criterion, which entails the M2N-shouldered models as better suited than the remaining ones. For instance, for composites involving a MEM, the reliability of fits remarkably improved when we switched from normal to logistic and then to a  (26) and (27)) already displays a fair spread. Similarly, Figure 9(f) showing the Q-Q M2N diagram corresponding to a TAMA-M2N fit (Equations (30) and (31)) indicates a spread unswerving to assumed M2N distribution. Tables 1-6 and 8-10 reveal that the relative standard errors of the estimates of the parameters α and β of the M2N model are smaller than those added to the normal or logistic counterparts. Moreover, reproducibility index values listed in Table 12 explain the advantage of the M2N distribution for a driving of good reproducibility strength of projections of averages of leaf biomass. This distribution enhanced the suitability of the DNLR-based methods even though the error structure implied by this protocol turned on unsuited. Indeed entries in Table 12 confirm that the DNLR-M2N-BP construct beats the TAMA-Poly(6) alternate in SEE, MPE, MPSE, and RD statistics. Additionally, comparison of CCC values in Table 12 unravels that the reproducibility power of a DNLR-BP-M2N rates equivalent to the one calculated for a TAMA-Poly(6) arrangement. Therefore, the M2N model outstrips its normal or logistic partners in the quality of the fit on the EHDS.
The statistics in Table 13 demonstrate that fitted reference values of the allometric parameters in the model of Equations (25) and (26) can be efficiently approximated through simulation procedures. Whatever bias among projections produced by mean response curves portrayed in Figure 12(b) and Figure 12(c) explains since the simulation procedure failed to reproduce a domed-like spread for the smaller leaf biomass values shown in the EHDS. In any event, differences in spreads among the EHDS and its simulated pool could provide further support to our judgment of a lack of systematization while gathering the former sample. Nevertheless, the closeness of reproducibility index values in Table 14 sustains a judgment that bias appointed by named differences in spreads bears irrelevant for practical purposes such as the allometric projection of leaf biomass values. Altogether, the result of addressed simulation runs allows Indeed, as we demonstrated here, it is possible to enhance the goodness-of-fit in allometric examination by holding up to Huxley's original theoretical perspective, that is, without recurring to the customary complex allometry and normal distribution approaches. It is also worth emphasizing that the presently offered analytical approach allows keeping the parsimonious structure of Huxley's model, thereby facilitating the exploration of a time invariance of the ensuing allometric relationship between eelgrass leaf biomass and area. This feature is paramount to genuinely nondestructive assessments of leaf biomass being relevant to eelgrass conservation.

Discussion
Huxley's power function model (cf. Equation (1)) beholds great relevance in biology. It offers a theoretical basis for static allometry that is conceived as the assessment of scaling relationships among individuals of a species, e.g., between one organ and total body size [55][56][57]). Huxley's model also provides a formal scheme for studies on evolutionary allometry that addresses scaling relationships between the sizes of organs of individuals across species [58][59][60]. And from a purely empirical perspective, whenever a time invariance of the parameters α and β upholds, Huxley's model can produce convenient nondestructive surrogates of the values of a given allometric response y by using previously fitted estimates and currently taken measurements of a descriptor variable x [61][62][63].
One relevant research subject circumscribing Huxley's model concerns the interpretation of the parameters α and β. For instance, Huxley stated that the intercept lnβ of the line deriving by log transformation on both sides of Equation (1) was of no particular importance. Still, the slope α was sufficient to mean static allometry on its own [2,64]. Furthermore, Huxley's perspective prevails in ongoing research as the valid theoretical standpoint for static allometry [59,60,65]. Another vital research subject around Huxley's model concerns the accuracy of estimates of the parameters α and β, which also relates to the suitability of the analytical scheme to get them. Supporters of the traditional perspective tie to a notion that an allometric response expressed through Equation (1) adapts to a multiplicative growth pattern in the direct scales of data (cf. Equation (2)) which also backs the TAMA (cf. Equation (6)) as the ordinary way to acquire estimates of α and β [40,43,60,[66][67][68][69][70][71][72][73]. But, other opinions sustain that a TAMA approach produces biased results, claiming that since the power function fostering Huxley's model is intrinsically nonlinear, the analysis should instead rely on a DNLR scheme (cf. Equation (8)) [74][75][76][77]. Moreover, there are settings where Huxley's model fails to produce a consistent fit. Exploring the reasons undermining the predictive power of this paradigm also endures a prime research subject. One attempt to address a lack of fit in Huxley's model adopts a standing that steers away from covariation among different traits, conceiving allometry as aligned on the covariation between size and shape [78,79]. From this standpoint, the analysis must rely upon regression schemes that include a systematic term expressed through modifications of Huxley's power function. So conceived variants engender constructs recognized as multiple parameter complex allometry forms (e.g., [75,[80][81][82][83][84]). Yet, since power series offer convenient representations of analytic functions, it turns reasonable to extend the complexity of Huxley's model to adapt polynomial regression schemes in geometrical space [37,71,[85][86][87][88]. But, embracing complex allometry attempts nurtures one of the most irreconcilable disagreements among schools of allometric examination. Indeed, for supporters of a TAMA way, examinations based on complex allometry constructs identified in the direct scales of data lose interpretation of biological theory to honor statistical appropriateness [59,60,73]. Besides defenders of traditional allometry claim that Huxley himself offered an approach to extending complexity while maintaining the original theoretical essence of static allometry. Indeed, when exploring the spread in the log-log plot of chela mass vs. body mass of fiddler crabs (Uca pugnax), Huxley acknowledged a breakpoint that was explained by an abrupt change in relative growth of the chela and assumed to take place around the time crabs reached sexual maturity [1,2,89]. The idea of a biphasic breakpointdetermined biological scaling conveyed the notion of nonlog-linear allometry [90][91][92][93][94][95]. Extension of Huxley's original idea of a biphasic scaling led to considering multiple breakpoints which in turn spawned the notion of polyphasic log- Table 15: Akaike information criterion values for comparison of the different protocols in direct scales contemplated in the present examination. AIC values appear in the intersection of rows for models and columns for the distribution of fundamental error random variable ϵ. For the TAMA, TAMA-logistic, TAMA-M2N, and TAMA-Poly(6) models fitted in geometrical space, we use the corresponding AICA values listed in Table 2, Table 5, Table 9, and Table 11 one to one. The ΔAIC entries to each AIC represent the difference in the AIC of each method relative to the smallest AIC value, that is, the one associated with the TAMA-M2N composite.  [96][97][98][99][100][101][102]. Broken-line regression techniques [37,[103][104][105][106][107][108] could deliver identification of breakpoints in polyphasic log-linear allometry schemes. Nevertheless, by relying upon nonlinear regression, broken-line schemes require starting values for the break-point estimation, so local maxima and inferences on estimates could make implementation difficult [107,109]. Then, even though polyphasic log-linear allometry adapts complexity by composing submodels, each one bearing Huxley's original theoretical envisioning, parameter estimation difficulties could make this approach problematic. Likewise, addressing complex multiparametric allometric forms could return overfitting, i.e., the statistical error of overestimation of the covariate's influence on determining the variability of the allometric response. So, in the tradeoffs to gaining reproducibility strength, overfitting could only offer an approach that overemphasizes empirical relevance, useful in reference only to the current data but not necessarily`to any other data sets. The present results illustrate an approach that enhances the reproducibility strength of Huxley's model while getting around overfitting associating with multiparameter complex allometry forms.
Maintaining Huxley's model within allometric examination offers relevant advantages for eelgrass conservation. Hamburg and Homann [110] used an allometric relationship to express leaf dry weight in eelgrass as a function of length and corresponding width. Solana-Arellano et al. [111] demonstrated that such a two-dimensional allometric dependence derives by assuming the validity of Huxley's laws of proportional growth between the leaf dry weight and the corresponding length and width one by one. Moreover, the leaf architecture in eelgrass approximates the related area as the product of the length times the associated width. Then, it is possible to call in Huxley's original setup and express eelgrass leaf dry weight as an allometric response of the linked area [112]. Likewise, the belt-shaped leaf expansion in eelgrass makes the corresponding length a good allometric descriptor of the dry weight [36]. Substantiation of these approaches allows the adaptation of allometric methods for the nondestructive estimation of aerial eelgrass biomass. These assessments are essential in remediation endeavors given current anthropogenic influences that threaten eelgrass permanence. But, given the general outline above, eelgrass allometric assessment methods based on Huxley's model are the subject of influences that affect their accuracy, precision, and pertinence. Remarkably, a quest for improving fit suitability of the sustaining Huxley's power function may suggest addressing multiparametric complex allometry forms pertinent, thereby drifting away from the theoretical perspective of Huxley's model of simple allometry [37]. It is worth emphasizing that focusing on Huxley's power function embraces parsimony that facilitates a concise exploration of an invariability property of the allometric parameters α and β, a condition upon which nondestructive assessments of eelgrass leaf biomass hinge in an essential manner.
Furthermore, the spread of present eelgrass leaf biomass to area data suggests strengthening the view by Montesinos-López et al. [36] of focusing on the suitability of error distri-bution to amend a lack of fit of Huxley's model. Indeed, it is worth recalling that the present data (EHDS) conforms to an extensive sample of (n = 10410) pairs of similar measurements to Montesinos-López et al. [36] data (MLDS) ðn = 537Þ, both collected at the Punta Banda estuary. They were, moreover, sampling to acquire the EHDS spanning over a complete year cycle, while MLDS stands for only one month. Therefore, as Figure 1(a) portrays, the Echavarria-Heras et al. [37] data set conforms to a remarkably more significant variability than that corresponding to the Montesinos-López et al. [36] sample (Figure 2(a)). Moreover, considerable variability in the EHDS could explain a lack of fit of a TAMA try reported by Echavarria-Heras et al. [37] and which we corroborated here. According to what these authors pointed out, our analysis of the EHDS also lets us become aware of apparent curvature in geometrical space. Likewise, we observed heterogeneity among dispersion patterns associated with the pools of smaller and larger leaves in the sample by analyzing the residual plot. Therefore, agreeing to a customary slant, in researching the EHDS, looking for complex allometry and a normal distribution approach seemed reasonable. But instead of attempting to amend TAMA's inconsistency based on a broken-line protocol or, more specifically, through a polynomial regression scheme as suggested in Echavarria-Heras et al. [37], we considered it adequate to keep Huxley's form of the systematic component while embedding the error term to nonnormal distribution, just as Montesinos-López et al. [36] devised when advocating for a logistic distribution-borne error term. But, for present data, residual dispersion accompanying a logistically distributed error term in geometrical scales could not improve a heavy tails pattern conforming to the normality of residues. We have to stress two points because the logistic error adaptation of TAMA by Montesinos-López et al. [36] failed to fit the present data well. First of all, as we have already pointed out, the EHDS is considerably larger in sampling dates, number of replicates, and number of processing participants than the MLDS. As a result, we stress that the present analysis deals with data expressing a significantly greater variability than that Montesinos-López et al. [36] dealt with. Second, heterogeneity of spread patterns among smaller and larger leaf sizes is more pronounced in the EHDS than it is in the MLDS. Perhaps this fact could explain the unsuitability of a unimodal distribution as a model for the random variable ϵ shaping the error term in present examination. Accordingly, at analyzing present data, a mixture of two normal distributions having a common zero mean but different standard deviations turned on suitable. The present review demonstrates that resultant allometric proxies of monthly averages of eelgrass leaf biomass exhibited similar reproducibility strength to those derived from the TAMA-Poly(6) scheme endorsed by Echavarría-Heras et al. [37].
Echavarría-Heras et al. [37] also advocated the advantages of a median absolute deviation-based data cleaning procedure to remove inconsistent replicates in the present data set. Indeed, we have stressed the significant heterogeneity in the variability of present leaf biomass replicates above. Perhaps such a spread ties to the participation of multiple 24 BioMed Research International data processing agents. Compared to Figure 1(a), the spread plots once data cleaning presented in Echavarría-Heras et al. [37] and Figure 8 here undoubtedly display a heteroscedastic pattern consistent with Huxley's model of simple allometry (cf. Equation (2)). Nevertheless, adjusting to a coherent spread to such a model was achieved only upon removing about 25% of the total number of replicates. Deletion of apparently undue replicates seems excessive, thereby questioning the pertinence of data cleaning at all as a suitable model improvement device. Then in analyzing the EHDS, data processing deems dubious even though Huxley's power function stands for the sound biological paradigm, masked by an overwhelming variability of replicates. It is worth adding that the markedly different spreads among smaller and larger leaf sizes in Figure 2(b) for the MLDS could have made it reasonable to call in data cleaning procedures. Nevertheless, as Montesinos-López et al. [36] demonstrated, a fit of a TAMA-logistic scheme in crude data turned on reliable. Furthermore, our results show that a revision of the complexity of error structure was enough to unmask the actual Huxley's power function-driven systematic trend despite a noticeable variability in crude data. Indeed, adapting the complexity of the error term distribution proposed here makes it unnecessary to call in data processing approaches to amend the reproducibility strength of a failing TAMA try. Then present results strengthen a perception of the appropriateness of the procedure suggested by Montesinos-López et al. [36]. Thus, it is possible to amend a lack of fit of a TAMA attempt without recurring to a non-log-linear complex allometry setup.

Conclusions
We have demonstrated that adapting the complexity of error structure could get around relying on multiple-parameter complex allometry forms as a mechanism to amend a lack of fit of Huxley's model. Maintaining this paradigm at the core of regression schemes aiming to identify scaling relationships offers a way to avoid the overfitting that could associate with using complex allometry-based amendments. For instance, inconveniencies could impair the corroboration of the time invariance of the parameters of the allometric relationship that substantiates nondestructive eelgrass aerial biomass assessments. Even though a lack of normality of residuals in the allometric examination is profusely acknowledged (e.g., [92,[113][114][115][116]), it is perhaps a drive to take advantage of the well-established assortment of statistical methods based on the normal distribution that mainly explains refraining from developing analytical tools conforming to alternative residual dispersion models. Given the present findings, it is pertinent to emphasize the advantages of the approach suggested by Montesinos-López et al. [36]. Certainly, allometric proxies for eelgrass leaf biomass produced by Huxley's model and normal mixture distribution composite delivered a reproducibility strength derived from complex allometry and normal distribution alternates [37]. But, intending to achieve empirical pertinence while keeping the parsimony entailed by Huxley's model, the present approach relies upon an error distribution assumption that requires elaboration from a theoretical perspective of allometry itself. Indeed, the problem of selecting the suitable regression model in allometry entails a statement about error structure. This problem remains unaddressed because, in allometry, error structure essentially depends on the inherent biological model and could not be resolved from statistical criteria alone . Perhaps the differences in spread patterns among smaller and larger leaf sizes that display in Figure 1(b) and Figure 2(b) could be the entrance to the path to providing a meaning to the involvement of a normal mixture distribution. To get substance to such an adaptation, we should rely on Huxley's quotation of breakpoint allometry itself. Indeed, the spread of present data in geometrical scales and corresponding to the MLDS suggests such a pattern. Indeed, it is reasonable to extend Huxley's original breakpoint idea so that we can justify the envisioning of a similar complexity of an error structure. Concisely for present data, we could assume that each of the pools of leaf sizes splitting by an identified breakpoint fits a parametrized normal distribution differently. The spreads of the collections before and after the breakpoint identified in the present data [37] corroborate suggested heterogeneity. Nevertheless, such a justification could spell a too simplistic rationale for embedding the mixture of two normal distributions in the allometric realm. But given the advantages that the present approach entails in identifying and validating traditional static allometry schemes for eelgrass conservation, it is worth exploring further substantiation and also the possibility of expanding its applicability elsewhere. Therefore, addressing further research on the subject is encouraged, and we particularly plan to undertake these tasks in upcoming contributions. BioMed Research International The AIC index permits evaluating the operation of different entrant models that fit a given data set. The model with the smallest AIC value is judged the best among competitors. The AIC index proves a settlement between the goodness-of-fit of a model and its complexity, which express through linked log-likelihood and number of parameters as a way to penalize inclusion of unnecessary ones. Since it bases on information entropy, an AIC index is frequently understood as an estimate of the information lost when a model is used to represent the process that generates the data. Lin's concordance correlation coefficient (CCC) evaluates how well one variable (Y) reproduces another (X); in other words, it represents a measure of the similarity (or agreement) between the two variables. The standard error of estimation (SEE) bears a global evaluation of goodnessof-fit of a model to observed data, since it measures the accuracy of (ŷ i ) predictions produced by a fitted regression model. This index takes nonnegative values. When SEE attains its minimum value, of zero, the observed values of the response coincide with the fitted mean response function, meaning that the model displays exact reproducibility of observed values. The MPE, which is now used to determine the goodness-of-fit of a model, is a standardized version of the coefficient of variation CV = ðSEE/ yÞ × 100 expressed as a percentage, as proposed by Schlaegen [120]. The MPSE backs a measure of the average absolute relative error, expressed as a percentage. This model assessment index was recommended by Schlaegen [120] and Meyer [121] as a measure of the absolute deviation of the expected and predicted responses, relative to the size of the prediction (jy i −ŷ i jŷ i ) expressed as a percentage average. For additional information on the use of the aforesaid statistics in allometric examination, the reader is referred to [122].

B. Comparison of Models Based on the AIC Index
In this appendix, we provide an explanation to the AICA entries in Tables 2, 5, 7, 9, and 11. To do so, we abide by the notation convention ðu, vÞ = ðln x, ln yÞ to understand the transformation that carries a point ðx, yÞ in arithmetical space into one ðu, vÞ in logscales.
In allometry, we ought to compare different models based on the AIC index, where some models are fitted in the arithmetic space and others in log scales. Whenever a model was fitted in log scales, we denote its log-likelihood through the symbol log G ð b θÞ. If we want to compare such a model with another one fitted in the direct scales of data, we are required to transport the likelihoods and loglikelihoods from log scales to arithmetical ones, because allometric data intrinsically associate there. Particularly, we denote the statistics resulting from transferring log G ð b θÞ to the arithmetical space by means of the symbol log A ð b θÞ. The conversion of concurring log-likelihood log G ð b θÞ into its arithmetical space equivalence log A ð b θÞ achieves through the relationship: where b θ stands for the set of estimated parameters for a given model. Equivalently, once a model fitted in geometrical space attains an AIC value, we make it correspond to its equivalent statistics in direct scales which we denote through AICA. Necessarily, the AICA derives from the log-likelihood expressed in direct scales log A ð b θÞ, that is, We explain how Equation (A.7) derives for the case in which the y as a random variable is lognormally distributed with log-mean μ and log-deviation σ, that is, y~lognormð μ, σÞ .Under such an assumption, v is a normally distributed having mean μ and deviation σ, that is, v~Nðμ, σÞ. We begin by establishing the relationship between the loglikelihoods in the arithmetic and geometric scales. Since we have assumed that y~lognormðμ, σÞ, the expression for the likelihood function in arithmetic scales log A ðθÞ becomes the following: Since v i = log ðy i Þ, i = 1, 2, ⋯, n, then we have the following: or equivalently, ðB:6Þ Therefore, To obtain the relation between the log-likelihoods log A ðθÞ and log G ðθÞ, we take logarithms on both sides of Equation (B.6) which yields the following:

Data Availability
Data will be available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.

Supplementary Materials
We provide MATLAB and R codes for the related identification tasks (Computer_codes.doc). Each code description includes comments referring to the paper section it relates to and what it generates (table of fitted parameters, figure, or reproducibility index). MATLAB 2016a or later version and RStudio are required. The set of observed data (referred by means of an EHDS acronym in the manuscript) could be available from the corresponding author-provided acceptance of a fair use agreement. However, we readily include the set of processed data resulting from the MAD procedure referred on Section 3.6 (dpro.text) and the simulated data in Section 3.11 (dsim_10410.text) that allow testing the performance of the codes. (Supplementary Materials)