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

Imagine a fixed-area circular inventory plot (with radius r ) with the full sample of trees measured for their DBH ( x _i , i = 1, ..., n _p ). If we want to fit a diameter distribution model with probability density f(x _i ; f ) where f is the parameter vector to be estimated we would maximize the likelihood (assuming independence):

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:

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

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. 5 ² /25 ² for trees with DBH in the interval (75 mm - 125 mm]

ii. 10 ² /25 ² for trees with DBH in the interval (125 mm - 225 mm]

iii. 15 ² /25 ² 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 ).

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,s ² ) ( 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 s ² . 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 ² to 0.015*8=0.12 stems/m ² (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 ² =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; Diggle, 2003 ). The within plot intensity ? of the point pattern varied according to the following intensity function ? (x ₁ ,x ₂ ) of the spatial coordinates x ₁ and x ₂ (considering the plot centre as the coordinate center):

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.