# Height-Diameter Modeling ofCinnamomum tamala Grown in Natural Forest in Mid-Hill of Nepal

Authors: Mukti Ram Subedi; Bishwa Nath Oli; Surendra Shrestha; Sophan Chhin
Journal: International Journal of Forestry Research
 (2018)

---

## Abstract

Tree height (H) and diameter at breast height (D) are key variables to calculate tree volume and biomass. We developed a height-diameter (H-D) model forCinnamomum tamala by evaluating 18 nonlinear models. Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Root Mean Square Error (RMSE), mean bias, Mean Absolute Error (MAE), graphical appearance, and biological logic were the criteria used to evaluate the predictive performance of the models. Gompertz model (M14) performed the best for predicting the total height ofC. tamala trees with the least RMSE (1.742 m), mean bias (0.012 m), and MAE (1.342 m) and satisfied model assumptions and biological logic. Validation data ranked the Gompertz model as the best model with RMSE (1.546 m), mean bias (-0.106 m), and MAE (1.149 m). Despite the consistent performance of the Gompertz model, it tended to underestimate the height prediction for taller (dominant crown class) and larger trees. Further work on refitting and validation of the proposed model with data from a larger geographic area, wider-ranging sites, and stand conditions is recommended.


---

## Body

## 1. Introduction

Tree height (H) and diameter at breast height (D) are fundamental variables in most forest inventories which are required to calculate tree volume, biomass, carbon storage, and survival analysis [1–4]. Accurate in situ measurement of D is easy and cost-effective. However, height measurement is labor-intensive, time-consuming, expensive, and prone to observational and measurement errors [5, 6]. Predicting forest dynamics through growth and yield simulation also requires individual tree level information, such as H and D. A height-diameter (H-D) relationship model can be easily built when both H and D variables are measured. The model then can be used to estimate missing tree height, biomass production, and stand dynamics.The allometric equation [7] serves as a tool to relate easily measurable morphometric variables (e.g., D) to the total height of the tree. Conventionally, H-D relationship models have been developed and applied to pure even-aged stands or plantations using D as a predictor variable [1, 2]. Recent studies employ other stand attributes (e.g., site quality, stand age, stand density, and dominant height) in mixed-effect H-D relationship models [8–13]. In mixed effects H-D relationship models, population-averaged (fixed parameters common to population), as well as subject-specific effects as random effects, are allowed. Mixed-effect H-D relationship models improve accuracy over nonlinear models that are based on the minimization of sums of squares.Generally, the site-specific H-D relationship models that are developed as H-D relationship can vary with differences in age, site quality, competition, stand age, and silvicultural treatment applied [14–17]. Although stand specific H-D relationship models are labor-intensive, costly, and time-intensive, these models have been proved to be more accurate in the realistic description of forest structure, growth simulations, and plot-level volume [18]. Due to a limited geographical range of naturalC. tamala stands, cross-sectional nature of data, and for simplicity of model use, we employed ordinary nonlinear least square equations. The model developed should be considered as a site-specific model as it was built considering fewer stands.Cinnamomum tamala T Nees and Eberm (family Lauraceae) is a moderate-sized, evergreen, monoecious tree species distributed at 900–2000 m above sea level in the tropical and subtropical Himalayas [19]. Natural stands ofC. tamala, being a moderate shade-tolerant species [20], are mostly found in shady-moist sites and grow well under full illumination.C. tamala is native to Nepal, China, India, and Bhutan [21].C. tamala trees are extensively managed, in situ and ex situ, for leaf and bark production due to dwindling natural population and high economic potentials [22]. The bark is traded to India [19] and is often used as a substitute forCinnamomum zeylanicum[23]. Department of forests, Nepal [24], has prioritizedC. tamala for research and management becauseC. tamala bark forms a significant part of national nontimber forest product trade, by both volume and economic value.Despite receiving a priority for research and management in Nepal and conservation in India, studies onC. tamala forests are limited. Since bark is a highly traded product, H-D relationship model is an essential decision-making tool both for management institutions and government for estimating forest volume and above-ground biomass, growth and yield modeling, and modeling ecophysiological process-based models of forests. However, research in this direction is limited within the distribution range ofC. tamala and has mostly focused on conservation, ethnopharmacological properties, dispersion pattern, cultivation, and harvesting practices [25–29]. In Nepal, Tree Improvement and Silviculture Component, TISC [30], and Poudel et al. [31] developed bark biomass production models in farm-grown trees and natural forests [32]. Equations for the H-D relationship forCinnamomum tamala grown in the natural forest have not been reported in past studies. Tree height increases nonlinearly with stem diameter [33, 34]. Therefore, linear models are inadequate for predicting tree height [35]. In this paper, we evaluated 18 nonlinear models to develop the H-D relationship model forC. tamala trees. To the best of our knowledge, this is probably the first attempt to establish an H-D relationship model for naturalC. tamala forests. Our model would be beneficial to researchers, managers, and academicians within the distributional range ofC. tamala. However, precaution should be used for accurate prediction of height outside the size and height ranges, site, and stand condition that differs from our study.


## 2. Materials and Methods

### 2.1. Study Area

The study was conducted at Mijure Danda Village Development Committee (VDC) of Kaski district, Nepal, which extends from 28°13′57′′ to 28°20′57′′N latitude and 84°08′53′′ to 84°12′42′′E longitude. The forest in the study area is managed by the Sikles unit of Annapurna Conservation Area Project (ACAP) and constitutesC. tamala as a dominant species having contiguous distribution pattern and with an Important Value Index (IVI) of 158.0 [27]. An IVI value was calculated by summing relative density, relative frequency, and relative dominance. The most frequent common associate species, with their IVI values, areSchima wallichii (77.89),Castanopsis indica (26.35,Madhuca latifolia (5.29), andMacaranga postulata (4.30) [27]. The study site falls under the humid subtropical monsoon climatic zone with hot and wet summers and relatively cold winters. The average monthly climate normal (1981-2010) temperature (minimum, maximum) to the nearest weather station was 21.1°C (13.4°C, 26.1°C) and precipitation was 325 mm (17.8 mm, 940.3 mm). Average annual precipitation was about 3900 mm year−1 [47].


### 2.2. Data Collection and Processing

We collected the data through the inventory of natural stands ofC. tamala, which includes 30 destructively sampled trees covering a wider range of H and D. Trees were selected for felling covering all the site conditions and tree size (H and D) distributions across the study sites in and outside of sampling plots. Nondestructive sample data came from a phytosociological survey of same forests [27]. We used systematic sampling plots with a size of 10 m × 10 m. Each plot has 5 m × 5 m subplots at northeast and southwest corners and 1 m × 1 m subplots at four corners. Information on species such as H and D were collected from 10 m × 10 m and 5 m × 5 m plots. The combined dataset size (n = 179) is smaller compared to other studies (e.g., [12, 13, 18, 48, 49]), but it can be used to construct H-D relationship models [50–52]. The combined dataset was split randomly into fitting (80%) and validation (20%) using R [53]. Summary statistics for H and D for both fitting and validation dataset are presented in Table 1.Table 1
Descriptive statistics of the variables for both fitting and validation dataset.
Statistic Fitting (n = 143) Validation (n = 36) Diameter (cm) Height (m) Diameter (cm) Height (m) Minimum 1.00 2.10 2.00 2.40 Mean 13.26 8.07 11.14 7.09 Maximum 32.00 14.20 30.00 13.58 Standard Deviation 7.72 3.26 7.47 3.15


### 2.3. Model Development and Evaluation

Initially, we selected 18 candidate nonlinear H-D relationship models (11 biparametric, and seven triparametric) based on a literature review (TableSA) to fit the model. We assumed that the total tree height is a function of diameter. In the later stage, we included only the models with all significant parameter estimates at 5% level of significance for further evaluation to compare their performance (Table 2). These models take the following functional form:(1)Hi=1.3+fDi,b+εiwhere Hi is the ith observation of the dependent variable tree height (m), Di is the ith observation of the independent variable diameter at breast height (cm), b is the vector of parameter to be estimated, εi is the random error term, f (.) is a nonlinear function, and i is the ith observation with i = 1, 2, 3, …, n. A constant term 1.3 is the height of stem above ground at which D was measured. We added this term in the function to avoid prediction of zero H when D approaches zero.Table 2
Models selected for further evaluation, their designation, and parameters of the models.
Model Designation Parameter Reference H i = 1.3 + b 0 D b 1 + ε i M1 2 [7] H i = 1.3 + b 0 1 - e x p - b 1 D b 2 + ε i M2 3 [33] H i = 1.3 + b 0 1 - e x p - b 1 D b 2 + ε i M3 3 [36] H i = 1.3 + D b 0 + b 1 D + ε i M4 2 [37] H i = 1.3 + exp ⁡ b 0 D b 1 + ε i M5 2 [38] H i = 1.3 + b 0 exp ⁡ - b 1 D + ε i M6 2 [39] H i = 1.3 + D 1 + b 0 D b 1 + ε i M7 2 [1, 2] H i = 1.3 + b 0 D 1 + D ) + ( b 1 D + ε i M8 2 [40] H i = 1.3 + D b 0 + b 1 D 3 + ε i M9 2 [41] cited in [13] H i = 1.3 + b 0 1 + 1 / b 1 D b 2 + ε i M10 3 [42] H i = 1.3 + D 2 b 0 + b 1 D 2 + ε i M11 2 [43] H i = 1.3 + b 0 D b 1 + D + ε i M12 2 [40, 44] H i = 1.3 + b 0 1 + b 1 exp ⁡ D b 2 + ε i M13 3 [45] H i = 1.3 + b 0 ⁡ exp ⁡ ( - b 1 × e x p - b 2 D ) ⁡ + ε i M14 3 [46] H i = 1.3 + b 0 ⁡ exp ⁡ ( - b 1 × D 2 ) + ε i M15 2 [2] H i is total height (m) of tree i = 1, 2, 3,…,n, D-diameter at breast height (cm), b0, b1, and b2 are parameters to be estimated, 1.3 is added to avoid prediction of zero height when D approaches zero, and εi is a random error term which is assumed to be normally and identically distributed with mean 0 and variance σ2 [NID (0, σ2)].The nonlinear function “nls2”, an improved version of the function nls, in statistical software R [54, 55] was used to estimate model parameters using the ordinary least square method. There are several criteria to assess model performance [56]. We evaluated model performance based on numeric analysis of the following: (1) asymptotic t-statistics of parameters, (2) Root Mean Square Error (RMSE: see (2)), also known as a measure of model’s precision, (3) mean bias (see (3)), and Mean Absolute Error (MAE: see (4)), and (4) information criteria—Akaike Information Criteria (AIC: see (5)) [56, 57] and Bayesian Information Criteria [58] (BIC: see (6)), followed by graphical presentation such as histogram, probability plots of residuals, and observed versus model predicted values. Criteria RMSE, AIC, BIC, Mean bias, and MAE are collectively called performance criteria. A coefficient of determination is not advisable to employ for assessing nonlinear regression. A statistic analogous to the coefficient of determination also called a fit index or modeling efficiency [14, 59] could be used. However, in this study, we did not use the fit index to evaluate our models as a fit index also suffers all the weaknesses of the coefficient of determination [14].(2)RMSE=1n-p-1∑i=1nyi-y^i2(3)Meanbias=∑i=1nyi-y^in(4)MAE=∑i=1nyi-y^in(5)AIC=n×ln⁡RSSn+2k(6)BIC=n×ln⁡RSSn+ln⁡n×kwhere yi,y^i are the observed and predicted values of height, respectively; n is the total number of observations used in model fitting; ln is natural logarithm; RSS is the residual sums of squares; and k is the number of parameters in the model.Models with least RMSE, AIC, BIC, mean bias, and MAE are considered to perform the best [60]. Reliability and validity of the model increase with its performance outside the test data [61–63]. However, due to limitation of dataset, model validation with an independent dataset covering more extensive geographical regions and management condition across the distributional ranges of species can be a potential topic for future research.We evaluated models based on the significance of model parameters, performance criteria, and practical application. After initial performance screening, models were subjected to further comparison by ranking the performance criteria and considering the practical application of the model. The models with the least value of performance criteria (RMSE, AIC, BIC, mean bias, and MAE) received the highest ranking and vice versa. The rank value of each criterion was added to get the total sum of ranking, which was then ranked in ascending order—the lower the value, the better the model. The model with the smallest rank value was considered the best H-D relationship model ofC. tamala naturally grown in mid-hill of Nepal.


## 2.1. Study Area

The study was conducted at Mijure Danda Village Development Committee (VDC) of Kaski district, Nepal, which extends from 28°13′57′′ to 28°20′57′′N latitude and 84°08′53′′ to 84°12′42′′E longitude. The forest in the study area is managed by the Sikles unit of Annapurna Conservation Area Project (ACAP) and constitutesC. tamala as a dominant species having contiguous distribution pattern and with an Important Value Index (IVI) of 158.0 [27]. An IVI value was calculated by summing relative density, relative frequency, and relative dominance. The most frequent common associate species, with their IVI values, areSchima wallichii (77.89),Castanopsis indica (26.35,Madhuca latifolia (5.29), andMacaranga postulata (4.30) [27]. The study site falls under the humid subtropical monsoon climatic zone with hot and wet summers and relatively cold winters. The average monthly climate normal (1981-2010) temperature (minimum, maximum) to the nearest weather station was 21.1°C (13.4°C, 26.1°C) and precipitation was 325 mm (17.8 mm, 940.3 mm). Average annual precipitation was about 3900 mm year−1 [47].


## 2.2. Data Collection and Processing

We collected the data through the inventory of natural stands ofC. tamala, which includes 30 destructively sampled trees covering a wider range of H and D. Trees were selected for felling covering all the site conditions and tree size (H and D) distributions across the study sites in and outside of sampling plots. Nondestructive sample data came from a phytosociological survey of same forests [27]. We used systematic sampling plots with a size of 10 m × 10 m. Each plot has 5 m × 5 m subplots at northeast and southwest corners and 1 m × 1 m subplots at four corners. Information on species such as H and D were collected from 10 m × 10 m and 5 m × 5 m plots. The combined dataset size (n = 179) is smaller compared to other studies (e.g., [12, 13, 18, 48, 49]), but it can be used to construct H-D relationship models [50–52]. The combined dataset was split randomly into fitting (80%) and validation (20%) using R [53]. Summary statistics for H and D for both fitting and validation dataset are presented in Table 1.Table 1
Descriptive statistics of the variables for both fitting and validation dataset.
Statistic Fitting (n = 143) Validation (n = 36) Diameter (cm) Height (m) Diameter (cm) Height (m) Minimum 1.00 2.10 2.00 2.40 Mean 13.26 8.07 11.14 7.09 Maximum 32.00 14.20 30.00 13.58 Standard Deviation 7.72 3.26 7.47 3.15


## 2.3. Model Development and Evaluation

Initially, we selected 18 candidate nonlinear H-D relationship models (11 biparametric, and seven triparametric) based on a literature review (TableSA) to fit the model. We assumed that the total tree height is a function of diameter. In the later stage, we included only the models with all significant parameter estimates at 5% level of significance for further evaluation to compare their performance (Table 2). These models take the following functional form:(1)Hi=1.3+fDi,b+εiwhere Hi is the ith observation of the dependent variable tree height (m), Di is the ith observation of the independent variable diameter at breast height (cm), b is the vector of parameter to be estimated, εi is the random error term, f (.) is a nonlinear function, and i is the ith observation with i = 1, 2, 3, …, n. A constant term 1.3 is the height of stem above ground at which D was measured. We added this term in the function to avoid prediction of zero H when D approaches zero.Table 2
Models selected for further evaluation, their designation, and parameters of the models.
Model Designation Parameter Reference H i = 1.3 + b 0 D b 1 + ε i M1 2 [7] H i = 1.3 + b 0 1 - e x p - b 1 D b 2 + ε i M2 3 [33] H i = 1.3 + b 0 1 - e x p - b 1 D b 2 + ε i M3 3 [36] H i = 1.3 + D b 0 + b 1 D + ε i M4 2 [37] H i = 1.3 + exp ⁡ b 0 D b 1 + ε i M5 2 [38] H i = 1.3 + b 0 exp ⁡ - b 1 D + ε i M6 2 [39] H i = 1.3 + D 1 + b 0 D b 1 + ε i M7 2 [1, 2] H i = 1.3 + b 0 D 1 + D ) + ( b 1 D + ε i M8 2 [40] H i = 1.3 + D b 0 + b 1 D 3 + ε i M9 2 [41] cited in [13] H i = 1.3 + b 0 1 + 1 / b 1 D b 2 + ε i M10 3 [42] H i = 1.3 + D 2 b 0 + b 1 D 2 + ε i M11 2 [43] H i = 1.3 + b 0 D b 1 + D + ε i M12 2 [40, 44] H i = 1.3 + b 0 1 + b 1 exp ⁡ D b 2 + ε i M13 3 [45] H i = 1.3 + b 0 ⁡ exp ⁡ ( - b 1 × e x p - b 2 D ) ⁡ + ε i M14 3 [46] H i = 1.3 + b 0 ⁡ exp ⁡ ( - b 1 × D 2 ) + ε i M15 2 [2] H i is total height (m) of tree i = 1, 2, 3,…,n, D-diameter at breast height (cm), b0, b1, and b2 are parameters to be estimated, 1.3 is added to avoid prediction of zero height when D approaches zero, and εi is a random error term which is assumed to be normally and identically distributed with mean 0 and variance σ2 [NID (0, σ2)].The nonlinear function “nls2”, an improved version of the function nls, in statistical software R [54, 55] was used to estimate model parameters using the ordinary least square method. There are several criteria to assess model performance [56]. We evaluated model performance based on numeric analysis of the following: (1) asymptotic t-statistics of parameters, (2) Root Mean Square Error (RMSE: see (2)), also known as a measure of model’s precision, (3) mean bias (see (3)), and Mean Absolute Error (MAE: see (4)), and (4) information criteria—Akaike Information Criteria (AIC: see (5)) [56, 57] and Bayesian Information Criteria [58] (BIC: see (6)), followed by graphical presentation such as histogram, probability plots of residuals, and observed versus model predicted values. Criteria RMSE, AIC, BIC, Mean bias, and MAE are collectively called performance criteria. A coefficient of determination is not advisable to employ for assessing nonlinear regression. A statistic analogous to the coefficient of determination also called a fit index or modeling efficiency [14, 59] could be used. However, in this study, we did not use the fit index to evaluate our models as a fit index also suffers all the weaknesses of the coefficient of determination [14].(2)RMSE=1n-p-1∑i=1nyi-y^i2(3)Meanbias=∑i=1nyi-y^in(4)MAE=∑i=1nyi-y^in(5)AIC=n×ln⁡RSSn+2k(6)BIC=n×ln⁡RSSn+ln⁡n×kwhere yi,y^i are the observed and predicted values of height, respectively; n is the total number of observations used in model fitting; ln is natural logarithm; RSS is the residual sums of squares; and k is the number of parameters in the model.Models with least RMSE, AIC, BIC, mean bias, and MAE are considered to perform the best [60]. Reliability and validity of the model increase with its performance outside the test data [61–63]. However, due to limitation of dataset, model validation with an independent dataset covering more extensive geographical regions and management condition across the distributional ranges of species can be a potential topic for future research.We evaluated models based on the significance of model parameters, performance criteria, and practical application. After initial performance screening, models were subjected to further comparison by ranking the performance criteria and considering the practical application of the model. The models with the least value of performance criteria (RMSE, AIC, BIC, mean bias, and MAE) received the highest ranking and vice versa. The rank value of each criterion was added to get the total sum of ranking, which was then ranked in ascending order—the lower the value, the better the model. The model with the smallest rank value was considered the best H-D relationship model ofC. tamala naturally grown in mid-hill of Nepal.


## 3. Results

### 3.1. Descriptive Statistics and Parameter Estimates

Model parameters were estimated using nonlinear least squares techniques for the model fitting dataset. Descriptive statistics includes mean ± st.dev (minimum, maximum) for both diameter and height [D (cm): 13.26 ± 7.72(1.0,32.0) and H (m): 8.07 ± 3.26(2.1,14.2)] for fitting dataset; and [D (cm): 11.14 ± 7.47(2,30) and H (m): 8.07 ± 3.15(2.4,13.58) for validation dataset]. All models were evaluated against multiple model performance criteria. Almost all (15) models (Table 3) showed significant parameter estimates (p<0.05, unless otherwise mentioned) and were considered for further evaluation.Table 3
Parameter estimates of the selected models.
Model Parameter Coefficient Std. error t p-value M1 b0 1.3491 0.1629 8.2812 <0.001 b1 0.6413 0.042 15.262 <0.001 M2 b0 13.7193 2.3975 5.7223 <0.001 b1 0.0544 0.0109 5.0132 <0.001 b2 1.0229 0.1495 6.84 <0.001 M3 b0 14.1312 2.3256 6.0765 <0.001 b1 0.0546 0.0242 2.2579 0.0255 b2 0.9905 0.1974 5.0186 <0.001 M4 b0 1.201 0.1003 11.9712 <0.001 b1 0.0467 0.0057 8.2513 <0.001 M5 b0 0.8545 0.0583 14.6447 <0.001 b1 0.3158 0.0228 13.844 <0.001 M6 b0 13.9604 0.6916 20.1857 <0.001 b1 7.9528 0.6879 11.5609 <0.001 M7 b0 14.3409 0.716 20.0279 <0.001 b1 8.6811 0.7223 12.0191 <0.001 M8 b0 0.8326 0.0696 11.9712 <0.001 b1 -0.9611 0.0078 -123.26 <0.001 M9 b0 1.4599 0.123 11.8697 <0.001 b1 0.4002 0.0081 49.4891 <0.001 M10 b0 20.4373 6.8756 2.9724 0.0035 b1 0.0388 0.0075 5.1557 <0.001 b2 1.0301 0.2058 5.0051 <0.001 M11 b0 1.6013 0.1337 11.9776 <0.001 b1 0.2454 0.0084 29.0816 <0.001 M12 b0 21.4115 2.5949 8.2514 <0.001 b2 25.715 5.1563 4.9871 <0.001 M13 b0 11.3241 0.498 22.7401 <0.001 b1 6.3863 1.1227 5.6884 <0.001 b2 0.0169 0.0019 8.9606 <0.001 M14 b0 12.1112 0.7606 15.9236 <0.001 b1 2.3069 0.2262 10.1997 <0.001 b2 0.115 0.0163 7.0766 <0.001 M15 b0 5.2026 0.2344 22.1983 <0.001 b1 -0.001 0.0001 -11.628 <0.001


### 3.2. Model Performance and Selection

Parameter estimates of models indicate that the best model cannot be determined solely based on significant t-statistics. We additionally considered performance criteria values to assess the model’s performance. Table4 presents the parameter estimates and fit statistics associated with the selected models for fitting and validation dataset.Table 4
Performance statistics of selected models for fitting and validation dataset.
Model Performance Criteria Fitting Validation Bias RMSE MAE AIC BIC Bias RMSE MAE (m) (m) (m) (m) M1 -0.016 1.801 1.381 576.972 585.861 -0.206 1.481 1.133 M2 0.033 1.777 1.382 574.167 586.019 -0.078 1.573 1.195 M3 0.025 1.777 1.381 574.185 586.036 -0.094 1.565 1.188 M4 0.020 1.776 1.381 573.093 581.981 -0.109 1.530 1.179 M5 -0.073 1.863 1.430 586.684 595.573 -0.341 1.523 1.195 M6 0.163 1.896 1.521 591.819 600.708 0.157 1.848 1.441 M7 0.136 1.869 1.494 587.717 596.606 0.112 1.794 1.393 M8 0.012 1.796 1.380 574.093 582.981 -0.129 1.639 1.179 M9 0.115 1.827 1.452 581.215 590.104 0.075 1.703 1.322 M10 0.025 1.782 1.383 575.072 586.923 -0.098 1.560 1.185 M11 0.090 1.804 1.427 577.579 586.467 0.029 1.644 1.277 M12 0.020 1.776 1.381 573.093 581.981 -0.109 1.530 1.179 M13 0.018 1.752 1.352 569.226 580.077 -0.190 1.546 1.138 M14 0.012 1.742 1.342 568.541 580.392 -0.106 1.546 1.149 M15 -0.050 2.419 1.972 661.457 670.346 -0.606 2.176 1.743Among 15 models, M14 produced the least bias (0.012 m), followed by M1 (0.016 m). The values of RMSE ranged from 1.742 m (M14) to 2.42 m (M15). M14 yielded the least MAE (1.342 m) and AIC(568.541), while M15 produced the highest AIC (661.457) and BIC (670.346). M14 performed significantly better than other models, produced similar parameter estimates for validation dataset (Table 6), and similar fit statistics for fitting and validation dataset (Table 4).We ranked models based on values of performance criteria and calculated average rank on sum ranks of performance criteria. M14 performed the best (rank = 1) for the fitting dataset, while M14, M12, and M4 performed the best with equal rank (rank = 1) for the validation dataset (Table5). Although M14, M12, and M4 have identical rank for validation dataset, we considered M14 the best model because it performed the best (rank = 1) for fitting dataset too, while M4 and M12 were the third best models for fitting dataset. Parameter estimates of M4 and M12 for validation dataset (Table SB), graph of the model, and residual distribution (Figures SA and SB) confirm that the models M4 and M12 are the best alternatives to model M14. Interestingly, M12 an M4 yielded strikingly similar performance statistics (identical to the three-decimal place) for model fitting, validation, and combined dataset (Tables 4 and SB).Table 5
Rank of selected models based on performance criteria for fitting and validation dataset.
Model Performance Criteria Fitting Rank Validation Rank Bias RMSE MAE AIC BIC Sum Bias RMSE MAE Sum (Rank) (Rank) M1 3 9 4 9 6 31 ( 6 ) 13 1 1 15 ( 4 ) M2 9 5 8 6 7 35 ( 8 ) 3 9 9 21 ( 8 ) M3 7 5 4 7 8 31 ( 6 ) 4 8 8 20 ( 7 ) M4 5 3 4 3 3 18 ( 3 ) 7 3 4 14 ( 1 ) M5 11 12 11 12 12 58 ( 11 ) 14 2 9 25 ( 11 ) M6 15 14 14 14 14 71 ( 15 ) 11 14 14 39 ( 14 ) M7 14 13 13 13 13 66 ( 13 ) 9 13 13 35 ( 13 ) M8 1 8 3 5 5 22 ( 5 ) 10 10 4 24 ( 10 ) M9 13 11 12 11 11 58 ( 11 ) 2 12 12 26 ( 12 ) M10 7 7 9 8 10 41 ( 9 ) 5 7 7 19 ( 5 ) M11 12 10 10 10 9 51 ( 10 ) 1 11 11 23 ( 9 ) M12 5 3 4 3 3 18 ( 3 ) 7 3 4 14 ( 1 ) M13 4 2 2 2 1 11 ( 2 ) 12 5 2 19 ( 5 ) M14 1 1 1 1 2 6 ( 1 ) 6 5 3 14 ( 1 ) M15 10 15 15 15 15 70 ( 14 ) 15 15 15 45 ( 15 )Table 6
Prediction statistics of the best model (M14) for validation dataset.
Dataset Parameter Coefficient Std. error t p-value RMSE Mean Bias Validation Data b0 17.63702 5.484934 3.21554 0.00291 1.439 0.0008 n = 36 b1 2.266303 0.218373 10.3781 <0.0001 b2 0.063602 0.023551 2.70063 0.01084


### 3.3. Residual Analysis and Shapiro-Wilk Test

To come up with the best model, we used test for residual normality. Residuals graphs (histograms, probability plots) and plot of fitted line overlaid on the observed heights data were produced for a visual check. Figure1 presents the plot of residuals against the observed diameters of fitting data. Most of the models showed no systematic departure from random pattern.Figure 1
Residuals analysis plots of 15 H-D relationship models fitted forCinnamomum tamala.The Shapiro-Wilk test of residuals of the best three models [M14- (W = 0.990, p = 0.4143) for fitting and (W = 0.978, p = 0.6858) for validation; M12- (W = 0.984, p = 0.627) for fitting and (W = 0.984, p = 0.871) for validation; M4- (W = 0.984, p = 0.627) for fitting and (W = 0.984, p = 0.871)] showed a p value greater than 0.05, which suggests that residuals do not violate the assumption of normal random error. Figure2 shows curves of the best fitting model M14 for fitting and validation dataset and histogram of residuals with a superimposed normal curve. Normally, the model tended to overestimate for shorter trees with larger diameter size.Figure 2
Scatter plot of total height (H) against D superimposed with M14 for fitting data (a) and validation data (c); a histogram of residuals of M14 for fitting data (b) and validation dataset (d).
(a) (b) (c) (d)


### 3.4. Paired t-Tests

The paired t-test of M14 with observed height values revealed insignificant statistics [t (142) = 0.084, p = 0.933]. Similarly, the paired t-test of observed height values with predicted heights of M12 [t (142) = 0.1734, p = 0.8625] and M4 [t (142) = 0.1368, p = 0.8914] for fitting dataset also resulted in statistically insignificant statistics. To determine the better model between M12 and M4, we compared the predicted height from both models between them and with the predicted height from the best model, i.e., M14. Our results of paired t-test to compare the predicted heights showed significant statistics between M4 and M12 [t (142) = 3.1143, p = 0.002] and between M12 and M14 [t (142) = -2.8968, p = 0.0044]. The M4 showed no significant difference in mean predicted heights from the M14 [t (142) = -0.8743, p = 0.3983]. This suggests that the M14 has more accurate prediction power followed by M4.


## 3.1. Descriptive Statistics and Parameter Estimates

Model parameters were estimated using nonlinear least squares techniques for the model fitting dataset. Descriptive statistics includes mean ± st.dev (minimum, maximum) for both diameter and height [D (cm): 13.26 ± 7.72(1.0,32.0) and H (m): 8.07 ± 3.26(2.1,14.2)] for fitting dataset; and [D (cm): 11.14 ± 7.47(2,30) and H (m): 8.07 ± 3.15(2.4,13.58) for validation dataset]. All models were evaluated against multiple model performance criteria. Almost all (15) models (Table 3) showed significant parameter estimates (p<0.05, unless otherwise mentioned) and were considered for further evaluation.Table 3
Parameter estimates of the selected models.
Model Parameter Coefficient Std. error t p-value M1 b0 1.3491 0.1629 8.2812 <0.001 b1 0.6413 0.042 15.262 <0.001 M2 b0 13.7193 2.3975 5.7223 <0.001 b1 0.0544 0.0109 5.0132 <0.001 b2 1.0229 0.1495 6.84 <0.001 M3 b0 14.1312 2.3256 6.0765 <0.001 b1 0.0546 0.0242 2.2579 0.0255 b2 0.9905 0.1974 5.0186 <0.001 M4 b0 1.201 0.1003 11.9712 <0.001 b1 0.0467 0.0057 8.2513 <0.001 M5 b0 0.8545 0.0583 14.6447 <0.001 b1 0.3158 0.0228 13.844 <0.001 M6 b0 13.9604 0.6916 20.1857 <0.001 b1 7.9528 0.6879 11.5609 <0.001 M7 b0 14.3409 0.716 20.0279 <0.001 b1 8.6811 0.7223 12.0191 <0.001 M8 b0 0.8326 0.0696 11.9712 <0.001 b1 -0.9611 0.0078 -123.26 <0.001 M9 b0 1.4599 0.123 11.8697 <0.001 b1 0.4002 0.0081 49.4891 <0.001 M10 b0 20.4373 6.8756 2.9724 0.0035 b1 0.0388 0.0075 5.1557 <0.001 b2 1.0301 0.2058 5.0051 <0.001 M11 b0 1.6013 0.1337 11.9776 <0.001 b1 0.2454 0.0084 29.0816 <0.001 M12 b0 21.4115 2.5949 8.2514 <0.001 b2 25.715 5.1563 4.9871 <0.001 M13 b0 11.3241 0.498 22.7401 <0.001 b1 6.3863 1.1227 5.6884 <0.001 b2 0.0169 0.0019 8.9606 <0.001 M14 b0 12.1112 0.7606 15.9236 <0.001 b1 2.3069 0.2262 10.1997 <0.001 b2 0.115 0.0163 7.0766 <0.001 M15 b0 5.2026 0.2344 22.1983 <0.001 b1 -0.001 0.0001 -11.628 <0.001


## 3.2. Model Performance and Selection

Parameter estimates of models indicate that the best model cannot be determined solely based on significant t-statistics. We additionally considered performance criteria values to assess the model’s performance. Table4 presents the parameter estimates and fit statistics associated with the selected models for fitting and validation dataset.Table 4
Performance statistics of selected models for fitting and validation dataset.
Model Performance Criteria Fitting Validation Bias RMSE MAE AIC BIC Bias RMSE MAE (m) (m) (m) (m) M1 -0.016 1.801 1.381 576.972 585.861 -0.206 1.481 1.133 M2 0.033 1.777 1.382 574.167 586.019 -0.078 1.573 1.195 M3 0.025 1.777 1.381 574.185 586.036 -0.094 1.565 1.188 M4 0.020 1.776 1.381 573.093 581.981 -0.109 1.530 1.179 M5 -0.073 1.863 1.430 586.684 595.573 -0.341 1.523 1.195 M6 0.163 1.896 1.521 591.819 600.708 0.157 1.848 1.441 M7 0.136 1.869 1.494 587.717 596.606 0.112 1.794 1.393 M8 0.012 1.796 1.380 574.093 582.981 -0.129 1.639 1.179 M9 0.115 1.827 1.452 581.215 590.104 0.075 1.703 1.322 M10 0.025 1.782 1.383 575.072 586.923 -0.098 1.560 1.185 M11 0.090 1.804 1.427 577.579 586.467 0.029 1.644 1.277 M12 0.020 1.776 1.381 573.093 581.981 -0.109 1.530 1.179 M13 0.018 1.752 1.352 569.226 580.077 -0.190 1.546 1.138 M14 0.012 1.742 1.342 568.541 580.392 -0.106 1.546 1.149 M15 -0.050 2.419 1.972 661.457 670.346 -0.606 2.176 1.743Among 15 models, M14 produced the least bias (0.012 m), followed by M1 (0.016 m). The values of RMSE ranged from 1.742 m (M14) to 2.42 m (M15). M14 yielded the least MAE (1.342 m) and AIC(568.541), while M15 produced the highest AIC (661.457) and BIC (670.346). M14 performed significantly better than other models, produced similar parameter estimates for validation dataset (Table 6), and similar fit statistics for fitting and validation dataset (Table 4).We ranked models based on values of performance criteria and calculated average rank on sum ranks of performance criteria. M14 performed the best (rank = 1) for the fitting dataset, while M14, M12, and M4 performed the best with equal rank (rank = 1) for the validation dataset (Table5). Although M14, M12, and M4 have identical rank for validation dataset, we considered M14 the best model because it performed the best (rank = 1) for fitting dataset too, while M4 and M12 were the third best models for fitting dataset. Parameter estimates of M4 and M12 for validation dataset (Table SB), graph of the model, and residual distribution (Figures SA and SB) confirm that the models M4 and M12 are the best alternatives to model M14. Interestingly, M12 an M4 yielded strikingly similar performance statistics (identical to the three-decimal place) for model fitting, validation, and combined dataset (Tables 4 and SB).Table 5
Rank of selected models based on performance criteria for fitting and validation dataset.
Model Performance Criteria Fitting Rank Validation Rank Bias RMSE MAE AIC BIC Sum Bias RMSE MAE Sum (Rank) (Rank) M1 3 9 4 9 6 31 ( 6 ) 13 1 1 15 ( 4 ) M2 9 5 8 6 7 35 ( 8 ) 3 9 9 21 ( 8 ) M3 7 5 4 7 8 31 ( 6 ) 4 8 8 20 ( 7 ) M4 5 3 4 3 3 18 ( 3 ) 7 3 4 14 ( 1 ) M5 11 12 11 12 12 58 ( 11 ) 14 2 9 25 ( 11 ) M6 15 14 14 14 14 71 ( 15 ) 11 14 14 39 ( 14 ) M7 14 13 13 13 13 66 ( 13 ) 9 13 13 35 ( 13 ) M8 1 8 3 5 5 22 ( 5 ) 10 10 4 24 ( 10 ) M9 13 11 12 11 11 58 ( 11 ) 2 12 12 26 ( 12 ) M10 7 7 9 8 10 41 ( 9 ) 5 7 7 19 ( 5 ) M11 12 10 10 10 9 51 ( 10 ) 1 11 11 23 ( 9 ) M12 5 3 4 3 3 18 ( 3 ) 7 3 4 14 ( 1 ) M13 4 2 2 2 1 11 ( 2 ) 12 5 2 19 ( 5 ) M14 1 1 1 1 2 6 ( 1 ) 6 5 3 14 ( 1 ) M15 10 15 15 15 15 70 ( 14 ) 15 15 15 45 ( 15 )Table 6
Prediction statistics of the best model (M14) for validation dataset.
Dataset Parameter Coefficient Std. error t p-value RMSE Mean Bias Validation Data b0 17.63702 5.484934 3.21554 0.00291 1.439 0.0008 n = 36 b1 2.266303 0.218373 10.3781 <0.0001 b2 0.063602 0.023551 2.70063 0.01084


## 3.3. Residual Analysis and Shapiro-Wilk Test

To come up with the best model, we used test for residual normality. Residuals graphs (histograms, probability plots) and plot of fitted line overlaid on the observed heights data were produced for a visual check. Figure1 presents the plot of residuals against the observed diameters of fitting data. Most of the models showed no systematic departure from random pattern.Figure 1
Residuals analysis plots of 15 H-D relationship models fitted forCinnamomum tamala.The Shapiro-Wilk test of residuals of the best three models [M14- (W = 0.990, p = 0.4143) for fitting and (W = 0.978, p = 0.6858) for validation; M12- (W = 0.984, p = 0.627) for fitting and (W = 0.984, p = 0.871) for validation; M4- (W = 0.984, p = 0.627) for fitting and (W = 0.984, p = 0.871)] showed a p value greater than 0.05, which suggests that residuals do not violate the assumption of normal random error. Figure2 shows curves of the best fitting model M14 for fitting and validation dataset and histogram of residuals with a superimposed normal curve. Normally, the model tended to overestimate for shorter trees with larger diameter size.Figure 2
Scatter plot of total height (H) against D superimposed with M14 for fitting data (a) and validation data (c); a histogram of residuals of M14 for fitting data (b) and validation dataset (d).
(a) (b) (c) (d)


## 3.4. Paired t-Tests

The paired t-test of M14 with observed height values revealed insignificant statistics [t (142) = 0.084, p = 0.933]. Similarly, the paired t-test of observed height values with predicted heights of M12 [t (142) = 0.1734, p = 0.8625] and M4 [t (142) = 0.1368, p = 0.8914] for fitting dataset also resulted in statistically insignificant statistics. To determine the better model between M12 and M4, we compared the predicted height from both models between them and with the predicted height from the best model, i.e., M14. Our results of paired t-test to compare the predicted heights showed significant statistics between M4 and M12 [t (142) = 3.1143, p = 0.002] and between M12 and M14 [t (142) = -2.8968, p = 0.0044]. The M4 showed no significant difference in mean predicted heights from the M14 [t (142) = -0.8743, p = 0.3983]. This suggests that the M14 has more accurate prediction power followed by M4.


## 4. Discussion

Almost all the candidate models performed well in explaining H-D relationship for C.tamala trees grown under natural forest. Model selection based on multiple criteria suggested that M14, M12, and M4 were more effective in predictingC. tamala tree height compared to other candidate models. The Gompertz model (M14) demonstrated competitive performance statistics for both fitting and validation dataset. Although Gompertz model follows growth theory ([8], Nord-Larsen and Cao 2006; Schmidt et al. 2011), it generally underestimates tree heights [10, 52]. Several studies have used the Gompertz model to describe the height-diameter relationship (Rupšys et al. 2011; Özçelik et al. 2013; [52, 64]). Zhang (1997) selected six models for evaluating H-D relationship and found that the Gompertz model underestimates tree heights for larger sizes trees (D>100 cm). Similarly, in our case, the Gompertz model tended to underestimate heights (10.84%) for taller trees with a larger diameter (H ≥12 m and D ≥17 cm). C.tamala trees from dominant crown category grow taller, larger, and faster. Trees under shade allocate more biomass to diameter growth than height as height growth requires more competitive advantage [65, 66]. This indicates thatC. tamala H-D relationship at our study site might have been mediated by shade tolerance or crown dominancy. Studies have shown that shade tolerance and climate jointly effect on allometric variation [66–68]. However, assessing the possible variation inC. tamala H-D relationship is beyond the scope of this study due to the small size of our dataset.H-D relationship models are often developed based on sampling data from permanent sample plots that are monitored for a relatively extended period [59, 64, 69, 70]. Studies that employed longitudinal data along with other plot-level covariates in nonlinear mixed H-D relationship models show improved model performance compared to ordinary nonlinear least square models [18, 52, 59, 71, 72]. Although the inclusion of stand variables may improve the prediction power of the model, this often requires greater sampling efforts especially when data collection involves an invasive sampling approach. Our study did not include such covariates because of the cross-sectional nature of data that were collected from both inside and outside of sample plots. Future researchers can test the validity of our model and improve the predictive performance utilizing stand and plot-level covariates when available.Predictive performance of the models beyond the observed range was not possible as the fitting dataset has larger sizes and taller trees. Independent datasets from different geographical locations and management systems with larger and taller trees than those used in this study could confirm the predictive power of the model [14]. Motallebi and Kangur [73], using the annual height and diameter increment data of 114 trees (observation periods spanned from 40 to 115 years) from three European countries, indicated a variable allometric relationship (significant coefficients across sites) between tree height and diameter which did not follow elasticity or geometric scaling rules across sites. H-D relationship varies among ecoregions due to differences in bioclimatic conditions [48, 69, 74] and among management regimes [32]. Therefore, H-D relationship of species should be evaluated across climatic and management systems to understand the changes in the stem allometry of species. AsC. tamala prefers moisture rich areas, water stress could change the H-D relationship. Therefore, the inclusion of water availability in original H-D relationship models and a mixed-effect modeling approach could further optimize the model’s performance. Since this study’s fitting and validation data came from same climatic and management conditions, validation of these models with independent data from different geographic, climatic, and management practices could be a potential topic for future research.


## 5. Conclusion

We calibrated and tested 18 nonlinear models onCinnamomum tamala trees grown under the natural multilayer tree structure in mid-hills of Nepal. The best model was selected based on multiple performance criteria, error diagnostics, growth theory, and practical application. Results showed that the Gompertz model (M14) best fitted the fitting dataset and performed consistently on the validation dataset and can be used reliably in biomass estimation and management forC. tamala trees in the mid-hills of Nepal.The existing studies onC. tamala height-diameter (H-D) modeling are limited to even-aged stands or plantation and have not attempted to develop H-D relationship models under natural forests. The methodology presented in this paper utilizes ordinary nonlinear models to generate H-D relationship in the uneven-aged natural forest ofC. tamala. Since we considered diameter as the only predictor variable, caution should be applied while using our model to other forests stands that differ in site and stand conditions to the basis of this study. In future research, our model can be extended to other regions and management regimes (e.g., the tree outside forestry), following updates through refitting and validation against independent data from the broadest possible ranges of size, site, and stand conditions, including stand and plot attributes across the distributional ranges ofC. tamala.



---

*Source: 6583948-2018-11-28.xml*