Investigation of parametric and artificial neural network modeling approaches for total tree height prediction in cedar plantations by

: Background: Height-diameter relationships are of critical importance in tree and stand volume estimation. Stand description, site quality determination and proper forest management decisions originate from reliable stem height predictions. Methods: In the context of this work, the prediction ability of the developed height-diameter models was investigated for cedar ( Cedrus libani A. Rich.) plantations in Western Mediterranean Region of Turkey. Towards this direction, parametric modeling methods such as fixed-effects, generalized models, and mixed-effects were evaluated. Furthermore, in an effort to come up with the construction of more reliable stem-height prediction models, artificial neural networks were employed using two different modeling algorithms: the Levenberg-Marquardt and the resilient back-propagation. Results: Taking into account the prediction behavior of each respective modelling strategy while using a new validation data set, the mixed-effects model with calibration using 3 trees for each plot seems to be a reliable alternative to the rest standard modelling approaches given the evaluation statistics regarding the predictions of tree heights. Conclusion : Finally, as for providing the most reliable results as compared to the remaining, the resilient propagation algorithm showed its capability of providing more accurate predictions of the tree stem height and thus it can be a reliable alternative to pre-existing modelling methodologies.


Background
Cedar (Cedrus libani A. Rich.) is an important tree species both ecologically and economically and also it is presently found primarily in the Taurus Mountains of Turkey (Boydak 2004). The most recent inventory data show that cedar covers an area of around 463,521 hectares in Turkey, of which a total growing stock is 27.4 million m 3 (GDF 2015). It is not only an important raw material resource of forest products industry, but also fulfils critical ecological tasks such as preserving soil and water resources, mitigation of and adaptation to the adverse impacts of climate change, and provides a vital for maintaining biological diversity in Turkey. It has moderate soil requirements and is tolerant to temperature and drought in summer and to the cold temperatures in winter (Saatçioğlu 1979); therefore, it has a big afforestation potential especially in areas with semi-arid climate. Given that 35% of Turkey's territory has semi-arid climate, cedar is a critical tree species that can be used for afforestation. As of 2000, an area of around 110.000 ha was afforested with cedar (Konukçu 2001).
Globally forest areas continue to decline in the last 30 years while plantations have increased in the same timeframe and they now comprise 7% of the total global forest area.
These plantations are of great importance in economic and ecological terms. Thus, it is essential to have robust and reliable information about the growth and development characteristics of these species for sustainable management of plantations and to measure the success of the afforestation efforts. Diameter at breast height (d, 1.3 m above ground level) and total stem height (h) are the two basic parameters for the inventory, planning and management of plantations. These parameters are used as the fundamental variables in several forest practices such as prediction of standing yield, growth projection, carbon accounting, and identification of stand structural diversity, determination of various damages in forest or stands and to predict missing tree heights (Curtis 1967;Peng et al. 2001;Newton and Amponsah 2007;Adame et al. 2008).
Although, there are published papers dealing with the tree height estimation and prediction models using mainly standard modeling techniques (Colbert et al. 2002;Newton and Amponsah 2007;Zhang et al. 2014;Kearsley et al. 2017), there is still space for further research about topic when evaluating more recent modeling methods such as artificial intelligence. Existing research on height-diameter (h-d) models for cedar plantations in Turkey was contacted by Çatal (2012). However, the aforementioned research was based on limited number of traditional h-d models, without any further investigation of the efficiency of other modeling approaches in the same topic.
In forest management applications or field inventories, d in trees can be measured rapidly, easily, and precisely while h is more time consuming, challenging and costly to measure tree height (Sharma and Parton 2007;Meng and Huang 2009;Corral-Rivas et al. 2014;Özçelik et al. 2018). It is also known that there is a strong correlation between d and h in the stem of trees. Because of this situation, models can be created to discover the relationship between d and h based on the tree diameter at breast height and the total stem height of trees that are measured. Finally, these created models are used to predict the height of any tree in the stand with unknown height (Huang et al. 1992;Robinson and Wykoff 2004;Lei et al. 2009).
Generally, h-d models can be defined as simple and generalized. Simple models express height as a function of tree diameter at breast height (Bronisz and Mehtätalo 2020). This modeling approach represents the best equation that can be obtained for specific stand and point in time. The main problem with this approach is the large sampling effort required (Gómez-García et al. 2014). On the other hand, as stated by Lhotka (2012), research evaluating the correlation between tree diameter and total height has found that this relationship varies among species and can be influenced by stand-specific variables such as stand structure and site productivity. Models that include measures of stand structure, site productivity and tree social position are commonly referred to as generalized h-d models. In several studies, the inclusion of additional stand-level variables improved tree height estimates (López Sánchez et al. 2003;Sharma and Zhang 2004;Castedo-Dorado et al. 2006;Vargas-Larreta et al. 2009). For example, Newton and Amponsah (2007) utilized stand dominant height as a measure of standlevel competition. Temesgen et al. (2008) stated that stand-level basal area per hectare and crown competition was useful variables for predicting tree height.
Nowadays, there are alternative modeling approaches such as the generalized, fixed and mixed-effects models, artificial neural network techniques, etc that can help to this direction by modeling the relationship between h-d and finally reducing the measurement effort on the ground.
In Turkey, h-d models were developed for some tree species at regional scale (Diamantopoulou and Özçelik 2012;Özçelik et al. 2014;Çatal and Carus 2018). However, research on the development of h-d model using mixed-effects approach for cedar plantations are very limited (Ercanlı 2015;Özçelik et al. 2013;Özçelik et al. 2018).
Until now, in the field of forest modeling, the implemented research led to the conclusion that artificial neural network models (ANNs) can be considered as a significant alternative modeling technique for many characteristics of trees against to standard modeling methods (Liu et al. 2003;Özçelik et al. 2008;Soares et al. 2013;Imada 2014;Ercanli et al. 2018;Socha et al. 2020). Moreover, according to the same research field, ANN models have gained publicity because they have successfully been applied for classification (Schmoldt et al. 1997;Sarigul et al. 2003;Cosenza et al. 2017) and estimation and prediction problems (Pertsen et al. 2010;Leite et al. 2011;Reis et al. 2016;Özçelik et al. 2017; Monteiro da Zhou et al. 2018;Socha et al. 2020). Despite the fact that ANNs need effort for avoiding over or under-fitting the data during their training phase, artificial intelligence has been successfully used for total tree height models construction (Li and Jiang 2010;Diamantopoulou and Özçelik 2012;Özçelik et al. 2013;Vieira et al. 2018;Thanh et al. 2019; Ercanli 2020), as well.
The majority of the research efforts of the literature employ multilayer perceptron architecture combined by the standard backpropagation algorithm for determining the optimal weights. The objective of this study was to investigate the effectiveness of different modeling approaches in terms of providing accurate predictions regarding the total stem height of the cedar tree plantations in Western Mediterranean Region of Turkey. Upon fitting the constructed models and in an effort to enhance the feasibility of our approach, we use easy to obtain measurements from the forest environment, such as the breast height tree stem diameter.
Parametric modeling methods such as the fixed-effects, generalized, adjusted fixed-effects and mixed-effects models after being localized using calibration data of one to three sample cedar pine trees were compared against the prediction ability derived by the constructed artificial neural network models. Two different learning algorithms, the Levenberg-Marquardt (LMANN) and the resilient propagation (RPANN) were examined towards optimizing the learning procedure and thus better harness the information residing in real data measurements originating from the forest ecosystem.

Data
This study was carried out in the Taurus cedar plantations in Western Mediterranean Region of Turkey. The sampling plots are within the Isparta Regional Directorate of Forestry, which includes Isparta and Burdur provinces. The elevation of the sampling plots from the sea level ranges from 1000 to 1600 m. The climate of the study site is characterized as transition climate from the Mediterranean climate to continental climate. It has a variable precipitation regime, while the annual precipitation ranges from 426 to 814 mm, the annual mean temperature in distribution areas of ranges from 11.9 to 13.2 °C (MGM 2013). The bedrock of the study site is composed of sedimentary rock such as limestone, sandstone and clay stone.
The data was collected from 30 sampling plots randomly distributed in the site in order to identify the diameter-height associations of the plantations in the Western Mediterranean Region. The sample plots were selected to represent the existing range of ages, stand densities and sites for the tree species in Turkey. Plot size ranged from 270 to 540 m 2 , depending on stand density, in order to achieve a minimum of 30 trees per plot. The age of the sampling plots varied from 17 to 35. For each tree, two perpendicular diameters (outside-bark 1.3 m above ground level) were measured to the nearest 0.1 cm and were then averaged to obtain diameter at breast height (d, cm). Total heights of these trees were measured to the nearest 0.5 m with a Blume-Leiss hypsometer. The dominant height (H0) and dominant diameter (D0) were defined as the average height of the 100 largest-diameter trees per hectare, depending on plot size.
The data were randomly divided into two groups; totally of 21 sample plots (70% of all sample plots) were randomly selected as the model fitting dataset. The remaining 9 sample plots were used as the model validation dataset for evaluating model performance when predictions were made (Table 1). Summary statistics for the fitting and the validation data sets are given in Table 1. Moreover, both the data sets are graphically depicted in Figure 1.

Generalized h-d models
One practical alternative is to develop a generalized h-d relationship that describes the relationship between d and h and includes stand variables as predictors (e.g. dominant height, quadratic mean diameter, dominant diameter, number of trees per hectare, stand basal area, etc.).
In this study, 21 generalized h-d models selected and evaluated from previous studies (López- Sánchez et al. 2003;Sharma and Zhang 2004;Castedo-Dorado et al. 2006;Sharma and Parton 2007;Crecente-Campo et al. 2010) were fitted using non-linear least squares. These models include generalizations of the most flexible equations for h-d relationships (Zhang 1997;Peng et al. 2001).

Fixed-effects model
A large number of nonlinear model forms were evaluated for cedar plantations; including those reported by Huang et al. (1992) and Lei et al. (2009) as Weibull, Exponential, Chapman-Richards, Gompertz, Schnute, and Korf-Lundgvist. The nonlinear models were fitted by the use of ordinary nonlinear least-squares (ONLS) to the fitting data set. Among the candidate models, the Chapman-Richards produced the best performance to model the h-d relationship for the cedar plantations.
where hij and dij are, respectively, total height (m) and diameter at breast height (cm) of the jth tree in the ith plot, 1 -3 are model parameters, and ij is random error. Temesgen et al. (2008) showed that when the heights of a subsample of nim trees from the i th stand are known, the predicted heights of the remaining trees from the same stand can be calibrated by use of the following correction factor:

Calibration of the fixed-effects model
where k * is correction factor, ĥ ℎ is predicted height from equation (1), and ℎ is observed height. Then the predicted height for a tree from the same stand can be calibrated as follows: The effect of the correction factor on the predictive performance of a nonlinear fixedeffects model was also investigated in this study.

Nonlinear mixed-effects model (NLME)
In the mixed-effects framework, all parameters of Eq.
(1) can be expressed as fixedeffects parameters (common to all trees), and some or all parameters contain additional random components, which are specific to individual plots. Eq. (1)  where R and D are diagonal matrices, assuming that the and are independent. Procedure NLMIXED from SAS (SAS Institute Inc. 2010) was used to obtain fixed-and random-effects parameters of Eq. (4).
The random parameters for plot can be predicted by the use of the first-order Taylor series expansion : where ̂ is estimate of the random parameters for tree at the kth iteration, ̂ is estimate of D, the variance-covariance matrix for , = ( , , ) | ,̂, R is the variance-covariance matrix for , is the × 1 vector of observed heights, and is number of tree measurements used in localizing the height growth model.
An iterative procedure was needed to estimate . Using a null starting value (̂= ), Eq.(5) was repeatedly updated until the absolute difference between two successive iterations was smaller than a predetermined tolerance limit. The result was an approximation of the empirical best linear unbiased predictor (EBLUP) for random effects.

ANN models
The most significant advantage of an artificial neural network model is the fact that it is trained from the measured data in hand without any need to add further information. After the system fitting with the real data, the ANN model discovers automatically the existing dependencies, and finally produces the trained model. For this purpose, it is very important the net topology and the learning algorithms, under specific rules, to be specified by the modeler.

Levenberg-Marquardt artificial neural network model (LMANN)
As an intermediate algorithm between the steepest descent method and the Gauss-Newton algorithm, the Levenberg-Marquardt (LM) algorithm (Levenberg 1944;Marquardt 1963) was selected to be applied in the multilayer perceptron learning of ANN. Past studies showed that the (LM) algorithm has the ability on the one hand to overcome the known The original description of the Levenberg-Marquardt algorithm is given in Levenberg (1944) and Marquardt (1963) papers. A detailed description of the application of Levenberg-Marquardt algorithm to neural network training is given in Hagan and Menhaj (1994) paper. In a few words LM algorithm introduces the approximation to Hessian matrix (Hm) as (Özçelik et al. 2019): where J: is the Jacobian matrix that contains first derivatives of the network errors with respect to the weights and biases, μ: is the combination coefficient that is always positive and I: the identity matrix.
The update rule of the LM algorithm is expressed by the function wt+1=wt -(Hm where w: are the weights, and e: are the biases. The effectiveness and convergence of the LMANN models are very sensitive to the adjustment of the combination coefficient (μ) of equation (6). As the value becomes larger, the more weight is given to gradient descent learning with a small step size. The smaller it is the more weight is given to large step sizes. For μ equal to zero, the algorithm turns to the Gauss-Newton method. Given the fact that the Gauss-Newton method's more accurate behavior near an error minimum, the aim was to shift toward the Gauss-Newton method as fast as possible (Diamantopoulou et al. 2015).

Resilient propagation artificial neural network model (RPANN)
The resilient propagation, known as Rprob, is an interesting neural network algorithm, which has been completely described (Riedmiller and Braun 1993) and successfully used (Prasad et al. 2013;Saputra et al. 2017;Florescu and Igel 2018) in order for known disadvantages of the standard back-propagation algorithm to be overcome. Back propagation algorithm has known drawbacks. As the learning of the standard back propagation is more or less elementary, the main problems of its learning face are the slow convergence because of small changes in weights, and the neural network weights trapping around local optima.
The resilient propagation algorithm mainly refers to the direction of the gradient. It calculates a delta value (Δij) for each weight that increases when the gradient doesn't change sign under the prerequisite that the system works on the right direction or decreases when the gradient does change sign. The following learning rule is applied to calculate delta value which improves the change of the network weights (Riedmiller and Braun 1993;Prasad et al. 2013): where, 0 < − < 1 < + Finally, the new weight value between i and j neurons in two consecutive layers on the (t-1) is given by:

Some aspects of ANN models' training
In geotechnical engineering, the selection of the proper input variables for ANN model building is usually based on a priori knowledge of the physical problem (Maier and Dandy 2000). In our case, in order to construct a reliable model for total tree stem height prediction, the choice of the input variables introduced to the input layer of the ANNs was based to sensitivity analysis (Cariboni et al. 2007). Sensitivity analysis was performed to the available variables were selected in advance due to the physical problem to be measured in the field.
In an effort to avoid overfitting and undertraining of the network reaching the best possible training efficiency of the model, the number of hidden nodes in the hidden layer was investigated by starting the training with the smallest possible network, without any hidden nodes. Then the construction of the model continued by adding any necessary hidden node in order the network learning to be deemed complete and the architecture finalized when satisfactory agreement was obtained between network predictions and target outputs.
Overfitting avoidance of the ANN models led to acceptable generalization ability of the models. For this reason, in order to have both training and testing data sets (Leahy 1994), due to its efficiency and easiness to apply, the k-fold cross validation resampling method with k=10, was used (Olson and Delen 2008;May et al. 2010). As described in "Data" section, the whole data set was randomly divided into fitting and validation data sets. The fitting data set was used for the construction of all different type models, while the models in their construction phase never saw the validation data set. Validation data set used only for the access of the performance of all models. Especially for the ANN models construction, the 10-fold cross validation division was used for the fitting data set that was further repeatedly divided into training (90% of the fitting data) and testing (10% of the fitting data) data sets as follows ( Figure 2).
Due to the fact that the number of weights is strongly related to the ANN model error value, the reliability of the ANN model and its ability to generalize was measured through the error that was estimated as the average error rate on these cross-validation examples (Fig. 2).
As mentioned, the effectiveness and convergence of the training of the LMANN models significantly depends on the proper value selection for the combination coefficient (μ) (Eq. 6).
Following the procedure described in Diamantopoulou et al. (2015), the initial value of (μ) was set to 0.001 and the final value of (μ) that led to the lowest sum of square errors model value was found equal to 1.0000e-07.
Similarly to the LMANN models, the training of the RPANN models is very sensitive to the value of delta rule ( ( ) ) (Eq. 7). The value of this rule is responsible for the sign (positive or negative) of the gradient to show the direction of the adjustment weight. Through this process the iteration of learning stops when the error target is reached. The initial delta value was 0.07, with increment step 1.2 and decrease step 0.5.
The architecture of both different algorithms LMANN and RPANN, finalized when the proper number of hidden nodes in the hidden layer was obtained by the trial and error process, starting from one hidden node until the desired target error value was reached. During the application of trial and error methodology, the selection of the proper combination of the transfer functions was made. The combination that boosted the models' learning to the most trustworthy behavior was a) the hyperbolic tangent transfer function (tansig) that transfers the information from the input layer to the hidden layer and b) the linear activation function (purelin) that transfers the information from the hidden layer to the output layer (Fausett 1994;Beale et al. 2014).
The neural network toolbox of the Matlab software (Beale et al. 2014) developed by The Math Works TM Inc, (R2011b version) that provided by the Aristotle University of Thessaloniki, Greece, was used for the artificial neural network models building.

Evaluation
Different modeling approaches evaluated in this study: (1) fixed-effects, (2) generalized, (3) ANN, (4) adjusted fixed-effects, and (5) calibrated mixed-effects techniques. According to the last two methods, parameters of the adjusted fixed-effects and mixed-effects models were "localized" by use of sampled heights in each plot, and then applied to predict all tree heights in the plots. As recommended by Trincado et al. (2007), three sampling scenarios were considered, corresponding to number of tree heights (ranging from one to three) measured in each plot. In order the generalization ability of the constructed models by the different modeling approaches to be compared and evaluated, the following statistics computed for the validation data set (Özçelik et al. 2018;Bronisz and Mehtätalo 2020): Root mean square error : where n is number of plots, ni is number of tree measurements for plot i, hij and ĥ are observed and predicted values of tree height, respectively, and ℎ ̅ is average value of hij for plot i.

Nonlinear regression models (Generalized models, fixed-and mixed-effects models)
As stated earlier, twenty-one generalized h-d models selected from previous studies and then the models were fitted using ONLS method using the fitting data set while the predictive capability of the models were evaluated using the validation dataset. The final generalized-h-d model (Krumland and Wensel 1988) including dominant diameter and dominant height as additional independent variables, resulting in: Different combinations of mixed parameters from equation (1) were implemented ( Table 2).
The results of Table 2 shows that model (1) with random components 1 and 3 produced the lowest values of Akaike's information criterion (AIC) (Akaike 1974) and Bayesian information criterion (BIC) among various combinations of mixed parameters.
The final mixed-effects model was: where 1 and 2 are random parameters.
Parameter estimates for the best models (generalized fixed-effects, and mixed-effects) were obtained using the fitting data set (Table 3).
Compared to fixed-effects models, mixed-effects models require prior height information, which is not necessary for the case of the fixed-effects models. The latest models consider all observations as independent and the residuals are measured against the population average. On the contrary, the residuals of the model including random parameters take into consideration tree-specific predictions (Bronisz and Mehtätalo 2020).

Artificial neural network modeling approaches
The aspect of the ANN models was the three-layer architecture that includes a) one with two hidden nodes as resulted by the error and trial procedure that has been followed. The value of the root mean square error was found equal to 0.5387. For both ANN algorithms, the optimum number of hidden nodes in their hidden layers was determined after the examination of different number of nodes. The best values were selected as the values that led to the best learning of the ANN models, with the smallest root mean square error value (Figure 3).
Moreover, the k=10 cross validation method applied in the fitting data set for both ANN algorithms for the models training and testing in their construction phase using the fitting data set, ensures the absence of overfitting and their ability to generalize. Both networks were trained for as many epochs as it needed in order the models to reach to the minimum mean square error value ( Figure 4). As can be seen (Figure 4), the estimation errors for the training and the test data sets were not significantly different for both ANN models, fact that indicated the absence of overfitting in the modeling systems. The LMANN was stopped after 200 training interactions (epochs) as there was no influential reduction of the mean square error value after the 28 epochs, while the RPANN was trained for 500 interactions as there was no influential reduction of the mean square error value after the 167 epochs (Fig. 4).

Evaluation of the different modeling approaches
According to the fitting and the validation data set, the evaluation statistics calculated for the different modeling approaches and their values are shown in Table 4. The best results among the generalized h-d models were obtained from Krumland and Wensel (1988) are presented in Table 4. Comparing the MD, MAD, RMSE, and FI values (Table 4) On top of the results of Table 4, the accuracy of the estimation results of the ANN models was examined by the correlation coefficient (r) between the h values and their estimations by both ANN models. According to the estimations obtained by the LMANN model, the r value was equal to 0.9223, while the r value for the estimations obtained by the RPANN model was equal to 0.9259. The variability of the estimation errors in relation to the mean of the total tree heights (CV%) in the fitting data set was 9.00% and 8.80% for the LMANN and the RPANN models, respectively. The reliability and the predictive capabilities of both ANN models was also examined based on the new validation data set consisting of 270 trees in 9 plots in the same area, never seen by the constructed models before. In order to further investigate the ability of our models to generalize, r, and the coefficient of variation (CV%) between the h values and their predictions by both ANN models were calculated.
According to the predictions obtained by the LMANN model, the r value was equal to 0.9318 and the CV% was 10.27%, while the r value for the RPANN model was equal to 0.9336, with CV% equals to 10.02%.

Standard modeling techniques assessment
The adjusted fixed-effects model had a much lower mean deviation (MD) for all calibration alternatives compared to the fixed-effects model. Considering all calibration alternatives, the fit index increased at variable rates from 27% to 39%. The mean increase was 35%. Similar findings were also reported by Lhotka (2012), VanderSchaaf (2013) and Özçelik et al. (2018). In this study, the most remarkable prediction improvement among the tested calibration alternatives was obtained through the use of three trees for calibration in each sampling plot. Different sampling alternatives tested for both adjusted fixed-effects model and NLME model, corresponding to number of trees measured in each plot. These tree heights were used to estimate random parameters of NLME model and also to calibrate heights and calculate the evaluation statistics for the validation data set (Table 4).
The mixed-effects model had lower mean deviation values obtained through the calibration with 3 trees compared to the adjusted fixed-effects model. Considering all calibration alternatives, there was a decrease of 32% and 39% on average in the mean deviation and mean absolute difference with the use of mixed-effects model, respectively. Previous studies (Trincado et al. 2007;Temesgen et al. 2008;Huang et al. 2009;VanderSchaaf 2013;Gómez-García et al. 2014;Özçelik et al. 2018) also demonstrated that the calibration improved height predictions significantly.  (2013) and Özçelik et al. (2018). However, the success of prediction started decreasing substantially after a certain number of trees. Therefore, the number of trees to be used for calibration should be selected in a way to strike a balance between the inventory-taking cost and the prediction accuracy and success. Calama and Montero (2004) suggested that calibration should be performed with additional trees to increase the success of height predictions. Trincado et al. (2007) argued that the accuracy of height predictions increased with the increased number of trees to be used for calibration but that increase was not remarkable. We also found a similar result in our study and the prediction success in the calibration with the use of three. Similar results were also reported by Temesgen et al. (2008) and Huang et al. (2009) used a series of 1 to 9 trees for calibration but the best result in height predictions was achieved in the calibration performed with 1 tree. Adame et al.
(2008) tested different calibration alternatives and reported that it was necessary and adequate to use 2 or 3 trees for calibration. Generally, this study reported that there was a minimal improvement in prediction capability of the model depending on the increasing number of trees beyond the number of three trees used for calibration.

Artificial neural network modeling techniques assessment
Due to their ability to overcome known disadvantages and limitations of the standard back-propagation algorithm (Rumelhart et al. 1986;Fausett 1994;Patterson 1996;Wilamowski and Yu 2010) which is the most used optimisation technique applied for the multilayerperceptron training of ANNs, the Levenberg-Marquardt (LM) and the Resilient Propagation (RP) algorithms were used. According to the proximity of each point to the 45-degree line throughout the range of the measured total heights (Fig. 5b) indicates that both constructed ANN models were showed similar adjustment to the measured data. In addition, the RPANN model appears to provide slightly more accurate adaptation as compared to the LMANN model ( Fig. 5a).
Moreover, taking into account the statistical exploration of the calculated estimation errors derived by the ANN models, having the errors frequency peak at zero, quickly falling off at the two sides of smaller error values, both ANN models can be considered as well-trained healthy networks.
Finally, the results of Table 4 Figure 6, the prediction errors histogram and the relative errors box plots showed more or less similar effectiveness of the predictive ability of the LMANN and the RPANN models, with slightly less error variability of the stem height predictions as derived from the constructed RPANN model.
Finally, according to the results of Table 4, upon comparing all modeling approaches, the resilient propagation artificial neural network model showed a clear estimation and prediction superiority over fixed, generalized, fixed-calibrated and mixed calibrated models.

Conclusion
In the context of this work, height prediction models were developed for cedar pine plantations in Western Mediterranean Region of Turkey, which involve the evaluation of five alternative methods; (1) fixed-effects model, (2) generalized h-d models, (3) adjusted fixedmodel, (4) mixed-effects model, and (5) ANN models. Except for the modelling process itself, this methodology in the cases of the adjusted and mixed-effects models also employs the calibration of the parameters of the fixed-effects and mixed-effects models using a subset of height measurements, ranging from 1 to 3 sample trees per plot. As compared to the fixedeffects, all constructed models (generalized, adjusted fixed-effects, mixed-effects, LMANN and BPANN models) achieve better performance given the fact that they provide more accurate and reliable tree height predictions. Given the results, it is obvious that the performance of the adjusted and mixed-effects models improved with increasing sample size. However, increasing the number of trees used in the calibration had a minimal effect on improving the predictive ability of the models.
As for the constructed ANN models, it seems that both alternatives exhibit similar performance in their construction phase as well as in terms of their ability to generalize. In addition, based on the cumulative results regarding all evaluated approaches, it is obvious that the resilient propagation algorithm is more effective in terms of predicting the total stem height of the cedar tress originating from Western Mediterranean Region of Turkey. Conclusively, it is worth noting that the effectiveness and efficiency of ANN modelling in general and in specific the resilient propagation algorithm along with its fit to the current dataset provides raises the potential of its applicability to further modelling cases that lie in the area of forest management.
Such directions constitute interesting future research endeavours.