Determining Ideal Timing of Row Thinning for a Cryptomeria japonica Plantation Using Event History Analysis

Effective time of thinning is essential for determining a silvicultural operation schedule. One of the most commonly used methods is the percentage of radial increase to assess the effect of thinning. However, it is difficult to determine the ideal time point due to variation in tree growth rates. Event history analysis was used to quantify the optimal timings for different row thinning types for a 45-year-old Cryptomeria japonica plantation in the mountainous region of Taiwan. The increase in tree-ring size was measured and converted to the basal area increment (BAI) to estimate annual tree growth; derived time-series data were entered into event history analysis to calculate the time to 50% probability of survival. Additionally, an accelerated failure time regression was applied to test the effects of thinning and its timing; model validation was carried out to examine the influence of thinning time variation on plant growth through time. Results showed that thinning modified the temporal dynamics of the BAI, and, in general, a positive trend was observed between strip-width and time of thinning. Simulated tree growth in the model validation corroborated that accurate timing may optimize thinning effects. Combining tree-ring measurement and event history analysis may facilitate determining the timing of row thinning, which can improve carbon sequestration of forest stands.


Introduction
Thinning is a forest management strategy to enhance tree growth by increasing spacing and adjusting the density of a stand.It is also commonly applied to improve the residual stems growing condition [1], to enhance biodiversity [2], and to reduce the risk of stand damage [3].One major challenge of thinning is to determine the effective time of thinning, which varies among biotic (e.g., species types, stand density), abiotic (mainly substrate fertility), and management-related factors.Sophisticated methods such as dynamic programming, numerical optimization, or stand growth models can be used to derive the information [4][5][6][7][8][9][10][11][12].However, most of these approaches are complex and require a wide array of long-term monitoring data including thin age, site index, initial stocking, residual stocking and/or rotation length that make the implementations challenging [13][14][15][16].Therefore, a straightforward and robust method is required.
Event history analysis (also known as survival analysis) is a statistical method to determine the probability of timing and duration to the occurrence of an event (time-to-event).It is a means in which to examine the distribution of event times often using a hazard function.The hazard function can be described in terms of the cumulative hazard function, so that we can transpose signs and fit in an exponential form to determine the connection between the hazard and survivor functions.A survival function can be used as the proportion of a given amount of time and characterizes the increase or decrease in the survival probability [17].Event history analysis is commonly used in ecology to examine tree and seedling survival in relation to various environmental factors [18][19][20][21].
Dendroecology is the science of utilizing dendrochronological data (e.g., historical (~100-1000-year time spans) climatic records derived via dating and analyzing annual rings of trees or old timber) to analyze historic ecological processes [22].Variations in ring widths reflect the influences of stand hierarchy, local site factors, climate, and anthropogenic activities.Dendroecological analysis can also precisely document perturbations such as drought [23], fire history [24], and/or canopy gap dynamics [25].Dates of apparent growth releases recorded on surviving trees may indicate the occurrence of disturbance events [26].From this, the effect of thinning should decrease over time, which can be treated as the time of failure in event history analysis.Therefore, we may be able to determine the optimal time of thinning by including tree-ring records in the analysis.Additionally, tree growth through time can also be retrieved by using this dendroecological data, which may facilitate assessing the effects of thinning.In this study, we aim to investigate the feasibility of event history analysis for determining optimal time of thinning with dendroecological data of Cryptomeria japonica (L.) D. Don (Japanese cedar).

Study Site and Thinning Treatments
The study site is a 2-ha C. japonica plantation stand located in the mountainous region (Chilan Mountain) of northeastern Taiwan (24 • 37 4.60 N, 121 • 29 23.37 E) (Figure 1) at an elevation of 1100 m, with a mean annual precipitation and air temperature of 3600 mm pre year and 13.3 • C, respectively.Ultisol and inceptisol are the main substrate types.C. japonica seedlings were planted at a density of 3300 trees ha -1 in 1966 right after a clear-cut; mean diameter at breast height (DBH, defined as stem diameter 130 cm above the forest floor) was approximately 20 cm in 1990 when thinning treatments were applied [27].Four thinning types (2-5 row (R) cuttings) with identical 2-m row widths were applied, and six rows were retained on either side of the cut strips (Figure 2).Each one was assigned to three 0.04-ha (20 m × 20 m) sub-plots (total of 12 plots), and three control plots (CTRL) of an identical size, without any perturbations since 1966, were preserved.C. japonica is the most widely planted species in Taiwan, occupying over 40,000 ha of the forested area [28].

Tree-Ring Data Collection and Analysis
Trees on the edges were sampled to assess the contribution of thinning to growth.In the summer of 2011, 20 or 21 pith-to-bark increment core specimens from each treatment and the CTRL group were collected (n = 104) by referring to the size distribution of DBH for each group (Figure 2).Each core was sampled at breast height (130 cm above ground).Trees may invest more on the canopy layer to compete for space and light, and slow the radial growth at breast height [29], which could potentially influence the estimation.However, the impact could be negligible since all sampled trees were along the edges and may have had sufficient space to grow (Figure 2b-e) during the period 1990-2011.The visual cross-dating was carried out for the individual core, and each series was compared to get master chronology and correctly dated.Tree-ring widths were then measured using the LINTAB system (RINNTECH ® , Heidelberg, Germany), and the performance of cross-dating was assessed by COFECHA [30].Since the study site was situated on the east-facing slope, trees were cored along the north-south direction to minimize the reaction wood effect [31].For the same reason, we also avoided taking samples from the upslope or downslope side.

Tree-Ring Data Pre-Processing
Tree-ring data are continuous measures of annual tree diameter increments, which are often influenced by climate.In our case, this would obscure the signal of failure events in event history analysis.Therefore, as suggested by Peltola et al. [32], a three-year moving average window was applied to smooth the time-series data to mitigate effects of short-term climate fluctuations on tree-ring widths.This should be appropriate since the region was wet (mean annual precipitation greater than 3000 mm) [33] and annual temperature fluctuation for this near-tropical humid island was much less pronounced in comparison to other regions [34].To assess the effects of thinning, the raw tree-ring series was divided into two periods, pre-(1966-1990) and post-(1991-2011) thinning.We note that tree-ring data have often been detrended using Spline, ModNegExp, Ar, Friedman, or other approaches before further analyses.However, it was not the case for this study since event history analysis was applied to assess the gradually declined effect of row thinning.A metric basal area increment (BAI), which is essential information in forest management, was derived for the analysis (Equation ( 1)): where rd and a are the radius measured from ring width and tree age, respectively.The BAI positively culminates in a linear phase based on the maturity stage in general [35] and a negative trend in the BAI is a clear indication of decline in tree growth [36].To calculate the effective length (accelerated growth), these two periods of tree-ring measures must be transformed to time-to-event binary data by setting a threshold as the following: effective (1): smoothed BAI > threshold, or failure (0): smoothed BAI < threshold.The pre-thinning threshold was the mean of the smooth BAI from 1966 to 1990, and for the post-thinning period, it was the mean of smoothed BAI (baseline) right before the treatment (average of [1988][1989][1990].When the growth rate was no longer greater than the rate prior to thinning, the time of failure was then assigned.Once the time-series data failed to reach the thresholds, failure (0s) were assigned thereafter.

Statistical Analysis
The Kaplan-Meyer method is a type of event history analysis that is specifically designed to estimate the following non-parametric survival function [37] (Equation ( 2)): where di is the number of event absences (0) within a threshold setting for each individual in a given time interval i, and ni is the number of event occurrences at risk (the probability of failing to reach the threshold).We examined the time effects between treatments by fitting an accelerated failure time model regression [19] (Equation ( 3)): where T is time to the event; X and β are vectors of explanatory variables and coefficients, respectively; σ is an error term; and W is a linking function assumed to have a specified distribution (e.g., Weibull, exponential, or Gaussian).Akaike's Information Criterion (AIC) [38] was used to assess the goodness of fit of model estimates [21].All of these aforementioned statistical analyses were conducted in R survival package (v.3.0.2;Stanford University, Stanford, CA, USA; http://www.r-project.org/).

Model Validation
Model validation was conducted to assess the performance of event history analysis.The polynomials were derived to fit the observed growth trends (1991-2011) of different treatments (2-5R), which were used to simulate tree-ring increments for the study period.Wide but reasonable ranges for repeat thinning (5-25 years) [39] and observation (100 years from the initiation of thinning in 1990) periods were set for the comparisons (%), which were the relative percentages defined as ratios of mean annual growth gains (cm•tree -1 •year -1 ) of different repeat thinning periods and treatment, and the corresponding event history analysis-derived time validities with the multiplication of 100.If the modeled time of thinning was ideal, the ratios should not be greater than 100%.

Field Assessment
The annual BAI among groups were very similar and temporal dynamics were synchronized before thinning based on visual inspection (Figure 3a).After the treatments (Figure 2) had been applied on the forest stands in 1990, the effect on tree growth through time was pronounced.Mean (± standard deviation, sd) basal area of the treatment groups (2-5R, n = 20-21) was very similar, from 712.5 ± 274.6 cm 2 (2R) to 782.7 ± 401.9 cm 2 (3R) (Table 1) for more than two decades after thinning.Basal area of the CTRL group was relatively small (604.9 ± 321.2 cm 2 ) and not different from treated stands (Tukey-Kramer test, p > 0.05).Dynamics of post-thinning relative to growth varied not only according to thinning but also the treatment type (Figure 3b).As expected, thinning enhanced the growth of C. japonica, whereas the growth rate of the CTRL group declined year after year.The temporal patterns of treated groups had a spherical trajectory peaking at 5 (2R)-14 (3R) years after thinning.Mean (± sd) relative growths ranged from 77.9% ± 14.3% for CTRL and from 109.4% ± 5.9% (2R) to 175.1% ± 33.2% (3R) for the treatment groups.Differences between CTRL and 2R (ANOVA, p = 0.04) and 3-5R (p ≤ 0.03) were observed 11 and 5-6 years after thinning, respectively.However, this became negligible (ANOVA, p > 0.05) within thinning groups before the tree-ring sampling.2 for the photographic illustrations.

Event History Analysis
Estimates of the pre-thinning median (time to 50% survival probability by comparing to the threshold, p = 0.974, log-rank statistics) and mean (ANOVA, p = 0.938) survival time derived from the Kaplan-Meyer analysis of the smoothed time-series data were similar for all groups (Table 2, Figure 4a).Differences between treatment groups for the post-thinning periods were quite pronounced (median: p < 0.001, log-rank statistics; mean: ANOVA, p < 0.001) (Table 2, Figure 4b).In the accelerated failure time regression analysis, Weibull and Gaussian distributions were selected to fit pre-and post-thinning time-to-event data, respectively, because these yielded the smallest AIC (626.6 and 687.5, respectively).There was no difference (Table 3) among the pre-thinning groups (Figure 4c) and accelerated growth of thinned groups (above the CTRL curve, Figure 4d).Overall, performance was satisfactory (high fitness, close to a 1:1 relationship and small offset) for the thinned groups before (r 2 = 0.951 (5R)-0.985(2R), slope = 0.998 (5R)-1.121(CTRL), offset = −0.051   2 for the photographic illustrations.

Model Validation
The four polynomial regression models were derived to fit the post-thinning (1991-2011) growth trends of different treatments (2-5R): y 2R = 14.568 + 0.7919 x -0.0509 x 2 + 0.0009 x 3 (4) y 4R = 9.7535 + 5.1983 x -0.4215 x 2 + 0.0100 y 5R = 10.6000+ 3.6923 x -0.2061 x 2 + 0.0033 x 3 (7) where y, x, and the subscripts of y are the BAI, year, and treatments, respectively.Performances of the models were satisfactory (R 2 = 0.80-0.96,p < 0.0001).The simulated results revealed that none of the thinning cycles (5-25 years) can facilitate the growth of C. japonica plantation more than the proposed approach within 100 years since no relative percentages went above 100% (Figure 6).Given the base condition (forest cleared in 1966 and thinned in 1990), implementation of event history analysis can enhance the relative BAI for a C. japonica tree by, on average, 20.9% more in the simulation.

Temporal Effects of Thinning
Thinning enlarges canopy gaps, which increases the penetration of solar radiation to the forest floor.It may also increase the amount of water that reaches the ground.Additionally, more space is provided for the expansion of tree crowns.Together, these would ameliorate the physical environment for plant growth [40][41][42].We studied the temporal effects of thinning in a 45-year-old C. japonica plantation in the mountainous region of northern Taiwan.According to our dendroecological analysis, the growth rate of C. japonica declined ~24 years after regeneration in 1966.The trends of three years of smoothing growth of the BAI showed that it would reduce gradually if there were no thinning.The effects of thinning were pronounced, especially with larger gaps (3-5R, Figure 3).However, only the trees on edges received full benefit from thinning and trees were much less affected without the direct contact of opened space [43].When a wider space is created, more resources (e.g., sunlight, nutrients, and water) are available for the dominant and predominant trees to enhance growth, especially those near the edge [32,44].In time, canopy closure and understory vegetation coverage increased, and the accelerated growth effects of thinning were gradually reduced (Figures 3b and 4).In the future, the growth of the interior trees should also be further assessed to better understand the effects of row thinning on the stand scale for the region.

Temporal Effects of Thinning
Thinning enlarges canopy gaps, which increases the penetration of solar radiation to the forest floor.It may also increase the amount of water that reaches the ground.Additionally, more space is provided for the expansion of tree crowns.Together, these would ameliorate the physical environment for plant growth [40][41][42].We studied the temporal effects of thinning in a 45-year-old C. japonica plantation in the mountainous region of northern Taiwan.According to our dendroecological analysis, the growth rate of C. japonica declined ~24 years after regeneration in 1966.The trends of three years of smoothing growth of the BAI showed that it would reduce gradually if there were no thinning.The effects of thinning were pronounced, especially with larger gaps (3-5R, Figure 3).However, only the trees on edges received full benefit from thinning and trees were much less affected without the direct contact of opened space [43].When a wider space is created, more resources (e.g., sunlight, nutrients, and water) are available for the dominant and predominant trees to enhance growth, especially those near the edge [32,44].In time, canopy closure and understory vegetation coverage increased, and the accelerated growth effects of thinning were gradually reduced (Figures 3b and 4).In the future, the growth of the interior trees should also be further assessed to better understand the effects of row thinning on the stand scale for the region.

Event History Analysis
We propose a novel approach to assess the timing in thinning-induced accelerated tree-ring growth of C. japonica, which has rarely been implemented in forest thinning assessment.The results showed that the forest stand reached the point of thinning at approximately 16 years after the initial clear-cut was applied in 1966 (Table 2).This indicates that plant competition for resources was conspicuous during the observation period, which reached a point where thinning was necessary.This finding coincides with typical pre-commercial thinning of 10-20-year-old stands [45][46][47] and should be delayed to about 15 years for trees on sites with poor physical conditions [48].Event history analysis provides more precise timing for management practices.After the thinning in 1990, the median values (survival probability = 0.5) derived from the dendroecological records using the Kaplan-Meyer method for 2R were relatively shorter than those for other treatments (11 years vs. 14-15 years, Figure 4b).As expected, the thinning time for the CTRL group was insignificant because 16 years had gone by after the clear-cut (Table 2).Thinning creates space that will be filled by lateral crown growth along edges.According to field observations in 2015, the canopy opening produced with the thinning might have been close to the stage of canopy closure in 2R, but not the other treatment groups (Figure 2), which may explain the shorter temporal validities for this treatment (Table 2 and Figure 3b), and imply that gaps wider than 3R might yield similar production (Figure 4b,d).
This study demonstrates the feasibility of accelerated failure time regression for assessing the effects (acceleration of tree growth) of thinning (Table 3 and Figure 4d) on a C. japonica stand in a mountainous region of northeastern Taiwan.The Kaplan-Meyer method can precisely derive temporal validities [19].Use of the accelerated failure time regression allows us to formally compare treatment effects and draw valid inference [18,21].It can also provide a reliable prediction of survival probability of any given point in time following different thinning treatments.We should note that accelerated failure time regression was significantly overestimated in the early stages of CTRL (Figure 5e) because the survival rate was less than 50% at this stage.This is an unavoidable limitation of the proposed approach when the initial probability is not equal to 1. Notwithstanding, this study shows that the Kaplan-Meyer method coupled with accelerated failure time regression curve fitting may provide a precise and robust estimate of ideal forest thinning temporal intervals, which was verified by model validation (Figure 6).The data show that the performance, mean annual BAI, which could be synonymous with annual carbon accumulation [49] for individual trees, is superior to those of shorter rotations with higher operational costs.

Potential for Future Research
The timing of thinning is key for determining silvicultural management, although it is very challenging to determine with conventional approaches [32,44,50,51].We demonstrate that it is possible to determine the exact ideal thinning time under different practices by analyzing dendroecological data and using event history analysis.The proposed approach would help make forest management decisions for scheduling thinning operations with high flexibility since the threshold (in our case, the average BAI of 1988-1990) can be changed to meet certain management plans.Additionally, it can be used to assess the benefits of different thinning methods (e.g., single, row, or block) [43,51] for forest carbon sequestration.Event history analysis can be further developed to better understand natural and anthropogenic alteration-induced canopy gap dynamics and to decipher the complex interactions among each tree by investigating the temporal variation of survival/mortality probabilities at different tree densities.

Conclusions
We used event history analysis to determine the ideal time for row thinning for a conifer planation in the mountainous region of Taiwan.The proposed method only requires tree-ring data to determining the ideal thinning time.We found that there was a positive trend between strip-widths and the timings of thinning.According to the model validation, event history analysis -derived temporal intervals of forest thinning may optimize carbon sequestration.

Figure 1 .
Figure 1.The 2-ha experimental plantation site (encompassed by a white-dashed polygon) near an unpaved road (solid white line) (a) on Chilan Mountain in northern Taiwan (b).

Figure 1 .
Figure 1.The 2-ha experimental plantation site (encompassed by a white-dashed polygon) near an unpaved road (solid white line) (a) on Chilan Mountain in northern Taiwan (b).

Figure 2 .
Figure 2. (a) Thinning treatment schematic depicting 2-5 row (R) cuttings with 2-m row widths.Six rows of trees were retained on either side of the cut strips (b-e).The selection of sample was based on the size distribution of diameter at breast height (DBH) for each group.Photographs show current conditions of row thinning from 2-5R, respectively, taken by C. Chung on April 16, 2015.

Figure 2 .
Figure 2. (a) Thinning treatment schematic depicting 2-5 row (R) cuttings with 2-m row widths.Six rows of trees were retained on either side of the cut strips (b-e).The selection of sample was based on the size distribution of diameter at breast height (DBH) for each group.Photographs show current conditions of row thinning from 2 to 5R, respectively, taken by C. Chung on 16 April 2015.

Forests 2017, 8 , 77 6 of 14 Figure 3 .
Figure 3. (a) The temporal dynamics of annual basal area increment (BAI) under different thinning treatments (2-5R and control plots (CTRL), Figure 2) from 1966-2011.The data were smoothed using a 3-year moving average window.(b) Post-treatment relative growth compared to average growth from 1988-1990.The dotted gray lines in (a) show the pre-(the mean BAI of 1966-1990) and post-thinning (the mean BAI of 1988-1990) thresholds, respectively.The dashed gray line in (b) indicates the 100% relative growth.

Figure 3 .
Figure 3. (a) The temporal dynamics of annual basal area increment (BAI) under different thinning treatments (2-5R and control plots (CTRL), Figure 2) from 1966 to 2011.The data were smoothed using a 3-year moving average window.(b) Post-treatment relative growth compared to average growth from 1988 to 1990.The dotted gray lines in (a) show the pre-(the mean BAI of 1966-1990) and post-thinning (the mean BAI of 1988-1990) thresholds, respectively.The dashed gray line in (b) indicates the 100% relative growth.

Figure 4 .
Figure 4.Estimated survival curves for control and thinning treatments using the Kaplan-Meyer method and an accelerated failure time regression model for pre-thinning (1966-1990) (a and c) and post-thinning (1991-2011) periods (b and d), respectively.

Figure 4 .
Figure 4.Estimated survival curves for control and thinning treatments using the Kaplan-Meyer method and an accelerated failure time regression model for pre-thinning (1966-1990) (a and c) and post-thinning (1991-2011) periods (b and d), respectively.

Figure 5 .
Figure 5. Accelerated failure time regression model probability estimation for pre-and post-thinning conditions based on the Kaplan-Meyer method for different treatments: (a) 2R, (b) 3R, (c) 4R, (d) 5R, and (e) CTRL.The gray lines indicate the 1:1 relationship.

Figure 5 .
Figure 5. Accelerated failure time regression model probability estimation for pre-and post-thinning conditions based on the Kaplan-Meyer method for different treatments: (a) 2R, (b) 3R, (c) 4R, (d) 5R, and (e) CTRL.The gray lines indicate the 1:1 relationship.

Figure 6 .
Figure 6.Model validation of event history analysis application on forest thinning by comparing a simulated 100-year mean annual basal area increment under different thinning cycles (5-25 years) in different treatments (2-5R) with those derived from event history analysis.The gray dashed line indicates relative percentage out of 100%.

Figure 6 .
Figure 6.Model validation of event history analysis application on forest thinning by comparing a simulated 100-year mean annual basal area increment under different thinning cycles (5-25 years) in different treatments (2-5R) with those derived from event history analysis.The gray dashed line indicates relative percentage out of 100%.

Table 1 .
The summary of stand characteristics and sampled tree measured in 2011.

Table 1 .
The summary of stand characteristics and sampled tree measured in 2011.
1Please see Figure2for the photographic illustrations.

Table 3 .
Accelerated failure time regression models for pre-and post-(values in parentheses) thinning periods (n = 104).