LiDAR-Based Modeling of Individual Tree Height to Crown Base in Picea crassifolia Kom. in Northern China: Comparing Bayesian, Gaussian Process, and Random Forest Approaches

Abstract

This study compared hierarchical Bayesian, mixed-effects Gaussian process regression, and random forest models for predicting height to crown base (HCB) in Qinghai spruce (Picea crassifolia Kom.) forests using LiDAR-derived data. Both modeling approaches were applied to a dataset of 510 trees from 16 plots in northern China. The models incorporated tree-level variables (height, diameter at breast height, crown projection area) and plot-level spatial competition indices. Model performance was evaluated using leave-one-plot-out cross-validation. The Gaussian mixed-effects process model (with an RMSE of 1.59 and MAE of 1.25) slightly outperformed the hierarchical Bayesian model and the random forest model. Both models identified LiDAR-derived tree height, DBH, and LiDAR-derived crown projection area as primary factors influencing HCB. The spatial competition index (SCI) emerged as the most effective random effect, with the lowest AIC and BIC values, highlighting the importance of local competition dynamics in HCB formation. Uncertainty analysis revealed consistent patterns across the predicted values, with an average relative uncertainty of 33.89% for the Gaussian process model. These findings provide valuable insights for forest management and suggest that incorporating spatial competition indices can enhance HCB predictions.

1. Introduction

Height to crown base (HCB) is an important tree attribute that determines the quality class of sawn logs and the value of timber [1]. At the stand level, mean HCB is also an indicator of wood quality [2,3]. Theoretically, HCB is a critical parameter in constructing allometric growth models of trees [4] and plays an important role in the allometric equations of biomass components, such as branches, bark, and foliage; it holds considerable value in elucidating the intrinsic connections between crown and trunk growth in various processes, such as carbon allocation and nutrient transfer [5,6]. Practically, HCB serves as a crucial basis for designing harvesting plans and determining appropriate harvesting intensity [7,8]. The accurate prediction of HCB can maximize the usable timber volume, reduce harvesting residues, and enhance both wood utilization and economic benefits [9,10].

Moreover, HCB is closely linked to stand tending and fire management. The proper regulation of HCB can improve light conditions within the stand and promote tree growth through tending measures, such as pruning [11,12]. On the other hand, HCB significantly influences the risk and intensity of wildfires’ spread [13,14]. Therefore, effective management of HCB is crucial for the prevention and control of forest fires [15,16].

Despite the significant importance of height to crown base (HCB) in forestry science and management, current research on HCB remains insufficient in several areas. Firstly, many existing HCB models are based on traditional statistical methods, such as linear regression, which often fail to capture the complexity and spatial heterogeneity of the processes involved in HCB formation. Advanced modeling approaches, including nonlinear simultaneous regression and the use of remote sensing data, such as LiDAR data, have shown promise in improving the accuracy of HCB predictions by accounting for these complexities [17]. Maltamo [18] developed HCB models by applying the linear mixed-effects (LME) modeling approach using LiDAR-retrieved percentiles of tree height, diameter at breast height (DBH), and field-measured tree height. However, given that HCB data often exhibit distinct nonlinear patterns, the use of linear functions to describe such relationships may be inadequate. Recognizing this limitation, Yang [17] established an HCB model based on LiDAR data using a nonlinear mixed-effects approach. Nevertheless, these studies have not fully addressed the inherent uncertainties associated with HCB modeling. Furthermore, the impact of spatial structure on HCB modeling has not been adequately accounted for in existing research. Spatial structure, including tree distribution patterns and competitive relationships, significantly influences trees’ growth and morphological development. Incorporating these spatial elements into the model could potentially enhance the accuracy of HCB predictions substantially.

Enhanced predictive accuracy of HCB models is crucial for both practical forestry operations and theoretical research [19]. Accurate HCB models can not only optimize forest management practices, such as thinning and pruning, to improve stand productivity and health, but they can also deepen our understanding of the processes governing trees’ growth and development, thereby providing important theoretical support for sustainable forest management.

Another crucial factor is data precision. This study employs data acquired through UAV-LiDAR technology, which provides more accurate measurements of tree height (H) and crown area compared to traditional ground-based techniques. UAV-LiDAR not only covers larger areas but also offers higher resolution and more precise three-dimensional forest structure information, thereby furnishing more reliable input data for HCB modeling. The use of high-precision UAV-LiDAR data is one of the major innovations of this study.

This study aims to systematically assess the uncertainties in HCB predictions using hierarchal Bayesian methods, mixed-effects Gaussian process regression models, and the random forest model, which has rarely been addressed in previous studies. A comprehensive exploration of uncertainty is crucial for understanding the reliability and limitations of models’ predictions, especially when dealing with complex ecological processes, such as tree growth and crown development [20]. Failing to account for various sources of uncertainty, such as measurement errors, model misspecification, and parameter variability, may lead to overconfident or misleading inferences [21,22]. Uncertainty analysis not only helps to understand the reliability and limitations of models’ predictions but also serves as an important evaluation metric for synthesis science, deepening our understanding of ecosystem processes and functions [23].

Bayesian inference is particularly well suited for this purpose, as it is based on probability theory and provides estimates of parameter uncertainty in the form of probability distributions, explicitly considering prior knowledge [24]. This contrasts with the frequentist approach, which treats model parameters as fixed and provides point estimates [25,26]. Bayesian methods have been increasingly applied in ecological and forestry research, providing new perspectives for modeling complex ecological data.

In addition to hierarchical Bayesian models, Gaussian process models have gained increasing attention in ecological and environmental modeling due to their flexibility and ability to capture complex nonlinear relationships [27]. Gaussian process models are nonparametric Bayesian models that define a prior distribution over functions, which can be updated based on observed data to obtain a posterior distribution [28]. The key advantage of Gaussian process models lies in their ability to provide a natural measure of uncertainty in predictions, taking into account both the inherent variability in the data and the uncertainty in the model itself [29]. Moreover, Gaussian process models can seamlessly incorporate various types of input data, such as remote sensing data and spatial information, making them particularly suitable for modeling ecological processes at different scales [30].

In addition to Bayesian and Gaussian process models, random forest models offer a powerful alternative for ecological modeling. Random forest models are an ensemble learning technique that builds multiple decision trees using different subsets of the training data and combines their predictions to improve robustness and reduce overfitting. This method does not rely on strict parametric assumptions, making it a versatile choice for a wide range of ecological applications, including the analysis of tree growth, species distribution, and environmental influences.

Hierarchical models, such as hierarchical Bayesian models and Gaussian process mixed-effects models, offer several advantages in addressing the challenges associated with modeling complex ecological data [31,32]. These advantages include (1) improving estimates for repeated measurements by avoiding underfitting or overfitting, (2) handling imbalanced sampling by automatically accounting for varying uncertainty across groups, (3) explicitly estimating variation at different levels, and (4) retaining variation by avoiding the pitfalls of pre-averaging data [33,34].

The main objective of this research is to develop a suitable LiDAR-based individual HCB prediction model for Picea crassifolia Kom. with canopy cover around 0.6, specifically in the Dayekou Basin. Therefore, we strive to construct a relatively complex model with strong predictive capability while accounting for plot-level variability. Compared to existing HCB models, the innovations of the proposed model in this study include the following: (1) focusing on the natural secondary forests of Picea crassifolia, a regionally important forest type, to fill the gap in the current growth model system; (2) introducing novel spatial structure indices to reveal the intrinsic relationship between spatial patterns and tree growth; and (3) comparing the performance of Bayesian models, Gaussian process regression models, and random forest models to provide a reference for model selection in future HCB research.

2. Methods

2.1. Data

2.1.1. Study Area

The study area was the Dayekou Observation Station of the Xishui forest farm in Sunan Yuguzu Autonomous County, Gansu Qilian Mountain National Nature Reserve (38°32′ N, 100°15′ E) (Figure 1). This area is a natural secondary forest in the shallow mountain area of the Qilian Mountains, belonging to the water conservation forest of Heihe River Basin, with an altitude of 2700–3000 m. The climate of this area is an alpine, semi-arid, and semi-humid mountain forest steppe climate. The slope exposed to the Sun is mountain grassland, while the northern slope is the forest landscape. The dominant species in the forest is Picea crassifolia Kom. The ground is covered with a moss layer ~10 cm thick, and the vegetation is in good condition.

2.1.2. Plot Setting

A 100 m × 100 m grid was laid out along the slope with an azimuth angle of 122 degrees. The grid was divided into sixteen 25 m × 25 m plots, and the stakes were set up at corners and at the approximate center points of each sample plot (shown as Figure 2). During the plot layout, a compass and measuring tape were used, with the compass for determining direction and the tape for measuring distances. During individual tree measurements, all trees with a DBH greater than 2.5 cm within the plot were surveyed, including leaning trees, standing dead trees, and snags. DBH was measured using a diameter tape, while HCB was measured using a laser rangefinder (TruPulse 200, Laser Technology, Inc., Centennial, CO, USA). For shorter trees, a measuring pole was used to determine the tree height and the height to the base of the crown. Crown width was measured using a tape measure along two main directions. Additionally, leaf area index (LAI) measurements were conducted using the LAI2000 (LI-COR Biosciences, Lincoln, NE, USA) and a fisheye camera.

2.1.3. ALS Data

The full-waveform LiDAR data were collected using a Riegel LMS-Q560 laser scanner (RIEGL Laser Measurement Systems GmbH, Horn, Austria) mounted on an airplane, with a 50 kHz pulse repetition frequency, a 49 kHz scanning frequency, and a 30° maximum scanning angle. The ground-measured data were also collected. The wavelength was 1550 nm, and the pulse length and laser beam divergence were 3.5 ns and 0.5 mrad, respectively. The average flight height was 3699 m, and the average flight speed was 230 $\mathrm{k}\mathrm{m}/{\mathrm{h}}^{-1}$. The mean density of point clouds was 4.34 m⁻². The spatial distribution of 510 neighbor-smoothed Picea crassifolia Kom. trees is shown in Figure 2. The data are summarized in Table 1.

The point cloud was generated from LiDAR waveforms provided by Pang et al. [35]. Ground points were classified using the Ground algorithm in TerraScan 005 (Terrasolid, Helsinki, Finland), following Axelsson’s method (1999) [36]. These points were used to create a digital elevation model (DEM) at a 0.25 m resolution using a triangulation algorithm. Similarly, a digital surface model (DSM) at a 0.25 m resolution was generated using ground and vegetation points with the highest hit z algorithm, assigning the highest z-value within each cell. The canopy height model (CHM) at a 0.25 m resolution was obtained by subtracting the DEM from the DSM [37]. To smooth pits in the CHM, a neighbor smoothing algorithm was applied [38].

Crown tops were identified from the CHM using a local maxima algorithm with a 2.0 m window size, and LiDAR-derived tree height (hL) was estimated as the height of these crown tops. Crown edges were delineated using a region-growing algorithm, and the LiDAR-derived crown area (caL) was calculated as the average horizontal diameter in the west–east and north–south directions [37].

2.2. Models

2.2.1. HCB Base Model

We applied the base model in logistic form, similar to Water and Hanns’s research [39], and the logistic form base model was introduced with variables, such as individual hL (LiDAR-based tree height), caL (LiDAR-based crown projection area), DBH (diameter at breast height), as follows:

(1)$HCB=hL/(1.0+ exp\mathit{(}{\beta }_0+{\beta }_{1}DBH+{\beta }_{2}caL\mathit{\left)}\right)$

2.2.2. Included Spatial Covariates

We calculated four spatial indices to characterize the competitive environment and spatial structure of the forest stand, namely Hegyi’s competition index (CI), the aggregation index (R), and the angular competition index (aCI), as follows:

(2)$\mathit{hcb} = \mathit{hL}/\mathit{(}1.0+\mathit{exp}({\beta }_0+{\beta }_1\mathit{DBH}+{\beta }_2\mathit{caL}+{\beta }_{3}S\left)\right)$

where S represents the spatial covariate included in the model.

Hegyi’s CI [40] quantifies the competitive pressure on a target tree from its neighbors. For each target tree i, ${CI}_i$ was calculated as follows:

(3)${CI}_i=\sum {\displaystyle \frac{D_j}{D_i\ast {dist}_{ij}}}$

where $D_i$ and $D_j$ are the diameter at breast height of the target tree i and competitor tree j, respectively, and ${dist}_{ij}$ is the distance between trees i and j. Only trees within a 7 m radius were considered competitors.

The aggregation index (R) [41] measures the spatial distribution pattern of trees within a plot as follows:

(4)$R=\left({\displaystyle \frac{1}{n}}\right)\ast \sum r_i$

where n is the number of trees in the plot and $r_i$ is the distance from tree i to its nearest neighbor.

The aCI [42] considers both the size and spatial arrangement of neighboring trees as follows:

(5)$aCI=\sum {\displaystyle \frac{(a_1+a_2\ast c_{ij})}{\pi }}$

where $a_1=arctan\left({\displaystyle \frac{H_j}{d_{ij}}}\right)$, $a_2=arctan\left({\displaystyle \frac{H_j-H_i}{d_{ij}}}\right)$, and $C_{ij}$ is 1 if $H_j$ > $H_i$ and 0 otherwise. $H_i$ and $H_j$ are the heights of the target and competitor trees, respectively, and $d_{ij}$ is the distance between them. The four nearest neighbors were considered for each target tree.

All indices were calculated using R, with tree locations specified by X and Y coordinates. The results provide a comprehensive characterization of the competitive environment and spatial structure of the forest stand.

2.2.3. HCB Hierarchical Bayesian Model

• (1). Bayes’ rule

Bayes’ rule is usually expressed as follows:

(6)$p(y,\theta )=p\left(y\right|\theta \left)p\right(\theta )=p(\theta \left|y\right)p\left(y\right)$

where $y=(y_1,y_2,y_3,...)$ is a vector of data and $\theta =({\theta }_1,{\theta }_2,{\theta }_3,...)$ is a vector of the parameters. p represents the density function or probability distribution, $p(y,\theta )$ is the joint probability of y and $\theta$ and, $p\left(\theta \right)$ is the marginal density of $\theta$.

Unlike a statistical reference, a Bayesian reference uses conditional probability distribution to describe uncertainty in the estimated parameters. In order to estimate parameters, posterior distribution should be obtained as follows:

(7)$p(\theta |y)={\displaystyle \frac{P\left(y\right|\theta \left)p\right(\theta )}{p\left(y\right)}}$

where $p\left(\theta \right|y)$ is the posterior distribution of parameters, $p\left(\theta \right)$ is the prior distribution of parameters, $P\left(y\right|\theta )$ is the likelihood function, and p(y) is the average likelihood, which is used for normalization to guarantee that the integral value of posterior distribution of the defined domain is 1.

Therefore, the posterior distribution is proportional to the product of likelihood and prior distribution, as follows:

(8)$p\left(\theta \right|y)\propto P(y\left|\theta \right)p\left(\theta \right)$

In our study, the relationship between HCB, DBH, and hL was described as follows:

(9)$\mathit{HCB}=\mathit{hL}/(1.0+\exp ({\beta }_0+{\beta }_{1}DBH+{\beta }_{2}caL\left)\right)$

(10)$HCB\sim N\left(f\right(hL,DBH,caL:{\beta }_0,{\beta }_1,{\beta }_2)$

where $f(hL,DBH,caL:{\beta }_0,{\beta }_1,{\beta }_2)=hL/(1.0+exp({\beta }_0+{\beta }_{1}DBH+{\beta }_{2}caL\left)\right)$ gives the mean distribution of HCB. Therefore, the likelihood function can be expressed as follows:

(11)$p(y|{\beta }_0,{\beta }_1,{\beta }_2)={\mathit{\Pi }}_i{\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}exp({\displaystyle \frac{-({HCB}_i-f(h{L}_i,DB{H}_i,ca{L}_i){)}^2}{2{\sigma }^2}})$

where i is the number of observations.

• (2). Hierarchical Bayesian model development

We developed a hierarchical Bayesian model with random effects to predict the height to crown base (HCB) of trees. The model was formulated as follows:

(12)$HCB=LH/(1+exp({\displaystyle \sum }\left({\beta }_i{X}_i\right)+{\displaystyle \sum }\left(u_j{X}_j\right)\left)\right)$

where $hL$ is the tree height, $X_i$ represents fixed-effect predictors (DBH, caL, and spatial competition index), $X_j$ represents predictors with random effects, ${\beta }_i$ represents fixed-effect coefficients, and $u_j$ represents random-effects coefficients varying by plot.

We assumed normal priors for fixed effects (${\beta }_i~N(0, {10}^2)$) and random effects ${(u}_j~N(0, {\sigma }_u^2)),$ with a half-normal prior on the random-effects standard deviation $\left({\sigma }_u~HalfNormal\left(1\right)\right).$ The model error was assumed to be normally distributed, with a half-normal prior on its standard deviation $(\sigma ~HalfNormal\left(1\right))$.

• (3). Parameter estimation

We implemented the model using PyMC3 running on Python version 3.9 and estimated the parameters via Markov chain Monte Carlo (MCMC) sampling, with 2000 tuning steps and 2000 sampling steps. We explored all possible combinations of random effects on predictors (excluding tree height) and compared the models using the root-mean-square error (RMSE), Akaike information criterion (AIC), Bayesian information criterion (BIC), and log-likelihood.

• (4). Uncertainty analysis of Hierarchical Bayesian Model

In Bayesian modeling using MCMC sampling, an uncertainty analysis was conducted for comprehensive decision-making based on full datasets. The uncertainty in HCB estimation was analyzed based on three main aspects, as follows:

a. Posterior predictive intervals

A total of 1000 parameter samples were drawn from the posterior distribution, and 90% posterior predictive intervals were calculated for all observations. The proportion of observed values falling within their respective intervals was assessed to confirm the reliability of the model and the uncertainty estimates.

b. Marginal posterior distributions

The marginal posterior distributions for each model parameter were summarized and visualized. The posterior distributions of regression coefficients ${\beta }_1$ to ${\beta }_3$ were examined for their significance in relation to zero, in order to determine the effects of crown diameter, live crown length, and total height on HCB. The posterior distributions of the random-effect variance parameters ${\sigma }_1$ and ${\sigma }_2$ were analyzed to assess how well both inter-plot and intra-plot HCB variations were characterized. The pseudocode of hierarchical Bayesian model was provided in Appendix A.

2.3. Mixed-Effects Gaussian Process Regression Models

2.3.1. Gaussian Process Regression (GPR)

GPR is a nonparametric Bayesian machine learning method used to model the nonlinear relationships between input features and target variables [27]. Given a training dataset $D=\left\{\right(x_i,y_i\left)\right|i=\mathrm{1,2},\cdots , n\}$, where $x_i\in R^d$ represents the input feature vector (DBH, hL, caL, and one of the four spatial variables) and $y_i\in R$ represents the corresponding target value (HCB), the goal of GPR is to learn a mapping function $f: R^d\to R$ from the input space to the target space.

In GPR, the mapping function f is assumed to be a Gaussian process with a prior distribution defined as follows:

(13)$f\left(x\right)~GO(m\left(x\right), k\left(x,x^{\prime }\right))$

where $m\left(x\right)$ is the mean function and $k\left(x,x^{\prime }\right)$ is the covariance function, also known as the kernel function. Commonly used kernel functions include the squared exponential kernel, the Matérn kernel, and the periodic kernel [43]. Given the training dataset D, the posterior distribution of GPR can be obtained using Bayes’ theorem, as follows:

(14)$p\left(f|D\right)={\displaystyle \frac{p\left(D|f\right)p\left(f\right)}{p\left(D\right)}}$

where $p\left(D|f\right)$ is the likelihood function, $p\left(f\right)$ is the prior distribution, and $p\left(D\right)$ is the marginal likelihood. By maximizing the log-marginal likelihood, the hyperparameters of the Gaussian process, such as the parameters of the kernel function and the noise variance, can be estimated.

For a new input $x_{\ast }$, GPR can predict the distribution of the corresponding target value as follows:

(15)$p\left(y_{\ast }|x_{\ast }, D\right)=\int p\left(y_{\ast }|x_{\ast }, f\right)p\left(f|D\right)df$

The mean and variance of this predictive distribution can be analytically derived, representing the prediction and the uncertainty estimation, respectively.

In this study, GPR was employed to develop a model for predicting HCB. The input features included caL, hL, DBH, and one of the four spatial pattern indices, while the target variable was HCB. We selected the squared exponential kernel as the covariance function and optimized the kernel hyperparameters using maximum likelihood estimation.

2.3.2. Mixed-Effects Gaussian Process Regression

To account for plot-level random effects, we extended the GPR model to a mixed-effects Gaussian process regression (ME-GPR) [27]. Consider m plots, each containing n_j tree samples, where $j=\mathrm{1,2},\cdots ,m$. Let $x_{ij}\in R^d$ denote the input features of the $y_{ij}\in R$ representing the corresponding height to crown base. The ME-GPR model can be expressed as follows:

(16)$y_{ij}=f\left(x_{ij}\right)+b_j+{\epsilon }_{ij}$

where $f\left(x_{ij}\right)$ is the fixed-effects component from the GPR model, b_j is the random effect for the j-th plot following a Gaussian distribution $b_j~N(0, {\sigma }_b^2)$, and ${ϵ}_{ij}$ is independent and identically distributed Gaussian noise, ${\epsilon }_{ij}~N(0, {\sigma }_n^2)$.

In ME-GPR, we assume that the fixed-effects component f(x_i) has a Gaussian process prior, as follows:

(17)$f\left(x\right)~GP(0, k\left(x,x^{\prime }\right))$

where $k\left(x,x^{\prime }\right)$ is the covariance function. The prior distribution of the random effects (b_j) is as follows:

(18)$b~N(0,{\sigma }_b^{2}I)$

where $b={[b_1,b_2,\cdots ,b_m]}^T$ is the vector of all plot random effects, and I is the identity matrix.

Given the training dataset $D=\{(x_{ij},y_{ij})|i=\mathrm{1,2},\cdots ,n_j;j=\mathrm{1,2},\cdots ,m\}$, the posterior distribution of ME-GPR can be obtained through Bayes’ theorem, as follows:

(19)$p\left(f,b|D\right)={\displaystyle \frac{p\left(D|f,b\right)p\left(f\right)p\left(b\right)}{p\left(D\right)}}$

where p(D|f, b) is the likelihood function, p(f) and p(b) are the priors of fixed and random effects, respectively, and p(D) is the marginal likelihood. Hyperparameters of the Gaussian process and the variance of random effects ${\sigma }_b^2$ can be simultaneously estimated by maximizing the log-marginal likelihood.

For a new input x*, ME-GPR predicts the distribution of its corresponding target value as follows:

(20)$p\left(y_{\ast }|x_{\ast }, D\right)=\uparrow \iint p\left(y_{\ast }|x_{\ast }, f, b_{\ast }\right)p\left(f, b_{\ast }|D\right)df d{b}_{\ast }$

where $b_{\ast }$ is the plot random effect associated with the new input $x_{\ast }$. The mean and variance of the predictive distribution can be obtained through numerical integration or approximate inference methods [44].

In this study, we employed ME-GPR to construct a prediction model for height to crown base, incorporating tree-level input features (hL, DBH, caL, and spatial pattern indices) and plot-level random effects. We selected the squared exponential kernel as the covariance function and used maximum likelihood estimation to optimize the model’s hyperparameters. The pseudocode of mixed-effects Gaussian process regression model was provided in Appendix A.

2.3.3. Uncertainty Analysis of Mixed-Effects Gaussian Process Regression Models

To complement the ME-GPR model, a comprehensive uncertainty analysis was conducted based on the full dataset. This analysis aimed to quantify and characterize the various sources of uncertainty in the model predictions, focusing on the following three main aspects:

a. Posterior predictive distributions

For each new input x*, the posterior predictive distribution p(y*|x*, D) was computed using Equation (21). From these distributions, 95% credible intervals were calculated for each prediction. These intervals represent the ranges within which the true HCB value is expected to lie with 95% probability, given the model and data. The proportion of observed values falling within their respective credible intervals was calculated in order to assess the reliability of the model’s uncertainty estimates.

b. Marginal posterior distributions

The posterior distributions of the model parameters were analyzed to assess parameter uncertainty. This included examining the posterior distributions of the fixed-effects coefficients (associated with caL, hL, DBH, and spatial pattern indices) to determine their significance and variability, analyzing the posterior distributions of the random-effects variance components to characterize the inter-plot variations in HCB, investigating the posterior distributions of the kernel function hyperparameters to understand the learned structure of the Gaussian process, and computing credible intervals and summary statistics for each parameter to quantify their uncertainty.

2.4. Random Forest Model

2.4.1. Random Forest HCB Model

The random forest technique is an ensemble learning technique that constructs multiple decision trees, each trained on a random subset of the data and features, to perform regression or classification tasks. The final prediction is obtained by averaging the results from individual trees, thereby reducing the risk of overfitting—a common issue with single decision trees [45]. This approach is particularly suitable for predicting height to crown base (HCB), a key parameter in forestry that influences light penetration, understory dynamics, and timber quality.

In this study, the random forest model was applied to predict HCB using key predictor variables, including LiDAR-derived tree height (hL), diameter at breast height (DBH), crown projection area (caL), and the spatial competition index (SCI). Model parameters, such as the number of trees (n_estimators = 500) and the number of features considered at each split (mtry), were optimized through 10-fold cross-validation to achieve a balance between accuracy and computational efficiency.

2.4.2. Bootstrap Method for Uncertainty Estimation

To quantify the uncertainty of the random forest model’s predictions, a bootstrap method was applied to each cross-validation training set. For each training set in a cross-validation fold, 100 bootstrap samples were generated by resampling the training data with replacement. For each bootstrap sample, a random forest model was trained, and the trained model was used to predict the response variable for the corresponding test data (the held-out fold).

For each test data point, the mean and standard deviation (SD) of predictions across all 100 bootstrap iterations were calculated. The mean prediction represents the central tendency of the RF predictions for each test sample. The standard error (SE) was calculated as the standard deviation of the predictions, which provides an estimate of the uncertainty associated with each prediction.

2.5. Model Validation

For model validation, new data could be applied; however, as our data were limited in number, it was not practical to perform such a validation. In our study, we applied a cross-validation approach to validate the Bayes models and Gaussian models, similar to the research of Fu [46]. Each time, we removed one plot and used the remaining plots to fit the HCB model and predict the removed one (leave-one-plot-out cross-validation (LOPCV)). This process was repeated N − 1 times, where N is the total number of plots. This method was chosen because it allowed us to assess the model’s ability to generalize to new, unseen plots, which is critical given the inclusion of random effects related to plot variability. For the machine learning algorithm (random forest model), we used a 10-fold cross-validation method. This approach is particularly suited for machine learning models that do not explicitly account for hierarchical structure, allowing for a robust assessment of the model’s overall predictive accuracy. The predictive performance of the model was assessed using RMSE and MAE. The model with the smallest RMSE and MAE was chosen as the final model to predict HCB. We performed all computations using R software, version 4.3.2 [47]. Equations (21) and (22) are as follows:

(21)$RMSE=\sqrt{{\overline{e}}^2+{\sigma }_e^2}$

(22)$\mathit{MAE} =|{\displaystyle \sum }_{i=1}^N{e}_i|/N=|{\displaystyle \sum }_{i=1}^N({HCB}_i-{\widehat{HCB}}_i)|/N$

where N is the number of observations, RMSE is the root-mean-square error, and MAE is the mean absolute error.

3. Results

3.1. Bayes Model Development

3.1.1. Hierarchical Bayes Base Model

The base model yielded parameter estimates (Table 2) for the two main coefficients (b1 and b2) and the error term (sigma). The first coefficient, b1, had a mean of 0.0113 (SD = 0.0032), with a 95% confidence interval of [0.0050, 0.0175]. The second coefficient, b2, showed a mean of 0.0262 (SD = 0.0058), with a 95% confidence interval of [0.0149, 0.0377]. Both coefficients were statistically significant, as their confidence intervals did not include zero. These results suggest that both predictors have a significant effect on the outcome variable, with the model accounting for some of the variability in the data.

The Bayesian base model yielded the following equation:

(23)$HCB={\displaystyle \frac{X_0}{1+\mathit{exp}\left(0.0113X_1+0.0262X_2\right)}}+\epsilon$

where:

HCB is the dependent variable;

X₀, X₁, and X₂ are the independent variables;

$\epsilon$ is the error term, normally distributed with a mean of 0 and a standard deviation of 1.7840.

Both coefficients were statistically significant (β₁: 95% CI [0.0050, 0.0175]; β₂: 95% CI [0.0149, 0.0377]). The model demonstrated modest predictive power (MAE = 3.17, RMSE = 1.78), suggesting a nonlinear relationship between the predictors and HCB. This structure may reflect complex underlying processes, although the moderate fit indicates that additional factors could improve the model’s performance.

3.1.2. Spatial Variables Introduced

The results show a performance ranking of five indices based on LOPCV (Table 3):

SCI performed best, with the lowest LOO value (−1135.88) and p loo (effective number of parameters) (4.61), along with a very low Dse (difference in standard error), because lower LOO values represent better performance, low p loo values represent a simple model, and lower Dse indicates proximity to the best model. The aCI was a close second, with slightly higher LOO (−1138.26) and p loo (4.70). The Hegyi indices performed significantly worse, with much higher LOO values (both above −1356). SCI and aCI demonstrated superior predictive accuracy and model parsimony. The aggregation index received a warning flag, suggesting potential issues. All of these models use a logarithmic LOO scale.

In summary, SCI and aCI outperformed the other indices across all metrics, indicating that they are likely the most suitable for accurate and reliable predictions in this context.

3.1.3. Hierarchical Bayes Model

• (1). Random effect selection

To investigate the influence of plot-level random effects on the prediction of HCB, we fitted a series of mixed-effects logistic regression models with different combinations of random effects. The model performance metrics, including the Akaike information criterion (AIC), Bayesian information criterion (BIC), and log-likelihood, were calculated for each model (Table 4).

Although the results showed that the model with random effects on all three predictor variables had the lowest AIC and BIC values, the difference in AIC between random effects on all three variables (1676.60) and DBH and SCI (1796.77) was not very big. Considering the convergence cost, the model with random effects on DBH and SCI was screened out as the most effective HCB hierarchical Bayes model, indicating the best balance between model fit and complexity.

• (2). Hierarchical Bayes model

We implemented a hierarchical Bayes model on the dataset using LOPCV. The model achieved an average MSE of 2.81 and an RMSE of 1.68, indicating moderate predictive accuracy.

The hierarchical Bayesian model yielded the following equation:

(24)$HCB={\displaystyle \frac{hL}{1+\mathit{exp}\left((0.0166+u_1)DBH+0.0272caL-(0.0902+u_2)SCI\right)}}+\epsilon$

where:

u₁ and u₂ are the random-effect variables;

$\epsilon$ is the error term, normally distributed, with a mean of 0 and a standard deviation of 1.7840.

• (3). Uncertainty analysis of hierarchical Bayes model

All parameters’ posterior distributions exhibit distinct bell-shaped curves, indicating some uncertainty in the model’s parameter estimation but relatively stable estimates. The estimated values of b₁ and b₂ are positive, with 95% HDI not including zero, suggesting potential statistical significance for these parameters. The estimated value of b₃ is negative, but its 95% HDI includes zero, indicating that b_3′s effect may not be statistically significant. The large estimated value of sigma indicates potentially high residual variance in the model. Most of the parameters have relatively narrow HDI intervals, suggesting high precision in estimation. However, b_3′s wider HDI interval indicates greater uncertainty in its estimation (Figure 3).

The distribution of the model’s residuals shows an approximately normal shape, centered near zero (Figure 4). Most residuals fall within the range of −2 to 2, indicating a reasonable model fit. There is a slight positive skew, suggesting that the model may underestimate some higher values. The peak of the distribution occurs around zero, with the majority of residuals clustered within ±1 units. While outliers exist beyond ±3 units, they are infrequent, pointing to generally good predictive accuracy, with some exceptions.

3.2. Gaussian Process Model

3.2.1. Gaussian Process Base Model

For the Gaussian process model, the kernel is a combination of a radial basis function (RBF) and a white kernel, as follows:

(25)$k\left(x, x^{\prime }\right)={0.716}^{2}RBF\left(length scale=1.85\right)+WhiteKernel(noise level=0.649)$

The Gaussian process base model demonstrated moderate predictive performance. The model’s accuracy is reflected in its MAE of 2.70 and RMSE of 1.64. The model’s complexity and fit are captured by an AIC of 415.15 and BIC of 449.27, with a log-likelihood of −194.57.

The kernel function, which defines the covariance structure, combines an RBF term with a length scale of 1.85 and a white noise term. This structure suggests that the model both captures smooth trends in the data and accounts for noise. The mean prediction interval of 6.66 and mean confidence interval of 3.33 provide insights into the model’s uncertainty estimates. The mean total variance of 2.90 indicates the overall variability in the model’s predictions.

3.2.2. Random Effect Selection

The comparison of GPR with different random effects revealed that the model incorporating random effects on SCI outperformed the other models (Table 5). This model demonstrated the lowest Akaike information criterion (AIC) of 401.03, indicating the best balance between model fit and complexity. The kernel for this model, a rational quadratic function with alpha = 0.0332 and length scale = 0.0736, captures both large-scale and small-scale variations in the data. These results indicate that accounting for random effects in SCI provides the most effective model for describing the underlying patterns in the dataset, potentially reflecting important spatial competition dynamics in the system under study.

3.2.3. Gaussian Mixed-Effects Process Model

We implemented a GPR with random effects on the dataset using LOPCV. Figure 5 illustrates the true versus predicted values with uncertainty, demonstrating the model’s performance across the range of observed values. The model achieved an average MAE of 1.25 and an RMSE of 1.59, indicating moderate predictive accuracy.

3.2.4. Uncertainty Analysis of Gaussian Mixed-Effects Process Model

The model’s uncertainty estimation was robust, with a 95% prediction interval coverage of 0.95, closely aligning with the expected 95% coverage. Figure 6 depicts the predictions with uncertainty (a) and the relative uncertainty versus true values (b), revealing a consistent pattern of uncertainty across the range of true values. The average relative uncertainty was 0.3389, indicating that the model’s predictions had an average uncertainty of ~33.89% relative to the predicted values. The average error variance was 2.61, reflecting the model’s estimate of the inherent variability in the predictions.

3.3. Random Forest Model

3.3.1. Random Forest HCB Model Hyperparameter Turning and Model Validation

In the process of tuning the random forest model, various values of the hyperparameters mtry and min.node.size were evaluated using a 10-fold cross-validation approach. Given the relatively small sample size, the number of trees (ntree) was set to 500. The results (Table 6) indicated that the optimal parameter combination was mtry = 3 and min.node.size = 5, which yielded the lowest RMSE of 1.639, suggesting the best predictive performance among the tested configurations. Additionally, the mean absolute error (MAE) for the optimal configuration was 1.300, indicating a relatively low average prediction error. Therefore, the best-performing HCB prediction model based on the random forest approach was the one configured with mtry = 3, min.node.size = 5, and ntree = 500.

3.3.2. Uncertainty Analysis of Random Forest Model

The random forest model’s predictions were analyzed for uncertainty using a bootstrap-based approach combined with 10-fold cross-validation. The resulting plot (Figure 7) displays the predicted values against the true values for the test data, including 95% confidence intervals for each prediction. The results showed the following: (1) predictions for lower true values (below 3.0) generally show wider intervals, indicating higher uncertainty in this region; (2) for intermediate true values (between 3.0 and 6.0), the intervals are narrower, suggesting more confident predictions in this range; (3) predictions for higher true values (above 6.0) tend to have slightly increased uncertainty, as shown by the wider error bars, which may reflect the model’s reduced ability to generalize in this range due to fewer training samples or increased data variability.

4. Discussion

This study employed hierarchical Bayesian, Gaussian process regression, and random forest models to predict HCB in Picea crassifolia forests. All models fitted the data well, explaining the majority of variation in HCB. The hierarchical Bayesian model achieved an RMSE of 1.6773 m, the Gaussian mixed-effects process model achieved an RMSE of 1.5902 m, and the random forest model achieved an RMSE of 1.6394 m. These results suggest that the Gaussian mixed-effects process model performs best, offering a robust alternative for HCB prediction despite the differences in methodology.

Regarding the influence of variables, the hierarchical Bayesian, Gaussian process, and random forest models consistently identified LiDAR-derived tree height, DBH, and LiDAR-derived crown projection area as the primary factors affecting HCB. Tree height exhibited a significant positive correlation with HCB, indicating that taller trees tend to have higher HCB. This relationship is consistent with the findings of Sharma and Parton [48] in their study on boreal forest species, suggesting a consistent pattern across different forest ecosystems. DBH also showed a positive correlation with HCB, which may be attributable to larger-diameter trees generally being older and having undergone a longer natural pruning process. Rijal et al. [49] observed a similar relationship in their study on Scots pine, further corroborating our findings. Interestingly, the live crown projection area was negatively correlated with HCB, implying that trees with larger crowns have relatively lower HCB. This inverse relationship might be explained by the allocation balance of carbohydrates between stem and crown organs during tree growth, as discussed by Mäkelä and Valentine [50]. The consistency of these relationships across models strengthens the reliability of our findings.

From an ecological perspective, factors influencing HCB can be categorized into two types: (1) tree growth and development characteristics, such as genetic traits, age, and growth rate; and (2) environmental conditions, including stand density and site index. These factors interact, creating spatial heterogeneity in HCB. Our study used plot-level random effects to account for some of this spatial variation. Notably, the superiority of the SCI in the Gaussian process model suggests that local competition dynamics play a crucial role in shaping HCB patterns. This highlights the importance of considering spatial relationships in the structural modeling of forests. Inter-plot differences in HCB may be related to the microenvironment, stand structure, soil fertility, and disturbance history. These factors, while not explicitly modeled in our study, likely contributed to the unexplained variance and represent an important avenue for future research. Understanding these spatial dynamics is crucial for developing more accurate predictive models and informing forest management practices aimed at optimizing stand structure and timber quality.

The findings of this study are generally consistent with those of previous research on other tree species and regions. For instance, Hasenauer and Monserud [51] found that HCB in Austrian spruce plantations was primarily governed by tree height and DBH. Similarly, Rijal [49] identified tree height as the most important variable for predicting HCB in Finnish pine forests. This suggests that the factors influencing HCB identified in this study may have some universality.

Both hierarchical Bayesian and Gaussian process models demonstrated strong performance in HCB prediction, each with distinct strengths. The hierarchical Bayesian model excels in its flexibility to incorporate prior knowledge, comprehensive quantification of uncertainty sources, and applicability to complex data structures, including spatial data and hyperparameters [32]. This makes it particularly suitable for capturing the multi-level nature of forest ecosystems. On the other hand, the Gaussian process model showed superior ability in fitting nonlinear relationships, providing probabilistic prediction intervals, and exhibiting strong adaptability [27]; its performance in our study, particularly with the SCI as a random effect, highlights its capacity to capture complex spatial competition dynamics.

The random forest model, while not outperforming the GPR model, still provided a competitive alternative for HCB prediction. The random forest approach is particularly advantageous due to its ability to handle large datasets, identify variable importance, and generate robust predictions without extensive assumptions about the underlying data distribution.

The HCB models developed in this study provide important references for the precision management of Picea crassifolia forests. Firstly, HCB is a crucial indicator for assessing usable timber volume, directly affecting the wood’s grade and value. Accurate HCB prediction aids in formulating rational harvesting plans to maximize volume utilization. Secondly, HCB significantly influences stand canopy closure, understory light conditions, and canopy structure, which are critical for stand regeneration and mixed-species regulation. Optimizing HCB can create favorable light and spatial conditions for seedling regeneration. Furthermore, changes in HCB are closely related to the spread of forest fires and pest outbreak processes. Accurate HCB estimation can inform forest disaster early warning and control measures [16].

Future research could consider the following directions: (1) expanding the scale of plot surveys to include more tree species and forest types, enhancing the models’ universality; (2) integrating remote sensing techniques to obtain high-spatiotemporal-resolution image data for monitoring and predicting dynamic HCB changes [14]; (3) developing model-based HCB management decision systems to directly serve production practices; and (4) incorporating more physiological and ecological processes to enhance models’ mechanistic understanding of HCB formation processes. Through multidisciplinary integration, the breadth and depth of HCB research can be continuously expanded.

The uncertainty analysis provides a more comprehensive understanding of the reliability of HCB predictions. Based on the current training data, the hierarchical Bayesian model and the Gaussian process regression model can provide reasonable and credible interval estimates and uncertainty measures, which serve as important references for the application of models’ results in practice. However, when using these models for estimation, it is crucial to fully consider the interval nature of the predictions and not to overly rely on point estimates [52]. The posterior distribution of the parameters in the hierarchical Bayesian model suggests that tree size has a stable promoting effect on HCB, while the influence of crown structure is weaker and potentially negative. The variation in HCB within and between plots is generally controllable. These findings are consistent with the results reported by Smith [53], who found that tree diameter and height were the primary drivers of HCB variation in a mixed coniferous forest.

While this study provides valuable insights into the factors influencing HCB in Qinghai spruce forests, its focus on natural spruce stands within a specific region may limit the generalizability of the findings to other forest types or geographic areas. Variations in ecological conditions, species composition, and forest management practices across different regions could influence the relationships between tree height, DBH, crown area, and HCB. For instance, the environmental conditions in boreal forests, temperate hardwood forests, or mixed-species stands might result in different dynamics of tree growth and competition [54]. Therefore, caution should be exercised when applying these findings beyond the specific context of Qinghai spruce forests. Expanding plot surveys to encompass a broader range of forest types could enhance the models’ general applicability [55]. Additionally, investigating how the identified relationships manifest under varying climatic and ecological conditions could provide a more comprehensive understanding of HCB dynamics.

In summary, due to the complexity of the HCB formation process, the existing models still struggle to exhaustively consider all influencing factors. Genetic, environmental, and disturbance variables not included in these models may be important sources of prediction errors [56]. This limitation presents an opportunity for future research to explore more comprehensive influencing factors through mechanistic studies, optimize the model structure and parameter space, and further improve the accuracy and robustness of the predictions [57]. Additionally, future studies could benefit from incorporating the spatial indices (such as SCI and aCI) that showed promise in our analysis, potentially leading to a more nuanced understanding of spatial competition dynamics in HCB formation.

5. Conclusions

This study applied hierarchical Bayesian, Gaussian process regression, and random forest models to predict LiDAR-based individual tree height to crown base (HCB) in Picea crassifolia Kom. forests. The Gaussian process model showed the best performance (RMSE = 1.5902 m), followed closely by the random forest and hierarchical Bayesian models. LiDAR-derived tree height and crown projection area and DBH were consistently identified as key factors influencing HCB across all models. The integration of spatial competition indices, such as the SCI, highlighted the role of local competition dynamics in shaping HCB patterns. While the models provide insights for precision forest management, future studies should aim to include other factors, such as genetic variability and environmental factors, to further enhance the models’ accuracy and applicability.

Author Contributions

Formal analysis, Z.Y., H.Y. and Z.Z.; methodology, G.D.; software, G.D., X.W. and H.Y.; supervision, H.Z. and Z.Y.; validation, X.W.; visualization, Z.Z.; writing—original draft, Z.Y., H.Y., X.W., G.D. and H.Z.; writing—review and editing, Z.Y. and H.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are available from the authors upon reasonable request and with permission from the Research Institute of Forest Resource Information Techniques, Chinese Academy of Forestry.

Acknowledgments

We thank the Research Institute of Forest Resource Information Techniques, the Chinese Academy of Forestry for data support. We sincerely thank the four anonymous reviewers for their valuable comments and suggestions.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Location of the sample plot.

Figure 2. Distribution of sampled trees within the sample plots.

Figure 3. Figures of hierarchical Bayesian uncertainty quantification.

Figure 4. Distribution of residuals of the HCB hierarchical Bayes model.

Figure 5. True vs. predicted values of the Gaussian mixed-effects process model, with uncertainty.

Figure 6. The predictions with uncertainty (a) and the relative uncertainty versus true values (b).

Figure 7. Predicted vs. true values, with 95% confidence intervals.

Table 1

Descriptive statistics of measurements of tree variables (hL, LiDAR-derived total tree height; caL, LiDAR-derived crown projection area; HCB, height to crown base; DBH, diameter at breast height; SD, standard deviation).

Variable	Min	Max	Mean	SD
hL (m)	4.62	22.15	13.84	3.20
caL (m2)	3.19	38.63	12.21	5.47
HCB (m)	0.90	10.20	4.80	3.52
DBH (cm)	3.60	81.10	22.57	8.54
aCI	1.36	3.10	1.93	0.34
Hegyi	0.24	11.86	2.13	1.43
Aggregation index	2.34	3.28	2.71	0.24
SCI	0.00	2.64	1.13	0.60

Table 2

Parameter estimates from the Bayesian base model.

Parameter	Mean	Std. Dev	95% CI
b 1	0.0113	0.0032	(0.0050, 0.0175)
b 2	0.0262	0.0058	(0.0149, 0.0377)
log(σ)	0.5784	0.0313	(0.5195, 0.6410)
σ	1.7840	0.0559	(1.6813, 1.8984)

Table 3

Comparisons of five spatial indices.

Rank	LOO	P Loo	Dse	Loo Scale
SCI	−1135.88	4.61	0.00	log
aCI	−1138.261	4.70	3.80	log
Aggregation index	−1352.79	145.21	8.31	log
Hegyi	−1356.54	142.81	7.55	log

Note: LOO represents the LOPCV information criterion (LOOIC), p loo represents the effective number of parameters, and Dse represents the difference in standard error.

Table 4

Performance of different random effects.

Random Effects	AIC	BIC	Log-Likelihood
(DBH)	1924.37	1949.78	−956.19
(c a L)	1934.86	1960.26	−961.43
(SCI)	2009.87	2035.27	−998.93
(DBH, caL)	1810.43	1840.07	−898.21
(DBH, SCI)	1796.77	1826.41	−891.38
(caL, SCI)	1853.03	1882.67	−919.51
(DBH, caL, SCI)	1676.60	1710.47	−830.30

Table 5

Results of random selection of the mixed-effects Gaussian process model.

Random Effects	AIC	BIC	Log-Likelihood
SCI	401.03	429.91	−189.52
caL	405.53	434.40	−191.76
caL, SCI	406.41	435.29	−1928.21
DBH, LCA	406.74	435.62	−192.37
DBH, caL, SCI	407.16	436.03	−192.58
DBH, SCI	407.29	436.16	−192.64
DBH	408.40	437.27	−193.20

Table 6

Results of hyperparameter tuning of the random forest HCB model.

mtry	Min Node Size	RMSE	MAE
2	1	1.651	1.309
2	3	1.650	1.308
2	5	1.653	1.312
3	1	1.652	1.309
3	3	1.642	1.301
3	5	1.639	1.300
4	1	1.643	1.301
4	3	1.642	1.301
4	5	1.643	1.308

Appendix A

1.

1. Hierarchical Bayesian model

• (1)

  Set up environment:

• (2)

  Define Functions:

  • calculate_metrics(y_true, y_pred): Calculate MAE and RMSE between true and predicted values.

  • mixed_effects_model(X_scaled, y, plot_index, random_effects):

    •     Define the Bayesian model:

    •     Return the fitted model and posterior samples.

  • cross_validation_by_plot(X_scaled, y, plot_index, random_effects):

    •     Loop over each unique plot:

    •     Store predictions for each plot.

    •     Return all predictions.

• (3)

  Main Analysis:

  • Load Data:

  • Define Random Effects Options: Create combinations of random effects terms for model testing.

  • Initialize an empty list to store results.

  • Loop through each random effects combination:

    •     Train and Evaluate Model:

    •     Handle any errors encountered during model fitting.

    •     Save intermediate results to an external file (optional for tracking progress).

Sort and Display Results:

• 2.

  Gaussian mixed-effects model

• (1)

  Set up environment:

• (2)

  Define Functions:

  • calculate_metrics(y_true, y_pred): Calculate MAE and RMSE between true and predicted values.

  • create_features(X, y):

    •     Generate polynomial features.

    •     Select relevant features using SelectKBest.

  • process_fold(plot, X_selected, y_scaled, plot_effects, affected_params, scaler_y):

    •     Remove one plot: Exclude all samples from the current plot as the test set.

    •     Train the model: Use data from the remaining plots to fit the model.

    •     Apply random effects to the selected parameters in the training data.

    •     Define a set of Gaussian Process kernels.

    •     Perform grid search to find the best kernel.

    •     Predict outcomes for the excluded plot using the best model.

    •     Return evaluation metrics and predictions.

  • gaussian_process_model(X, y, plot_effects, affected_params):

    •     Standardize features and response variable.

    •     Generate polynomial features and select the best.

    •     Perform cross-validation by removing one plot at a time

    •     Calculate mean evaluation metrics (MAE, RMSE) across all plots.

    •     Return the model performance metrics and predictions.

• (3)

  Main Analysis:

  • Load Data

  • Define Parameters: Specify parameter names and initialize a list for results.

  • Base Model: Fit the base Gaussian Process model without random effects and store the results.

  • Random Effects Models:

    •       Loop through all combinations of random effects:

      •          Apply selected random effects to affected parameters.

      •          Fit the Gaussian Process model using gaussian_process_model with leave-one-plot-out cross-validation.

      •          Store model name and results.

  • Sort Results: Sort models based on mean RMSE and MAE.

References

1. Maltamo, M.; Bollandsås, O.M.; Vauhkonen, J.; Breidenbach, J.; Gobakken, T.; Næsset, E. Comparing different methods for prediction of mean crown height in Norway spruce stands using airborne laser scanner data. Forestry; 2010; 83, pp. 257-268. [DOI: https://dx.doi.org/10.1093/forestry/cpq008]

2. Rais, A.; Jacobs, M.; Van de Kuilen, J.W.; Pretzsch, H. Crown structure of European beech (Fagus sylvatica): A noncausal proxy for mechanical–physical wood properties. Can. J. For. Res; 2021; 51, pp. 834-841. [DOI: https://dx.doi.org/10.1139/cjfr-2020-0382]

3. Krajnc, L.; Farrelly, N.; Harte, A.M. The influence of crown and stem characteristics on timber quality in softwoods. For. Ecol. Manag.; 2019; 435, pp. 8-17. [DOI: https://dx.doi.org/10.1016/j.foreco.2018.12.043]

4. Vanninen, P.; Mäkelä, A. Needle and stem wood production in Scots pine (Pinus sylvestris) trees of different age, size and competitive status. Tree Physiol.; 2000; 20, pp. 527-533. [DOI: https://dx.doi.org/10.1093/treephys/20.8.527] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12651433]

5. Valentine, H.T.; Gregoire, T.G.; Burkhart, H.E.; Hollinger, D.Y. A stand-level model of carbon allocation and growth, calibrated for loblolly pine. Can. J. For. Res.; 1994; 24, pp. 1575-1585. [DOI: https://dx.doi.org/10.1139/x97-033]

6. Mäkelä, A. Derivation of stem taper from the pipe theory in a carbon balance framework. Tree Physiol.; 2002; 22, pp. 891-905. [DOI: https://dx.doi.org/10.1093/treephys/22.13.891]

7. Hyytiäinen, K.; Hari, P.; Kokkila, T.; Mäkelä, A.; Tahvonen, O.; Taipale, J. Connecting a process-based forest growth model to stand-level economic optimization. Can. J. For. Res.; 2004; 34, pp. 2060-2073. [DOI: https://dx.doi.org/10.1139/x04-056]

8. Maguire, D.A.; Kershaw, J.A.; Hann, D.W. Predicting the effects of silvicultural regime on branch size and crown wood core in Douglas-fir. For. Sci.; 2005; 51, pp. 503-514. [DOI: https://dx.doi.org/10.1093/forestscience/37.5.1409]

9. Kershaw, J.A.; Maguire, D.A.; Hann, D.W. Longevity and duration of radial growth in Douglas-fir branches. Can. J. For. Res.; 2008; 20, pp. 1690-1695. [DOI: https://dx.doi.org/10.1139/x90-225]

10. Jiang, L.; Liu, R. Segmented taper equations with crown ratio and stand density for Dahurian Larch (Larix gmelinii) in northeastern China. J. For. Res.; 2011; 22, pp. 347-352. [DOI: https://dx.doi.org/10.1007/s11676-011-0178-4]

11. Nyland, R.D. Silviculture: Concepts and Applications; Waveland Press: Long Grove, IL, USA, 2016.

12. Mäkinen, H.; Isomäki, A. Thinning intensity and growth of Scots pine stands in Finland. For. Ecol. Manag.; 2004; 201, pp. 311-325. [DOI: https://dx.doi.org/10.1016/j.foreco.2004.07.016]

13. Fernández-Alonso, J.M.; Alberdi, I.; Álvarez-González, J.G.; Vega, J.A.; Cañellas, I.; Ruiz-González, A.D. Canopy fuel characteristics in relation to crown fire potential in pine stands: Analysis, modelling and classification. Eur. J. For. Res.; 2013; 132, pp. 363-377. [DOI: https://dx.doi.org/10.1007/s10342-012-0680-z]

14. Andersen, H.E.; McGaughey, R.J.; Reutebuch, S.E. Estimating forest canopy fuel parameters using LIDAR data. Remote Sens. Environ.; 2005; 94, pp. 441-449. [DOI: https://dx.doi.org/10.1016/j.rse.2004.10.013]

15. Cruz, M.G.; Alexander, M.E.; Wakimoto, R.H. Assessing canopy fuel stratum characteristics in crown fire prone fuel types of western North America. Int. J. Wildland Fire; 2003; 12, pp. 39-50. [DOI: https://dx.doi.org/10.1071/WF02024]

16. Agee, J.K.; Skinner, C.N. Basic principles of forest fuel reduction treatments. For. Ecol. Manag.; 2005; 211, pp. 83-96. [DOI: https://dx.doi.org/10.1016/j.foreco.2005.01.034]

17. Yang, Z.; Liu, Q.; Luo, P.; Ye, Q.; Duan, G.; Sharma, R.P.; Fu, L. Prediction of individual tree diameter and height to crown base using nonlinear simultaneous regression and airborne LiDAR data. Remote Sens.; 2020; 12, 2238. [DOI: https://dx.doi.org/10.3390/rs12142238]

18. Maltamo, M.; Karjalainen, T.; Repola, J.; Vauhkonen, J. Incorporating tree- and stand-level information on crown base height into multivariate forest management inventories based on airborne laser scanning. Silva Fenn.; 2018; 52, 10006. [DOI: https://dx.doi.org/10.14214/sf.10006]

19. Yan, Y.; Wang, J.; Liu, S.; Gaire, D.; Jiang, L. Modeling the influence of competition, climate, soil, and their interaction on height to crown base for Korean pine plantations in Northeast China. Eur. J. For. Res.; 2024; 9, pp. 1-4. [DOI: https://dx.doi.org/10.1007/s10342-024-01710-z]

20. Clark, J.S. Uncertainty and variability in demography and population growth: A hierarchical approach. Ecology; 2003; 84, pp. 1370-1381. [DOI: https://dx.doi.org/10.1890/0012-9658(2003)084[1370:UAVIDA]2.0.CO;2]

21. Regan, H.M.; Colyvan, M.; Burgman, M.A. A taxonomy and treatment of uncertainty for ecology and conservation biology. Ecol. Appl.; 2002; 12, pp. 618-628. [DOI: https://dx.doi.org/10.1890/1051-0761(2002)012[0618:ATATOU]2.0.CO;2]

22. Xu, C.; Gertner, G.Z. Uncertainty and sensitivity analysis for models of complex systems. Ecol. Model.; 2008; 213, pp. 385-397. [DOI: https://dx.doi.org/10.1007/978-3-540-77362-7_9]

23. Harmon, M.E.; Fasth, B.; Halpern, C.B.; Lutz, J.A. Uncertainty analysis: An evaluation metric for synthesis science. Ecosphere; 2015; 6, pp. 1-12. [DOI: https://dx.doi.org/10.1890/ES14-00235.1]

24. Reichert, P.; Omlin, M. On the usefulness of overparameterized ecological models. Ecol. Model; 1997; 95, pp. 289-299. [DOI: https://dx.doi.org/10.1016/S0304-3800(96)00043-9]

25. Ellison, A.M. Bayesian inference in ecology. Ecol. Lett.; 2004; 7, pp. 509-520. [DOI: https://dx.doi.org/10.1111/j.1461-0248.2004.00603.x]

26. McCarthy, M.A. Bayesian Methods for Ecology; Cambridge University Press: Cambridge, UK, 2007; [DOI: https://dx.doi.org/10.1017/CBO9780511802454.002]

27. Rasmussen, C.E.; Williams, C.K.I. Gaussian Processes for Machine Learning; MIT Press: Cambridge, MA, USA, 2006; pp. 79-104.

28. Chu, W.; Ghahramani, Z. Gaussian processes for ordinal regression. J. Mach. Learn. Res.; 2005; 6, pp. 1019-1041.

29. Quiñonero-Candela, J.; Rasmussen, C.E.; Williams, C.K. Approximation methods for gaussian process regression. Large-Scale Kernel Machines; MIT Press: Cambridge, MA, USA, 2007; pp. 203-224. [DOI: https://dx.doi.org/10.7551/mitpress/7496.003.0011]

30. Camps-Valls, G.; Verrelst, J.; Munoz-Mari, J.; Laparra, V.; Mateo-Jimenez, F.; Gomez-Dans, J. A survey on Gaussian processes for earth-observation data analysis: A comprehensive investigation. IEEE Geosci. Remote Sens. Mag.; 2016; 4, pp. 58-78. [DOI: https://dx.doi.org/10.1109/MGRS.2015.2510084]

31. Cressie, N.; Calder, C.A.; Clark, J.S.; Hoef, J.M.V.; Wikle, C.K. Accounting for uncertainty in ecological analysis: The strengths and limitations of hierarchical statistical modeling. Ecol. Appl.; 2009; 19, pp. 553-570. [DOI: https://dx.doi.org/10.1890/07-0744.1]

32. Gelman, A.; Carlin, J.B.; Stern, H.S.; Dunson, D.B.; Vehtari, A.; Rubin, D.B. Bayesian Data Analysis; Chapman and Hall/CRC: London, UK, 2013; [DOI: https://dx.doi.org/10.1002/wcs.72]

33. Gelman, A. Data Analysis Using Regression and Multilevel/Hierarchical Models; Cambridge University Press: Cambridge, UK, 2007; Available online: https://www.jstor.org/stable/20461880 (accessed on 12 October 2024).

34. Qian, S.S.; Cuffney, T.F.; Alameddine, I.; McMahon, G.; Reckhow, K.H. On the application of multilevel modeling in environmental and ecological studies. Ecology; 2010; 91, pp. 355-361. [DOI: https://dx.doi.org/10.1890/09-1043.1]

35. Pang, Y.; Chen, E.; Liu, Q.; Xiao, Q.; Zhong, K.; Li, X.; Ma, M. WATER: Dataset of airborne LiDAR mission at the super site in the Dayekou watershed flight zone on Jun. 23, 2008. Chinese Academy of Forestry; Institute of Remote Sensing Applications, Chinese Academy of Sciences; Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences; Heihe Plan Science Data Center: Lanzhou, China, 2008.

36. Axelsson, P. Processing of laser scanner data—Algorithms and applications. ISPRS J. Photogramm.; 1999; 54, pp. 138-147. [DOI: https://dx.doi.org/10.1016/S0924-2716(99)00008-8]

37. Liu, Q.; Fu, L.; Wang, G.; Li, S.; Li, Z.; Chen, E.; Pang, Y.; Hu, K. Improving estimation of forest canopy cover by introducing loss ratio of laser pulses using airborne LiDAR. IEEE Trans. Geosci. Remote; 2020; 58, pp. 567-584. [DOI: https://dx.doi.org/10.1109/TGRS.2019.2938017]

38. Liu, Q. Study on the Estimation Method of Forest Parameters Using Airborne Lidar. Ph.D. Thesis; Chinese Academy of Forestry: Beijing, China, 2009; (In Chinese)

39. Walters, D.K.; Hann, D.W. Taper equations for six conifer species in southwest Oregon. For. Res. Lab. Or. State Univ. Corvallis Res. Bull.; 1986; 56, 41.

40. Hegyi, F. A simulation model for managing jack-pine stands Growth Model Tree Stand Simul. R. Coll. Res. Notes; 1974; 30, pp. 74-90.

41. Clark, P.h.J.; Evans, F.C. Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology; 1954; 35, pp. 445-453. [DOI: https://dx.doi.org/10.2307/1931034]

42. Hui, G.; Hu, Y.; Zhao, Z.; Yuan, S.; Liu, W. A Forest Competition Index Based on Intersection Angle. Ann. For. Sci.; 2013; 49, pp. 68-73. (In Chinese)

43. Duvenaud, D. Automatic Model Construction with Gaussian Processes. Ph.D. Dissertation; Cambridge University Press: Cambridge, UK, 2014; [DOI: https://dx.doi.org/10.17863/CAM.14087]

44. Hensman, J.; Matthews, A.; Ghahramani, Z. Scalable variational Gaussian process classification. Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics; San Diego, CA, USA, 9–12 May 2015; Volume 38, pp. 351-360.

45. Breiman, L. Random Forests. Mach. Learn.; 2001; 45, pp. 5-32. [DOI: https://dx.doi.org/10.1023/A:1010933404324]

46. Fu, L.; Liu, Q.; Sun, H.; Wang, Q.; Li, Z.; Chen, E.; Pang, Y.; Song, X.; Wang, G. Development of a system of compatible Individual tree diameter and aboveground biomass prediction models using error-in-variable regression and airborne LiDAR Dat-a. Remote Sens.; 2018; 10, 325. [DOI: https://dx.doi.org/10.3390/rs10020325]

47. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018; Available online: https://www.R-project.org/ (accessed on 12 October 2024).

48. Sharma, M.; Parton, J. Modeling stand density effects on taper for jack pine and black spruce plantations using dimensional analysis. For. Sci.; 2009; 55, pp. 268-282. [DOI: https://dx.doi.org/10.1093/forestscience/55.3.268]

49. Rijal, B.; Weiskittel, A.R.; Kershaw, J.A. Development of height to crown base models for thirteen tree species of the North American Acadian Region. For. Chron.; 2012; 88, pp. 60-73. [DOI: https://dx.doi.org/10.5558/tfc2012-011]

50. Mäkelä, A.; Valentine, H.T. Crown ratio influences allometric scaling in trees. Ecology; 2006; 87, pp. 2967-2972. [DOI: https://dx.doi.org/10.1890/0012-9658(2006)87[2967:CRIASI]2.0.CO;2]

51. Hasenauer, H.; Monserud, R.A. A crown ratio model for Austrian forests. For. Ecol. Manag.; 1996; 84, pp. 49-60. [DOI: https://dx.doi.org/10.1016/0378-1127(96)03768-1]

52. Hooten, M.B.; Hobbs, N.T. A guide to Bayesian model selection for ecologists. Ecol. Monogr.; 2015; 85, pp. 3-28. [DOI: https://dx.doi.org/10.1890/14-0661.1]

53. Smith, J.P.; Hoffman, C.M.; Landhäusser, S.M. Modelling crown base height and crown ratio in mixed-species stands of Douglas-fir and western hemlock. For. Ecol. Manag.; 2019; 447, pp. 115-124.

54. Kneeshaw, D.; Kobe, R.K.; Coates, K.D.; Messier, C. Sapling size influences shade tolerance ranking among southern boreal tree species. J. Ecol.; 2002; 90, pp. 17-24. [DOI: https://dx.doi.org/10.1111/j.1365-2745.2005.01070.x]

55. Temesgen, H.; LeMay, V.M.; Mitchell, S.J. Tree crown ratio models for multi-species and multi-layered stands of southeastern British Columbia. For. Chron.; 2005; 81, pp. 133-141. [DOI: https://dx.doi.org/10.5558/tfc81133-1]

56. Lines, E.R.; Zavala, M.A.; Purves, D.W.; Coomes, D.A. Predictable changes in aboveground allometry of trees along gradients of temperature, aridity and competition. Glob. Ecol. Biogeogr.; 2012; 21, pp. 1017-1028. [DOI: https://dx.doi.org/10.1111/j.1466-8238.2011.00746.x]

57. Xing, Z.; Bourque, C.P.A.; Swift, D.E.; Clowater, C.W.; Krasowski, M.; Meng, F.R. Carbon and biomass partitioning in balsam fir (Abies balsamea). Tree Physiol.; 2005; 25, pp. 1207-1217. [DOI: https://dx.doi.org/10.1093/treephys/25.9.1207]