Performance of Kernel Estimator and Johnson S_B Function for Modeling Diameter Distribution of Black Alder (Alnus glutinosa (L.) Gaertn.) Stands

Abstract

We compare the usefulness of nonparametric and parametric methods of diameter distribution modeling. The nonparametric method was represented by the new tool—kernel estimator of cumulative distribution function with bandwidths of 1 cm (KE1), 2 cm (KE2), and bandwidth obtained automatically (KEA). Johnson S_B (JSB) function was used for the parametric method. The data set consisted of 7867 measurements made at breast height in 360 sample plots established in 36 managed black alder (Alnus glutinosa (L.) Gaertn.) stands located in southeastern Poland. The model performance was assessed using leave-one-plot-out cross-validation and goodness-of-fit measures: mean error, root mean squared error, Kolmogorov–Smirnov, and Anderson–Darling statistics. The model based on KE1 revealed a good fit to diameters forming training sets. A poor fit was observed for KEA. Frequency of diameters forming test sets were properly fitted by KEA and poorly by KE1. KEA develops more general models that can be used for the approximation of independent data sets. Models based on KE1 adequately fit local irregularities in diameter frequency, which may be considered as an advantageous in some situations and as a drawback in other conditions due to the risk of model overfitting. The application of the JSB function to training sets resulted in the worst fit among the developed models. The performance of the parametric method used to test sets varied depending on the criterion used. Similar to KEA, the JSB function gives more general models that emphasize the rough shape of the approximated distribution. Site type and stand age do not affect the fit of nonparametric models. The JSB function show slightly better fit in older stands. The differences between the average values of Kolmogorov–Smirnov (KS), Anderson–Darling (AD), and root mean squared error (RMSE) statistics calculated for models developed with test sets were statistically nonsignificant, which indicates the similar usefulness of the investigated methods for modeling diameter distribution.

1. Introduction

Diameter at breast height (dbh) is a basic tree attribute that is important to both forestry practice and science. It is a simple to measure, and closely related to other difficult to obtain tree features. Diameter as an explanatory variable is used to predict tree attributes such as height [1,2], crown length [3,4], bark thickness [5], volume [6,7], biomass [8,9], and carbon storage [10].

The mean diameter provides only a general view of a tree stand, and it is required to know dbh distribution in order to obtain detailed information about the stand [11]. Research focused on mathematical description of dbh distribution has long history. Initially, diameter frequency models were mainly used as part of yield models. Over time, it was obvious that dbh distribution models can essentially support other areas of forestry practice and science. Examples include the optimization of silviculture treatments, preharvest prediction of the kind and amount of wood assortments, tree bucking control, forest resource evaluation, preparation of timber trade contracts, and financial planning undertaken by forest holdings. A close relationship between diameter distribution and stand volume from one side and between stand volume and financial income from the second side justify the inclusion of diameter distribution models as part of decision support systems [12,13,14,15]. In addition to the abovementioned strictly quantitative approach, mathematically described diameter distribution can provide valuable information about forest dynamics [16] and support the management of uneven-aged and mixed forests [17,18].

Forests as a multipurpose site provide wood and other goods and services that are important to society. Forest managers requires tools for evaluating different management practices and their effects on stand structure [19]. In this context, dbh distribution models can provide very valuable information that constitute the first link of the management chain—the so-called three Ds: data, dynamics, and decisions [20].

Dbh distribution modeling is an example of pattern recognition [21]. It combines a set of different issues related first to the choice of method that is adequate to approximate the number of trees belonging to distinguished diameter classes. Other issues related to dbh distribution modeling include inventory of forest resources [22], geostatistics [23,24], and remote data acquisition [25,26,27,28] combined with the application of advanced machine learning methods [10,29].

One of the most popular approaches that allows a mathematical description of dbh distribution relies on the use of available theoretical probability density functions and is commonly referred to as “parametric”. As diameters are thought to be realizations of some kind of continuous distribution [30], the research scope is to find adequate function. To date, there are no biological reasons to justify the choice of any particular function [31,32]. The basic premise justifying the use of theoretical distribution is usually the ease of estimating and interpreting its parameters as well as the goodness-of-fit to empirical data. Examples of theoretical functions applied to date to model frequency of diameter include the following distributions: normal [33,34], double normal [35], Weibull [36], Johnson S_B [37], beta and gamma [37]. The usefulness of the aforementioned models for describing dbh distribution of black alder stands was assessed in the authors’ previous paper [38]. The best results were obtained using Johnson S_B function, while the worst result was obtained for normal distribution. In the frame of parametric approach are also used finite mixture models. They combine two or more single functions (e.g., Weibull or gamma) and are mainly used to model irregularly shaped diameter frequencies occurring in mixed-species and uneven-aged stands [39,40,41,42].

In addition to the parametric approach, there are a number of alternative methods for modeling dbh distribution that do not rely on any predefined function and are jointly termed as “nonparametric”. The concept of nonparametric methods is based on the optimal use of information contained in the collected data [43]. Previous studies have reported the following examples of nonparametric approaches used in the context of dbh distribution modeling: percentile method [44,45,46], k-nearest neighbor method [25,47,48,49,50], decimal curves [51], neural networks [52,53], and kernel estimators [54,55,56,57,58,59].

The application of kernel estimator to model dbh frequency does not require any assumptions on the form of the underlying density, which is a characteristic of nonparametric methods [60]. On the other hand, this method requires the selection of the type of kernel function and adequate value of bandwidth, which is also known as a smoothing parameter or window. Kernel function ensures continuity and is responsible for weighting observation [61]. Examples of the available kernel types include Gaussian, Epanechnikov, uniform, biweight, and triweight. The type of kernel has minor influence on the performance of estimator. However, the properties of kernel function, including, for example, the order of its derivatives, play an important role in the optimization process used to select smoothing parameter. Bandwidth determines the neighborhood of the point (diameter record) used during the estimation of the probability value [62]. Taubert et al. [63] focused on the importance of the chosen diameter class width to the final effects of parametric diameter distribution modeling. Similar to kernel estimator, the selected bandwidth has major influence on the results of diameter frequency estimation. Kernel estimator with too small bandwidth provides a dbh distribution model with excessively highlighted local extrema. On the other hand, an estimator with high values of bandwidth causes a very flat final distribution shape and neglects local extrema. This shows the crucial role of bandwidth in determining a bias–variance tradeoff [64].

Because of few studies on the use of kernel estimator to describe dbh distribution, it is actually a poorly recognized tool in forest resource modeling. Results of previous research by Pogoda et al. [65] provide only preliminary information about the size and variability of the smoothing parameter depending on the age of the alder stand, the size of used sample and the chosen method of bandwidth selection. However, a detailed analysis focused on the assessment of the direct use of kernel estimator to model dbh distribution is lacking.

The present study aimed to assess the effects of diameter distribution description for black alder stands by using kernel estimator with different bandwidth values and then performing a comparative analysis of the usefulness of kernel estimator and Johnson S_B function for modeling frequency of diameters.

2. Materials and Methods

2.1. Data Collection

Data from 36 managed black alder stands located in southeastern Poland in the Dąbrowa Tarnowska (10 stands) and Niepołomice (26 stands) Forest District were used in this study. The stands represent four age groups, i.e., 20, 40, 60, and 80 years old and three distinguished site types: moist mixed broadleaved forest (MMB), moist broadleaved forest (MB), and alder swamp forest (AS). Each age group was represented by three stands, while site type by 12 stands (4 age groups of 3 stands). All stands were characterized by the share of black alder not less than 0.9 and the ratio of actual to tabular growing stock of above 0.8. The basic attributes characterizing the studied stands are summarized in

Table 1. Basic attributes of the studied black alder stands.

Attribute	Minimum	First Quartile	Median	Third Quartile	Maximum	Mean	Standard Deviation
Number of measured trees	138	181	206	260	359	219	51
Density (stems ha−1)	348	424	570	914	1505	725	377
Reineke stand density index (stems ha−1)	249	414	506	584	692	504	103
Quadratic mean diameter (cm)	11.1	15.7	24.1	30.0	35.0	23.1	8.2
Basal area (m2 ha−1)	9.83	17.57	25.17	30.95	39.49	24.78	8.01

.

In each stand, 10 circular sample plots were established that were located systematically in the nodes of a grid of squares. The side of a square depended on the stand area and ranged from 20 to 70 m. The size of sample plots was 0.02, 0.03, 0.04, and 0.05 ha for stands in age 20, 40, 60, and 80 years old, respectively. For each tree in the sample plot, the dbh (1.3 m above ground) was measured using a caliper in two perpendicular directions with the accuracy of 0.1 cm. The average value from the two measurements was taken for further analysis. In total, 7867 black alders were measured, whose diameters ranged from 6.7 to 61.7 cm.

2.2. Diameter Distribution Models

In the first stage of analysis, the cumulative frequency of diameters in the considered stands was approximated using kernel estimator of cumulative distribution function (CDF), proposed in 1964 by Nadaraya, which is expressed by the following formula [66]:

$${\hat{F}}_h\left(x\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}K\left(\frac{x-X_i}{h}\right)},$$

(1)

where n is sample size, x is the value of dbh (cm) for which the CDF was estimated, $X_i$ is individual dbh belonging to the sample, h is a bandwidth (h > 0, cm), $K\left(·\right)$ is an integral of kernel function k expressed by the following formula:

$$K\left(x\right)={\displaystyle \underset{-\infty }{\overset{x}{\int }}k\left(t\right)dt},$$

(2)

In the presented analysis, the estimator with Gaussian kernel was applied as follows:

$$k\left(x\right)=\frac{1}{\sqrt{2\mathsf{\pi }}}\exp \left(-\frac{x^2}{2}\right),$$

(3)

Cumulative frequency of diameters was approximated by kernel estimator of CDF with three different smoothing parameter values. Kernel estimator denoted as “KEA” use bandwidth obtained automatically in the process of optimizing the discrepancy between empirical and predicted probability. Optimization is carried out with the selected criterion, which may include mean integrated squared error (MISE) or asymptotic mean integrated squared error (AMISE) [67]. The research problem is to find an adequate method to estimate the value of predicted probability. In the present study, a plug-in method proposed by Altman and Legér [68] was used for this purpose. Pilot bandwidth required by the above method during the estimation of predicted probability value was selected by Silverman’s rule of thumb [67,69]. Kernel estimators denoted as “KE1” and “KE2” use bandwidth 1 cm and 2 cm, respectively. The mentioned values correspond to the span of diameter classes commonly used in practical presentation of building stand tables.

In the second stage of analysis, frequency of diameters was described by Johnson S_B (“JSB”) function, which was applied for the first time in dbh distribution modeling context by Hafley and Schreuder [37]. The JSB probability density function is given by the formula:

$$f\left(x\right)=\frac{\mathsf{\delta }}{\sqrt{2\mathsf{\pi }}}\cdot \frac{\mathsf{\lambda }}{\left(x-\mathsf{\xi }\right)\left(\mathsf{\xi }+\mathsf{\lambda }-x\right)}\exp \left\{-\frac{1}{2}{\left[\mathsf{\gamma }+\mathsf{\delta }\ln \left(\frac{x-\mathsf{\xi }}{\mathsf{\xi }+\mathsf{\lambda }-x}\right)\right]}^2\right\},$$

(4)

and the CDF takes the form as:

$$F\left(x\right)=\mathrm{\Phi }\left(\mathsf{\gamma }+\mathsf{\delta }\ln \left(\frac{z}{1-z}\right)\right),z\equiv \frac{x-\mathsf{\xi }}{\mathsf{\lambda }},$$

(5)

where x is dbh (cm), γ (–∞ < γ < +∞) and δ (δ > 0) are shape parameters denoting asymmetry and kurtosis, respectively, λ (λ > 0) is a scale parameter, ξ (–∞ < ξ < +∞) is location parameter, and Φ is Laplace integral.

The parameters of JSB function were estimated from detailed diameter records using the maximum likelihood method (MLE).

2.3. Model Evaluation and Validation

The performance of KE1, KE2, KEA, and JSB function was assessed by the cross-validation method [70,71,72]. It allows to avoid some too optimistic results (“optimism principle”) [73], which can occur when model assessment is performed using the same data that was used to construct the model. Considering the hierarchical structure of the collected data (levels: single tree and sample plot), it was assumed that the sample plot will be the basic unit to group the collected data for resampling. This corresponds to the fact that in forestry research, the data representing sample plot are most often used for model development. It refers both to the traditional process of model development based on field inventories and to the construction of models using remotely sensed data in the frame of area-based approach (ABA) [27]. In the present study, cross-validation was performed using the leave-one-plot-out cross-validation method (LOPOCV) [74]. From the set of 10 sample plots established in stand, 9 sample plots were used as a training set, and the remaining one plot constituted the test set. Detailed diameter records forming the training set were used to estimate the parameters of JSB function and to the obtain bandwidth for KEA. The LOPOCV procedure was repeated 10 times, so that each sample plot established in the stand was used exactly once as a test set. Goodness-of-fit was assessed using the following statistics:

-

mean error (ME):

$$ME=\frac{{\displaystyle {\sum }_{i=1}^n\left(F\left(x_i\right)-\widehat{F}\left(x_i\right)\right)}}{n},$$

(6)

-

root mean squared error (RMSE):

$$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^n{\left(F\left(x_i\right)-\widehat{F}\left(x_i\right)\right)}^2}}{n}},$$

(7)

-

the Kolmogorov–Smirnov statistic (KS) [73]:

$$KS=\underset{x}{\sup }\left(F\left(x_i\right)-\widehat{F}\left(x_i\right)\right),$$

(8)

-

the Anderson–Darling statistic (AD) [75]:

$$AD=-n-{\displaystyle {\sum }_{j=1}^n\left(2j-1\right)}\left[\ln \left({\widehat{F}}_j+\ln \left(1-{\widehat{F}}_{n-j-1}\right)\right)\right]/n,$$

(9)

where $F\left(x_i\right)$ is the value of empirical CDF for tree with dbh $x_i$, $\hat{F}\left(x_i\right)$ is the value of CDF predicted using the JSB function or kernel estimator for tree with dbh $x_i$, n is the number of measured trees, j are the ranks assigned to dbh’s set in ascending order, and ${\hat{F}}_j$ and ${\hat{F}}_{n-j-1}$ are the CDF values predicted by JSB function or kernel estimator to dbh ranked as j and n − j − 1, respectively.

Final comparison of kernel estimator and JSB function was based on the analysis of mean values and variability of statistics obtained with Equations (6)–(9) for training and test sets. The analyses were performed for all 36 stands and within the distinguished age groups and site types. By comparing calculated KS and AD statistics with their critical values, it was assessed in how many cases the CDF approximated by kernel estimator or JSB function differed significantly from the empirical CDF at the assumed level of significance α = 0.05. Because of the lack of normality, the Kruskal–Wallis test was used to validate the significance of differences between the average values of goodness-of-fit statistics (Equations (6)–(9)) for assessing the results of dbh distribution modeling by the considered methods.

Calculations and analyses were performed using the EasyFit 5.5 Professional Program and packages sROC [76] and SuppDists [77] available in R environment [78].

3. Results

By using kernel estimator with three bandwidth variants, namely KE1, KE2, and KEA, the nonparametric models of dbh distribution for all (360) training sets were obtained. The KS test showed significant (α = 0.05) differences between empirical and estimated dbh frequencies for 34 training sets for the KE2 model. A similar analysis (KS test, α = 0.05) performed for test sets showed significant differences for 100, 94, and 95 sets for KE1, KE2, and KEA models, respectively. The outcomes of the AD test showed significant (α = 0.05) differences for 68 training sets for only the KE2 model and for 134, 131, and 124 test sets for KE1, KE2, and KEA, respectively. The MLE used for parameter estimation of JSB function failed for 48 (13.3%) training sets due to lack of convergence. Consequently, the parametric model was used for only 312 training and test sets. From the results of the KS test, significant differences were found between empirical dbh distribution and estimated by JSB function for 25 training and 82 test sets, while the AD test showed significant differences for 25 training and 118 test sets.

For training sets, the best fit for nonparametric dbh distribution models (KS = 0.0293, AD = 0.2594, and RMSE = 0.0121) was observed for the KE1 model and the worst fit (KS = 0.0569, AD = 1.1252, and RMSE = 0.0268) for the KEA model (

Table 2. Mean and standard deviation (in brackets) of measures of goodness-of-fit; the best values are in bold.

Method	ME	RMSE	KS	AD	ME	RMSE	KS	AD
	Training Sets				Test Sets
KEA	0.0059 (0.001)	0.0268 (0.004)	0.0569 (0.009)	1.1252 (0.138)	0.0128 (0.016)	0.1426 (0.046)	0.2440 (0.068)	3.0245 (1.745)
KE1	0.0059 (0.001)	0.0121 (0.004)	0.0293 (0.008)	0.2594 (0.303)	0.0121 (0.016)	0.1450 (0.048)	0.2458 (0.069)	3.4159 (2.071)
KE2	0.0059 (0.001)	0.0231 (0.011)	0.0514 (0.022)	1.1853 (1.480)	0.0126 (0.016)	0.1441 (0.047)	0.2446 (0.068)	3.2046 (1.918)
JSB	0.0039 (0.046)	0.0326 (0.059)	0.0666 (0.103)	4.2440 (16.432)	0.0079 (0.030)	0.1431 (0.049)	0.2427 (0.069)	3.4290 (2.128)

KE1, KE2, and KEA: kernel estimator with bandwidth 1 cm, 2 cm, and obtained automatically, respectively; JSB: Johnson SB function; ME: mean error; RMSE: root mean squared error; KS: Kolmogorov–Smirnov statistic; and AD: Anderson–Darling statistic.

). In terms of fitting to training sets, the parametric model developed on the basis of JSB function was the worst (KS = 0.0666, AD = 4.2440, and RMSE = 0.0326), in which the largest variability of KS, AD, and RMSE was also observed (Table 2). The differences between the average values of KS, AD, and RMSE observed for parametric and nonparametric models for training sets were statistically significant (Kruskal–Wallis test, α = 0.05).

The goodness-of-fit assessment carried out for test sets showed that for the nonparametric method, the smallest differences between the empirical and modeled CDF (KS = 0.2440, AD = 3.0245, and RMSE = 0.1426) were observed for KEA and the largest for KE1 (Table 2). The assessment of the effectiveness of the parametric JSB model strongly depended on the fitting criterion used. In terms of KS statistics (KS = 0.2427), it had the best fit, while according to AD statistics (AD = 3.4290), it had the worst fit. The value of RMSE statistics (0.1431) placed JSB function in the second spot immediately below the KEA. The differences between the average values of the goodness-of-fit statistics observed for test sets approximated by KE1, KE2, KEA, and JSB function were statistically nonsignificant (Kruskal–Wallis test, α = 0.05).

An example of diameter distribution modeling by the nonparametric method represented by KE1, KE2, KEA, and the parametric method based on JSB function is shown in Figure 1. The frequency of diameters belonging to the training set (Figure 1a), which was used for model development, was very well fitted by KE1. The application of KEA led to a much more general model that does not consider local maxima and minima in diameter frequency, which were observed in sample distribution. The parametric model based on JSB function is also a more general model. The comparison of estimated models for an independent data—the test set (Figure 1b)—showed the advantage of general models reflecting the rough shape of distribution over models that excessively fit the shape of the empirical distribution from the sample used for estimation.

The analysis of the average values of goodness-of-fit measures (KS, AD, and RMSE;

Table 3. Mean and standard deviation (in brackets) for KS statistic evaluated goodness-of-fit in stands of distinguished site and age groups; the best and worst values are in bold and italics, respectively.

Method	KEA	KE1	KE2	JSB	KEA	KE1	KE2	JSB
Group of Stands	Training Sets				Test Sets
MMB	0.0585 (0.008)	0.0304 (0.009)	0.0523 (0.027)	0.0517 (0.038)	0.2607 (0.089)	0.2625 (0.092)	0.2617 (0.090)	0.2529 (0.095)
MB	0.0544 (0.010)	0.0286 (0.007)	0.0552 (0.022)	0.0876 (0.106)	0.2204 (0.051)	0.2214 (0.051)	0.2209 (0.050)	0.2213 (0.052)
AS	0.0577 (0.009)	0.0289 (0.009)	0.0466 (0.018)	0.0841 (0.124)	0.2510 (0.056)	0.2534 (0.057)	0.2511 (0.056)	0.2508 (0.055)
20	0.0561 (0.011)	0.0405 (0.009)	0.0847 (0.019)	0.0393 (0.011)	0.2000 (0.041)	0.2008 (0.041)	0.2011 (0.039)	0.2041 (0.041)
40	0.0525 (0.008)	0.0274 (0.003)	0.0459 (0.006)	0.0963 (0.131)	0.2108 (0.047)	0.2120 (0.048)	0.2105 (0.047)	0.2123 (0.052)
60	0.0596 (0.009)	0.0255 (0.003)	0.0384 (0.006)	0.1303 (0.159)	0.2949 (0.078)	0.2977 (0.080)	0.2958 (0.079)	0.2915 (0.089)
80	0.0592 (0.006)	0.0239 (0.002)	0.0364 (0.004)	0.0422 (0.003)	0.2705 (0.055)	0.2726 (0.057)	0.2709 (0.057)	0.2587 (0.059)

MMB: moist mixed broadleaved forest; MB: moist broadleaved forest; AS: alder swamp forest; 20, 40, 60, and 80: stands in age of 20, 40, 60, and 80 years old, respectively; other abbreviations as denoted in Table 2.

,

Table 4. Mean and standard deviation (in brackets) for AD statistic evaluated goodness-of-fit in stands of distinguished site and age groups; the best and worst values are in bold and italics, respectively.

Method	KEA	KE1	KE2	JSB	KEA	KE1	KE2	JSB
Group of Stands	Training Sets				Test Sets
MMB	1.1730 (0.196)	0.2875 (0.362)	1.1527 (1.257)	1.0390 (2.355)	3.3959 (2.313)	3.9156 (2.782)	3.6743 (2.605)	3.5222 (2.858)
MB	1.1247 (0.071)	0.3409 (0.358)	1.6664 (1.895)	7.9067 (19.147)	2.5227 (1.004)	2.7436 (1.142)	2.5922 (1.049)	2.8625 (1.292)
AS	1.0780 (0.110)	0.1498 (0.102)	0.7368 (0.760)	7.5225 (21.186)	3.1548 (1.700)	3.5886 (1.951)	3.3474 (1.773)	3.8371 (2.091)
20	1.2658 (0.124)	0.6775 (0.355)	3.2760 (1.598)	0.4296 (0.288)	2.1975 (0.936)	2.3119 (1.099)	2.2198 (0.935)	2.4966 (1.152)
40	1.0813 (0.096)	0.1690 (0.080)	0.8240 (0.500)	11.1023 (23.287)	2.5412 (1.895)	2.8149 (2.130)	2.5891 (1.913)	3.1441 (2.539)
60	1.1104 (0.143)	0.1074 (0.026)	0.3755 (0.146)	11.7151 (24.138)	3.9786 (2.016)	4.5411 (2.317)	4.2709 (2.148)	4.4073 (2.423)
80	1.0434 (0.075)	0.0836 (0.022)	0.2657 (0.102)	0.2946 (0.071)	3.3805 (1.611)	3.9957 (2.018)	3.7389 (1.954)	3.6024 (2.131)

Abbreviations as denoted in Table 3.

and

Table 5. Mean and standard deviation (in brackets) for root mean squared error (RMSE) evaluated goodness-of-fit in stands of distinguished site and age groups; the best and worst values are in bold and italics, respectively.

Method	KEA	KE1	KE2	JSB	KEA	KE1	KE2	JSB
Group of Stands	Training Sets				Test Sets
MMB	0.0271 (0.004)	0.0120 (0.004)	0.0223 (0.011)	0.0242 (0.024)	0.1525 (0.062)	0.1554 (0.065)	0.1547 (0.063)	0.1478 (0.068)
MB	0.0253 (0.004)	0.0129 (0.004)	0.0257 (0.013)	0.0449 (0.061)	0.1289 (0.035)	0.1302 (0.036)	0.1297 (0.035)	0.1300 (0.037)
AS	0.0281 (0.003)	0.0113 (0.003)	0.0213 (0.010)	0.0419 (0.080)	0.1466 (0.038)	0.1496 (0.039)	0.1480 (0.038)	0.1482 (0.038)
20	0.0253 (0.005)	0.0175 (0.003)	0.0398 (0.009)	0.0173 (0.006)	0.1136 (0.028)	0.1142 (0.028)	0.1152 (0.026)	0.1163 (0.029)
40	0.0259 (0.004)	0.0114 (0.001)	0.0214 (0.003)	0.0502 (0.075)	0.1191 (0.034)	0.1203 (0.034)	0.1192 (0.034)	0.1212 (0.038)
60	0.0290 (0.003)	0.0103 (0.001)	0.0169 (0.003)	0.0705 (0.089)	0.1784 (0.052)	0.1825 (0.055)	0.1804 (0.0504)	0.1781 (0.060)
80	0.0272 (0.002)	0.0091 (0.001)	0.0143 (0.002)	0.0161 (0.002)	0.1595 (0.038)	0.1632 (0.040)	0.1618 (0.040)	0.1528 (0.044)

Abbreviations as denoted in Table 3.

) within the distinguished site and age groups confirmed the previously observed relationships. Among the nonparametric models: KE1, KE2, and KEA, in the majority of the distinguished group of stands, the best fit to empirical diameter frequencies of training sets was obtained using KE1, whereas the worst fit was observed for KEA (Table 3, Table 4 and Table 5). Compared to the nonparametric method of kernel estimation, the parametric model based on JSB function showed in most cases a much worse fit to data forming training sets. The evaluation of dbh distribution models performed on independent samples of diameter forming test sets indicated that among the nonparametric models, the best fit in the majority of site and age groups of stands was found for KEA, and the worst was for KE1 (Table 3, Table 4 and Table 5).

The average value of bandwidth to KEA determined using samples of diameters representing stands growing in MMB, MB, and AS forest was 2.91, 2.33, and 2.95 cm, respectively, and for particular age groups, namely 20, 40, 60, and 80 years old, the value was 1.46, 2.41, 3.34, and 3.71 cm, respectively. Thus, in every case, the bandwidth was greater than 1 cm, and in most cases, it was more than 2 cm, which resulted in a much worse fit of KEA to the training sets compared to that of KE1 and KE2. At the same time, a more general shape of the distribution obtained by KEA allows better fit to independent data forming test sets. The evaluation of JSB function carried out on test sets for the distinguished groups of stands, especially for age groups, showed that this model was inferior to the nonparametric models in stands of age 20 and 40 years old. In older stands, i.e., 60 and 80 years old, the model based on JSB function proved to be better than KE1, KE2, and KEA (Table 3, Table 4 and Table 5).

The Kruskal–Wallis test (α = 0.05) conducted independently for separate site and age groups of stands showed significant differences between the average values of RMSE, KS, and AD statistics calculated for the analyzed models on the basis of data of training sets. For independent test sets, differences in the average values of the abovementioned statistics were nonsignificant, except for the KS statistic determined for MB forest (

Table 6. Summary of Kruskal–Wallis test; bold and italic values refer to training and test sets, respectively.

Group of Stands	Statistic	χ2	p
MMB	KS	22.5214/0.2970	0.0001/0.9606
	AD	20.9270/0.4461	0.0001/0.9306
	RMSE	28.8302/0.2932	0.0000/0.9613
	ME	0.2846/0.0900	0.9629/0.9930
MB	KS	22.4515/18.9881	0.0001/0.0003
	AD	15.7228/0.5850	0.0013/0.8998
	RMSE	19.5009/0.0264	0.0002/0.9989
	ME	26.4907/0.2372	0.0000/0.9714
AS	KS	23.3006/0.1037	0.0000/0.9914
	AD	26.2651/1.0981	0.0000/0.7775
	RMSE	26.8973/0.2249	0.0000/0.9735
	ME	2.4891/2.9956	0.4773/0.3923
20	KS	22.0410/0.1732	0.0001/0.9818
	AD	25.3904/0.7377	0.0000/0.8643
	RMSE	24.3173/0.3894	0.0000/0.9424
	ME	2.2152/0.1612	0.5290/0.9836
40	KS	19.9688/0.2954	0.0002/0.9609
	AD	18.3556/0.3430	0.0004/0.9518
	RMSE	19.8372/0.0541	0.0002/0.9967
	ME	1.1652/0.7631	0.7614/0.8583
60	KS	25.8482/0.0746	0.0000/0.9947
	AD	24.7631/0.7852	0.0000/0.8530
	RMSE	26.1388/0.2393	0.0000/0.9710
	ME	0.0704/3.0270	0.9951/0.3875
80	KS	31.5405/0.4695	0.0000/0.9255
	AD	29.6987/1.6026	0.0000/0.6588
	RMSE	30.4795/0.7417	0.0000/0.8633
	ME	0.0831/0.4294	0.9938/0.9341
All data	KS	66.9480/0.1140	<0.0001/0.9901
	AD	60.2600/0.9725	<0.0001/0.8079
	RMSE	69.6550/0.1521	<0.0001/0.9849
	ME	0.0429/1.0486	0.9977/0.7895

χ²: value of test statistic; p: probability value; and other abbreviations as denoted in Table 2 and Table 3.

). Models were also evaluated in terms of the value of systematic error (ME). For both training and test sets, a slightly larger ME error was found for KE1, KE2, and KEA models than that noted for the parametric model of JSB function (Table 2). The observed differences of ME statistics for the developed models were nonsignificant, except for training sets representing MB forest (Kruskal–Wallis test, Table 6).

4. Discussion

The present study is a continuation of a previous investigation focused on modeling dbh distribution for managed black alder stands located in southeastern Poland [38,46,65]. Diameter distribution as a basic attribute of stand structure [79,80] is modeled using different methods that can be divided into two distinct categories: parametric and nonparametric. The parametric method used in the present study is represented by the JSB function. It was first used to describe the frequency of diameters for pine stands growing on the east coast of the USA [37]. Because of its good performance, the JSB function was widely used to model dbh structure of different stands growing in various geographical zones [30,81,82,83,84,85]. Previous studies on the application of JSB function in Poland showed contradictory findings. Zasada [86] reported satisfactory performance of the JSB function when applied to modeling diameter frequency of fir stands. The same author found low usefulness of this function to model dbh distribution of birch stands [87]. Jagiełło et al. [88] confirmed the usefulness of the JSB function for beech stands, while Rymer-Dudzińska and Dudzińska [89,90] indicated high instability of the results obtained using this model. Interesting results were also reported by Siekierski [91], who, as the first researcher in Poland, paid attention to the problem associated with the estimation of parameters to JSB function.

Our study results confirmed that the performance of JSB function used to modeling dbh distribution for black alder stands was in line with the earlier works listed above. The JSB function shows the worst fit to data forming training sets. Additionally, problems occurred for parameter estimation for 48 training sets. Frequencies of diameters belonging to test sets were, however, better approximated by the JSB function. While analyzing the performance of the JSB function, it should be noted that it is a four-parameter function. In some situations, a larger set of parameters cause that model to be more flexible and therefore, more adequate to describe dbh distributions with different shapes. Simultaneously, the application of the abovementioned model may cause the problem of lack of convergence to iterative procedures used for parameter estimation and additionally hamper their biological interpretation [92]. It should also be noted that JSB function is a member of the Johnson family of distribution [93]. Functions belonging to Johnson’s family have different normalizing transformations and by this cover other distribution shapes. It is possible that the estimation of JSB function parameters failed in our study because the collected diameter records did not meet the conditions for the given transformation.

Some studies [32,50,94] have emphasized that the method applied during parameters estimation has substantial influence on the relative performance of the chosen probability density function. In this situation, the application of nonparametric methods of dbh distribution modeling and using information from the observed stand tables seems to be a good support to the parametric approach. The performance of kernel estimator used in the present study strictly depends on the selected bandwidth. KEA, by using bandwidth obtained automatically as a result of the optimization process, emphasizes the general shape of distribution by ignoring local irregularities observed in training set diameter frequency. Different results were achieved using KE1 that can model extrema (minima and maxima) occurring in empirical diameter frequency. The good performance of KE1 can be very helpful in modeling irregular diameter distributions, although the well fit carries the risk that the final model is over-fitted, as observed in the present study. KE1 properly describes the distribution of the training set (Figure 1a), but when transferred to independent data (Figure 1b), it presents much worst fit. KEA, through its general shape, approximate very well to the distribution of diameter records forming test set. The performance of KE2 can be considered as intermediate between KEA and KE1. For training sets characterized by irregular shape of diameter distribution, KE2 depicted local extrema less excessively than KE1.

The observed good performance of KEA justifies the use of bandwidths obtained automatically on the basis of data used for model development instead of arbitrarily adopted values. In the process of bandwidth selection, it is worth additionally considering the hint given by Sheather [69] and to analyze simultaneously few values of bandwidth.

The practical application of kernel estimator requires some diameter records; thus, this approach allows only the interpolation of probability values without the possibility of their extrapolation [58]. By considering two aspects of dbh distribution modeling, first: prediction of the number of trees belonging to distinguished diameter classes solely on the basis of assessed stand attributes, and second: the shape description of dbh distribution, it seems that kernel estimator can effectively support the smoothing of empirical diameter frequencies. This ability may be important especially in the stage of preliminary approximations [56], as it influences the follow up selection of theoretical probability density function for dbh distribution modeling.

Uutera and Maltamo [54] reported that the use of kernel estimator with a relatively small amount of data can be seen as an advantage to this method for dbh distribution modeling, considering that traditional field inventories are often laborious and expensive. Remote-sensing techniques and remotely sensed materials, which according to some researchers drastically improve the outcome of inventories [27], also have some limitations related to dbh distribution modeling [20] and still need to be validated by field measurements. As forest inventories and planning require cost reduction [95], it seems that the application of kernel estimator can fit in this scope. Siipilehto et al. [96] and Siipilehto and Rajala [97] reported promising examples of dbh distribution smoothing by using a pre-harvest measurement tool termed as EMO. The practical application of EMO requires the measurement of several diameters in the stand and the subsequent application of kernel estimator, to obtain the diameter frequency.

The results obtained during the approximation of dbh frequency confirmed that KE1, KE2, KEA, and JSB function are slightly different models. They differ in the sensitivity to the structure of the data and behave differently depending on the data provided for their development. KE1 and KE2 to some extent, highlight local extrema (minima and maxima) in diameter frequency. KEA and JSB are susceptible to changes in kurtosis and skewness. The similar performance of KE1, KE2, KEA, and JSB function in the case of test sets can be partially explained by the number of trees that constitute training sets (several dozen or even several hundred trees) used in the stage of model development. The performance of the investigated models developed with training sets containing fewer numbers of trees seems to be an interesting topic for further analysis.

Our investigation confirmed that the final evaluation of dbh distribution models requires adequate procedure of data resampling. As shown in Table 2 and Figure 1, the evaluation conducted on the data set used for model construction and on independent data sample yields completely different results. From the practical point of view, more valuable information about model performance can be obtained with test sets. Diameter records gathered in the test sets include new information that should be properly presented by the dbh distribution model.

Despite the different methodological grounds, both kernel estimator and JSB function are suitable to describe empirical diameter frequencies of black alder stands. A slight disadvantage of the JSB function is the lack of possibility to estimate the parameters for 48 training sets. It is difficult to unambiguously define a situation in which kernel estimator or JSB function is preferred. Considering the importance of information about diameter frequency for forestry practice and science, it seems that both these methods should be treated as complementary to each other unless the analysis purpose, or data set characteristics (size, variability, etc.), discriminate against one of them.

5. Conclusions

Kernel estimator with three bandwidth variants, namely KE1, KE2, and KEA, and JSB function represent the nonparametric and parametric approaches of density estimation and are useful tools for dbh distribution modeling in alder stands. The KE1 model is characterized by a very high degree of mapping the shape of empirical diameter distributions. KE1 can accurately represent the irregularities in the shape of distribution observed in the sample used for estimation; this can be considered as an advantage for the estimation of multimodal irregularly shaped distributions. However, the use of KE1 carries the risk of excessive adjustment of the model to the distribution observed in the sample, which implies that the obtained model may differ from the actual dbh distribution occurring in a given population of trees. Models developed using KEA and JSB function provide a more general shape of distribution and do not adjust to local irregularities observed in the sample; thus, ensuring greater compatibility with real-world distribution. Reliable assessment of diameter distribution models requires consideration of the data structure. Such an assessment should be performed on an independent set of measurement data, and conducting it only on the data used in the development of the model can lead to overly optimistic conclusions.