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

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).

Models

Huxley’s formula of Simple Allometry

A characterization of $w\left(x,\mathit{p}\right)$ 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 (11)${\displaystyle y=w_s\left(x,\mathit{p}\right)+ϵ}$ with $w_s\left(x,\mathit{p}\right)=\beta x^{\alpha }$ and ϵ assumed to be normally distributed with zero mean and variance σ², that is, $ϵ\sim \psi \left(0,{\sigma }^2\right)$. According to our notation convention Eq. (11) yields the mean response function $E_{as}\left(y|x\right)=\beta x^{\alpha }$.

Similarly, the logtransformation method produces the TAMA’s regression model, that is, (12)${\displaystyle v=v_s\left(u,\pi \right)+ϵ}$ with (13)${\displaystyle v_s\left(u,\pi \right)=ln\beta +\alpha u}$ and ϵ as specified above. Equations (12) and (13) determine $E_s\left(v|u\right)=v_s\left(u,\pi \right)$. Accordingly, back transformation of Eq. (12) to arithmetical space brings up a mean response $E_{gs}\left(y|x\right)$ given by (14)${\displaystyle E_{gs}\left(y|x\right)=\beta x^{\alpha }\delta ,}$ where δ stands for necessary CF.

It often occurs that even after contemplation of proper form for δ this TAMA’s derived proxy for $w\left(x,\mathit{p}\right)$ produces biased projections of observed values of the response y. This means, that complexity of Huxley’s formula of simple allometry $w_s\left(x,\mathit{p}\right)$ becomes inappropriate to identify the true $w\left(x,\mathit{p}\right)$ 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\left(x,\mathit{p}\right)$ could depend on MPCA forms. Corresponding logtransformmed expression $v\left(u,\pi \right)$ 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.

Data

Plots in 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\left(a,\mathit{p}\right)=\beta a^{\alpha }$ that includes scaling parameters α and β depending in the covariate, that is, (23)${\displaystyle y=\beta \left(x\right)x^{\alpha \left(x\right)}}$ where $\beta \left(x\right)$ and α(x) are continuous functions and with $\beta \left(x\right)$ assumed to be positive. This sets (24)${\displaystyle w\left(x,\mathit{p}\right)=\beta \left(x\right)x^{\alpha \left(x\right)}.}$ Thus, formally whenever the scaling functions $\beta \left(x\right)$ and α(x) are not simultaneously constant $w\left(x,\mathit{p}\right)$ 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 (25)${\displaystyle v_{TSK}\left(u,\pi \right)=ln{\beta }_{TSK}\left(u\right)+{\alpha }_{TSK}\left(u\right)u}$ where (26)${\displaystyle {\beta }_{TSK}\left(u\right)=e^{\sum _1^q{\vartheta }^i\left(u\right)ln{\beta }_i}}$ and (27)${\displaystyle {\alpha }_{TSK}\left(u\right)=\sum _1^q{\vartheta }^i\left(u\right){\alpha }_i.}$ Thus, Eq. (17) takes on the equivalent representation, (28)${\displaystyle v=ln{\beta }_{TSK}\left(u\right)+{\alpha }_{TSK}\left(u\right)u+{ϵ}_{TSK}.}$

The functions $ln{\beta }_{TSK}\left(u\right)$ and ${\alpha }_{TSK}\left(u\right)$ 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 (29)${\displaystyle E_{gTSK}\left(y|x\right)={\beta }_{TSK}\left(x\right)x^{{\alpha }_{TSK}\left(x\right)}\delta .}$

This sets $exp\left(v_{TSK}\left(u,\pi \right)\right)={\beta }_{TSK}\left(x\right)x^{{\alpha }_{TSK}\left(x\right)}.$ But, from Eq. (5) we have $v_{TSK}\left(x,\mathit{p}\right)=ln\left(w_{TSK}\left(x,\mathit{p}\right)\right)$ then $w\left(x,\mathit{p}\right)$ as identified by retransformation of $v_{TSK}\left(x,\mathit{p}\right)$ admits the form specified by Eq. (24).

Fitting results of the TAMA protocol: Zostera marina

For comparison aims, we present fitting results of the TAMA on the Echavarría-Heras et al. (2019a) data. Estimates for the allometric parameters α and β derive from linear regression on log-transformed data (v, u) (cf. Eq. (12)). Table 1 summarizes the results of the analysis. Corresponding, mean response $E_s\left(v|u\right)$ 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\left(v|u\right)$. Spread may lead on first sight to the impression that the distribution of v around the mean response line $E_s\left(v|u\right)$ 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., (30)${\displaystyle ϵ{\left(u\right)}^2=b+du+\zeta }$ 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 ${\vartheta }^{\mathit{i}}\left(\mathit{u}\right)$

In order to identify the required firing strength factors ${\vartheta }^{\mathit{i}}\left(\mathit{u}\right)$ 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 ${\mu }_{{\Phi }_k}\left(u\right)$ 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 ${\vartheta }^1\left(u\right)$ and ${\vartheta }^2\left(u\right)$ turn out to be (31)${\displaystyle {\vartheta }^1\left(u\right)=\frac{1}{1+exp\left\{-\frac{1}{2}\left[{\left(\frac{u-{\theta }_2}{{\lambda }_2}\right)}^2-{\left(\frac{u-{\theta }_1}{{\lambda }_1}\right)}^2\right]\right\}}}$ (32)${\displaystyle {\vartheta }^2\left(u\right)=\frac{exp\left\{-\frac{1}{2}\left[{\left(\frac{u-{\theta }_2}{{\lambda }_2}\right)}^2-{\left(\frac{u-{\theta }_1}{{\lambda }_1}\right)}^2\right]\right\}}{1+exp\left\{-\frac{1}{2}\left[{\left(\frac{u-{\theta }_2}{{\lambda }_2}\right)}^2-{\left(\frac{u-{\theta }_1}{{\lambda }_1}\right)}^2\right]\right\}}.}$

Plots of the estimated membership functions μ_Φ1(u) and μ_Φ2(u) and normalized firing strength factors ${\vartheta }^1\left(u\right)$ and ${\vartheta }^2\left(u\right)$ 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 ${\mathit{f}}^{\mathbf{i}}\left(\mathit{u}\right)$

A second step in the procedure to get v_TSK(u, π) concerns acquiring the consequent functions $f^{\mathrm{i}}\left(u\right)$ in Eq. (22). Since, for this data, we obtained q = 2, the code ought to establish consequent functions $f^1\left(u\right)$ and $f^2\left(u\right),$each one assumed to be linear, that is, (33)${\displaystyle f^1\left(u\right)=ln{\beta }_1+{\alpha }_{1}u}$ and (34)${\displaystyle f^2\left(u\right)=ln{\beta }_2+{\alpha }_{2}u.}$

With the aim of assessing heteroscedasticity, we replaced the forms of ${\vartheta }^1\left(u\right)$ and ${\vartheta }^2\left(u\right)$ identified by SC technique in regression Eq. (17). In turn the involved consequent functions $f^1\left(u\right)$ and $f^2\left(u\right)$ 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 (35)${\displaystyle {\sigma }_{TSK}\left(u\right)=log\left(\sigma +ku\right),}$ 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 ϑ¹(u)f¹(u) and ϑ²(u)f²(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\left(u,\pi \right)$, 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 ϑⁱ(u)fⁱ(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 ϑ¹(u)f¹(u). Conversely, the larger the distance between u and u_bfor leaves in the large phase u ≥ u_b the closer v_TSK(u, π) will be to ϑ²(u)f²(u). Relating to E_s(v|u) shown in Fig. 2A, we can assess from Table 4 and Fig. 4A that the reproducibility strength of E_TSK(v|u) is higher. Besides compared with Fig. 2B, the plot in Fig. 4B shows that distribution of residuals about the zero line for the TSK fit improved. Also compared to Fig. 2C, normal QQ plot in 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⁻¹⁶, 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\left(u\right)$ and $f^2\left(u\right)$ as given by Eqs. (33) and (34) in Eq. (20) direct algebraic manipulation yields (36)${\displaystyle v_{TSK}\left(u,\pi \right)=ln\left({\beta }_1^{\left(1-{\vartheta }^2\left(u\right)\right)}{\beta }_2^{{\vartheta }^2\left(u\right)}\right)+\left({\alpha }_1\left(1-{\vartheta }^2\left(u\right)\right)+{\alpha }_2{\vartheta }^2\left(u\right)\right)u.}$

Similarly, it can be directly verified that firing strengths ${\vartheta }^1\left(u\right)$ and ${\vartheta }^2\left(u\right)$ as given by Eqs. (31) and (32) can be also expressed in the form (37)${\displaystyle {\vartheta }^1\left(u\right)=\frac{1}{1+e^{\tau \left(u,\theta ,\lambda \right)}}}$ and (38)${\displaystyle {\vartheta }^2\left(u\right)=\frac{e^{\tau \left(u,\theta ,\lambda \right)}}{1+e^{\tau \left(u,\theta ,\lambda \right)}},}$ where (39)${\displaystyle \tau \left(u,\theta ,\lambda \right)=\psi \left(\lambda \right)\varphi \left(u,\theta ,\lambda \right)+\xi \left(\theta ,\lambda \right),}$ (40)${\displaystyle \psi \left(\lambda \right)=\frac{\left({\lambda }_2^2-{\lambda }_1^2\right)}{2{\left({\lambda }_1{\lambda }_2\right)}^2},}$ (41)${\displaystyle \varphi \left(u,\theta ,\lambda \right)={\left[u^{}+\left(\frac{{\lambda }_2^2{\theta }_1-{\lambda }_1^2{\theta }_2}{{\lambda }_1^2-{\lambda }_2^2}\right)\right]}^2}$ (42)${\displaystyle \xi \left(\theta ,\lambda \right)=\frac{{\left({\lambda }_2^2{\theta }_1-{\lambda }_1^2{\theta }_2\right)}^2-\left({\lambda }_1^2-{\lambda }_2^2\right)\left[{\left({\lambda }_1{\theta }_2\right)}^2-{\left({\lambda }_2{\theta }_1\right)}^2\right]}{2\left({\lambda }_1^2-{\lambda }_2^2\right){\left({\lambda }_1{\lambda }_2\right)}^2}.}$ with θ, λ standing for parameter vectors $\left({\theta }_1,{\theta }_2\right)$ and (λ₁, λ₂) one to one. We can then ascertain that $\varphi \left(u,\theta ,\lambda \right)$ remains positive for all values of u. Also, the term $\xi \left(\theta ,\lambda \right)$,does not depend on u. Consequently, whenever the factor $\psi \left(\lambda \right)$ in Eq. (40) is positive, $\tau \left(u,\theta ,\lambda \right)$ will approach infinity as u approaches infinity. Then, Eq. (37) implies ${\vartheta }^1\left(u\right)$ asymptotically vanishing as u approaches infinity. Reversely, whenever $\psi \left(\lambda \right)$is negative, the firing strength factor ${\vartheta }^1\left(u\right)$ will asymptotically approach one as u approaches infinity. For the Echavarría-Heras et al. (2019a) data set we obtained $\psi \left(\lambda \right)=0.2756$, then we must have (43)${\displaystyle \underset{u\to \mathrm{\infty }}{lim}{\vartheta }^1\left(u\right)=0}$ and since Eq. (A13) implies ${\vartheta }^2\left(u\right)=1-{\vartheta }^1\left(u\right)$, we also have (44)${\displaystyle \underset{u\to \infty }{lim}{\vartheta }^2\left(u\right)=1.}$

Agreeing with Eqs. (19) and (20) back-transformation e^v produces (45)${\displaystyle E_{gTSK}\left(y|x\right)={\beta }_1^{\left(1-{\vartheta }^2\left(u\left(x\right)\right)\right)}{\beta }_2^{{\vartheta }^2\left(u\left(x\right)\right)}{a}^{\left({\alpha }_1\left(1-{\vartheta }^2\left(u\left(x\right)\right)\right)+{\alpha }_2{\vartheta }^2\left(u\left(x\right)\right)\right)}\delta .}$

We denote by means of the symbol, $E_{gTSK}^{\infty }\left(y|x\right)$ the limit of $E_{gTSK}\left(y|x\right)$ as x approaches infinity, that is, (46)${\displaystyle E_{gTSK}^{\infty }\left(y|x\right)\underset{x\to \infty }{=lim}{E}_{gTSK}\left(y|x\right).}$ Then, Eqs. (31), (32) and (44) through (46) imply (47)${\displaystyle E_{\text{gTSK}}^{\infty }\left(y|x\right)={\beta }_2{x}^{{\alpha }_2}\delta .}$

Then, the asymptotic mode $E_{gTSK}^{\infty }\left(y|x\right)$ identified for the Echavarría-Heras et al. (2019a) data set, is attains a form like Huxley’s formula of simple allometry $w_s\left(x,p\right)$. Estimated parameters are α₂ = 1.1126 and β₂ = 7.8398 × 10⁻⁶. Figure 5C displays observed leaf biomass values y and their projections through the $E_{gTSK}^{\infty }\left(y|x\right)$ proxy. We can learn of a remarkable correspondence between the power function $w_s\left(x,p\right)=\beta x^{\alpha }$ of Eq. (11) fitted by direct nonlinear regression methods and the asymptotic mean response $E_{gTSK}^{\infty }\left(y|x\right).$ Besides as established by Eq. (45) we can directly asses from Fig. 5C that for sufficiently large values of x in the Echavarría-Heras et al. (2019a) data set, the mean response E_gTSK(y|x) behaves as the power function $E_{gTSK}^{\infty }\left(y|x\right)$. Moreover, the order relationship u ≥ u_b holds for about 86% of analyzed data. This explains why corresponding phase of the TSK output can be considered dominant. This by the way elucidates the apparent benefit of fitting Huxley’s formula of simple allometry by means of nonlinear regression model directly in arithmetical scale for the considered data. Indeed, such a fitting could deliver reasonable model adequacy results. But, as the present results show direct nonlinear examination based on Huxley’s formula of simple allometry will fail to detect the different allometrical phases conforming the real variation pattern in the data. Then, as we have elaborated a log transformation step followed by nonlinear model identification in geometrical space could overcome the reproducibility deficiency of the TAMA approach.

Fitting results of the TSK-PLA assembly: Metrosideros polymorpha

For trying the main_fun_tsk_pla_model_fit.m function on the Mascaro et al. (2011) data we set r_a = 0.80. This returned q = 2, heterogeneity. Figure 6 displays the plots of membership functions μ_Φi(u), firing strength factors ${\vartheta }^i\left(u\right)$, consequent linear segments $f^i\left(u\right)$ and component products ${\vartheta }^i\left(u\right)f^i\left(u\right)$ identified by the fit of the TSK fuzzy model to the Mascaro et al. (2011) data set. Membership functions are shown in Fig. 6A. Fit suggests heterogeneity determined by a break point u_b = 1.575 as shown in Fig. 6B displaying firing strength factors. 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. Break point suggest a growth phase 0 < u ≤ u_b and a complementary u > u_b. We can interpreted these regions as dominance realms for the component product functions ${\vartheta }^1\left(u\right)f^1\left(u\right)$ and ${\vartheta }^2\left(u\right)f^2\left(u\right)$ one to one (Figs. 6C and 6D). Correspondingly, Fig. 7A shows spread about mean response function regions in geometrical space matching identified phases. Moreover, residual plot in Fig. 7B displays a 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.2676 and α₂ = 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 ${\vartheta }^1\left(u\right)f^1\left(u\right)$ and ${\vartheta }^2\left(u\right)f^2\left(u\right)$. 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}\left(u,\pi \right)$. Accordingly, setting r_a = 0.2 brought about q = 7 sub models, inducing a maximum reproducibility strength of interpolation function $v_{TSK}\left(u,\pi \right)$ 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}\left(u,\pi \right)$ suggests a pattern accommodating two linear segments that alternate about u_bI. Moreover, using the interpolation points $\left(u,v_{TSK}\left(u,\pi \right)\right)$ 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}\left(u,\pi \right)$, 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 ${\vartheta }^1\left(u\right)$ set to one (Fig. 10B), and a conforming single TAMA’s form linear consequent $f^1\left(u\right)$ (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\left(x,\mathit{p}\right)$ as intended for MPCA can be modeled by $w_{TSK}\left(x,\mathit{p}\right)$ as expressed by means of Eq. (A14), in arithmetical space, namely (48)${\displaystyle w_{TSK}\left(x,\mathit{p}\right)=\sum _1^q{\vartheta }^i\left(x\right)f^i\left(x\right)}$ with firing strengths ${\vartheta }^i\left(x\right)$ given by (49)${\displaystyle {\vartheta }^i\left(x\right)=\frac{{\mu }_{{\Phi }_i}\left(x\right)}{\sum _1^q{\mu }_{{\Phi }_k}\left(x\right)}}$ being μ_Φi(x) for i = 1, 2……..q the involved membership functions. We also undertake that both $w\left(x,\mathit{p}\right)$ and x remain positive, and that (50)${\displaystyle \underset{x\to 0^{+}}{lim}w\left(x,\mathit{p}\right)=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 (51)${\displaystyle {\mu }_{{\Phi }_i}\left(x\right)=c\left(e^{h_i\left(x\right)}-1\right)}$ where c = (1 − e)⁻¹ and (52)${\displaystyle h_i\left(a\right)=e^{\left\{-\frac{1}{2}\left[{\left(\frac{lna-{\theta }_i}{{\lambda }_i}\right)}^2\right]\right\}}}$ with θ_i and λ_i for i = 1, 2, …., q parameters. Correspondingly we consider the consequents $f^i\left(x\right)$ to be linear functions, that is, (53)${\displaystyle f^i\left(x\right)=c_{i1}+c_{i2}x.}$

It is worth recalling that Eq. (15) provides the form of linked regression protocol.

Identification of $w_{TSK}\left(x,\mathit{p}\right)$ as given by Eq. (48) through Eq. (53) from data pairs (y, x) in direct scale is performed by means of the Matlab function main_fun_tsk_mpca_model_fit.m supplied in the supplemental files section. Heterogeneity and reproducibility strength features of $w_{TSK}\left(x,\mathit{p}\right)$ can be explored in an interactive way through different characterizations of the clustering radius parameter r_a as specified by Eq. (B7) through Eq. (B9).

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}\left(x,\mathit{p}\right).$ 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}\left(y|x\right)$ appear in Fig. 12B. This plot allows comparison to its counterpart $E_{gTSK}\left(y|x\right)$ 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}\left(x,\mathit{p}\right)$ of reliable reproducibility. Identified firing strengths, display in Fig. 13A. Again analysis in direct scale detected by the TSK-MPCA approach 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 display residual and QQ plots for TSK-MPCA fit to Metrosideros polymorpha. Equivalent plots associating to the retransformed TSK- PLA are displayed in Figs. 13E and 13F, correspondingly. Exploring interpolation capabilities of the TSK-MPCA 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\left(a_b\right)$. 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 function. This consequently associates to a single firing strength ${\vartheta }^1\left(x\right)=1$. As a result, we have to contemplate a single component product of a linear form in Eq. (48). No break point composed heterogeneity in direct scales can be verified for this data. Moreover, implied linear form of Eq. (48) does not fit required complexity in direct scales. Nevertheless, previous knowledge on consistency of corresponding TSK-PLA fit suggest using the interpolation features of TSK-MPCA to grant adequacy. For this empirical aim, for instance taking the clustering parameter r_a = 0.22 in Eq. (B7) we can manage to obtain an heterogeneity index of q = 3. This entails three sub models composing Eq. (48). Figures 15A through 15C display spread about resulting form of interpolation mean response function 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.