Examination of the Effects of Curvature in Geometrical Space on Accuracy of Scaling Derived Projections of Plant Biomass Units: Applications to the Assessment of Average Leaf Biomass in Eelgrass Shoots

Conservation of eelgrass relies on transplants and evaluation of success depends on nondestructive measurements of average leaf biomass in shoots among other variables. Allometric proxies offer a convenient way to assessments. Identifying surrogates via log transformation and linear regression can set biased results. Views conceive this approach to be meaningful, asserting that curvature in geometrical space explains bias. Inappropriateness of correction factor of retransformation bias could also explain inconsistencies. Accounting for nonlinearity of the log transformed response relied on a generalized allometric model. Scaling parameters depend continuously on the descriptor. Joining correction factor is conceived as the partial sum of series expansion of mean retransformed residuals leading to highest reproducibility strength. Fits of particular characterizations of the generalized curvature model conveyed outstanding reproducibility of average eelgrass leaf biomass in shoots. Although nonlinear heteroscedastic regression resulted also to be suitable, only log transformation approaches can unmask a size related differentiation in growth form of the leaf. Generally, whenever structure of regression error is undetermined, choosing a suitable form of retransformation correction factor becomes elusive. Compared to customary nonparametric characterizations of this correction factor, present form proved more efficient. We expect that offered generalized allometric model along with proposed correction factor form provides a suitable analytical arrangement for the general settings of allometric examination.


Introduction
The model of relative growth of Huxley [1] is formally stated by means of a scaling relationship of the form where and are measurable traits and the parameter is designated as the allometric exponent, while is identified as the normalization constant. This model, also termed equation of simple allometry, has been extensively used in research problems in biology [1][2][3][4][5], physics [6], economics [7], earth sciences [8], resource management, and conservation [9,10], among other fields. Eelgrass provides nursery for waterfowl and fish species. By trapping sediment and stumping wave energy, this seagrass promotes shoreline stabilization. Eelgrass services also include nutrient recycling, water filtration, and carbon dioxide removal. Current anthropogenic influences threaten eelgrass permanence. Conservation efforts rely on plot transplanting in a fundamental way. Monitoring effectiveness 2 BioMed Research International depends on measurements of standing stock and productivity through time. This makes the assessment of average leaf biomass in shoots a necessary input. But traditional estimation of eelgrass leaf biomass units relies on destructive methods. This could alter shoot density in a developing transplant. Thus evaluation renders indirect assessment methods necessary [10,11]. Results show that an allometric scaling of the form of (1) for eelgrass leaf biomass and associated area is consistent [11]. Derived projections of individual leaf biomass convey useful surrogates for mean leaf biomass in shoots. Moreover, estimates of the parameters and are invariant within a given geographical region [10][11][12]. Hence, estimates fitted at site can endow suitable projections of leaf biomass values currently observed in other places of the region. This bears the referred allometric projections of a convenient nondestructive feature ( [11,13]).
Simplicity of (1) makes allometric projection of eelgrass leaf biomass units convenient. But there are caveats on dependability. For instance, even though and are invariant, environmental influences can induce a relative extent of variability on local estimates ( [12,14]). Besides, the response in the power function-like scaling of (1) is very sensitive to variation of parameter estimates. Then, accuracy of derived proxies is subject of error propagation of estimates. In addition, there are factors of biological scaling which can influence precision of estimates (e.g., [10,11,13,15]). Packard [16] questioned results of Mascaro et al. [17] on allometric examination, and Mascaro et al. [18] responded to criticism. Going over this exchange highlights the relevance of procedural factors in determining precision of parameter estimates of allometric scaling. It also offers a convenient framework for the aims of the present research.
An important factor influencing precision of estimates of allometric parameters is analysis method. A widespread approach is the traditional analysis method of allometry (TAMA hereafter). It relies on a log transformation of data in arithmetical scale in order to contemplate a linear regression model in geometrical scale. Then, the fitted line is backtransformed to yield a two-parameter power function in the original scale. But embracing this procedure fuels a vivid unresolved debate. Views assert that this protocol can lead to biased results (e.g., [16,[19][20][21][22][23][24][25][26][27][28][29]). And other practitioners consistently wed to the idea that logarithmic transformations are necessary (e.g., [18,[30][31][32][33][34][35][36][37][38][39]). An alternative to the TAMA approach is using nonlinear regression methods in the direct scale of the data [16]. Echavarria-Heras et al. [11] concluded that producing allometric projections of average leaf biomass in eelgrass shoots must rely on this protocol. Yet a direct nonlinear regression approach in allometry is also not unfailing. For instance, inadequate identification of the inherent error structure can lead to significant bias [18]. Besides, Lai et al. [30] found that estimates of allometric parameters fitted by nonlinear regression can exhibit a high sensitivity to the largest values of the covariate. Therefore, evaluation of analysis method suitability in acquiring consistent eelgrass leaf biomass proxies needs revision.
The adoption of methods of curvature in geometrical space could offer a way to overcome inadequacy of the TAMA procedure ( [27,[40][41][42]). In particular, it is pertinent to examine if taking curvature into account leads to improved accuracy of eelgrass leaf biomass proxies. But, according to Mascaro et al [18], curvature could manifest because of methodological factors of data gathering. Thus, an examination of the effect of curvature in eelgrass leaf biomass allometry must also take into account a possible participation of data quality effects. Mascaro et al. [18] reminds on three ways of handling curvilinearity in geometrical space. One is by separating data to contemplate different local linear models to account for heterogeneity of effects of the covariate [43][44][45]. A second one is by fitting a polynomial model [46][47][48][49]. A third approach endorses direct nonlinear regression assuming a heteroscedastic error structure as contemplated by Mascaro et al. [17]. Either approach above bears complexity beyond the linear model in geometrical space that associates with the customary bivariate power function of allometry. This suggested putting forward a generalized allometric model intended to deal with curvature in geometrical space. This paradigm incorporates parameters that change as continuous functions of the log transformed covariate ( [27,39,46]). As we explain further on, the curvature arrangements recommended by Mascaro et al. [18] can be all derived from the offered formalization. Moreover, a nonzero intercept power function that Packard [50] recommends to handle curvilinearity in geometrical space also derives from the presented generalized scaling model.
But any scheme addressing curvature in geometrical space depends on a factor for correction of bias of retransformation of the regression error. In the general settings, if stands for the regression error, then the said correction factor, through denoted using the symbol ( ), is given by the mean of the exponentiated error random variable; that is, ( ) = ( ∈ ) [51][52][53]. Furthermore, the TAMA approach relies on the essential assumption that is additive, normally distributed, and homoscedastic [33]. When this happens, ( ) takes its lognormal-mean form [51][52][53]. But if fails to be normally distributed, there are two possibilities. If the distribution of is known, we could derive a closed form for ( ). In turn, if the error distribution is not identified a priori, a widespread approach is taking ( ) in the nonparametric form given by the smearing estimator of bias of Duan [54]. Still, there are provisions on this. A smearing estimate form can fail to compensate the downward retransformation bias of logged data ( [53,55,56]). Thus, in a circumstance where is unspecified, characterizing ( ) seems elusive. Here, we put forward an arrangement for ( ) aimed to get around this circumstance. Zeng and Tang [57] proposed a nonparametric alternate to the smearing form. It matches the first three terms' partial sum of a power series expression of ( ∈ ), assuming ( ) = 0. Form suggested here corresponds to a generalization of this construct. It does not abide the restriction ( ) = 0 and matches an −terms partial sum approximation of the exponential series representation of ( ∈ ). The partial sum that maximizes reproducibility strength of retransformed mean response sets criterion to choose .
Present results show that a consideration of curvature in geometrical space, as well as a suitable characterization of the correction factor of retransformation bias, offers consistent allometric proxies of observed mean leaf biomass in eelgrass shoots. Hence, contrary to views asserting direct nonlinear regression as mandatory in allometric examination, our findings validate a parallel reliability of log transformation based methods. This is well in line with claims of Mascaro et al. [18] and many others about blaming the use of logarithms of incongruent results in allometric analysis. Moreover, keeping the analysis in geometrical space unraveled heterogeneity in the inherent leaf biomass scaling pattern. This could not be achieved by clinging to direct nonlinear regression in arithmetic space as the only valid approach of allometric examination. Offered analytical arrangement is expected to be applicable in the general settings of allometry.

Data.
For the aims of the present research, we relied on an extensive eelgrass data set collected in San Quintin Bay, a coastal lagoon on the Pacific side of the Baja California Peninsula, México (30 ∘ 30' N -116 ∘ 00 W), and through a 13 months' long sampling period covering a whole-year cycle. Data composes measurements of length (mm), width (mm), and dry weight (g) of a total of 10412 individual eelgrass leaves taken from 20 randomly thrown 400cm 2 quadrats every monthly visit to the site. A sampling visit will be further referred to as "sampling time" in the text. The length times width proxy [11] provided estimations of leaf area (mm 2 ). In order to test for methodological influences of data gathering, we processed raw data set according to a mean plus or minus two standard deviations outlier's removal procedure [58]. Appendix A presents results of an exploratory analysis of data.

Models.
As above specified symbols, and stand for the biomass of an individual eelgrass leaf and its respective area one to one. Echavarría-Heras et al. [11] assert that these variables can be related through the bivariate allometric model of (1). One procedure to acquire estimates for the parameters and is fitting directly in arithmetical scale a nonlinear homoscedastic regression model. Besides, we can use a TAMA approach, that is, fitting the linear regression model where V = ln , = ln , and is a random error term assumed to be normally distributed with zero mean and variance 2 , that is, ∼ (0, 2 ).
We conceive curvature in geometrical space as a circumstance where fitting results of regression model of (2) are inconsistent. Dealing with this situation amounts for considering complexity beyond incorporated by (1). One possible approach to address curvature is assuming that scaling parameters and in (2) depend continuously on the covariate ( [27,39,46]). This is consistent with the generalized allometric model, with ( ) and ( ) intended to be continuous and differentiable functions defined on + and with ( ) being positive.
Certainly a log transformation V = ln , = ln of (3) establishes the regression model where V ( , ) = ln ( ( )) + ( ( )) , where , a residual error term, is conceived as a random variable that in the general settings is -distributed with mean and variance set by a function 2 ( ) of the covariate , that is, ∼ ( , 2 ( )). Setting ( ( )) = and ( ( )) = with and constants reduces (4) to the regression model of (2). In Appendix B, we explain that (3) accommodates all curvature paradigms suggested by Mascaro et al. [18]. These include a biphasic and a polynomial model in geometrical space, as well as the nonlinear heteroscedastic model referred to by Mascaro et al. [17] in direct arithmetical space. Moreover, as shown in Appendix B, the three-parameter power function chosen as an alternate standard for curvature [59] also derives from (3).

Biphasic Model in Geometrical Space.
In order to characterize the model of (3) in a biphasic mode, we let ( ) = ( ) and ( ) = ( ), such that including parameters and and the function ( ( )) given by for = 1, 2. ( − ) is a Heaviside function ( ) [60], evaluated at = − and correspondingly = ln , with ≤ ≤ being a point separating growth phases ≤ and > . Then, denoting by means of V ( , ) the resulting form of V ( , ) from (4), we get the biphasic regression model where is a random error term as defined in (4) and with ( ( )) = ln + ( ) , where and for = 1, 2 parameters are to be estimated from data.
and for 1 ≤ ≤ are coefficients; one can acquire a polynomial representation V ( , ) for the generalized mean response function in geometrical space V ( , ). This way the polynomial form of regression (4) becomes where is a random error term as defined in (4), with and , for = 0, 2, . . . , parameters.

Nonlinear Heteroscedastic Model in Arithmetical Space.
As we explain ahead, direct algebraic manipulation of (3) leads to the consideration of the nonlinear heteroscedastic regression model addressed by Mascaro et al. [17]; namely, with , , and being parameters and being a zero mean normally distributed error term with covariate dependent variance 2 ( ) = 2 2 , that is, ∼ (0, 2 2 ).
A nonlinear homoscedastic form derives from (14) by setting = 0; that is, And, again is an additive error term assumed to be normally distributed with zero mean and variance 2 , that is, Appendix A deals with exploratory analysis of data Appendix B presents notation convention and also explains how all addressed paradigms derive from the generalized model of (3). Appendix C explains the addressed forms of correction factor for bias of retransformation of the regression error. Fitting results of the geometrical space based models appear in Appendix D. Those corresponding to the nonlinear heteroscedastic and homoscedastic models pertain to Appendix E. Agreement between observed and projected values is commonly evaluated by analyzing values of Lin's Concordance Correlation Coefficient (CCC) [61]. This correspondence index is commonly denoted by means of the symbol . Agreement will be defined as poor whenever < 0.90, moderate for 0.90 ≤ < 0.95, good for 0.95 ≤ < 0.99, or excellent for ≥ 0.99 [62]. Besides CCC values, we assessed reproducibility by comparing goodnessof-fit statistics, such as the coefficient of determination (CD), standard error of estimate (SEE), mean prediction error (MPE), total relative error (TRE), average systematic error (ASE), and mean percent standard error (MPSE) ( [63][64][65]). For statistical tasks, we relied on the R package release 3.5.

Results
Exploratory analysis in Appendix A identifies maximum, minimum, and sample mean values for observed leaf area values and associated dry weights . We also explain distribution of variables in terms of quantiles of probability 0.1, 0.25, 0.50, 0.75, and 0.90, for both crude and processed data. Statistical exploration extends to log transformed values of these variables. We present Q-Q plots (quantile-quantile) for comparison of distribution patterns, as well as boxplots for the 13 months' long sampling scheme for both crude and processed data. We can learn that, from month 2 to month 6, a reduction in the values of weight and area occurred; this perhaps is explained by an increase in temperature during the period. We can be also aware that a similar variation pattern over time is shown in both raw and processed data sets.

Fitting Results of Geometrical Space Models.
In order to validate curvature in geometrical space, we compared the linear model derived from (1) as well as biphasic and polynomial alternates derived from the generalized model of (3). Appendix B explains formal matters. Tables 1 and 2 summarize notation convention. Equations numbered beyond (15) belong to the appendices. Appendix D explains corresponding regression protocols.
Fitting results of the TAMA arrangement of (2) appear in Appendix D. Figure 1 shows the spread about TAMA's linear mean response function (V | ). We can visually ascertain that deviations from the linear mean response function (V | ) suggest curvature (red dots). Thus, data processing removed inconsistent replicates but shown spread still deviates from a linear mean response. This suggests that curvature in geometrical space could not be explained by methodological factors related to data gathering.
Fitting results of the biphasic protocol of (8) are summarized in Appendix D. Figure 2 displays the spread about mean response function (V | ) in geometrical space. Compared with Figure 1, we can ascertain that the biphasic fit provides a consistent account of different variation patterns among smaller and larger leaves. We can visually ascertain that fit produced consistent results. This confirms a judgement that identified curvature might be due to intrinsic factors of leaf growth rather than methodological influences related to data gathering.
Appendix D presents fitting results of the polynomial model ( = 6). Figure 3 displays dispersion about the polynomial mean response function in geometrical space (V | ). A polynomial representation also exhibits higher consistency than the TAMA arrangement. Recalling the biphasic scheme, the polynomial suggests a smooth transition between two growing phases.

Model Selection in Geometrical
Space. Assessment of models fitted on geometrical space relied on goodness-of-fit statistics, that is, the coefficient of determination, standard error of estimate, mean prediction error, total relative error, average systematic error, and mean percent standard error ( [63][64][65]). Besides, we took into account concordance correlation coefficient [61] and Akaike's information index [66]. Table 1: Summary of notation convention of the bivariate scaling model of (1) and its generalization to account for curvature as given by (3). GS stands for geometrical space, AS stands for arithmetical space, and CF means correction factor of retransformation bias.

Typical bivariate
Eq. Generalized for curvature Eq. Table 2: Notation convention for curvature models in geometrical space. We include the biphasic and polynomial characterizations of the generalized allometric model of (3).GS stands for geometrical space, AS stands for arithmetical space, and CF means correction factor of retransformation bias of the regression error.

CF Zeng and Tang
Biphasic model Eq.
Polynomial model Eq. Table 3 presents results. Goodness-of-fit statistics and and AIC values disfavored the TAMA protocol. On the contrary, comparison indices favored the biphasic model. Moreover, differences among indices but TRE and ASE for this scheme and the polynomial ( = 6) are slight. Particularly, the highest AIC is associated with the TAMA protocol (AIC = 15069.9, ûAIC = 2491.9). Therefore, this model bears the less support. The biphasic choice delivered the smallest AIC's value (AIC = 12578.0, ûAIC = 0). Nevertheless, difference in AIC is just barely relative to the = 6 polynomial model, since this choice conveyed (AIC = 12615.0 and ûAIC = 37). Thus, model confrontation shows that the TAMA protocol is unsuited, thus backing the assertion that whatever model aims to be consistent with the present data, it ought to be nonlinear in geometrical space.

Retransformation
Results. The TAMA protocol was not supported by the model selection criteria. Anyway, for comparison, corresponding retransformation results are included in Appendix D. Related with the TAMA protocol, fitting results of the biphasic model display a relatively improved distribution of residuals about the zero line. Nevertheless, normal Q-Q plot still shows heavier tails than those expected for a normal distribution. And, again both test statistics and Table 3: Assessment of geometrical space fitted models. Comparison took into account goodness-of-fit statistics, that is, the coefficient of determination (R 2 ) standard error of estimate (SEE), mean prediction error (MPE), total relative error (TRE), average systematic error (ASE), and mean percent standard error (MPSE) ( [63][64][65]  p values of an Anderson-Darling test [67] provide evidence against normality of residuals. This justifies choosing the nonparametric forms ( ), ( ), or ( ) for compensation of downward bias induced by retransformation of the regression error ( [11,52,53]). Table 4 displays comparison statistics for the reproducibility strength of the biphasic mean response ( | ) as shaped by the different forms of ( ). We can learn that agreement between the biphasic mean response and leaf biomass data is best for ( ). Figure 4 shows spread of processed leaf biomass values about the biphasic mean response function ( | ) as shaped by the considered forms of ( ). We can observe that both ( ) and ( ) overcompensate the bias correction by ( ). Moreover results show that as opposed to the TAMA a biphasic protocol along with the ( ) form offers consistent proxies of individual leaf biomass. But it is worth mentioning that in spite of the fact that model selection favored the biphasic scheme, examination of the polynomial model output reveals similar predictive strength to the biphasic alternate.

Assessing Curvature by Direct Nonlinear Regression.
As suggested by Mascaro et al. [18], effects of curvature in geometrical space can be analyzed by means of the direct nonlinear heteroscedastic regression model of (14). In Appendix B, we explain that such a protocol also derives from the generalized bivariate allometric model of (3). Table 5 presents pertinent notation convention. For comparison, we also present results for the associated homoscedastic case.
Fitting results of the heteroscedastic and homoscedastic models appear in Appendix E. We can learn that estimates for the normalization constant and scaling exponent parameters are very similar. Certainly, corresponding 95% confidence intervals display some overlap. As a result, we can expect similar reproducibility features for both models. Table 6 presents comparison statistics.
We can be aware that model assessment backs the heteroscedastic model. But selection here is mainly on qualitative grounds. It actually concerns the ability of the heteroscedastic model to identify an expected dependence of variance in the covariate. Certainly, the reproducibility strengths of both paradigms are equivalent. Indeed, Figure 5 shows that mean response curves ( | ) and ( | ) differ just barely.
Results show that as it occurred for models fitted in geometrical space, data cleaning failed to correct a heavy tails problem for the nonlinear fits. This can be ascertained from the normal Q-Q plot of residuals. This strengthens our point on the consideration of a different error structure from the one assumed here. Exploring the effects of error structure in the fitting of models for curvature addressed here will be a matter of further research. Interestingly, both the homoscedastic and heteroscedastic models seem to induce the same reproducibility strengths. Table 3 favored the biphasic protocol. Correspondingly, statistics in Table 6 support the nonlinear heteroscedastic model. Results of Table 4 endure ( ) as required for largest reproducibility of retransformation output. Table 7 allows assessment of these models. We can learn that half the number of comparison indices coincide ( , R 2 , SEE, and MPE). In addition, the biphasic model is favored by AIC, ASE, and MPSE. This sets criterion for selection of curvature in geometric space as a consistent paradigm for the present data. Accordingly, the biphasic model bears adequate.

Implications for Allometric Proxies of Mean Leaf Biomass in Eelgrass Shoots.
We in turn consider allometric proxies for average leaf biomass in eelgrass shoots. In getting these surrogates, we aggregate allometric projections of individual leaf biomass conforming a shoot. For comparison, we consider individual leaf biomass surrogates produced by the   Figure 1, we can be aware that, as opposite to the TAMA protocol, the polynomial scheme offers a consistent account of curvature.   (14) and (15), respectively.
Nonlinear model Heteroscedastic Eq. Homoscedastic Eq.   [11] stablished that proxies derived from the TAMA protocol are inconsistent with observed values. This endorsed nonlinear regression in the direct scale as a requirement for reliability of allometric projections of mean leaf biomass in eelgrass shoots. But Table 8 shows that a curvature model fitted in geometrical space can offer proxies entailing similar predictive power to a nonlinear regression protocol. Plots in Figure 6 allow getting a glimpse of this assertion.

Discussion
The customary bivariate allometric model of (1) offers nondestructive surrogates for average leaf biomass in eelgrass shoots [11]. But there are methodological factors that could influence dependability. Views assert that parameter identification based on logarithmic transformations leads to biased projections [20][21][22][23][24][25][26][27][28][29]. But other practitioners clung to this approach as meaningful and necessary in allometric examination [30][31][32][33][34][35][36][37][38][39]. This going over suggests that surpassing this controversy amounts to considering curvature in geometrical space. For that aim, we proposed the generalized model of (3). Approaches such as direct nonlinear heteroscedastic regression, as well as biphasic and polynomial protocols in geometrical space [18], became logical resultants from this construct. For present data model selection validated maintenance of the analysis in geometrical space. Nevertheless, at an empirical level, addressed protocols produced allometric projections of individual leaf biomass of correspondent precision. This was also verified for concomitant projections of average leaf biomass in shoots. But, from a qualitative standpoint, the nonlinear regression protocol mainly contributed by identifying expected dependence of the variance on covariate. Moreover. Figure 7(a) depicts manifest differences in mean response trends between the polynomial fit and the nonlinear heteroscedastic model. Nonetheless, those in Figure 7 the biphasic models differ but only barely. Then, a nonlinear regression scheme at best shaped a reasonable approximation of the mean response function resultant from curvature methods. Differences in patterns of the biphasic and polynomial mean response functions relative to the nonlinear protocol exhibit that clinging to this last paradigm could impair detection of the true allometric relationship. Moreover, relying on direct nonlinear regression impairs identification of  heterogeneity in the log transformed response as covariate changes. This further stresses on limitations of this device as a tool for allometric examination ( [18,30]). Oppositely, output of the selected biphasic model shown in Figure 2 suggests differentiation of growth patterns among smaller and larger leaves. Besides, the polynomial mean response in Figure 3 suggests a gradual transition between different growth phases. Thus, as opposed to direct nonlinear regression, a consideration of curvature in geometrical space could elucidate an inherent leaf growth pattern. This strengthens a judgement that the log transformation step, essential to traditional allometric examination, cannot be thrown away without losing relevant information ( [33,37]).
Mascaro et al. [18] conceived curvature in geometrical space as related to methodological factors of data gathering. But present examination corroborated consistency of curvature models for processed data. This suggests that manifestation of curvature is rather explained by intrinsic factors in leaf growth. Additionally, data processing failed to amend the heavy tails problem detected on Q-Q plots. This indicates departure of residuals from an assumed error structure. As a result, numerical values of the addressed correction factor forms turned out to be different, thus conveying ambiguity in selection of suitable mean response of models fitted in geometrical space. This could entail the only advantage of nonlinear regression method over log-transformationcurvature paradigms. But, also for this analysis method, an inadequate postulation of inherent error structure can lead to significant bias [18]. It seems then reasonable considering that a suitable characterization of error structure could lead to robustness of built allometric proxies, even when they derive from crude data. Steering to an error structure different from what is assumed here is a worthwhile subject of further research.
When dealing with similar data we suggest taking into account recommendations that come up from this examination. First, it is highly advisable to perform a preliminary examination of the spread around the straight line in geometrical space resulting from the model of (1). If further statistical exploration confirms that linearity and assumed error structure are consistent with data, Huxley's bivariate allometric model could suit. Otherwise, the arrangement of curvature, error structure, and correction factor form such as proposed here could be called into account for the analysis. The use of data cleaning procedures in order to achieve a better fit is controversial [11]. Instead of performing data processing a posteriori, it is highly advisable to rely on standardized data gathering procedures. This will prevent proliferation of inconsistent replicates that could exacerbate a heavy tails problem on Q-Q plots.

Conclusion
Failure to perform both a preliminary exploration of spread of log transformed allometric data and a sound evaluation of model adequacy could impair detecting a possible manifestation of curvature. As a consequence, the output of a traditional analysis method could set biased predictions of observed values. This circumstance could result in dismissal of a log transformation step in the analysis, giving way to contemplation of direct nonlinear regression as the only protocol to acquire reliable parameter estimates [11]. Results of this examination suggest that consideration of curvature in geometrical space as set by the model of (3) could offer dependable allometric proxies of average leaf biomass in eelgrass shoots.
From a general perspective, complexity as encompassed by the model of (3) can stand for curvature as conceived in allometric examination. Particularly, biphasic or polynomial protocols in geometrical space, as well as a direct nonlinear heteroscedastic regression model, derive as particular characterizations of this paradigm. Moreover all statistical models for accurate estimates of relative growth contemplated by Bervian et, al., [68] can be also accomodated by the present generalization of the model of simple allometry of Huxley. But empirical convenience on its own does not validate adoption of this paradigm as a general tool. Certainly, the Weierstrass approximation theorem [69] backs a polynomial regression model as a reasonable identification device for the generalized allometric model expressed in geometrical space. But suitability of retransformation results will sensibly depend on correction factor form. And a mean function resulting from a polynomial fitted in geometrical space will not enable characterization of functions ( ) and ( ) one to one. Furthermore, complexity of (3) could pose significant difficulties while attempting its identification through direct nonlinear regression methods. Needless to say, biological interpretation of the scaling functions ( ) and ( ) is also pending. A quest for efficient tools of nondestructive assessment of plant biomass units justifies addressing these examinations in a further research.

A. Data Exploratory Analysis
This appendix presents an exploratory analysis of raw and processed data sets contemplated in this examination. Table 9 describes the distribution pattern, in terms of quantiles, for a sample of 10412 measurements of eelgrass leaf weights and related areas taken over 13 months conforming  Figure 7: Trends of mean response functions calculated from the retransformed outputs of polynomial (a) and biphasic models (b) involving the ( ) form of ( ). We can observe a differentiation in trends relative to that of a power function fitted in directly by means of the nonlinear heteroscedastic protocol of (14). Although relative deviations manifest for the biphasic model, differences in mean response patterns are more clearly depicted for the polynomial choice.
the present raw data. The first four columns in the uppermost row label the minimum followed by quantiles of probability 0.1, 0.25, and 0.50. Correspondingly, the fifth column in the first row presents the sample mean followed by quantiles of probability 0.75 and 0.90 before the maximum. Second row presents leaf dry weight values ( ) and third row shows corresponding areas . Third and fourth rows present transformed ln and ln values one to one. Similarly, Table 10 shows the variation pattern of the 10023 observations resulting after applying to raw data a mean plus or minus two standard deviations outlier's removal procedure [58] designed to remove replicates considered significantly discrepant from the mean response function of a simple allometric model for a leaf dry weight response in terms of a leaf area covariate.
Comparing quantile values for the leaf dry weight ( ) and corresponding area ( ) variables reported in Tables 9 and 10 for crude and processed data respectively, we can ascertain that in spite of removing a large number of discrepant observations by data cleaning, both original and remnant sets show an equivalent distribution pattern. This similarity in distribution before and after data processing of the data can be better perceived in the Q-Q (quantile-quantile) graphs observed in Figure 8 for raw (a) and processed data (b), respectively.
Q-Q plots in Figure 8 compare the distributions (regardless what they are) of the leaf area and associated dry weight variables. The linear relationship observed between the quantiles of both variables underlines that both variables have similar distributions, although with different parameters. This similarity between distributions of both variables is observed for both sets of observations. The great similarity of the fitted straight lines for both sets of observations in Table 11 confirms that the referred distribution pattern does not change after data processing. Figures 9 and 10 display boxplots for the 13 months' long sampling scheme. We can learn that, for raw data, from month 2 to month 6, a reduction in the values of both leaf dry weight and linked area occurred. This is perhaps explained by an increase in temperature during those months. Processed data exhibits similar dynamics through time for these variables. Moreover, the overall variation patterns of raw and processed data, throughout the 13 months of sampling, are similar

B. Derivation of Regression Protocols
In this appendix, we first present notation convention for statistics related to the generalized bivariate model of (3). We then explain how the different regression protocols addressed in this examination derive from this paradigm. We also include notation convention for related statics.

B.1. Generalized Bivariate Allometric Model.
A transformation V = ln and ( ) = ln of (3) establishes the generalized regression model in geometrical space given by (4). The term V ( , ) stands for the corresponding mean response function. Using the customary notation convention, we represent V ( , ) through symbol (V | ); that is,    with ( ) = ( ) interpreted as a correction factor of bias of retransformation of the regression error .

B.2. Traditional Analysis Method of Allometry (TAMA)
. The bivariate allometric model of (1) derives from the model of (3) setting the scaling parameters constant; that is, ( ) = and ( ) = . The resultant regression model in geometrical space is given by (2). Corresponding mean response function will be denoted here by means of the symbol (V | ). Hence, from (2)  This explains (8). Since as given by (3) is assumed to be a continuous function of , this property ought to be satisfied by the biphasic model; that is, we require to set a condition, 1 ( ) = 2 ( ). This leads to the equation Moreover, V ( , ) as given by (8) can be equivalently represented by with the continuity condition of (B.15).
Correspondingly, the regression model of (4) in its biphasic form becomes As we explained in (4), is a residual error term conceived as a random variable distributed according to an unknown distribution having mean and variance set by a function 2 ( ) of the covariate ; that is, ∼ ( , 2 ( )).
The form of the generalized mean response (V | ) for the biphasic model is denoted through (V | ) and becomes It follows from (B.16) that the expression for (V | ) can be equivalently represented by

B.4. Polynomial Model in Geometrical Space.
In this section, we explain how (12) can be derived from the generalized structure of model of (3). We also set the notation convention for related statistics. For these aims, we begin by considering polynomials ( ) and ( ) given by where, for 1 ≤ ≤ , and stand for coefficients. Now (B.31) This explains the result of (12). Correspondingly, a polynomial characterization of the regression model of (4) takes the form with being a residual error as described in (4). The form of the generalized mean response (V | ) for the polynomial characterization is denoted through (V | ). It becomes

B.5. Nonlinear Regression Models.
In this section, we explain how the nonlinear heteroscedastic model of (14) can be also associated with (3). We also explain the notation convention for related statistics. For these aims, we begin by noticing that V ( , ) in (5)  which suggests the nonlinear heteroscedastic regression model of (14). That is, which can be considered as a tool to analyze curvature in geometrical space. Particularly, by setting = 0, we can consider the homoscedastic model where in (B.43) and (B.44) the error term is assumed to be normally distributed random variable having zero mean and homogeneous variance; that is, ∼ (0, 2 ). It turns out that the mean response function ( | ) associated with the heteroscedastic model becomes The mean response function in arithmetical space becomes Packard [59] asserts that, commonly, any data set on arithmetical scale that is consistently described by a threeparameter power function will track a curved path when transformed to the logarithmic scale. Equations (B.45) through (B.47) provide a formal set-up accommodating such a statement.

C. Forms of the Correction Factor of Bias of Retransformation of the Regression Error
In this appendix, we explain the different characterizations of ( ), defined as a correction factor for bias of retransformation of the regression error introduced by (B.3) for ( | ).
If is normally distributed with zero mean and constant variance 2 , that is ∼ (0, 2 ), resulting form of ( ), denoted here by means of the symbol ( ), is given by ( ) will be referred to as the Baskerville [51] form of ( ).
In case of a nonnormally distributed residual error term , Newman [52] asserts that ( ) must take the form provided by Duan's smearing estimate of bias [54]. Here, this form is correspondingly represented by means of the symbol ( ), and it is calculated by means of with standing for the − ℎ residual of the contemplated regression model. Nonetheless, this form of ( ) could produce bias overcompensation ( [53,55]). Moreover, ( ) corresponds to the sample mean of retransformed residuals. Therefore, whenever outliers occur, the actual central tendency of retransformed residuals data could not be represented by ( ). Under such a circumstance, a suitable form of the correction factor ( ) seems indefinite. Moreover, Zeng and Tang [57] suggest a distribution-free form of ( ) represented here by means of ( ) and given by We notice that ( ) is actually an approximation to ( ) as it corresponds to a three terms' partial sum of the power series expression of ( ) assuming ( ) = 0. By the same token, we can consider an alternative approximation for ( ); this is represented through the symbol ( ) and given by anterms partial sum of the series representation of ( ); that is, The value of leading to the highest concordance correlation coefficient value between ( | ) projections and observed values sets the explicit form of ( ) in (C.4).

D. Fitting Results of Regresion Models in Geometrical Space
This appendix presents fitting results of the geometrical space protocols derived from the generalized model of (3). This includes TAMA and biphasic and polynomial schemes. Identification method contemplated here relies on finding parameter estimates that maximize the log-likelihood function ( , , ). This is given by Fitting results of the linear model of (D.1) through (D.3) to raw data appear in Table 12. Figure 11(a) displays a biased distribution of residuals around the zero line. The normal Q-Q plot in Figure 11(b) can be considered as a visual test of goodness of fit. In this case, this provides primary evidence against normality of residuals. Indeed, we can learn of a heavy-tails pattern in this plot. Moreover, an Anderson-Darling [67] goodness-offit test resulted in a test statistic of 246.7 and in a p value of <2.2e-16, which confirms lack of normality of residuals. Figure 12 shows spread about TAMA's mean response ( | ) as produced by different forms of the correction factor ( ). Plot suggests moderate bias about ( | ) for leaves of small-to-medium-sized leaves. Nevertheless, deviations between projected and processed data are notorious for larger area values. Thus, spread in Figure 12 suggests that the TAMA protocol is unsuited for the analysis of present data. and variance 2 . Therefore, the associated log-likelihood function ( 1 , 2 , 1 , 2 , , ) becomes

D.2. Fitting
(D.7) Table 13 Presents fitting results of the biphasic regression model of (D.4) through (D.7). It is worth mentioning that in this case we are performing a nonlinear fit in geometrical space.
Comparing the behavior of the residuals of the TAMA and biphasic models from Figure 13, we can learn of a relative stabilization of variability of the residues. Nevertheless, the biphasic normal Q-Q plot still shows heavier tails than those expected for a normal distribution. This amounts to primary evidence against normality of the residuals. An Anderson-Darling goodness-of-fit test [67] that resulted in a test statistic of 254.5 and a p value <2.2e-16 confirms a lack of normality for the involved residuals. Nevertheless, compared with the TAMA fit, the normal Q-Q plot for the biphasic model reveals a larger region where data track a normal distribution pattern.  (Table 1). Black lines associated with ( | ) green lines go with ( | ), and red ones go with ( | ) and in yellow we show those for ( | ). We can learn of a biased spread of observed values about the mean response functions ( | ) even after data processing. where the additive errors 1 , 2 , . . . , are independent and identically distributed, with distribution (0, 2 ). Hence, the response V above is normally distributed having mean ( , ) with standing for the parameter vector ( 0 , 1 , . . . , ); namely, (D.9)     Table 14.
Reviewing fitting results of the polynomial model, we observe similar results to those we found earlier for the fit of the biphasic model. Moreover, Figure 14 displays similar patterns to those shown in Figure 13. Also for this fit, an Anderson-Darling normality test [67] delivered a test statistic 257.1 and a p value <2.2e-16 that yields remarkable  Figure 14: Residuals and normal probability plots ((a) and (b) one to one) for the fitting of the polynomial regression model of (D.8) through (D.10). Compared with TAMA results, distribution of residuals about the zero line improved. Normal Q-Q plot reveals a large plateau, where residuals conform to a normal distribution, but anyhow heavier tails than expected for such a pattern remain. evidence against normality of residuals. Again, comparing with TAMA results, distribution of residuals about the zero line improved. Moreover, the Q-Q plot shows a larger plateau, where residuals agree to a normal distribution pattern.

E. Fitting Results of Direct Nonlinear Regression Models
This appendix presents the exploration of curvature by means of direct nonlinear regression methods. For that aim, we present fitting results of the nonlinear heteroscedastic model of (14). For comparison, we include fitting results of the homoscedastic case of (15).
Since we deal with data pairs ( , ) for the heteroscedastic case, the regression equation to be considered turns out to be where is a normally distributed zero mean random variable with standard deviation . Then the response is normally distributed having mean function ( , , ) given by And a variance set as a function 2 ( ) of covariate is, namely, Again, fitting this regression model relied on a likelihood approach, that is, obtaining estimates for the parameters , , , and , which result in maximizing the loglikelihood function ( , , , ) expressed by It can be ascertained from Table 15 that estimates for the parameters and for raw data are very similar to , as well as , for the heteroscedastic regression model of (E.1) and (E.4), respectively. Table 15 show some overlap. Therefore, it can be established that parameter estimates differ but just barely. This statement is supported by the estimated value of , Lin's concordance correlation coefficient between projected and observed leaf biomass values. For both regression models, homoscedastic and heteroscedastic, the estimate was = 0.972. Likewise, plotting the mean response curves = and = on a dispersion diagram of weight and area values reveals that curves corresponding to parameter estimates fitted considering heteroscedasticity or not are very similar ( Figure 5). Nevertheless, the heteroscedastic fit clearly identifies for the variance a variation pattern that grows in a power function as leaf area increases. Indeed, Figure 15 shows remarkable differences in scatter patterns of dry weight residuals against leaf area for the heteroscedastic and homoscedastic models ((a) and (b), resp.). This can be ascertained from values in Table 6, unveiling that the heteroscedastic model is selected by most agreement indices.

Confidence intervals in
Nevertheless, normal Q-Q plots of residuals of the fits of the heteroscedastic model shown in Figure 16 reveal that in spite of data cleaning we can be aware of a severe problem of heavy tails. This points to consideration of a different error structure from the normal one assumed here.

Data Availability
Data will be provided by the corresponding author upon usage agreement.

Disclosure
This research was completed while Enrique Villa-Diharce (CVU:208964) was on sabbatical at CICESE.

Conflicts of Interest
The authors declare that they have no conflicts of interest.  Figure 16: Residuals and normal probability plots ((a) and (b) one to one) for the fitting of the nonlinear heteroscedastic regression model of (E.1). We can observe a slightly biased distribution of residuals around the zero line. Also, normal Q-Q plot in (b) displays heavier tails than those expected for a normal distribution.