Modeling Height—Diameter Relationship Using Artificial Neural Networks for Durango Pine Species in Mexico

1. Introduction

Artificial Intelligence (AI) refers to the simulation of human intelligence in machines, wherein they are programmed to think and learn in a manner similar to humans. The AI has been used in forest modelling for different purposes and objectives. The machine learning (ML) is a subset of AI that focuses on the development of algorithms and models that enable computer to learns from a specific dataset and make predictions or take actions without being explicitly programmed [1,2]. One of the most ML techniques is the artificial neural network (ANN) and the resilient backpropagation artificial neural network (RBPANN) perform supervised ML in multi-layer perceptron, and the main principle is to eliminate the harmful influence of the size of the partial derivative on the weight step [3,4,5]. The ANNs are computational models inspired in the natural neurons and they represent a generalization of mathematical models of human cognition or neural biology [1,6,7]. In ANNs, the training and testing datasets are used to train and evaluate the performance of the network for a specific randomly selected dataset. The training dataset is used to train the neural network. It consists of a set of input data points and their corresponding target output values, while, the testing dataset is a separate dataset that is used to evaluate the performance of the trained neural network [1,8,9].

One of the most important relationships in forest modeling is the total tree height and diameter at breast height (h-dbh) relationship, and this relationship is usually applied in forest inventory or for height estimation in forest management and planning [10]. The knowledge of h-dbh relationship is a fundamental both developing and applying many growth and yield models [11,12]. This relationship has mainly studied with nonlinear mixed effect modeling (NLMEM) with fix and random parameters for several species and grouping level or ecological conditions [10,11,13,14,15,16]. Lately, this relationship has been studied with AI, and the ML thought ANN has been used [7,17,18]. Also, other variables as basal area [19], crown width [20], biomass [21], volume [22], forest fire [23], and annual radial growth with competition indices [24] have been studied with different ML algorithms. Occasionally, the clustering analysis based in unsupervised ML has been included in to group similar data point together based on their inherent characteristics or similarities [1,25,26,27]. The unsupervised clustering analysis could identify patterns or structures in dataset to improve the fitted models in forest modeling.

Specifically, in the Mexican Forestry the h-dbh relationship has been extensively studied with NLMEM for local and generalized models and occasionally the unsupervised cluster analysis was included in modelling [12,28,29]. The NLMEM are better than fitted models by ordinary least squares method and those use random parameter to explain the variability between groups, sites, or ecological regions. Lately, the ML algorithms are taken attention in forestry research and the results outperform the NLMEM approach for the h-dbh relationship. In ANN analysis is convenient to separate the dataset in two parts, one for training and other one for testing or validation [7,17]. The main used model for NLMEM has been th Chapman-Richards model [30], which is base in a sigmoid relationship growth based on age [31].

Considering the above schemes and the context of AI in forestry research, this study aim the h-dbh relationship for Durango Pine species (Pinus durangensis Martínez) by NLMEM and ANN for unsupervised clustered dataset for training and testing sets. The algorithms were compared in both training and testing phases and some conventional statistics like root mean square error, coefficient of determination, Akaike’s information criterion, Bayesian information criterion, and loglikelihood were uses to perform the approaches. The resilience backpropagation of ANN (RBPANN) was employed, and three activation functions were computed and evaluated. The activation functions were tangent hyperbolicus (RBPANN-tanh), softplus (RBPANN-softplus), and logistic (RBPANN-logistic), and those were trained by resilience backpropagation and maximum likelihood was used.

2. Materials and Methods

2.1. Study area

The study was developed in a forest community in Norther Mexico, specifically in Durango state. The forest community is called San Diego de Tezains, and the total area is around 61,098 ha, which 30,000 ha are used for forest management and timber harvesting. The main applied silvicultural treatments are based on continuous cover forestry (CCF) and rotation forest management (RFM) [32]. The silvicultural treatments for CCF area base on selection, while three thinning and shelterwood cutting treatments for RFM with 15-years of forest cycle cutting [33]. The location of study area is showed in Figure 1. The mean annual temperature ranges from 5 to 18 °C, and the lowest temperature occurs in January (− 6 °C) and the hottest in May (28 °C). The altitude varies from 1,500 to 3,033 m. The mixed-species stands are represented by seven genera: Pinus, Quercus, Juniperus, Cupressus, Pseudotsuga, Arbutus, and Alnus. The main species are Pinus durangensis Martínez and Quercus sideroxyla Bonpl. Lately, the Improve forest management combine the forest management and credit carbon offsets according to the Mexican Protocol developed by The Climate Action Reserve [34].

2.2. Dataset description

The dataset came from temporary forest inventory plots with a random sampling design for framework of 30,000 ha. A total of 14,390 temporary forest inventory plots were considered and Durango pine species was selected. A random sample of 1,000 plots was selected in sampling R package [35], and 11,472 measured trees were considered. Firstly, the unsupervised clustering analysis was used for grouping the dataset [1,25] according to the k-Means Clustering of kmeans R package [35]. Ten clusters were generated according with density (N, trees per hectare), basal area (BA, m²), mean diameter (Dm, cm), mean total tree (Hm, m), quadratic mean diameter (QMD, cm), altitudes (A, m), slope (S, %), and aspect (As, categorical variable). All variables were standardized, which perform values bounded 0.0 and 1.0 [21,36,37]. The standardization was performed according to Milligan and Cooper [38] and the Equation 1.

$$Z=(x-Min(x\left)\right)/(Max\left(x\right)-Min(x\left)\right)$$

(1)

where Z is the standardized variable, x is the variable, Min and Max are the minimum and maximum values of x.

In the Clustering analysis, the cluster (k = 10) were generated, and the explained variance was 61.8812%. The data set is shown in Figure 2 for the total dataset and the clusters.

The site-specific variables for clustering analysis are recorded as a descriptive statistic in Table 1. All variables were standardized with Equation 1 to improve the clustering analysis [38].

The total dataset (11,472 pairs of h-dbh relationship) was randomly divided in two sets; 50% for training and 50% for testing or validation. The main statistics for both training and testing dataset are shown in Table 2 for the total tree height and diameter at breast height (h-dbh).

2.3. Nonlinear mixed effect modeling (NLMEM)

The base growth model developed by Richards [30] was used to model the nonlinear h-dbh relationship. This model is based on a sigmoid curve, and it is represented by Equation 2. This model had been extensively used in this kind of relationship [7,11,17,28].

$$h_{ij}=A_0+{\alpha }_0{\left(1-e^{-{\alpha }_{1}db{h}_{ij}}\right)}^{{\alpha }_2}+e_{ij}$$

(2)

where $h_{ij}$ = total tree height j in the plot i; $A_0$ = lower asymptote parameter; ${\alpha }_0$, ${\alpha }_1$ and ${\alpha }_1$ = upper asymptote, growth rate, and slope of growth parameters; e = exponential function; ${dbh}_{ij}$= diameter at breast height j in the plot i; $e_{ij}$ = residual j in the plot i. In this case for $A_0$ parameter, the value 1.3 was fixed. This warranted the total tree height equal to 1.3 m when diameter at breast height is equal to 0 mentioned by Fang and Bailey [39].

For NLMEM, the parameter vector of nonlinear model was defined according to Pinheiro, et al. [40] and summarize as follows [7,11,13] (Equation 3):

$${\mathsf{\Phi }}_{\mathit{k}}={\mathit{A}}_{\mathit{k}}\lambda +{\mathit{B}}_{\mathit{k}}{\mathit{b}}_{\mathit{k}}$$

(3)

where $\lambda$ is the $p\times 1$ vector of fixed parameters ($p$ is the number of fixed parameters in the model), ${\mathit{b}}_{\mathit{k}}$ is the $q\times 1$ vector of random effect associated with $lth$ cluster ($q$ is the number of random parameters in the model), and ${\mathit{A}}_{\mathit{k}}$ and ${\mathit{B}}_{\mathit{k}}$ are the design matrices of size $r\times p$ and $r\times q$ ($r$ is the total number of parameters in the model) for the fixed and random effects specific to each cluster. The residual vector (${\mathit{e}}_{\mathit{i}\mathit{j}}$), and the random effect vector (${\mathit{b}}_{\mathit{k}}$) are frequently assumed uncorrelated and normally distributed with mean zero and variance-covariance matrices ${\mathit{R}}_{\mathit{k}}$ and $\mathit{D}$, respectively.

The upper asymptote parameter (${\alpha }_o$) was treated as a random parameter in the analysis for each cluster (${\alpha }_{ol}$ $l=1,2,\dots ,10$), which explain about the maximum relationship between h and dbh. The random effect vector represent the variability between clusters for the asymptote parameter.

2.4. Artificial neural network (ANN)

The ANN are inspired by the early models of sensory processing by the brain. An ANN can easy be created by simulating a network of model neurons in a computer or specific programming language. Also, by applying mathematical algorithms that mimic the process of real neurons, we can make the network “learn” to solve many types of problems [41]. The ANNs can learn by themselves. Because of they have similarities with the information processing features of the human brain (nonlinearity, high parallelism, capability to generalize), this modeling technique has the potential ability to solve problems that are difficult to formalize, such as problems of biological nature [7,41]. The resilient backpropagation artificial neural network (RBPANN) is a logical adaptative learning scheme, which perform supervised batch learning in multi-layer perceptron. The basic principle of RBPANN is to eliminate the harmful influence of the size of the partial derivative on the weight steps [3,4,19].

According to Barbosa, Costa, Schons, Finger, Liesenberg and Bispo [19], and Anastasiadis, et al. [42], the RBP for ANN employs a sign-based scheme to update the weights in order to eliminate harmful influences of the derivatives´ magnitude on the weight updates. The size of the update step along a weight direction is exclusively determined by a weight-specific “update-value” as follows:

$$∆w_{ij}^k=\left\{\begin{array}{c}-{∆}_{ij}^{k}if\frac{\partial E\left(w^k\right)}{\partial w_{ij}}>0\\ +{∆}_{ij}^{k}if\frac{\partial E\left(w^k\right)}{\partial w_{ij}}<0\\ 0otherwise\end{array}\right.$$

(4)

where $\partial E\left(w^k\right)$ denotes the partial derivative of bathc error with respect to weight $w_{ij}$ at the kth iteration.

The second step of RBP learning is to determine the new update values [19], as follows:

$${∆}_{ij}^k=\left\{\begin{array}{c}{\eta }^{+}{∆}_{ij}^{k}if\frac{\partial E\left(w^{k-1}\right)}{\partial w_{ij}}\frac{\partial E\left(w^{k-1}\right)}{\partial w_{ij}}>0\\ {\eta }^{-}{∆}_{ij}^{k}if\frac{\partial E\left(w^{k-1}\right)}{\partial w_{ij}}\frac{\partial E\left(w^{k-1}\right)}{\partial w_{ij}}<0\\ {∆}_{ij}^{k-1}otherwise\end{array}\right.$$

(5)

where $0<{\eta }^{-}<1<{\eta }^{+}$.

The total number of parameters in RBPANN is five; (i) the increase factor isa set to ${\eta }^{+}=1.2$; (ii) the decrease facto is set to ${\eta }^{-}=0.5$; (iii) the initial updates-value is set to ${∆}_o=0.1$: (iv) the maximum step; (v) the minimum steps.

According to Cartwright [43], the first step in using ANNs is to determine a suitable topology, optimal if possible (number of inputs and outputs, number of hidden layers, and neurons in each layer) and optimal (weights, biases, and activation functions). The process of ANNs begins by setting up the weights as small random variables. Then, each input pattern undergoes a feedforward phase, where the input signal is received and transmitted to all nodes in the hidden layer. In ANNs, every hidden node calculates the sum of its weighted input signals, applies an activation function to determine its output signal, and transmits this signal to the output node. At the output node, the final output signal is computed using the received signals from the hidden nodes [7]. Within the context of RBPANN, the associated error (${\delta }_k$) is computed, and this error is utilized to adjust the weights. The weights correction term is determined based on the error, and it is subsequently employed to update the corresponding weights. Additionally, the ${\delta }_k$ is transmitted to each hidden node. Each hidden node then calculates its error information term by summing the inputs received from the output node, multiplied by the derivative of its activation function [7,43,44]. According to Fausett [45] and Cartwright [43], the general formulation for RBPANN could be as follows (Equation 6):

$$∆w_{ij}\left(t+1\right)=\alpha {\delta }_k{Z}_j+\mu ∆w_{ij}$$

(6)

where $w_{ij}$ is the bias on output unit $k$, $\alpha$ is the learning rate, ${\delta }_k$ is the ratio of error correction weight fitted for $w_{ij}$ that is due to an error at output $O_k$, also the information about the error at unit $O_k$ that is propagated back to the hidden units that feed into unit $O_k$, $Z_j$ is the output activation of the hidden unit $j$ and $\mu$ is the momentum parameter (refers to the contribution of the gradient calculated at the previous time step to the correction of the weights) [43].

The used activation functions for smoothing the h-dbh relationship through RBPANN were tangent hyperbolicus (RBPANN-tanh), softplus (RBPANN-softplus), and logistic (RBPANN-logistic) functions [4,42], these activation functions occur between the hidden layers or between the input layer and hidden layer [17]. These functions were defined for RBPANN-tanh, RBPANN-softplus, and RBPANN-logistic in Equations 7, 8, and 9, respectively. Also, the derivatives are in Equations 10, 11 and 12, respectively.

$$f\left(x\right)=\frac{e^x-e^{-s}}{e^x+e^{-s}}$$

(7)

$$f\left(x\right)=\frac{1}{1+e^{-s}}$$

(8)

$$f\left(x\right)=\mathrm{l}\mathrm{o}\mathrm{g}(1+e^s)$$

(9)

$$f^{\text{'}\left(x\right)}=1-{\left(\frac{e^s-e^{-s}}{e^s+e^{-s}}\right)}^2$$

(10)

$$f\text{'}\left(x\right)=\frac{e^{-s}}{{\left(1+e^{-s}\right)}^2}$$

(11)

$$f\text{'}\left(x\right)=\frac{e^s}{1+e^s}$$

(12)

where $s=\sum w_i{x}_i$ is the information of the node transmits, in which $w_i$ are the weights and $x_i$ are the input values with $s\in \left[-1,1\right]$, $s\in \left[0,\infty \right]$, and $s\in \left[0,1\right]$ for RBPANN-tanh, RBPANN-softplus, and RBPANN-logistic, respectively.

In the ANN learning process, a different vector from 1 to 10 for each hidden layer was performed in a preliminary analysis, and the best results were obtained when the vector was 10 nodes for each hidden layer c(10, 10, 10). In Figure 3, the ANNs plot are presented for vectors c(3, 3, 3), c(5, 5, 5), and c(10, 10, 10) in the hidden layer.

For both RBPANN-tanh and RBPANN-logistic functions, the topology was as follows: (i) two inputs (dbh and cluster); (ii) one output (h); (iii) a vector c(10, 10, 10) hidden layers; and (iv) two nodes for the first layer, 11 nodes for the second, third and fourth layer layers (Bias node is included) and two nodes for the fifth layer. The ANNs for RBPANN-tanh, RBPANN-softplus, and RBPANN-logistic activation functions are presented in Figure 4. Input variables are represented by “I” in nodes, hidden nodes are represented by “H”, input variable is represented by "O”, and Bias nodes are represented by “B”.

For RBPANN-tanh, RBPANN-softplus, and RBPANN-logistic activation functions. The number of repetitions was 10, the maximum steps for training of the NN was 10⁷ , and the threshold was 0.1, which is similar that arguments used by Özçelik, Diamantopoulou, Crecente-Campo and Eler [7] and Shen, Hu, Sharma, Wang, Meng, Wang, Wang and Fu [17]. Also, the training algorithm for ANNs was the resilient backpropagation with weight backtracking [4,42,43]

2.5. Fitting modeling

For NLME, the total tree height and diameter at breast height (h-dbh) relationship was fitted in “nlme” R package [35] and used maximum likelihood estimation method [40] for fixed and random parameters withing cluster groups. While for ANN models, the “neuralnet” R package [35] was used. For ANN, the resilient backpropagation (RPROP) for tangent hyperbolicus (tanh), softplus (softplus), and logistic (logistic) functions were programed for smoothing the result of the cross product of the covariate or neurons and the weights [3,4,45]. All functions about fitting statistics for both training and testing were programmed in R environment [35].

2.6. Models performance criteria

For both training and testing steps, the fitting statistics were obtained in two levels, first one for the entire dataset and second one for each cluster. The statistics were the root mean square error (RMSE), standard error of estimate (SEE), relative SEE (RSEE), fitting index (FI), mean error (E), relative E (RE), Akaike information criterion (AIC), Bayesian information criterion (BIC) and the log-likelihood (logLik). The statistics were computed as follows:

$$RMSE=\sqrt{\frac{\sum _{i=1}^n{(O_i-P_i)}^2}{n}}$$

(13)

$$SEE=\sqrt{\frac{\sum _{i=1}^n{(O_i-P_i)}^2}{n-p}}$$

(14)

$$RSEE=\frac{SEE}{\overline{O}}$$

(15)

$$FI=1-\frac{\sum _{i=1}^n{(O_i-P_i)}^2}{\sum _{i=1}^n{(O_i-\overline{O})}^2}$$

(16)

$$E=\frac{\sum _{i=1}^n(O_i-P_i)}{n}$$

(17)

$$RE=\frac{E}{\overline{O}}$$

(18)

$$AIC=nlog\left(\frac{\sum _{i=1}^n{(O_i-P_i)}^2}{n}\right)+2p$$

(19)

$$BIC=nlog\left(\frac{\sum _{i=1}^n{(O_i-P_i)}^2}{n}\right)+plogn$$

(20)

$$logLik=-nlog\left(\frac{\sum _{i=1}^n{(O_i-P_i)}^2}{n}\right)$$

(21)

where $O_i$, $P_i$ and $\overline{O}$ are observed, predicted and average values of h variable; n = observations; p = number of parameters estimated; and log = logarithm function.

In all cases, the residual values were obtained with the implementation of NLMEM or RBPANN models and the statistics were programmed in R environment [35]. Lately, the NLMEM and ANNs models were ranked based on the overall dataset and cluster-group for all fitting statistics. A ranking system of Kozak and Smith [46] was used. All fitted statistics were equally weighted and Rank 1 was used for the best model and 4 for the poorest.

3. Results

3.1. Training phase

3.1.1. NLMEM

The fitted growth equation for h-dbh relationship by NLMEM performed well and all parameters were significantly different to zero at 5% of significance level. The relationship between total tree height and diameter at breast height can be explained with fixed and random parameters. In Table 3, the estimated parameters and their statistical properties can be found for the entire training dataset. Furthermore, the confidence interval for each parameter is recorded at a 95% confidence level.

The training phase's fitting statistics are presented in Table 4, which includes the overall training dataset as well as the cluster-groups individually. For both overall training dataset, and cluster-groups, the fitting statistics were accurate and showed the potential to offer the NLMEM approach for the h-dbh relationship. The RMSE value for the overall training dataset was 3.1085 m, the best value was 2.4735 for cluster-group 1 (C1) and the worst value for cluster-group 4 (C4). Additionally, the overall training dataset exhibited an E value of -0.0005. The highest value was observed in C3, indicating the poorer performance, while the lowest value was found in C4, indicating the best performance. In terms of AIC, the C4 demonstrated the best performance, whereas C6 exhibited relatively poorer performance.

3.1.1. RBPANN

The results about performed ANNs for RBPANN-tanh, RBPANN-softplus, and RBPANN-logistic activation functions are shown in Table 5. The statistics are in standardized variables and observed in 10 repetitions in the learning process. In this scenario, all three activation functions exhibit favorable outcomes. Specifically, both RBPANN-tanh and RBPANN-softplus deliver comparable performance. In contrast, RBPANN-logistic exhibited the lowest performance among the three activation functions. Additionally, the RBPANN-logistic achieved the minimum number of steps (88) required for convergence.

The fitting statistics for three ANNs applied to examine the h-dbh relationship in Durango pine are presented for both the overall dataset and each cluster-group in Table 6 for training phase. The nine fitting statistics illustrate the accuracy of RBPANN-tanh, RBPANN-softplus, and RBPANN-logistic activation functions in the ANNs models. The topology of each ANN, as depicted in Figure 4, exhibited satisfactory results in predicting h based on dbh and an unsupervised clustering analysis for ten groups. Overall, the estimations demonstrate similar characteristics across the three activation models. However, RBPANN-tanh exhibits certain advantages that are comparable to the other activation functions in the training phase. All ANNs were trained using the resilience backpropagation learning algorithm, and the likelihood function was employed.

The residual and predicted values are presented in Figure 5 for each cluster-group by RBPANN-tanh, which was the best approach in training phase to model the h-dbh relationship. In general, the residuals ranged between -6 m and 6 m. Enhancing the training phase could involve increasing the number of repetitions or epochs, however, the computational process ought to be significantly time-consuming.

Finally, The ranks and sums of ranks for hierarchy of NLMEN and ANNs are presented in Table 7. The statistics for both the overall dataset and the cluster-group were ranked from 1 to 4. In terms of the overall dataset, the RBPANN-softplus exhibited the best performance during the training phase, which was comparable to the NLMEN model and other ANNs. The RBPANN-tanh activation function ranked second, while the RBPANN-logistic ranked third. On the other hand, the NLMEN approach had the lowest rank sum, indicating poorer performance compared to the other models. A similar pattern was observed within the cluster-groups, and the sum of ranks resulted in the following rankings: 176, 241, 269, and 304 for RBPANN-tanh, RBPANN-softplus, NLMEN, and RBPANN-logistic, respectively. It is worth noting that only the RBPANN-logistic demonstrated lower performance compared to the NLMEN approach. The number in parenthesis indicate the ranking for models in the overall dataset.

3.1. Testing or Validation phase

3.1.1. NLMEM

During the testing phase, 5,736 pairs of heights and diameters from 50% of the dataset were utilized. Height estimations were performed using fixed and random parameters for each cluster-group provided by NLMEM approach. The nine testing statistics were computed at two levels: for the overall dataset and for each cluster-group. The Table 8 presents the statistics for testing the advantages of NLMEN for overall dataset and for each cluster-group. All the statistics displayed satisfactory performance, which depended on the number of observations. Additionally, the cluster-groups with limited information exhibited the lowest values. Among the cluster-groups, C4 had the maximum number of observations, whereas C9 had the minimum number of observations.

3.1.1. RBPANN

During the testing phase, the results for the ANNs were similar to the training phase. The ANN utilizing the tangent hyperbolicus activation function (RBPANN-tanh) exhibited the highest performance, followed by RBPANN-logistic, and finally RBPANN-softplus. The Table 9 records the statistics for the testing dataset, both in the overall dataset and within each cluster-group. The FIs values for the overall dataset were higher than 0.7029, when RBPANN-tanh demonstrating the best performance and RBPANN-logistic exhibiting the poorest performance. A similar pattern was observed for other statistics such as AIC BIC and logLik. Furthermore, in this instance, the ANNs demonstrate superior performance compared to the NLMEM approach.

The Figure 6 displays the representation of residual and predicted values for each cluster-group obtained through RBPANN-tanh. In this scenario, the residual dispersion appears to be larger compared to the training phase. However, this can be attributed to the utilization of a new dataset, where the predictions are made under different training conditions. Cluster-groups 4 and 10 exhibited higher levels of dispersion compared to the remaining cluster-groups.

Finally, the ranks and sums of ranks for hierarchy of NLMEN and ANNs are presented in Table 10. During the testing phase, the RBPANN-tanh achieved the highest rank of 1 with a sum of ranks of 9, while the NLMEM approach performed the poorest, ranking 4 with sum of ranks of 36. The RBPANN-softplus and RBPANN-logistic exhibited relatively similar performance conditions. Among the proven models of ANNs, the RBPANN-logistic exhibited the poorest performance. In terms of the sum of ranks for the combined overall dataset and cluster-groups, the RBPANN-tanh demonstrated the best performance with a sum of ranks of 181. Following that, the RBPANN-logistic, RBPANN-softplus, and NLMEM were ranked 2nd, 3rd, and 4th, with sum of ranks of 201, 240, and 368, respectively.

4. Discussion

Having knowledge about the total tree height and diameter at breast height is essential for both the development and application of many growth and yield models. Models focusing on the h-dbh relationship serve as valuable tools for accurately predicting tree height based on dbh measurements. Because of the dbh can be conducted quickly, easily, and accurately, but the measurement of total tree height is comparatively complex, time consuming, and expensive [11]. The NLME had been a capable approach to generate models in h-dbh relationship for different species and assumed fixed and random parameters for specific-groups or covariables to study the variability inter-and intra-plots, ecological regions or cluster-groups [10,16,39]. Also, these models have been studied for local and generalized formulations with NLMEM approach [12,13,16,28]. In this case of study the NLMEM performance was accurately strong to model the h-dbh relationship for Durango pine and the inclusion of unsupervised clustering analysis improve the estimated parameters and its statistics properties [36,47], which involve fixed parameter for the overall dataset in training phase and random parameter for each cluster-group, also parameter to give information about general variability and variability within cluster-group.

The NLMEM demonstrated outstanding performance during the training phase, with the fitting process converging quickly and effortlessly. Additionally, the maximum likelihood approach yielded favorable and suitable results particularly when expressing the asymptote parameter with mixed effects (Table 3 and Table 4). All parameters in fitting process were significantly different to zero at 5% of significance level and the random parameters allow suitable estimations in training phase and those were used for cluster-groups in testing phase. The application of the NLMEM approach on the testing dataset resulted in successful outcomes that aligned with the expected results (as shown in Table 8), accompanied by the utilization of appropriate statistical measures. As an illustration, the root mean square error (RMSE) for the overall dataset during the testing phase was determined to be 3.1438 m, with an average value of 3.3773 m observed within the cluster-groups (refer to Table 8). By employing a mixed-effect model and incorporating cluster-group inclusion, the Chapman-Richards growth equation [30] (Equation 2) proves to be a highly effective model for predicting the height of Durango pine trees. Similar results have been conducted for several species an different conditions [11,16,28]. Even though the NLMEM method is accurate for height prediction based on diameter measurements, it is worth considering that ANNs could be a suitable alternative for modeling the h-dbh relationship under several dataset conditions and the incorporation of grouped strategies [7,14,48]. In recent times, there has been a growing application of AI and ML techniques in the fields of biology and forestry. These advanced approaches have proven valuable in addressing challenges that require substantial computational resources and unsupervised learning methods [1,41], Several of these approaches have been employed in studying the height-diameter at breast height (h-dbh) relationship, leading to notable outcomes and reported successes for various species and under diverse forest management conditions, demonstrating their versatility and effectiveness [7,14,15,17,48]. In this context the ANNs model outperformed the NLMEM approach.

In this study, the ANNs were evaluated and compared with the traditional NLMEM method. The ANNs utilized the RBP learning algorithm along with three activation functions. In most cases, the ANNs employing RBPANN-tanh, RBPANN, and RBPANN-logistic (Equations 7, 8, an 9, respectively) exhibited superior performance compared to the results obtained by NLMEM, both during the training and testing phases. The training statistics for three ANNs, as presented in Table 6, exhibited enhanced fitting performance compared to the statistics obtained by NLMEM (see Table 4). This improvement was observed in both the overall dataset and cluster-group analyses. These findings provide evidence that the clustering analysis using the k-means algorithm effectively grouped the dataset utilized in this study [36,47]. The RBPANN-tanh model, employing a tangent hyperbolic activation function, demonstrated the highest performance in predicting height measurements during both the training and testing phases (it can see in Table 6 and 9). Furthermore, the ranks and sum of ranks, based on the ranking system proposed by Kozak and Smith [46], provided evidence of the advantages of the ANNs models over the NLMEM approach. Models such as RBPANN-logistic were reported by Özçelik, Diamantopoulou, Crecente-Campo and Eler [7] revealed that models such as RBPANN-logistic exhibited advantages over NLMEM models when predicting the growth of Crimean juniper in the southern and southwestern regions of Turkey. Conversely Shen, Hu, Sharma, Wang, Meng, Wang, Wang and Fu [17] developed ANNs models utilizing RBPANN-logistic and RBPANN-tanh transfer or activation functions for Populus spp. L. in China, where the RBPANN-logistic model outperformed both the NLMEN and the RBPANN-tanh model. In our case, the best ANN was the RBPANN-tanh and this outperformed other tested ANNs and NLMEM approach. Similar results have been reported the advantages of ANNs or deep learning algorithms over the ordinary least square and NLMEM models in both training and testing or validation phases [14,15,48,49]. In all cases, the implementation of ANNs exhibited significant advantages over traditional approaches when modeling the h-dbh relationship.

In this study, based on the implemented ranking system, the RBPANN-tanh model emerged as the top performer (residual and predicted values are showed in Figure 5 and Figure 6). It achieved a sum of ranks of 176 for the training phase and 81 for the testing phase. These sums of ranks account for both the overall dataset and cluster-groups, as illustrated in Table 7 and Table 10, respectively. In terms of training, the RBPANN-softplus model ranked second, whereas during the testing phase, the RBPANN-logistic model exhibited the second-best performance. On the other hand, the RBPANN-logistic model performed least effectively in the training phase, while the NLMEM model demonstrated comparatively lower performance during the testing phase. The ANNs developed in this study, as depicted in Figure 4, were trained using the RBP algorithm. The ANNs were then evaluated using three different activation functions: RBPANN-tanh, RBPANN-softplus, and RBPANN-logistic. These models comprised a total of five layers, including three hidden layers. The training process involved ten repetitions to ensure robustness and accuracy. Even though the RBPANN-logistic converging in 88 steps, it exhibited relatively poorer performance compared to the RBPANN-tanh, which achieved better results within 301 steps. Interestingly, the RBPANN-logistic required a longer convergence time of 1885 steps, indicating its comparatively poorer performance in this aspect. As a result, the developed ANNs model showcased a high capability for predicting total tree height measurements. This highlights the potential application of AI in modeling the h-dbh relationship, not only for Durango pine but also for general forest modeling purposes or other variables [6,19,21,50]. The ANNs could be used to improve the estimations in forest inventory and forest management and planning in mixed-species forest in Durango, Mexico.

5. Conclusions

The nonlinear mixed effect modeling (NLMEM) and Artificial Neural Networks (ANN) with resilience backpropagation (RBP) were employed to model the height-diameter at breast height (h-dbh) relationship for Durango pine species. Unsupervised clustering analysis was conducted to enhance the capability of the trained and tested models.. Three activation functions, namely tangent hyperbolicus (RBPANN-tanh), softplus (RBPANN-softplus), and logistic (RBPANN-logistic), were utilized in RBPANN. Those activation functions were trained and tested on both overall dataset and each cluster group. In general, the ANNs outperformed the NLMEM for predictions of heights in training and testing phases. The best model in both training and testing phases was the RBPANN-tanh, which assumed five layers in the ANN and three of them were hidden. The use of ANNs proves to be a suitable and effective approach for estimating the total tree height of Durango pine species. Additionally, incorporating unsupervised clustering analysis enhances the estimation accuracy in ANN models, highlighting the capabilities of artificial intelligence (AI) in this context. In conclusion, AI techniques such as ANNs prove to be suitable and modern statistical tools for forest modeling.