Modelling Diameter at Breast Height Distribution for Eight Commercial Species in Natural-Origin Mixed Forests of Ontario, Canada

Abstract

Diameter at breast height (DBH) is a unique attribute used to characterize forest growth and development for forest management planning and to understand forest ecology. Forest managers require an array of DBHs of forest stands, which can be reconstructed using selected probability distribution functions (PDFs). However, there is a lack of practices that fit PDFs of sub-dominating species grown in natural mixed forests. This study aimed to fit PDFs and develop predictive models for PDF parameters, so that the predicted distribution would represent dynamic forest structures and compositions in mixed forest stands. We fitted three of the simplest forms of PDFs, log-normal, gamma, and Weibull, for the DBH of eight tree species, namely balsam fir (Abies balsamea [L.] Mill.), eastern white pine (Pinus strobus L.), paper birch (Betula papyrifera Marshall), red maple (Acer rubrum L.), red pine (Pinus resinosa Aiton), sugar maple (Acer saccharum Marshall), trembling aspen (Populus tremuloides Michx), and white spruce (Picea glauca [Moench] Voss), all grown in natural-origin mixed forests in Ontario province, Canada. We estimated the parameters of the PDFs as a function of DBH mean and standard deviation for these species. Our results showed that log-normal fit the best among the three PDFs. We demonstrated that the predictive model could estimate the recovered parameters unbiasedly for all species, which can be used to reconstruct the DBH distributions of these tree species. In addition to prediction, the cross-validated R² for the DBH mean ranged between 0.76 for red maple and 0.92 for red pine. However, the R² for the regression of the standard deviation ranged between 0.00 for red pine and 0.69 for sugar maple, although it produced unbiased predictions and a small mean absolute bias. As these mean and standard deviations are regressed with dynamic covariates (such as stem density and stand basal area), in addition to climate and static geographic variables, the predicted DBH distribution can reflect change over time in response to management or any type of disturbance in the regime of the given geography. The predictive model-based DBH distributions can be applied to the design of appropriate silviculture systems for forest management planning.

1. Introduction

Tree diameter at breast height (DBH) is the fundamental attribute used to characterize individual trees and stand-level forest growth and development. It has always been used in forest management planning [1,2] to predict other tree- and stand-level attributes, such as tree height, tree forms, stem density, basal area, stand structure and composition, and wood biomass and volume [3]. Model-based prediction systems for assessing forest resources require DBH at the foremost, as it is an input variable used in growth model simulators, which are used in management planning and implementation. The Forest Vegetation Simulator (FVS) [4], Mixedwood Growth Model [5], NATURA [6], and the Tree and Stand Simulator [7] are a few examples of growth simulators used in various parts of North America that rely on DBH. These simulators are used to characterize forest stand growth and predict long-term timber volumes, stand density, and basal area [8,9] as needed for forest management.

Silviculture practices in commercially managed forest largely depend on stand-level basal area and stem density to regulate and direct the succession trajectory throughout the rotation cycle of forest stands. The stand density management diagram [10,11] and the stocking diagram [12,13] are examples of decision tools for regulating silvicultural systems. The proper maintenance of DBH (class) distribution, stand age, structure, composition, relative frequency (or proportion) [14], and gap dynamics [15] are the basis for regulating forests over rotation cycles.

The DBH distribution of naturally grown forests can serve as a reference for management. For example, it helps minimize the differences of age classes in frequency distributions, while applying silviculture systems and harvest schedules. Management requires accounting for gap dynamics [15] for competition-related mortality, natural senescence, and other forms of natural disturbance like wildfire, insect epidemics, and disease [16,17], in order to mimic the natural process of stand development. An appropriate silvicultural strategy can convert evenly aged forest stands into close-to-naturally grown unevenly aged structures [18] using stem density manipulation over time. Growth models provide a useful decision support tool in long-term forest management, so more precise predictive systems are beneficial [19].

An array of individual tree DBHs may follow some (standard) probability distribution functions (PDFs) depending on various environmental variables, such as the following dendrometric variables: stem density, basal area, dominant height, and species composition, including age structure. Geo- and topographic variables like elevation, slope, aspect, geographic position, and climate, particularly precipitation and temperature, among others, also affect species’ DBH distribution [20,21]. Some of these variables are static (e.g., geographic position) but many are dynamic (e.g., vegetation structure and bioclimate). Therefore, the parameter of the PDFs, once fitted, may change over time, which may be predicted by developing some predictive models [22,23,24]. Alternatively, the parameters of the probability density function can be mathematically calculated by solving the system of equations of moments and percentiles [22,25,26], which can be regressed in terms of dynamic variables.

Several tree species in the province of Ontario, Canada are commercially important. Although black spruce (Picea mariana (Mill) B.S.P) and jack pine (Pinus banksiana Lamb) are dominating species in terms of presence and abundancy across the province, 82 species are listed in provincial inventory data, of which 42 occur in our current study site in various proportions [27]. Among these, a few have a substantial presence, depending on location and ecological distribution. These species have two sources of origin: naturally grown, mainly by wildfire disturbance, which encompasses about 75% of the sample plots or plantation origin. Naturally originated stands are richer in diversity in DBH, height, and age than the plantation origin stands [28,29].

There have been a few studies on the development of diameter distribution models for jack pine and black spruce grown in Ontario [22,30,31,32]. However, these models are not available for many other commercially important tree species which are less common or sub-dominating in terms of abundancy in Ontario. The study scoped that the developed model would be used in need-based forest management planning, along with being valued for use in the scientific literature for DBH simulations using simple statistics-based parameter recovery techniques. The goal of this study was to develop dynamic models for DBH distributions of eight commercial tree species grown in natural-origin mixed stands across Ontario, namely balsam fir (Abies balsamea (L.) Mill.), eastern white pine (Pinus strobus L.), paper birch (Betula papyrifera Marshall), red maple (Acer rubrum L.), red pine (Pinus resinosa Aiton), sugar maple (Acer saccharum Marshall), trembling aspen (Populus tremuloides Michx), and white spruce (Picea glauca [Moench] Voss). The specific objectives were to: (a) fit selected three PDFs, gamma, log-normal, and Weibull, for the DBHs of these species, (b) compare DBH distributions reconstructed based on the parameters obtained using maximum likelihood and parameter recovery methods, and (c) construct predictive models for the estimated parameters as functions of stand-level dynamic variables, using data collected from permanent sample plots across Ontario.

2. Methods

2.1. Description of Study Area

Our study area is spread across Ontario, Canada. It expands between 75° W and 95° W longitude and 43° N and 53° N latitude. All sample plots taken in this study are located in the Ontario Shield and Mixedwood Plains ecozones (Figure 1).

The area features wildfire-origin coniferous pure or mixed forest stands in the north and windthrow and insect epidemics (especially, spruce budworm—Choristoneura fumiferana [Clemens])-origin mixed species in the south [33,34]. Black spruce, jack pine, balsam fir, and tamarack (Larix laricina) are dominating coniferous species in the northern region, whereas deciduous mixtures of sugar maple and American beech (Fagus grandifolia) are found in the southern region of this ecozone. Although the structural compositions of the vegetation vary, black spruce and jack pine are universal across this ecozone [35]. Other species can be found in various proportions.

2.2. Data Source

Data collected from permanent sample plots, which were established and repeatedly measured by the Ontario Ministry of Natural Resources and Forestry, were used in our study. These plots were established and maintained in the natural and planted forest stands, comprising different forest types in terms of ecology, geography, and species composition [36]. Most of the sample plots have repeated measurements, but the number of repetitions varied between two and five. However, a few sample plots had one time measurement. Among these, we used DBH data selected from live trees grown only in natural mixed stands.

The tree species (other than black spruce and jack pine) that accounted for at least 20% of the basal area in each sample plot was selected for this study. That threshold resulted in eight species with varying numbers of sample plots, repeatedly measured plots, and individual trees (

Table 1. Species included in this study and sample numbers for two plot types and individual trees. The species codes in this table are used throughout the paper to refer to the respective species.

Species Code	Species	Botanical Name	Sample Number
			Spatial Plot	Repeated Plot	Individual Trees
BF	Balsam fir	Abies balsamea (L.) Mill.	857	1433	47,597
EWP	Eastern white pine	Pinus strobus L.	139	472	8409
PB	Paper birch	Betula papyrifera Marshall	783	1149	23,033
RM	Red maple	Acer rubrum L.	176	454	9153
RP	Red pine	Pinus resinosa Aiton	56	165	2865
SM	Sugar maple	Acer saccharum Marshall	136	578	19,274
TA	Trembling aspen	Populus tremuloides Michx.	881	1211	43,953
WS	White spruce	Picea glauca (Moench) Voss	278	423	4812

). The mean proportions of the selected species varied between 27% and 71%, but in some cases were 100%, signifying a single-species sample plot.

To avoid the non-convergence problem when fitting probability distribution functions, we selected plots for the species that included at least five individual trees’ DBH measurements of the same species individually. Summary statistics of tree-level DBHs are presented in

Table 2. Summary statistics of diameters at breast height of selected eight species from the selected sample plots.

Species	Min. (cm)	Median (cm)	Mean (cm)	Max (cm)	Std. Dev (cm)
Balsam fir	1.00	5.00	5.46	57.50	4.20
Eastern white pine	2.50	11.00	16.34	88.50	13.98
Paper birch	2.50	6.30	8.16	58.90	5.69
Red pine	2.50	22.60	22.04	57.30	2.31
Red maple	1.30	4.40	5.95	51.40	4.49
Sugar maple	2.50	7.60	12.09	83.50	10.96
Trembling aspen	2.50	7.30	9.68	64.80	7.15
White spruce	2.50	8.00	11.03	57.20	8.43

.

2.3. Model Development

We fit and evaluated two types of models. The first was to fit the probability distribution functions (PDFs) of the DBHs of individual trees of eight selected species. To fit the PDFs, we used the entire individual trees for global fit and, for each repeated plot, DBH for a plot (stand)-level fit, assuming the repeated data were independent (Table 1). The second was to construct two regression models; one for the DBH’s mean and another for standard deviation, from which we could recover parameters of the PDFs. We constructed the regression fitting at the spatial plot level by considering a hierarchical random effect for the repeated measurements within the spatial plot.

2.3.1. Fitting Probability Distribution Functions

The parameters of any PDFs may be estimated using distribution fitting that involves the maximum likelihood method. Alternatively, these parameters can be estimated or calculated using simple statistics like the mean and standard deviation, which is often referred to as the parameter recovery method [25,37]. This method has been used in studies by several researchers such as Liu et al. [22,23] and Zhang et al. [38], among others. In this method, the parameters of PDFs are functionally linked to the stand-level mean and standard deviation of DBH [22,23]. We reconstructed DBH distributions for eight species (Table 1) based on the parameters estimated using both methods and constructed 5 cm DBH classes. The mathematical forms of selected PDFs, their relationships with the mean and standard deviation of the data (DBH), and their respective parameters are described in the following subsections.

We fit three PDFs, namely log-normal [39] adopted from [24], gamma [40], and Weibull [41,42,43], which are suitable for fitting positive and continuous random variables. Although there are several types of PDFs suitable for these types of variables that have been used and tested for selected tree species (see [24,26,44,45]), we selected the simplest that could readily be applicable for forest management for all species considered in this study.

Log-Normal Probability Distribution Function

Let us suppose X is a random variable (DBH in our case) with positive values (non-negative). The log-transformed variable (log(X) follows a normal distribution with a mean (µ > 0) and standard deviation (σ > 0), provided the variable X follows a log-normal distribution. Here, µ and σ are the mean and standard deviation of the log-transformed random variable log(X). The PDF (f) of X is as given in Equation (1) below. For simplicity, Equations (1)–(9) are reproduced here from Rijal and Sharma [32].

$$f\left(X|\mu ,{\sigma }^2\right)=\frac{1}{x\sigma \sqrt{2\pi }}\exp [\frac{{\left(\mathit{log}\left(x\right)-\mu \right)}^2}{2{\sigma }^2}]$$

(1)

We estimate the parameters µ and σ using the maximum likelihood method. The sample mean ($\overline{x}$) and standard deviation (s) of the random variable (X) can be used to calculate the parameters (µ: population mean and σ: standard deviation), as they are unbiased estimators, using the following simultaneous Equations (2) and (3):

$$\mu =\log (\frac{\overline{x}}{\sqrt{\left(\frac{s^2}{{\overline{x}}^2}+1\right)}})$$

(2)

$$\sigma =\sqrt{log\left(\frac{s^2}{{\overline{x}}^2}+1\right)}$$

(3)

The “fitdistrplus” package [46] in R [47] was used to fit Equation (1) and estimate the parameters for the available DBH data values. In the parameter recovery method, these parameters were calculated by solving two simultaneous equations (Equations (2) and (3)), relating the parameters with the sample mean ($\overline{x}$) and sample variance s². The “nleqslv” package [48] in R was used to solve the system of simultaneous equations (Equations (2) and (3)).

Weibull Probability Distribution Function

Assume X is a positive and continuous random variable. The PDF of the Weibull distribution with scale parameter α (>0) and shape parameter β (>0) is given in Equation (4), as follows:

$$f\left(X,\alpha ,\beta \right)=\frac{\beta }{\alpha }{\left(\frac{x}{\alpha }\right)}^{\beta -1}\ast \mathrm{e}\mathrm{x}\mathrm{p}{(-\frac{x}{\alpha })}^{\beta }$$

(4)

The parameters α and β are estimated using the maximum likelihood method. We used the “fitdistrplus” package to estimate these parameters, as above in Equation (1), for an array of the DBH. Alternatively, these parameters can be calculated using the sample mean ($\overline{x}$) and sample standard deviation (s) of an array of DBHs by solving simultaneous Equations (5) and (6), as given below:

$$\mu =\alpha \mathsf{\Gamma }(1+\frac{1}{\beta })$$

(5)

$${\sigma }^2={\alpha }^2\mathsf{\Gamma }\left(1+\frac{2}{\beta }\right)-{\mu }^2$$

(6)

where Γ is the gamma function. Please note that the sample mean and standard deviation are the unbiased estimators of population mean (µ) and standard deviation ($\sigma$).

As described in the section “Log-Normal Probability Distribution Function”, parameters α and β were found as the solution values of Equation (4) using the maximum likelihood method. The parameters in Equations (5) and (6) were calculated as the solution of two simultaneous equations. Like above, the “nleqslv” package was used to solve the system of simultaneous equations (Equations (5) and (6)).

Gamma Probability Distribution Function

Let us suppose X is a positive continuous random variable (>0 to $\infty$). The PDF of gamma distribution of random variable X with the (shape) parameter β (>0) and (rate) parameter λ (>0) is given by Equation (7), as follows:

$$f\left(x\right)=\frac{{\lambda }^{\beta }}{\Gamma \beta }{x}^{(\beta -1)}{e}^{-\lambda x}$$

(7)

$$\mu =\frac{\beta }{\lambda }$$

(8)

$${\sigma }^2=\frac{\beta }{{\lambda }^2}$$

(9)

Once again, as described in the sections “Log-normal probability distribution function” and “Weibull probability distribution function” above, parameters β and λ were estimated by solving Equation (7) using the “fitdistrplus” package. The simultaneous Equations (8) and (9) were solved to find parameters β and λ, using the sample mean and standard deviation as an unbiased estimator of the population mean (μ) and standard deviation (σ). As above, the simultaneous equations were solved using the “nleqslv” package.

2.3.2. Development of Regression Models

The distribution of DBH depends on various environmental variables, including dendrometric variables such as competition-related attributes like basal area, stem density, dominant height, species composition [22,49], and age structure; geographic variables such as elevation, slope, aspect, geographic position, and climate [20,21,49], chiefly precipitation and temperature, among others. These are either static or dynamic variables. The dendrometric variables are dynamic as they change over time in response to management interventions or natural growth and disturbances, whereas the geographic variables are considered static as they do not change within the span of a forest rotation. It indicates that the parameters of the PDFs may change over time and can be predicted by developing the models based on the dendrometric variables like site index, stand density, age, basal area, or quadratic mean diameter (QMD) [22,23,24], in addition to static variables like geographic attributes.

To construct predictive models, we require model covariates. From an array of measured DBHs of species-wise individual trees listed in each sample plot, we computed the following stand-level dendrometric variables: species-wise and total basal area per ha (BAPH), stem density per ha, quadratic mean diameter, and dominant height for each sample plot (

Table 3. Summary statistics of stand-level derived variables for eight selected species found in sample plots of this study.

Species	Min.	Median	Mean	Max	Std. Dev
Species-wise density (No./ha)
Balsam fir	55.56	500.00	793.17	9975.00	872.06
Eastern white pine	125.00	300.00	445.39	3800.00	436.48
Paper birch	61.77	300.00	502.23	4150.00	515.04
Red pine	125.00	325.00	434.09	2375.00	342.42
Red maple	125.00	375.00	504.02	2975.00	425.59
Sugar maple	125.00	775.00	833.65	2900.00	456.37
Trembling aspen	58.59	450.00	890.21	9475.00	1250.45
White spruce	48.83	200.00	269.90	1800.00	188.68
Total stand density (No./ha)
Balsam fir	200.00	2100.00	2459.68	11,250.00	1435.08
Eastern white pine	175.00	1800.00	1932.20	6800.00	986.06
Paper birch	361.45	2350.00	2763.22	12,375.00	1689.96
Red pine	425.00	2000.00	2078.33	6075.00	824.39
Red maple	425.00	1900.00	2107.27	7775.00	1069.05
Sugar maple	175.00	1200.00	1329.24	4000.00	632.46
Trembling aspen	175.00	2500.00	2853.84	17,625.00	1837.25
White spruce	225.00	2025.00	2343.60	11,250.00	1445.39
Species-wise basal area (m2/ha)
Balsam fir	0.05	2.06	3.67	33.71	4.49
Eastern white pine	0.09	15.32	16.17	50.91	10.15
Paper birch	0.07	1.96	3.83	31.05	4.82
Red pine	0.29	21.67	21.73	50.36	12.69
Red maple	0.10	1.44	2.20	16.71	2.37
Sugar maple	0.09	18.77	17.44	43.52	9.81
Trembling aspen	0.09	6.98	9.95	45.13	9.29
White spruce	0.11	2.51	3.89	39.07	4.61
Stand level basal area (m2/ha)
Basam fir	0.42	26.65	26.83	67.95	11.37
Eastern white pine	0.93	32.42	32.07	67.95	11.49
Paper birch	0.87	24.62	24.95	67.95	10.27
Red pine	1.32	38.20	37.56	67.95	12.08
Red maple	4.23	29.70	30.07	61.08	10.10
Sugar maple	5.48	26.86	27.01	51.52	6.94
Trembling aspen	0.44	25.04	25.29	61.79	9.95
White spruce	6.13	28.86	29.15	65.23	10.37

). Additionally, we received simulated climatic attributes such as annual and monthly mean precipitation, temperature, number of growing days, and moisture index from Natural Resources Canada, Great Lake Forestry Centre, Ontario, for each sample plot. We calculated the overall annual mean and standard deviation of these variables and used them as model covariates. Some variables were log-transformed to eliminate the possible effects of multicollinearity among the covariates considered in our full models, especially basal area, stem density, and total and species-wise attributes for the given stands. The formulation of the regressions is given in Equation (10) and detailed with the examined covariates in Equation (11).

We fit the two models, one for the DBH mean and another for standard deviation, as a function of stand-level covariates. These response variables (mean and standard deviation) were used to compute the parameters of each PDF using the equations aforementioned (Equations (2), (3), (5), (6), (8) and (9)) and then simulate the DBH distributions using selected three PDFs.

Let $Y$ be a random variable (DBH) representing one of our response variables (the mean of DBH and standard deviation of DBH of (each) sample plot i).

The expected value of $Y_i$ is subjected to the random effect of the sample plot i. The expected values (E(Y_i)) are regressed with model covariates using a generalized linear modelling technique [47] with m (=1, 2, 3, …) covariates using link function g(.) [50], as follows:

$$g\left(E\right(Y_i\left)\right)={(\beta }_0+{\phi }_i+{\phi }_{ij})+{\beta }_m{x}_{im}$$

(10)

where β₀ is the intercept parameter for fixed effects, φ_i is the plot level random effect, and φ_ij is the year-level random effect nested within the plot. The random effects were considered only in the intercept portion of the parameters model for simplicity, and β_m (m = 1, 2, 3, …) are parameters associated with other model covariates, as presented in Equation (11).

Our model response variables, the mean and standard deviation, are positive continuous numbers. Therefore, the parameters of these models were estimated assuming the gamma-distributed random errors, e.g., [22], instead of the traditional practice of assuming the normality of the errors.

Foremost, we fit full models that consisted of all dendrometric, geographic, and climatic attributes as model covariates, as mentioned above, for these two models, one for DBH mean and another for DBH standard deviation, for all eight species. We sorted out the model covariates to construct the reduced models using Bayesian information criteria (BIC) [51] and removed the covariates that had the highest non-significant (>0.05) p-value systematically one at a time. We were cautious of the possible impacts of multicollinearity on the estimated parameters. We calculated the pairwise variance inflation factor (VIF) among the covariates that we considered in our full model. We kept the threshold value of VIF less than 5 (or Pearson’s correlation coefficient of about 0.90), following the thumb rule [52]. In addition to a high VIF, the impact of multicollinearity could be observed through the inconsistent signs of estimated parameters and/or high t-values associated with model variables [53].

Our preliminary analyses showed that adding a hierarchal level for repeated measures, i.e., random effect (φ_ij), did not improve the model. Therefore, these repeated sample plots were also considered to be independent samples in our study. The mathematical formulation of the regression model is as follows:

$$\begin{array}{cc}g\left(E\left(Y_{ik}|{\mu }_i\right)\right)=& {(\beta }_0+{\phi }_i)+{\beta }_{1}Tot\_BA+{\beta }_{2}Species\_BA+{\beta }_3\log (Tot\_Density)\hfill \\ & +{\beta }_{4}Species\_Density+{\beta }_{5}Tot\_QMD+{\beta }_{6}Species\_QMD\hfill \\ & +{\beta }_{7j}ClimateVariables\hfill \end{array}$$

(11)

where Species_BA and Tot_BA are species-wise and total basal area (m²ha⁻¹), respectively. Likewise, density is the stand-level stem density (no. ha⁻¹) and QMD is the quadratic mean diameter (cm). The model includes a set of climate variables (j). The random effects were applied only to the plot level in our reduced model. The glmmTMB package [54] in R was used to fit the equations. We fit Equation (11) for the DBH mean and standard deviation separately. When we fit for the mean, the µ represents the mean of the DBH mean, whereas when we fit for standard deviation, the notation µ represents the mean of standard deviation.

2.4. Model Performances

The foremost statistical tool we used to compare the PDF was the BIC [51]. We simulated DBH distribution using the parameters estimated, using both approaches for the same number of individual trees for all species separately, as we had in the original dataset for each sample plot (Table 1). To predict consistent DBH distributions, we repeated the simulations 1000 times and calculated the mean of the model-predicted DBH distributions of these 1000 repetitions. We used the same random seed at the start of each species and PDF to make the results comparable while reconstructing the simulation-based DBH distributions. We compared the simulated DBH with the observed data. In addition, DBH class intervals of five centimeters (<5 cm, 5–10 cm, …, 45–50 cm) were constructed for the observed data array and the simulated values. We used the terms “fitted parameter” for the solutions obtained from distribution fitting and “recovered parameters” for those obtained from solving simultaneous equations. A Chi-squared test was performed to evaluate the hypothesis of whether the simulated and observed DBH were significantly different or not. In addition, we calculated the mean absolute bias (MAB), mean bias (MB), and R² from the mean of 1000 repeated simulations and observed DBHs. We relied on the differences in the graphical presentation, as well as on other statistics.

We applied a cross-validation technique to compare the model performance statistics. For this, a one-fold 20% leave-out method was applied, and the process was repeated 200 times [55,56]. These 200 iterations produced the distribution of MB, MAB, and R², but, for simplicity, we present the median and 90% confidence band of these statistics.

3. Results

3.1. Distribution Fitting

Firstly, we selected three probability distribution functions, which were fitted for eight species separately, using all trees from all sample plots to obtain a landscape level (Table 1) model for each species. Comparatively, these PDFs fit well with the observed data across the data range, as shown in Figure 2.

However, the BIC values showed that the log-normal distribution fit best for seven species, excluding red pine, which seemed to have a bimodal or multimodal distribution that differed from the other seven species, which had reverse j-shaped curves (Figure 2). Using the sample plot-wise large dataset, the model-estimated parameters were used to reconstruct the DBH distributions at the forest landscape level (

Table 4. Fitted and recovered parameters of species-wise three probability distribution functions (Equations (1)–(9)) at landscape level (all sample plots) for eight species. The species are balsam fir (BF), eastern white pine (EWP), paper birch (PB), red maple (RM), red pine (RP), sugar maple (SM), trembling aspen (TA), and white spruce (WS).

Probability Distribution	Estimation Method	Parameter	Species Code
			BF	EWP	PB	RM	RP	SM	TA	WS
Gamma	Fitted	Shape	3.232	1.411	2.611	2.947	2.382	1.540	2.387	2.011
		Rate	0.501	0.086	0.320	0.495	0.108	0.127	0.247	0.182
	Recovered	Shape	2.362	1.367	2.056	1.759	3.204	1.217	1.830	1.710
		Rate	0.366	0.084	0.252	0.296	0.145	0.101	0.189	0.155
Log-normal	Fitted	Mean	1.703	2.399	1.895	1.604	2.869	2.134	2.046	2.132
		Standard deviation	0.545	0.922	0.625	0.552	0.755	0.838	0.655	0.732
	Recovered	Mean	1.689	2.519	1.901	1.559	2.957	2.193	2.052	2.170
		Standard deviation	0.542	0.676	0.575	0.612	0.476	0.707	0.603	0.619
Weibull	Fitted	Shape	1.701	1.200	1.580	1.542	1.816	1.212	1.503	1.417
		Scale	7.308	17.437	9.172	6.693	24.729	12.984	10.830	12.225
	Recovered	Shape	1.595	1.185	1.479	1.359	1.882	1.112	1.388	1.338
		Scale	7.200	17.313	9.018	6.497	24.832	12.578	10.604	12.007

) using any of the three PDFs.

We then fitted these three PDFs of the individual tree DBHs for each species for each sample plot. The number of sample plots varied from 165 (for red pine) to 1433 (for balsam fir). The sample plot-wise distributions fits for all species using these three PDFs seemed to be alike in graphical representation (Figure 3a,b). However, the values of BIC showed that the log-normal distribution fit best in most of the sample plots (for all except red pine;

Table 5. Number of sample plots with best fit for three probability distribution fittings based on Bayesian information criteria for eight species taken in this study.

Species	Total No. of Sample Plot	Number of Superior Plots
		Gamma	Log-Normal	Weibull
Balsam fir	1433	168	1005 (70%)	260
Eastern white pine	472	45	231 (49%)	196
Paper birch	1149	158	701 (61%)	290
Red pine	165	11	38	116 (70%)
Red maple	454	26	352 (77%)	76
Sugar maple	578	55	408 (71%)	115
Trembling aspen	1211	233	514 (42%)	464
White spruce	423	37	236 (56%)	150

).

We compared the recovered and fitted parameters of these three PDFs. The differences between these two for the landscape level fit were not substantial (Table 3). Our analysis results for sample plot-wise model fits and recoveries for all cases (

Table 6. The median and 90% confidence band of fitted and recovery parameter distributions of the diameters at breast height of eight species for all sample plots fitted individually. The species are balsam fir (BF), eastern white pine (EWP), paper birch (PB), red maple (RM), red pine (RP), sugar maple (SM), trembling aspen (TA), and white spruce (WS).

Species (Code)	Statistics	Gamma				Log-Normal				Weibull
		Fitted		Recovered		Fitted		Recovered		Fitted		Recovered
		Shape	Rate	Shape	Rate	Mean	Std.Dev	Mean	Std.Dev	Shape	Scale	Shape	Scale
BF	Median	5.145	0.817	4.493	0.685	1.666	0.444	1.660	0.436	5.145	0.817	2.261	6.661
	5%	2.246	0.232	1.583	0.197	1.263	0.225	1.248	0.225	2.246	0.232	1.283	4.022
	95%	19.821	4.759	17.336	4.254	2.357	0.694	2.389	0.660	19.821	4.759	4.707	13.692
EWP	Median	3.458	0.210	3.689	0.191	2.797	0.567	2.888	0.470	3.484	0.210	2.031	22.488
	5%	1.059	0.058	0.686	0.053	1.582	0.240	1.542	0.225	1.056	0.058	0.815	6.020
	95%	18.062	1.182	16.426	1.047	3.560	1.057	3.565	0.933	19.032	1.190	4.571	42.043
PB	Median	7.170	0.923	6.564	0.823	1.872	0.381	1.884	0.363	7.170	0.923	2.778	8.211
	5%	2.306	0.228	1.713	0.185	1.233	0.175	1.226	0.178	2.306	0.228	1.340	3.899
	95%	31.933	5.979	26.977	5.248	2.888	0.676	2.897	0.648	31.933	5.979	5.985	21.192
RM	Median	5.840	0.963	4.755	0.772	1.576	0.414	1.564	0.419	5.840	0.963	2.332	6.082
	5%	1.905	0.210	1.146	0.151	1.187	0.196	1.165	0.202	1.905	0.210	1.077	3.736
	95%	25.406	6.100	20.223	5.073	2.504	0.714	2.529	0.759	25.406	6.100	5.118	15.704
RP	Median	8.852	0.374	8.723	0.366	3.167	0.358	3.176	0.306	8.852	0.374	19.501	25.871
	5%	1.920	0.105	1.906	0.114	2.109	0.130	2.141	0.121	1.920	0.105	6.268	9.518
	95%	62.278	2.444	57.955	2.191	3.678	0.804	3.677	0.631	62.278	2.444	32.341	40.794
SM	Median	2.076	0.166	1.788	0.133	2.179	0.716	2.234	0.653	2.076	0.166	1.371	12.643
	5%	1.145	0.075	0.665	0.061	1.553	0.316	1.492	0.291	1.145	0.075	0.801	5.756
	95%	11.206	1.218	10.621	1.076	3.029	1.001	3.059	0.945	11.206	1.218	3.607	26.087
TA	Median	10.177	1.003	9.985	0.964	2.284	0.320	2.298	0.301	10.177	1.003	3.489	11.728
	5%	2.188	0.157	2.040	0.158	1.357	0.142	1.340	0.141	2.188	0.157	1.473	4.523
	95%	49.690	4.007	43.706	3.700	3.359	0.718	3.359	0.608	49.690	4.007	7.777	31.540
WS	Median	11.728	0.435	4.310	0.397	2.114	0.479	2.129	0.431	4.811	0.435	2.211	11.077
	5%	4.523	0.109	1.172	0.095	1.352	0.231	1.334	0.231	1.502	0.109	1.090	4.592
	95%	31.540	3.461	17.105	3.105	3.059	0.863	3.066	0.741	18.689	3.461	4.672	25.452

) showed that we can precisely estimate the parameters of these distributions using simple and widely available variables, namely the computed mean and standard deviation of DBH. Depending on individual sample plots and the size of the sample, DBH distributions varied across species based on the median and 90% percentiles (Table 6).

Thereafter, we grouped the continuous values of observed and simulated DBHs into class intervals of 5 cm and reconstructed the DBH class-wise frequency distribution. The class interval of simulated DBH and observed DBH applying the three PDFs (log-normal, gamma, and Weibull) and recovered and fitted parameters showed that they were alike in general. The p values of the Chi-squared test ranged from 0.22 to 0.99 for all PDFs for all species, except for the Weibull distribution for red pine (p = 0.01). The p values greater than our threshold, 0.05, suggested that the distribution which we selected can represent the observed values that correspond to the ones in Figure 3.

3.2. Fitting Regression Models

We constructed predictive models for the mean and standard deviation DBH (Equation (11)) in terms of dendrometric and environmental covariates. The covariates had different levels of impact and significance in predicting the mean (

Table 7. (a) Estimated parameters of species-wise regression fit for (response variable) stand-level mean with environmental variables. The species are balsam fir (BF), eastern white pine (EWP), paper birch (PB), red maple (RM), red pine (RP), sugar maple (SM), trembling aspen (TA), and white spruce (WS). The numbers presented in the table are significant values at a significance level of 0.05. Parameter details with t and p values are presented in Table A1. (b) Estimated parameters of species-wise regression fit for (response variable) stand-level standard deviation with environmental variables. The species are balsam fir (BF), eastern white pine (EWP), paper birch (PB), red maple (RM), red pine (RP), sugar maple (SM), trembling aspen (TA), and white spruce (WS). Numbers presented in the table are significant values at the probability significance level of 0.05. The details of the parameters with t and p values are presented in Table A2.

(a) Estimates for Species—Mean
Parameters	BF	EWP	PB	RM	RP	SM	TA	WS
Intercept	-	-	-	-	-	-	-	-
BAPH Total	-	-	−0.0023	−0.0041	−0.0120	−0.0081	0.0038	-
BAPH Species	-	-	-	-	0.0095	-	-	−0.0102
log (Density Total)	0.0291	0.0753	0.0728	0.0819	0.2445	0.0558	0.0934	0.0741
Density Species (Respective)	-	−0.0001	0.0000	0.0000	−0.0003	-	0.0000	0.0504
log (QMD Total)	-	0.1769	0.2710	0.1807	0.7400	0.2809	0.4125	0.1623
QMD Species	0.1143	0.0515	0.0897	0.0998	0.0352	0.0645	0.0554	0.0805
Dominant Height	-	-	-	0.0039	-	-	-	0.0042
Longitude	0.0070	0.0569	-	-	0.0193	-	-	-
Latitude	0.0396	-	-	-	-	-	-	-
Elevation	-	-	-	-	-	-	-	-
Temperature: Mean	0.0675	0.0832	-	-	-	-	0.0215	-
Temperature: Standard Deviation	-	3.7628	-	-	-	-	0.4068	-
No. of Growing Days: Mean	−0.0048	-	-	-	-	-	-	-
Precipitation: Standard Deviation	−0.0004	0.0075	-	-	0.0036	-	-	-
No. of Growing Days: Standard Deviation	0.0130	-	-	-	-	-	−0.0501	-
Moisture Index: Mean	-	0.2077	-	-	-	0.0025	-	-
Random effect	0.0608	0.1588	0.1124	0.1575	0.0905	0.1148	0.1510	0.1086
Residual standard deviation	0.0606	0.0819	0.0661	0.0817	0.0553	0.0790	0.0776	0.0755
(b) Estimates for Species—Standard Deviation
Parameters	BF	EWP	PB	RM	RP	SM	TA	WS
Intercept	-	−14.2566	-	19.6753	−0.0226	−5.0884	−5.1396	−1.8670
BAPH Total	-	−0.0211	−0.0092	-	-	-	−0.0183	-
BAPH Species	−0.0210	-	−0.0227	-	-	-	-	−0.0489
log (Density Total)	-	0.7612	0.1621	-	0.9856	-	0.4854	0.2585
Density Species (Respective)	0.0001	-	0.0002	0.0708	−0.1343	0.0002	-	0.1443
log (QMD Total)	-	1.4041	0.6792	0.2721	1.7694	−0.2352	1.0895	0.3235
QMD Species	0.1743	0.0283	0.1049	0.1335	-	0.0844	0.0457	0.1198
Dominant Height	0.0037	-	−0.0081	-	-	-	0.0113	-
Longitude	-	-	-	−0.0447	-	-	-	-
Latitude	−0.0167	-	-	−0.4636	−0.1843	0.2268	-	-
Elevation	-	-	-	−0.0030	-	-	-	-
Temperature: Mean	−0.0326	-	-	−0.4820	-	0.5496	-	-
Temperature: Standard Deviation	-	3.2186	-	-	-	-	0.8155	−1.5078
Precipitation: Mean	-	-	-	-	-	0.0007	-	-
No. of Growing Days: Mean	-	0.0221	−0.0082	-	-	−0.0380	-	-
Precipitation: Standard Deviation	-	-	-	-	-	-	−0.0026	-
No. of Growing Days: Standard Deviation	-	-	-	-	-	-	−0.0774	-
Moisture Index: Mean	0.1264	-	−0.2333	-	-	-	-	-
Random effect	0.1718	0.3392	0.3237	0.2811	0.2814	0.2620	0.3525	0.3174
Residual standard deviation	0.2700	0.2474	0.2990	0.3096	0.2793	0.2171	0.2482	0.2374

a) and standard deviation (Table 7b).

We initially considered two levels of random effects, namely sample plot and repeated measurement. Our analyses indicated that the impact of repeated-measurement-related random effects were nonsignificant. They only increased computational complexity and caused problems with model convergence, along with producing inconsistent parameter estimates.

The sample plot-level random effects were higher than the residual standard deviation for the models of the mean for all species (Table 7a). Similarly, the random effects were higher for the model of standard deviation for all species, except for paper birch and red maple. The comparative numbers of random effects and residuals standard deviation could infer that the random effects were substantial in explaining the model variability.

We examined models’ performances based on cross-validated statistics, namely MB, MAB, and R². The models yielded unbiased prediction at 90% confidence intervals for the mean (

Table 8. Species-wise regression model performance statistics, using regression for the (a) mean and (b) standard deviation. The species are balsam fir (BF), eastern white pine (EWP), paper birch (PB), red maple (RM), red pine (RP), sugar maple (SM), trembling aspen (TA), and white spruce (WS).

Species	Mean Bias			Mean Absolute Bias			R2
	Median	5%	95%	Median	5%	95%	Median	5%	95%
(a) Model performance for DBH mean (cm)
BF	0.046	−0.488	0.439	0.481	0.420	0.564	0.906	0.828	0.947
EWP	−0.213	−5.204	3.773	2.742	2.005	4.354	0.844	0.000	0.927
PB	−0.046	−1.081	0.884	1.034	0.821	1.389	0.786	0.564	0.933
RM	−0.427	−4.388	1.143	0.885	0.567	4.453	0.758	0.000	0.947
RP	1.134	−7.437	5.492	1.807	1.276	3.869	0.915	0.631	0.968
SM	−0.071	−2.526	1.669	1.394	1.039	1.837	0.852	0.722	0.923
TA	−0.183	−1.915	1.452	1.794	1.422	2.311	0.776	0.578	0.913
WS	0.008	−1.978	1.540	1.236	0.987	1.620	0.909	0.823	0.950
(b) Model performance for DBH standard deviation (cm)
BF	−0.044	−0.560	0.310	0.871	0.722	1.098	0.349	0.000	0.680
EWP	0.803	−0.424	1.785	3.092	2.555	3.662	0.281	0.000	0.500
PB	−0.233	−0.862	0.315	1.338	1.096	1.745	0.182	0.000	0.621
RM	−0.222	−7.427	0.923	1.298	0.845	8.229	0.188	0.000	0.652
RP	0.495	−0.441	1.426	2.724	2.097	3.652	0.000	0.000	0.351
SM	0.126	−1.977	1.482	2.353	1.821	3.192	0.408	0.000	0.687
TA	−0.128	−0.648	0.372	1.503	1.254	1.769	0.461	0.249	0.599
WS	0.106	−1.016	0.966	1.798	1.471	2.294	0.388	0.000	0.673

a) and the standard deviation (Table 8b) for all species. The median of mean absolute bias varied between 0.48 and 2.74, and between 0.72 and 2.56 for mean and standard deviation, respectively. Likewise, the R² varied between 0.76 and 0.92 for the mean and 0 and 0.46 for the standard deviation.

Corresponding to unbiased parameters (Table 8), our simulated DBH distributions applying the recovered parameters (using the parameter recovery technique) from the model-predicted mean and standard deviation (Equation (11)) also produced similar patterns of DBH class distribution, as presented in Figure 4a,b.

4. Discussion

DBH is a uniquely important attribute in forest research and management planning. It can be measured or estimated at the individual tree level or at a stand level. For management and harvest schedule planning, stand-level attributes are needed. Although the measurement of individual tree DBHs is common, stand-level estimates may often be more precise than individual tree level measurements, due to higher chances of error accumulation [57]. On the other hand, stand-level DBH distribution can be used to reconstruct individual tree DBHs [30] or size class distributions [58].

Fitting DBH probability distributions and undertaking the reconstruction of distributions are not new techniques, but they are rarely used in practice and only for dominating species in localized geography [24,26]. Nonetheless, these still lack models for multiple species growing in mixed forests. We evaluated three commonly tested PDFs for eight commercial tree species in Ontario grown in natural mixed stands. The use of Weibull PDF for DBH distribution fitting is very common in a wide range of the literature [59], with a few studies reflecting comparative analyses among competing PDFs for different species globally, e.g., [60,61].

We compared three PDFs: Weibull, gamma, and log-normal. Log-normal is the simplest form of distribution, after negative exponential, to fit positive continuous random variables like DBH. Our preliminary analysis indicated that the one-parametric negative exponential did not fit at all for any species. But all three PDFs generally fit well for all species except red pine. Among these PDFs, log-normal fitted the best at the landscape and sample plot levels (Figure 2, Table 5) for all species. In a similar study, Rijal and Sharma [30] found that gamma PDF was the best among these three PDFs for two dominant species, jack pine and black spruce, which are also grown in Ontario. This difference could be because of the number of individual trees in the sample plots and the range of DBH. The less dominating eight species consisted of a lower number of trees within the sample plots, but the range of DBH was wider. This resulted in higher standard deviations and, hence, a longer right tail (Figure 3 and Figure 4) for these species, compared to the ones for jack pine and black spruce. Highly right-skewed data favor log-normal distribution, which is consistent with the previous works, e.g., [60,61].

We confined this study to the naturally grown forest stands. The DBH distribution of species in naturally grown forests features opposite J-shaped curves [62]. This characteristic reverse J-shape is approximately met based on the graphical presentation for all species, with a very high negative slope at the beginning of the curves (Figure 3). The DBH class distributions show a similar pattern for all species, except for trembling aspen and white spruce in the DBH class (0–5] and the DBH class up to (25–30] for red pine (Figure 4). A reason for such a lower density in the DBH class below five cm might have been attributed to the undercounting of seedlings and saplings. To our knowledge, there were no recent disturbances in the area. Nevertheless, the selected PDF could fit well for trembling aspen and white spruce. In the case of red pine, it seemed that there could have been more intensive forest management intervention in the recent decades that resulted roughly uniform distribution or may require multimodal PDF. The selected eight species were sub-dominant species after jack pine and black spruce in the study area. There could be a small number of individual trees that could escape from suppression, resulting in a DBH with a long right tail, which supports log-normal distribution.

The first part of this study examined two approaches to parameter estimations, namely maximum likelihood and parameter recovery methods for the selected three PDFs, but both require solving non-linear functions. The parameter recovery technique involves stand-level attributes [37], belonging to simple statistics, such as the mean, the standard deviation, a percentile of an array of the data, or the dominant diameter [22,26], among others. Among them, we used the mean and standard deviation, which are the simplest statistics and most common variables in forest management planning [3]. Using these variables in determining distribution functions worked well compared to the maximum likelihood method (Table 3 and Table 5, Figure 4). Mathematically, the maximum likelihood or moment-based methods (Equations (1)–(9)) are the same, and we should have the same results, provided there are no computational problems. Our results with the use of the selected three PDFs have demonstrated that DBH distribution can be simulated using any parametric PDFs for any range of DBHs for any species and can evaluate their performance from such simple statistics (mean and standard deviation).

From the perspective of forest management planning, our study would help the development of DBH distributions in response to various silvicultural treatments with the use of dynamic variables. The statistics of DBH change over time in response to various natural and anthropogenic disturbance factors. Consequently, either only the parameters of the PDF or the form of distribution may vary though time. More importantly, when there is intensive forest management, the growth trajectories of different species can vary substantially in a way that can be addressed with the use of simple statistics. Fitting regression models with dendrometric, geographic, and climatic variables in the second component of our work helped to update the parameters of PDF. Among the evaluated model covariates, the dendrometric variables could explain more than the selected geographic and climatic variables for most of the examined species. The parameters associated with geographic variables (longitude and latitude), as well as the annual mean and variability (standard deviation) of bioclimatic variables, were also significant for some species like balsam fir and eastern white pine when we fit of an array of DBH for eight species for mean and standard deviation (Table 7a-b). These results correspond to previous pieces of work that showed that geography and climate can substantially impact DBH growth [20,49]. When dendrometric variables have a strong relationship with the response variables, like DBH, the bioclimatic variables may come to be insignificant [32], even though these variables are always important in forest growth and development.

Stand dendrometric attributes like basal area, stem density, QMD, dominant height, and age are common as covariates in predictive modelling, e.g., [22,26,30,31]. Based on the data availability, we limited dendrometric covariates (Table 1), and the regression evaluations demonstrated that total and species-wise stand attributes, especially QMD, stem density, and basal area, were significant covariates for predicting the mean and standard deviation for the selected species, but showed different magnitudes and signs of parameters (Table 5). For example, total basal area as a model covariate had a negative parameter sign whenever it was significant for all species in regressions for DBH mean and standard deviation, except for the mean of trembling aspen. Stem density had a positive relationship with the DBH mean. Likewise, the QMD was positively related to the mean and the standard deviation. A reason for the opposite behaviors of the stem density and basal area could be due to the fact that the stands were dominated by small-sized trees [63]. Such a parametric relationship with the response variable, DBH mean, was consistent with an earlier study for two dominant species: black spruce and jack pine, in the same study area [32].

Our regression models for the DBH mean and standard deviation for eight species using various model covariates performed differently in terms of model performance statistics. In all cases, they yielded unbiased predictions. The magnitude of cross-validated other statistics, namely MAB and R², varied by model for each species, but the median and 90 percentile range (confidence interval) values inferred that these models’ fit had precise predictions. However, there was a wide range of variations in R², especially for the model of standard deviation. The model for the mean had a range of R² between 0.76 and 0.91, whereas the values for the standard deviation were between zero and 0.46. The mean values are generally related linearly to the examined dendrometric covariates, resulting in a high R². On the other hand, the standard deviation, a form of transformation of the DBH, is not linearly related to the model covariates, resulting in low R² values despite seeing unbiased predictions.

The simulation of DBH distribution for multiple species is challenging for sub-dominating species because of their limited presence in various proportions or absence in many sample plots. There were more than 80 individual tree species listed. Modelling for all individual trees’ DBHs is neither desirable for increased complexity, nor do we achieve better estimates for very small number of sample trees. For better consistency of model parameters, we set a minimum threshold of five individual trees, and, even with this, we often encountered problems with convergence or inconsistent parameters and required several iterations for better results. Despite the possibility that other forms of PDF may still be better than ours for some species we covered, we had to limit ourselves to only these three PDFs. Nevertheless, there can be some extended avenue of improvement, leading to covering a wide spectrum by developing multi-modal distributions to capture the characteristics of the growth responses of different types of silvicultural operations in managed forests. For example, distribution fitting for red pine, in our case, may have warranted exploring multimodal PDFs, but we limited ourselves to unimodal PDFs for some resource constraints and because we hoped to keep this method operationally simple and feasible for other practitioners. A general model can be developed for minor species, in terms of their presence, by taking them together. A comparative study of the model parameters of individual trees and combined parameters for all minor species may help us to understand the species’ ecology and growth dynamics in the given forest stands. These parameters, more specific to general levels of PDF and regression fit, would help with selecting better silviculture systems for directing and maintaining forest growth and development trajectory, as desired by forest management planners and managers.

5. Conclusions

This study comprised a comparative analysis of three probability distribution functions for the individual tree diameters at breast height (DBHs) of eight sub-dominating tree species. Past studies have shown that the Weibull probability distribution function (PDF) is the best fit for dominant species, whereas our results indicated that log-normal PDF was the best fit for sub-dominant species. While several PDFs have been evaluated in research, only a few have been used in practice. The parameter estimation of individual PDFs using the mean and standard deviation can be simpler than the maximum likelihood method. Moreover, we illustrated that these simple statistics, mean and standard deviation, can be used simultaneously to estimate the parameters of three PDFs and to reconstruct DBH distributions using these PDFs. As all statistical PDFs are mathematically related to the mean and standard deviation, albeit with varying computation complexity, our study has opened an avenue of reconstruction and evaluation of any PDF using just the mean and standard deviation. Another advantage of the use of the mean is that it is usually modelled for many purposes. Our construction of predictive models using dynamic variables shows that we can fit dynamic probability distribution functions in response to forest management and or any type of forest growth or disturbances over time.