Estimating Crimean juniper tree height using nonlinear regression and artificial neural network models

Forest Ecology and Management 306 (2013) 52–60

Contents lists available at SciVerse ScienceDirect

Forest Ecology and Management journal homepage: www.elsevier.com/locate/foreco

Full length article

Estimating Crimean juniper tree height using nonlinear regression and artificial neural network models Ramazan Özçelik a,⇑, Maria J. Diamantopoulou b, Felipe Crecente-Campo c, Unal Eler a a

Faculty of Forestry, Süleyman Demirel University, East Campus, 32260 Isparta, Turkey Faculty of Forestry and Natural Environment, Aristotle University of Thessaloniki, GR-54124 Thessaloniki, Greece c Departamento de Ingeniería Agroforestal, Universidad de Santiago de Compostela, Escuela Politécnica Superior, R Benigno Ledo, Campus universitario, 27002 Lugo, Spain b

a r t i c l e

i n f o

Article history: Received 22 April 2013 Received in revised form 3 June 2013 Accepted 10 June 2013

Keywords: Mixed-effects model Generalized h–d model Back-propagation neural network model Tree height estimation

a b s t r a c t Artificial neural network models offer a number of advantages including the ability to implicitly detect complex nonlinear relationships between input and output variables, which are very helpful in tree height modeling. Back-propagation artificial neural network models were produced for individual-tree height estimation and the results were compared with the most used tree height estimation methods. Height diameter (h–d) measurements of 1163 Crimean juniper trees in 63 sample plots from southwestern region of Turkey were used. A calibrated basic h–d mixed model, a generalized h–d model and backpropagation artificial neural network h–d models were constructed and compared. When the variability of the h–d relationship from stand to stand can be incorporated into the model, then both mixed-effects nonlinear regression and back-propagation neural network modeling approaches can produce accurate results, reducing the root mean squared error by more than 20% as compared to a basic nonlinear regression model. The use of a generalized h–d model also showed reliable results (reduction of 13% in root mean squared error as compared to a nonlinear regression model). The back-propagation artificial neural network model seems a reliable alternative to the other methods examined possessing the best generalization ability. Further, from a practical point of view it has the advantage that no height measurements are needed for its implementation. On the contrary prior information is required for the mixed-effects model calibration which is a limiting factor according to its use. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Currently, Juniper species are economically and ecologically ones of the most common and important forest tree species in Turkey, where they constitute large forest areas and occupy approximately 450,000 ha, with a current standing volume of approximately 8 million cubic meters (GDF, 2006). These species are more common in the South and North sides of the Taurus Mountains, a mountain complex in southern Turkey, which divides the Mediterranean coastal region from the central Anatolian Plateau. Six species grow naturally in Turkey, among which Crimean juniper (Juniperus excels Bieb) and Foetid odur juniper (Juniperus foetidissima Willd.) form pure stands and are a valuable timber resource. The other species are less common. These species in the Taurus Mountains occur generally between 500 and 2500 m elevations. Especially, Crimean juniper forests have a great economic potential and have a key role in providing important benefits and environmental services such as protection of soil and water resources, conservation of biological diversity. ⇑ Corresponding author. Tel.: +90 246 2113956; fax: +90 246 2371810. E-mail address: [email protected] (R. Özçelik). 0378-1127/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.foreco.2013.06.009

In recent years, Turkey has adopted the principles of multipurpose and ecologically based forest management. Therefore, the General Directorate of Forests (GDF) needs to develop and evaluate growth and yield prediction models for management of forest resources. Furthermore, height–diameter models are needed to better understand the nature of various relationships that characterize, differentiate, and influence the development of forest ecosystems (Peng et al., 2001). However, the available information about height-diameter relationships concerning these species is very limited. More accurate and versatile information about important characteristics of forest resources is necessary for evaluating the numerous management and utilization alternatives for timber resources. Height–diameter (h–d) models are very useful for yield estimation (Curtis, 1967; Parresol, 1992), site index and dominant height estimation (Curtis, 1967; Calama and Montero, 2004), stand structural analysis (Spies and Cohen, 1992; Morrison et al., 1992; Gadow et al., 2001), damage appraisal and stand stability (Parresol, 1992; Vospernik et al., 2010), stand growth dynamics (Curtis, 1967; Burkhart and Strub, 1974; Wykoff et al., 1982), and product recovery and carbon budgeting models (Newton and Amponsah, 2007).

53

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

Most h–d models have been applied to pure even-aged stands or plantations (e.g. Soares and Tomé, 2002; López-Sánchez et al., 2003; Diéguez-Aranda et al., 2005). Despite the homogeneous characteristics of this type of forest, a single h–d model is not usually adequate for all possible situations, since the h–d relationship varies from stand to stand, and is not even constant within the same stand over time (Curtis, 1967). The most widely used method for minimizing this level of variance is to fit a local h–d equation (h as a function of only d and related parameters) for each plot and measurement occasion, which is the most accurate model that can be developed for a specific stand. The main problem of this approach is that it requires a high sampling effort. Alternative approaches which reduce this effort are generalized h–d models, mixed-effects h–d models and, more recently, the use of artificial neural network models (ANNs). Generalized h–d models use d and stand-specific variables as regressors, accounting for the differences in the h–d relationship within stands and with time. Several generalized and region-wide equations have been developed for many tree species all over the world (e.g., Curtis, 1967; Larsen and Hann, 1987; López-Sánchez et al., 2003; Temesgen and Gadow, 2004). Mixed-effects models allow for both population-averaged and cluster (or subject)-specific models. The first considers only fixedparameters, common to the population, while the second considers both fixed- and random-parameters, common to each subject. The inclusion of random parameters, specific for each stand, allows modeling the variability of the h–d relationship among different locations and time, after defining a common fixed functional structure (Lindstrom and Bates, 1990). If a prediction for a new stand is required and prior information (a small sample of trees measured for h and d) is available, the h–d curve can be calibrated to obtain a stand-specific response. Many studies have used mixed-models to develop h–d relationships (e.g., Lappi, 1997; Calama and Montero, 2004; Trincado et al., 2007; Sharma and Parton, 2007). Artificial neural network models (ANNs) have been applied successfully in the field of forest modeling. Specifically, ANNs have been used for internal log defect classification (Schmoldt et al., 1997), for the prediction of diameter distribution (Leduc et al., 2001), probability of juvenile trees and tree species per unit area (Hasenauer and Kinderman, 2002), forest attributes (Corne et al., 2004), bark volume of standing trees (Diamantopoulou, 2005), dielectric properties of wood (Avramidis et al., 2006), inside-bark diameter and heartwood diameter (Leite et al., 2011), and as a modeling technique for site index prediction of homogenous stands (Aertsen et al., 2010), among others. While ANNs have been applied for the prediction of tree volume (Nelson et al., 2009; Diamantopoulou and Milios, 2010; Özçelik et al., 2010), none of the above studies have incorporated h–d relationships into ANN models. Further, according to the knowledge of the authors, no work has been reported in the literature that has conducted a comparative study between the most used tree height estimation methods and the back-propagation ANN modeling technique in total tree-height estimation. The objective of this work was to conduct a comparative study between several methods for obtaining height predictions. For this purpose, several widely used nonlinear growth functions were fitted using nonlinear mixed-effects model techniques; several generalized h–d models were fitted by ordinary nonlinear least squares, and artificial neural network models were developed using sample plots of Crimean juniper located in the Southern of Turkey. The fitted models represent new management tools for forest managers and could be implemented in whole-stand or individual tree level models for Crimean juniper stands in the region.

2. Materials and methods 2.1. Study area and data description The data used in this study was obtained from 63 sample plots established in pure (more than 90% of juniper trees), natural, evenaged stands of Crimean juniper in south and southwestern region of Turkey. The plots were located throughout the area of study, and were subjectively selected to represent the existing range of ages, stand densities and sites for the species in Turkey. Plot size ranged from 373 to 2400 m2, depending on stand density, in order to achieve a minimum of 30 trees per plot. 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 arithmetically averaged to obtain diameter at breast height (d, cm). Total tree height (h, m) was measured to the nearest 0.5 m, using a Blume– Leiss hypsometer, in a randomized sample of approximately onethird of the trees in the plot. The trees were subjectively selected to ensure a representative distribution by diameter and height classes. In order to assess the developed h–d model performance, the available data set was randomly divided into a construction data set (85% of the plots) and an evaluation data set (the remaining 15% of the plots). In the model fitting process, the evaluation data was not used. Summary statistics for the two mentioned tree variables, plus some stand variables, are provided in Table 1. 2.2. Model development and comparison 2.2.1. Non-linear mixed-effects modeling Six basic h–d models, selected from previous studies (Peng et al., 2001; Fang and Bailey, 1998; Huang et al., 2000), were fitted to the dataset. Some modifications to the models were also tested (i.e. adding 1.3 to the equation). These models include the most flexible equations for h–d relationships (Zhang, 1997; Peng et al., 2001), the Bertalanffy–Richards (Bertalanffy, 1949, 1957; Richards, 1959), Weibull (Yang et al., 1978) and Schnute (Schnute, 1981) models: a2

Bertalanffy—Richards hij ¼ 1:3 þ a0  ð1  expða1  dij ÞÞ a

Weibull hij ¼ 1:3 þ a0  ð1  expða1  dij2 ÞÞ Exponential hij ¼ 1:3 þ a0  exp



a1 dij þ a2

ð1Þ ð2Þ

 ð3Þ a

Korf—Lundgvist hij ¼ 1:3 þ a0  expða1  dij 2 Þ

ð4Þ

Gompertz hij ¼ 1:3 þ a0  expða1  expða2  dij ÞÞ

ð5Þ

 1 1  expða2  dij Þ a0 a0 a0 a0 Schnute hij ¼ 1:3 þ ða1  1:3 Þ 1  expða2  100Þ

ð6Þ

where hij and dij are, respectively, the total height and diameter at breast height of the jth tree in the ith plot, and ai are the parameters to be estimated. These models were firstly fitted by the ordinary nonlinear least squares (ONLS) method with the MODEL procedure of SAS/ETSÒ (SAS Institute Inc (2002)). Statistical and graphical analyses were used to compare the performance of the models, as explained later. A nonlinear mixed-effects modeling framework was then used to fit the above-mentioned models (e.g. Lappi, 1997; Calama and Montero, 2004; Castado-Dorado et al., 2006) in an attempt to find the best possible fit to the database. Basically the parameter vector of the nonlinear model can be defined (Pinheiro and Bates, 1998) as:

54

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

Table 1 Summary statistics for the fitting and the evaluation data sets. Variable

d (cm) h (m) t(years) SI (m) BA (m2 ha1) N (trees ha1)

Construction data set (967 trees in 54 plots)

Evaluation data set (196 trees in 9 plots)

Mean

Min.

Max.

S.D.

Mean

Min.

Max.

S.D.

24.1 8.6 95 10.9 31.6 900

4.0 2.0 44 7.0 16.7 267

110.0 19.0 162.0 15.0 60.2 2905

13.9 3.2 37 1.9 10.8 550

21.7 7.9 85 10.7 30.1 880

5.0 3.0 48 7.0 21.0 373

83.0 19.0 158 14.0 42.0 1650

12.7 3.3 34 1.9 7.5 325

d, diameter at breast height (1.3 m above ground level); h, total tree height; t, individual tree age; SI, site index (calculated with the model developed by Eler, 1986); BA, stand basal area; N, number of trees per hectare.

Ui ¼ Ai k þ Bi bi

ð7Þ

where k is the p  1 vector of fixed population parameters (where p is the number of fixed parameters in the model), bi is the q  1 vector of random effects associated with the ith plot (where q is the number of random parameters in the model), and Ai and Bi are design matrices of size r  p and r  q (where r is the total number of parameters in the model) for the fixed and random effects specific to each plot. The SAS macro NLINMIX (Litell et al., 2006) was used to fit the models. Maximization of the marginal likelihood function was achieved using the best linear unbiased predictor (BLUP) estimation (Beal and Sheiner, 1982). All the possible linear combinations of fixed and random effects in one, two and tree parameters were tested and compared with the criteria explained later. 2.2.2. Calibrated response An advantage of mixed effects models is that if a subsample of k tree heights is available, such data can be used to predict the random effects vector bi, with the following expression (Vonesh and Chinchilli, 1997):

^ D bZ b TðR bi þ Z b i DZ b T Þ1 ^ei b i i i

ð8Þ

b is a q  q unstructured (in this case) variance–covariance where D matrix for the among-plot variability, common to all plots and estib i is the k  k variance– mated in the general fitting of the model; R ei is the residual covariance matrix for within-plot variability; ^ vector k  1, the components of which are given by the difference between the observed height value for each tree included in the subsample, and the value predicted by the model including only b i is the k  q matrix of partial derivatives evalufixed effects; and Z ^ i . Once b ^ i is predicted, the value of the vector of heights, i.e. ated at b the calibrated response vector, can be calculated with the general expression of the parameter vector of the nonlinear model (Eq. (7)). If we wish to predict heights of a particular stand with no prior h–d observations, the fixed-effects typical response should be used, considering that the best linear unbiased predictor of bi is the null ^ i ¼ Ai ^ q  1 vector, i.e., U k. For the best mixed-effects model, the calibrated response was evaluated for different height sampling designs and sampling sizes within each plot, using the evaluation dataset (9 plots, 15% of the total number of plots used to fit the models) and using the conb The calibration alternastruction dataset for estimating ^ k and D. tives selected were: (i) (ii) (iii) (iv) (v)

Total height of 1–10 randomly selected trees per plot. Total height of 1–10 largest trees per the plot. Total height of 1–10 smallest trees per the plot. Total height of 1–10 medium-size trees per plot. Total height of 3, 6 and 9 trees per plot, in the largest, small^ was then calculated for each est and medium-size classes.b i of the 9 plots. The calibrated h–d model was then applied to

all the trees in each of these plots. The five alternatives were evaluated with the criteria explained later, and compared with the estimations obtained with ordinary nonlinear least squares (ONLS) in the individual fit of the selected model to each of the calibration plots. For the randomly selected trees, mean values of the statistics after 100 simulations were obtained. 2.2.3. Generalized h–d models In a first step, 25 generalized h–d models selected from previous studies (Krumland and Wensel, 1988; Tomé, 1989; López-Sánchez et al., 2003; Sharma and Zhang, 2004; Castado-Dorado et al., 2006; Sharma and Parton, 2007; Crecente-Campo et al., 2010) were fitted to the dataset using non-linear least squares. These models include generalizations of the most flexible equations for h–d relationships (Zhang, 1997; Peng et al., 2001). In a second step some generalizations of Eqs. (1)–(6) were tested in this study by relating the equations parameters to stand variables. For this purpose, and individual fitting of Eqs. (1)–(6) to each plot was carried out using ONLS. Graphical analyses were performed to know which stand variables were more related with local parameters and to observe the type of relationship (e.g., linear, allometric or exponential). All the developed models were evaluated with the criteria explained later. 2.2.4. Artificial neural network models Artificial neural network models (ANNs) are computational models inspired in the natural neurons. They have been developed as generalizations of mathematical models of human cognition or neural biology. A neural network is characterized by its pattern of connections between nodes (its architecture); its method of determining the weights of the connections (its training or learning); and its activation function (Fausett, 1994). ANNs can learn from experience. Because of its remarkable similarities with the information processing features of the human brain (nonlinearity, high parallelism, capability to generalize), this modeling approach has the ability to solve problems that are difficult to formalize, such as problems of biological nature. The fact that they require no assumptions about the form of a fitting function makes this approach a very interesting estimation tool. As a result, this modeling approach can liberate the modeler from reliance on parametric approximating functions that may fit the observed data. The back-propagation algorithm has been used as the central algorithm to the training of multilayer perceptrons, which is the layout of the most popular neural network. This algorithm has been thoroughly described (Rumelhart et al., 1986; Fausett, 1994; Haykin, 1994; Patterson, 1996; Gurney, 1999). In brief, the network is first initialized by setting up all its weights to be small random values. The feedforward face is being applied by each input pattern to receive input signal and broadcasts this signal to all

55

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

nodes in the hidden layer. Each hidden node sums its weighted input signals, applies its activation function to compute its output signal, and sends this signal to the output node where its output signal is computed. In the back-propagation face, the associated error (dk) is calculated. This error is used to change the weights by calculating its weight correction term that is later used to update the relative weights and sends dk to the hidden nodes. Each hidden node sums its input from the output node multiplied by the derivative of its activation function to calculate its error information term. Convergence is sometimes faster if a momentum term is added to the weight update formula (Fausett, 1994):

Dwjk ðt þ 1Þ ¼ adk zj þ lDwjk

ð9Þ

where wjk is the bias on output unit k, a is the learning rate, dk is the portion of error correction weight adjustment for wjk that is due to an error at output unit Yk, also the information about the error at unit Yk that is propagated back to the hidden units that feed into unit Yk, zj is the output activation of the hidden unit j and l is the momentum parameter. The new weights for training step (t + 1) are based on the weights at training step (t). Momentum also reduces the likelihood that the neural network will find weights that are a local minimum (Fausett, 1994). Another significant parameter that is used in the training procedure is the learning rate (a) that places an upper limit on the amount by which a weight can be changed. Consequently the learning rate influences the speed of the learning. The network architecture, the number of epochs, the learning rate and momentum factor values were selected after evaluating several combinations. Due to its important characteristics (continuous, differentiable, and computational efficient), the activation function that has been used for the back-propagation artificial neural network (BPANN) models, training is the logistic function, which is a sigmoid function with a range of (0, 1) and is defined (Swingler, 2001) as:

f ðxÞ ¼

1 1 þ es

ð10Þ

having derivative:

1

.h

es ð1 þ eðsÞ Þ P

2

i

ð11Þ

where s ¼ wi xi . That is the information the node transmits, in which wi are the weights and xi are the input values with s e [1, +1]. Two ANN models were fit: (1) the BPANN model which uses the diameter at breast height (d), and (2) the BPANN_V model which uses the diameter at breast height (d) and the diameter at breast height taking into account the d variance for each plot, introducing in the ANN model the transformed variable (d  dvar0.5)/dmean. It is worth noted that the BPANN model is not taking into account the d variance for each plot and it is comparable to any ordinary nonlinear regression model that could be fitted to our data, while the BPANN_V model is comparable to the mixed-effect and the generalized h–d models. The above models were built using the neural network toolbox in Matlab software (Demuth and Beale, 2000; Beale et al., 2012) developed by The Mathworks™. To avoid overfitting in ANN model building, the K = 10-fold cross-validation resampling technique was used. The advantage of this method is that all of the examples in the fitting data set are eventually used for both training and testing. Following this method, the fitting data set was randomly partitioned into training (the 70% of the fitting data) and testing (the remaining 30% of the fitting data) data sets. After the selection of the training parameters, the models were trained using the training data set. The inner model evaluation was achieved using the test data set. The above procedure was repeated for each K = 10 experiments, using K  1 = 9 folds for training and the remaining one for testing. The entire

procedure was repeated about 2000 times using the combination of different training parameters (number of epochs, momentum and learning rate values, and number of hidden nodes). Finally, training was stopped when the generalization ability reached a maximum. This ability was measured by the average error rate of the cross-validation examples. 2.2.5. Model performance criteria Four statistical criteria obtained from the residuals were examined to compare the performance of all the developed models: the root mean square error (RMSE), which indicates the accuracy of the estimates; the correlation coefficient (R), which is a measure used to determine the relative correlation and the goodness-of-fit between the estimated and the measured data; the mean bias (E), which indicates the precision of the model; and Akaike’s information criterion (AIC), which is an index used to select the best model from a group of candidate models (Akaike, 1974), summarized as follows:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn ^ 2 i¼1 ðyij  yij Þ RMSE ¼ n

ð12Þ

Pn  ^  i¼1 ½ðyij  yÞ  ðyij  yest Þ ffi R ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn P n 2 ^  2 i¼1 ðyij  yÞ  i¼1 ðyij  yest Þ Pn E¼

i¼1 ðyij

ð13Þ

^ij Þ y

ð14Þ

n

AIC ¼ n lnðRMSEÞ þ 2p

ð15Þ

^ij and y  are the measured, estimated and average values where yij, y est is the average value of of the dependent variable, respectively; y the estimated values; n is the total number of observations used for fitting the model; p is the number of model parameters to be estimated; and ln is the natural logarithm. Graphical analysis, specifically plots of residuals versus predicted values, and plots of observed versus predicted values, was used as a helpful tool to identify lack of fit. 3. Results The best ONLS model, according to the fitting statistics (Table 2) was the model by Gompertz (Eq. (5)), with all parameter estimates significant at the 0.01 level (Table 3). The best mixed-effects model, according to the fitting statistics (Table 2) was the model by Gompertz (Eq. (5)) with parameters a0 and a1 expanded to be mixed, resulting in:

hij ¼ 1:3 þ ða0 þ ui Þ  expðða1 þ mi Þ  expða2  dij ÞÞ

ð16Þ

where ui and vi are random parameters, associated to the ith plot, and the rest of variables are as previously defined. All parameter estimates for Eq. (16) were significant at the 0.01 level (Table 3). The mixed-effects model led to a decrease in RMSE of 25.3% compared with the ONLS model. The best generalized h–d model, according to the fitting statistics (Table 2) was a generalization of Eq. (1), i.e., the Bertalanffy– Richards model, including dominant diameter and dominant height as additional independent variables, resulting in: a

a

a4

hij ¼ 1:3 þ a0  h01  ð1  expða2  d03  dij ÞÞ

ð17Þ

where d0i and h0i are the dominant diameter and dominant height, respectively, of the ith plot, ai are the parameters of the model, and the rest of variables are as previously defined. All parameter estimates for Eq. (17) were significant at the 0.01 level (Table 3). The

56

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

Table 2 Goodness of fit statistics for: a basic h–d model fitted by ONLS, a basic h–d mixed-model, a generalized h–d model, a BPANN model and a BPANN_V model, for the fitting and the evaluation data sets (A calibrated model, using 3, 6 and 9 trees in the calibration process, is shown for the mixed model when using the evaluation dataset). Model

Type

RMSE

R

E

AIC

Fitting data set (967 trees in 54 plots) Gompertz (Eq. (5)) BPANN Gompertz (Eq. (16)) Bertalanffy–Richards (Eq. (17)) BPANN_V

ONLS ANN Mixed Generalized ANN

1.738 1.713 1.298 1.507 1.326

0.8316 0.8368 0.9108 0.8768 0.9062

0.00165 0.01467 0.02115 0.00778 0.06094

541 524 266 407 281

Evaluation data set (196 trees in 9 plots) Gompertz (Eq. (5)) BPANN Gompertz (Eq. (16)) Gompertz (Eq. (16)) Gompertz (Eq. (16)) Bertalanffy–Richards (Eq. (17)) BPANN_V

ONLS ANN Calibrated(3) Calibrated(6) Calibrated(9) Generalized ANN

1.883 1.728 1.434 1.371 1.277 1.386 1.343

0.8342 0.8598 0.9068 0.9172 0.9260 0.9112 0.9157

0.33926 0.72429 0.16095 0.17183 0.16175 0.00571 0.15256

130 107 85 76 62 74 57

RMSE, root mean square error; R, correlation coefficient; E, mean bias; AIC, Akaike’s information criterion; ANN, artificial neural network model.

Table 3 Parameter estimates for the basic h–d ONLS model, the h–d mixed model, and the h–d generalized model selected (standard errors in brackets). Parameter estimates

Model Eq. (5)

Eq. (16)

Eq. (17)

a0 a1 a2 a3 a4

13.40 (0.41) 2.114 (0.069) 0.05505 (0.00341)

11.17 (0.41) 1.894 (0.110) 0.06603 (0.00387)

2.080 (0.264) 0.7463 (0.0465) 0.2362 (0.0616) 0.5046 (0.0608) 0.9616 (0.0783)

r2

1.805 (0.087) 4.408 (1.112)

r2u r2m

0.2998 (0.0916)

ruv

0.7573 (0.2581)

2

2 u,

r , residual variance of the model; r variance of the random-effect ui, r2m , variance of the random-effect vi, ruv, covariance of the random-effects ui and vi.

generalized h–d model led to a decrease in RMSE of 13.3% compared with the basic ONLS model (Eq. (1)). Sensitivity analysis was performed for the BPANN_V model in order to obtain information about the relative importance of the two variables used in the network. The ratio between the error with an omission of each variable and the error obtained using all variables was calculated, and it was found to be much higher than one, for both variables, indicating that both variables contribute significantly to the estimation of the total height values. The best values for the learning rate and the momentum for both models are given in Table 4. The number of nodes in the hidden layer of the neural network that was used for various training models were 1, 2, 3, . . ., 7. The networks were trained for 1000 epochs for both models, as there was very negligible reduction of the root mean square error (RMSE) values after 1000 epochs (Table 4). One and two hidden nodes were chosen in the hidden layer for the BPANN and the BPANN_V trained models, respectively, because for the above values, the

average lowest residual error (RMSE) for the training set were obtained, for both models. The fit statistics for the constructed BPANN and BPANN_V models are given in Table 2. The BPANN model, which does not take into account the d variance of each plot, is comparable to the ONLS regression model (1.4% of reduction in RMSE compared with the ONLS model), while the BPANN_V model which takes the d variance of each plot into account is comparable to the mixed-effects regression model (23.7% of reduction in RMSE compared with the ONLS model). For the fitting data set, the fitness capability of the BPANN_V model as compared to the other examined regression models, was found to be almost as reliable as the mixed-effects model, which showed slightly better results. However, when using the evaluation data set, the BPANN_V model gave the most precise results, and more than 6 trees should be measured for calibration of the mixed-model to obtain better results than the BPANN_V model (Table 2) and the generalized h–d model (Eq. (17)). The measured total height values of the juniper trees were graphically compared to the corresponding values estimated by the Gompertz’s mixed-model (Eq. (16)), the Bertalanffy–Richards’ generalized h–d model (Eq. (17)) and the artificial neural network models (BPANN and BPANN_V) for the fitting and the evaluation data sets (Fig. 1). All models appear to be a valid solution for the tree height estimation, with Eq. (16) showing the best results as compared to the other approaches, following very closely by the BPANN_V model. The ONLS model is not compared as it showed the worse results over all the tested models. Residuals of Eqs. (16) and (17) and the BPANN and BPANN_V models showed no lack of fit, with a homogenous distribution around zero (Fig. 2). Specifically, the residuals of the mixed-effect model gave the smallest standard deviation (std. dev. = 1.271 m) for the fitting data set, while the BPANN_V model gave the smallest (std. dev. = 1.337 m) for the evaluation data set. Results for the calibration response of the mixed-effects model (Fig. 3) showed that the largest values of RMSE were obtained when applying the fixed effects response model, without predicting random parameters. In contrast, individual fitting to each plot

Table 4 Combinations of parameters that conducted to the best learning of the BPANN and BPANN_V models. Model

BPANN BPANN_V

Number of nodes Input layer

Hidden layer

Output layer

1 2

1 2

1 1

Number of Epochs

Learning rate

Momentum factor

1000 1000

0.10 0.10

0.30 0.30

57

12 8 4 fitting data

0

4

8

12

16

20

12 8 4 evaluation data

0

0

20 16 12 8 4 fitting data

0 0

4

8

12

4

8

12

16

16 12 8 4 fitting data

0

0

20

16

20 16 12 8 4 evaluation data

0

20

0

4

Measured tree height (m)

8

12

16

4

8

12

16

20 16 12 8 4 evaluation data

0

20

0

Measured tree height (m)

Measured tree height (m) BPANN model estimations (m)

BPANN model estimations (m)

Measured tree height (m)

20

BPANN_V model estimations (m)

0

16

20

Mixed model estimations (m)

16

20

20 16 12 8 4 fitting data

0 0

4

Measured tree height (m)

8

12

16

4

8

12

16

20

Measured tree height (m)

20

BPANN_V model estimations (m)

20

Mixed model estimations (m)

Generalized model estimations (m)

Generalized model estimations (m)

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

20 16 12 8 4 evaluation data

0 0

4

Measured tree height (m)

8

12

16

20

Measured tree height (m)

4 2 0 -2 -4 -6 -8

BPANN model residuals (m)

8

6

6 4 2 0 -2 -4 -6 -8

0

4

8

12

16

20

0

Generalized model estimations (m)

4

8

12

16

BPANN_V model residuals (m)

8

Mixed model residuals (m)

Generalized model residuals (m)

Fig. 1. Measured total height values of the juniper trees versus to the corresponding values estimated by mixed-models using the Gompertz’s model (Eq. (16)), by generalized h–d models using the Bertalanffy–Richards’ model (Eq. (17)) and by the BPANN and BPANN_V models, for the fitting and the evaluation data sets (for the evaluation dataset and the Gompertz’s mixed model, the residuals for the calibrated response using 6 trees of all the size classes are shown).

8 6 4 2 0 -2 -4 -6 -8

20

0

Mixed model estimations (m)

4

8

12

16

8 6 4 2 0 -2 -4 -6 -8

20

0

BPANN model estimations (m)

4

8

12

16

20

BPANN_V model estimations (m)

Fig. 2. Residuals versus predicted values of the juniper trees for the Gompertz’s mixed model (Eq. (16)), for the Bertalanffy–Richards’ generalized h–d model (Eq. (17)), and for the BPANN and BPANN_V models.

large

medium

small

all

random

fixed

0,95

0,50

0,90

0,25

ONLS

E

R

RMSE

2,00

1,50

0,00

0,85

1,00

-0,25

0,80 1

2

3

4

5

6

7

8

9

No. of trees for calibration

10

1

2

3

4

5

6

7

8

9

No. of trees for calibration

10

1

2

3

4

5

6

7

8

9

10

No. of trees for calibration

Fig. 3. Root Mean Squared Error (RMSE), correlation coefficient (R) and mean bias (E) for the fixed-effects model (fixed), the individual fit for each plot with ONLS (ONLS), and the calibrated model with different height sampling designs and sampling sizes within each plot, for calculating the random parameters (large = largest trees; medium = medium size trees, small = smallest trees; all = a mixture of large, medium and small trees; random = randomly selected trees).

58

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

with nonlinear least squares produced the smallest value, because it is the best possible fit of the function for each plot. The greatest reduction in RMSE in the calibration response was obtained using data from a mixture of trees selected from the smallest, largest and medium-size classes. Use of one tree in each class (3 trees) led to reduction in RMSE by 29%, and reduction of bias by 61% relative to the fixed-effects response. Use of the smallest trees in the plot is not a recommendable option since it showed the worse results. Finally, with randomly selected trees, medium-size trees, or the largest trees, the reduction in RMSE was smaller than in the case of selection of a mixture of trees from different classes.

4. Discussion A single h–d model is not usually adequate for all possible situations for a single species, because the h–d relationship usually varies within each location and also with time (Curtis, 1967). To minimize this level of variance, h–d relationships can be improved by taking into account stand variables that introduce the dynamics of each stand into the model (Curtis, 1967; Larsen and Hann, 1987; Temesgen and Gadow, 2004). Such models are the so-called generalized h–d models. Further, these relationships have been shown to be applicable to both even-aged (e.g. Castado-Dorado et al., 2006; Trincado et al., 2007) and uneven-aged stands (e.g. Sharma and Zhang, 2004; Sharma and Parton, 2007), and even to a mixture of both types of stands (e.g. Crecente-Campo et al., 2010). However, as Trincado et al. (2007) stated, the incorporation of additional predictor variables (i.e. stand variables) in a basic h–d model has a major effect on the ability of the fixed-effects model to explain between-individual variability, but certainly not on the mixed-effects model. Therefore, mixed-effects basic h–d models show a good alternative to generalized h–d models. It is also usually argued (Calama and Montero, 2004; Castado-Dorado et al., 2006) that mixed-models take into account the correlations between observations that belong to the same sampling unit, thus dealing with the dependence of error terms, which should be avoided in regression models. A major problem with mixed-effects models is that the predictive performance assessed from the residuals is not reasonable from an application point of view (Temesgen et al., 2008; Meng et al., 2009), since the model would never be calibrated with all the height measurements used in model fitting. So, the calibrated response should be taken into account when comparing the mixed model with other fitted models or fitting techniques, as it gives an idea of the error that would be obtained in practical applications. Among several widely used generalized h–d models tested in this study, Eq. (17) was found to be the most accurate, with the lowest AIC value. This model has been previously used by Crecente-Campo et al. (2010), who found it to be the most accurate for modeling the h–d relationship of blue gum eucalyptus (Eucalyptus globulus L.) stands in NW Spain. According to López-Sánchez et al. (2003), it is a medium sampling effort model, which needs measurements of diameter and of a sample of tree heights (those heights are necessary to calculate dominant height). Compared to the basic mixed-effects h–d model (Eq. (16)), it requires a similar sampling effort since the measurement of 6 heights in the stand allows calibration of Eq. (16), to obtain similar results than with Eq. (17). It should be noted that this is valid for a plot of approximately 600 m2, in which 6 trees should be measured to calculate dominant height. On plots more than 650 m2, at least 7 trees should be measured to calculate dominant height. This fact increase the sampling effort when the Bertalanffy–Richards generalized model (Eq. (17)) is used. Calibration of a mixed h–d model with information of only one tree (e.g., Trincado et al. 2007) has shown poor results in many studies (e.g., Calama and Montero, 2004; Castado-Dorado et al.,

2006; Crecente-Campo et al., 2010) including the results in this study, as well. The best results for the calibration of the basic h–d mixed-model (Eq. (16)) were obtained by selecting a mixture of trees from the smallest, largest and medium-size classes. The greater the number of measurements included in the subsample, the greater the decrease in RMSE (Fig. 3) and increase in R. However, a large sample is not often justifiable because of the increased cost of sampling (Castado-Dorado et al., 2006). The model bias follows a similar pattern, although the reduction is slower than for the previous statistics. Forest management requires reliable models that have the potential to describe the state and the dynamics of a forest. Following this certainty, effort has been made for the construction of growth models that could effectively specify the forest growth system; which is a complicated nonlinear system. In this system, as a comparable methodology to statistical modeling, artificial neural network models have received considerable attention. Specifically, according to their fitness capability that has been examined in this paper, they seem to be a reliable alternative approach for tree height prediction, which is a significant part of the forest growth system. Their most attractive advantage over non-linear regression models (generalized or mixed) is their ability to learn from the data without assumptions about the form of a fitting function, which is a very difficult procedure especially in the case of modeling the complex behavior of the tree bole growth and its dimensions (Diamantopoulou, 2010). The developed BPANN and BPANN_V models, can be easily applied because the only information that has to be provided to the model are the values of the diameter at breast height variable, and the most common descriptive statistics of this variable, i.e. its average value and its variance. Further, a BPANN model based on artificial networks associative ability is generally well fitted to missing or inaccurate data, once it has been developed (Diamantopoulou, 2010). However, because of the times of required repetitions and the need for the best combinations of the ANN attributes to be specified, big modeling effort is needed in order to optimize the accuracy and the generalization ability of an ANN model (Aertsen et al., 2010). From the results illustrated in this study, the estimation accuracy of the BPANN_V model was almost as high as that obtained with the mixed-effects models, and can be considered as a valid alternative for estimating the tree bole height for Crimean juniper in the study area. On the other hand, as mentioned a mixed model needs to be calibrated when used to take advantage of the random effects calculation (Trincado et al., 2007); for this purpose a sample of tree heights is necessary, thus increasing the sampling effort as regarding the BPANN_V model, which only needs diameter information. More than 6 trees should be measured for calibration of the mixed model to obtain as good results as that obtained from the BPANN_V model in this study.

5. Conclusions The comparative analysis between mixed-effects models, generalized models and back-propagation artificial network models conducted in this study, showed that when the variability of the h–d relationship from stand to stand can be incorporated into the model, both mixed-effects nonlinear regression and back-propagation neural networks modeling approaches are useful tools for tree bole height prediction in forest resource management. Specifically, the BPANN_V model was found to be superior to other models according to its generalization ability (Table 2). The Gompertz’s mixed-effects model (Eq. (16)) gave the most accurate results in terms of RMSE for the fitting data set. However,

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

according to the practical use of the examined models, there are limitations of each one of the above methodologies that should be taken into account, in order to select the most suitable model. The mixed-effects model (Eq. (16)) requires the measurement of at least 6 tree heights to be calibrated. The generalized h–d model (Eq. (17)) needs additional information, particularly the dominant height and the dominant diameter of each plot as input variables, meaning extra effort in the field. On the contrary, the BPANN_V model requires only the breast height diameter values as input variable. Finally, the BPANN_V model seems to be an effective alternative to regression analysis, providing an accurate tree height estimation approach which can supply the forest growth modeling system with better information, increasing the ability to properly manage forest ecosystems. Acknowledgements The authors thank the Turkish Forest Research Institute for its contribution to Field works. We would like to thank Dr. John R. Brooks and Dr. Harry V. Wiant Jr. for his valuable comments on the earlier drafts of the manuscript and for revising the English grammar of the text. We appreciate the valuable comments and constructive suggestions from the reviewers and the Editor. References Aertsen, W., Kint, V., Van Orshoven, J., Özkan, K., Muys, B., 2010. Comparison and ranking of different modelling techniques for prediction of site index in Mediterranean mountain forest. Ecological Modelling 221, 1119–1130. Akaike, H., 1974. A New Look at the Statistical Model Identification. IEEE Transactions on Automatic Control. AC-19, 716–723. Avramidis, S., Iliadis, L., Mansfield, S.D., 2006. Wood dielectric loss factor prediction with artificial neural networks. Wood Science and Technology 40, 563–574. Beal, S.L., Sheiner, L.B., 1982. Estimating population kinetics. CRC Critical Reviews in Biomedical Engineering 8, 195–222. Beale, M., Hagan, M., Demuth, H., 2012. Neural Network Toolbox™ User’s Guide. The MathWorks Inc., Noatic, MA, 404p.. Bertalanffy, L.V., 1949. Problems of organic growth. Nature 163, 156–158. Bertalanffy, L.V., 1957. Quantitative laws in metabolism and growth. Quarterly Review of Biology 32, 217–231. Burkhart, H.E., Strub, M.R., 1974. A model for simulation of planted loblolly pine stands. In: Growth models for tree and stand simulation. Royal College of Forestry, Stockholm, Sweden. pp. 128–135. Calama, R., Montero, G., 2004. Interregional nonlinear height–diameter model with random coefficients for stone pine in Spain. Canadian Journal of Forest Research 34, 150–163. Castado-Dorado, F., Diéguez-Aranda, U., Barrio Anta, M., Sánchez Rodríguez, M., Gadow, K.V., 2006. A generalized height–diameter model including random components for radiata pine plantations in northwestern Spain. Forest Ecology and Management 229, 202–213. Corne, S.A., Carver, S.J., Kunin, W.E., Lennon, J.J., Van Hees, W.W.S., 2004. Predicting forest attributes in southeast Alaska using artificial neural networks. Forest Science 50, 259–276. Crecente-Campo, F., Tomé, M., Soares, P., Diéguez-Aranda, U., 2010. A generalized nonlinear mixed-effects height-diameter model for Eucalyptus globulus L. in northwestern Spain. Forest Ecology and Management 259, 943–952. Curtis, R.O., 1967. Height–diameter and height–diameter–age equations for secondgrowth Douglas-fir. Forest Science 13, 365–375. Demuth, H., Beale, M., 2000. Neural network toolbox for use with Matlab: Computation, Visualization, Programming. User’s Guide. The MathWorks Inc., Natic, MA, 504p. Diamantopoulou, M.J., 2005. Artificial neural networks as an alternative tool in pine bark volume estimation. Computers and Electronics in Agriculture 48, 235–244. Diamantopoulou, M.J., 2010. Filling gaps in diameter measurements on standing tree boles in the Urban Forest of Thessaloniki, Greece. Environmental Modelling and Software 25, 1857–1865. Diamantopoulou, M.J., Milios, E., 2010. Modelling total volume of dominant pine trees in reforestations via multivariate analysis and artificial neural network models. Biosystems Engineering 105, 306–315. Diéguez-Aranda, U., Barrio Anta, M., Castedo Dorado, F., Álvarez González, J.G., 2005. Relación altura-diámetro generalizada para masas de Pinus sylvestris L. procedentes de repoblación en el noroeste de España. Forest Systems (Formerly: Inv Agrar: Sist Rec For) 14 (2), 229–241. Eler, Ü., 1986. Yield studies for Crimean Juniper in Turkey. Turkish Forest Research Institute, Technical Bulletin No. 192, Ankara, 40p. Fang, Z., Bailey, R.L., 1998. Height–diameter models for tropical forest on Hainan Island in southern China. Forest Ecology and Management 110, 315–327.

59

Fausett, L., 1994. Fundamentals of neural networks architectures. Algorithms and applications. Prentice Hall, USA, 461p.. Gadow, K.V., Real, P., Álvarez González, J.G., 2001. Modelización del Crecimiento y la Evolución de los Bosques. IUFRO World Series, vol. 12, Vienna, 242p. GDF, 2006. Forest Resources. The General Directorate of Forests, Ankara, 159p. Gurney, K., 1999. An Introduction to Neural Networks. UCL Press, London, 234p. Hasenauer, H., Kinderman, G., 2002. Methods for assessing regeneration establishment and height growth in uneven-aged mixed species stands. Forestry 75 (4), 385–394. Haykin, S., 1994. Neural Networks: A Comprehensive Foundation. Macmillan, N. York, USA. Huang, S., Price, D., Titus, S.J., 2000. Development of ecoregion-based height– diameter models for white spruce in boreal forests. Forest Ecology and Management 12, 125–141. Krumland, B.E., Wensel, L.C., 1988. A generalized height–diameter equation for coastal California species. Western Journal of Applied Forestry 3, 113–115. Lappi, J., 1997. A longitudinal analysis of height/diameter curves. Forest Science 43, 555–570. Larsen, D.R., Hann, D.W., 1987. Height–Diameter Equations for Seventeen Tree Species in Southwest Oregon. Research paper 49. Forest Research Laboratory, Oregon State University, Corvallis, OR, 16p. Leduc, D.J., Matney, T.G., Belli, K.L., Baldwin, V.C., 2001. Predicting Diameter Distributions of Longleaf Pine Plantations: A Comparison Between Artificial Neural Networks and Other Accepted Methodologies. USDA Forest Service, SRS025. Leite, H.G., Marques da Silva, M.L., Binoti, D.H.B., Fardin, L., Takizawa, F.H., 2011. Estimation of inside-bark diameter and heartwood diameter for Tectona grandis Linn. Trees using artificial neural networks. European Journal of Forest and Research 130, 263–269. Lindstrom, M.J., Bates, D.M., 1990. Nonlinear mixed effects for repeated measures data. Biometrics 46, 673–687. Litell, R.C., Milliken, G.A., Stroup, W.W., Wolfinger, R.D., Schabenberger, O., 2006. SAS for Mixed Models, second ed. SAS Institute Inc., Cary, NC, 814p.. López-Sánchez, C.A., Gorgoso Varela, J., Castedo Dorado, F., Rojo Alboreca, A., Rodríguez Soalleiro, R., Alvarez González, J.G., Sánchez Rodríguez, F., 2003. A height-diameter model for Pinus radiata D. Don in Galicia (Northwest Spain). Annals of Forest Science 60, 237–245. Meng, S.X., Huang, S., Yang, Y., Trincado, G., VanderSchaaf, C., 2009. Evaluation of population-averaged and subject-specific approaches for modeling the dominant/co-dominant height of lodgepole pine trees. Canadian Journal of Forest Research 39, 1148–1158. Morrison, M.L., Marcot, B.G., Mannan, R.W., 1992. Wildlife Habitat Relationships: Concepts and Applications. Univ. Wisconsin Press, Madison, 343p.. Nelson, R., Ranson, K.J., Sun, G., Kimes, D., Kharuk, V., Montesano, P., 2009. Estimating Siberian timber volume using MODIS and ICESat/GLAS. Remote Sensing of Environment 113, 691–701. Newton, P.F., Amponsah, I.G., 2007. Comparative evaluation of five height–diameter models developed for black spruce and jack pine stand-types in terms of goodness-of-fit, lack-of-fit and predictive ability. Forest Ecology and Management 247, 149–166. Özçelik, R., Diamantopoulou, M.J., Wiant, H.V., Brooks, J.R., 2010. Estimating tree bole volume using artificial neural network models for four species in Turkey. Journal of Environmental Management 91 (3), 742–753. Parresol, B.R., 1992. Baldcypress height–diameter equations and their prediction confidence interval. Canadian Journal of Forest Research 22, 1429–1434. Patterson, D., 1996. Artificial Neural Networks: Theory and Application. Prentice Hall, Singapore, 477p. Peng, C., Zhang, L., Huang, S., Zhou, X., Parton, J., Woods, M., 2001. Developing Ecoregion-Based Height–Diameter Models for Jack Pine and Black Spruce in Ontario. Forest Research Report 159. Ministry of Natural Resources. Ontario Forest Research Institute, Ontario. 10p. Pinheiro, J.C., Bates, D.M., 1998. Model Building for Nonlinear Mixed Effects Model. Department of Biostatistics and Department of Statistics, University of Wisconsin, Madison, Wis, 11p. Richards, F.J., 1959. A flexible growth function for empirical use. Journal of Experimental Botany 10, 290–300. Rumelhart, D.E., Hinton, G.E., Williams, R.J., 1986. Learning internal representations by error propagation. In: Rumelhart, D.E., McClelland, J.L. (Eds.), PDP Research Group, Parallel Distributed Processing, Foundations, vol. 1. MIT Press, pp. 318– 362. SAS Institute Inc., 2002. SAS/ETS User’s Guide, Version 9.0, SAS Institute Inc., Cary, NC. Schmoldt, D.L., Pei, L., Lynn Abbott, A., 1997. Machine vision using artificial neural networks with local 3D neighborhoods. Computers and Electronics in Agriculture 16, 255–271. Schnute, J., 1981. A versatile growth model with statistically stable parameters. Canadian Journal of Fisheries and Aquatic Sciences 38, 1128–1140. Sharma, M., Parton, J., 2007. Height–diameter equations for boreal tree species in Ontario using a mixed-effects modeling approach. Forest Ecology and Management 249, 187–198. Sharma, M., Zhang, S.Y., 2004. Height–diameter models using stand characteristics for Pinus banksiana and Picea mariana. Scandinavian Journal of Forest Research 19, 442–451. Soares, P., Tomé, M., 2002. Height–diameter equation for first rotation eucalypt plantations in Portugal. Forest Ecology and Management 166, 99–109.

60

R. Özçelik et al. / Forest Ecology and Management 306 (2013) 52–60

Spies, T.A., Cohen, W.B., 1992. An index of canopy height diversity. Coastal Oregon Productivity Enhancement Program Report 5. Oregon State University, Corvallis, OR, pp. 5–7. Swingler, K., 2001. Applying Neural Networks. A Practical Guide. Morgan Kaufman Publishers Inc., Great Britain, 303p. Temesgen, H., Gadow, K.V., 2004. Generalized height–diameter models – an application for major tree species in complex stands of interior British Columbia. European Journal of Forest Research 123, 45–51. Temesgen, H., Monleon, V.J., Hann, D.W., 2008. Analysis of nonlinear tree height prediction strategies for Douglas-fir forests. Canadian Journal of Forest Research 38, 553–565. Tomé, M., 1989. Modelação do crescimento da árvore individual em povoamentos de Eucalyptus globulus Labill. (1ª rotação) na região centro de Portugal. Ph.D. Thesis, Instituto Superior de Agronomía, Universidade Técnica de Lisboa, Lisbon, Portugal. 256 pp.

Trincado, G., VanderSchaaf, C.L., Burkhart, H.E., 2007. Regional mixed-effects height–diameter models for loblolly pine (Pinus taeda L.) plantations. European Journal of Forest Research 126, 253–262. Vonesh, E.F., Chinchilli, V.M., 1997. Linear and Nonlinear Models for the Analysis of Repeated Measurements. Marcel Dekker Inc., New York, 560 pp. Vospernik, S., Monserud, R.A., Sterba, H., 2010. Do individual-tree growth models correctly represent height-diameter ratios of Norway spruce and Scots pine. Forest Ecology and Management 260, 1735–1753. Wykoff, W.F., Crookston, N.L., Stage, A.R., 1982. User’s Guide to the Stand Prognosis Model. USDA Forest Service. General Technical Report INT-133, Intermountain Forest and Range Experimental Station, Ogden, UT, 113p. Yang, R.C., Kozak, A., Smith, J.H., 1978. The potential of Weibull-type functions as a flexible growth curve. Canadian Journal of Forest Research 8, 424–431. Zhang, L., 1997. Cross-validation of nonlinear growth functions for modeling tree height–diameter distributions. Annals of Botany 79, 251–257.