Assessment of a Takagi–Sugeno-Kang fuzzy model assembly for examination of polyphasic loglinear allometry

Background The traditional allometric analysis relies on log- transformation to contemplate linear regression in geometrical space then retransforming to get Huxley’s model of simple allometry. Views assert this induces bias endorsing multi-parameter complex allometry forms and nonlinear regression in arithmetical scales. Defenders of traditional approach deem it necessary since generally organismal growth is essentially multiplicative. Then keeping allometry as originally envisioned by Huxley requires a paradigm of polyphasic loglinear allometry. A Takagi-Sugeno-Kang fuzzy model assembles a mixture of weighted sub models. This allows direct identification of break points for transition between phases. Then, this paradigm is seamlessly appropriate for efficient allometric examination of polyphasic loglinear allometry patterns. Here, we explore its suitability. Methods Present fuzzy model embraces firing strength weights from Gaussian membership functions and linear consequents. Weights are identified by subtractive clustering and consequents through recursive least squares or maximum likelihood. Intersection of firing strength factors set criterion to estimate breakpoints. A multi-parameter complex allometry model follows by adapting firing strengths by composite membership functions and linear consequents in arithmetical space. Results Takagi-Sugeno-Kang surrogates adapted complexity depending on analyzed data set. Retransformation results conveyed reproducibility strength of similar proxies identified in arithmetical space. Breakpoints were straightforwardly identified. Retransformed form implies complex allometry as a generalization of Huxley’s power model involving covariate depending parameters. Huxley reported a breakpoint in the log–log plot of chela mass vs. body mass of fiddler crabs (Uca pugnax), attributed to a sudden change in relative growth of the chela approximately when crabs reach sexual maturity. G.C. Packard implied this breakpoint as putative. However, according to present fuzzy methods existence of a break point in Huxley’s data could be validated. Conclusions Offered scheme bears reliable analysis of zero intercept allometries based on geometrical space protocols. Endorsed affine structure accommodates either polyphasic or simple allometry if whatever turns required. Interpretation of break points characterizing heterogeneity is intuitive. Analysis can be achieved in an interactive way. This could not have been obtained by relying on customary approaches. Besides, identification of break points in arithmetical scale is straightforward. Present Takagi-Sugeno-Kang arrangement offers a way to overcome the controversy between a school considering a log-transformation necessary and their critics claiming that consistent results can be only obtained through complex allometry models fitted by direct nonlinear regression in the original scales.


INTRODUCTION
Julian Huxley introduced the theory of constant relative growth between a trait y and overall body size x (Huxley, 1924;Huxley, 1932;Strauss & Huxley, 1993). This paradigm is commonly refereed as Huxley's model of simple allometry and is essentially formulated through the power law y = βx α with α identified as the allometric exponent and β as the normalization constant. In biology allometric relationships are to within species, as well as, between species (evolutionary allometry) (Houle et al., 2011;Marquet et al., 2005;West & Brown, 2005;Pélabon et al., 2014). Power function models are also extensively used in other research fields, e.g., physics (Newman, 2007), ecology (Harris, Duarte & Nixon, 2006;Hood, 2007) earth and atmospheric sciences (Hills, 2013), and economics (Li et al., 2015). This has encouraged many research endeavors addressing interpretation of involved parameters, as well as, suitability of analysis method for getting estimates. Indeed, concomitant to Huxley's theory of relative growth is the Traditional Analysis Method of Allometry (TAMA hereafter). This is, a widespread device to acquire estimates of the parameters α and β. It contemplates a logarithmic transformation of the original bivariate data in arithmetical scale in order to consider a linear regression model in geometrical space, and then retransforming to acquire Huxley's model of simple allometry in the original scale. This approach implicitly embraces a notion that variability of the response conforms to a pattern of multiplicative growth. On Huxley's elucidation (Huxley, 1932) the intercept lnβ of TAMA's line was of no specific biological importance, but the slope b was significant enough as to mean allometry itself. This interpretation has permeated contemporary research to such an extent that many practitioners still consider it to be the valid theoretical perspective for static and ontogenetic allometry (Eberhard, 2009;Houle et al., 2011;Pélabon et al., 2018). However, views assert that a TAMA approach produces inconsistent results, thereby recommending allometric examination by relying instead on nonlinear regression in the direct scales of data (Packard, 2017a;Packard, 2013;Packard, 2009;Packard & Birchard, 2008). This is, debatable for defenders of the traditional approach that claim that, as it is conceived in the original theoretical context of allometry, a logarithmic transformation deems necessary in the analysis (Lai et al., 2013;Klingenberg, 1998;Nevill, Bate & Holder, 2005;Kerkhoff & Enquist, 2009;Xiao et al., 2011;White et al., 2012;Ballantyne, 2013;Glazier, 2013;Niklas & Hammond, 2014;Lemaître et al., 2015;Pélabon et al., 2018). Yet steering further away from Huxley's perspective on covariation among different traits, other views conceive allometry centered on the covariation between size and shape (Mosimann, 1970;Klingenberg, 2016).
From this standpoint, analysis must rely in Multiple Parameter Complex Allometry (MPCA after this) formalizations through all varieties of nonlinear or discontinuous relationships (e.g., Frankino, Emlen & Shingleton, 2010;MacLeod, 2014;Bervian, Fontoura & Haimovici, 2006;Lovett & Felder, 1989;Packard, 2013). However, adoption of MPCA approaches nourishes one of the most fundamental discrepancies among schools of allometric examination. Indeed, for advocates of the traditional approach, examination based on MPCA models fitted in arithmetical scale sacrifices appreciation of biological theory in order to privilege statistical correctness (Houle et al., 2011;Lemaître et al., 2015;Pélabon et al., 2018). A way to keep the analysis in geometrical space while amending unreliability of a linearity assumption is conceiving the notion of non-log linear allometry (Packard, 2012b;Strauss & Huxley, 1993;Echavarría-Heras et al., 2019a). As every analytic function can be expanded as a power series, curvature in geometrical space has been addressed through polynomial regression schemes (Kolokotrones et al., 2010;Lemaître et al., 2014;MacLeod, 2010;Glazier, Powell & Deptola, 2013;Tidière et al., 2017;Echavarría-Heras et al., 2019a). But besides difficulties related to biological interpretation of a polynomial mean response, this approach maintains a single functional form of the response over the whole covariate range. This could not account for inherent heterogeneity in the logtransformmed response as contemplated in Huxley's theoretical perspective. Certainly, Huxley reported a breakpoint in the log-log plot of chela mass vs. body mass of fiddler crabs (Uca pugnax). It was attributed to a sudden change in relative growth of the chela approximately when crabs reach sexual maturity (Huxley, 1924;Huxley, 1927;Huxley, 1932). This suggests a slant aimed at adding complexity in geometrical space while keeping the theoretical essence of traditional allometry in the analytical set up. Is this conception that hosts polyphasic loglinear allometry approaches (PLA afterwards) (Packard, 2016;Gerber, Eble & Neige, 2008;Strauss & Huxley, 1993;Hartnoll, 1978). PLA characterizes heterogeneity of the logtransformmed response by composing covariate range into sectors separated by break points. Each subdivision associates to a linear sub model. Brokenline regression (Beckman & Cook, 1979;Ertel & Fowlkes, 1976;Tsuboi et al., 2018;Ramírez-Ramírez et al., 2019;Muggeo, 2003;Echavarria-Heras et al., 2019b). Forbes & López (1989) furnish an empirical approach to identification of PLA patterns. Nevertheless, by relying in nonlinear regression this technique requires starting values for the break-point estimation. Therefore, complications set by local maxima, as well as, inferences on estimates could make implementation difficult (Julious, 2001;Muggeo, 2003).
An operating regime based modeling approach offers a structure supporting model adaptation amid an empirical and mechanistic standpoint. Local models valid over restricted domains are combined by smooth interpolation into an overall general output (Johansen & Foss, 1997). Therefore, this structure naturally hosts heterogeneity as conceived by PLA (Echavarria-Heras et al., 2019b). One example of a hybrid-operating regime based modeling is the Takagi-Sugeno-Kang fuzzy model (Sugeno & Kang, 1988;Takagi & Sugeno, 1985) (TSK in what follows). This construct composes a fuzzy logic step intended to characterize smoothing weight factors. Then, conventional statistical methods are used to acquire estimates of parameters characterizing sub models. It turns out that the general output of a first order TSK fuzzy model can uniformly approximate any continuous function to arbitrarily high precision (Ying, 1998;Zeng, Nai-Yao & Wen-Li, 2000). As we show in this examination, an advantage of TSK over conventional PLA, is that it can offer convenient non-statistical proxies of break points for transition among phases. Moreover, consideration of sub models of a TSK scheme as TAMA's linear functions in geometrical space not only offers a congruent PLA model, but it could also entail a highly biologically meaningful model of allometry, because it can model the breakpoints while keeping the meanings of allometric exponents as in Huxley's original formulation. A comprehensive exploration of suitability of the TSK scheme to examine PLA patterns has not been undertaken so here we attempted to fill this gap. In what follows a formulation of PLA by means of a the TSK fuzzy model will be referred as TSK-PLA for short.
The outstanding approximation capabilities of a TSK fuzzy model entail reliable identification of whatever MPCA functional form renders necessary in arithmetical space (Echavarría-Heras et al., 2018a;Echavarria-Heras et al., 2019b). Adaptation of the TSK fuzzy model for that aim will be forward designated by means of the TSK-MPCA abbreviation. As a criterion to evaluate the performance of the TSK-PLA proxy we verified the dependability of linked retransformation results, including break point placement and reproducibility strength of mean response function against corresponding estimations produced via TSK-MPCA. It turns out that proposed TSK-PLA analysis method endorsed reliable identification of heterogeneity of examined allometries. Furthermore, the affine structure of the present fuzzy protocol can accommodate either complex or simple allometry as required to analyzing the data. Thus, the presented TSK-PLA model can be considered as a general tool for examination of zero intercept allometries. Moreover, from a theoretical standpoint a TSK-PLA representation implies an allometric model in arithmetical space that seemingly fits MPCA. This expresses the response as a generalized power function including scaling parameters expressed as functions of the covariate (Bervian, Fontoura & Haimovici, 2006;Echavarría-Heras et al., 2019a;Echavarria-Heras et al., 2019b). But, above all, present fuzzy approach contributes by offering a way of overcoming the controversy between a school considering analysis in geometrical space as a must in allometry, and critics claiming that consistent results can only come along by using a MPCA formulation followed by nonlinear regression protocol in the original scale of data. Interestingly, present TSK-PLA arrangement also contributed on qualitative grounds. Certainly, Huxley reported a breakpoint in the log-log plot of chela mass vs. body mass of fiddler crabs (Uca pugnax). Packard (2012a) inferred this point was only putative. In his own interpretation, perhaps due to combined effects of a log transformation itself and the format of graphical display of Huxley's data. However, application of present Takagi-Sugeno-Kang protocol supports existence of a break point in Huxley's Uca pugnax log-log plot.
This article is organized as follows: In the Materials and Methods section, we formally explain the steps backing the identification of the offered TSK-PLA scheme. There, we explain why this construct can be considered as a generalized protocol for allometric analysis in geometrical space. We also clarify why the offered TSK model implies a MPCA scheme in arithmetical space. The presentation includes an elucidation of sufficient conditions under what the asymptotic mode of the acquired TSK proxy behaves as the power function in Huxley's model of simple allometry. There, we also suggest a correction factor (CF here after) for bias of retransformation of the regression error that grants highest reproducibility for derived mean response function in arithmetical space. The Results section highlights on the advantages of the present approach over conventional counterparts. A Discussion section elaborates on the contribution that our approach bears for the general subject of suitability of analysis method in bivariate allometry. An Appendix includes a detailed explanation of the steps involved in the construction and identification of the general form of the addressed TSK models.

Data
Allometric examination here mainly relied on a primary data set exhibiting curvature in geometrical space. This composes 10,412 measurements of Zostera marina (Eelgrass) leaf biomasses y and corresponding leaf areas x as reported in Echavarría-Heras et al. (2019a) and Echavarría-Heras et al. (2018b). For comparison, we also considered data reported in Mascaro et al. (2011) comprising 30 Biomass-Diameter at Breast Height measurements on Metrosideros polymorpha. Analisis also extended to data reported in De Robertis & Williams (2008) including 29,363 Length-Weight measurements on Gadus chalcogrammu. This last data set allowed illustration of the performance of the TSK paradigm in a circumstance where the TAMA protocol is consistent. Finally, we analysed the fitness of the TSK in detecting break points in the log-log plot of chela mass vs. body mass of fiddler crabs (Uca pugnax) (Huxley, 1924;Huxley, 1932).

General formula of allometry
The methods engaged here aim to identification of the suitable form of the allometric function representing the variation of a trait y depending on a descriptor x. For that purpose, we firstly introduce the formal framework and the notation convention used through. We assume that a response y and its covariate x belong to domains Y and X of positive real numbers one to one and with y having a zero limit when x approaches zero. We also consider that there exist a function w x,p : X → Y where p = (p 1 ,...,p n ) is a parameter set, and a concomitant approximation error function (x) : X → Y that combine to model whatever form, the linkage between x and y acquires. Moreover, we take on, that such a relationship can be expressed through an additive error description or else through the multiplicative error alternate In order to get w x,p , we can consider the error term (x) as a random variable . Then, specifications above offer two commonly addressed analysis protocols in allometry. A regression model with additive error in arithmetical scale with taken as ψ−distributed with zero mean and variance generally expressed as a function σ 2 (x) of covariate, that is, ∼ ψ 0, σ 2 (x) . Fitting Eq. (3) generally requires direct nonlinear regression protocols. This returns a mean response function E aw y|x = w(x,p). For the sake of facilitating comparison aims in further developments, this a subscript will be maintained to typify a mean response function gotten by means of identification protocols applied in arithmetical space. A second procedure circumscribes to the multiplicative error model of Eq.
(2) and relies in a logtransformation procedure in order to consider a parallel regression model in geometrical space. Formally, we contemplate a mapping (y,x) → (v,u) such that u = lnx and v = lny. This sets variation domains U and V for u and v to one. We concomitantly have the regression model with additive error in geometrical space where formally v (u,π) = ln w(x,p) and is random variable as specified above. It follows that back-transforming Eq. (4) to arithmetical space yields, Then, concomitant mean response function is symbolized through E gw (y|x) and becomes where δ = E (e ). Notice that in E gw (y|x) we have used the notation convention of a subscript g referring to identification of w(x,p) based on the regression model of Eqs. (4) and (5). The CF, δ above provides the necessary adjustment for bias of retransformation of the regression error (Mascaro et al., 2011;Baskerville, 1972;Newman, 1993). Assuming ∼ N 0,σ 2 sets e to be lognormally distributed. Then, CF becomes But, Newman (1993) asserts that whenever is not normally distributed, δis given by the smearing estimate of bias of Duan (1983). Nevertheless, in some settings this nonparametric form can produce bias overcompensation (Manning, 1998;Smith, 1993;Koch & Smillie, 1986). Zeng & Tang (2011a) propose an alternate nonparametric form of δ namely Actually, δ given this way corresponds to a three terms partial sum approximation of the power series expression of E(e ) assuming E ( ) = 0. By the same token, Echavarría-Heras et al. (2019a) suggest a representation for δ given by a n-terms partial sum of series representation of E (e ), that is, Maximization of Lin's Concordance Correlation Coefficient (CCC) (Lin, 1989) between observed values and mean response E g (y|x) resulting using this form of δ sets criterion to choose n.

Huxley's formula of Simple Allometry
A characterization of w x,p as a power function βx α has been traditionally referred as Huxley's formula of simple allometry (Strauss & Huxley, 1993). This model will be ahead epitomized by a subscript s as a mnemonic device for ''simple''. Equation (3) becomes with w s x,p = βx α and assumed to be normally distributed with zero mean and variance σ 2 , that is, ∼ ψ 0, σ 2 . According to our notation convention Eq. (11) yields the mean response function E as y|x = βx α . Similarly, the logtransformation method produces the TAMA's regression model, that with v s (u,π) = lnβ + αu (13) and as specified above. Equations (12) and (13) determine E s (v|u) = v s (u,π). Accordingly, back transformation of Eq. (12) to arithmetical space brings up a mean response E gs y|x given by where δ stands for necessary CF. It often occurs that even after contemplation of proper form for δ this TAMA's derived proxy for w x,p produces biased projections of observed values of the response y. This means, that complexity of Huxley's formula of simple allometry w s x,p becomes inappropriate to identify the true w x,p form (Gould, 1966;Huxley, 1932;Bervian, Fontoura & Haimovici, 2006;MacLeod, 2014;Echavarría-Heras et al., 2019a). From the settings of Eq. (1) it is reasonable assuming that adapting complexity as it is needed to identify w x,p could depend on MPCA forms. Corresponding logtransformmed expression v (u,π) is inferred to be a nonlinear function of covariate u. This rears PLA as a likely device to acquire complexity for identification of MPCA through geometrical space methods. We adopt the affine structure of a first order TSK fuzzy model as a device for identification both MPCA or PLA alternates.

The TSK account of w x,p
The general output of a first order TSK fuzzy model bears a fuzzy alternate to a statistical mixture regression model (Cohn, Ghahramani & Jordan, 1997). It is then reasonable to assume that such an structure could efficiently address the problem of identifying w x,p expressed as a MPCA model in arithmetical scale or its assumed PLA forms in geometrical space. The symbol w TSK x,p will stand for the TSK-MPCA surrogate for w x,p . Accordingly, adaptation of Eq. (3) becomes with TSK a ψ−distributed residual random variable with zero mean and variance expressed as a function σ 2 TSK (x) of x, that is, TSK ∼ ψ(0,σ 2 TSK (x)). Denoting through the symbol E aTSK y|x the mean response function determined by Eq. (15), we have Since, the general output of a first order TSK fuzzy model can uniformly approximate any continuous function to arbitrarily high precision (Ying, 1998;Zeng, Nai-Yao & Wen-Li, 2000) then whatever MPCA form w x,p embraces, this can be accurately projected through a consistent identification of E aTSK y|x .
By assumption, we take v TSK (u,π) in the form given by Eq. (A14), that is, with ϑ i (u) and f i (u) a one to one the firing strengths and consequent functions for i = 1,2,...,q. Since, the domain U of covariate is R, we can assume membership functions µ k (u) to have a Gaussian form i.e., being θ k and λ k for k = 1,2,...,q, parameters. We also consider that consequent local models f i (u) have a description, that is,  (u,π) can be interactively explored through different characterizations of the clustering radiustraining parameter r a as specified by Eqs. (B7) through (B9).
As described in Appendix A, acquiring v TSK (u,π) requires on first stage a fuzzy partition L u of the input domain U (cf. Eq. (A3)). Achieving this relies on a Subtractive Clustering (SC after this) technique to establish the value of the parameter q (Castro et al., 2016;Chiu, 1994). A brief description of the SC method is provided in Appendix B. This stage also sets the number of inference rules R i specified by Eq. (A10) and concomitant number local models in v TSK (u,π). The SC step also produces estimates for the parameters θ k and λ k in characterizing the membership functions µ k (u). Then, the normalized firing strength functions, ϑ i (u) follows from (Eqs. A11) and (A12). A second step targets at characterization of the linear consequents f i (u) as given by Eq. (22). This is achieved by replacing the identified factors ϑ i (u) and the assumed form of the consequents f i (u) into Eq. (20) to characterize the regression model of Eq. (17). Then, the parameters in the consequents f i (u) are estimated by implementing a Recursive Least Squares (RLS) routine (Jang, Sun & Mizutani, 1997;Wang & Mendel, 1992). This identification step could be also performed through a maximum likelihood approach (Kalbfleisch, 1985). The whole identification algorithm is explained in Appendix B. Codes are included in the Supplementary Information section. Following Echavarria-Heras et al. (2019b) reproducibility will be primarily estimated by comparing values of Lin's concordance correlation coefficient, symbolized by means of ρ C (Lin, 1989). Agreement will be defined as poor whenever ρ C < 0.90, moderate for 0.90 ≤ ρ C < 0.95, good for 0.95 ≤ ρ C < 0.99, or excellent for ρ C ≥ 0.99 (McBride, 2005). Assessment of reproducibility will also rely on model performance metrics, such as the Coefficient of Determination (CD), Standard Error of Estimate (SEE), Mean Prediction Error (MPE), and Mean Percent Standard Error (MPSE) (Gupta, Sorooshian & Yapo, 1998;Hauduc et al., 2011;Zeng et al., 2017;Zeng & Tang, 2011b;Parresol, 1999;Meyer, 1938;Schlaegen, 1982). Related statistics are included in Appendix C. Matlab and R codes are provided in the supplemental files section. Fig. 1 display the spread response-covariate in geometrical space for data sets included in this examination. Figure 1A relates to the Echavarría- Heras et al. (2019a) data. Figure 1B is for Mascaro et al. (2014) data. Figure 1C shows Huxley's plot of chela mass vs. body mass of fiddler crabs (Uca pugnax) in log-log scale (Huxley, 1924;Huxley, 1932). Figure 1D displays spread for the De Robertis & Williams (2008) data. Assessment of curvature will be performed for all data sets by analyzing fitting results of the TSK-PLA and TSK-MPCA models. For easy of presentation detailed results on a TAMA-TSK model comparison will only circumscribe to the Echavarría- Heras et al. (2019a) data.

Representation of the back-transformed TSK-PLA proxy as a MPCA formula
This section explains that assuming TSK-PLA implies a multiple parameter complex allometry form in direct arithmetical scales. Indeed following Bervian, Fontoura & Haimovici (2006), Echavarría-Heras et al. (2019a) proposed an extension of Huxley's formula of simple allometry w s a,p = βa α that includes scaling parameters α and β depending in the covariate, that is, where β (x) and α(x) are continuous functions and with β (x) assumed to be positive. This sets Thus, formally whenever the scaling functions β (x) and α(x) are not simultaneously constant w x,p as given by Eq. (24) entails MPCA (Echavarría-Heras et al., 2019a;Echavarria-Heras et al., 2019b). We now demonstrate that the mean response function E gTSK (y|x) in arithmetical space derived from a TSK-PLA arrangement implies MPCA in the form set by Eq. (24). For that drive, we notice that replacing Eq. (22) into Eq. (20) and then rearranging leads to v TSK (u,π) = lnβ TSK (u) + α TSK (u)u where Thus, Eq. (17) takes on the equivalent representation, The functions lnβ TSK (u) and α TSK (u) above, suggest u − dependent forms of the parameters lnβ and α involved in the regression model of the TAMA approach. Then, under the assumption of Eq. (22) a TSK-PLA arrangement interprets as generalization of the TAMA scheme (Echavarría-Heras et al., 2019a). Applying the back-transformation e v of Eq. (28) and recalling Eq. (19) yields This sets exp(v TSK (u,π)) = β TSK (x)x α TSK (x) . But, from Eq. (5) we have v TSK x,p = ln(w TSK x,p ) then w x,p as identified by retransformation of v TSK x,p admits the form specified by Eq. (24).  (12)). Table 1 summarizes the results of the analysis. Corresponding, mean response E s (v|u) in geometrical scale is shown in Fig. 2A. Log-transformation is a mechanism aimed to reduce variability of data (Feng et al., 2014). Nevertheless, Fig. 2A still displays a significant dispersion of v values about E s (v|u). Spread may lead on first sight to the impression that the distribution of v around the mean response line E s (v|u) for small values of u is fair. Agreeing with (Packard, 2017b), on the importance of graphs in allometry, led to a careful revision of the spread which suggested curvature. Moreover, the assessment of dispersion of residuals of Eq. (12) suggested lack of normality, as well as, heteroscedasticity (Fig. 2B). Further, QQ plot shows heavier tails than expected for a normal distribution (Fig. 2C). Indeed, an Anderson & Darling (1952) test to ascertain normality of residuals produced a test statistics value of A = 310.848 and a : p-value < 2.2e-16. In turn a Lilliefors (Kolmogorov-Smirnov-type) test, delivered D = 0.1305, as well as, a relatively small p-value < 2.2e-16. Therefore, the hypothesis of normally distributed errors in the analysis should be rejected since obtained p-values are extremely small (< 2.2e-16). What is more, a lack of normality of errors in the linear regression analysis of Eq. (12) can be also ascertained from the normal QQ plot shown in Fig. 2C. It can be perceived that the distribution of residuals exhibits heavier tails than those expected for a normal distribution. Besides, a Breush-Pagan statistic (Breusch & Pagan, 1979), provided a way to assess heteroscedasticity of the residuals. In order to perform this test, the squared errors in the linear model of Eq. (12) were assumed to depend linearly on the independent variable i.e., where b and d are parameters and ζ the error term. The null hypothesis is that the parameter d in Eq. (30) vanishes. Rejection of the null hypothesis not only corroborates heteroscedasticity but also provides information on variability. The test statistics turned out to be BP = 808.8119 with one degree of freedom with a p − value (< 2.2e-16), that is sufficiently small as to provide strong evidence against homoscedasticity, while undoubtedly favoring heteroscedasticity. Thus, the presently fitted straight line does not comply the essential assumptions of normality and homoscedasticity of the analysis. Therefore, the TAMA protocol turned out unsuited for analyzing the allometric relationship in the Echavarría-Heras et al. (2019a) data. Consequently, characterization of the general function w(x,p) entailed by Huxley's model of simple allometry (cf. Eq. (11)) does not fit required complexity. Thus data spread shown in Fig. 2A, submits curvature in geometrical space. We now turn to explore the capabilities of the TSK-PLA construct to identify this pattern.
Fitting results of the TSK-PLA assembly: Zostera marina

Identification of firing strength factors ϑ i (u)
In order to identify the required firing strength factors ϑ i (u) for i = 1,2,...,q. We executed the main_fun_tsk_pla_model_fit.m code on log-transformed values (v,u) from the Echavarría-Heras et al. (2019a) data set. This try assumed membership functions µ k (u) having a form given by Eq. (21) for k = 1,2,...,q. Setting r a = 0.47 returned a value of q = 2. Then, we have to consider two membership functions characterizing the fuzzy partition of imput space U . Moreover, in compliance with Eq. (A12) normalized firing strength factors ϑ 1 (u) and ϑ 2 (u) turn out to be Plots of the estimated membership functions µ 1 (u) and µ 2 (u) and normalized firing strength factors ϑ 1 (u) and ϑ 2 (u) appear in Fig. 3A and Fig. 3B respectively. We observe that agreeing curves intersect at a point u b = 3.998.

Identification of consequent functions f i (u)
A second step in the procedure to get v TSK (u,π) concerns acquiring the consequent functions f i (u) in Eq. (22). Since, for this data, we obtained q = 2, the code ought to establish consequent functions f 1 (u) and f 2 (u),each one assumed to be linear, that is, and With the aim of assessing heteroscedasticity, we replaced the forms of ϑ 1 (u) and ϑ 2 (u) identified by SC technique in regression Eq. (17). In turn the involved consequent functions f 1 (u) and f 2 (u) were assumed to have both the form given by Eqs. (33) and (34) correspondingly. Then, we assumed the involved TSK residuals to be normally distributed with zero mean, but with a standard deviation set as a function σ TSK (u) of the covariate u, Namely where σ and k are parameters to be estimated and such that σ + ku > 1. Table 2 presents maximum-likelihood parameter estimates and associated uncertainties for the related fit. We can ascertain a highly significant fit, since, in all cases the standard error is extremely small, this mainly explained by the large amount of data in the analysis. To judge heteroscedasticity of residuals we study the confidence interval of parameter k. This parameter determines the change in residual variability per unit change in u. It turns out that figures in Table 2 show that confidence interval of parameter k does not include a zero value. Therefore, we may conclude that with high probability the variability of the residuals changes as u changes. Meanwhile, setting k = 0 in Eq. (35) allowed consideration of a parallel maximum likelihood fit of homoscedastic case of the TSK regression model of Eq. (17). Table 3 provides fitting results. Model performance metrics allow assessment of the fits of the heteroscedastic and homoscedastic versions of the TSK-PLA protocol. Accordingly, we can assert that the heteroscedastic model stands a better fit than its homoscedastic counterpart. Particularly, in the heteroscedastic case we have a negative log-likelihood  value of −logL = 6304.60, which turns out to be notably smaller than −logL = 8111.49 obtained for the homoscedastic model. These figures bear a notable difference that backs the selection of the heteroscedastic model. This difference in fit quality favoring the heteroscedastic model is mainly due to the fact that the latter models the pattern of variation of the errors in a more reliable way. Plots of identified consequents appear in Fig.  3C, component products ϑ 1 (u)f 1 (u) and ϑ 2 (u)f 2 (u) appear in Fig. 3D. As it occurs for the membership functions and firing strength factors for this data the component products also intersect at value of u b = 3.98. This estimate of u b is consistent with value previously reported by Echavarría-Heras et al. (2019a) for this data and deriving from conventional maximum likelihood biphasic regression. Figure 4A displays dispersion about resulting mean response function v TSK (u,π). Figure 4B shows residual scattering about the zero line. Region bounded by red lines determine (95%) confidence intervals. Figure 4C shows corresponding QQ plot.

Comparison TAMA vs. TSK-PLA
Compared with corresponding fitting results for the TAMA protocol ( Fig. 2A) we can verify that plots in Fig. 4 show fairer residual spread patterns. Nevertheless, the QQ plot in Fig. 4C, still suggest deviation of TSK residuals from a normal distribution pattern. Table 4 allows further comparison of performances of the TAMA and TSK proxies. This undoubtedly favor selection of the TSK scheme. Therefore, opposed to the linear model v s (u,π), the affine characterization of variability granted by v TSK (u,π) can better refer to inherent non-log linear allometry for the Echavarría-Heras et al. (2019a) data set. Certainly, the point u b = 3.98 shown in Fig. 3B can be interpreted as a point separating two phases describing the variation pattern of the log transformed response v. One for small leaves valid over the region u < u c and another for large leaves over u ≥ u b . The form of the component products ϑ i (u)f i (u) shows that while u drifts away from u b taking smaller and smaller values the closer the TSK output v TSK (u,π) will be to the component product ϑ 1 (u)f 1 (u). Conversely, the larger the distance between u and u b for leaves in the large phase u ≥ u b the closer v TSK (u,π) will be to ϑ 2 (u)f 2 (u). Relating to E s (v|u) shown in Fig. 2A, we can assess from  Fig. 4C shows a larger plateau where TSK residuals track a normal distribution pattern. Nevertheless, application of an Anderson & Darling (1952) test to the residuals of regression Eq. (17) resulted in AD = 370.17. This associates a p-value < 2.2 × 10 −16 , that provides evidence against a normality assumption for the TSK residuals. This is, in agreement with the observation that he normal Q-Q plot shown in Fig. 4C showing heavier tails than those expected for a normal distribution. It is worth pointing out that the break point u b identified by the fuzzy proxy v TSK (u,π) coincides with corresponding value obtained by Echavarría-Heras et al. (2019a) using conventional biphasic regression methods. Correspondingly, Fig. 5A displays the plot of the estimated form of the mean response function E gTSK (x|y) of Eq. (19). Since, residuals TSK are not normally distributed, Eq. (10) provided CF form. Figure 5B allows a visual assessment of the extent of biased projections in arithmetical scale tied to the TAMA surrogate E gs (x|y) calculated with Duan's form of δ. Compared with spread deriving from the TSK model, TAMA's bias is significant. Besides, Table 4 allows assessment of differences in associated predictive strengths. All indices favor the TSK-PLA scheme. As suggested by perceptible bias shown in Fig. 5B, CCC value for TAMA's projections point to poor reproducibility of observed values. Besides, relevance of accounting for curvature, this assessment highlights on the importance of choosing a proper form of δ for assuring consistency or retransformation results.

Asymptotic analysis of TSK-PLA assembly
In this section we explain that adoption of a TSK-PLA approach allows exploration of asymptotic behavior of the allometric mean response function in arithmetical space. After replacing f 1 (u) and f 2 (u) as given by Eqs. (33) and (34) in Eq. (20) Similarly, it can be directly verified that firing strengths ϑ 1 (u) and ϑ 2 (u) as given by Eqs. (31) and (32) can be also expressed in the form and where with θ,λ standing for parameter vectors (θ 1 ,θ 2 ) and (λ 1 ,λ 2 ) one to one. We can then ascertain that φ (u,θ,λ) remains positive for all values of u. Also, the term ξ (θ,λ),does not depend on u. Consequently, whenever the factor ψ (λ) in Eq. (40) is positive, τ (u,θ,λ) will approach infinity as u approaches infinity. Then, Eq. (37) implies ϑ 1 (u) asymptotically vanishing as u approaches infinity. Reversely, whenever ψ (λ)is negative, the firing strength factor ϑ 1 (u) will asymptotically approach one as u approaches infinity.
We denote by means of the symbol, E ∞ gTSK y|x the limit of E gTSK y|x as x approaches infinity, that is, Then, Eqs. (31), (32) and (44) through (46) imply Then, the asymptotic mode E ∞ gTSK y|x identified for the Echavarría-Heras et al. (2019a) data set, is attains a form like Huxley's formula of simple allometry w s x,p . Estimated parameters are α 2 = 1.1126 and β 2 = 7.8398 × 10 −6 . Figure 5C displays observed leaf biomass values y and their projections through the E ∞ gTSK y|x proxy. We can learn of a remarkable correspondence between the power function w s x,p = βx α of Eq. (11) fitted by direct nonlinear regression methods and the asymptotic mean response E ∞ gTSK y|x . Besides as established by Eq. (45)  fair distribution about the zero line. And, normal QQ-plot in Fig. 7C shows a large plateau where residuals track a normal distribution pattern. We can also ascertain from goodness of fit statistics in Table 5, that compared to the linear regression scheme of the TAMA protocol, the affine modeling approach composing the TSK-PLA scheme entails consistent identification of curvature in geometrical space.

Identification of the TSK-PLA proxy: Uca pugnax
Huxley conceived a breakpoint in the log-log plot of chela mass vs. body mass of fiddler crabs (Uca pugnax) (Huxley, 1924;Huxley, 1932). Huxley situated this point between the 15th and 16th observations and assumed it meant a to a sudden change in relative growth of the chela approximately when crabs reach sexual maturity. Examination of Huxley's data by Packard (2012a) implied such a break point to be only putative and in Packard's own interpretation, perhaps explained by the fact that Huxley could have been misled by the effects of the log transformation itself, along with the format of graphical display that might have exaggerated the slopes of segments before and after the change point. In order to test the performance of the TSK-PLA protocol in analyzing Huxley's Uca pugnax data, we took averages of both body mass and chela mass form Table 1 in Huxley's report (Huxley, 1932). Concurrent log transformed values appear in Fig. 1C.
For easy of presentation a break point as conceived by Huxley's will be denoted here through the symbol u bH . One substantial advantage of the fuzzy logic approach over conventional probabilistic slants is that the former facilitates knowledge based modeling. In order to incorporate previous knowledge, we we abided by Huxley's assertion of biphasic allometry in Uca pugnax. Then, we examined heterogeneity patterns predicted by the TSK-PLA system for different values of clustering radius r a . Particularly, setting r a = 0.8 returned q = 2, arranging biphasic allometry. Acquired firing strengths appear in Fig. 8A, exhibiting a break point at u b = 5.831. Figure 8B display consequent linear functions with estimated slopes α 1 = 1.2676 and α 2 = 1.4708 one to one respectively. In the settings of performed TSK-PLA analysis these correspond to exponents characterizing allometric phases as conceived in Huxley's original theoretical standpoint. Correspondingly, Fig. 8C portrays consequent component products ϑ 1 (u)f 1 (u) and ϑ 2 (u)f 2 (u). Similarly, Fig. 8D shows spread about mean response v TSK (u,π) including placement of u b in a display in compliance with that in Fig. 3 in Huxley (1932). We can be aware that location of u b is shifted back relative to u bH . Figure 8E displays placement of u b and spread about v TSK (u,π) in the original scale of data (cf. Fig. 1C). But instead, we may integrate previous knowledge by considering for that the break point u bH actually exists. Then, we can search among different values of r a , the one for what the TSK-PLA arrangement compromises a break point u b placed as u bH and also a maximum reproducibility strength of interpolation by v TSK (u,π). Accordingly, setting r a = 0.2 brought about q = 7 sub models, inducing a maximum reproducibility strength of interpolation function v TSK (u,π) and where u bI , one of six detected break points is placing as u bH . (Fig. 8E and Table 6). Interestingly, visual examination of plot showing v TSK (u,π) suggests a pattern accommodating two linear segments that alternate about u bI . Moreover, using the interpolation points (u,v TSK (u,π)) we fit two linear segments of slopes α 1I = 1.626 and α 2I = 1.274 before and after u bI one to one (Fig. 8F). Since u bI . can be taken as a proxy for u bH the TSK-PLA interpolation mode could suggest Huxley's reasoning of biphasic allometry in in Uca pugnax as consistent. In the meantime acquired q = 7, interpolation confirms the outstanding capabilities of the TSK-PLA device to approximate the dynamics of the logtransformmed allometric response. This can be better ascertained from Fig. 9A through Fig. 9C presenting spread about the high order interpolation function v TSK (u,π), as well as, concomitant residual and QQ plots in conforming order. Moreover, Fig. 9D through Fig. 9F allow visual comparison of parallel results by a TAMA's fit. Additionally, Table 6 compares model performance metrics for the TSK-PLA interpolation and TAMA's output fits. We can ascertain that the TSK-PLA interpolation stands a better fit. In any event, the non-probabilistic interpretation of uncertainty backing the TSK-PLA approach seems to advocate biphasic heterogeneity in geometrical space for Huxley's Uca pugnax data.

Fitting results of the TSK-PLA assembly: Gadus chalcogrammu
A fit of the TSK-PLA protocol to Gadus chalcogrammu data reported in De Robertis & Williams (2008), can exhibit reliability of this paradigm in further way. Visual examination of spread in geometrical space may suggest curvature. But, setting r a = 0.5 led to highest reproducibility of v TSK (u,π) characterized in a linear form. Indeed, plots in Fig. 10 show   that identification of the TSK-PLA model for this data, produced only one membership function µ 1 (u) (Fig. 10A). This corresponds to a firing strength ϑ 1 (u) set to one (Fig.  10B), and a conforming single TAMA's form linear consequent f 1 (u) (Fig. 10C). This matched the single linear component product function shown in Fig. 10D. As a result, no heterogeneity as determined by breaking points u b was detected for this data. Consequently, the TSK arrangement suggests a fit equivalent to the TAMA approach. Moreover, spread abut mean response, residual and Normal QQ-plots for a TAMA fit performed in this data ( Fig. 11A through Fig. 11C respectively) seem to faithfully agree to corresponding plots (Fig. 11D, through Fig. 11F) associating to the TSK-PLA fit. In turn model performance metrics in Table 7 corroborate these alternate fits as equivalent. Therefore, the TSK-PLA assembly seemingly adapts required complexity. This supports judgement on this paradigm being considered as a generalized tool for allometric examination in geometrical space.

Assembly of the TSK-MPCA proxy
We assume that w x,p as intended for MPCA can be modeled by w TSK x,p as expressed by means of Eq. (A14), in arithmetical space, namely with firing strengths ϑ i (x) given by being µ i (x) for i = 1,2........q the involved membership functions. We also undertake that both w x,p and x remain positive, and that lim x→0 + w x,p = 0.  This sets the imput space X to be R + . We then contemplate that membership functions can be expressed through a composite log normal form that satisfies the constrain by Eq. (50), namely where c = (1 − e) −1 and with θ i and λ i for i = 1,2,....,q parameters. Correspondingly we consider the consequents f i (x) to be linear functions, that is, It is worth recalling that Eq. (15)

Identification of the TSK-MPCA proxy: Zostera marina
For the Echavarría-Heras et al. (2019a) data a try of the main_fun_tsk_mpca_model_fit.m code setting r a = 0.5416 returned q = 2 for a biphasic mode and a maximum reproducibility of w TSK x,p . Figure 12A displays acquired firing strength functions. The estimated break point was x b = 49.632 consistent with the retransformed value of u b = 3.9 for this data set. This means that variability of the response y indeed conforms to a MPCA pattern in the direct scale of data. Corresponding spread about fitted mean response function E aTSK y|x appear in Fig. 12B. This plot allows comparison to its counterpart E gTSK y|x produced by retransformation of mean function v TSK (u,π) to arithmetical space. We can be aware of remarkable correspondence through x values. This validates adequacy of a TSK-PLA analysis for this data. Figures 12C and 12D show residual spread and QQ-plot for the TSK-MPCA fit. Figures 12E and 12F show corresponding plots for retransformed TSK-PLA fit. Besides Table 8 allows assessment of addressed proxies through model performance metrics. This could endure a judgement that concurrent MPCA pattern in arithmetical space can be consistently characterized by retransformation of PLA results.

Identification of the TSK-MPCA proxy: Metrosideros polymorpha
Correspondingly, taking as previous knowledge a manifestation of biphasic allometry as detected by the TSK-PLA scheme, for the Mascaro et al. (2011), we examined the possibility of the TSK-MPCA arrangement identifying a similar pattern in direct scales. Indeed, by setting r a = 0.855 the main_fun_tsk_mpca_model_fit.m function returned q = 2 for a biphasic mode and a w TSK x,p of reliable reproducibility. Identified firing strengths, display in Fig. 13A. Again analysis in direct scale detected by the TSK-MPCA approach spread about fitted mean response function E aTSK (y|x) compared to E gTSK (y|x) derived from retransformation of the TSK-PLA output. We can be aware that reproducibility strengths are equivalent. This can be stressed by performance metrics in Table 8. (C) through (D) presenting residual and QQ-plots confirm equivalence of E aTSK (y|x) and E gTSK (y|x).
Full-size DOI: 10.7717/peerj.8173/ fig-12 corroborates the consistency of break point allometry assumption for this data. We can learn of a break point estimated at a b = 8.8662. This estimate is consistent with resulting from a two linear segment mixture regression model performed by present authors. Spreads about fitted mean functions shown in Fig. 13B reveal remarkable correspondence of projections by E aTSK (y|x) and E gTSK (y|x). This can be stressed by performance metrics in Table 9. In turn this demonstrates adequacy of a TSK-PLA approach for the analysis of this data. Figures 13C and 13D   for this data led to considering an alternate clustering radius set at a value r a = 0.52. Resulting heterogeneity index was q = 3 that resulted in good model assessment metrics (Table 9) and a break point placed at x b = 4.832. This is in agreement with retransformed TSK-PLA estimation for this data. Nevertheless, forcing an interpolation mode of the TSK-MPCA to achieve a break point placed in agreement with previous estimation brings about complexity that renders biological interpretation difficult.

Identification of the TSK-MPCA proxy: Uca pugnax
Firing strengths, of a r a = 0.668, q = 2, TSK-MPCA fit to Huxley's Uca pugnax data set are displayed in Fig. 14A. We can be aware of heterogeneity as corresponds to dominance of sub models before and after the break point placed at x b = 340.7, matching exp(u b ) with u b = 5.83 the break point determined by a TSK-PLA fit to this data. Figure 14B show spread about resulting mean function E aTSK (y|x) and compares to E gTSK (y|x) gotten by retransformation of fitted TSK-PLA. Plot suggest equivalent reproducibility strengths. Nevertheless as it can be made certain by model performance metrics in Table 10 the E gTSK (y|x) proxy entails relatively higher reliability. Figures 14C and 14D one to one show residual and QQ plots corresponding to the TSK-MPCA fit. Similarly, residual and QQ plots for the back-transformed TSK-PLA fit appear in Figs. 14E and 14F. Table 10 also includes performance metrics for E gTSK (y|x) gotten by retransforming the output of the r a = 0.2 and q = 7, fit of the TSK-PLA. We can assert that resulting interpolation E gTSK (y|x) yields a relatively better fit. Therefore, results of the retransformed form of a TSK-PLA approach entails consistent results in direct scales. In other words, logtransformation based procedures do not lead to biased results for this data. But above all, results of a TSK-MPCA fit could provide a clue clearing an apparent misinterpretation of Huxley about existence of a break point in his analysis of Uca pugnax data. Indeed, as stated above we have u b = ln(a b ). This implies u b being the image of a b under logtransformation. Then, claiming existence of u b attributable to distortion set by a logtransformation itself is inappropriate. Agreeing with Packard (2012a), we have no doubt that conventional statistical methods do not put up with existence of u b as detected by the present fuzzy inference system. But, this fact cannot be exhibited to question fuzzy methods. These relying in non-probabilistic approaches have provided reliable interpretation of uncertainty as it can be inferred by fuzzy approach solutions to many problems of identification and control of nonlinear systems.

Identification of the TSK-MPCA proxy: Gadus chalcogrammu
When we assessed the performance of the TSK-MPCA device on the De Robertis & Williams (2008) data we found results consistent to the TSK-PLA fit reported above. Indeed, a TSK-MPCA analysis based on r a = 0.50 for this data returned q = 1, for a single membership  Table 10, the E gTSK (y|x) proxy entails relatively higher reliability. (C) and (D) display residual and QQ plots one to one for E aTSK (y|x). (E) and (F) show corresponding plots for E gTSK (y|x). We may be aware that log-transformation based procedures do not lead to biased results for this data. Full-size DOI: 10.7717/peerj.8173/ fig-14 E aTSK (y|x), as well as, residual and normal QQ plots in that order. Figure 15A also allows visual appraisal of a better reproducibility by E gTSK (y|x). Figure 15D presents spread about E as (y|x) and compares to E gTSK (y|x). We can observe that both proxies entitle similar reproducibilities. Figures 15E and 15F present residual spread and QQ plot accompanying E as (y|x). Table 11 compares reproducibility metrics for the E aTSK (y|x), E gTSK (y|x) and E as (y|x) proxies for this data. Again confrontation of model performance metrics shows that retransformation of the TSK-PLA output stands reliable results in direct scales.

DISCUSSION
A logarithmic transformation in allometry is often vindicated as a natural way to lodge a variation pattern resulting from multiplicative growth in plants and animals. Indeed, Gingerich (2000) and Kerkhoff & Enquist (2009) state that a number of biological processes, (i.e., growth, reproduction, metabolism and perception), are essentially multiplicative and are therefore prone to fit in better to a geometric error model. Beyond biological arguments supporting the traditional approach, Kerkhoff & Enquist (2009) underline that fitting models to log-transformed data is seamlessly adequate, since taking into account proportional rather than absolute variation is more significant. Therefore, from this standpoint, the fact that log-transformation places numbers into a geometric domain could bestow advantages beyond a purely statistical convenience. Nevertheless, there are remarks that a logtransformation approach procedure produces biased results, and that direct nonlinear regression methods in arithmetical scale, should be preferred in parameter identification tasks (e.g., Packard, 2013;Packard, 2009;Packard & Birchard, 2008;Packard & Boardman, 2008). But, these views are debatable for a school of defenders of the traditional protocol. For instance, Mascaro et al. (2014), stress on an important drawback in findings in Packard (2013) that refuted the traditional analysis method of allometry. This concerns the apparent significant bias linked to small values of the explanatory variable, that result from a use of nonlinear regression with the assumption of homoscedastic errors. Besides, Mascaro et al. (2014), underline that a lack of a CF misled Packard (2013), thereby explaining his assertion of biased results attributed to the logarithmic transformation protocol. Other practitioners have also placed a vigorous defense of this procedure, (e.g., Ballantyne, 2013;Glazier, 2013;Lai et al., 2013;White et al., 2012;Xiao et al., 2011). This is reasonably understood since inferences of many allometric studies could be invalidated by a substantiated rebuttal of this analysis method. But, Packard (2017a) asserts for instance, that adherence to a TAMA approach has been maintained even in situations when the resulting bivariate distribution is curvilinear in geometrical scale. Consequent pattern is generally referred as non-log linear allometry (Packard, 2012b;Strauss & Huxley, 1993;Echavarría-Heras et al., 2019a). Moreover, G.C. Packard has considered deviations from linearity in log-log plots of allometry as mainly attributable to a logtransformation itself (Packard & Boardman, 2008;Packard, 2012a;Packard, 2012b;Packard, 2013). From this perspective, overcoming the bias due to curvature in log scale requires extending complexity of Huxley's model of simple allometry in direct scales, which bears a paradigm of multiple parameter complex allometry (Gould, 1966;Lovett & Felder, 1989;MacLeod, 2014;Bervian, Fontoura & Haimovici, 2006;Packard, 2012a). Again, for promoters of the traditional approach this viewpoint sacrifices appreciation of biological theory in order to privilege statistical correctness (Houle et al., 2011;Lemaître et al., 2015;Pélabon et al., 2018). The approach underwent here demonstrates that a merging of points above can be achieved by evoking Huxley's report on the existence of a breakpoint in the log-log plot of chela mass vs. body mass of fiddler crabs (Uca pugnax) (Huxley, 1924;Huxley, 1927;Huxley, 1932). A generalization of this perspective explains adoption of a polyphasic loglinear allometry paradigm (Packard, 2016;Gerber, Eble & Neige, 2008;Strauss & Huxley, 1993;Hartnoll, 1978). This notion bestows curvature in geometrical space as determined by breakpoints interpreted as thresholds for transition among successive growth phases. Formally, this conception adds complexity for improving statistical consistency while keeping the meanings of allometric exponents as Huxley's original formulation. Conventional approaches have handled curvature in geometrical space by means of polynomial regression (Kolokotrones et al., 2010;Lemaître et al., 2014;MacLeod, 2010;Glazier, Powell & Deptola, 2013;Tidière et al., 2017;Echavarría-Heras et al., 2019a). Nevertheless, complexity underneath precludes accounting for heterogeneity as determined by break-point allometry. Conventional identification procedures also offer refined broken-line regression protocols (Beckman & Cook, 1979;Ertel & Fowlkes, 1976;Muggeo, 2003;Tsuboi et al., 2018;Ramírez-Ramírez et al., 2019;Echavarría-Heras et al., 2019a). Nevertheless, this slant relies on nonlinear regression that requires starting values for break-point estimation. Besides, crucial profile log likelihood could be log-concave so local maxima problems may exist. Surpassing this inconvenience may depend on using smoothed scatter plots to get candidate break points and consider several additional trials for estimation sensitivity to different starting points. Also necessary inferences on estimates could make implementation difficult (Julious, 2001;Muggeo, 2003).
The approaches in Bitar, Campos & Freitas (2016), Echavarría-Heras et al. (2018a) and in Echavarria-Heras et al. (2019b) typify fuzzy logic based hybrid paradigms aimed to allometric examination. Present TSK constructs can be placed in this framework. Moreover, as our results demonstrate offered fuzzy paradigm can naturally host complexity as intended in a break point assimilation of allometry. Moreover, conceived TSK arrangements offer direct-intuitive and starting value free identification of breakpoints. Certainly, beak points as conceived here correspond to points of intersection of TSK-firing strength factors. Besides, intervals in between break points can be interpreted as dominance realms of corresponding sub models. The TSK break point identification in geometrical space for the Echavarría-Heras et al. (2019a) andthe Mascaro et al. (2011) was paralleled by conventional broken-line regression. This confirms consistent capabilities by the fuzzy paradigm to identify heterogeneity in of the logtransformmed response. Thus, the offered TSK fuzzy model can be considered a tool entailing efficient automatic detection of weighted polyphasic log linear allometry patterns. And, the fact that the TSK model identified linearity in geometrical space for the De Robertis & Williams (2008) data demonstrates this approach can adapt complexity as necessary in an efficient way. But, we must emphasize that consistency of results in arithmetical space hinges on suitability of CF form. We suggest contemplating the optimal reproducibility criterion around Eq. (10) for this matter.
Motivation for the present research mainly stirred from the idea that identification based on a logarithmic transformation is suited for allometric examination. Visual inspection of TSK proxies fitted in geometrical space, as well as, included model performance metrics provides partial validation of our assertion. But, from the perspective of MPCA proponents, validation of detected heterogeneity should be made on the original arithmetic scales. Moreover, the addressed TSK-MPCA proxy corresponds to an expression of the general output of the TSK fuzzy model involving linear consequents in arithmetical scales. This arrangement is consistent with a MPCA approach as conceived in Lovett & Felder (1989). Furthermore, identification of a TSK-MPCA arrangement allows examination of break point allometry in arithmetical scales. Existence of break points in direct scales of data, confirms that a corresponding structure detected in geometrical space was not induced by effects of a logtransformation itself. And, using the Weierstrass approximation theorem, it can be demonstrated that the general output of a TSK fuzzy model can uniformly approximate any continuous function to arbitrarily high precision (Ying, 1998;Zeng, Nai-Yao & Wen-Li, 2000). Therefore, the high order interpolation capabilities of the TSK-MPCA scheme sets criterion to evaluate performance of retransformed TSK-PLA output E gTSK y|x . Certainly, as our results demonstrate this can be achieved by comparing the reproducibility strength of E gTSK y|x against that of E aTSK (w|a) for a given data set. And in the present settings the offered TSK-PLA or TSK-MPCA approaches were equally suited. This demonstrates that it is possible to maintain a logtransformation as part of a consistent allometric examination arrangement. This is a controversial subject whose clarification seems to be overcome by adopting presently offered analytical arrangement.
Packard (2012a) applied conventional statistical methods to conclude that a break point in Huxley's Uca pugnax data (Huxley, 1932) was inexistent. Nevertheless, application of present fuzzy methods detected a break point shifted back from the locus Huxley conceived. Corroboration of existence of this point seems to endure a biologically meaningful interpretation by Huxley of existence of a threshold for a sudden change in relative growth of the chela at about the time crabs reach sexual maturity. Likewise, detected break point in Zostera marina could be interpreted as a threshold beyond which plant assigns to leaves a relatively greater amount of tissue to resist damage and separation from shoots induced by drag forces. This implies different scaling parameters among small and large leaves (Echavarría-Heras et al., 2019a;Echavarría-Heras et al., 2018b;Echavarria-Heras et al., 2019b). Similarly, a detected break point in Metrosideros polymorpha may suggest different allometric scaling depending on tree size. Certainly, resource allocation to different tree traits like diameter or height could vary through growth in response to different environmental-biotic settings, and also to changes in resource availability. In this perspective, a risk of suppression by competitors may drive small trees to assign more resources to increase height (Echavarria-Heras et al., 2019b). Then, past a threshold height at which suppression risk is at a minimum resource may be apportioned to horizontal growth parameters such as diameter, crown and root cover (Weiner, 2004;Ramírez-Ramírez et al., 2019;Echavarria-Heras et al., 2019b). Therefore since the aim of allometric examination is understanding the biological processes that bring about covariance among traits, analytical approaches entailing break point identification must be preferred over conventional complex multi-parameter approaches (Echavarria-Heras et al., 2019b). Indeed, on spite of any gains in statistical fit attributable to the latter, characterization of inherent heterogeneity by the former could enhance biological insight. Particularly, a TSK-PLA slant could be a highly biologically significant model of allometry, because it can model the breakpoints while keeping the meanings of allometric exponents as Huxley's original formulation (Echavarria-Heras et al., 2019b).
As it is demonstrated by the steps in the derivation of Eq. (29), an imbedding of the TSK-PLA in the original theoretical perspective of allometry makes MPCA in arithmetical scale its logical consequence. By the same token the TSK-PLA approach grants direct-intuitive and starting values-free estimation of break points for transition among growth phases. We can also refer to benefits derived from the outstanding high order interpolation capabilities by this device. This functional mode of the TSK paradigm can be achieved by adjusting the value of the clustering parameter r a in Eq. (B7) (radii in supplied code) as to let the identification algorithm increase the number q of interpolation sub models in Eq. (20) (Table  9). We can be aware for instance that for membership functions in the form set by Eq. (51) consistent break point transference among geometrical and arithmetical scales is only achieved when r a = 0.52 which implies heterogeneity set by q = 3. Nevertheless, this by the way leads to a penalization in reproducibility strength relative to a fit by r a = 0.855. Setting a compromise between both fits depends on integration of previous knowledge into the analysis. This could help for instance by suggesting ad hoc forms of membership function with the aim of achieving high reproducibility and consistent break point placement relative to that previously estimated on geometrical space. In any event present digression on integration of subjective knowledge in the analysis of Huxley's data illustrates the extent on which a fuzzy logic approach can elucidate issues in allometric examination.

CONCLUSIONS
The offered TSK-PLA as formalized by the v TSK (u,π) paradigm can be interpreted as a generalized tool for the analysis of log transformed allometric data, that allows to contemplate: (1) the regression arrangement of the TAMA way (the case q = 1 and ϑ 1 (u) = 1), (2) a generalized nonlinear model for identification of weighted polyphasic nonlinear allometry (the case q > 1). (3) A direct-intuitive identification of concomitant break points for transition among successive growth phases.
On spite of a seemingly complicated formal set up of the v TSK (u,π) scheme, this can be conveniently identified by loading logtransformmed data into the provided code. Analysis of model performance metrics show that the mean response function E gTSK y|x deriving from retransformation of v TSK (u,π) to arithmetical space produces similar reproducibility strength as its counterpart E aTSK y|x following from identification of its arithmetical space TSK-MPCA counterpart w TSK x,p . Available conventional like broken line or weighted linear segment mixture regression approaches could offer reasonable analytical paradigms. Nevertheless, the offered TSK approach bears a flexible computational assembly for previous knowledge integration in an intuitive-interactive way. The present digression on Huxley's break point illustrates this advantage in a more proper way.
Present results confirm the pertinence of the quotation of Kerkhoff & Enquist (2009), on the uselessness of a distinction between logarithmic transformations and nonlinearity in many instances of allometric examination. Moreover, in our view, whenever we can manage to exhibit a suitable CF form proposed Takagi Sugeno Kang generalization can elucidate a glowing controversy. Surely, this paradigm allows the coexistence of the log transformation step claimed by practitioners as a must in allometry, and the unbiasedness of parameter estimates attributed to alternate direct nonlinear regression approaches in the original scale defended by many others.
However, the fact that TSK-PLA modeling provided meaningful interpretation in present settings does not rear this paradigm as a general tool of allometric examination. In the elucidating around Eq. (1) we established a condition on the response being positive and having a zero limit as covariate approaches zero. Therefore, the TSK-PLA slant essentially aims to analyse zero intercept allometries. And, there are instances where the initial timing of development of the trait itself and overall size are different. This situation will lead to consideration of a negative intercept in direct scales, ruling out transference of the examination into geometrical space. Then, modeling should be necessarily kept in direct scales and relying in MPCA turns out to be biologically reasonable. There are also situations where the error structure can be additive while the biological process underlying allometry is multiplicative. Again, this requests keeping the analyses on the arithmetic scales or modeling heteroscedastic errors in geometrical space. Certainly, we briefly addressed this approach while analyzing the eelgrass data. However, offering a heteroscedastic TSK-PLA protocol suited for the general settings requires further exploration. output of a Generalized Fuzzy Model (GFM) are mathematically equivalent. From this it follows that the probability density function of a Gaussian mixture can be reduced to accommodate the TSK fuzzy model. Likewise, by using the Weierstrass approximation theorem, Ying (1998) demonstrated that the general output of a TSK fuzzy can uniformly approximate any continuous function to arbitrarily high precision. In what follows we describe steps in the derivation of a w TSK x,p proxy for w x,p when the former is expressed as the general output of a first order TSK fuzzy model. In what follows, we review the steps in order to conceive a TSK surrogate w TSK x,p for w (x,p). For that aim, we recall our definition of w(x,p) as a continuous function w x,p : X → Y with both domain and range as sets of real numbers and p a set of parameters, and such that the allometric response y admits the representation y = w x,p .
In what follows, we use the symbol L x to denote the collection of membership functions describing the input variable x, that is, L x = µ k (x)|k = 1,2,....,q .
(A3) L x is known as a fuzzy partition of the input variable x in the domain X (Mendel, 2001;Bodjanova, 1993;Bezdek, 1981).
In the same way, we designate a collection T y of linguistic terms j associating to the output variable y. Namely T y = j |j = 1,2,...,r .
This way, the L y sets a fuzzy partition for the output variable y in the domain Y (Mendel, 2001;Bodjanova, 1993;Bezdek, 1981). Additionally, for i = 1,2,...,q, we advance correspondences i → A k(i) and i → B j(i) , therefore, we can contemplate antecedents P i (x) of the form and consequents Q i y backing inferential rules R i R i : We may now think about a single input-single output fuzzy inference system (Mendel, 2001). This is conceived as an application F : X → Y incorporating (1) a fuzzification module that characterizes the fuzzy partitions L x and L y , (2) an inference engine that uses the rules R = p 1 R i to convert a fuzzy input into a fuzzy output, and (3) a defuzzification operator D that transforms the fuzzy set obtained by the inference engine into a crisp value y in Y .
A first order single input-single output Takagi-Sugeno-Kang fuzzy inference system (Sugeno & Kang, 1988;Takagi & Sugeno, 1985) considers decision rules R i having an antecedent P i (x) of the form given by Eq. (A7) but with the consequent Q i y in expression (A8) taking a crisp functional form f i (x). That is, in the TSK fuzzy inference system we consider R i rules of the form for i = 1,2,...,q. We notice also that being the consequent a real number the use of a defuzzification operator is not necessary. An important component in a TSK fuzzy model is the firing strength ϕ i (x) of the antecedent P i (x) of a rule R i . For a first order single input-single output TSK fuzzy model we take A normalized firing strength ϑ i (x) takes a form (Mendel, 2001) It follows that The final output w TSK x,p of the Takagi-Sugeno-Kang inference system is the normalized firing strength weighted average of all rule outputs (Sugeno & Kang, 1988;Takagi & Sugeno, 1985), that is, Where p stands for the set of parameters identifying the membership and consequent functions in Eqs. (A3) and (A10) one to one.

APPENDIX B. IDENTIFICATION PROCEDURES FOR W TSK (X ,P)
Description of structure and parameter estimation of the TSK fuzzy model interrelate (Echavarría-Heras et al., 2018a). A first stage relies on Subtractive Clustering (Castro et al., 2016;Chiu, 1994). This technique sets decision rules R i and produces parameter estimates for the µ k (x) membership functions. This acquires as well the forms of the normalized firing strength factors ϑ i (x). A second stage of the identification task is achieved by placing weight factors ϑ i (x) in Eq. (A14) in order to obtain parameter estimates for the consequents of the rules f i (x). Regularly, this is achieved by means of recursive least squares methods (Jang, Sun & Mizutani, 1997;Wang & Mendel, 1992).

The subtractive clustering method
The subtractive clustering method (SC) is an extension of the Mountain Function Clustering method (MFC) proposed originally by Yager & Filev (1994). This procedure estimates cluster centers based on the notion of a density function. For clarifying aims, before we explain the SC procedure, it is convenient to describe the MFC method. Following Yager & Filev (1994) and Chiu (1994), the MFC procedure assembles the following steps: 1. The set of n data to be analyzed is arranged as a vector X namely: X = {x 1 ,x 2 ,...,x n }.
(B1) 2. Generation a nth dimensional space grid on which the data is located. Intersections of grid lines provide nodes N i . i = 1,...,m. Cluster centers to be are restricted to grid nodes. The set of cluster centers to be denotes through C.
Based on the distribution of the data devise a mountain function (MF). This represents a data density measure. The height of the mountain function at a node point N i is where α is a positive constant and d (x k ,N i ) is a measure of the distance between data point x k and grid node N i . Equation (B2) enunciates that the data density measure at a point N i is influenced by all data points x k . Such a density measure varies inversely proportional to distance between data points and node N i . The parameter α not only influences the maximum value of MF (N i ) but also its smoothness. It can be ascertained that the values of the mountain function on nodes are closely dependent on density of data points in the neighborhood of N i . Mountain function values can be also interpreted as the potential suitability of a grid node to become a cluster center estimate. A node with many neighboring data points will have a large mountain function value. Cluster center C 1 is chosen as node N i such that MF (N i ) attains the maximum value among all remaining nodes, formally Since the nodes close to C 1 will also have high mountain-function values, it deems necessary to remove the effect of C 1 before obtaining the next cluster center C 2 . A new mountain function MF 2 (N i ) is shaped by taking off a scaled Gaussian function centered at C 1 this eliminates the effect of the first cluster. Iteratively the mountain function, after eliminating the effects of the cluster center that was previously identified becomes 3. The previous step is repeated until the number of desired cluster centers is found or until a stopping condition is met. This expresses in terms of the ratio of first maximum value of mountain function found MF 1 (C 1 ) to corresponding penultimate maximum value found MF k−1 (C k−1 ). Iteration stops when this ratio attains a value less than a certain positive constant δ, Since, the mountain function has to be evaluated at each grid point the processing time of the mountain clustering method rises exponentially with dimension of task. A SC scheme amends these difficulties by taking data points as potential cluster centers. This implies processing time becoming proportional to problem size instead of problem dimension. Since this method does not take into account any grid intersections execution time is reduced. Nevertheless, the real clusters centers do not necessarily place at data points, but in most cases this approach offers a good approximation, to cluster center identification. The SC approach also bases on a function representing the distribution of the data. Actually, the SC method consists of equations very similar to those used in the MF method, steps in the later are: 1. The set of n data to be analyzed is defined as follows: X = {x 1 ,x 2 ,...,x n } (B6) 2. Since, each data point is an aspirant for a cluster center a density measure at data point x i is defined as where r a is a positive constant, this constant acts as the radius that defines the area of proximity to the potential cluster center. Once all density assessments are obtained, the one with the highest value is taken. Formally ) this will identify a first cluster center C 1 , 3. Next cluster center is acquired by subtracting a scaled Gaussian function centered at C 1 , that is where r b < r a , is a radius that defines the proximity area of the cluster center C 1 . 4. The maximum value of rescaled densities D i of Eq. (B9) is chosen as cluster center C 2 .
This procedure is repeated until a desired number m of cluster centers C 1 ,C 2 ,....C m is determined.

The recursive least squares method
The general least-squares problem establishes the output of a model y as given by a linearly parameterized expression, namely y = γ 1 h 1 (x) + γ 2 h 2 (x) + ··· + γ n h n (x), where x = x 1 ,··· ,x p T is the model's input values vector, h 1 (x),··· ,h n (x) are known functions of x, and γ 1 ,··· ,γ n called regression coefficients are to be fitted. Without loss of generality, we address the case q = 2 assuming consequent linear functions in the form given by Eq. (53), so that the general output of the TKS of Eq. (A14) is w TSK x,p = p 1 1 ϑ 1 (x)x + p 2 1 ϑ 2 (x)x + p 1 2 ϑ 1 (x) + p 2 2 ϑ 2 (x).
with ϑ i (x) given by Eq. (A12). Then, w TSK x,p as given by Eq. (B11) becomes a particular characterization of Eq. (B10) by taking model's input values x = [x] T , and p 1 1 , p 1 2 , p 2 1 and p 2 2 unknown parameters in the consequent functions.
To obtain parameter estimates, we take into account that in the present settings the target system to be modeled involves an input-output relationship x → w TSK x,p being x the descriptor variable and w TSK x,p standing for the response y. Therefore, we have a training data composing pairs x k : y k , for k = 1,··· ,m that stand for replicates of the considered input-output relationship. Therefore, in order to identify the unknown parameters p 1 1 , p 1 2 , p 2 1 and p 2 2 , we must fill in for each data pair x k : y k , into Eq. (B11) in order to obtain the set of m linear equations: 1 ϑ 1 (x 1 )x 1 + p 2 1 ϑ 2 (x 1 )x 1 + p 1 2 ϑ 1 (x 1 ) + p 2 2 ϑ 2 (x 1 ) = y 1 ϑ 1 (x 2 )p 1 1 x 2 + p 2 1 ϑ 2 (x 2 )x 2 + p 1 2 ϑ 1 (x 2 ) + p 2 2 ϑ 2 (x 2 ) = y 2 . . . . . . . . .
This system of equations can be equivalently written in a concise form BP = y, where B is the m × n matrix, ϑ 1 (x 1 )x 1 ϑ 1 (x 1 ) ϑ 2 (x 1 )x 1 ϑ 2 (x 1 ) ϑ 1 (x 2 )x 2 ϑ 1 (x 2 ) ϑ 2 (x 2 )x 2 ϑ 2 (x 2 ) . . . and being y the m × 1 output values vector: The i-th row of the data matrix B . . .y is denoted by b T i ,y i and formally represented by, Then, Eq. (B12) modifies to include an error vector e that accounts for random noise or modeling error, that is, Since e = y − BP then e T e = y − BP T y − BP , and if we let E (P) = e T ewe will have We call E (P) the sum of squared errors. Then, we need to search for a characterization P of the vector P, which minimizes E (P). Furthermore, the vectorP is known as the least-squares estimator (LSE) of P. Since E (P) is in quadratic form,P is unique. It turns out thatP satisfies the normal equation Furthermore,P is given bŷ A n-order least squares estimatorP n ofP defined by means of the expression is a description ofP that associates to n data pairs taken out of the training data set x i : y i . Once we have gottenP n we can acquire the following estimatorP n+1 with a minimum of effort, through a recursive least-squares estimator (RLSE) technique, a procedure where the nth row of B . . .y , with (1 ≤ n ≤ m) denoted by b T n . . .y n is recursively obtained. We now explain the procedure behind the RLSE method.
A new pair b T n+1 ;y n+1 becomes available as the (n + 1) th entry in the data set, producing theP n+1 estimate, Further, in order to simplify the notation, the pair b T n+1 ;y n+1 will be symbolized by (b T ;y) and we also introduce the p × p matrices H n and H n+1 defined by means of and ThusP n+1 can be recursively identified in terms of the preceding estimateP n and the new data pairs b T ;y . Furthermore, the current estimateP n+1 is expressed as the previous oneP n plus a correcting term based on the new data b T ;y ; this adjusting term can be understood as an adaptation gain vector H n+1 multiplied by a prediction error (y − b TP n ) linked to the previous estimatorP n .
The descending gradient method is one of the oldest techniques to minimize, a given function defined in a multidimensional input space. This method bases many direct optimization methods for restricted and unrestricted problems. On spite of its slow convergence, its simplicity makes it the most used non-linear optimization technique. Formally ϕ next = ϕ now − ηG. (B32) A slightly different formulation of Eq. (B32) results from a gradient normalization, that is being κ the actual size of step ,which interprets as the Euclidean transition distance from ϕ now to ϕ next , namely In order to typify (Eqs. B32) and (B33) the former one refers as simple descendent gradient and the later as its normalized version. The term ηG in Eq. (B32) stands for the extent of step. With a fixed η, step magnitude changes automatically in each iteration due to different gradients of G. If the minimum point is on a flat surface or plateau, G tends to be infinitesimally small. Consequently, the simple descendent gradient in Eq. (B32) exhibits a slow convergence. On the other hand, for a fixed κ the normalized simple descending gradient in Eq. (B33) always does the same steps, neglecting how steep the slope is Chan & Fallside (1987) and Rumlhart, Hinton & Williams (1986). It is then necessary to actualize step size κ for efficiency. data and, therefore, more rules. The value of λ in a membership function is related to the radii value and the range of entries through > λ = (radii.*Range(X))/ sqrt(8.0) (B37) (Chiu, 1994). The κ in the descendent gradient method adapts the θ and λ values of the membership functions (cf. Eqs. (21) and (52)) in the antecedents of each iteration or epoch while minimizing the objective function.

APPENDIX C. MODEL PERFORMANCE METRICS
Besides AIC and ρ indices, model assessment in this examination relies on the SEE, MPE and MPSE indices based on statistics of squared and absolute deviations of observed to predicted values. According to Parresol (1999), SEE, MPE and MPSE statistics as model performance metrics were first recommended by Meyer (1938), then by Schlaegen (1982) and have subsequently been used by Zeng & Tang (2011a); Zeng & Tang (2011b). We provide ahead related formulae and explanation.

Akaike information criterion (AIC)
Lin's Concordance Correlation Coefficient (ρ C ) with ρ standing for Pearson's correlation coefficient. The ρ C index estimates througĥ