Point-level density.

We estimated stem density at each point with a Morisita plotless density estimator that uses the measured distances from each survey point to the nearest trees at the point location [28]. The standardized approach for the Morisita method is well-validated [28]. However, over time the survey design used by PLS surveyors changed as protocols were updated, which affects how we estimate density from the information at each point. S2 Appendix summarizes the changes in the information recorded and how we developed and applied spatially-varying correction factors to the Morisita estimator [7] to account for these changes when estimating stem density at a point. Distances from the tree to the survey point were taken to be the distance from the survey notes plus one-half the diameter of the tree.

In S1 Appendix, we detail the steps taken to either exclude points or adjust the density estimates at points where direct estimation of density was impossible or posed a risk of bias. After all removals we estimated stem density at 66,648 Illinois points, 67,072 Indiana points, 113,801 Michigan points, 226,047 Minnesota points, and 159,058 Wisconsin points (Fig 1).

We limited estimates of density to trees greater than or equal to 8 inches (20.32 cm) dbh because that is approximately the size below which surveyors tended to avoid sampling small trees. However, at many points smaller trees were reported by surveyors. Smaller trees reported by surveyors were a censored sample to some degree. We included all trees that were surveyed in our initial density estimate (including those with missing diameters), giving a raw stem density estimate whose meaning (in terms of the implicit diameter threshold) varies spatially based on how surveyors selected for tree size in a given area. We then used a spatially-varying correction factor [7] to scale the raw density estimates to a corrected stem density estimate for trees greater than or equal to 8 inches dbh.

Estimation of individual tree biomass.

We use the allometric scalings from [29] to estimate aboveground biomass (AGB) from dbh. The assignment of allometric coefficients (for simple linear regressions of log biomass (kg) on log dbh (cm)) to taxa is provided in our GitHub repository (https://github.com/PalEON-Project/PLS_products). Note that some of the 21 taxa use the same allometric equations.

Our original goal was to make use of the full set of allometric information in [29,35] to incorporate uncertainty in scaling dbh to tree biomass, using the Bayesian statistical methods provided in the PEcAn software [36] allometry module. However, even at the taxonomic aggregation inherent in our 21 taxa, there are often few allometries available for a given taxonomic group, and, in many cases, the allometries come from locations outside of our midwestern US spatial domain. Furthermore, although there are more allometries for stem biomass (component 6 in the nomenclature of [35]; note that this excludes branches) than for aboveground (component 2) or total biomass (component 1), most research focuses on aboveground biomass rather than stem biomass. As a result, it was not feasible to robustly estimate the aboveground biomass allometries with uncertainty and therefore we have omitted incorporating allometric uncertainty.

Estimation of point-level biomass and basal area.

Here we describe how we calculate biomass at each PLS point. Calculations for basal area follow an analogous process.

In the usual case of having two trees, we calculated the point-level biomass as one-half the stem density multiplied by the estimated biomass of each tree. When the two trees were from different taxa, this calculation produces point-level biomass for two taxa that were added to estimate total biomass. When both trees are from the same taxon, this is equivalent to averaging the tree-level biomass for the two trees and multiplying by stem density.

For simplicity we excluded all 3,221 points with any tree-level missing biomass values (i.e., missing diameters), although we note that it is possible to estimate (1) total biomass based on having one of two trees with available biomass and (2) taxon-level biomass from the available tree. Since one-tree points with missing biomass cannot be used for estimation, excluding two tree points with missing biomass data treats one- and two-tree points similarly, with the goal of limiting bias at the grid cell level.

To estimate biomass, we used the original density (i.e., without using the correction factors that account for size-biased sampling to estimate density of trees greater than or equal to 8 inches dbh) combined with biomass estimates for all individual trees (including those less than 8 inches dbh) to produce an unbiased biomass estimate without an explicit size threshold. We recognize that the different surveyor behavior regarding the diameter threshold introduces some imprecision, but the effect should be small given the limited contribution of smaller trees to total biomass. In contrast, for density estimation it is critical to define a diameter threshold in order to have a meaningful quantity.

Grid-level estimation.

Before statistically modeling at the 8-km grid scale, we aggregated the point-level data to the 8-km grid by averaging over point-level stem density, basal area, and biomass values for all points in a grid cell. In addition, for our statistical modeling to be able to account for the high abundance of points with either no trees or (for taxon-specific analyses) no trees of a given taxon, we also calculated the proportion of points in each grid cell with no trees. For taxon-specific analysis, we calculated the proportion of points with no trees of the taxon of interest.

Traditionally, basal area at an aggregated level has been calculated as the product of the mean density and mean tree basal area, but because of their negative correlation, this traditional approach overestimates the average values [37–39]. Our estimates of biomass and basal area in a grid cell avoid this problem by instead calculating the mean of the point-level multiplication of density and tree size [28]. Similarly, our estimate of density for a given taxon is calculated as the average of the point-level density estimates for that taxon. Furthermore, the biomass of each taxon is the mean of the point-level biomass estimates.

Statistical smoothing.

The major challenge of modeling stem density, basal area, and biomass data is that these quantities are both non-negative and continuous, with a discrete spike at zero (i.e., “zero inflated”); few statistical distributions are available for this type of data. The description below is specifically for biomass for concreteness and clarity of presentation, but the modeling structure is the same for basal area and stem density.

There are many zero-inflated models in the statistical literature, most focusing on count or proportional data [40,41]. In early efforts we considered a Tweedie model [42,43]. However, computational difficulties affected model convergence, and the Tweedie model produced poor fits. Given this, we developed a two-stage model to address the challenge of zero inflation in non-negatively valued distributions. Our model was motivated by the biological insight that local conditions may prevent a taxon from occurring locally even though the taxon may be present at high density nearby. Thus, we combine a model for “potential biomass”, which reflects the large-spatial-scale patterns in biomass with a model for “occupancy”, which reflects the propensity for a given forest stand to contain the taxon. This model allows for zero inflation because a low probability of occupancy can easily produce observations that are zero at the grid cell aggregation.

Let N(s) be the number of PLS sample points in grid cell s. Let n_p(s) be the number of points in grid cell s that have one or more trees of taxon p. Let be the average biomass for taxon p calculated only from the n_p(s) points at which the taxon is present. In other words, where i indexes the N(s) sample points within cell s. Let m_p(s) be the potential biomass process (we consider it both on the log scale and the original scale) and θ_p(s) be the occupancy process, both evaluated at grid cell s. The biomass in a grid cell can then be calculated as b_ps = θ_p(s) exp (m_p(s)), (for the case when m_p(s) is modelled on the log scale) namely weighting the average biomass in “occupied points” by the proportion of points that contain the taxon.

First consider the occupancy model. The likelihood is binomial, n_p(s)~Bin(N(s), θ_p(s)). Note that the occupancy model represents the occupancy of points within a grid cell for taxon p and that ∑_pθ_p(s)>1 because two taxa will often “occupy” the same point, since most PLS points have two trees. Next consider the (log) biomass process. We considered modeling potential biomass both on the original scale, , and on the log scale, , where the scaling of the variance by n_p(s) is the usual variance of an average. This likelihood accounts for heteroscedasticity related to the number of points at which the taxon is observed (not the number of PLS points in the cell). Finally, for total (non-taxon-specific) biomass, n_p(s) is simply the number of points with any trees.

We found that working on the log scale produced more accurate point and uncertainty estimates based on cross-validation (S3 Appendix). This improved performance likely results from 1) downweighting the influence of outliers and 2) the log-scale model inherently having its variance scale with the mean (when both are considered on the original scale), which we observe empirically in the raw grid-level data. A disadvantage of using the log scale rather than the original scale (akin to using the geometric rather than arithmetic mean) is that by smoothing on the log scale, our estimates have a downward bias when transformed back to the original scale. This occurs because outlying values are discounted on the log scale, whereas on the original scale, outlying values have a strong influence through the squared error loss inherent in the likelihood. We also note that, based on the delta method [44], the correct approximate distribution when working on the log scale is . There is no clear means of accounting for the extra m_p(s) in the denominator when fitting the potential biomass on the log scale using generalized additive modeling (GAM) software (see below). We briefly considered using gamma regression (where the potential biomass is taken to be distributed according to a gamma distribution) in order to remain on the original scale, but we encountered computational difficulties in doing so with this large a dataset when using the GAM software discussed below. Hence, we employed the log-scale model while accepting some downward bias in the resulting estimates.

This two-stage model is able to account for the zero inflation produced by structural zeros (the taxon is not present because local conditions prevent it) through the use of the occupancy model. Through the potential model, it is also able to capture the smooth larger-scale variation in biomass. And by having both component models, we can account for the differential amounts of information caused by the large number of zeros and different numbers of sampling points in each grid cell.

Note that m_p(s) is likely to be quite smooth spatially, at least for the PLS data, because when a point is occupied by a given taxon, the tree is likely to be of adult size, regardless of whether the tree is common in the grid cell. Thus, most of the spatial variation in biomass may be determined by variability in occupancy. The potential biomass is meant to correct for the fact that density and tree size may vary somewhat, but probably not drastically, across the domain.

We fit the two-component models using penalized splines to model the spatial variation, with the fitting done by the numerically-robust GAM methodology implemented in the R package mgcv [45]. We use the GAM implementation intended for large datasets encoded in the bam() function [46] in place of the usual gam() function.

The GAM methodology determines a data-driven compromise between simply using the noisy grid-level estimates (in which much of the spatial variation is statistical noise) and averaging over large spatial regions (which prevents seeing real finer-scale spatial variation). Internally during the fitting, the methodology uses a form of cross-validation (separate from the cross-validation results reported in S3 Appendix) to determine the optimal degree of spatial averaging that smooths over noise while retaining spatial signal. As can be seen in Results, this spatial averaging smooths out sharp spatial variations (likely from statistical noise) at the scale of 1–5 grid cells (several townships) while resolving larger-scale spatial features.

We accounted for the heterogeneity in the number of occupied points per grid cell by setting the ‘weights’ argument in the bam() function to equal to n_p(s). We also considered scaling all weights by dividing by 70, where 70 is the approximate number of points in a cell that was fully surveyed. This treats a fully-covered cell as having one ’unit’ of information and scales the contribution to the likelihood from cells with a different number of points relative to that. However, the results with and without the division by 70 were identical for the point estimates and very similar for the uncertainty estimates, so our final results omit this scaling.

We did not use environmental covariates as predictors in our statistical model for several reasons. First, we have fairly complete coverage (see Fig 1), such that the use of covariates is expected to provide limited additional information. Second, covariates such as climate for the settlement time period are not available, and we were reluctant to assume that present-day values are sufficiently similar to values in the past. Finally, without developing complicated statistical models that allow the effect of covariates to vary spatially (so-called varying-coefficient models), using regression coefficient estimates that are constant spatially can cause biases, such as inferring the presence of a taxon outside of its range boundary. For these reasons, we chose to rely only on spatial smoothing of the raw data. In the future, researchers could use our raw data products in combination with covariates.

Finally, to estimate total biomass, we fit the model above to raw total biomass values from the survey points, aggregating in the same fashion as described above for individual taxa, but including data from all trees. We note that the estimated total biomass at individual grid points is systematically higher than summing over the taxon-specific estimated values; this results from smoothing on the log scale, which more strongly discounts the more extreme taxon-specific values in individual grid cells than would be the case if working on the original scale.

Quasi-Bayesian uncertainty estimates.

As discussed in [45], one can derive a quasi-Bayesian approach and simulate draws from an approximate Bayesian posterior by drawing values of the spline coefficients based on the approximate Bayesian posterior covariance provided by gam() or bam() and, for each draw, calculating a draw of θ_p(s) and similarly a draw for m_p(s) for the biomass process. We combined 250 draws from the occupancy and potential biomass processes (assuming independence between the processes) to produce biomass draws for each taxon and for total biomass. The procedure for stem density and basal area was analogous.

Note that one major drawback of our methodology of fitting the taxon-specific models separately is that individual taxon estimates are not constrained to add to the total biomass values estimated from using our model on raw total biomass values because the taxa are fit individually. Further, as in our related modeling of composition [22], we do not capture correlations between taxa, in part to reduce computational bottlenecks and in part to avoid inferring the value of one taxon based on the value of another. While there are real correlations, the correlation structure likely varies substantially over space (e.g., two taxa that positively covary can have different range boundaries such that the presence of one beyond the boundary of another does not indicate the presence of the second taxa). Since information is available for all taxa at any location with data, there is little need to borrow strength across taxa to improve our point estimates (unlike the need to borrow strength across space to fill in missing areas and smooth over noise caused by limited data in each grid cell). However, since the taxa are fit independently of each other, but are not truly independent, one cannot estimate uncertainty in any quantities that are functions of two or more taxon-specific estimates. Also note that one might scale the taxon-level point estimates to add to the total estimates, but there is no clear way to do this at the level of the posterior draws, because the draws are computed independently between the total and taxon-specific fits. In summary, we are able to calculate uncertainties in the total biomass and in the biomasses for individual taxa but unable to calculate uncertainty for any variable that combines two or more of the taxon-specific variables.

In the GAM fits, we noticed some anomalies in the quasi-posterior draws for the occupancy model that were likely caused by numerical issues. In particular, for some taxa, there were draws of the occupancy probability that were more than five times as large as the corresponding point estimate of the probability. Most of these occurred for very low occupancy probabilities in areas outside the apparent range boundary for the taxa. There were also cases where draws of total (non-taxon-specific) occupancy probability were near zero even though the probability point estimate was essentially one. As an ad hoc, but seemingly effective solution, we made the following adjustments to the draws:

1. set draws where the taxon-specific occupancy probability is greater than five times the point estimate to be equal to the point estimate, and
2. set all draws where the total point estimate was greater than 0.999 to be equal to 1.

Choice of smoothing and scale of averaging.

We used cross-validation for total biomass, basal area, and stem density, as well as on a per-taxon basis, at the grid scale to:

1. choose between estimating potential biomass on the log-scale or original scale, and
2. determine the maximum number of spline basis functions, denoted as k.

With regard to the maximum number of basis functions, while the generalized additive modeling methods of [45] choose the amount of smoothing based on the data, using a large number of basis functions can result in slow computation. We hence chose to limit the number of basis functions and thereby impose an upper limit on the effective degrees of freedom of the spatial smoothing estimated from the data. The choice of that limit was informed by cross-validation. However, the imposed upper limit to the number of possible basis functions was large enough to have little effect on the amount of smoothing, although possibly imposing slightly more smoothing than without the limitation.

We used 10-fold cross-validation, randomly dividing the grid cells into 10 sets and holding out each set in turn. This allows us to assess the ability of the model to estimate biomass for cells with no data (and also gives us a good sense of performance for cells with very few points). Note that even with our incomplete sampling in Indiana and Illinois (Fig 1, 65% and 42% of the total area of each state, respectively), most unsampled grid cells are near to other grid cells with data. We considered values of k in the set {100, 250, 500, 1000, 1500, 2000, 2500, 3000, 3500}.

The metrics used in cross-validation were absolute error loss for the point predictions relative to the grid-level raw data and statistical coverage of prediction intervals of the grid-level raw data.

We calculated absolute error weighted by the number of PLS points in the held-out cell and truncated both held-out values and predictions to maximum values of 600 Mg/ha to avoid having very large values overly influence the assessment. This also allowed us to work on the original (not log) scale in our evaluation, as we didn’t want to accentuate small differences at low biomass values when choosing amongst modeling approaches, although we did end up choosing to fit the model on the log scale based on this empirical evaluation.

We calculated 90% prediction interval coverage using a modified version of the quasi-Bayesian uncertainty procedure described previously. The modification to the sampling procedure involves drawing a random binomial value based on each draw of the occupancy probability, multiplied in turn by a random normal draw (exponentiated when fitting the model on the log scale) centered on the draw of the potential surface with variance equal to the residual variance from the potential model. The addition of the binomial draw and the residual variance produces a prediction interval for the data rather than the unknown process and allows us to assess coverage relative to an observed quantity. We calculated 90% prediction intervals using the 5th and 95th percentiles of the 250 draws for each held-out cell. Coverage was determined as the proportion of cells for which the observation fell into the interval, considering only grid cells with at least 60 PLS points. We also calculated the median length of intervals (and median log-length) to assess the sharpness of the intervals, as high coverage can always be trivially obtained from overly-wide intervals.

Unfortunately, the coverage results cannot directly assess the uncertainty estimates provided in the gridded data product. This is because the true biomass is unknown and thus cannot be used to judge coverage. We can only judge coverage of prediction intervals for the data. Thus, under- or over-estimation of uncertainty for the true quantities may be masked by compensating over- and under-estimation of the residual error of the data around the truth.

Based on the cross-validation, we chose to use the models where the potential biomass (or stem density or basal area) is fit on the log scale, as this approach produced lower absolute error loss and gave a much better tradeoff between coverage and interval length (S3 Appendix).

Data product.

We provide the following data products via the LTER Network Data Portal.

• Raw gridded stem density, aboveground biomass and basal area at https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=26&revision=1 (DOI: 10.6073/pasta/801601af769fa5acade1ef07f6892bdd).
• Gridded statistically-smoothed stem density at https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=24&revision=0 (DOI: 10.6073/pasta/1b2632d48fc79b370740a7c20a70b4b0).
• Gridded statistically-smoothed aboveground biomass at https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=23&revision=0 (DOI: 10.6073/pasta/b246e05afb25dbe06b3006c5d18a4a2b).
• Gridded statistically-smoothed basal area at https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=25&revision=0 (DOI: 10.6073/pasta/c3ae2363e4ae2e0f42a7c02b6f12b50a).

We also provide point-level raw data:

• Indiana - https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=27&revision=0 (DOI:10.6073/pasta/c3e2404f5b34204b5871a743ebce3c51)
• Illinois - https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=28&revision=0 (DOI: 10.6073/pasta/b8fdabd7cfba3a2b3d55fe4c1dc5383f)
• Michigan -
  1. ○. Southeastern - https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=29&revision=0 (DOI: 10.6073/pasta/409ec6dfb218b6a3e98022916d2b4438)
  2. ○. Southern - https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=30&revision=0 (DOI: 10.6073/pasta/8d033c1cfadca42bf060f9f38940c81e)
  3. ○. Northern - https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=31&revision=0 (DOI: 10.6073/pasta/3760eec82562e0a8b7cd493c0a3e3ef4)
• Wisconsin - https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=32&revision=0 (DOI: 10.6073/pasta/c3e680e51026e74a103663ffa16cb95d)
• Minnesota - https://portal.lternet.edu/nis/mapbrowse?scope=msb-paleon&identifier=33&revision=0 (DOI: 10.6073/pasta/f55f6b7f4060a9b4f07374e7db8443cd)

The project GitHub repository (https://github.com/PalEON-Project/PLS_products) provides code for processing the point-level data and producing the data products above in the subdirectory named ‘R’. In the subdirectory ‘data/conversions’, we provide:

• our translation tables for translating surveyor taxon abbreviations to modern common names, including aggregation for the raw gridded values and statistical modeling done in this work,
• correction factors for the subregions of the domain for estimating point-level tree density, and,
• our assignments of allometric relationships for the PalEON taxa, based on [29].