A Nonlinear Mixed-Effects Height-to-Diameter Ratio Model for Several Tree Species Based on Czech National Forest Inventory Data

Abstract

Height-to-diameter at breast height (DBH) ratio (HDR) is an important tree and stand stability measure. Several factors such as stand dynamics, natural and anthropogenic disturbances, and silvicultural tending significantly affect HDR, and, therefore, in-depth investigation of HDR is essential for better understanding of ecological processes in a forest. A nonlinear mixed-effects HDR model applicable to several tree species was developed using the Czech national forest inventory data comprising 13,875 sample plots and 348,980 trees. The predictive performance of this model was evaluated using the independent dataset which was originated from 25,146 trees on 220 research sample plots. Among various tree- and stand-level variables describing tree size, site quality, stand development stage, stand density, inter-tree spacing, and competition evaluated, dominant height (HDOM), dominant diameter (DDOM), relative spacing index (RS), and DBH-to-quadratic mean DBH ratio (dq) were identified as the most important predictors of HDR variations. A random component describing sample plot-specific HDR variations was included through mixed-effects modelling, and dummy variables describing species-specific HDR variations and canopy layer-specific HDR variations were also included into the HDR model through dummy variable modelling. The mixed-effects HDR model explained 79% of HDR variations without any significant trends in the residuals. Simulation results showed that HDR for each canopy layer increased with increasing site quality and stand development stage (increased HDOM) and increasing competition (increased RS, decreased DDOM and dq). Testing the HDR model on the independent data revealed that more than 85% of HDR variations were described for each individual species (Norway spruce, Scots pine, European larch, and European beech) and group of species (fir species, oak species, birch and alder species) without significant trends in the prediction errors. The HDR can be predicted with a higher accuracy using the calibrated mixed-effects HDR model from measurements of its predictors that can be obtained from routine forest inventories. To improve the prediction accuracy, a model needs to be calibrated with the random effects estimated using one to four randomly selected trees of a particular species or group of species depending on the availability of their numbers per sample plot. The HDR model can be applied for stand stability assessment and stand density regulation. The HDR information is also useful for designing a stand density management diagram. Brief implications of the HDR model for designing silviculture strategies and forest management planning are presented in the article.

1. Introduction

Height-to-diameter ratio (HDR) is total height divided either by diameter at breast height or by diameter at the root collar of a tree. The HDR is a tree-level index of slenderness and is used for the evaluation of tree and stand stability. Large values of HDR indicate that either trees have grown in a crowded stand with mutual support from neighboring trees or they have grown in an extremely open stand where no significant competition exists [1,2,3,4]. Smaller values of HDR indicate longer crown length, higher crown projection area, better developed root system, lower position of the center of gravity, and higher stability of the trees. Many studies involving stand stability assessments have used HDR for characterizing tree and stand stability and their vulnerability to natural disasters [3,5,6,7,8,9]. Trees with larger HDR are always associated with increased risk of vulnerability to uprooting by windstorm and breakage by windstorm and icing [3,4,10,11,12,13]. Trees with large HDR may be more vulnerable because their stems may not be adapted to a condition of higher mechanical perturbation [3,4,14,15]. The mechanical properties of stem wood may be evaluated using HDR, for example, trees with small HDR may have higher bending movement than trees with larger HDR of similar heights [15,16,17,18].

Various HDR applications are possible for the effective silvicultural tending and forest management. By understanding the degree to which trees and stands are more susceptible to snow, icing, and wind damages, forest managers may better design silviculture prescriptions based on the range of HDR [1,3,5,9]. The HDR may be used for the evaluation of effectiveness and efficiency of thinning, as thinning substantially influences HDR, both mean stems, and even the upper layer of the trees in a stand [19,20]. Furthermore, HDR can be used as a prominent predictor to describe competition effects in various forest models [19,21,22,23,24,25,26]. The HDR is used for assessing tree vigor and health [19,27], for example, a tree with a higher HDR may have a lower vigor in general. The HDR derived from any empirical HDR model may be considered as a benchmark value that may be compared against HDR derived from tree height and diameter increment models, and models of height-diameter relationships. Since HDR is significantly correlated to crown ratio of the trees, trees with unfavorable properties in both HDR and crown ratio can be removed [1,20,28].

Tree and stand characteristics change over time, and thus, various empirical forest models including HDR models are necessary to update the information of these characteristics. Updating is done either from fresh measurements of HDR for all trees on each sample plot or empirical HDR models are used to predict HDR. In general, diameter and height increment models or height–diameter models are used to update the HDR instead of using HDR prediction models. However, the former models are often based on different datasets and the parameters of these models might not have been estimated using the simultaneous regression; therefore, HDR derived from these models is significantly biased [29]. Small bias associated with these models may lead to a large bias in the derived HDR. If a diameter increment model over-predicts and a height increment model under-predicts or vice versa, bias associated with HDR derived from these models could be substantial. Therefore, a species-specific HDR model is necessary for more precise prediction of HDR at tree-level. This can be possible with the HDR model developed using the comprehensive data that cover a wide variability in terms of tree- and stand-level characteristics (tree size and vigor, site quality, stand development stage, stand age, stand density, stand structure, silvicultural tending, and management regime). The methods employed to estimate model parameters also determine the degree of accuracy of the developed model.

The HDR information is essential for better understanding forest ecological processes in any forest as HDR is substantially influenced by tree species, stand age, stand structure, stand density, site quality (soil type, soil moisture, litter depth, slope, elevation, and exposition), stand development stage, climatic and natural disturbance (light, wind, snow, and icing), silvicultural tending, and species provenance [14,20,30,31,32,33,34,35,36,37,38,39]. Stand location, such as a stand on the forest edge, a stand in large gap or an unstable stand fragment also influence HDR [40,41,42]. The HDR varies with inter-tree spacing, even for a single species within the same stand, and extreme HDR may reach the extremely spare stands and the crowded stands [2,9]. The HDR also significantly varies on the trees with top or lower canopy layers within a stand [43]. The HDR also varies with tree root system, crown width and crown length [9]. Among the above-mentioned factors, inter-tree spacing, competition and stand density have the biggest influences on HDR [2,14,40,42,44,45,46,47]. Thus, the inclusion of variables describing the effects inter-tree spacing and competition, stand density, site quality and stand development stage into the HDR model is necessary to increase its prediction accuracy. A composite model developed using a large dataset that includes the interactive effects of the above-mentioned factors affecting HDR could provide a significantly high prediction accuracy.

The prediction accuracy of a composite HDR model can be further increased through integration of an unstructured random component by applying the mixed-effects modelling, which describes HDR variabilities caused by various stochastic factors. The hierarchically-structured dataset often leads to dependency among the measurements, as observations from the same sampling unit are most likely to be correlated. The assumptions of independent errors are largely violated while ordinary least square regression is used for estimating models [48]. Consequently, estimated parameters and variances would be significantly biased, resulting in the invalidation of the hypothesis tests [49,50]. The mixed-effects modelling allows for inclusion of all types of variabilities into the model and significantly increases the prediction accuracy [48,51,52,53,54,55].

To develop the robust mixed-effects HDR model (model with stable parameter and variance estimates), a large dataset is required, which can be supplied by a national forest inventory (NFI) program. This dataset covers a wide variability of growth conditions and silviculture tending on the forest stands across the country. In recent decades, several studies have used NFI data to develop forest models in many European countries [25,56,57,58,59,60,61,62,63]. Developing models using NFI data and testing against research sample plot data has been frequently practiced [25,58,63]. Testing the models against research sample plot data would increase the credibility of the models. A large dataset allows the interactive effects of various factors that affect HDR to be integrated into the model using complex statistical analyses.

Several studies have been carried out to show the relationships between HDR and the vulnerability of trees and stands; however, only a few have developed the statistical HDR models [2,14,43]. All existing HDR studies were based on smaller datasets acquired from few species compared to the datasets used in the current study, and HDR models were developed with a conventional modelling approach, i.e., least square regression approach. HDR modelling may be based on the data acquired from several tree species in the monospecific and mixed stands. Developing separate species-specific models for several species and their separate applications would be tedious and effort demanding. To increase the application efficiency of the model, this study applied the dummy variable modelling approach, which allows making a single HDR model applicable to several tree species. Therefore, this study aims to (1) develop a nonlinear mixed-effects HDR model (composite model) using a large NFI dataset through integration of the predictor variables describing the effects of inter-tree spacing and competition, site quality and stand development stage, canopy height differentiation and social position of individual trees, and unstructured component describing sample plot-level HDR random variations; and (2) evaluate the HDR model’s prediction performance using a large independent dataset acquired from research sample plots distributed in various parts of the Czech Republic. The parameters of a nonlinear mixed-effects HDR model are estimated using maximum likelihood procedure. The presented HDR model, which requires only a few predictors that can be derived from a routine forest inventory database, may serve as a useful tool in assessing tree and stand stability. This model will be used for various forest management implications that are briefly discussed in the article. Since conifer tree species are much more susceptible to ice- and wind-related damages than broadleaved species, the presented HDR model may be more useful for conifer species.

2. Materials and Methods

2.1. Data Materials

Our study is based on two large datasets: One originated from the Czech NFI program (model fitting dataset, hereafter termed as a training dataset) and another originated from permanent research sample plots (validation dataset). These datasets were largely different from each other in terms of sampling design and measurement method used, measurement accuracy, sample plot coverage in the tree population, and statistical characteristics. A brief description of each dataset is presented below.

2.2. Training Dataset

The Czech NFI program provided the training data, which were used to fit the model and were collected from sample plots distributed across the forested parts of the country (Figure 1a). The NFI sample plots were laid out following the protocols devised by Field-Map technology of the Institute of Forest Ecosystem Research (IFER) Monitoring and Mapping Solutions Ltd. [64]. Measurements were carried out for all trees with diameter at breast height (DBH, 1.3 m above-ground) ≥7 cm on each sample plot. We used measurements from two inventory cycles, i.e., the first measurement was done between 2001 and 2004 and the second measurement between 2011 and 2015. All measurements in both inventory cycles comprise of 13,875 sample plots and 348,980 trees, representing a large variability in site quality, stand development stage, stand density, species mixture, stand structure, silvicultural tending and other management interventions, and natural disturbances. Tree- and stand-level variables were measured by considering the inventory protocols of the Forest Management Institute [65]. Details of sampling design, NFI sample plot establishment and measurement methods are found in the literature [59,66,67,68]. The training dataset comprised of 30% monospecific sample plots and 70% mixed species plots. The definition of the monospecific sample plots considered the inclusion of all individuals other than a species of interest if they had over-bark DBH <4 cm [43]. In terms of species composition, coniferous tree species predominate with a share of 58.9% (Norway spruce 44.1%) and broadleaved tree species occupy a share of 41.1% (European beech 10.3%) [59]. The number of trees by species in this dataset is given in Table S1.

2.3. Validation Dataset

This dataset was obtained from the measurements of 220 research sample plots. These plots, hereafter, are termed as sample plots. Sample plots are square or rectangle, and size varies from 2500 m² to 4900 m². Sample plots are laid out in different parts of the Czech Republic (Figure 1b). Several stand characteristics (stand structure, natural regeneration, dead wood stock, and stand development stage and site quality) were considered as main criteria while laying out sample plots. Sample plot locations vary from 240 m to 1,370 m a.s.l., where mean climate characteristics (i.e., annual temperature (4–9.5° C), mean annual precipitation (500–1450 mm), and growing season length (45–180 days) also largely vary. In-depth descriptions of this dataset are also available in the literature [43,54,59,69,70,71,72,73,74]. This dataset comprises 24% monospecific sample plots and 76% mixed species sample plots. The definition of the monospecific sample plots considered the inclusion of all individuals other than a species of interest if they had over-bark DBH <4 cm [43]. Measurements were carried out between April 2007 and March 2016. However, no repeated measurement was involved, meaning that there was no temporal variation in this dataset. Total height and DBH were measured with precisions of 0.1 m and 1 mm, respectively. The number of trees by species in the validation data is given in Table S2. Graphs of HDR plotted against DBH in both training and validation datasets are presented in Figure 2.

2.4. Data Analysis

2.4.1. Deriving Stand-Level Variables

We derived several stand-level variables, which describe site quality and stand development stage, stand density and competition, and inter-tree spacing, and they were assessed for their contributions to the HDR variations. We then selected easily accessible variables that contributed significantly to the HDR model. Even though site quality is precisely described by site index, we could not use this as our data lacked this information. Instead, we included dominant height (HDOM) into the HDR model. We identified dominant trees and computed dominant height for each sample plot by applying the same methods as in Sharma et al. [63]. The HDOM has been frequently used to describe the combined effect of site quality and stand development stage in modelling various tree characteristics [43,51,52,54,55,58,70,71,75]. For mixed species sample plots, where the number of dominant trees of a particular species or group of species was inadequate according to the definition of dominant height (i.e., average of height of 250 largest trees per hectare), measured heights of the largest trees per sample plot regardless of tree species were used to calculate HDOM [59]. We also calculated the dominant diameter (DDOM) for each sample plot.

We evaluated the potential contributions of other stand-level measures to the HDR model such as stem per hectare (N, ha⁻¹), stand basal area (BA, m² ha⁻¹), arithmetic mean DBH (AMD, cm), and quadratic mean DBH (QMD, cm) per sample plot. We also evaluated two tree-centered competition measures: DBH-to-QMD ratio (dq) and basal area of the trees larger in diameters than a subject tree (BAL, m² ha⁻¹) for their potential contributions to the HDR model. We evaluated the relative spacing index (RS), which is defined by $\mathrm{RS}=\sqrt{10000/N}/\mathrm{HDOM}$ [55,76,77] for its potential contribution to the HDR model. The RS is assumed as a more effective density measure than others, because it incorporates number of stems per hectare, site quality and stand development through HDOM. All measures were calculated using the dendrometric measurements of all trees regardless of species per sample plot. Summary statistics of important variables in both training and validation datasets are presented in

Table 1. Summary statistics of data.

Variables	Statistics (Mean ± Standard Deviation (Range))
	Training Data	Validation Data
Number of sample plots	13,875	220
Number of HDR sample trees	348,980	25,146
Number of HDR sample trees per sample plot	36.5 ± 15.9 (4–90)	232.3 ± 167.8 (8–664)
Number of stems (N ha−1)	1514 ± 671 (40–5460)	940 ± 711 (32–4700)
Stand basal area (BA, m2 ha−1)	38.3 ± 13.5 (0.07–85.4)	41.3 ± 14.4 (0.1–81.1)
BA of trees lager than a subject tree (BAL, m2 ha−1)	28.5 ± 14.4 (0–83.3)	32.2 ± 17.6 (0–79.5)
Quadratic mean DBH per sample plot (QMD, cm)	27.3 ± 6.9 (8.7–84.2)	28.6 ± 10.3 (9.6–87.3)
DBH-to-QMD ratio (dq)	0.95 ± 0.31 (0.13–4.5)	0.86 ± 0.52 (0.04–7)
Arithmetic mean DBH per sample plot (AMD, cm)	25.8 ± 6.6 (8.5–78.4)	24.9 ± 10.6 (9.1–84.4)
Dominant height per sample plot (HDOM, m)	23.2 ± 6.2 (4.4–42.4)	27.5 ± 6.9 (8–42.8)
Dominant diameter (DDOM, cm)	31.5 ± 8.1 (8.6–78.4)	49.7 ± 14.2 (13.3–84.4)
Total height (H, m)	21.1 ± 7.1 (1.5–54.3)	17.2 ± 9.3 (1.4–50.6)
Diameter at breast height (DBH, cm)	26.4 ± 11.6 (7–117.3)	25.2 ± 17.4 (2–118.1)
Height-to-DBH ratio (HDR, m cm−1)	0.85 ± 0.19 (0.08–2.15)	0.81 ± 0.28 (0.11–2.3)

. The number of HDR sample trees per sample plot in the validation dataset is higher compared to that in the training dataset, because of larger sample plot size (2500–4900 m²) in the former dataset compared to the latter plot size (500 m²). However, the number of stems per hectare in the training dataset is higher compared to that in the validation dataset. It is because that training dataset includes all age and DBH classes of the forest stands, whereas the majority of the validation data came from middle-aged and mature forest stands.

2.4.2. Model Construction

Considering the nonlinear pattern on the HDR plotted against DBH for each species or group of species (Figure 2), which is an exponential decrease of HDR with increasing DBH, we used a simple exponential decay function, hereafter termed as a base function (Equation (1)), to fit data. This base function was then expanded through the introduction of additional predictors that describe the effects of stand density, inter-tree spacing, competition, site quality and stand development stage, and sample plot-specific random variability.

$${\mathrm{HDR}}_{ij}= b_{1}exp\left(-b_2{\mathrm{DBH}}_{ij}^{b_3} \right)+{\epsilon }_{ij}$$

(1)

where HDR_ij and DBH_ij are the height-to-diameter ratio (m cm⁻¹) and diameter at breast height of the j^th tree (cm) on the i^th sample plot, respectively; b₃ = 0.5; b₁ and b₂ are the parameters to be estimated; and ${\epsilon }_{ij}$ is an error term.

All stand-level measures and tree-centered competition measures (Table 1) were evaluated for their contributions to the HDR variations. The two-stage variable selection approach [78,79] was used to identify significantly contributing variables to the HDR model. This involved two steps: Firstly, the base function was fitted to the pooled sample-plot data from a particular species or group of species by sample plot for all sample plots. Secondly, the matrix-scattered graphs of the sample plot-specific parameter estimates (i.e., estimates of b₁, b₂) plotted against each variable and its transformation and interactions with other predictors were also examined, and those showing stronger correlations with b₁ or b₂ were identified. Because of its biological logic and interpretation, the two-stage approach has frequently been used for modelling tree characteristics [53,70,80,81]. Among various predictors evaluated using the two-stage approach, HDOM, DDOM, RS, and dq showed stronger relationships only with b₁ of Equation (1). Following the principles of modelling categorical variables [62,82], we formed two dummy variables (CC₁, CC₂) to account for HDR variations caused by three canopy layers denoted by CHC, as below:

Canopy Height Class	CC1	CC2
CHC1	0	0
CHC2	1	0
CHC3	0	1

The inventory crew did not record CHC information (i.e., suppressed, intermediate and dominant) in the inventory. However, we formed CHC on each sample plot in our data by assuming CHC1: height >66% height of the tallest tree (dominant canopy layer); CHC2: 33% height of the tallest tree < height < 66% height of the tallest tree (intermediate canopy layer); and CHC3: height <33% height of the tallest tree (suppressed canopy layer). This classification was possible due to measurements of total height of all trees regardless of species on each sample plot. This classification was also used in other studies [43,71].

For simplicity (i.e., easy fitting and efficient application of HDR model), only four individual species and three species group were considered (Table S1). We formed six dummy variables (TS₁, TS₂, TS₃, TS₄, TS₅, TS₆) for seven species to describe species-specific or group-specic HDR variaitons, as below:

Tree Species	TS1	TS2	TS3	TS4	TS5	TS6
Norway spruce	0	0	0	0	0	0
Scots pine	1	0	0	0	0	0
European larch	0	1	0	0	0	0
Fir species	0	0	1	0	0	0
European beech	0	0	0	1	0	0
Oak species	0	0	0	0	1	0
Birch and alder speices	0	0	0	0	0	1

Four measurable variables (HDOM, DDOM, RS, dq) and eight dummy variables (CC₁, CC₂, TS₁, TS₂, TS₃, TS₄, TS₅, TS₆) were included into the base function (Equation (1)) through expression of b₁ as a function of all these variables as in Equation (2). This expression was chosen as this provided the convergence with the smallest residual variations (smallest sum of squared errors) in fitting the model.

$$b_1=f\left(\mathrm{HDOM}, \mathrm{DDOM}, \mathrm{RS},\mathrm{dq}, C{C}_1, C{C}_2, T{S}_1, T{S}_2, T{S}_3, T{S}_4, T{S}_5, T{S}_6\right)$$

(2)

where HDOM = dominant height (m); DDOM = dominant diameter (cm); RS = relative spacing index; dq = DBH-to-quadratic mean DBH ratio; and dummy variables are the same as defined above.

Our database contained hierarchical sampling units and, therefore, it was necessary to include variabilities among the sampling units. The mixed-effects HDR model was formulated using Equation (1) (with Equation (2) included) through introduction of a random component to describe the stochastic variabilities of HDR at the sample plot-level [48]. Among the three mixed-effects model alternatives formulated, convergence with a global minimum was achieved only with the random effect parameter (u_i) added to b₁. Details of the mixed-effects model formulation are found in the statistical textbooks [48,83]. We present only a final form of a nonlinear mixed-effects model (Equation (3)) that describes HDR variations with the smallest residual variations in our data.

$${\mathrm{HDR}}_{ij}= \left(b_1+u_i\right)exp\left(-b_2{\mathrm{DBH}}_{ij}^{b_3}\right)+{\epsilon }_{ij}$$

(3)

with $b_1={\alpha }_{1}C{C}_1+{\alpha }_{2}C{C}_2+{\alpha }_3{\mathrm{HDOM}}_i^{{\alpha }_4}+{\alpha }_5{\mathrm{DDOM}}_i+{\alpha }_6{\mathrm{RS}}_i+{\alpha }_7{\mathrm{dq}}_{ij}+{\alpha }_{8}T{S}_1+{\alpha }_{9}T{S}_2+{\alpha }_{10}T{S}_3+{\alpha }_{11}T{S}_4+{\alpha }_{12}T{S}_5+{\alpha }_{13}T{S}_6$, where CC_k is a dummy variable (k = 1,2); TS_k is dummy variable (k = 1, 2, …, 6); b₃ = 0.5; b₂ and α_k (k = 1,2, …, 13) are parameters to be estimated; and ${\epsilon }_{ij}$ is an error term; the other symbols and abbreviations are the same as in Equations (1) and (2). When CC belongs to k, then CC_k = 1, and otherwise is 0; similarly, when TS belongs to k, then TS_k = 1 and otherwise is 0 (see above two tables of dummy variable formation). When b₃ in Equations (1) and (3) was tried to be estimated along with other parameters by optimization, the convergence with a global minimum was not achieved. We then compared the sum of squared errors (SSE) produced from the fitting with several alternative values of b₃ (0.1 to 2 by 0.1 increment) iteratively, and b₃ = 0.5 resulted in the smallest SSE. This is the reason that b₃ = 0.5 was chosen in both Equations (1) and (3).

In Equation (3), the vectors of errors (${\epsilon }_{ij}$) and sample plot-level random effects (u_i) are defined by ${\mathit{\epsilon }}_{\mathit{i}}~\mathit{N}\left(0,\mathit{R}\right)$ and ${\mathit{u}}_{\mathit{i}}~\mathit{N}\left(0, \mathit{D}\right)$ respectively. It was assumed that error vector ${\mathit{\epsilon }}_{\mathit{i}}$ would have a normal distribution with zero mean and within-plot variance–covariance matrix ${\mathit{R}}_{\mathit{i}}$ defined by:

$${\mathit{R}}_{\mathit{i}}={\sigma }^2{\mathit{G}}_{\mathit{i}}^{0.5}{\mathsf{\Gamma }}_{\mathit{i}}{\mathit{G}}_{\mathit{i}}^{0.5}$$

(4)

where σ² is a scaling factor and is equivalent to a residual variance of the estimated HDR model and common to all plots [50]. A matrix ${\mathsf{\Gamma }}_{\mathit{i}}$ accounts for within-plot residual autocorrelations. However, we assumed this matrix as an identity matrix, as our data lacked significant autocorrelations. The matrix G_i is a diagonal matrix that accounted for a within-sample plot heteroscedasticity. The random effect vector u_i was assumed to have a normal distribution with zero mean and plot variance–covariance matrix D defined by:

$$\mathit{D}={\sigma }_{ui}^2$$

(5)

2.4.3. Estimating Model Parameters and Evaluation

Parameters of the mixed-effects HDR model (Equation (3)) were estimated using the SAS macro NLINMIX version 9.4 [84] with the expansion-around-zero method [85] and restricted maximum likelihood method. This is a linearization method of Sheiner and Beal [86], which expands the likelihood function with a first-order Taylor series approximation. The fitting performance of the model was evaluated using root mean square error (RMSE), adjusted coefficient of determination (R²_adj), and Akaike information criterion (AIC). The AIC is based on minimizing the Kull-back-Lieber distance and it imposes a penalty for the number of parameters in the model [87,88]. Formulae of all these statistical measures are found in the statistical textbooks [89]. The graphs of residuals, prediction errors and simulated HDR curves produced with fitted HDR model were also examined. Unless otherwise specified, we used 1% level of significance in our analyses. Realizing the fact that validation with an external independent dataset increases the credibility of the model, we tested our HDR model using validation data. Sample plot-specific HDR curves generated using the calibrated mixed-effects HDR model after adjustment of the random effects predicted with the empirical best linear unbiased prediction (EBLUP) theory [48] were examined for all sample plots in the validation data.

2.4.4. Calibrating Mixed-Effects Model and Predicting Sample Plot-Specific HDR

Predicting the response variable (HDR in our case) with the random effects estimated from its prior information is termed as calibration (localization) of the mixed-effects model [48,90]. Prediction of the random effect and its adjustment to the fixed part of the mixed-effects HDR model might be possible using prior information of HDR of any number of trees per sample plot. We calibrated the mixed-effects HDR model using the random effects predicted from one to four randomly selected trees depending on the available number of trees of a particular species or group of species per sample plot. This selection strategy was employed because many sample plots had less than four trees of a particular species or group of species in the validation data [59]. The EBLUP theory (Equation (6)) was applied to predict the random effect for each sample plot using the procedure of calculating matrices in the PROC IML in SAS version 9.4 [84].

$$u_i=\mathit{D}{\mathit{Z}}_{\mathit{i}}^{\mathit{T}}{\left({\mathit{Z}}_{\mathit{i}}\mathit{D}{\mathit{Z}}_{\mathit{i}}^{\mathit{T}}+{\mathit{R}}_{\mathit{i}}\right)}^{-1}{\mathit{\epsilon }}_{\mathit{i}}$$

(6)

where u_i is a random effect vector that describes plot-level HDR variations for the i^th plot. Equations (4) and (5) provide the elements of matrices ${\mathit{R}}_{\mathit{i}}$ and $\mathit{D}$ respectively. Vector ${\mathit{\epsilon }}_{\mathit{i}}$ of the prediction errors was derived from the mean model (only fixed part of the mixed-effects HDR model). The element (z_ij) of the designed matrix Z_i for each tree and sample plot was determined from partial derivatives of a nonlinear mixed-effects model (Equation (3)) with respect to its fixed-effect parameter (b₁) to which random effect (u_i) has been added [48,53,58,91]. The following equation was obtained for this purpose.

$$z_{ij }\left(b_1\right)=exp\left(-b_2{\mathrm{DBH}}_{ij}^{b_3}\right)$$

(7)

where abbreviations and parameters are the same as in Equation (1).

3. Results

The fixed-effect parameter estimates of all predictors including dummy variables are highly significant (

Table 2. Parameter estimates, variance components and fit statistics of the mixed-effects HDR model (Equation (3)).

	Estimate	Standard Error	t-Value	Pr > |t|
Fixed
α1	−0.4005	0.001765	−226.98	<0.0001
α2	−1.0526	0.004407	−238.81	<0.0001
α3	0.7998	0.004939	161.95	<0.0001
α4	0.4568	0.002249	203.06	<0.0001
α5	−0.02188	0.000223	−97.98	<0.0001
α6	0.1236	0.01565	7.90	<0.0001
α7	−0.4599	0.002771	−165.97	<0.0001
α8	0.01049	0.001457	7.20	<0.0001
α9	0.07762	0.002047	37.91	<0.0001
α10	−0.04291	0.003544	−12.11	<0.0001
α11	0.01252	0.002087	6.00	<0.0001
α12	−0.04845	0.002019	−23.99	<0.0001
α13	0.04806	0.001804	26.65	<0.0001
b2	0.1823	0.000747	244.09	<0.0001
Variance
σ2ui	0.01109
σ2	0.006923
Fit statistics
RMSE	0.0923
R2adj	0.7863
AIC	−729778

R²_adj: adjusted coefficient of determination; RMSE: root mean squared errors; AIC: Akaike’s information criterion; α₁, …, α₁₃, b₂: fixed-effect parameters; u_i = random effect for i^th sample plot; σ²: residual variance; other symbols are the same as in Equations (1)–(3).

). The estimated values and signs of all parameters are biologically plausible and interpretable. The model described most parts of the HDR variations (R²_adj = 0.7863) in the training data without significant trends in the residuals.

A box plot of the residuals plotted against each of the potential predictors (Table 2) was examined for each canopy height class; however, for the brevity of space, only residuals plotted against DBH, height and estimated HDR are presented here (Figure 3). In each tree species or group of species, there was no significant trend in the residuals across the observed ranges of the predictors and fitted values of HDR. Only a few larger deviations could be seen for some extremely large or small trees. The variance heteroskedasticity was not observed in the residuals either. The Gaussian distribution patterns (bell-shaped patterns) of the residuals indicated the absence of skewness for each individual species or group of species.

The effects of measurable predictors on HDR were simulated using only fixed parts of the mixed-effects HDR model (i.e., u_i = 0) (Figure 4, Figure 5, Figure 6 and Figure 7). The same scale in the y-axis for all these figures was used for better comparison of the effects of predictors on the HDR. A wide range of the observed data of two predictors HDOM and DDOM in each of the canopy height classes in the training data resulted in several HDR curves for each species or group of species (Figure 4 and Figure 5). However, a relatively narrower range of the observed data of two predictors RS and dq in the suppressed canopy heights (CHC = 3) produced fewer curves for some individual species or groups of species (Figure 6 and Figure 7). The combined effects of site quality and stand development stage described by HDOM emerged as the biggest effect, followed by the effect of tree crowding or competition described by DDOM and dq. For a given condition of stand density and inter-tree spacing, HDR increased with increasing site quality and stand development stage (increased HDOM), and increasing competition (decreased DDOM and dq) and increased RS). The magnitudes of the effects by each predictor on HDR for each canopy height class significantly differed.

Our mixed-effects HDR model was tested using the validation data acquired from research sample plots (Figure 1b). Overall, our calibrated mixed-effects HDR model described 94.4% of HDR variations (R² = 0.9443; RMSE = 0.06745) and more than 85% of HDR variations for individual species or groups of species in this dataset, with a description of the largest part for European beech and the smallest part for fir species (

Table 3. Prediction statistics in the validation data.

Species	RMSE	R2
Norway spruce	0.0673	0.8922
Scots pine	0.0657	0.9039
European larch	0.0667	0.8821
Fir species	0.0654	0.8574
European beech	0.0676	0.9605
Oak species	0.0681	0.9176
Birch and alder species	0.0668	0.9328

RMSE: root mean squared errors; R²: coefficient of determination.

). The mixed-effects HDR model was calibrated with the random effect predicted from measured HDR of one to four randomly selected trees of a particular species or group of species depending on the availability of their numbers on each sample plot. The model test confirmed that a single HDR model developed with the pooled dataset from several species in the NFI sample plots was precise enough for predicting the HDR of the corresponding tree species in the research sample plot data. Graphs of the prediction errors of the calibrated mixed-effects HDR model for each species or group of species showed a promising appearance, i.e., absence of significant trends in the prediction errors for each species and group of species (Figure 8). The prediction errors were found falling within the ±0.25 range for almost all sample plots, indicating that the calibrated mixed-effects HDR model was precise for the validation data.

The simulated curves of the calibrated mixed-effects HDR model for all sample plots in the validation data were examined (Figure 9). Except for some sample plots where outlier measurements existed, curves generated with the calibrated mixed-effects HDR model showed complete covering of the measured data for each species or group of species. For mixed species sample plots, curves were clearly differentiated and each curve was found to be passing through the middle of the data points.

4. Discussion

With four covariate predictors (HDOM, DDOM, RS, and dq) and one sample plot-specific random effect parameter included into the HDR model, the fitting improvement achieved was the highest relative to that of the HDR model fitted with DBH as a single predictor. The model describes most of the HDR variations without systematic residual trends (Table 2, Figure 3). The absence of residual trends across the observed range of the predictors for each species or group of species suggests that the chosen base function and covariate predictors and the assumed canopy height classes were all adequate. Fitting improvement could be marginally increased through the inclusion of additional predictors. However, the inclusion of several predictors not only impedes convergence, but also results in biased estimation due to over-parameterization [89]. Introducing several predictors into the model increases inventory cost and reduces the model use efficiency [92,93].

Site quality and stand development stage are important characteristics of a stand and, therefore, forest managers need to assess the stand characteristics (e.g., HDOM) that describe them in relationship with HDR. Thus, including HDOM into the HDR model that describes the combined effect of site quality and development stage would be an appealing approach. For a given stand density, the effects of site quality and stand development stage on the HDR model for each species and group of species emerged the highest, but substantially varied with canopy layers (Figure 4). The stand dynamics and individual tree characteristics are tightly linked to HDOM and yield-capacity of the stands [94]. This is the reason that HDOM is commonly used as a predictor in different forest models such as height-diameter models [58,91,95,96,97], crown models [52,55,70,75,98,99], HDR models [43], height to crown base models [51,54], crown-to-bole diameter ratio models [71], and crown height increment model [94]. Unlike the site index, HDOM cannot be considered as an effective measure of site quality, because the same HDOM may be resulted in by younger stands on better sites and older stands on poorer sites. However, in the absence of site index information, HDOM may be used to explain the combined effects of site quality and stand development stage in various forest models.

The other three covariate predictors related to stand density and competition (DDOM, RS and dq) significantly contributed to the fitting improvement of the HDR model; however, each of their contributions is less than that of HDOM (Figure 4, Figure 5, Figure 6 and Figure 7). For a given site quality, HDR increases with increasing competition as depicted by decreased DDOM and dq (Figure 5 and Figure 7) and increased RS (Figure 6). The negative signs of the parameter estimate of DDOM and dq suggest that HDR increases with decreasing DDOM and dq. This is due to the influence of tree crowding that results in taller height, smaller crown and thinner bole. A tree growing in a very dense stand under the effect of mutual support of its neighbors would have extremely large HDR [2,3,4]. Denser stands not only make the trees become taller, but also become thinner. Diameter growth is more sensitive to the effect of competition than that of the height growth. The studies carried out in different stand densities show that height growth is almost stable within a broad range of tree crowding, whereas growth of the bole diameter and crown size of a subject tree are sensitive to the competitive impact of their neighbors [100].

The RS is also related to stand density and depicts a significant contribution to the HDR variations, but less than that of other predictors. The positive sign of the parameter estimates of RS (Table 2) indicate that HDR increases with increasing RS (Figure 6), which is biologically logical and interpretable. For a given site quality and development stage of a stand, RS increases with decreasing stem number per hectare and consequently HDR increases, because trees attain a higher rate of height growth than that of DBH under the normal stand density level. However, this trend might not exist in very sparse stands or trees grown in extremely open stands, because trees have attained a higher diameter growth rate relative to that of the height growth to increase their strength against windstorm and ice load. Compared to the suppressed trees, trees with dominant canopy positions can allocate more resources to diameter growth relative to height growth, which may result in smaller HDR and are more stable against the external forces. The inter-tree spacing determines crown shape, crown size, height-to-crown base, and increments of height and diameter attained by trees [55,60,77,101]. The RS incorporates the number of stems per hectare, site quality and stand development stage through HDOM, and, therefore, this can be a very effective measure of stand density [55,77,102].

The prediction is done without calibration of the mixed-effects model; however, prediction accuracy can be very low [51,53,58,91,103]. Therefore, the calibrated mixed-effects model is always recommended for more accurate prediction of a response variable (HDR in our case). However, prediction accuracy of the calibrated mixed-effects HDR model may not always be high, as this depends on the application conditions of the stands, such as vertical structure, number of sample trees selected, and the range of HDR values used to calibrate the mixed-effects HDR model. For stands with homogenous canopies, the measured HDR of a single tree per sample plot may provide a precise result [104], but more trees are needed for the stands with heterogeneous canopies [53,54,58,105,106]. A localization with the smallest sample size, i.e., one tree can be most efficient for application of the mixed-effects model [96]. The mixed-effects HDR model calibrated using HDR measurements of one to four randomly selected trees of a particular species or group of species depending on the availability of their numbers per sample plot resulted in a promising prediction accuracy (Table 3, Figure 8). However, some large prediction error can be seen in this figure, which might be caused by extreme HDR observations in the stands that have heterogeneous vertical structures. The prediction accuracy for such sample plots might be increased with increasing number of trees per sample plot for calibration. However, using many sample trees in calibrating the mixed-effects model leads to a higher inventory cost with a little gain of accuracy. Using one to four trees per sample plot may compromise between sampling cost, model use efficiency, and prediction accuracy [59]. Sample plot-specific HDR curves generated with the calibrated mixed-effects HDR model for each species or group of species show complete coverages of the measured HDR values in the validation data (Figure 9), indicating that our model is precise enough for the tree species recorded in the research sample plot inventory. The reason HDR is not covered by the model curves for some sample plots is due to the presence of extremely small or large heights, which is obvious in the stands that have heterogeneous vertical structures.

As mentioned in the introduction section, HDR has various forest management implications. The HDR can be used for designing density management diagrams. The vulnerability condition of trees and stands to the external forces (wind, snow and icing) may be assessed using HDR. The HDR is significantly correlated to the crown ratio (CR) (Figure 10), and this relationship might be used for their removal if any unfavorable properties are detected in the relationship [1,20,28]. For example, trees with small CR and large HDR may be removed, as these trees may be less productive. Similarly, trees with large CR and large HDR need to be removed, as these trees may be more prone to snow-, icing-, and wind-related damages. This may be true that trees with a combination of both small CR and large HDR or large CR and large HDR might be susceptible to icing and windstorm. However, it would be difficult to predict this combination without a joint modelling of HDR and CR. Due to a lack of CR information in the training data, we were unable to include CR into the HDR model.

Based on the results presented in the article, forest managers may develop the guidelines for minimizing the devastating effects of ice- and wind-related damages of forest stands, especially conifer stands. With better understanding of the susceptibility to the potential damages, forest managers may design silvicultural treatments based on a certain range of HDR [43,107]. A stand stability strategy can be developed through manipulation of factors affecting HDR. Keeping the HDR of all trees in any forest stand under a certain threshold, for example, below 1.0 (young spruce stands), may reduce the overall vulnerability to the devastating ice- and wind-related damages, as trees with a smaller HDR than this threshold could have higher stability. The HDR information (HDR model) may allow forest managers to identify practical stand structural targets and design thinning strategies accordingly. The HDR model may be used to assess the quality and efficiency of thinning as HDR is significantly affected by thinning and other selection felling. The HDR is used as an effective predictor into various forest models to describe competition effects [24,25,26,108]. The HDR differs from other competition indices in that HDR is a direct measure of competition whereas others are indirect measures [20]. The HDR is significantly affected by several factors such as stand dynamics, natural disturbances and silvicultural tending, and, therefore, knowledge of HDR is essential for better understanding of forest ecological processes. The frequency of icing and wind disturbances in the European countries has recently increased, intensifying the need to advisee forest practitioners on how to manage their forests, especially conifer forests, to minimize such consequences. The HDR is one, but not only one, of the possible indicators that can be used to assess the potential risk of a forest stand to snow-, ice- and wind-related damages. In this context, in-depth analyses of HDR using large datasets that include the interactive effects of the several factors affecting HDR could be timely and appealing. The HDR model can be included as a sub-model into the individual tree or stand simulators, which are very useful tools for decision making in forestry. The forest risk model, which includes all kinds of potential risks caused by factors of climate, anthropogenic and natural disturbances, and tree and stand dynamics, can be applied to evaluate the overall risk on the stability of forest stands. Thus, predicting HDR with our model may be one of the intermediate steps for developing stand stability strategies and increasing the productivity of forest stands. Our HDR model not only predicts HDR, but factors that modify HDR are of paramount interest in the overall assessment of site productivity, preparing thinning strategies, and density management diagrams. The HDR model may be used as a benchmark model for comparing HDR derived from diameter and height increment models and height–diameter models.

Manipulation of the factors affecting HDR may be possible to reduce the vulnerability of stands to the devastating ice- and wind-related damages. With the help of results presented in this article (Table 2, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 9), one can modify the HDR of a stand. For example, as magnitudes of the parameter estimates of HDOM and dq seem to be much more influential on HDR, forest managers need to focus more on the modification of HDOM and dq. For a given stand density condition, maintaining very high and low HDOM on the younger stands may result in a large HDR, which may not be favorable for these stands to be more stable. On the stands with very high site quality (e.g., high HDOM), forest managers may reduce the stand density (by increasing dq and DDOM) so that the HDR of individual trees can decrease through thickening of their boles relative to height growth. When a strong wind approaches the stands with very high HDOM, it may destroy the top canopy trees while most of the suppressed trees usually manage to survive, and this is the season that the selection-type forest management is considered to be more ecologically stable. Understanding the links between risk assessment and forest management can be important in the context of climate change, whereby snow-, ice- and wind-related damages are common in the European forests. Only HDR prediction may be less important for forest managers, but what is really important for them is how they can manage their forests through modification of factors affecting HDR. Thus, in-depth exploration of the factors affecting HDR, which our study has done, would certainly be appealing. However, the information of many more factors rather than those analyzed in our study is also necessary for making the effective stand stability plans, as HDR is not the only indicator that perfectly gauges the vulnerability of a stand to icing and windstorm damages. The HDR model presented in this article is static one, which was evolved from a static relationship associating HDR with other variables. The dynamic HDR models (i.e., ΔHDR/Δt) can be more important than the static ones as HDR changes over time are included in the dynamic models. However, developing a dynamic HDR model needs longer-time series data. When these data are available, we will develop dynamic HDR models in the future.

5. Conclusions

This article presents a nonlinear mixed-effects HDR model applicable for several tree species that was developed using a large Czech NFI dataset through integration of covariate predictors describing the effects of inter-tree spacing and competition, site quality and stand development stage, canopy height differentiation and social position of individual trees, and an unstructured random component describing sample plot-level random variations. The model was tested using a large independent dataset acquired from research sample plots. The model described most of the HDR variations in the training data without significant residual trends. The model test also confirmed its high precision. The HDR can be predicted with reasonable accuracy by a properly calibrated mixed-effects model using measurements of DBH and other predictors that can be easily derived from the routine forest inventories. For higher prediction precision, the mixed-effects HDR model should be calibrated using HDR measurements of one to four randomly selected trees of a particular species or group of species depending on the availability of their numbers per sample plot.

The main uses of the HDR model can be for tree and stand stability assessments and regulation of stand density. The model presented in this article can be used for the formulation and implementation of silvicultural strategies for the species of interest. Silviculture strategies leading to increased stability have become increasingly important in terms of advancing global changes. Silviculture can play an important role in mitigating the climate change effects on the growth and stability of the forest stands, adapting the strategies that would increase the structural diversity of forests. Many other implications of HDR are also presented in the introduction and discussion sections of this article.