Introduction Top

The average height of dominant or co-dominant trees at a specified index age, defined as the site index (SI), is used as an indicator of forest site quality in both pure and mixed stands (Carmean, 1972). The prediction of SI is required for efficient forest management planning since it is an imperative starting point for forest-level plans and silvicultural strategies. In particular, forest managers require site quality predictions to assess the amount of wood volume to be produced for a specific species within a time period (Skovsgaard & Vanclay, 2008). In the forestry literature, many statistical techniques have been proposed to evaluate site qualities with SI predictions. Importantly, Bailey & Clutter (1974) pioneered the algebraic difference approach (ADA) in the forestry literature for formalizing the base-age invariance feature in SI models. ADA is independent of the initial selected base age and produces polymorphic SI curves with single asymptotes for a given SI class. Thereafter, Cieszewski & Bailey (2000) introduced the generalized algebraic difference approach (GADA) as a generalization of ADA. It can produce polymorphic site index curves with multiple asymptotes.

Many nonlinear growth models such as Chapman-Richards, Bertanlaffy, Gompertz, Soloboda, & Korf models have been intensively used for developing SI models and constructing SI curves. These growth models have provided desirable characteristics including polymorphism, multiple horizontal asymptotes, sigmoid trend, base-age invariance, and inflection point presence in SI predictions. In spite of the attractive performance of these growth models in predicting height-age relations, some data obtained from permanent sample plots over time or stem analyses of dominant trees cause them to violate independence among observations due to highly correlated data features in spatial and time terms. Ignoring basic statistical assumptions often results in biased standard errors in parameter predictions with these growth models (Grégoire et al., 1995). Historically, Clutter (1961) and Lappi & Bailey (1988) are the earliest pioneering studies that have reported that both spatial and temporal correlation structures are possible with data from permanent sample plots or stem analysis of individual trees. For dealing with statistical issues arising from temporal and spatial correlations in height-age relations of dominant trees, autoregressive models (ARMs) such as the first autoregressive structure AR (1), the first-order continuous autoregressive error structure (CAR (1)), and the first-order autoregressive and moving average structures (ARMA (1, 1)) have often been implemented (Weiskittel et al., 2009; Wang et al., 2011). Similarly, linear and nonlinear mixed-effects models (NLMEs) (Lappi & Bailey, 1988; Wang et al., 2007; 2008) and dummy variable models (DVM) (Martín-Benito et al., 2008; Wang et al., 2008) have been applied to minimize impacts of spatial autocorrelations resulted from variation between sample trees or plots. Many studies have shown an improvement ranging from moderate to high upgrading fitting statistics accompanied by NLMEs or ARMs. However, the above-mentioned statistical issues limit the use of regression models for SI predictions and require new methods for satisfactory predictions of site quality.

Figure 1. Location of the study areas.

Artificial neural networks (ANNs), as an application and a type of artificial intelligence system, differ from many innovative prediction techniques since they need no statistical assumptions. Furthermore, ANNs can provide efficient predictions in complex and highly nonlinear data from dynamic forest systems. Therefore, forest modelers have considered ANNs as an effective tool to obtain sound and competent predictions of some individual trees as well as stand and forest attributes. Many studies have compared ANNs and regression models (e.g. NLMEs) using fit statistics (e.g. bias and root mean square error (RMSE)). ANNs need to be tested in the context of desirable characteristics of dominant height growth, including sigmoidal growth, polymorphism, and multiple asymptotes for a thorough assessment of its prediction performance.

The objectives of this study were to: (1) evaluate the performance of ANNs to predict dominant height growth obtained from stem analysis; (2) compare ANNs and regression models (i.e., NLME, DVM, and ARM) in terms of fit statistics; and (3) evaluate the ability of ANNs in the context of desirable characteristics of dominant height growth.

Material and methods Top

Study site

This study was conducted in mixed Oriental spruce (Picea orientalis L.) and Scotch pine (Pinus sylvestris L.) stands in the northeast of Türkiye. These mixed stands are located at Torul State Forest Enterprise (SFE) and Maçka SFE, in the Trabzon Forest District Directorate (FDD); Tirebolu SFE and Espiye SFE in the Giresun FDD (Fig. 1). The altitudes and slopes of the studied locations vary from 650 to 1450 m and 15 to 50%, respectively. These areas are characterized geomorphologically as high mountainous land with moderate to steep slopes. The annual mean temperature is 10.9 ºC and mean daily maximum temperature ranges from 11.8 ºC to 28.6 ºC; mean daily minimum temperature ranges from -3.8 ºC to 17.4 ºC, respectively. This climate is a typical Black Sea climate, characterized by a mild winter and a cool summer. The mean annual rainfall varies from 900 to 1350 mm.

Data

The sample plots were carefully investigated for traces of intensive silvicultural treatments or clear-cutting within naturally regenerated and regularly stocked stands (60–90% of tree layer cover). The stands inventoried, that were thinned moderately and had no evidence of historical damage such as fire or storms, were sampled in this study. The Oriental spruce and Scotch pine mixed stands were selected to have uniform stratification of these species when both of them had been in the upper stratum, such that there were site trees of both species in the plot. A total of 161 temporary circular plots were sampled subjectively by Ercanli (2010) to represent the existing range of stand ages, density indexes and site qualities. The sampling plot sizes ranged from 600 to 1200 m² according to stand density to ensure a minimum of 30 trees per sample plot. In each sampling stand, one dominant Oriental spruce tree in accordance with 100 trees of the greatest height per hectare was sampled for stem analysis. The area where the study was carried out is under the management and administration of the Turkish National Forestry Directorate, but an authorization for only one felled tree was enabled in these plots. Therefore, the study was limited to only one felled tree at each sample point among the 100 trees of the greatest height per hectare.

These trees were felled, leaving stumps of a height of 0.30 m, and cross-sectional discs were cut at 3, 4, or 5 m intervals above ground surface to avoid distorting economic appraisal for tree stems in harvest applications of Turkish forestry. The total height of these felled trees was measured to the nearest 0.01 m. The number of rings was counted at each cross-sectioned point, and then these values were converted to stump age, which can be considered equal to tree age. As cross section lengths do not coincide with periodic height growth, it was necessary to adjust height and age data from stem analysis. Fabbio et al. (1994) presented the ISSA method that uses second differences in ring counts to have smoother height/age curves. Height at each age of the felled trees was calculated by the ISSA method based on height and age data at cross-sectional cuts. The data set in this study consisted of a total of 161 stem analyses with 873 height/age pairs. The mean, standard deviations, minimum, and maximum values for the data used are presented in Table 1. Also, profile plots of stem analyses for the felled trees are presented in Fig. 2.

Figure 2. Line plots of 873 dominant height data obtained from 161 stem analysis trees.

Methods

Base nonlinear model

In the literature on forests, a variety of statistical growth functions have been employed to describe the dominant height-age relationship and to predict site index of forest stands. Besides the fitting abilities of these statistical models, the model used to obtain site index predictions has some desirable growth properties, including: polymorphism, base-age invariance, sigmoid growth pattern with an inflection point, passing a zero point, multiple asymptotes with function of site index (increasing with increasing site index), parallel asymptote at older ages and a non-decreasing curve over time (Corral-Rivas et al., 2004; Diéguez-Aranda et al., 2006; Martín-Benito et al., 2008). In the present study, the GADA expansion of Chapman-Richards growth model (Chapman, 1961; Richars, 1959) was selected as the base function due to its better predictive results compare to other functions in our preliminary analyses, in which some preliminary analyses were performed based on both fitting criteria and the desirable growth properties to determine the predictive activation function of choice from these statistical functions. The Chapman-Richards growth model

was solved by a and assuming

ResultsTop

From the results for these ANN models, a network architecture with a tan- or log-sigmoid transfer function for the input layer improved model accuracy for both FF and CC network architectures (Table 2). In contrast, selecting the purelin transfer function for the input layer considerably reduces model accuracy and efficiency. For instance, the RMSE value for CC-tansig-purelin (cascade correlation learning algorithms with tan-sig function between input layer and hidden layer and Pure-Lin function between hidden layer and output layer, abbreviated as tan-pure) with 46 neurons was 1.14, while it was 2.40 for CC-purelin-tansig with 21 neurons. Also, FI for CC-tansig-purelin with 46 neurons was 0.97, while it was 0.86 for CC-purelin-tansig with 21 neurons. In summary, the selection of this transfer function, purelin, for the input layer for dominant height growth was inappropriate, and the selection of the tansig or logsig transfer function for the input layer was more suitable for our data set. Note that our network architecture with a logsig transfer function for the output layer for both ANN architectures did not converge.

In general, ANNs resulted in moderately worse predictive ability with higher AAE, max. AAE, RMSE, RMSE%, AIC, BIC, Bias, and Bias%, and the lowest r and FI compared to given regression techniques. While the RMSE differences between ARMA (1,1) and ANNs with tan-or-log-sigmoid transfer functions for the input layer were negligible (0.2 m), the differences between the ANNs and NLME-f were significant (> 0.5 m), implying a 40% reduction in RMSE (Table 2). These results suggested that spatial autocorrelations were more influential than temporal autocorrelations in the dominant height growth trends of Oriental spruce located in these studied stands.

Figs. 3, 4, and 5 depict graphical analyses of RMSE% and lagged-residuals, as well as observed vs. predicted values. Fig. 3 shows that NLME-f had a smaller RMSE% across all ages than other models. CC-tansig-tansig with 53 neurons and the ARMA (1,1) displayed a very similar pattern in terms of RMSE% and had poorer estimates across ages over 100 years. For lower ages than 100 years, the base model and NLME-f-AR1 yielded worse predictions. While NLME-f, CC-tansig-tansig with 53 neurons and ARMA (1, 1) behaved consistently above 50 years, the base model and NLME-f-AR1 showed a steady decrease in terms of RMSE%. All models produced the largest RMSE% for young trees below 50 years. While RMSE% for NLME-f, CC-tansig-tansig with 53 neurons, and ARMA (1, 1) was less than 10% for ages older than 50 years, it was less than 10% and RMSE% for the base model and NLME-f-AR1 over 90 years, indicating that NLME-f, and ARMA (1, 1) models for over 50 years, and NLME-f-AR1 over 90 years can be used reliably in prediction of dominant height growth in the studied stands.

Figure 3. Line plots of regression techniques, the base nonlinear model, ARMA (1, 1), NLME-f, NLME-f-AR1 and cascade-forward network architecture with tan-sigmoid transfer function for both input and output layers (CF-tan-tan-53) in terms of RMSE.

The one-lagged residuals of all models, except for CC-tansig-tansig with 53 neurons, suggested no temporal autocorrelations, and were identically distributed around zero (Fig. 4). CC-tansig-tansig with 53 neurons showed a slight temporal autocorrelation, but it was not clear evidence of temporal autocorrelation. Correcting for temporal autocorrelations in residuals of the base model did not improve prediction accuracy. While the RMSE for the base model was 1.32 m, it was 1.09 m for ARMA (1, 1). Likewise, Bias was -0.01 and -0.02 for the base model and ARMA (1, 1), respectively. Graphics of observed vs. predicted dominant heights showed NLME-f to be better than other models (Fig. 5). Patterns of the base model, NLME-f-AR1, and CF-tan-tan with 53 neurons were highly similar, but were slightly different from ARMA (1, 1).

Figure 4. The one lagged residuals of regression techniques, the base nonlinear model, ARMA (1, 1), NLME-f, NLME-f-AR1 and cascade-forward network architecture with tan-sigmoid transfer function for both input and output layers (CF-tan-tan-53).

Some graphs of dominant height growth for five site index values (i.e. 5, 10, 15, 20, 25, and 30 m at base age 100) for given regression techniques, feed-forward and cascade-forward network architectures are shown in Fig. 6, Fig. 7 and Fig. 8, respectively. All regression techniques provided expected behaviors of dominant height growth (e.g. polymorphism, multiple asymptotes, and sigmoid pattern). It is interesting that while FF-purelin-tansig with 28 neurons and CC-purelin-tansig with 21 neurons performed worse than other ANN architectures in terms of given fit statistics, these network architectures provided most of the desirable properties of dominant height growth. In other words, ANN architectures that were fitted well, failed to assure the desirable properties of dominant height growth in this study. In detail, dominant height growth for the regression techniques, FF-purelin-tansig with 28 neurons and CC-purelin-tansig with 21 neurons showed a sigmoidal pattern (known as the “S curve”). Site index curves for these models are polymorphic for varying site indexes, and follow multiple asymptotes at older ages (i.e. beyond 150 years). While the regression techniques had strong asymptotic behaviors for five site index curves, that is a desirable property for dominant height growth and provide the biological interpretation, FF-purelin-tansig with 28 neurons and CC-purelin-tansig with 21 neurons, had fairly weak asymptotic behaviors. For the predictions from these ANN models, it is obvious that h is not 0 at age 0 and extrapolating beyond the measurements range indicate the SI curves go to a common asymptote. In the same way, whereas dominant heights for the regression techniques were equal to zero at age zero, those for all ANNs were not at age zero. Also, whereas dominant heights for the regression techniques were equal to site index values at base age 100, those for all ANNs architecture were different from site index values at a given base age. All these findings on site index curves show that Chapman-Richards site index model, which is an n-degree polynomial regression model to a certain extent, is more effective than the different ANN model structures in obtaining site index predictions.

DiscussionTop

In this study, dominant height data were obtained from stem analysis and these data, which were obtained from the same tree, were serially correlated, and this feature may cause a larger error variance as reported by many studies (Lappi & Bailey, 1988; Wang et al., 2007; Martín-Benito et al., 2008; Wang et al., 2008). We implemented ARMA (1, 1) to remove the effects of the covariance structure on predictions and achieved a significant improvement in predictions (Table 2). Also, mixed-effects models have been recommended to account for variation among observations driven by unobserved influences. This unobservable variability is taken into account by adding subject-specific parameters at the tree or plot level (even region level) to a given linear or nonlinear model (Wang et al., 2008). In context of site index, Calegario et al. (2005) stated that subject-specific parameters are needed to achieve desirable site curves that vary with a site quality. In this study, the nonlinear models with a single random parameter generally outperformed those with an autoregressive structure (Table 2). NLME with two and three random parameters performed generally worse than those with one random parameter. Conversely, Wang et al. (2011) reported that NLME with two or more random parameters performed worse relative to those with a single random parameter, which probably resulted from predicting two or more random parameters using a single height measurement. Wang et al. (2008) also reported that in the case of using multiple measurements, NLME with two random parameters outperformed the dummy variable method since NLME considers both the subject-specific parameter and the error parameter.

Fortin et al. (2007) demonstrated that combining correlated error terms that appeared in repeated measurements with random effects provided the smallest fit statistics in predicting merchantable basal area. However, they also stated that considering correlated error terms only may provide more accurate estimates than combining correlations in residuals and random effects in the case of less variability among subjects (e.g., trees). Likewise, Martín-Benito et al. (2008) used dummy variables as subject-specific responses and could not achieve higher accuracy because of less variability resulting from compensating effects among regions. Conversely, in some experiments, as evidenced by our study, the correlated errors seemed to be negligible, whereas the variability among subject trees was significant (Wang et al., 2011). One possible reason is that correlations between closer measurements in time were insignificant (<1%) (Nothdurft et al., 2006). Additionally, since data sets were collected from different ecological regions in this study, adding tree-specific parameters as random effects to the base nonlinear model enhanced the prediction accuracy significantly (Gregorie, 1987; Wang et al., 2008).

Most of the dominant height models developed resulted in erratic estimates at especially younger ages below 50 years. This could be due to changes in the status of dominant trees over time, which could be influenced by factors such as a lack of knowledge about the selection of dominant trees at these ages, as well as temperature and rainfall (García, 2005; Martn-Benito et al., 2008; Weiskittel et al., 2009). Therefore, forest managers should treat young stands with care before applying this model.

Figure 5. Scatter plots relating to observations vs. predictions using regression techniques, the base nonlinear model, ARMA (1, 1), NLME-f, NLME-f-AR1 and cascade-forward network architecture with tan-sigmoid transfer function for both input and output layers (CF-tan-tan-53).

The discussions associated with ANN have focused on prediction accuracy, but few studies have focused on biological reality. Since several factors affect the growth of a tree (or a forest stand), relations between variables and the significance of the variables are needed to biologically interpret parameter estimates of the models. The black-box structure into the internal workings of an ANN used for training, on the other hand, complicates biological interpretations of model parameters, inferences on the population of interest, and comparisons of differences among subjects (e.g., trees in this study). Yee et al. (1993), Aertsen et al. (2010) and Jevšenak & Levanič (2016) evaluated some modeling techniques such as ANN and generalized additive model (GAM) in terms of statistical, ecological interpretability, and user-friendliness, and suggested that other modeling techniques (e.g. GAM) may be preferred over ANN to predict site index. Also, Yee et al. (1993) reported the lack of biological interpretation as one of the disadvantages of ANN, and stated that assessment of an ecological model should be based on both the fit statistics and the ability to describe growth curves. Brosofske et al. (2014) commented that if a researcher’s aim is to improve the accuracy of predictions, ANN may be a good choice in the case of a large dataset, but not a small dataset. In a case study (Castaño-Santamaría et al., 2013), ANN resulted in poor predictions in an uneven-aged forest due to very high height-diameter diversity, which resulted from different site qualities that obscured the training of ANN. They also explained that including diameter at breast height variance to each diameter class resulted in ANN being able to provide similar results with a local-mixed model.

The shape of the fitted curve is highly important for a reliable dominant height model (Guan & Gertner, 1991). However, the overfitting issue, as reported in previous studies by Aertsen et al. (2010), Özçelik et al. (2013) and Jevšenak & Levanič (2016) turns out to be decisive in assuring desirable properties of dominant height growth. To overcome this problem, an alternative way is to use the Levenberg-Marquardt training algorithm, as recommended by Jevšenak & Levanič (2016). However, in our study, we failed to avoid overfitting in the validation data set. The required properties of dominant height growth were reported by Goelz & Burk (1992). Accordingly, our results showed that most ANN architectures could not provide well-suited fitted curves in terms of dominant height growth (Figs. 4 and 5). Only FF-purelin-tansig with 28 neurons and CC-purelin-tansig with 21 neurons enabled a few required properties, including polymorphism and sigmoidal pattern, suggesting that the ANN needs improving to meet biological presumptions in terms of individual dominant height growth.

Figure 6. Site index curves for five site index values for Oriental spruce stands using regression techniques, the base nonlinear model, ARMA (1, 1), NLME-f, and NLME-f-AR1.

Figure 7. Site index curves for five site index values for Oriental spruce stands using feed-forward network architectures, FF-log-tan (52), FF-tan-tan (100), FF-tan-pure (87), FF-log-pure (75), FF-pure-pure (77), and FF-pure-tan (28).

Figure 8. Site index curves for five site index values for Oriental spruce stands using cascade-forward network architectures, CF-log-tan (24), CF-tan-tan (53), CF-tan-pure (46), CF-log-pure (28), CF-pure-pure (90), and CF-pure-tan (21).

ConclusionsTop

NLMEs with a single random parameter were better at fitting the stem analysis data when compared with NLMEs with two or three parameters. The worst estimates among the NLMEs were obtained by those with three parameters having approximately four times greater RMSE when compared to the base nonlinear model. In the data obtained by the stem analysis, the spatial autocorrelation was more evident than the temporal autocorrelation. These regression techniques had a homogeneous residual variance and had greater correlations with the measurements. The desirable properties for dominant height growth were provided by these regression techniques, such as multiple asymptotes, base age invariance, and sigmoidal pattern.

ANNs had worse estimates than the regression techniques and did not fulfill the desirable properties for dominant height growth. ANNs suggested very different curves depending on the transfer functions. This study showed that while selecting tan or log sigmoid transfer functions for the input layer was more appropriate in terms of the fit statistics, the purelin transfer function was more suitable in terms of ensuring the desirable properties. Selecting a purelin transfer function for the input layer and a tan-sigmoid transfer function for the output layer provided some desirable properties, including polymorphism and sigmoidal pattern. In the output layer, choosing a tan-sigmoid other than the purelin function, which presents a linear trend, will be more likely to provide a nonlinear trend compatible with polymorphism and sigmoidal pattern. Also, this study showed that determining an accurate transfer function was essential for the biological assumptions but was not sufficient. Considering the learning rate factor may increase the model efficiency and provide opportunities for ensuring the assumptions of interest.