Evaluation of Individual Distance-Independent Diameter Growth Models for Japanese Cedar ( Cryptomeria japonica ) Trees under Multiple Thinning Treatments

: Understanding the tree growth process is essential for sustainable forest management. Future yields are a ﬀ ected by various forest management regimes such as thinning; therefore, accurate predictions of tree growth are needed under various thinning intensities. This study compared the accuracy of individual-level distance-independent diameter growth models constructed for di ﬀ erent thinning intensities (thinning intensity-dependent multiple models: TDM model) against the model designed to include all thinning intensities (thinning intensity-independent single model: TIS model) to understand how model accuracy is a ﬀ ected by thinning intensity. We used long-term permanent plot data of Japanese cedar ( Cryptomeria japonica ) stands in Japan, which was gathered from four plots where thinning was conducted at di ﬀ erent thinning intensities: (1) intensive (41% and 38% of trees removed at 25 and 37 years old, respectively), (2) moderate (38% and 34%), (3) light (32% and 34%), and (4) no thinning. First, we speciﬁed high interpretability distance-independent competition indices, and we compared the model accuracy both in TDM and TIS models. The results show that the relative spacing index was the best competition index both in TDM and TIS models across all thinning intensities, and the di ﬀ erences in the RMSE (Root mean square error) and rRMSE (relative RMSE) in both TDM and TIS models were 0.001–0.01 cm and 0.2–2%, respectively. In the TIS model, rRMSE varied with thinning intensity; the rRMSE was the lowest for moderate thinning intensity (45.8%) and the highest for no thinning (59.4%). In addition, bias values were negative for the TIS model for all thinning intensities. These results suggest that the TIS model could express diameter growth regardless of thinning intensities. However, the rRMSE had varied with thinning intensity and bias had negative values in the TIS model. Therefore, more model improvements are required for accurate predictions of long-term growth of actual Japanese cedar stands.


Introduction
Yield predictions are essential for sustainable forest management. Generally, future yields are affected by different forest management regimes, site conditions, other climate factors, and species [1,2]. Such growth processes are described by growth models [3]. Growth models can be classified into stand-level and individual-level models. The stand-level model, as the basis of growth prediction, has been studied for a long time [2]. This model can characterize yield and growth with a little information about the stand; therefore, many researchers have adapted this model for a variety of

Data Collection
The data used in this study were obtained from Kudarukawa national forest, located in Kochi prefecture, western Japan (33 • 23' N, 133 • 10' E) ( Figure 1). The study site is located in areas with an inclination of 5 • -10 • at an elevation of 500 m. In 1972, four permanent plots were established in planted forests with different thinning intensities, which were: (1) intensive thinning, (2) moderate thinning, (3) light thinning, and (4) no thinning [26]. Each plot included approximately 200 Japanese cedar trees. The first thinning was conducted at 25 years old, and the second thinning was conducted at 37 years old (Table 1). Trees' heights and diameters at breast height (DBH) in each plot were measured every 5-6 years from 14 years old to 61 years old. The DBH of each tree was measured at breast height position (1.2 m) from two directions using a caliper, whereas tree height was measured for approximately 30 selected trees with various DBH classes in each plot ( Table 2). A tree number and breast height line were painted for each tree when tree measurements were conducted. Before our analysis, we compared the mean annual diameter growth for 47 years among four plots using a multiple-comparison test. The mean annual diameter growth of the unthinned plot was significantly smaller than that of the other thinned plots (Tukey-Kramer test, p < 0.05). Moreover, we estimated the heights of unmeasured trees using the diameter-height equation, created using the data of measured trees (R 2 values in each thinning intensity were 0.88-0.99). To calculate the competition indices, we assumed that thinning was conducted in prior measurement periods in this study. In total, we used nine growing seasons of 1077 trees' growth data.

Data Collection
The data used in this study were obtained from Kudarukawa national forest, located in Kochi prefecture, western Japan (33°23' N, 133°10' E) ( Figure 1). The study site is located in areas with an inclination of 5°-10° at an elevation of 500 m. In 1972, four permanent plots were established in planted forests with different thinning intensities, which were: (1) intensive thinning, (2) moderate thinning, (3) light thinning, and (4) no thinning [26]. Each plot included approximately 200 Japanese cedar trees. The first thinning was conducted at 25 years old, and the second thinning was conducted at 37 years old (Table 1). Trees' heights and diameters at breast height (DBH) in each plot were measured every 5-6 years from 14 years old to 61 years old. The DBH of each tree was measured at breast height position (1.2 m) from two directions using a caliper, whereas tree height was measured for approximately 30 selected trees with various DBH classes in each plot ( Table 2). A tree number and breast height line were painted for each tree when tree measurements were conducted. Before our analysis, we compared the mean annual diameter growth for 47 years among four plots using a multiple-comparison test. The mean annual diameter growth of the unthinned plot was significantly smaller than that of the other thinned plots (Tukey-Kramer test, p < 0.05). Moreover, we estimated the heights of unmeasured trees using the diameter-height equation, created using the data of measured trees (R 2 values in each thinning intensity were 0.88-0.99). To calculate the competition indices, we assumed that thinning was conducted in prior measurement periods in this study. In total, we used nine growing seasons of 1077 trees' growth data.

Individual-Tree Diameter Growth Model
We used the linear-mixed effects model to characterize Japanese cedar diameter growth. We defined individual diameter growth of the TDM model as follows: where G I,i, j is the growth of the subject tree i in the I th plot in the period between the jth measurement and the subsequent measurement. DBH I,i,j and DBH 2 I,i,j are DBH and the DBH squared, respectively, and Age I,i,j is the stand's age. Generally, tree diameter growth is affected by individual diameter size [27]. Moreover, tree diameter growth follows a sigmoid curve, indicating the existence of restraint by aging and increasing DBH [3]. To express this restraint, we used a simple quadratic function of DBH. α 0 , . . . , α 3 are each parameters, and ϕ i is a random parameter for subject tree i. ε I,i,j is the error term produced by the normal distribution with mean 0 and variance σ 2 ( ε I,i,j ∼ Normal 0, σ 2 ). The annual diameter growth in the period between two successive measurements was expressed using DBH, DBH 2 , and the age at the beginning of the period. A full TDM model derived by adding each distance-independent competition index CI I,i,j is given by where α 0 , . . . , α 4 are fixed effects parameters, and CI I,i,j is the distance-independent competition index of the Ith plot in subject tree ith at jth measurements. We estimated the parameters of Equation (2) for each of the four plots, with varying thinning regimes. Thus, we obtained a total of 20 fixed effects parameters for the TDM model. Corresponding to Equations (1) and (2), the TIS model was defined as follows: where ϕ I,i is the random parameter for the Ith plot in subject tree i. The TIS model is a single model that expresses diameter growth for all plots because the model includes the plot-level random effect ϕ I . In this study, we used six distance-independent competition indices as described in the following section. The generalized liner mixed effects model was used in the lme4 package [28] in R version 3.6.0 [29]. Also, we used a Correction Factor (CF = exp((MSE)/2, where MSE is the mean square error) to predict expected diameter growth. We applied the CF to correct bias occurring in the log transformation [30].

Distance-Independent Competition Index
We used six distance-independent competition indices: basal area of trees larger than the subject tree (BAL), diameter ratio (DR), basal area ratio (BR), cumulative distribution function (CDF), relative spacing index (Sr), and SDI according to the definition of Sun et al [11].

Basal Area of Trees Larger than the Subject Tree (BAL)
BAL is defined as where BAL I,i,j is the sum of the basal area of competitor trees for subject tree i in Ith plot at the jth measurement, and DBH c is competitor trees that have larger diameter than subject tree i in the Ith plot at the jth measurement [14,31].

Diameter Ratio (DR)
DR is defined as where DR I,i,j is the DR of Ith plot in subject tree i at the jth measurement. DBH I,i,j is the DBH of subject tree i at the jth measurement, and QMD I,j is the quadratic mean diameter of the Ith plot at the jth measurement [32].

Basal Area Ratio (BR)
BR is defined as where BR I,i,j and BA I,i,j are the BR and BA of the Ith plot in subject tree i at the jth measurement, respectively. BA I,j is the mean BA of the Ith plot at the jth measurement. In this study, we used the square of Equation (5) [11].

Cumulative Distribution Function (CDF)
CDF is usually used to express the diameter distribution in the stand. Also, CDF expresses the relative position F I,i,j of the Ith plot in subject tree i at the jth measurements as follows: where k is the rank of the size of subject tree i ordered from smallest to largest in the Ith plot. N I,j is the number of trees per hectare in the Ith plot at the jth measurement.

Relative Spacing Index (Sr)
Sr is a stand density index that is calculated as the ratio of distance between stand trees and mean height as follows: where Sr I,j and H I,j are Sr and mean height in the Ith plot at the jth measurement, respectively [33].

Site Density Index (SDI)
Similarly, SDI is expressed by the number of trees in the stand and quadratic mean diameter as follows [34]: where c is a constant defined by a self-thinning model. In this study, we used c = 1.486, which was evaluated in a Japanese cedar stand [35,36].

Model Fitting and Evaluation
We evaluated the goodness of the fit of the model for each competition index based on the marginal R 2 , conditional R 2 , Akaike's information criteria (AIC), and ∆AIC. Marginal R 2 and conditional R 2 were calculated using observed and predicted diameter growth. Marginal R 2 was calculated with only the fixed effect, and conditional R 2 was calculated with both fixed and random effects [37]. ∆AIC is the difference between the AIC of the best model and each other model. Larger values of marginal R 2 and conditional R 2 indicate better model accuracy; conversely, a lower value of AIC indicates better model accuracy. Additionally, to evaluate the effectiveness of the competition indices under different thinning intensities, we calculated the root mean square error (RMSE), relative RMSE (rRMSE), bias (mean error), and variance error using the best model. rRMSE is calculated as RMSE divided by mean diameter of observed diameter growth, and a smaller rRMSE indicates better prediction accuracy [38]. Accuracy is higher when bias is closer to 0. Negative bias values indicate underestimation, and vice versa.

The Goodness of Fit for Distance-Independent Competition Index
In the TDM model, the best model included Sr for each thinning intensity based on marginal R 2 , conditional R 2 , AIC and ∆AIC (Table 3). Then, the marginal R 2 and conditional R 2 of each thinning intensity were 0.36-0.49 and 0.62-0.70, respectively ( Table 3). The second-best model included BAL in intensive and light thinning, and SDI in moderate thinning and no thinning ( Table 3). Coefficients of the TDM model in all thinning intensities showed that DBH, DBH 2 , age, and Sr were 0.174 to 0.213, −0.002 to −0.001, −0.056 to −0.030, and 0.161 to 0.209, respectively ( Table 4). The TIS model showed that the best model included Sr (Table 5). Then, for the best model, marginal R 2 and conditional R 2 were 0.42 and 0.67, respectively (Table 5). Also, the coefficients of the best model showed that DBH, DBH 2 , age, and Sr were 0.177, −0.002, −0.040, and 0.163, respectively (Table 6).

Comparison of Model Accuracy between TDM Model and TIS Model
Visual assessment of the relationships between actual and predicted DBH growth showed little difference in prediction performance between TDM and TIS models of each thinning intensity (Figure 2). RMSE and rRMSE also differed little for both TDM and TIS models of all thinning intensities ( Table 7). The differences in RMSE and rRMSE among TDM and TIS models were 0.001-0.01 and 0.2-2%, respectively ( Table 7). The bias of the TDM and TIS models specifically differed in no thinning, with a difference of 0.032. Variance decreased with decreasing thinning intensity, and the lowest variances for no thinning in both TDM and TIS models were 0.039 and 0.033, respectively. Table 3.
Evaluation of thinning intensity-dependent multiple (TDM) models using each distance-independent competition index.

Model Accuracy under Different Thinning Intensities
For the TDM model, the lowest RMSE value was 0.203 for no thinning, and the highest RMSE value was 0.249 for intensive thinning (Table 7). rRMSE was the lowest for moderate thinning (47.8%), and the highest for no thinning (60.6%). Bias values were negative for all thinning intensities, indicating that the predicted values were underestimated compared with the observed values. For the TIS model, RMSE values showed that the lowest value was 0.199 for no thinning, and the highest value was 0.250 for intensive thinning (Table 7). rRMSE values were lowest for moderate thinning (45.8%), and highest for no thinning (59.4%). Bias values were negative for all thinning intensities, indicating that the predicted values were underestimated compared to the observed values.   Table 4, and the best model of TIS model according to Table 6.

Model Selection for TDM and TIS Models
Several studies that used distance-independent competition indices showed that the best competition indices were BAL, SDI, and CDF [6,10,20]. Goodness of fit of competition indices may be species-specific. In this study, Sr was the best competition index both in the TDM and TIS models (Tables 3 and 5). This result indicates that individual level competition is expressed as other explanatory variables such as DBH and DBH squared because Sr is expressed as stand level competition. Moreover, a previous study showed that Sr was related to crown length [39]. Tree crown is widely used for predicting tree growth, and there is high correlation between crown and tree growth [9,40]. Sr may indirectly reflect information of the tree crown; therefore, the best models might perform well. Moreover, tree height is suited to express site productivity. Sr could reflect differences in productivity among each study plot, therefore increasing the accuracy of the model. In the TDM model, the second-best models were different for each thinning intensity (Table 3). BAL was selected in intensive and light thinning, and SDI was selected in moderate and no thinning. This result suggests that BAL and SDI may also be possible to express diameter growth of Japanese cedar under different thinning intensities.  Table 4, and the best model of TIS model according to Table 6.

Model Selection for TDM and TIS Models
Several studies that used distance-independent competition indices showed that the best competition indices were BAL, SDI, and CDF [6,10,20]. Goodness of fit of competition indices may be species-specific. In this study, Sr was the best competition index both in the TDM and TIS models (Tables 3 and 5). This result indicates that individual level competition is expressed as other explanatory variables such as DBH and DBH squared because Sr is expressed as stand level competition. Moreover, a previous study showed that Sr was related to crown length [39]. Tree crown is widely used for predicting tree growth, and there is high correlation between crown and tree growth [9,40]. Sr may indirectly reflect information of the tree crown; therefore, the best models might perform well. Moreover, tree height is suited to express site productivity. Sr could reflect differences in productivity among each study plot, therefore increasing the accuracy of the model. In the TDM model, the second-best models were different for each thinning intensity (Table 3). BAL was selected in intensive and light thinning, and SDI was selected in moderate and no thinning. This result suggests that BAL and SDI may also be possible to express diameter growth of Japanese cedar under different thinning intensities.

Diameter Growth Model Accuracy for TDM and TIS Models
Previous studies evaluated model accuracy using indexes such as RMSE, rRMSE, and AIC [9,10,41]. For example, Palahí et al. [41] constructed an individual-tree distance-independent diameter growth model for Scots pine (Pinus sylvestris L.) in northeast Spain using long-term growth data. They calculated model accuracy using rRMSE, and the rRMSE in their model was 64.1%. On the other hand, in our study, the rRMSE values of TDM and TIS models were 45.8-60.6% (Table 7). Thus, our model accuracy was not inferior in comparison with previous studies.
There was little difference in RMSE and rRMSE among TDM and TIS models (Table 7). Also, the best competition index was the same for both TDM and TIS models (Tables 3-6). These results indicate that the TIS model can explain diameter growth under different thinning intensities. If the TIS model can explain the response of diameter growth under different thinning intensities, even a single model could estimate tree growth in stands with various thinning intensities. rRMSE showed that there were differences in the TIS model under each thinning intensity (Table 7). Specifically, the rRMSE values for the model of intensive thinning and no thinning were higher than those of other thinning intensities, indicating low model accuracy. The model may not suit specific thinning intensities, because the TDM model had a similar tendency. On the other hand, both rRMSE and bias were smallest with moderate thinning; therefore, the model can roughly predict DBH growth for forests with moderate thinning intensities.
There were not significant differences in diameter growth among three thinned plots in our dataset, as shown in "2.1. Data Collection". Therefore, our analysis could not examine in detail the behavior of TDM and TIS models under conditions in which there are large differences in diameter growth. On the other hand, there were significant differences in diameter growth among thinned plots and no thinning plot in our dataset. In addition, the TIS model accurately predicted diameter growth in thinned and no thinning plots, suggesting that the TIS model is effective for the multiple plots with largely different diameter growth. However, this suggestion is not based on actual data from the multiple plots with largely different diameter growth. Thus, the discussion here requires further analysis.

Application for Forest Management
Our results suggest that the TIS model can express diameter growth under various thinning intensities. In this study, the model predicts individual-tree diameter growth. When the model is applied to actual stands, the diameter size is predicted based on tree diameter sizes in the prior period. The long-term growth is repeatedly predicted. In addition, stand-level growth is predicted by accumulating individual tree-level growth. If the error is large, such as in intensive and no thinning, error will propagate. The more we predict long-term growth, the more the error will increase between actual stands and predicted values. Also, bias values were negative for all thinning intensities in the TIS model, indicating the predicted diameter growth was underestimated ( Table 7). Underestimation of tree diameter is problematic for optimal operation planning, such as determining thinning intensity. Stands with low tree diameter usually require intensive thinning to allow residual trees to grow well; thus, underestimated tree diameter may lead to unnecessary thinning intensity. For those reasons, we need to improve models in which the error is constant with thinning intensity and bias values are not negative.

Conclusions
This study compared the model accuracy of TDM and TIS models with different thinning intensities. For both the TDM and TIS models, the selected competition index was the same and there was little difference in model accuracy. These results suggest that the TIS model could express diameter growth regardless of thinning intensity. This result contributes to yield prediction techniques for various thinning intensities. There are little differences in diameter growth among thinned plots in our study site. We need to consider model construction under the variety of thinning intensities. On the other hand, the TIS model had varying error under each thinning intensity, and the model accuracy was low in intensive thinning and no thinning. Also, bias values were negative; thus, predicted diameter size may be underestimated when the model is adapted to actual stands. These results suggest that our model has limitation regarding adaptation in unusual thinning intensities or other sites. Therefore, we need to improve the model to generate accurate predictions of the long-term growth of actual stands. We focused on individual tree growth in terms of tree competition in this study. Recently, some researchers have mentioned that site productivity is affected by climate change [42,43]. Thus, we need to consider factors such as temperature, precipitation, and carbon concentration when predicting long-term individual or forest stand growth.
Author Contributions: K.F., T.N. and F.K. conceived and designed the study. F.K. provided field measurements data. All authors contributed to data analysis. K.F and T.N. wrote the paper. All authors have read and agreed to the published of the manuscript.