Machine Learning for the Estimation of Diameter Increment in Mixed and Uneven-Aged Forests

Abstract

Estimating the diameter increment of forests is one of the most important relationships in forest management and planning. The aim of this study was to provide insight into the application of two machine learning methods, i.e., the multilayer perceptron artificial neural network (MLP) and adaptive neuro-fuzzy inference system (ANFIS), for developing diameter increment models for the Hyrcanian forests. For this purpose, the diameters at breast height (DBH) of seven tree species were recorded during two inventory periods. The trees were divided into four broad species groups, including beech (Fagus orientalis), chestnut-leaved oak (Quercus castaneifolia), hornbeam (Carpinus betulus), and other species. For each group, a separate model was developed. The k-fold strategy was used to evaluate these models. The Pearson correlation coefficient (r), coefficient of determination (R²), root mean square error (RMSE), Akaike information criterion (AIC), and Bayesian information criterion (BIC) were utilized to evaluate the models. RMSE and R² of the MLP and ANFIS models were estimated for the four groups of beech ((1.61 and 0.23) and (1.57 and 0.26)), hornbeam ((1.42 and 0.13) and (1.49 and 0.10)), chestnut-leaved oak ((1.55 and 0.28) and (1.47 and 0.39)), and other species ((1.44 and 0.32) and (1.5 and 0.24)), respectively. Despite the low coefficient of determination, the correlation test in both techniques was significant at a 0.01 level for all four groups. In this study, we also determined optimal network parameters such as number of nodes of one or multiple hidden layers and the type of membership functions for modeling the diameter increment in the Hyrcanian forests. Comparison of the results of the two techniques showed that for the groups of beech and chestnut-leaved oak, the ANFIS technique performed better and that the modeling techniques have a deep relationship with the nature of the tree species.

1. Introduction

Estimating the diameter increment of forests is one of the most important relationships in forest management and planning [1,2]. Therefore, choosing an accurate method to determine this relationship in the forests has great importance. A great variety of growth model systems exist, and they are usually grouped in different levels of resolution such as stand, diameter class, and individual tree [3].

Growth and yield models predict the dynamics of forests, including future forest growth and products that enable the study of alternative management options [4,5]. The diameter increment of trees is affected by internal factors such as physiology, species, age, and genetic characteristics and external factors such as climatic conditions, ground slope, soil type, and competition from nearby trees [6].

The revolution in computing technology has affected growth and yield modeling as it has all fields of science [7]. Computationally intensive methods have been applied in various studies, showing that the use of machine learning methods has many advantages [8,9,10,11]. Machine learning methods have been used to detect complex nonlinear relationships and work with qualitative data. The utility of these methods has been well demonstrated in many studies [12,13,14,15,16,17]. ANNs form a subset of machine learning that is used for a wide range of problem-solving tasks, including in optimization, prediction, diagnosis, and control [18]. Another method of machine learning that has many applications in various fields, including forestry, is the ANFIS method, and the results of various studies in this field are satisfactory [15,19]. This method has been widely used in recent years because it is not limited to the assumptions of experimental models and has the ability to solve problems related to nonlinear relationships [13]. The ANFIS is a fuzzy system with a parallel structure, and neural network learning algorithms are used to adjust the parameters of the fuzzy system [20].

Among the available methods in machine learning, ANNs are the most common method for modeling at the individual-tree level [21], while usage of other techniques such as the adaptive neuro-fuzzy inference system (ANFIS) and fuzzy logic are in their early stages [18]. The ANFIS as a nonlinear model was proposed by Jang [22] for the first time. In its structure, there are two different parametric groups that include premise and consequence. The training of the ANFIS involves determining these parameters using an optimization algorithm [23]. The ANFIS is a two-in-one package, as if you are utilizing a neural network and fuzzy logic advantages at the same time [22].

Different studies have been conducted in various fields of forest science using ANNs and ANFIS methods. Reis et al. [24] investigated the mortality and survival of individual-trees in selectively harvested forests using ANNs in Amazon rainforests. Training and testing of ANNs were performed for modeling, classification, mortality, and survival, using different input variables. Reis et al. [24] argued that the overall efficiencies of the classification in the training and testing phases were above 89% and 90%, respectively. In another study conducted by Vieira et al. [18], the predicted growth of diameter at breast height (DBH) and height of eucalyptus trees were determined using both techniques (ANN and ANFIS). These two methods have high predictive accuracy for DBH growth and tree height. Reis et al. [25], in their research, predicted a diameter growth model using an ANN with a correlation above 99% and root mean square error (RMSE) below 11%. They argued that the ANN could be used effectively to help in the management of rainforests and ensure the environmental and economic sustainability of the forest. Ashraf et al. [26], using machine learning, developed growth and yield models (basal area and volume increment) in forests. The coefficients of determination (R²) of the model using field data were 0.38 and 0.60 for the increments of the basal area and volume, respectively. Other studies have also highlighted the successful use of ANNs in forest model development, including predicting the growth and mortality of trees (e.g., [27,28,29]). In the Hyrcanian uneven-aged and mixed forests, Bayat et al. [15] estimated forest tree height using an ANN and ANFIS and concluded that these methods have good ability for tree height estimation.

Hyrcanian forests are of great importance as the only commercial forests in Iran that have temperate broad-leaved species. Their approximate area is 1.85 million hectares and about 15% of Iran’s forests. This region is of special importance with about 44% of the vascular plants in Iran. There are about 500 species of plants native to Iran in these forests. Therefore, sustainable and effective planning and management in this region have always been a major concern. Accurate estimation of the characteristics of these forests, including diameter increment, has great importance for planners [30].

However, the application of ANNs and the ANFIS to predict the diameter increment in uneven-aged forests has not been previously conducted in the world. Therefore, by providing diameter increment models in mixed and uneven-aged Caspian forests using these two techniques, we are presenting here an account of novel approaches to build a foundation for future studies.

Therefore, in this study, some specific issues were investigated: (1) estimation of the diameter increment in mixed and uneven-aged forests using machine learning (the main objective), (2) the potential use of nonparametric models including ANNs and the ANFIS for the estimation of the diameter increment, (3) a comparison between the ANFIS and ANN in estimating tree increment, and finally (4) optimization of ANN and ANFIS methods when used in predicting tree diameter increment.

2. Materials and Methods

2.1. Study Area and Data

Kheyroud educational research forest is located in watershed 45 in northern forests of Iran seven kilometers east of the city of Noshahr in Mazandaran Province. The forest plan categorizes this area into eight districts including the Gorazbon district. The Gorazbon district, with an area of 934.24 hectares, is the third largest district of Kheyroud forest (Figure 1). This district is located in the longitudinal range 51° 36′ 30″ to 51° 39′ 30″ and latitude 36° 32′ 0″ to 36° 34′ 0″. The general orientation of this district is north–south, and the minimum and maximum altitudes are 840 and 1350 m above sea level. Brown soil is typical in this district.

The data used in modeling were obtained from two sampling periods (2003 and 2012) with fixed-area sample plots. In 2003, using a 150 × 200 m grid, 256 permanent circular sample plots with area of 0.1ha were present in the Gorazbon district in the form of a systematic sample. Each sample plot was recorded with a GNSS receiver. The accuracy of the registration of the forest monitoring permanent sample plots’ center depends on several factors, including the density of the canopy of the study forest [31]. Inside the sample plots, the diameter at breast height (DBH) of all living trees with a diameter of more than 7.5 cm was measured using caliper, and their values were recorded in inventory forms and in one-centimeter classes. These operations were repeated after nine years. In total, seven tree species consisting of Fagus orientalis, Carpinus betulus, Quercus castaneifolia, Alnus subcordat, Acer cappadocicum, Acer velutinum, and Tilia begonifolia were tallied in sample plots. For the sake of our study and simplification, we put these species into four broad groups such as beech, hornbeam, chestnut-leaved oak, and other species [32].

2.2. Determination of the Model Inputs

The data we recorded were used as input to calculate models (

Table 1. General characteristics of the main tree groups.

Groups of Trees	Year	Maximum	Average	Standard Deviation
Beech	2003	178	39.1	29.1
	2012	180	41.8	36.1
Chestnut-leaved oak	2003	108	22.7	13.3
	2012	115	25.9	18.3
Hornbeam	2003	133	23.9	17.5
	2012	136	26.1	19.5
Other species	2003	186	34.3	25.6
	2012	188	37.5	26

). These inputs consisted of natural logarithm of DBH (cm), square of DBH (cm²), basal area of the largest tree (BAL) (m²/ha), Shannon–Weiner’s species diversity index (H_S) (Equation (1)) [33], size diversity index (H_d) (Equation (2)), and the mean of basal area (BA) in sample plot (m²) (

Table 2. Forest stand status in terms of different parameters.

Variable	Minimum	Maximum	Standard Deviation
The number of trees per hectare	20	1220	241
DBH (cm)	7.5	188	24.7
Average of BA in sample plot (m2)	0.02	0.633	0.1
BAL (m2/ha)	0	52	8.4
Hd in sample plot	0	2.468	0.314
Hs in sample plot	0	1.8	0.663

Wherein: diameters at breast height (DBH), basal area (BA), basal area of the largest tree (BAL) (m²/ha), size diversity index (H_d), and Shannon–Weiner’s species diversity index (H_S).

).

$$H_s=-{\displaystyle \sum _{i=1}^{n_s}\frac{B_i}{B}\ln \left(\frac{B_i}{B_{}}\right)}$$

(1)

$$H_d=-{\displaystyle \sum _{j=1}^{n_d}\frac{Bj}{B}\ln \left(\frac{Bj}{B}\right)}$$

(2)

2.3. Model Development

Almost 80% of the studied trees come from two species (beech and hornbeam), and nearly 20% fall into chestnut-leaved oak and other species groups in the study site. Consequently, the species were divided into four groups of beech (2215 trees), hornbeam (3131 trees), chestnut-leaved oak (340 trees), and other species (856 trees), and a separate model was developed for each group.

Machine learning has various algorithms including artificial neural network (ANN), genetic algorithm, fuzzy logic, and neuro-fuzzy [30,34,35]. For predicting the diameter increment of individual trees, two techniques of adaptive neuro-fuzzy inference system (ANFIS) and multi-layer perceptron (MLP) neural network were designed and implemented in MATLAB software version 2016. Prior to delving into the details of these techniques, it is better to describe the data stratification method. In general, the data used in the ANN model were divided into two subsets of training and testing data. Training section is larger and usually involves 70 to 80% of the data and is used for training the learning patterns to the network [36,37]. The test section is not utilized in the network learning process. Rather, it is used for comparing different models in terms of network power (goodness of fit) and generalization [38]. In the present study, the k-fold cross validation method was utilized. This method has been applied by Kozak and Kozak [39] and others. In our study, all observations after random distribution were divided into four-fold (i.e., k = 4), and for each species group every data set contained 25% of observations. Therefore, four data sets were created with different structures (

Table 3. The formation of data set.

Data Set	Validation	Training
A	1	2	3	4
B	2	4	3	1
C	3	3	2	1
D	4	4	2	1

). The data sets were utilized to construct four models for both MLP and ANFIS techniques. The three sets of data were used for model development (training) and one data set for error estimation in each model. For each species group, 4 models with different formation of the data sets were obtained. The final step included averaging the calculations of outputs of the four models.

2.3.1. Multilayer Perceptron Artificial Neural Network (MLP)

The typical architecture of MLP consists of three layers: an input, hidden, and output layer [40]. The architecture of these networks includes the number of layers, the number of neurons, the transfer function of each neuron, and how the layers are interconnected depending on the type of problem [26]. The number of independent variables of the problem determines the number of network inputs, and the number of dependent variables specifies the neurons’ number of output layer. The main issue in the design of the ANN is the determination of the number of suitable neurons in the hidden layer(s), which is achieved with the trial-and-error method. A single neuron is not very useful for modeling real and complex systems, whereas several neurons have a high power when they are linked with one or more layers. In general, the larger number of layers enables the system to represent more complexities. However, on the other hand, a high number of layers decreases the accuracy of the prediction model and may prevent the network from converging. In such a case, it has been proven that MLPs with one hidden layer, a sigmoid function (in the hidden layer), and a linear function in the output layer are able to provide an approximation of all the functions, provided that there are enough neurons in the hidden layer [41,42]. Hence, in this study, the number of neurons in each layer was determined as follows:

(i)

The number of neurons in input layer: The determination of inputs number for each species group was obtained by trial and error. Thus, the combination with the lowest mean squared error (MSE) values was selected as the final combination of the model. The tested combinations are shown in Table 4. Therefore, according to each species group, the suitable and normalized inputs were entered into the network. Then, these inputs, using hyperbolic tangent (tansig) function [26], were passed from hidden layer to the output layer.

(ii)

The number of neurons in hidden layer: In MLP, the nodes of hidden layer are a structural sensitive parameter. This means that very low nodes could lead to ill training and high nodes lead to over-fitting [40]. In this study, the structures of MLP with 1 to 30 neurons in the hidden layer were analyzed and their prediction accuracy was evaluated.

(iii)

The number of neurons in output layer: According to the aim of this study (diameter increment estimation) and number of parameters, the neuron’s number in the output layer is equal to 1. The outputs of hidden layer are eventually passed through the output layer with a linear transfer function (purelin) to provide the final output of the model.

Table 4. Combinations of model inputs.

No.	Model Inputs
1	Natural logarithm of DBH (cm), square of DBH (m), BAL
2	Natural logarithm of DBH (cm), square of DBH (m), BAL, the average of BA at sample plot
3	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hs
4	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hd
5	Natural logarithm of DBH (cm), square of DBH (m), BAL, average BA at sample plot, Hs
6	Natural logarithm of DBH (cm), square of DBH (m), BAL, average BA at sample plot, Hd
7	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hs, Hd
8	Natural logarithm of DBH (cm), square of DBH (m), BAL, average BA at sample plot, Hs, Hd

Table 4. Combinations of model inputs.

No.	Model Inputs
1	Natural logarithm of DBH (cm), square of DBH (m), BAL
2	Natural logarithm of DBH (cm), square of DBH (m), BAL, the average of BA at sample plot
3	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hs
4	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hd
5	Natural logarithm of DBH (cm), square of DBH (m), BAL, average BA at sample plot, Hs
6	Natural logarithm of DBH (cm), square of DBH (m), BAL, average BA at sample plot, Hd
7	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hs, Hd
8	Natural logarithm of DBH (cm), square of DBH (m), BAL, average BA at sample plot, Hs, Hd

The outputs of the model were compared with the target values, and the differences were calculated as suggested by earlier studies [43,44]. The process was continued until network’s weights were optimized. The prediction capability and performance of the developed ANN were determined using mathematical and statistical methods. In order to assess the network performance and terminate the training, the MSE was utilized as a measure of the accuracy of the ANN results. The Levenberg–Marquardt (LM) algorithm depicted the best performance in fitting functions. Therefore, this LM algorithm was used in this study.

2.3.2. Neuro-Fuzzy Inference System (ANFIS)

ANFIS was also used to predict diameter increment of individual trees. ANFIS structure consists of five layers: input nodes, rule nodes, average nodes, consequent nodes, and output nodes [45,46]. Each node has a function with an adjustable or fixed component [46]. The proper structure of the ANFIS is selected according to the input data, membership degree, the rules, and the membership functions of the input and output. The input values could become closer to the actual ones, if the membership degree components are adjusted during the training stage.

For our study, we used the same stratification method as was conducted for MLP (Table 4) during modeling and testing stages. Just as with MLP (Table 4), the inputs of ANFIS technique were determined by trial and error. ANFIS uses neural network learning algorithms and fuzzy logic to design nonlinear mapping between input and output space [22]. It also has the advantage of adaptation capability and rapid learning capacity [47]. In addition, it can adjust the conversion of human intelligence to fuzzy systems. The learning rule is based on the back-propagation algorithm with the aim of minimizing the MSE between the network output and the actual output. There are three types of systems included in the fuzzy inference system (FIS): Mamdani system, Sugeno system, and Tsukamoto system. In most cases, Mamdani system is used; however, the Sugeno system performs better than the former in computing and has a definite (crisp) output [46]. Therefore, for our study, modeling was performed by ANFIS technique using Sugeno system. The grid partitioning method was used to determine the degree and type of desirable membership functions such as triangular-shaped, trapezoidal-shaped, bell-shaped, Gaussian, pi-shaped, Gaussian 2, and sigmoid through trial and error. Further, the hybrid algorithm, which is a combination of the back-propagation algorithm and the least-squares error method, was used for training and adapting with the FIS.

2.3.3. Model Evaluation

Model evaluation is an important part of the modeling process and should be considered at different stages of modeling, including design, adjustment, and implementation. The evaluation of the model consists of two main stages, which include verification (qualitative) and validation (quantitative) [48]. In this study, five criteria including coefficient of determination (R²) (Equation (3)), root mean square error (RMSE) (Equation (4)), Akaike information criterion (AIC) (Equation (5)), Bayesian information criterion (BIC) (Equation (6)), and correlation coefficient (r) were used. The coefficient of determination (R²) is the proportion of the variance in the dependent variable that can be explained by the independent variables [49,50]. The root mean square error (RMSE) is the standard deviation of the residuals, which is the difference between the surveyed data and the fitted model [51,52].

$$R^2=\frac{{({\displaystyle {\sum }_{i=1}^N(O_i-\overline{O})(P_i-\overline{P})})}^2}{{\displaystyle {\sum }_{i=1}^N(O_i-\overline{O}}{)}^2{\displaystyle {\sum }_{i=1}^N{(P_i-\overline{P})}^2}}$$

(3)

$$\mathrm{RMSE}={{\displaystyle \left[{{\displaystyle N}}^{-1}{\displaystyle {\sum }_{i=1}^N{{\displaystyle \left({{\displaystyle P}}_i-{{\displaystyle O}}_i\right)}}^2}\right]}}^{0.5}$$

(4)

Akaike [53] proposed the Akaike information criterion (AIC) to estimate the expected Kullback–Leibler distance between the model generating the data and the fitted candidate model. AIC can be expressed as:

$$\mathrm{AIC}=N.\ln (SSE)+2n_{par}$$

(5)

In order to overcome the over-estimation problem of AIC, Schwarz [54] proposed the Bayesian information criterion (BIC). BIC was derived using a Bayesian perspective, and it selects the model with the maximum posterior probability for a given prior probability. BIC can be expressed as:

$$\mathrm{BIC}=N.\ln (SSE)+n_{par}\ln (N)$$

(6)

where O, p, $\overline{O}$, $\overline{P}$, N, n_par, and SSE are observed values, estimated values, the average of observed values, the average of estimated values, total number of observations, number of inputs used in the model, and sum of the squared error (SSE), respectively. The best model has the highest value of R² and the lowest of three other criteria (RMSE, AIC, and BIC).

3. Results

For the development of the ANN and ANFIS models, all networks examined in the MATLAB software environment (version 2016) were designed and implemented, and the results were derived separately for the groups of beech, hornbeam, chestnut-leaved oak, and other species. During the simulation process, the settings of the networks’ parameters were selected for both techniques (ANN and ANFIS) based on the minimum MSE via a trial-and-error procedure. The results indicated that for all groups of species in all data sets, the training of the network is more successful and has a higher accuracy (R² and Pearson correlation coefficient (r)) than the network validation. In addition, the correlation values of the two techniques (ANN and ANFIS) for all species groups in all data sets were significant at a 0.01 level. The best combination of input variables was the same for both techniques (

Table 5. The best combination of variables in relation to each species group.

Species Group	Input Variables
Beech	Natural logarithm of DBH (cm), square of DBH (m), BAL, the average of BA, Hd
Hornbeam	Natural logarithm of DBH (cm), square of DBH (m), BAL
Chestnut-leaved oak	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hs
Other species	Natural logarithm of DBH (cm), square of DBH (m), BAL, Hs

).

In the ANFIS technique, the k-fold strategy was used just as with the MLP. The results of this technique are presented separately for the four groups. The optimal ANFIS models were obtained by implementing a hybrid learning algorithm. The type and number of membership functions (in the form of trial and error) were determined using MSE criteria. The number of membership functions for all models was considered to be two, because the network was over-fitted with the increase in the number of membership functions. Moreover, triangular membership function was obtained as the best membership function. We used a linear transfer function in the output layer. The grid partitioning function was used to create the network architecture and the FIS. The number of iterations for optimization of the model was 100.

3.1. Beech Group

In the MLP model, data set D showed the best results in terms of the correlation and coefficient of determination (R²). However, by comparing the four data sets in terms of the three evaluation criteria of RMSE, AIC, and BIC, a different result is obtained, so that the lowest RMSE of training and validation were seen in groups C and B, respectively. Further, group C had the least value of AIC and BIC criteria and showed better performance (

Table 6. Evaluation criteria of MLP and ANFIS for beech species.

Data Set	Type	MLP					ANFIS
		r	R2	RMSE	BIC	AIC	r	R2	RMSE	BIC	AIC
A	Training	0.47	0.22	1.61	-	-	0.52	0.27	1.58	-	-
	Validation	0.42	0.18	1.68	-	-	0.35	0.12	1.71	-	-
	Total	0.44	0.20	1.64	19,345	19,278	0.48	0.23	1.62	19,217	19,189
B	Training	0.45	0.204	1.66	-	-	0.54	0.29	1.51	-	-
	Validation	0.41	0.17	1.59	-	-	0.35	0.12	1.84	-	-
	Total	0.45	0.204	1.64	19,332	19,265	0.49	0.24	1.6	19,150	19,122
C	Training	0.48	0.23	1.60	-	-	0.53	0.28	1.58	-	-
	Validation	0.40	0.16	1.83	-	-	0.40	0.16	1.62	-	-
	Total	0.47	0.22	1.63	19,331	19,234	0.50	0.25	1.59	19,136	19,108
D	Training	0.47	0.22	1.631	-	-	0.51	0.26	1.57	-	-
	Validation	0.42	0.18	1.619	-	-	0.41	0.17	1.7	-	-
	Total	0.46	0.21	1.632	19,309	19,242	0.49	0.24	1.6	19,191	19,163
Final model		0.48	0.23	1.613	19,256	19,189	0.51	0.26	1.57	19,117	19,089

).

In the ANFIS model, given the R² of both validation and training, data set C provided the best results. Moreover, by comparing four data sets in terms of the three evaluation criteria of RMSE, AIC, and BIC, data set C had the lowest value and presented the best result and performance (Table 6).

The error rates in the correlation and RMSE (MLP model) for the training and validation (four data sets) from the beech species are presented in Figure 2. The results of the models showed that for different data sets, correlation and RMSE values varied with increasing number of nodes, so that by increasing the number of nodes, the correlation value of the training of all data sets increases, while, for the validation data set, the trend is approximately unchanged. This trend is also seen in the case of RMSE criteria. The lowest correlation with two nodes and the highest correlation with 30 nodes were observed (Figure 2). Finally, according to the evaluation criteria of training and validation, it can be said that the best MLP network had 21 hidden nodes (i.e., a 5-21-1 structure).

3.2. Hornbeam Group

In the MLP model of the hornbeam group, data set B had the highest R2, data set A and C had the lowest RMSE, and data set B had the lowest AIC and BIC. Therefore, data set B generally gives the best result. In the ANFIS model, data set B had the highest R², and C had the lowest RMSE, AIC, and BIC. Thus, for the hornbeam species, data set B generally provides the best outcome (

Table 7. Evaluation criteria of MLP and ANFIS for hornbeam species.

Data Set	Type	MLP					ANFIS
		r	R2	RMSE	BIC	AIC	r	R2	RMSE	BIC	AIC
A	Training	0.35	0.12	1.45	-	-	0.32	0.10	1.57	-	-
	Validation	0.33	0.11	1.38	-	-	0.3	0.09	1.38	-	-
	Total	0.35	0.12	1.42	27,474	27,426	0.32	0.10	1.52	27,871	27,853
B	Training	0.37	0.14	1.42	-	-	0.35	0.12	1.57	-	-
	Validation	0.33	0.11	1.46	-	-	0.3	0.09	1.46	-	-
	Total	0.36	0.13	1.43	27,468	27,419	0.33	0.11	1.54	27,948	27,930
C	Training	0.34	0.11	1.41	-	-	0.3	0.09	1.51	-	-
	Validation	0.26	0.07	1.53	-	-	0.26	0.07	1.51	-	-
	Total	0.32	0.10	1.42	27,473	27,424	0.28	0.08	1.51	27,814	27,796
D	Training	0.33	0.11	1.42	-	-	0.33	0.11	1.58	-	-
	Validation	0.23	0.05	1.44	-	-	0.28	0.08	1.43	-	-
	Total	0.32	0.10	1.44	27,521	27,473	0.32	0.10	1.54	27,941	27,923
Final model		0.36	0.13	1.42	27,465	27,407	0.32	0.10	1.49	27,762	27,740

).

The results of the MLP model for hornbeam species indicated that by increasing the number of nodes, the error rates in the correlation of training and validation data were incremental and decreasing, respectively. However, in the case of RMSE criteria, this trend is reversed. The lowest correlation was observed with two nodes and the highest correlation occurred with 29 nodes (Figure 3). According to the description provided, the best MLP model for hornbeam species consists of 29 hidden nodes (i.e., a 5-29-1 structure).

3.3. Chestnut-Leaved Oak Group

In the MLP model of this species, the difference between correlation and RMSE of training and validation is more than other studied species. The R² of the training was approximately twice the R² of the validation data. In both MLP and ANFIS models, in terms of correlation, RMSE, BIC, and AIC, it can be said that the data sets of A, B, B, and B provide the best results, respectively. Additionally, data set B also had the best performance in both models (

Table 8. Evaluation criteria of MLP and ANFIS for chestnut-leaved oak species.

Data Set	Type	MLP					ANFIS
		r	R2	RMSE	BIC	AIC	r	R2	RMSE	BIC	AIC
A	Training	0.65	0.42	1.58	-	-	0.68	0.46	1.33	-	-
	Validation	0.34	0.12	2.24	-	-	0.52	0.27	1.7	-	-
	Total	0.59	0.35	1.76	2413	2375	0.63	0.40	1.52	2264	2260
B	Training	0.60	0.36	1.63	-	-	0.65	0.42	1.42	-	-
	Validation	0.29	0.08	1.56	-	-	0.57	0.33	1.62	-	-
	Total	0.56	0.31	1.59	2345	2306	0.61	0.37	1.48	2259	2255
C	Training	0.29	0.08	1.93	-	-	0.62	0.38	1.42	-	-
	Validation	0.31	0.10	1.82	-	-	0.45	0.20	1.75	-	-
	Total	0.29	0.08	1.84	2442	2403	0.56	0.32	1.51	2275	2271
D	Training	0.51	0.26	1.65	-	-	0.62	0.38	1.41	-	-
	Validation	0.18	0.03	2.23	-	-	0.45	0.20	1.74	-	-
	Total	0.47	0.22	1.94	2479	2440	0.57	0.33	1.5	2254	2250
Final model		0.53	0.28	1.55	2328	2289	0.62	0.39	1.47	2251	2247

).

The results of the changes in correlation and RMSE of chestnut-leaved oak species are slightly more complicated and do not show a tangible trend. It can only be said that by increasing the number of nodes, the value of correlation for training and validation data increases and decreases, respectively. The lowest correlation was observed with two nodes and the highest correlation with 21 nodes (Figure 4). The best MLP network of chestnut-leaved oak species has 21 hidden nodes (i.e., a 4-21-1 structure).

3.4. Other Species Group

In the MLP model of this group, the data set B had the highest R² in training. For the other three evaluation criteria (RMSE, AIC, and BIC), the best performance was for the data set D. Finally, we can say that data set D provides the best performance because the validation provides better R² than the training (

Table 9. Evaluation criteria of MLP and ANFIS for other species group.

Data Set	Type	MLP					ANFIS
		r	R2	RMSE	BIC	AIC	r	R2	RMSE	BIC	AIC
A	Training	0.52	0.27	1.53	-	-	0.50	0.25	1.52	-	-
	Validation	0.27	0.07	1.87	-	-	0.33	0.11	1.64	-	-
	Total	0.50	0.25	1.55	6555	6536	0.46	0.21	1.54	6557	6538
B	Training	0.65	0.42	1.42	-	-	0.53	0.28	1.43	-	-
	Validation	0.46	0.22	1.51	-	-	0.32	0.10	1.85	-	-
	Total	0.63	0.40	1.48	6474	6455	0.47	0.22	1.54	6507	6488
C	Training	0.59	0.35	1.48	-	-	0.54	0.29	1.57	-	-
	Validation	0.40	0.16	1.52	-	-	0.35	0.12	1.49	-	-
	Total	0.56	0.31	1.48	6482	6463	0.5	0.25	1.51	6499	6480
D	Training	0.59	0.35	1.48	-	-	0.53	0.28	1.49	-	-
	Validation	0.52	0.27	1.49	-	-	0.37	0.14	1.64	-	-
	Total	0.58	0.34	1.48	6466	6466	0.48	0.23	1.53	6511	6492
Final model		0.57	0.32	1.44	6432	6413	0.49	0.24	1.5	6498	6479

). In the ANFIS model, data set C presented the best result in terms of all evaluation criteria (Table 9).

The results of the MLP model for the other species group showed that by increasing nodes (up to 25 nodes), the correlation showed an increasing trend. This means that by increasing the nodes, the error rate of the model decreases. The lowest correlation was found with one node (Figure 5). Therefore, the best MLP model of the other species group consists of 25 hidden nodes (i.e., a 4-25-1 structure).

4. Discussion

Since the habitat conditions of the Hyrcanian forests change from west to east and in elevation, it is very difficult to develop increment models. Therefore, identifying the appropriate predictors (inputs of the model) on increment and incorporating them into increment models is key to making predictions. For this purpose, the best input variables in the diameter increment models of four groups of beech, hornbeam, chestnut-leaved oak, and other species were identified using MLP (a class of feedforward ANN) and ANFIS techniques. The results were then compared.

Our results indicated that in relation to each group of species, different parameters could be significant, such as a previous study concluded [55,56]. For example, in relation to beech species, factors such as DBH and natural logarithm of DBH, BAL, H_d, and average of BA were significant. These variables are associated with hornbeam (DBH, natural logarithm of DBH, and BAL) and chestnut-leaved oak, and other species (DBH and natural logarithm of DBH, BAL, and H_s) showed a different pattern. The effectiveness of the two indices, size diversity and species diversity of Shannon–Wiener, on the increment rate were consistent with the results of Liang [57]. In the research for oak forests, an individual-tree model was presented, and it was concluded that the indices affecting the diameter increment include tree size, basal area, tree diameter, volume inventory, and site index [58]. The results of Lhotka and Loewenstein [59] on the development of individual-tree diameter growth models for the Missouri forest stands in the United States showed that the BAL parameter was effective for all species, while the patterns of remaining predictors were different for other species. The results of Lhotka and Loewenstein [59] were in line with the findings of the present study, except in one case. Unlike their findings, in our study the species composition (species diversity index) was not effective for all species groups. In general, the BAL variable has been used as the most important competition variable in previous growth and increment modeling studies (e.g., [60,61,62,63,64]). In this study, we incorporated all the measured independent variables, i.e., DBH, BAL, H_S, H_d, BA, number of trees per hectare, slope, aspect, and altitude into the modeling process. However, only DBH, BAL, BA, H_S, and H_d were significant and considered in the final model. Other variables were removed from the final model because they were insignificant in terms of correlation because the study area was relatively consistent in some respects (e.g., slope). If a larger area is selected for study, these variables may become meaningful [30,40]. However, variables of size diversity and species diversity have not been closely evaluated, except for few studies (e.g., [1,57]). Given the fact that the Hyrcanian forests have high biodiversity (species diversity) and high size diversity [30], consideration of these two indices in new studies is most important.

Comparison of the MLP and ANFIS techniques showed that the latter performed better for of the beech and chestnut-leaved oak groups. On the other hand, MLP showed promising results for the hornbeam and other species groups. However, in the study of Vieira et al. [18], no significant difference was observed for diameter growth estimation between the ANFIS and ANN techniques. It seems that with increasing complexity of the relationship, the ANN technique brings better results [65]. Beech and chestnut-leaved oak clusters occupy more distinctive site conditions than the other two groups in the Hyrcanian forests. For instance, the beech class is mainly situated in the northern side [6], while the chestnut-leaved oak group is situated in the southern and southwestern directions. The hornbeam group appeared in almost all sample plots (high distribution). Therefore, the determination of its diameter increment was more challenging and complicated, and the coefficient of determination (R²) of the model showed the lowest value compared to the other tree groups. This means that the inputs of the model did not properly account for the variation in the diameter increment. It is recommended that hybrid models be fitted for this group. The low R² values can be attributed to a large number of heterogeneous data used in this study. However, the values were statistically significant and allowed for identifying major growth driving factors. For future works, other factors (e.g., soil properties) from a large-scale area might be incorporated into the models for measuring, and perhaps improving, their performance.

Using the k-fold strategy in evaluation showed that different combinations of data give different results. Therefore, the use of this strategy and the acquisition of the final model of the four-folds (where, k = 4) average provides the most realistic result. We found this strategy useful, as it has been considered promising in other recent studies (e.g., [26,61,63]). The correlation coefficient of models, despite not having high values, were significant at a 0.01 level. There could be several reasons for such a pattern of correlation coefficients. One possibility could be the existence of trees with the same diameter and significant differences in diameter increment values [66], which resulted from forest competition and succession [67]. The second reason might be the autocorrelation of data used in increment models due to both time and spatial components [26]. Although permanent sample plots are the most popular method for estimating forest growth and yield and have been used in numerous studies (e.g., [60,61,62,63]), they are not immune from the autocorrelation. The autocorrelation created within the sample plots was greater than the autocorrelation between the sample plots, because the sample plots are separated by more distance than the trees within a sample plot [26].

Since previous studies carried out using ANNs and the ANFIS are different from our study in terms of diverse aspects such as forest type (even-aged or uneven-aged) and research goal (diameter increment, diameter growth, basal area growth survival, etc.), comparisons between our results and those of most prior studies are not possible.

5. Conclusions

The aim of this study was to provide insight into the utilization of two machine learning methods (MLP and ANFIS) for developing a diameter increment model for the Hyrcanian forests. The findings of this study would for the first time provide general information on the number of optimal hidden neurons (in ANNs) and the type and number of membership functions (in ANFIS) for diameter increment modeling.

RMSE and R² of the MLP and ANFIS models were estimated for the four groups of beech ((1.61 and 0.23) and (1.57 and 0.26)), hornbeam ((1.42 and 0.13) and (1.49 and 0.10)), chestnut-leaved oak ((1.55 and 0.28) and (1.47 and 0.39)), and other species ((1.44 and 0.32) and (1.5 and 0.24)), respectively. Despite the low coefficient of determination, the correlation test in both techniques was significant at a 0.01 level for all four groups. In general, the ANFIS worked better when the data from the oak and beech groups were used, probably because these groups represent the dense-covered areas with trees. On the other hand, the ANN performed better with the hornbeam and other species groups that represent a wide but less-covered area with trees. Our study also provides information on optimizing and adjusting the parameters necessary for the application of machine learning in developing prediction models for the estimation of the diameter increment in the Hyrcanian forests. We found that modeling techniques have a deep relationship with the nature of the tree species. The results provide guidance for future studies in the same area (Hyrcanian forest) or elsewhere in the diverse forests of the world.