RESEARCH ARTICLE
Fitting diameter distribution models to data from forest inventories with concentric plot design
Universidad Politécnica de Madrid, ETSI Montes, Forestal y del Medio Natural, Ciudad Universitaria s/n, 28040, Madrid, Spain
Present address: ELGO Dimitra, Forest Research Institute, 57006, Vasilika-Thessaloniki, Greece
Umeå University, Dept. Mathematics and Mathematical Statistics, SE-901 87 Umeå, Sweden.
Abstract Aim of study: Several national forest inventories use a complex plot design based on multiple concentric subplots where smaller diameter trees are inventoried when lying in the smaller-radius subplots and ignored otherwise. Data from these plots are truncated with threshold (truncation) diameters varying according to the distance from the plot centre. In this paper we designed a maximum likelihood method to fit the Weibull diameter distribution to data from concentric plots. Material and methods: Our method (M1) was based on multiple truncated probability density functions to build the likelihood. In addition, we used an alternative method (M2) presented recently. We used methods M1 and M2 as well as two other reference methods to estimate the Weibull parameters in 40000 simulated plots. The spatial tree pattern of the simulated plots was generated using four models of spatial point patterns. Two error indices were used to assess the relative performance of M1 and M2 in estimating relevant stand-level variables. In addition, we estimated the Quadratic Mean plot Diameter (QMD) using Expansion Factors (EFs). Main results: Methods M1 and M2 produced comparable estimation errors in random and cluster tree spatial patterns. Method M2 produced biased parameter estimates in plots with inhomogeneous Poisson patterns. Estimation of QMD using EFs produced biased results in plots within inhomogeneous intensity Poisson patterns. Research highlights: We designed a new method to fit the Weibull distribution to forest inventory data from concentric plots that achieves high accuracy and precision in parameter estimates regardless of the within-plot spatial tree pattern. Additional keywords: expansion factors; forest growth and yield; National Forest Inventory; spatial point pattern; Weibull. Abbreviations used: DBH (Diameter at Breast Height); EF (Expansion Factor); MAPE (Mean Absolute Percentage Error); ML (Maximum Likelihood); MLF (Maximum Likelihood of the Full sample); MLT (Maximum Likelihood ignoring Truncation); NFI (National Forest Inventory); QMD (Quadratic Mean Diameter) Authors´ contributions: Simulations, drafting of the manuscript: NN. Mathematical and statistical methods, manuscript revision: SSdL. Supplementary material (Figures S1, S2, S3 and S4) accompanies the paper on FS’s website Citation: Nikos Nanos, N.; Sjöstedt de Luna, S. (2017). Fitting diameter distribution models to data from forest inventories with concentric plot design. Forest Systems, Volume 26, Issue 2, e01S. https://doi.org/10.5424/fs/2017262-10486 Received: 14 Oct 2016. Accepted: 05 Jul 2017. Copyright © 2017 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution (CC-by) Spain 3.0 License. Funding: The authors received no specific funding for this work. Competing interests: The authors have declared that no competing interests exist. Correspondence should be addressed to Niko Nanos: nikolaos.nanos@upm.es |
CONTENTS |
IntroductionTop
National Forest Inventories (NFIs) were designed and executed in Nordic countries and in the United States during the beginning of the 20th century (Kangas & Maltamo, 2006; Zeng et al., 2015). Initially, NFIs were perceived as a means to provide quantitative information on national timber resources whilst collected data were used primarily by state authorities, forest policy makers and timber companies, to forecast timber supply and design forest policies. Lately, however, NFI planners are adapting the inventory design to satisfy newly emerging societal requirements such as biodiversity conservation and planetary CO2 budget. Thus, not only are NFI data employed in monitoring forest carbon stocks and biodiversity (Montero et al., 2005; Winter et al., 2008) but also planning of NFIs is being coordinated on transnational scales to meet the challenges faced in modern globalized societies (McRoberts et al., 2010a). Finally, NFIs provide data-based evidence to a broad spectrum of scientific fields that need up-to-date information on broad spatial scales (i.e. species migration to climate change (Montoya et al., 2008), growth and yield modeling (Nanos et al., 2002; Bravo et al., 2011; Weiskittel et al., 2011; Sharma & Breidenbach, 2015), survival after forest fires (González et al., 2007) or forest pests and diseases outbreaks (Wulff et al., 2006).
In order to diminish the inventory´s budget and reduce efforts in field measurements, most countries with boreal and temperate forests use circular plots with multiple concentric fixed-area subplots (McRoberts et al., 2010b). This type of plot design is placing more emphasis on larger trees (the ones with the most timber) whilst smaller-diameter trees are inventoried when they fall within the smaller-radius subplots and ignored otherwise. In this sense plot designs using multiple, concentric and circular plots (concentric plots for convenience) provide a truncated sample of the trees of the inventory plot using a complex truncation mechanism that selects trees depending on both their distance to the plot centre and their size. Concentric plots are not the rule in NFIs but they are used in several countries, e.g. Switzerland (Böhl & Brändli, 2007), Japan (Kitahara et al., 2009) and Portugal (Barreiro et al., 2010). Plot layout and minimum (threshold) diameters may vary across countries according to budget constraints and country-specific necessities (see Henttonen & Kangas, (2015) for optimal design patterns). French authorities, for instance, are using three concentric subplots (6-, 9- and 15-m-radius plots) with minimum measured diameters of 7.5, 22.5 and 37.5 cm (Teissier du Cros & Lopez, 2009), while Spain is using four subplots with minimum measured diameter of 7.5, 12.5, 22.5 and 42.5 for subplot radius of 5, 10, 15 and 25 m respectively (Alberdi et al., 2010).
Forest inventories with concentric plot design are very useful in providing precise information on the larger trees of a forest subject to budget constraints but they provide incomplete information on the number of trees lying in low-diameter classes. This fact precludes the use of NFI data in a variety of projects and applications such as stand structural diversity studies (Del Río et al., 2003; Pommerening, 2006), the estimation of ingrowth (Adame et al., 2010) or growth and yield modeling through size-class distributions, i.e. models based on diameter distribution models (Cañadas et al., 2002; Nanos & Montero, 2002; Nanos et al., 2004a; Gorgoso et al., 2007).
Parametric models describing the distribution of Diameters at Breast Height (DBH) are very frequently used in growth and yield studies to model the actual distribution of timber into size classes and to predict its evolution under alternative silvicultural schemes (Pique-Nicolau et al., 2011). Although several parametric families have been employed to model the stand diameter distribution, the Weibull function is one of the most frequently used models due to its flexibility (Bailey & Dell, 1973; Rennolls et al., 1985; Diamantopoulou et al., 2015; de Lima et al., 2015; Sghaier et al., 2016).
Fitting the Weibull function to the full sample of trees from the plot is straightforwardly done, e.g. by the Maximum Likelihood (ML) method (Johnson et al., 1995; Casella & Berger, 2001). However, for partially unobserved data such as the ones used in concentric plots the likelihood function to be maximized need to be adjusted to take into account sample truncation. Methods of inference for censored data based on graphical methods (Balakrishnan & Kateri, 2008) or the EM-algorithm (Ng et al., 2002) do exist but they are not applicable in this particular truncated sampling design. To handle the problem of data truncation in concentric plots Mehtätalo et al. (2011) derived the likelihood under the assumption that the within-plot spatial tree pattern was completely random. However, the efficiency and accuracy in parameter estimation of this method has not been assessed previously while the hypothesis of randomness in the spatial tree pattern may limit its applicability since both random and non-random spatial tree patterns may be sampled in a NFI.
In this paper we designed a method to fit a Weibull function to diameter data from concentric plots making no assumptions about the within-plot spatial tree pattern. After deriving the likelihood, parameter estimation is performed through ML. We run a simulation study (based on 40000 simulated plots) to evaluate the new method in terms of precision and accuracy in parameter estimates. Simulated data were generated using four spatial point processes to distribute trees within the inventory plots (random Poisson, cluster process and two inhomogeneous Poisson). Our method was compared to the one proposed by Mehtätalo et al. (2011) and also to two other reference methods (one using the full data sample and another ignoring data truncation).
Material and methods Top
Method description
Imagine a fixed-area circular inventory plot (with radius r) with the full sample of trees measured for their DBH (xi, i = 1, ..., np). If we want to fit a diameter distribution model with probability density f(xi;φ) where φ is the parameter vector to be estimated we would maximize the likelihood (assuming independence):
with respect to φ, yielding the ML parameter estimate.
Forest inventory data from concentric plots, however, include diameter measurements only for a subsample of trees from the whole plot. The inventory sample is collected in the following way: First define S circular concentric subplots with radii r1< r2 < ... < rs-1 < rs= r that split the whole inventory plot into S regions (i.e. S - 1 annuli delimited by two successive subplots and the circular plot with the smallest radius, see Fig. 1). The diameter xis of the i-th tree in the s-th region (i = 1, 2, ... ns) is sampled (belongs to the inventory sample) if xis ≥ ysfor some pre-specified ys (ys is the threshold diameter of the s-th region). The diameters from region s are thus following a truncated distribution with density function fs(xis;φ,ys) = f(xis;φ )/(1 - F(ys,φ )), where F(ys,f) is the distribution function of the non-truncated diameter distribution. The likelihood is computed using the truncated density function fs (xi; φ, ys). Therefore, the likelihood to be maximized for the observed inventory sample is:
Figure 1. Concentric plot design showing subplots, regions and truncation diameters. Trees are represented by green or black circles (circle size is proportional to the tree´s diameter at breast height or DBH). Only green trees having a DBH larger than the truncation diameter were sampled. Black circles represent non-sampled trees. |
Consider, furthermore, one of the most frequently used diameter distribution models, the Weibull probability density function with low-truncation at ys:
where a and b are the shape and scale parameters, respectively. Assuming independence the likelihood is:
The ML parameter estimates are given by maximizing the likelihood or, equivalently, the log-likelihood function with respect to the unknown parameters φ = (a,b) (note that ys,s = 1,…,S are known).
Standard errors in parameter estimates
Standard errors of the ML estimates were approximated from the observed information matrix: the covariance matrix of the parameter vector is equal to the inverse of the expected information matrix (the Fisher information matrix) which has as its jt-th elements:
This matrix can be approximated by the observed information matrix Iwith elements:
An estimate of the standard error (SE) of the parameter estimate is thus given by
where Ijj is the j-th diagonal term of the matrix i.e. the inverse of I. Approximate 95% confidence intervals for parameter φj is constructed as
Application to the NFI of Spain
The inventory plot of the Spanish NFI consists of four concentric plots (s = 1,2,3,4) with minimum measured diameter equal to y1 = 75 mm, y2 = 125 mm, y3 = 225 mm and y4 = 425 mm. The log-likelihood for this plot design for the Weibull distribution is:
where is the total number of observed trees in the truncated sample. The ML estimator of b in terms of a is obtained by setting the derivative δlogL(a,b)/δ bequal to zero and solving for b yielding the following expression:
By substituting (9) into (8) we get the following expression as a function of a:
Maximizing logL(a) with respect to a using numerical algorithms (here we used the nlminb function of R (R Core Team, 2013)) yields the ML estimator â of a. The ML estimator of b is then
The method of Mehtätalo et al. (2011)
Mehtätalo et al. (2011) used weights according to tree DBH thresholds when estimating diameter distribution models from inventory samples from concentric plots. Their weights correspond to the inclusion probability (in the sample from concentric plots) of a tree given its diameter under the assumption that trees are randomly distributed (according to a homogeneous Poisson process). Specifically for the Spanish NFI plot-design the tree weights are as follows:
i. 52/252 for trees with DBH in the interval (75 mm - 125 mm]
ii. 102/252 for trees with DBH in the interval (125 mm - 225 mm]
iii. 152/252 for trees with DBH in the interval (225 mm - 425 mm]
iv. 1 for trees with DBH > 425 mm
Weights are then used to build the weighed density function of the observed data assuming that a Weibull function (with low truncation at 75 mm) has generated the sample of tree diameters. Finally, the likelihood is maximized using numerical optimization to obtain ML parameter estimates (for more details see Mehtätalo et al., 2011).
Simulation study
A simulation study was designed to assess the adequacy of the two methods, M1 and M2, in estimating parameters of the Weibull (diameter) distribution when data from concentric plots are available. The basic unit of our simulation is the -simulated- plot. This is a circular plot with 25 m radius (similar to the largest-radius subplot of the Spanish NFI).
Spatial tree pattern simulation. Within each plot we generated the spatial location of trees as a realization of a random point process X in R2 within the 25 m-radius circular plot. Four models of spatial point process were used to generate, accordingly, four simulated sets of plots (each simulated set consists of 10000 plots). Spatial point patterns were simulated in the R environment using the corresponding functions of the spatstat package (Baddeley & Turner, 2005).
- Set 1 (random tree patterns). The spatial location of trees in the circular plot was generated as a realization of a random Poisson process with constant intensity (λ) over the plot (Diggle, 2003). The intensity parameter was assigned to each plot using a uniform distribution in the interval [0.015, 0.25]. Under this parametric setting, the expected tree density within this set of plots will be larger than 150 stems/ha and smaller than 2500 stems/ha (Fig. S1 [suppl]).
- Set 2 (cluster tree patterns). Plots of Set 2 were generated using the Thomas process (Thomas, 1949). This is a cluster point process (a Cox process) resulting from the superposition of a random Poisson process of cluster centers with intensity λ (frequently called the parent process) and an effective clustering mechanism of µ points dispersed around the parent (i.e. the cluster center) with distance from the parent being radially symmetric Gaussian N(0,σ 2) (Baddeley, 2010). The Thomas-process model parameters were assigned to each plot using a uniform distribution in the interval [1, 8] for µ and [1, 16] for σ2. Note, additionally, that λ (the intensity for the parental process) was fixed to λ = 0.015. Under this parametric setting our plots will have an expected stem density varying between 0.015*1=0.015 stems/m2 to 0.015*8=0.12 stems/m2 (or between 150 stems/ha and 1200 stems/ha). Within clusters the stem density may reach its maximum expected value of 2820 stems/ha in plots with µ =8 and σ2 =1 (computed within radial distances of 3σ (=3 meters) around the cluster center that will receive the 99.7% of the events assigned to the cluster). Plots of this type may be present in coppice forest stands (see Fig. S2 [suppl]).
- Set 3 (patterns with non-constant intensity). Plots of this set were generated as realizations of a non-stationary Poisson process (called the inhomogeneous Poisson process; Diggle, 2003). The within plot intensity λ of the point pattern varied according to the following intensity function λ (x1 ,x2) of the spatial coordinates x1 and x2 (considering the plot centre as the coordinate center):
Note that this intensity function will generate a minimum stem density of 50 stems/ha in the plot center while the maximum stem density (expected value =1300 stems/ha) is produced on the plot edge (Fig. S3 [suppl]). Spatial tree patterns of this type may be found in inventory plots including a forest gap in their surface.
- Set 4 (patterns with non-constant intensity). The tree pattern in this set was generated by an inhomogeneous Poisson process with intensity function:
This intensity function generates plots with minimum stem density (0 stems/ha) in one semicircle of the inventory plot and maximum stem density (=2250 stems/ha) in the other semicircle (Fig. S4 [suppl]). This spatial tree distribution may be found in plots within mixed-species stands when one species occupies half of the plot´s area.
Assigning tree diameters. Diameters at breast height were assigned to trees according to the Weibull distribution (diameters were independently assigned to events of the point pattern). The shape and the scale parameter were assigned to plots using the uniform distribution within the interval [1.5, 4.5] for the former and [150, 600] for the latter (we do not consider plots in uneven-aged stands, having a value for the shape smaller than one, due to the high proportion of non-measured trees). After DBH assignment for each tree in each plot we computed the plot´s basal area (i.e. the sample statistic in Eq. 16). Plots with basal area outside the interval [15 m2/ha, 60 m2/ha] were removed from the set and a new simulated plot was generated until the desired set-size was reached (10,000 plots in each set).
Choosing trees to include in the sample. Simulated trees were included in the sample according to their DBH and distance from the plot centre: trees forming the sample were the ones having a minimum measured diameter of 75, 125, 225 and 425 mm for concentric circular plots of radius 5, 10, 15 and 25 m, respectively, following the plot design of the NFI of Spain.
Estimating plot attributes. Following plot simulation we computed some relevant plot-level statistics, which were subsequently used to facilitate the interpretation of results and error analysis. More specifically we computed the following quantities for each plot:
(i) The model-based estimation of the Quadratic Mean Diameter (QMD) of the simulated plot (mbQMD) was computed using the following equation:
where Γ is the gamma function and a and b is the real shape and scale parameter, respectively, of the Weibull distribution (i.e. the ones used to generate the simulated diameters).
(ii) The full-sample based estimation of the QMD (fsbQMD):
This was computed for all trees in the plot as well as for trees larger than 75 mm (the minimum measured diameter for the Spanish NFI).
(iii) The truncated-sample QMD estimation using expansion factors (efQMD):
where hi is the expansion factor of the i-th measured diameter (see also Husch et al., 2003):
(iv) The full-sample based plot´s basal area.
The plot´s basal area (in m2/ha) was computed using the full-sample DBH data:
(v) Tree density measures
The tree density for each plot was expressed in stems/ha and also in stems in the whole plot. In addition we computed the number of stems in the truncated sample. Also, we computed the number of trees with DBH>75 mm (the minimum measured diameter in the Spanish NFI).
Estimating the Weibull parameters. Parameters estimates were obtained in each plot using the two methods described previously (referred to as M1 and M2). Method M1 is the method presented in this paper and method M2 corresponds to the one presented by Mehtätalo et al. (2011). For method M1 we computed approximate standard errors for the estimated parameters using the Fisher information matrix.
Furthermore we used two additional methods to estimate the Weibull parameters. The first one (referred to as method MLF) uses the full sample of np trees in the inventory plot, which in practice is not observed. The parameter estimates are found by maximizing the likelihood, being the product of the (non-truncated) Weibull density function over the np diameter values (xi):
MLF parameter estimates provided a comparison reference to evaluate the magnitude of the errors committed when truncation of data is applied in conjunction to methods M1 and/or M2 (however, MLF parameter estimates are never known in practical applications).
Finally, we estimated the Weibull parameters based on the truncated inventory data sample, ignoring the truncation and thus maximize the likelihood of the product of the densities in (17) over the n diameter values (a parameter estimation method referred to as MLT). MLT parameter estimates were used to provide a rough estimate of errors committed when data truncation is not taken into account by the analyst.
Error analysis. The relative performance of methods M1 and M2 (but also of MLF and MLT) for estimating the parameters of the Weibull distribution and the QMD was assessed using two error indices. For error indices computation, we used only plots with more than n = 10 trees in the truncated-data sample (this is a practice employed frequently in diameter distribution modeling (see, for instance, Palahí et al., 2006) to avoid large estimation errors originating from plots with small sample-size). Thus, the number of plots used for error computation varies across sets and it is always smaller than 10000.
The first error index, the relative bias in the q-th set, was computed as:
where and Y_{kq} are the estimated and the real parameter values, respectively, of the k-th simulated plot of set q (q = 1,2,3,4). Here, Mq designates the total number of simulated plots in the q-th set (with more than 10 trees in the truncated-data sample).
The second error index was the Mean Absolute Percentage Error (MAPE(q)) of the q-th set:
Error indices for QMD for estimation methods M1, M2, MLF and MLT were computed by setting (see also Eq. 12 for mbQMD computation):
and
into Eq. 18 and Eq. 19. Symbolsand stand for the estimated shape and scale parameters, respectively, for all estimation methods. In addition, akq and bkq stand for the real shape and scale parameters, respectively, of the Weibull function of the k-th plot in the q-th set.
Error indices were also computed for QMD estimation through Expansion Factors (EFs) (i.e. using Eq. 14 for QMD estimation). In this case Eq. 20 takes the form:
ResultsTop
Random patterns (Set 1 plots)
Summary statistics for Set 1 plots are presented in Table 1. In addition, Fig. 2a shows the plot distribution according to stem density, QMD and basal area. The number of plots with more than 10 trees in the truncated-data sample was 9991 in this set.
Table 1. Summary statistics of the simulated plots
Figure 2. Stand density vs. quadratic mean diameter of the simulated plots showing, in addition, stand basal area (in m^{2}/ha). (a) Plots generated using a homogeneous Poisson process. The two line graphs in this plot show two stand density management guidelines suggested for Pinus sylvestris in Sistema Ibérico (Madrigal et al., 1999, p. 121) and for Betula pendula (continuous line) in northern Spain (Madrigal et al., 1999, p. 155); (b): Plots generated using the Thomas process; Panels (c) and (d) show plots generated with an inhomogeneous Poisson process using the intensity in Eq. 10 and Eq. 11, respectively. |
The relative bias in QMD estimation of M1 and M2 (-0.27% and 0.23%, respectively) is almost negligible and very similar to the relative bias of MLF (-0.18%), based on all trees in the whole plot. Estimation through EFs is more biased (3.17%) but still reasonably small compared to the relative bias of MLT (19.21%), which completely ignores the truncation. In addition, the bias in QMD estimation of M1, M2 and EF shows no evident trends when plotted against the real shape or real scale parameter of the Weibull function (Fig. 3). The MLT estimator, however, is shown to produce conditionally biased results with larger (positive) errors for right-skewed diameter distributions (Fig. 3).
Figure 3. Error in Quadratic Mean plot Diameter (QMD) against the real shape (upper panel) and scale (lower panel) parameter of the Weibull distribution in four sets of simulated plots and for four estimation methods. M1: method presented in this paper; M2: method of Mehtätalo et al. (2011); MLT: method ignoring the truncation of diameters in concentric plots; Expansion factors: QMD estimation using expansion factors (^{ef}QMD). The bold lines are loess curves (smoothing parameter=0.4). |
The relative performance of methods M1, M2 (but also MLF and MLT) in shape and scale parameter estimation is similar to the one already reported for QMD if we consider the relative bias and the MAPE (Table 2). However, both relative bias and MAPE for the shape parameter were much larger for all methods showing just how difficult the estimation of the shape parameter is when using truncated data from concentric plots to fit a diameter distribution model.
Table 2. Relative bias and Mean Absolute Percentage Error (MAPE) in parameter estimates across four sets of simulated plots.
Approximate standard errors in parameter estimates of M1 are in agreement with the uncertainty in parameter estimates. In 9505 plots (95.13%) the real shape parameter of the distribution model lied inside the 95% confidence interval constructed through the Fisher information matrix. Also in 9364 plots (93.72%) the 95% confidence interval included the real scale parameter.
Cluster patterns (Set 2 plots)
The QMD and the stand density are somewhat smaller than in the previous set as a result of the spatial distribution of the Thomas model (Table 1). In addition, the number of plots with more than 10 trees in the truncated-data sample was 9884. The distribution of plots according to plot stem density, QMD and basal area is shown in Fig. 2b.
The relative bias in QMD of this set shows a slightly worse performance of M2 with respect to M1 (Table 1). The QMD relative bias was -0.34% and 0.98% for M1 and M2, respectively, and -0.10% for MLF. Estimation through EFs is again more biased (3.29%) than either M1 or M2, but still reasonably more accurate than MLT (relative bias equal to 18.87%). Similar performance results among methods M1, M2, MLF and MLT were observed for the relative bias and MAPE in scale and shape parameter estimation (see Table 2). In addition, the bias in QMD seems to be reasonably independent of the real scale and shape parameter of the Weibull function for methods M1, M2 and EF (see Fig. 3 for Set 2 plots).
Regarding the standard error estimation in parameters of method M1, the results of the simulations showed that the 95% confidence interval included the real shape and scale parameter in 9386 (94.96%) and 9210 (93.18%) simulated plots, respectively, out of 9884 plots (having more than 10 trees in the concentric plots).
Inhomogeneous Poisson patterns (Set 3 and Set 4 plots)
Plots with inhomogeneous intensity of the spatial tree pattern present a remarkably smaller variation in stand density and in QMD (see Fig. 2c, Fig. 2d and Table 1) as a result of the functions used to generate the within-plot intensity of the spatial tree-pattern. The number of plots having more than 10 trees in the truncated-data sample was smaller than in the previous simulated sets and varied between 6924 (Set 3) and 9033 plots (for Set 4).
In these sets the relative bias in parameter estimates of the Weibull function and in QMD is very large for methods M2, MLT and EF. More specifically, the relative bias in QMD estimation for M2 was 16.31% for Set 3 and 9.99% for Set 4 (for EF 17.36% and 11.40% for Set 3 and Set 4, respectively). On the contrary, method M1 seems to provide unbiased results with relative bias equal to -0.10% (Set 3) and -0.28% (Set 4). Not surprisingly, the bias in QMD for methods M2, MLT and EF is conditional on the shape parameter of the Weibull distribution being larger for plots with positive skewness (Fig. 3). Method M1, on the contrary, seems to provide unbiased estimates of QMD regardless of the real shape (or scale) parameter.
Large relative bias was detected for the shape parameter across plots of Set 3 (=30.55%) and Set 4 (=23.23%) for M2 (much larger than the bias obtained for method M1, 10% approximately for both sets). In addition, the MAPE in QMD for M2 was 17.09% and 11.24% for Set 3 and Set 4, respectively, and 11.96% (Set 3) and 11.12% (Set 4) for method M1. Similar results for the MAPE in scale and shape parameter estimation are reported in Table 2. Noteworthy, using MLT in Set 3 and Set 4 plots resulted in the highest relative bias and MAPE (more than 55%).
When using the Fisher method for standard error approximation of parameter estimates the true parameter value lied inside the 95% estimated confidence interval in 6571 (94.92%) and 6285 (90.80%) plots for the shape and the scale parameter of the function, respectively, in Set 3 and in 8596 (95.16%) and 8312 (92.01%) for the shape and the scale parameter of the function, respectively,
DiscussionTop
The use of concentric plots in forest inventories is a practice that reduces substantially the inventory budget by precluding measurements of smaller-diameter trees (Tomppo et al., 2010). As we have shown in this paper, the truncated sample originating from concentric plot designs can be used in combination with methods M1 or M2 to recover the unknown diameter distribution, as long as the user is willing to assume a Weibull family of parametric distributions to have produced the sample diameters. However, the use of either method will produce larger relative bias and MAPE than the MLF method (using all the trees in the whole plot) given that they are using a truncated –and thus smaller- data sample. In addition, our results show that the –inadequate– use of method MLT will overestimate the QMD as much as 19% to 33% (depending on the spatial tree pattern considered).
The relative performance of methods M1 and M2 depends on the spatial distribution of trees in the inventory plot. Indeed, method M2 produced (slightly) smaller relative bias and MAPE when the within-plot spatial tree pattern was random (i.e. Poisson homogeneous). This result seems to stem from the fact that the method of Mehtätalo et al. (2011) is using additional information for parameter estimation by assuming uniformity in the within-plot spatial distribution of trees (a realization of a homogeneous Poisson process). Contrastingly, M1 resulted in smaller relative bias in QMD when the spatial tree pattern in the inventory plot was clustered or inhomogeneous.
In general, our results indicate that accuracy in QMD estimation is pattern-independent for M1. Indeed, relative bias in QMD (ranging from -0.1 to -0.34%) will be negligible for M1 regardless of the spatial tree pattern considered while method M2 will overestimate the QMD as much as 10-16% in plots with inhomogeneous intensity (observed in patterns of Set 3 and Set 4). Similar results were reported when estimating the QMD using EFs with relative bias ranging from 3% to 17%. Additionally, QMD overestimation via method M2 in plots of either Set 3 or Set 4 is dependent on the shape parameter of the Weibull distribution being larger for right-skewed distributions. A very similar skewness-dependent positive bias in QMD was observed when using EFs for point estimation of QMD in plots of Set 3 or Set 4. This pattern- and skewness-dependent bias is most probably related to the number of inventoried trees lying in the smallest-radii subplots (radii of 5 m and 10 m) that in Set 3 and Set 4 plots is dramatically reduced. However, lower-DBH trees are only inventoried in these subplots. If, in addition, the real diameter distribution is right-skewed then several small-diameter trees located in the outermost subplots are not inventoried due to sample truncation. Lack of small-diameter trees from the sample, however, is producing the observed overestimation of the shape parameter (and in consequence of the QMD). A similar skewness-dependent bias was also observed in the MLT estimators for all sets of plots; however, as we have shown in this paper, the use of the full-sample likelihood function (ignoring truncation) in conjunction with the truncated sample of concentric plot designs is erroneous and should be avoided.
Therefore, as long as the assumption of within plot homogeneity and randomness in the spatial tree pattern is questioned method M2 cannot be recommended for fitting a Weibull distribution to inventory data. However, the real tree pattern is never known to the inventory analyst since the sample is truncated; thus as long as auxiliary information (usually airborne) is not available to infer the spatial tree pattern in a forest (Packalen et al., 2013; Valbuena-Rabadán et al., 2016) the use of method M2 should be avoided. Instead, it seems more plausible to recommend the method presented in this paper for estimating parameters of the Weibull distribution model that will provide on average (i.e. across several spatial tree patterns) a standard accuracy and precision.
In this paper we design a method to fit the Weibull distribution model to DBH data from concentric plots and we run a simulation study based on the Spanish NFI plot design. Thus our results may not directly reflect the magnitude of the estimation errors for other NFIs. Indeed, NFI of other countries use different plot designs that may provide different estimation errors (Kitahara et al., 2009; Teissier du Cros & Lopez, 2009). In addition, our results may not represent the estimation errors to be expected when dealing with data from point patterns where the position and the marks (i.e. the diameters) may guard some relationship (Schlather et al., 2004). Instead, we have made the assumption of independence of diameters (almost always this assumption is being made to fit a diameter distribution model). However, this assumption may be violated in some cases due to tree-to-tree competition or due to microhabitat variability or spatial genetic structure (Nanos et al., 2004b; Dimov et al., 2005; Pommerening & Särkkä, 2013; Fraver et al., 2014). In addition, the method presented in this paper is building a likelihood function using the Weibull function. Thus, our results should only be considered valid for this specific distribution family while only future research on other parametric distribution models may further generalize our results.
The application of the method presented in this paper includes the possibility to be applied for NFI harmonization across different measurement occasions and across different countries (McRoberts et al., 2010 a,b; Tomppo et al., 2010; Fridman et al., 2014). The Swedish NFI, for instance, has been using different plot-designs across measurement occasions to keep the inventory up-to-date (i.e. response to newly societal requirements or because newly emerging technologies may reduce the effort -or increase the accuracy and the spatial extend- of the inventory (Fridman et al., 2014). The method proposed in this paper can be used to provide a unique framework to study the evolution of the diameter distribution of permanent inventory plots regardless of the plot design used insofar, as long as the user is willing to assume the reported errors in parameter estimates and other possible errors due to misspecification of the family of parametric density. However, the development –in this paper- of the equations to approximate the standard errors in parameter estimates may prove useful for comparisons among different measurement occasions.
Finally, the method presented in this paper can be used in a variety of projects and applications such as (i) building growth-and-yield models using the NFI data (González et al., 2007; Gorgoso et al., 2007; Palahí et al., 2007; Adame et al., 2010; Aunós et al., 2010) or (ii) assessing stand diversity indices (Chirici et al., 2012) from NFI data, (iii) estimating the number of small-diameter trees to estimate ingrowth (Mcgarrigle et al., 2011) and (iv) estimating total biomass and carbon stocks in forest vegetation (Ruiz-Peinado et al., 2011; Gobakken et al., 2012; Breidenbach et al., 2014; Ruiz-Peinado et al., 2016).
As conclusions, fitting the Weibull diameter distribution model to data from concentric plots can be achieved through the method presented in this paper. In addition, the method to approximate standard errors in parameter estimates can prove a trustful means for uncertainty assessment. Relative bias seems to be small (negligible) in estimating the QMD and also in the estimation of the scale parameter of the diameter distribution model. In contrast larger errors are reported in the shape parameter estimation. The method presented in this paper exhibits stable performance results across all spatial tree patterns studied (random, clustered and patterns with inhomogeneous intensity). The method of Mehtätalo et al. (2011) is very efficient in reducing the bias in parameter estimates in plots with random patterns. However, its efficiency is dramatically reduced when the spatial tree pattern in the inventory plot has an inhomogeneous intensity. Therefore, fitting a diameter distribution model (of the Weibull family) to data from concentric plots should be produced using the method presented in this paper, unless the user is willing to assume random placement of trees within the inventory plot.