National Forest Inventories (NFIs) were designed and executed in Nordic countries and in the United States during the beginning of the 20th century (
_{2}
budget. Thus, not only are NFI data employed in monitoring forest carbon stocks and biodiversity (

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 (

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 (

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 (

Fitting the Weibull function to the full sample of trees from the plot is straightforwardly done, e.g. by the Maximum Likelihood (ML) method (

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

Imagine a fixed-area circular inventory plot (with radius
_{i}
,
_{p}
). If we want to fit a diameter distribution model with probability density
_{i}

with respect to

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
_{1}
<
_{2}
< ... <
_{s-1}
<
_{s}
=
_{is}
of the
_{s}
) is sampled (belongs to the inventory sample) if
_{is}
=>
_{s}
for some pre-specified
_{s}
(
_{s}
is the threshold diameter of the

_{s}
_{i}
_{s}

Consider, furthermore, one of the most frequently used diameter distribution models, the Weibull probability density function with low-truncation at
_{s}
:

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 f = (a,b) (note that ys,s = 1,…,S are known).

Standard errors of the ML estimates f = (a,b) were approximated from the observed information matrix: the covariance matrix of the parameter vector f 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 I (f) with elements:

An estimate of the standard error (SE) of the parameter estimate is f
_{j}
thus given by

where I
^{jj}
is the
*j*
-th diagonal term of the matrix I(f)
^{-1}
,
_{j}
is constructed as

The inventory plot of the Spanish NFI consists of four concentric plots (s = 1,2,3,4) with minimum measured diameter equal to y
_{1}
= 75 mm, y
_{2}
= 125 mm, y
_{3}
= 225 mm and y
_{4}
= 425 mm. The log-likelihood for this plot design for the Weibull distribution is:

where n=?
^{4}
_{S=1}
n
_{s}
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
*dlogL(a,b)/db*
equal to zero and solving for

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 (

i. 5
^{2}
/25
^{2}
for trees with DBH in the interval (75 mm - 125 mm]

ii. 10
^{2}
/25
^{2}
for trees with DBH in the interval (125 mm - 225 mm]

iii. 15
^{2}
/25
^{2}
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

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 (

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 (

Set 2 (cluster tree patterns). Plots of Set 2 were generated using the Thomas process (
^{2}
) (
^{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/m
^{2}
to 0.015*8=0.12 stems/m
^{2}
(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 s
^{2}
=1 (computed within radial distances of 3s (=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;
_{1}
,x
_{2}
) of the spatial coordinates x
_{1}
and x
_{2}
(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
^{2}
/ha, 60 m
^{2}
/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 (
_{mb}
QMD) was computed using the following equation:

where G 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 (
^{fsb}
QMD):

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 (
^{ef}
QMD):

where
*
h
_{i}
*
is the expansion factor of the i-th measured diameter (see also

(iv) The full-sample based plot´s basal area.

The plot´s basal area (in m
^{2}
/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 (
*a,b*
) 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

Furthermore we used two additional methods to estimate the Weibull parameters. The first one (referred to as method MLF) uses the full sample of n
_{p}
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 n
_{p}
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,

The first error index, the relative bias in the q-th set, was computed as:

where Y
_{kq}
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, M
_{q}
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
^{mb}
QMD computation):

and

into
_{kq}
and b
_{kq}
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 (EF
_{s}
) (i.e. using

Summary statistics for Set 1 plots are presented in

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 (

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 (

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.

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 (

The relative bias in QMD of this set shows a slightly worse performance of M2 with respect to M1 (

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).

Plots with inhomogeneous intensity of the spatial tree pattern present a remarkably smaller variation in stand density and in QMD (see

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 (

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

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,

in Set 4.

The use of concentric plots in forest inventories is a practice that reduces substantially the inventory budget by precluding measurements of smaller-diameter trees (

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

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 (

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 (

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 (

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 (

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