Regional Modeling of Forest Fuels and Structural Attributes Using Airborne Laser Scanning Data in Oregon

: Airborne laser scanning (ALS) acquisitions provide piecemeal coverage across the western US, as collections are organized by local managers of individual project areas. In this study, we analyze different factors that can contribute to developing a regional strategy to use information from completed ALS data acquisitions and develop maps of multiple forest attributes in new ALS project areas in a rapid manner. This study is located in Oregon, USA, and analyzes six forest structural attributes for differences between: (1) synthetic (i.e., not-calibrated), and calibrated predictions, (2) parametric linear and semiparametric models, and (3) models developed with predictors computed for point clouds enclosed in the areas where ﬁeld measurements were taken, i.e., “point-cloud predictors”, and models developed using predictors extracted from pre-rasterized layers, i.e., “raster-ized predictors”. Forest structural attributes under consideration are aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load. Results from our study indicate that semiparametric models perform better than parametric models if no calibration is performed. However, the effect of the calibration is substantial in reducing the bias of parametric models but minimal for the semiparametric models and, once calibrations are performed, differences between parametric and semiparametric models become negligible for all responses. In addition, minimal differences between models using point-cloud predictors and models using rasterized predictors were found. We conclude that the approach that applies semiparametric models and rasterized predictors, which represents the easiest workﬂow and leads to the most rapid results, is justiﬁed with little loss in accuracy or precision even if no calibration is performed.


Introduction
Forest fires and carbon accounting are tightly interrelated areas of interest with critical importance for sustainable forest management [1,2]. Detailed spatial information about forest fuels and biomass accumulation in the forest are necessary to guide decision-making processes in these areas [3][4][5]. However, forested areas are typically large and remote, and obtaining this information using only field measurements is an expensive or inefficient option for many applications. For example, national forest inventories, such as the US 1.
Transferability and effect of calibration. Comparisons focused on analyzing differences in accuracy and precisions between synthetic predictions and predictions obtained using the same models but performing an additional calibration step with available ground and ALS data. We will refer to this factor in the following sections as "calibration". 2.
Differences between modeling techniques. Comparisons focused on analyzing differences between parametric linear mixed-effects models and semiparametric models. We will refer to this factor in the following sections as "modeling technique".
3. Differences between models using a different source of ALS metrics. Comparisons focused on analyzing differences between models using ALS predictors computed from point-clouds clipped around the training plot footprints, i.e., point-cloud predictors, and models using predictors extracted from raster layers, i.e., rasterized predictors. We will refer to this factor in the following sections as "source of predictors".

ALS Data Acquisitions and ALS, Climate, and Topographic Metrics
Eight ALS data acquisitions collected by different agencies in the state of Oregon during the period 2008-2016 were used in this study (Figure 1). These acquisitions covered the main forested areas in the state and included areas of temperate coastal coniferous forest, areas with Mediterranean influence in the south of the state, mountain areas on the Cascades range, and drier and more continental forest ecosystems east of the Cascades mountains. The area, completion year, flying altitude, sensor information, return density and number of field plots available for each ALS acquisition are indicated in Table  1. Figure 1. Study area and airborne laser scanning (ALS) data acquisitions used in the analyses. Individual ALS acquisitions can be identified using Table 1.  Individual ALS acquisitions can be identified using Table 1.
For all acquisitions, 30 m resolution rasters containing the ALS metrics indicated in Table 2 were available. The ALS predictors only included descriptors of the distribution of the ALS point cloud such as percentiles or moments and proportions of returns in height categories. For each FIA plot we obtained two sets of ALS predictors. The first set, "rasterized" predictors, was obtained by intersecting the FIA plot center location with the rasters containing the ALS metrics. The second set, "point-cloud" predictors, was obtained by first clipping the points inside the four macro-plots of each FIA plot. Then, for each FIA plot we normalized the point cloud with digital terrain models (DTMs) provided by the vendor of each ALS data acquisition using the R-package lidR [47]. Once normalized, metrics for each FIA plot were computed using FUSION [48]. An extended analysis for the Remote Sens. 2021, 13, 261 5 of 35 factors: (1) transferability and effect of calibration, and (2) modeling technique, including 12 additional ALS acquisitions where only rasterized predictors were used is presented in Appendix A. Total area and number of FIA plots 46,146.51 1133 In addition to the ALS predictors, we used topographic and climate predictors listed in Table 2 in the models for the selected responses. Topographic predictors were derived from a DTM derived from the Shuttle Radar Topography Mission and climate predictors were obtained from the Climate-FVS Ready Data Server [49]. Both topographic and climate indexes were rasterized at a 30 m resolution and a grid that aligned with the grids containing ALS metrics. For each FIA plot, topographic and climate predictors were obtained intersecting the FIA plot center with the corresponding raster layers. Table 2. Sets of candidate predictors used in the study.

Group Description Auxiliary Variables Acronym
ALS predictors (Derived from point clouds or extracted from raster files) Mean, standard deviation, skewness coefficient, kurtosis coefficient, coefficient of variation of the distribution of heights of the point cloud.

Ground Data and Response Variables
Ground data to train models was obtained from the FIA database [50]. Coordinates of FIA plots were obtained using mapping-grade GPS units. To the best of our knowledge no study has analyzed directly the reliability of the GPS coordinates of the FIA database in the study area, but based on previous experiences [51,52] the errors are expected to have accuracies in the range of a meter and maximum location errors are expected to be in the range of 5 to 10 m [53]. Following [26], only those FIA plots that were measured at most three years before or after each ALS data acquisition project was completed were used. In addition, any plots presenting any sign of disturbance between the plot visit date and the completion date of their corresponding ALS data acquisition were removed from the dataset. For each one of the remaining FIA plots, AGB per hectare were computed by aggregating tree-level AGB provided in the FIA database with their corresponding expansion factors [50] (Table 3). Estimates of DWB per hectare were obtained for each FIA plot condition class and weighed by the transect length of the appropriate condition class. Finally, values of CBD, CH, CBH and CFL were obtained for the FIA plots using the FireCalc program [54] and the tree-lists of the FIA plots as inputs (Table 3).

Parametric Models
Parametric models for AGB, DWB, CBD, CH, CBH, and CFL, were linear mixed-effects models with random intercepts and slopes for each ALS data acquisition. The general form of these models is where x ij is a p-dimensional vector with the first element a 1 and the p-1 remaining elements being the values of p − 1 covariates for the j th plot in the i th ALS data acquisition, β is a p-dimensional fixed-effect parameter vector, v i is a p-dimensional vector of random effects for the i th ALS data acquisition, and e ij is an additive model error. The random effects v i were assumed to be normally distributed with a variance-covariance matrix G(δ v ) with δ v a vector of variance-covariance parameters that is the same for all ALS data acquisitions; i.e., v i ∼ N(0, G(δ v )) ∀ i. Model errors were also assumed to be normally distributed, but we allowed for nonconstant error variances. The variance of the model errors was assumed to be proportional to a power α of the predictor most correlated with the response, mcp ij . That is, e ij ∼ N(0, σ 2 e mcp 2α ij ), where σ 2 e and α are model parameters that are the same for all ALS data acquisitions. For any pair of ALS data acquisitions, i and k, random effects v i and v k were assumed to be independent so cov(v i , v k ) = 0 p x p and model errors e ij and e kl for any pair of FIA plots were also assumed to be independent.
Hereafter we will use lower-case Greek letters to denote sets of FIA plots possibly distributed across several ALS data acquisitions. Four sets will be systematically considered in the following sections. The first one, τ, will denote a set of FIA plots used to fit a model, the second one, ρ, will be a set of FIA plots or prediction points in an ALS acquisition for which predictions are sought. The third and fourth sets that will be considered are ρ s and ρ u . The set ρ s is the subset of ρ that contains the sampled elements for which both ground and auxiliary information are available. Analogously, ρ u is the subset of ρ for which only auxiliary data is known; i.e., the complement of ρ s . Letting ξ denote an arbitrary set containing n ξ FIA plots, grouped by ALS acquisitions, the model in Equation (1) can be specified in matrix notation as where X ξ = col 1≤i≤m (X i ) is a matrix with n ξ rows and p columns obtained by stacking the matrices X i = (x i1 , . . . , x ij , . . . , x in i ) t associated to the n i elements from the i th ALS data acquisition in ξ and Z ξ is a matrix with n ξ rows and mp columns formed by mxm blocks, where all elements in off diagonal blocks are zeros and blocks in the diagonal are the matrices X i . The vector of random effects is and the vector of model errors is e ξ = (e t 1 , . . . , e t i , . . . , e t m ) t with e i = (e i1 , . . . , e n i ) t . The variance-covariance matrix of y ξ is denoted as , the variance-covariance matrix of e ξ , a diagonal matrix where the elements in the diagonal equal σ 2 e mcp 2α ij . Grouping all variance-covariance components into a single vector δ = (δ t v , σ 2 e , α) t , we will simplify the notation for V ξ (δ v , σ 2 e , α) as V ξ (δ).

Model Selection
A model selection process consisting of four steps was run separately for every response and type of ALS metrics (i.e., point-cloud and rasterized predictors). In the first step we obtained the four linear fixed-effects models with highest R 2 with one, two, three, and up to seven predictors using the R-package leaps [55]. For most models we observed that residuals tended to increase with the predicted value, thus, for each candidate we obtained a second fixed-effects model where the error variance was proportional to mcp 2α ij and compared it to the first model using a likelihood ratio test. When the p-value of the likelihood ratio test was smaller than 0.05, the candidate model with constant error variance was replaced by its counterpart with error variance proportional to mcp 2α ij . Finally, for each candidate we obtained a mixed-effects model having the same fixed effects and error variance but incorporating random effects as indicated in Equation (1). The significance of the fixed-effect coefficients associated with each predictor was tested for each candidate and those coefficients that were not different from zero at a 0.05 significance level were sequentially removed from the model until all fixed effects were significantly different from zero at a 0.05 confidence level. The result was a list of 28 models from which we selected a final one for the response variable under consideration.
To select the final model, we first removed from the list of 28 candidates all models where the maximum variance inflation factor, VIF, was larger than 5. To balance parsimony and predictive power, we observed at the increases of explained variance when increasing the number of predictors. We initially removed from the list all models where the R 2 was 2.5% less than the maximum. Then, we only kept the models with the smallest number of predictors and selected the model with lowest root mean square error if more than one remained in the list. All models in this process were obtained using maximum likelihood and functions from the R package nlme [56].

Prediction and Calibration with Parametric Models
Once models were fitted, we obtained synthetic predictions and calibrated predictions for each FIA plot. Lettingδ be the estimated variance-covariance parameters of the model, we obtained the estimated fixed-effect parameters aŝ Synthetic predictions for points in a new ALS acquisition, ρ, were obtained by a direct extrapolation using the fixed-effects parameters obtained in the model fitting stage aŝ To emphasize that these are synthetic predictions, we will use the superscript syn. It is important to note that synthetic predictions do not perform any calibration to the local conditions of a new ALS data acquisition and can be obtained without any new ground information.
If a set ρ s of ground observations with their corresponding values of the auxiliary variables is available for the new ALS data acquisition, then, following p. 314 [57], it is possible to obtain calibrated predictions aŝ where the superscript cal is used to denote that these are calibrated predictions and

Semiparametric Models
For each response variable we obtained semiparametric random-forest models using point-cloud and rasterized predictors. We followed the approach proposed by [40] but included some modifications to accommodate nonconstant error variances because increasing error variances are commonly observed when modeling forest attributes with ALS auxiliary information (e.g., [58,59]).
The form of the semiparametric models can be described as where f (x ij ) is a fixed and unknown function that will be approximated by the random forest algorithm and θ ij a random variable with zero mean that includes all the variation that is not explained by f (x ij ). We further assumed that the random component not explained by f (x ij ) was where u i is a random effect specific of the i th ALS data acquisition and ε ij is an additive model error for the j th plot in the i th ALS acquisition. Random effects and model errors are assumed to be independent of each other and normally distributed with , and κ are model parameters. Finally, u i was assumed to be independent of u k for any pair of ALS data acquisitions and ε ij independent of ε kl for any pair of FIA plots.
For an arbitrary set of units, model (8) can be expressed as where col ] is a n ξ dimensional column vector obtained stacking the values f (x ij ) of all units in ξ and U ξ is a matrix with n ξ rows and m columns, where the row corresponding to the j th element of the i th acquisition has zeros everywhere except for the i th position, which has a value of f (x ij ). The vector u ξ = (u 1 , . . . , u i , . . . , u m ) t groups the random effects for all ALS acquisitions in ξ and the vector ε ξ is a vector of model errors obtained in the same way as e ξ . The variance-covariance matrix of u ξ is J ξ (σ 2 u ) = σ 2 u I mxm with I mxm an identity matrix of dimension m and the variance-covariance matrix of Variance-covariance parameters σ 2 u , σ 2 ε , and κ were obtained after a random forest model was fitted to the training sample τ. Then we assumed that the random forest provided a close approximation,f (.), to the unknown function f (.). Usingf (.) we obtained θ ij = y ξ −f (x ij ), which were assumed to have a normal distribution with zero mean and variance covariance matrix as indicated in Equation (10). The estimated residuals from the random forest model,θ ij , were finally used to estimate σ 2 u , σ 2 ε , and κ using maximum likelihood.

Model Selection
A model selection consisting of three steps was developed for every response variable and source of ALS metrics. In the first step we followed the approach described by [28] and eliminated correlated predictors using the QR decomposition implemented in the multi.collinear function of the rf.Utilities R package [60]. This step was run only once for each combination of response variable and source of ALS metrics. In the second step, we used the function rf.modelSel [60] to identify the best combination of variables selected in the final random-forest model for each response variable and source of ALS metrics. Once predictors were selected, random-forest models were fit using the randomForest R package [61]. These models provided the nonparametric component of the semiparametric models and were used to compute values forθ ij . In the last step, we obtained the parametric part of the model that e×plains differences inθ ij due to ALS acquisition membership. Two models as indicated in Equation (10) were obtained for the variance-covariance of random effects and model errors using the R package nlme [56]. The first model had a fixed exponent κ = 0; i.e., homoscedastic errors, and the second model had an exponent κ that was free to vary during the model-fitting stage. Both models were compared using a likelihood ratio test. We selected the model with variable κ when the p-value of this test was smaller than 0.05; otherwise, we selected the model with κ = 0.

Prediction and Calibration with Semiparametric Models
A direct application of the random forest model to a new ALS data acquisition provides synthetic predictionsŷ Following [41] (pp. 307, 353), calibrated predictions can be obtained aŝ 2.5. Accuracy Assessment and Comparisons 2.5.1. Cross-Validation and Performance Metrics To assess the effect of calibration for new ALS data acquisitions, we used the following procedure for each response variable and modeling technique ( Figure 2). First, for each ALS data acquisition, we split the entire training dataset in two parts, τ and ρ. The set τ, the training subsample, contained all FIA plots not included in the selected ALS data acquisitions, and the set ρ contained the FIA plots within the ALS data acquisition under consideration. Models obtained in the model selection stage were re-fitted using only τ and synthetic predictions were obtained for all elements in ρ. Finally, random effects and calibrated predictions for parametric and semiparametric models were obtained using Equations (5) and (6) and Equations (12) and (13), respectively. To consider that only FIA plots measured prior to the acquisition date will be available for the calibration of new ALS data collections, we performed the calibration using ρ s , the subset of ρ that contained the FIA plots that were measured the year the ALS data acquisition was finished or the three previous years. It is important to note that while models were selected including also plots measured up to three years after the ALS data collection was finished, we tested the effects of the calibration by excluding the plots that were measured after the ALS data collection was finished when computing the random effects using Equations (6) and (13). The potential use of the models developed in this study is to perform calibrations in new ALS data acquisitions. Thus, excluding the plots that were measured after the acquisitions were finished provides a closer approximation to the sample sizes that could be used for calibration in a real case. The result of this process was (1) a set of synthetic predictions from models developed with data from different acquisitions and (2) a set of calibrated predictions where only the field plot data available at the end of the ALS acquisition project was used to compute random effects.
We used both sets of predictions to computeê ij =ŷ ij − y ij and obtained the following performance metrics: where n i represents the number of elements in the i th ALS data acquisition andŷ ij can represent synthetic or calibrated predictions. In addition to these metrics we also computed their values relative to the mean, y, of the response variable under consideration as: And the coefficient of determination as:

Differences between Calibrated and Synthetic Predictions, Modeling Techniques, and Sources of ALS Metrics
For each response variable, we analyzed changes in performance due to each of the three factors under analysis ( Figure 2) by computing changes in the relative root mean square error and relative bias. For a given factor, changes were computed as: and : where f , A, and B are generic superscripts that respectively indicate the factor under analysis and the two alternatives that are compared. For the transferability analysis we used the symbol t for f , and A = cal for calibrated predictions and B = syn for synthetic predictions. To analyze differences due to the modeling technique we used error metrics computed after performing the corresponding calibration. We used the symbol m to indicate the factor "modeling technique", and A = par for the parametric models and B = spar for the semiparametric models. Finally, for the comparisons of models using rasterized and point-cloud ALS predictors, error metrics were also computed after performing the corresponding calibration. The letter s was used to indicate the factor "source of predictors", and A = r for models using rasterized ALS metrics, and B = pc for models using point-cloud predictors.
To test differences in model performance due to: (1) the modeling technique and (2) the source of ALS metrics, we used a t-test on the differencesê ij A −ê ij B to test differences in accuracy and a t-test on the differences ê ij A − ê ij B to assess differences in precision.
Finally, in addition to the global values of SE, RBI AS, ∆RMSE f , and ∆RBI AS f , these metrics were also computed for each ALS acquisition separately.

Parametric and Semiparametric Models
Regardless of the source of ALS metrics and modeling technique, several patterns were observed with respect to the predictive performance of the models for different response variables. The largest R 2 values were obtained for CH, followed by AGB, with values ranging from 88.59% to 85.29% and from 80.03% to 76.45%, respectively. For the remaining variables the explanatory power of the models was largest for CFL with R 2 values above 50%. This variable was followed by CBH and CBD, which tended to have R 2 values around 30% and 25% respectively. Finally, the explanatory power of the models for DWB was very poor, and larger for the parametric models where R 2 reached 13.79% for the model using point-cloud metrics and 13.49% for the model using rasterized predictors. Patterns observed for RRMSE cal were similar to those observed for R 2 but in the opposite direction. Finally, RBI AS cal was positive for most models indicating a systematic overprediction with a magnitude that varied substantially between variables. For CH and AGB it was below 2.5% in absolute value. For CBD, CBH, and CFL, RBI AS cal was below 10% in absolute value and for DWB RBI AS cal exceeded 5% for the parametric models and 15% for the semiparametric models ( Figure 3).
for the model using point-cloud metrics and 13.49% for the model using rasterized predictors. Patterns observed for were similar to those observed for but in the opposite direction. Finally, was positive for most models indicating a systematic overprediction with a magnitude that varied substantially between variables. For CH and AGB it was below 2.5% in absolute value. For CBD, CBH, and CFL, was below 10% in absolute value and for DWB exceeded 5% for the parametric models and 15% for the semiparametric models ( Figure 3).

Figure 3.
Relative root mean squared error (RRMSE cal ), relative bias (RBIAS cal ), and coefficient of determination (R 2 ) for parametric models. Error metrics were calculated for models fitted with the entire dataset under consideration and predictions calibrated using all plots in the corresponding ALS acquisition, which were measured either in the year that the ALS data acquisition was finished, or in the two previous years. Model response variables are: aboveground biomass (AGB), downed woody biomass (DWB), canopy bulk density (CBD), canopy height (CH), canopy base height (CBH), and canopy fuel loading (CFL). Relative root mean squared error (RRMSE cal ), relative bias (RBIAS cal ), and coefficient of determination (R 2 ) for parametric models. Error metrics were calculated for models fitted with the entire dataset under consideration and predictions calibrated using all plots in the corresponding ALS acquisition, which were measured either in the year that the ALS data acquisition was finished, or in the two previous years. Model response variables are: aboveground biomass (AGB), downed woody biomass (DWB), canopy bulk density (CBD), canopy height (CH), canopy base height (CBH), and canopy fuel loading (CFL).
All models included random effects; however, for the parametric models for DWB and CBH using point-cloud predictors, the inclusion of random effects did not improve the model fit significantly at a 0.05 confidence level. For CBH the model selection algorithm for the parametric models described in Section 2.3.1 provided models with up to six predictors, followed by models with three predictors but with an R 2 3.2% points larger than the maximum R 2 . These simpler models were considered more appropriate and selected as final models for this variable because they allowed for a substantial model simplification without causing important performance losses. For most variables, sources of ALS metrics, and modeling techniques, models improved when including random effects for the ALS acquisitions and a parameter κ = 0 to account for nonconstant error variances (Tables A4, A5 and A7). The only exceptions to this trend were the parametric models for CFL and the parametric model for CH using rasterized ALS metrics. Important differences were observed between parametric and semiparametric models for the number of predictors included in the model. The median number of predictors for the semiparametric models was 33, and ALS, topographic, and climate metrics were present in all but the semiparametric model for DWB using rasterized predictors. The median numbers of ALS metrics, topographic, and climate metrics were 18, 10, and 4 variables respectively ( Figure A5). Parametric models always had less than five predictors and in most cases these predictors were derived from the ALS data (Tables A4 and A5).

Differences between Synthetic and Calibrated Predictions
All models included ALS acquisition random effects; however, the magnitude of these effects varied significantly between variables and modeling techniques. For semiparametric models, reductions in global RMSE and RBI AS were negligible, which indicates that with these models using synthetic predictions and omitting the calibration step will not have consequences of practical importance in accuracy and precision. While semiparametric models seem to capture differences between ALS acquisitions directly, leaving little room for improvements to the calibration stage, parametric models operate differently, and their performance improves substantially with the calibration. Parametric models captured the main relationship with auxiliary information in their fixed-effect component, and a part of the variability between ALS acquisition was accounted for in the calibration stage ( Figure 4). Focusing on the parametric models, for AGB, CBD, CBH, and CFL calibration resulted in consistent reductions of the global RMSE cal and RBI AS cal . For DWB and CH, three parametric models failed at reducing the global RBI AS cal . We further inspected this issue by looking at ALS-acquisition-specific values of RBI AS cal for the parametric models. the model fit significantly at a 0.05 confidence level. For CBH the model selection algorithm for the parametric models described in Section 2.3.1 provided models with up to six predictors, followed by models with three predictors but with an R 2 3.2% points larger than the maximum R 2 . These simpler models were considered more appropriate and selected as final models for this variable because they allowed for a substantial model simplification without causing important performance losses. For most variables, sources of ALS metrics, and modeling techniques, models improved when including random effects for the ALS acquisitions and a parameter ≠ 0 to account for nonconstant error variances (Tables A4, A5 and A7). The only exceptions to this trend were the parametric models for CFL and the parametric model for CH using rasterized ALS metrics.
Important differences were observed between parametric and semiparametric models for the number of predictors included in the model. The median number of predictors for the semiparametric models was 33, and ALS, topographic, and climate metrics were present in all but the semiparametric model for DWB using rasterized predictors. The median numbers of ALS metrics, topographic, and climate metrics were 18, 10, and 4 variables respectively ( Figure A5). Parametric models always had less than five predictors and in most cases these predictors were derived from the ALS data (Tables A4 and A5).

Differences between Synthetic and Calibrated Predictions
All models included ALS acquisition random effects; however, the magnitude of these effects varied significantly between variables and modeling techniques. For semiparametric models, reductions in global and were negligible, which indicates that with these models using synthetic predictions and omitting the calibration step will not have consequences of practical importance in accuracy and precision. While semiparametric models seem to capture differences between ALS acquisitions directly, leaving little room for improvements to the calibration stage, parametric models operate differently, and their performance improves substantially with the calibration. Parametric models captured the main relationship with auxiliary information in their fixed-effect component, and a part of the variability between ALS acquisition was accounted for in the calibration stage ( Figure 4). Focusing on the parametric models, for AGB, CBD, CBH, and CFL calibration resulted in consistent reductions of the global and . For DWB and CH, three parametric models failed at reducing the global . We further inspected this issue by looking at ALS-acquisition-specific values of for the parametric models. For all variables and source of ALS metrics the |BI AS| decreased in the majority of the ALS acquisitions after the calibration was performed ( Figure 5). These reductions in the ALS acquisition BI AS also resulted in a reduction of RMSE. Similar patterns were observed for the e×tended dataset using only rasterized metrics ( Figure A3). The effect of the calibration in DWB was very small and erratic, which, in conjunction with the poor performance of the models for this variable, indicate that the variability in DWB cannot be predicted well by the auxiliary information used in this study or by geographic factors such as membership in a particular ALS acquisition area. For the remaining variables the effect of the calibration varied in magnitude depending on the variable and source of predictors. Most important effects of calibration were observed for CH followed by CBH and AGB, and the smallest reductions of ALS-specific BI AS were observed for CBD. Overall, the calibration step had positive impact in reducing the ALS acquisition biases and that is the main advantage of the mixed-effects parametric models developed in this study.
Remote Sens. 2020, 12, × FOR PEER REVIEW 16 of 37 Figure 5. Absolute value of bias, | |, and root mean square error, , for synthetic and calibrated predictions of parametric models by ALS data acquisitions. AGB, DWB, CBD, CH, CBH, and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load.

Differences between Parametric and Semiparametric Models
Parametric models provided lower values of the global RBI AS cal in 11 of the 12 possible combinations of response variables and sources of ALS metrics, while semiparametric models provided lower values of the global RRMSE cal in 8 of the 12 possible comparisons. With the exception of DWB, for which parametric models provided consistently better results in terms of RBI AS cal and RRMSE cal , there was no modeling technique that can be considered as clearly superior than the other. For DWB, CBD, CH, CBH, and CFL, the RBI AS cal of parametric models was smaller than the bias of their semiparametric counterparts regardless of the source of ALS predictors (Table 4). For AGB there was no modeling technique better than the other in terms of RBI AS cal . Semiparametric models provided better results in terms of global RRMSE cal for AGB, CBD, and CH, while for DWB the parametric models had consistently smaller values of RRMSE cal . The magnitude of the differences in RBI AS cal and RRMSE cal between parametric and semiparametric models were of small magnitude for those responses with larger R 2 . For CH and AGB, differences between modeling techniques in RBI AS cal and in RRMSE cal were below 1.35% and 2.27% respectively, regardless of the source of ALS metrics. Excluding DWB, magnitude of the differences in RBI AS cal were most important for the models for CBH and CFL using rasterized predictors, while differences in RRMSE cal were below 4% for all variables (Table 4). Similar results were obtained when comparing parametric and semiparametric models with an extended dataset that included 12 additional ALS acquisitions (Table A3). The smaller RBI AS cal of the parametric models and smaller RRMSE cal in semiparametric models clearly indicates performance tradeoffs between modeling techniques; however, the magnitude of these tradeoffs is for most variables of very small magnitude. Table 4. Differences between parametric and semiparametric models. AGB, DWB, CBD, CH, CBH, and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load. Superscripts p and sp indicate error metrics for parametric and semiparametric models respectively and the superscript m denote changes in RBI AS and RRMSE due to using a different modeling technique.  (Figure 6). In terms of RMSE cal , semiparametric models provided slightly lower acquisition specific values of RMSE cal for CBD and CH and for the remaining variables no modeling technique consistently provided lower values of this performance metric. For both sources of ALS predictors, the magnitude of the differences between parametric and semiparametric models in RMSE cal was, in general, orders of magnitude smaller than the values of RMSE cal themselves. Similar results were obtained when comparing parametric and semiparametric models with the extended dataset that included 12 additional ALS acquisitions ( Figure A4) and only used rasterized predictors. These results show that multiple exceptions to the trends observed for the global BI AS cal and RMSE cal can be observed when focusing on specific ALS acquisitions and confirmed that differences in performance due to the modeling technique are very minor.

t-Test
Remote Sens. 2020, 12, × FOR PEER REVIEW 18 of 37 smaller than the values of themselves. Similar results were obtained when comparing parametric and semiparametric models with the extended dataset that included 12 additional ALS acquisitions ( Figure A4) and only used rasterized predictors. These results show that multiple exceptions to the trends observed for the global and can be observed when focusing on specific ALS acquisitions and confirmed that differences in performance due to the modeling technique are very minor. Figure 6. ALS-acquisition-specific values of absolute bias, BI AS cal , and root mean squared error, RMSE cal , for calibrated predictions for different modeling techniques by ALS data acquisitions. AGB, DWB, CBD, CH, CBH, and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load.

Differences between Sources of ALS Metrics
For the parametric models, differences in global RBI AS cal between models based on point-cloud and rasterized ALS predictors were significant only for CBH, and the magnitude of these differences only reached 1.43%. For the semiparametric models, using point-cloud predictors resulted in significantly more accurate models for AGB, CH, CBH, and CFL, but the magnitude of these differences was below 2% for AGB and CH and below 4% and 6% for CBH and CFL respectively. Except for the semiparametric models for DWB, models using point-cloud predictors had lower values of global RMSE cal . However, the magnitude of these differences only reached a 5% for the semiparametric models for CBH, and for the remaining variables differences were always below 3.3% (Table 5). Table 5. Differences between models developed using point-cloud metrics and rasterized predictors. AGB, DWB, CBD, CH, CBH, and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load. Superscripts pc and r indicate error metrics for models point-cloud and rasterized predictors and the superscript s denote changes in RBI AS and RRMSE due to using a different source of ALS metrics. Finally, side by side comparison of ALS-acquisition-specific values of BI AS cal and RMSE cal for models using point-cloud and rasterized metrics showed that models using point-cloud predictors were not always better than the models using rasterized metrics ( Figure 7). The only exceptions to this trend were the models for CH in which using point-cloud predictors resulted in RMSE cal about 1 m smaller for most acquisitions. More important, the magnitude of the change in acquisition specific BI AS cal and RMSE cal was generally negligible compared to the values of these metrics. This suggests that using rasterized metrics instead of point-cloud metrics has a limited impact on accuracy.

Discussion
We considered the prediction of a wide array of attributes of interest for forest managers and observed important differences between variables. Models for CH and AGB showed a high predictive performance, similar to that found in previous studies [20,26,62,63]. The percentages of explained variance for both parametric and semiparametric models for CBH and CFL were similar to those obtained by [20] for a single ALS acquisition. For CBD the percentage of explained variance of the semiparametric models was about 10% smaller than those reported by [20] for ALS acquisition-specific models. The low R 2 values obtained for DWB were also about 10% smaller than those reported by [20] for local models, indicating that regional models developed using the auxiliary information considered in this study only provide a very crude approximation to reality for this variable.
Important efforts have recently been made to analyze the potential transferability of ALS models to new acquisitions not considered in the training dataset [27][28][29] or to a new collection over the same study area [64]. These studies showed that the degradation in model performance due to temporal offsets, and differences in the configuration of the ALS data collection are small compared to the losses in performance due to (1) ecological differences between the regions where the models were fit and transferred, and (2) the effect of the modeling technique. In our study we focused on developing models that can be calibrated to new acquisitions as a way of developing maps of forest fuels and structural attributes in a rapid manner. Our results partially align with those obtained by [29] and for almost all response variables, synthetic predictions from semiparametric models had lower RRMSE than synthetic predictions from parametric models. The most important effect of the calibration in the parametric models was the reduction in the ALS-acquisition-specific BI AS. Such BI AS reduction cannot be obtained transferring fixed-effect models to new ALS acquisition as fixed-effect models can only generate synthetic predictions. Once the parametric models were calibrated, few differences were observed between both modeling techniques. One factor that should be considered in future applications is that models developed with more ALS acquisitions will result in more reliable estimates of the variance parameters of the random effects [57]. For cases with fewer ALS acquisitions of large size or where the calibration is not an issue because no new ALS acquisitions are expected, other modeling techniques using only fixed-effects models would be more appropriate. We developed this study using eight ALS acquisitions, but main results were confirmed in the analysis of the extended dataset with 20 ALS acquisitions (Appendix A).
The minor effect that calibration had in the semiparametric models can be explained by the ability of random forest to effectively use large numbers of predictors and model nonlinear patterns [37,38]. If the training dataset covers the main ecological gradients of the area under analysis, random forest captures the most important sources of variability that can be explained by the auxiliary information leaving little room for improvements in the calibration step. The linear mixed-effects models with random intercepts and slopes developed in this study operate in a very different manner. The fixed component of the model characterizes the main relationships between auxiliary information and response variables and complex effects, specific of a given ALS data collection, are accounted for in the calibration step.
Results from our study support the idea that both techniques produce comparable results if calibration is a possibility and suggest the following recommendations for future applications. When calibration is a possibility, linear mixed-effects models do not provide consistently worse results than semiparametric models and have clear advantages in terms of simplicity, interpretability, and ease of use. Apart from the corresponding ground data and auxiliary information, only parameters reported in Tables A4 and A5 are necessary to develop the calibrations of the parametric linear models. For random forest-based semiparametric models, external users need to have access to the original random forest data-structure developed by the modeler to obtain predictions in new areas. In addition, linear mixed-effects models allow mapping uncertainty of predictions using well-known methods [58,65], but reporting the same uncertainties for semiparametric models is still a field under development [38,66]. When calibration is not a possibility, or a fast and direct model transfer is necessary, semiparametric models are clearly a better option. For these models, synthetic predictions and calibrated predictions showed little differences for all variables, which indicates that omitting the calibration step should not imply significant losses in accuracy with respect to a scenario where calibration is possible.
Regardless of the modelling technique and source of ALS predictors, for all responses analyzed in this study, ALS-acquisition-specific biases did not completely disappear after the calibration step indicating that global models should not be a replacement for ALS-acquisition-specific models, but instead a solution to enable a fast mapping of forest attributes. This result is consistent with [30,31] for nonparametric random forest models. The calibration approach used in this study relies on the assumption of having different error structures due to membership in a particular ALS acquisition. This assumption is reasonable because ALS acquisitions usually cover relatively homogeneous areas and the specifications in the ALS data collection can differ among ALS data-collection projects. The inclusion of ALS-acquisition-specific random effects allows accounting for these differences between ALS acquisitions in the calibration step. Nevertheless, other approaches to calibration should be investigated in the future. Of special interest is the approach proposed by [31] that develops calibrated predictions by considering the variability of geographic areas that do not have to match a specific ALS acquisition.
Differences between models using point-cloud and rasterized metrics were negligible (Figures 7 and A7). Differences of such a small magnitude were not expected; however, this result can be explained by: (1) the presence of GPS positioning errors that even if expected to be in the range of a meter introduce some noise in the ALS metrics [67], and (2) by the averaging that occurs when aggregating point-cloud metrics at the FIA plot level. FIA plots are composed of four 17.95 m radius subplots, with one subplot at the center and the other three subplots 36.58 m from the central plot. Point clouds for the subplots are combined when computing ALS predictors for an FIA plot, resulting in an aggregation effect similar to the one that occurs when computing rasterized metrics. While the use of pre-rasterized layers to e × tract ALS predictors implies a departure from traditional workflows using ALS data, this process is equivalent to that used in many applications using optical sensors such as Landsat images (e.g., [26]) and has been used in previous studies to combine ALS, optical, and topographic predictors [68]. Additionally, and more importantly: (1) ALS metrics are computed in a standardized and structured way that for every point of the territory assigns a unique support area [69] for the computation of ALS metrics, and (2) empirical results from our analysis showed that, by using rasterized ALS metrics instead of point-cloud metrics, the predictive performance of the models does not worsen substantially and allows generating cartographic products useful for multiple planning tasks, with accuracies and precisions similar to those obtained using point-cloud predictors.
The small differences between models developed using point-clouds or rasterized predictors in combination with the fact that sharing rasterized products is far more operationally feasible than sharing large point-clouds, opens the possibility of developing models in a more efficient way. Future mapping applications that require rapid delivery could be developed extracting ALS metrics for the training units using interpolation over pre-rasterized products available from spatial data infrastructures serving raster layers with metrics derived using a standardized workflow. Operating in this way eliminates the need of intersecting ground measurements with large point clouds that are more difficult to share and manipulate, since a larger number of users will be more comfortable extracting auxiliary variables from raster files than clipping and normalizing point-cloud data. These two factors combined can reduce (1) the time needed to develop models and (2) the time that forest managers need to wait for maps of structural attributes once an ALS data are acquired and processed.

Conclusions
The increasing availability of ALS data in the western US states raises questions about methodological aspects to derive cartographic products of forest attributes such as AGB, DWB, CBD, CH, CBH, and CFL. This study analyzes three factors of importance for development of a regional strategy to map forest structural attributes that allows using available ALS data acquisitions and incorporating new ALS project areas in a rapid manner. The main conclusion regarding these factors are: Transferability and calibration. For both modeling techniques, calibration reduced bias problems when predicting to ALS data acquisitions not included in the training dataset; however, the effect of calibration was much more important for parametric linear mixedeffects models than for semiparametric models using random forest. This indicates that linear models developed over such a large region only account for main trends with respect to the auxiliary information and need to (1) account for the variability between ALS data acquisitions and (2) be calibrated to local conditions to eliminate bias problems.
Modeling technique. Once the calibration is performed, both modeling techniques have similar performance in terms of acquisition-specific RBI AS or RRMSE. Interpretability of results and simplicity make the linear mixed-effects models more appealing for situations where performing a calibration for new acquisitions is feasible, and the small effect of calibration in semiparametric models makes them a better choice for cases where calibration is not possible.
Source of ALS metrics. Differences in performance between models developed with predictors obtained by interpolation over rasterized ALS metrics and models developed using predictors computed by clipping the point clouds coincident with areas where field measurements were taken were minimal for the variables considered in this study. This clearly suggests that raster layers can be used to drive the entire modeling workflow with little loss in terms of performance, which can be important for time-sensitive applications because the pre-processing of the ALS data can be substantially simplified. in Table A2. For each variable we obtained parametric and semiparametric models using the FIA plots in the combined dataset using the model selection methods described in Sections 2.3.1 and 2.4.1. Error metrics and model comparisons described in Section 2.5 were computed for the parametric and semiparametric models fitted with the extended dataset.
Selected parametric models are reported in Table A6, and predictors and variance components of the semiparametric models are shown in Figure A6 and Table A8. Models fitted with the extended dataset were very similar to those obtained with the set of eight ALS acquisitions used in the main study and few differences were observed with respect to: (1) global model performance ( Figure A2), (2) predictors included in the model, Table A6 and Figure A6, and (3) structure of model errors and random effects (Tables A6 and A8).
The effect of the calibration was important for the parametric models and negligible for the semiparametric models. The most important effect of calibration in the parametric models was a reduction in the ALS-acquisition-specific BI AS ( Figure A3). The magnitude of this effect varied among variables from a negligible effect for CBD to substantial reductions in ALS-acquisition-specific BI AS for AGB, CBH, and CFL.
Once calibrated, the parametric models had smaller values of global BI AS except for AGB and larger values of global RMSEs. These results are similar to those obtained with the main dataset and confirm that, globally, parametric models are less prone to bias problems, and semiparametric models tend to have smaller RMSEs ( Figure A4). Difference in BI AS and RMSEs for specific ALS acquisitions also showed that trends observed for the global BI AS and RMSEs have multiple exceptions when observing ALS acquisitions specific performance metrics. ALS-specific BI AS and RMSEs showed that no modeling technique was consistently better than the other one in terms of these metrics ( Figure A4). This result indicates that factors such as model interpretability or the possibility of omitting the calibration are expected to be more relevant than differences in model performance when selecting a modeling technique for a similar problem. Overall, results from this extended analysis confirmed the main conclusions obtained for the factors transferability and effect of calibration and differences between modeling techniques.
Remote Sens. 2020, 12, × FOR PEER REVIEW 24 of 37 2.3.1 and 2.4.1. Error metrics and model comparisons described in Section 2.5 were computed for the parametric and semiparametric models fitted with the extended dataset. Selected parametric models are reported in Table A6, and predictors and variance components of the semiparametric models are shown in Figure A6 and Table A8. Models fitted with the extended dataset were very similar to those obtained with the set of eight ALS acquisitions used in the main study and few differences were observed with respect to: (1) global model performance ( Figure A2), (2) predictors included in the model, Table  A6 and Figure A6, and (3) structure of model errors and random effects (Tables A6 and  A8).
The effect of the calibration was important for the parametric models and negligible for the semiparametric models. The most important effect of calibration in the parametric models was a reduction in the ALS-acquisition-specific ( Figure A3). The magnitude of this effect varied among variables from a negligible effect for CBD to substantial reductions in ALS-acquisition-specific for AGB, CBH, and CFL. Once calibrated, the parametric models had smaller values of global except for AGB and larger values of global s. These results are similar to those obtained with the main dataset and confirm that, globally, parametric models are less prone to bias problems, and semiparametric models tend to have smaller s ( Figure A4). Difference in and s for specific ALS acquisitions also showed that trends observed for the global and s have multiple exceptions when observing ALS acquisitions specific performance metrics. ALS-specific and s showed that no modeling technique was consistently better than the other one in terms of these metrics ( Figure A4). This result indicates that factors such as model interpretability or the possibility of omitting the calibration are expected to be more relevant than differences in model performance when selecting a modeling technique for a similar problem. Overall, results from this extended analysis confirmed the main conclusions obtained for the factors transferability and effect of calibration and differences between modeling techniques.        Figure A3. ALS-acquisition-specific and for synthetic and calibrated predictions for parametric and models. AGB, DWB, CBD, CH, CBH and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load. Figure A3. ALS-acquisition-specific BI AS and RMSE for synthetic and calibrated predictions for parametric and models. AGB, DWB, CBD, CH, CBH and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load.

Set of Acquisitions t-Test Accuracy
Remote Sens. 2020, 12, × FOR PEER REVIEW 28 of 37 Figure A4. ALS-acquisition-specific values of absolute relative bias, , and relative root mean squared error, , for synthetic and calibrated predictions for parametric and models. AGB, DWB, CBD, CH, CBH and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load. Figure A4. ALS-acquisition-specific values of absolute relative bias, RBI AS cal , and relative root mean squared error, RRMSE cal , for synthetic and calibrated predictions for parametric and models. AGB, DWB, CBD, CH, CBH and CFL respectively indicate aboveground biomass, downed woody biomass, canopy bulk density, canopy height, canopy base height, and canopy fuel load.    Figure A5. Predictors included in the semiparametric models and variable importance as percent reduction in mean squared error. R and PC indicate models developed with rasterized and point-cloud metrics respectively. Model response variables are: aboveground biomass (AGB), downed woody biomass (DWB), canopy bulk density (CBD), canopy height (CH), canopy base height (CBH), and canopy fuel loading (CFL).   igure A6. Predictors included in the semiparametric models and variable importance as percent reduction in mean uared error for the extended set of 20 ALS acquisitions using rasterized predictors described in Appendix A. Model sponse variables are: aboveground biomass (AGB), downed woody biomass (DWB), canopy bulk density (CBD), canopy eight (CH), canopy base height (CBH), and canopy fuel loading (CFL).

Appendix B. Selected Parametric and Semiparametric Models
able A8. Variance parameters for the random effects and model errors for the semiparametric models obtained for each sponse variable for the extended set of 20 ALS acquisitions using rasterized predictors described in Appendix A. Model sponse variables are: aboveground biomass (AGB), downed woody biomass (DWB), canopy bulk density (CBD), canopy eight (CH), canopy base height (CBH), and canopy fuel loading (CFL).  Figure A6. Predictors included in the semiparametric models and variable importance as percent reduction in mean squared error for the extended set of 20 ALS acquisitions using rasterized predictors described in Appendix A. Model response variables are: aboveground biomass (AGB), downed woody biomass (DWB), canopy bulk density (CBD), canopy height (CH), canopy base height (CBH), and canopy fuel loading (CFL). Table A8. Variance parameters for the random effects and model errors for the semiparametric models obtained for each response variable for the extended set of 20 ALS acquisitions using rasterized predictors described in Appendix A. Model response variables are: aboveground biomass (AGB), downed woody biomass (DWB), canopy bulk density (CBD), canopy height (CH), canopy base height (CBH), and canopy fuel loading (CFL).