Estimation of Changes of Forest Structural Attributes at Three Different Spatial Aggregation Levels in Northern California using Multitemporal LiDAR

Abstract

Accurate estimates of growth and structural changes are key for forest management tasks such as determination of optimal rotation times, optimal rotation times, site indices and for identifying areas experiencing difficulties to regenerate. Estimation of structural changes, especially for biomass, is also key to quantify greenhouse gas (GHG) emissions/sequestration. We compared two different modeling strategies to estimate changes in V, BA and B, at three different spatial aggregation levels using auxiliary information from two light detection and ranging (LiDAR) flights. The study area is Blacks Mountains Experimental Forest, a ponderosa pine dominated forest in Northern California for which two LiDAR acquisitions separated by six years were available. Analyzed strategies consisted of (1) directly modeling the observed changes as a function of the LiDAR auxiliary information ($\delta$-modeling method) and (2) modeling V, BA and B at two different points in time, including a term to account for the temporal correlation, and then computing the changes as the difference between the predicted values of V, BA and B for time two and time one. We analyzed predictions and measures of uncertainty at three different level of aggregation (i.e., pixels, stands or compartments and the entire study area). Results showed that changes were very weakly correlated with the LiDAR auxiliary information. Both modeling alternatives provided similar results with a better performance of the $\delta$-modeling for the entire study area; however, this method also showed some inconsistencies and seemed to be very prone to extrapolation problems. The $y$-modeling method, which seems to be less prone to extrapolation problems, allows obtaining more outputs that are flexible and can outperform the $\delta$-modeling method at the stand level. The weak correlation between changes in structural attributes and LiDAR auxiliary information indicates that pixel-level maps have very large uncertainties and estimation of change clearly requires some degree of spatial aggregation; additionally, in similar environments, it might be necessary to increase the time lapse between LiDAR acquisitions to obtain reliable estimates of change.

1. Introduction

Light detection and ranging LiDAR data have been extensively used in forest inventories to provide auxiliary information that is highly correlated with multiple forest structural attributes [1,2,3]. This strong correlation allows estimating forest structural attributes more efficiently than if only field measurements are available [4]. In addition, the spatially explicit nature of LiDAR enables the mapping of forest attributes at fine resolutions (e.g., [2,5]). Accurate estimates of growth and structural changes are key for forest management as multiple management tasks such as determination of optimal rotation times, calculation of site indexes or the identification of areas experiencing difficulties in regeneration. Estimation of biomass is also key to quantifying greenhouse gas (GHG) emissions/sequestration, to comply with the International Panel on Climate Change (IPCC) reporting and good practice guidelines [6], and to develop a correct appraisal of forest resources for carbon markets. The extensively used area based approach (ABA) [1] provides a way to estimate forest attributes at multiple levels ranging from single pixels to large areas using LiDAR auxiliary information [7]. Availability of repeated LiDAR data acquisitions has opened the door to estimation of changes in forest structural attributes over time (e.g., [8,9]) using the ABA method.

In the ABA, the area under study is covered by a regular grid that will define a population of pixels or grid cells. In this approach, the field plots used to train models and the grid cells are of the same size, typically between 400 m² and 900 m². A direct application of predictive models will render predictions for grid units of size too small to be considered of interest for reporting in forest inventories. Areas of interest (AOIs) (i.e., the areas for which estimates are needed) are typically geographic units that can vary in size depending on the particular application. For worldwide inventories or inventories over continents or countries, AOIs are typically administrative or political units such as countries or municipalities. In forest management applications, AOIs are typically stands, compartments or even complete forests or landscapes. All these AOIs require spatial aggregation of grid units. However, validations of predictive models in the ABA literature are typically performed using global metrics of model fit, such as the sample-based root mean square error or bias, that provide average measures of uncertainty for predictions made for pixels or plots. These measures of uncertainty derived from the model fitting stage do not directly translate into measures of uncertainty for predictions for AOIs composed of multiple pixels (i.e., countries, municipalities, forests, stands, etc.). In addition, even when considering single pixels, they are not AOI-specific, as they only provide an average value, across the entire population, of the error that can be expected using a given model.

Thus, it is clear that uncertainty measures used as quality controls in forest inventories need to be made at the AOI-level and change estimation is not an exception. For large areas holding large sample sizes, AOI-specific estimates of means or totals and their measures of uncertainty can be obtained using direct estimators (e.g., [10,11,12]) that use only use sample data from the AOI under consideration. However, if the AOI sample sizes are not large enough to support direct estimates with reliable precision, then they must be regarded as small areas [13].

Small area estimation (SAE) techniques, especially empirical best linear unbiased predictors (EBLUPs) in combination with the ABA approach have been used to obtain estimates, and their corresponding measures of uncertainty for subpopulations such as municipalities [14], groups or management units [15] and stands [4,14,16,17]. SAE techniques allow correcting the potential bias problems of synthetic predictions (i.e., predictions developed assuming that a general model developed at the population level holds for all subpopulations) and also permit reducing the large variance problems of direct estimators when AOIs sample sizes are small [13]. In addition, while EBLUPs have been extensively used in SAE contexts, they can also be used to produce estimates for subpopulations or AOIs with large sample sizes and preserve important advantages over other methods. First, they allow obtaining model-unbiased estimates and their corresponding measures of uncertainty for all AOIs using a single model that explicitly considers potential variations between AOIs. This is a clear advantage over synthetic methods that assume that a certain relation derived for the entire population holds in all AOIs. A second advantage of EBLUPs is that it is possible to reduce the modeling effort required by direct model-based or model-assisted methods where a model is needed for each AOI. It is thus clear that SAE techniques in combination with LiDAR auxiliary information have potential applications in multiple forest inventories contexts. Unfortunately, to the best of our knowledge, all studies on SAE and forest inventories have focused on estimation of structural attributes at a given point in time, and little is known about: (1) their performance when applied to forest structure change estimation, and (2) about how these techniques compare to other methods used for estimation of changes in AOIs comprising entire populations [10,12,18,19] and especially subpopulations [20].

In this study, we analyzed the two most commonly used strategies to model changes in structural forest attributes using repeated LiDAR acquisitions, and analyzed their performance when used to obtain EBLUPs for AOIs of different size. The first strategy, referred hereafter as the $\delta$-modeling method, considers the change, $\delta$, over the time between LiDAR acquisitions as the model response. The second strategy, which we will call $y$-modeling method, focuses on modeling the structural attributes $y$, and their derived change over time. As a novelty, in the $y$-modeling method, the temporal correlation of both model errors and AOI random effects were taken into account. We considered changes in three structural variables, and AOIs at three different spatial aggregation levels in order to provide insights for future applications where estimates for an entire population and for subpopulations of different sizes are needed. Variables under study are standing volume (V), above ground biomass (B) and basal area (BA) and AOIs subject to analysis are (1) an entire forested area or landscape, (2) subpopulations that in this case are forest stands and (3) pixels as gridded maps are common output in mapping applications.

2. Materials and Methods

2.1. Study Area

The study area is Blacks Mountains Experimental Forest (BMEF), a 3715 ha forest managed by the United States Forest Service, located northeast of Lassen National Park in northern California, USA (Figure 1). Elevation ranges from 1700 m to 2100 m above sea level. Slopes are gentle (<10%) on the lower parts of the forest and moderate (10%–40%) at higher elevations. Climate is Mediterranean with a certain degree of continentality, with dry summers and wet and cold winters when precipitation is in the form of snow. Average precipitation is 460 mm per year with monthly average temperatures that range from −9 °C to 29 °C. Soils are developed over basalts with depths that range from 1 to 3 m. Ponderosa pine (Pinus ponderosa Lawson & C. Lawson) dominated forest occupies the majority of the area. Incense cedar (Calocedrus decurrens (Torr.) Florin), white fir (Abies concolor (Gordon & Glend) Hildebr) and Jeffrey Pine (Pinus jeffreyi Grev & Balf.) are abundant accompanying species. Forest structure is relatively open and the canopy cover varies greatly within the forest (see Figure 1). A more detailed description of the study area can be found in [21,22].

2.2. Sampling Design and Field Data

In total, 106 forested stands were delineated in BMEF. Small non-forested patches were masked in the study area and hence were not considered part of the population under study. Out of the 106 forested stands, 24 were selected and sampled in the field. Nine of the remaining 82 unsampled stands were subject to thinning during the period between two available LiDAR acquisitions (i.e., 2009–2015) and all thinning operations were finished by fall 2011. These nine stands are located on the southwestern edge of BMEF and were analyzed separately because the sample of field plots used to train the LiDAR models did not include any stand subject to similar silvicultural interventions (Figure 1). Sampled stands come from a long-term research project initiated in BMEF in 1991 and, excluding the nine thinned stands, were representative of the forest structures and forest management treatments applied in rest of BMEF.

Sampled stands were subject to six different types of treatments resulting from crossing two different factors. The first factor is the structural diversity. It has three levels referred hereafter as low structural diversity (LoD), high structural diversity (HiD) and research natural areas, RNA, or controls. Low structural diversity stands are subject to thinning operations aiming to generate simplified single-strata structures. High diversity stands are subject to thinning where all canopy layers and age groups are preserved, resulting in a multi-storied forest structure with trees of different sizes and ages. Neither the HiD stands nor the LoD stands were subject to thinnings during the period between the two available LiDAR flights. Finally, RNA stands are not subject to any thinning or harvest operation. In total, 10 LoD, 12 HiD and two RNA stands were measured in the field. The second factor under consideration was the presence or absence of prescribed forest fires. Half of the LoD, HiD and RNA stands sampled in the field had been subject to prescribed fires, but only one of the RNA stands was subject to prescribed fires during the period 2009–2015.

A sample of 151, 16 m radius plots (804 m²) were measured in the field during the summer of 2009 and then remeasured during the summer of 2016. All field plots were located on nodes of the 100 m by 100 m grid of monumented markers at BMEF. Coordinates of the makers were determined using traverse methods and survey grade GPS observations and have an accuracy of 15 cm or better (see [23]). For each of the 26 stands selected for sampling, a node of the BMEF 100 m grid was randomly selected and used as a starting node for a 282 m by 282 m grid formed by selecting every other plot of the 100 m grid moving in the diagonal directions. Field measurements were taken on the nodes of the 282 m by 282 m grid (see Figure 1).

Within each field plot all live trees with DBH larger than 9 cm, and all dead standing trees with DBH larger than 12 cm, were stem mapped and measured for DBH and height. Plot basal area (BA) was derived directly from the field measurements. Volume (V) and above ground biomass (B) were computed as the sum of the individual tree volumes and biomasses of all standing trees. Individual tree volumes and biomasses were estimated using species-specific allometric models included in the national volume estimation library (NVEL) and in the national biomass estimation library (NBEL). To account for the one-year difference between acquisition of field measurements in 2016, and the second LiDAR data acquisition obtained in 2015; plot-level values of the variables under analysis were computed for 2015 by linearly interpolating between the values obtained for 2009 and 2016. Finally, for each field plot we computed the change of V, B and BA on a per year basis, as the difference of the plot-level values in 2009 and 2015 divided by 6. For two plots close to the southeastern boundary of the forest, changes in V were extremely large, more than three standard deviations away from the mean value for the change in volume. These anomalous plots were removed from the analysis because such large changes seemed to be derived from edge effects. Plot-level values for 2009, 2016 and per year increments for the period, 2009-2015, for V, B and BA, in the remaining 149 plots are summarized in

Table 1. Minimum (Min), mean (Mean), standard deviation (Sd), and maximum (Max) of the plot-level values for 2009, 2015 and yearly increments for the period 2009–2015. Values of volume V, basal area BA and biomass B are expressed on a per-hectare basis.

Variable (Units)	Period	Min	Mean	Sd	Max
V(m3 ha−1)	2009	19.87	166.93	119.66	619.43
BA(m2 ha−1)		3.81	23.43	12.02	66.54
B(Mg ha−1)		8.31	83.65	61.55	323.30
V(m3 ha−1)	2015	17.20	175.52	117.04	644.30
BA(m2 ha−1)		3.42	25.45	12.01	67.47
B(Mg ha−1)		8.34	89.38	60.29	335.03
V(m3 ha−1year−1)	Increment 2009–2015	−10.89	1.43	3.88	11.19
BA(m2 ha−1year−1)		−0.91	0.34	0.45	1.74
B(Mg ha−1year−1)		−5.81	0.95	1.97	5.99

.

For the nine unsampled stands thinned during the period 2009–2015 all thinning operations were completed by fall 2011. In total 427.40 hectares were thinned with prescriptions that varied among stands. Approximately 80% of the area was thinned from below, leaving a residual basal area of 17.22 m² ha⁻¹ to 25.25 m² ha⁻¹. For the remaining 20% of the area, approximately one quarter was not thinned while the other three quarters were thinned to a residual BA that ranged from 6.89 m² ha⁻¹ to 13.77 m² ha⁻¹. Fresh weight of total extractions for the 427.40 hectares subject to thinning was 11,009.38 Mg of logs and 23,164.32 Mg of chipped material.

2.3. LiDAR Data Acquisitions

Two LiDAR acquisitions are available for BMEF. The first LiDAR dataset was acquired during the summer of 2009 using a Leica ALS 50 discrete return sensor. Flying altitude was 900 m, side-lap between adjacent flight lines was at least 50% and scanning angle was ±14°. The LiDAR data vendor generated digital terrain models (DTMs) with an accuracy of 15 cm at 95% confidence level. Additional details on the LiDAR data collection for the 2009 acquisition can be found in [24]. The same vendor in the study area performed a second LiDAR acquisition during the summer of 2015 using the same sensor, flying altitude and side-lap specifications. DTMs were also created for 2015 by the vendor.

Four sets of auxiliary variables were considered in this study. The first two sets are composed of 42 LiDAR predictors computed for each acquisition date. Set 1 will represent the predictors for 2009 and Set 2 the predictors for 2015. These predictors are descriptors of the point cloud height distributions and were all relative quantities to avoid introducing noise due to local differences in the point cloud densities of 2009 and 2015 [25]. The third set of predictors, Set 3, was computed as the differences between the 2009 and the 2015 LiDAR predictors. Finally, the fourth set of predictors, Set 4, included the incoming solar radiation computed using the Environmental Systems Research Institute (ESRI) ArcGIS Area Solar Radiation tool [26] with the 2009 digital surface model (DSM) as input; and two treatments: (1) single- or multi-story structural diversity and (2) presence or absence of prescribed fires. All predictors were computed for each field plot and for a grid with a cell size of 805 m² covering the entire BMEF. The cell size matched the field plot size and each cell of the grid was considered a population unit, equivalent to the field plots. Predictors and their corresponding acronyms used in further sections are summarized in

Table A1. Sets of candidate predictors used in the study. Predictors included in the models to predict structural changes are highlighted with a boldface font. HiD, LoD and RNA represent the high diversity, low diversity and research natural areas respectively.

Description Auxiliary Variables Sets 1, 2 and 3	Acronym			Description Auxiliary Variables Set 4	Acronym
	Set 1 Year: 2009	Set 2 Year: 2015	Set 3, Difference 2015-2009		Set 4
Minimum, maximum, mean, mode, standard deviation, variance, coefficient of variation and interquartile range of the distribution of heights of the point cloud.	Elev_min09	Elev_min15	$\delta$Elev_min15-09	Incoming solar radiation	Solar_radiation
	Elev_max09	Elev_max15	$\delta$Elev_max15-09
	Elev_mean09	Elev_mean15	$\delta$Elev_mean15-09	Structural diversity, factor with three levels HiD, LoD and RNA. Coded using two dummy variables. RNA reference level.	HiD
	Elev_mean209	Elev_mean215	$\delta$Elev_mean215-09
	Elev_mode09	Elev_mode15	$\delta$Elev_mode15-09
	Elev_stddv09	Elev_stddv15	$\delta$Elev_stddv15-09		LoD
	Elev_var09	Elev_var15	$\mathit{\delta }$Elev_var15-09
	Elev_CV09	Elev_CV15	$\delta$Elev_CV15-09	Presence absence of prescribed fires. Coded using a dummy variable taking value 1 for stands where prescribed fires are applied and 0 otherwise.	Burned
	Elev_IQ09	Elev_IQ15	$\delta$Elev_IQ15-09
	Elev_AAD09	Elev_AAD15	$\delta$Elev_AAD15-09
	Elev_MADmed09	Elev_MADmed15	$\delta$Elev_MADmed15-09
	Elev_MADmod09	Elev_MADmod15	$\delta$Elev_MADmod15-09
Percentiles of the distribution of heights of the point cloud.	Elev_P0109	Elev_P0115	$\delta$Elev_P0115-09
	Elev_P0509	Elev_P0515	$\delta$Elev_P0515-09
	Elev_P1009	Elev_P1015	$\delta$Elev_P1015-09
	Elev_P2009	Elev_P2015	$\delta$Elev_P2015-09
	Elev_P3009	Elev_P3015	$\delta$Elev_P3015-09
	Elev_P4009	Elev_P4015	$\delta$Elev_P4015-09
	Elev_P5009	Elev_P5015	$\mathit{\delta }$Elev_P5015-09
	Elev_P6009	Elev_P6015	$\delta$Elev_P6015-09
	Elev_P7009	Elev_P7015	$\delta$Elev_P7015-09
	Elev_P7509	Elev_P7515	$\delta$Elev_P7515-09
	Elev_P8009	Elev_P8015	$\delta$Elev_P8015-09
	Elev_P9009	Elev_P9015	$\delta$Elev_P9015-09
	Elev_P9509	Elev_P9515	$\delta$Elev_P9515-09
	Elev_P9909	Elev_P9915	$\delta$Elev_P9915-09
Canopy relief ratio	CRR09	CRR15	$\mathit{\delta }$CRR15-09
Percentage of first (Fst) and all (All) returns above 2 m	PcFstAbv209	PcFstAbv215	$\mathit{\delta }$PcFstAbv215-09
	PcAllAbv209	PcAllAbv215	$\mathit{\delta }$PcAllAbv215-09
Ratio all returns above 2 m to first returns	AllAbv2Fst09	AllAbv2Fst15	$\delta$AllAbv2Fst15-09
Percentage of first returns above the mean and mode	PcFstAbvMean09	PcFstAbvMean15	$\delta$PcFstAbvMean15-09
	PcFstAbvMode09	PcFstAbvMode15	$\delta$PcFstAbvMode15-09
Percentage of all returns above the mean and mode	PcAllAbvMean09	PcAllAbvMean15	$\delta$PcAllAbvMean15-09
	PcAllAbvMode09	PcAllAbvMode15	$\delta$PcAllAbvMode15-09
Ratio of all returns above the mean and mode to number of first returns	AllAbvMeanFst09	AllAbvMeanFst15	$\delta$AllAbvMeanFst15-09
	AllAbvModeFst09	AllAbvModeFst15	$\delta$AllAbvModeFst15-09
Proportion of points in the height intervals [0,0.5), [0.5,1), [1,2), [2,4), [4,8) and [8,16) meters.	Prop0_0509	Prop0_0515	$\delta$Prop0_0515-09
	Prop05_109	Prop05_115	$\delta$Prop05_115-09
	Prop1_209	Prop1_215	$\delta$Prop1_215-09
	Prop2_409	Prop2_415	$\delta$Prop2_415-09
	Prop4_809	Prop4_815	$\delta$Prop4_815-09
	Prop8_1609	Prop8_1615	$\delta$Prop8_1615-09

.

2.4. AOIs, Target Parameter and Overview of Modelling Strategies

Two different types of subsets of population units will be repeatedly used throughout the manuscript in remaining sections. These subsets and their corresponding notation are: the sample of plots measured in the field, denoted using sub-index $s$ and the target AOIs represented by pixels, denoted using sub-index $\alpha$.

Three different groups of AOIs representing different levels of spatial aggregation were analyzed. The first group represents the largest level of spatial aggregation and represents the entire population under study. Within this group, we considered the set of all sampled stands, SS, and the entire BMEF study area after removing the nine thinned and unsampled stands, SA (i.e., sampled and unsampled but not thinned stands). The second group consists of the 106 forested stands in BMEF. In this group, we considered separately the unsampled and thinned stands (nine stands), unsampled and not thinned stands (73 stands) and sampled and not thinned stands (24 stands). Finally, the third group is the set of all pixels of the LiDAR grid covering the forested area in BMEF.

The main objective of this study was to analyze AOI-specific estimates of the change between 2009 and 2015 for three different structural variables, (V, B and BA). We will use the generic term variable of interest and the letter $y$ to refer to the forest structural variables, and the term target parameter and Greek letter $\Delta$ to refer to the quantities that we seek to estimate. Hereafter, target parameters will always refer to changes over time for the totals of the variables of interest in the considered AOIs, and will be expressed in a per hectare and year basis.

Considering that all pixels have the same area, the target parameter ${\Delta }_{\alpha }$ for a generic AOI or subset of population units, $\alpha$, can be expressed as:

$${\Delta }_{\alpha }={\displaystyle \sum }_{i=1}^{N_{\alpha }}{K}_{y\alpha }(y_{i\alpha 15}-y_{i\alpha 09})={\displaystyle \sum }_{i=1}^{N_{\alpha }}{K}_{\delta \alpha }{\delta }_{i\alpha },$$

(1)

where $N_{\alpha }$ is the number of population units (i.e., pixels) in the AOI. The terms $y_{i\alpha 15}$ and $y_{i\alpha 09}$ respectively represent the value of the variable of interest for 2009 and 2015 for the $i^{th}$ population unit of $\alpha$, and ${\mathsf{\delta }}_{i\alpha }$ is the change, for the $i^{th}$ pixel of $\alpha$, in the variable of interest during the period 2009 to 2015. Finally, for comparability with previous studies, the variables of interest will be expressed in a per unit area basis, and the increments ${\mathsf{\delta }}_{i\alpha }$ will be expressed in a per unit area and year basis. Thus, to ensure that ${\Delta }_{\alpha }$ is expressed in the correct units, it is necessary to introduce the factors $K_{y\alpha }$ and $K_{\delta \alpha }$. When $y_{i\alpha 15}$ and $y_{i\alpha 09}$ are expressed in a per unit area basis $K_{y\alpha }=\frac{1}{6N_{\alpha }}$ and for ${\mathsf{\delta }}_{i\alpha }$ in a per unit area and year basis $K_{\delta \alpha }=\frac{1}{N_{\alpha }}$.

We calculated AOI estimates using two different methods. The first, $\delta$-modeling method, for estimation of change uses models similar to those in approach A5 of Poudel et al. [8]. In this approach, the change in a structural variable at the plot/pixel-level (${\delta }_{i\alpha }$) is directly modeled as a function of the LiDAR auxiliary variables available for the study area. The second, the $y$-modeling method, uses a modified version of approach A4 of Poudel et al. [8] to obtain AOI-specific estimates of change. Models in this approach jointly relate structural variables ($y_{\alpha 15}$ and $y_{\alpha 09}$) and LiDAR auxiliary information at a given point in time and account for the correlation between errors obtained for the same plot/pixel at different times. For both methods, variability between stands was accounted by considering them as small areas. Thus, stand-level random effects were included in the models.

2.5. $\delta$-Modeling Method

2.5.1. Model $\delta$-modeling Method

Models in the $\delta$-modeling method relate the change (per year) of the variable of interest in a population unit to the auxiliary variables for the population unit. To indicate that these models consider change in the variables of interest directly, model parameters, stand-level random effects and model errors will include the subscript $\delta$. Three different types of auxiliary variables were considered as potential predictors in the $\delta$-modeling method. First, changes in the LiDAR auxiliary variables for the period 2009-2015, Set 3, were considered following Poudel et al. [8] as changes in LiDAR predictors are expected to correlate with growth or changes in forest attributes. Forest structure relates to growth. Thus, the LiDAR auxiliary variables for 2009, Set 1, were also considered as potential predictors that act as proxies for forest structure at the beginning of the period 2009-2015. Finally, the incoming solar radiation and the structural diversity factors and presence of prescribed fires, Set 4, were also considered as potential predictors.

For the $j^{th}$ population unit in the $i^{th}$ stand, models of the $\delta$-modeling method have the form:

$${\delta }_{ij}={\mathit{x}}_{\delta ij}^t{\mathit{\beta }}_{\delta }+v_{\delta i}+{\epsilon }_{\delta ij},$$

(2)

where $t$ indicate the transpose operator and ${\mathit{x}}_{\delta ij}^t$ is a vector of auxiliary variables in which the first element takes the value 1 for the intercept. The term ${\mathit{\beta }}_{\delta }$ is a vector of model coefficients where the first element is the intercept of model (2). Selection of auxiliary variables included in the model was performed using the method described in [27]. Stand-level random effects $v_{\delta i}$ are assumed to be independently and identically distributed (i.i.d.) normal random variables $v_{\delta i}\sim N(0,{\sigma }_{\delta v}^2)$ for all $i=1,\dots ,D$, where $D$ is the total number of stands in the study area. Model errors are i.i.d. normal random variables ${\epsilon }_{\delta ij}\sim N(0,{\sigma }_{\delta \epsilon }^2)$ independent of the stand-level random effects (i.e., $Cov({\epsilon }_{\delta ij}{v}_{\delta ,k})=0$, for all $i,j$ and $k$). Models with spatially correlated errors and with non-constant error variances were initially considered but discarded in the model selection stage, as they were not found to be significant (see Section 2.5.3).

For a generic set of population units denoted by subscript $\xi$ (which can represent either $s$ or $\alpha$), the relation in matrix notation between vector of changes of structural variables ${\mathsf{\delta }}_{y\xi }$, and the auxiliary variables included in the model (${\mathit{X}}_{\delta \xi }$), is expressed as:

$${\mathsf{\delta }}_{\xi }={\mathit{X}}_{\delta \xi }{\mathit{\beta }}_{\delta }+{\mathit{Z}}_{\delta \xi }{\mathit{v}}_{\delta }+{\mathit{\epsilon }}_{\delta \xi },$$

(3)

where ${\mathsf{\delta }}_{\xi }={({\delta }_1,\dots ,{\delta }_{N_{\xi }})}^t$, with ${\delta }_k$ being the yearly change for the forest structural variable $y$, in the $k^{th}$ unit of $\xi$, and $N_{\xi }$ is the number of elements in the set $\xi$. The $k^{th}$ element of $\xi$ will be an element of a given stand. To explicitly indicate this membership we will use, when necessary, the sub-indexes $i^{th}$ and $j^{th}$ to respectively indicate the stand and index of the element within the stand. The $k^{th}$ row of the matrix ${\mathit{X}}_{\delta \xi }$ is ${\mathit{x}}_{\delta k}^t$. The vector ${\mathit{v}}_{\delta }={(v_{\delta 1},\dots ,v_{\delta D})}^t$ is a vector of stand-level random effects with variance covariance matrix ${\mathit{G}}_{\delta }={\sigma }_{\delta v}^2{\mathit{I}}_D$, where ${\mathit{I}}_{\mathit{D}}$ is the identity matrix of dimension $D$. The matrix ${\mathit{Z}}_{\xi \delta }$ is a $N_{\xi }$ x $D$ incidence matrix that describes stand membership for each population unit. The $r^{th}$ row of ${\mathit{Z}}_{\xi \delta }$ have zeros at all positions except at position $i$, where $i$ is the index of the stand to which the $k^{th}$ unit of $\xi$ belongs. Finally, ${\mathit{\epsilon }}_{\xi \delta }$ is a vector of $N_{\xi }$ model errors with diagonal variance covariance matrix ${\mathit{R}}_{\xi \delta }={\sigma }_{\delta \epsilon }^2{\mathit{I}}_{N_{\xi }}$.

To simplify the notation, hereafter ${\mathit{\theta }}_{\mathsf{\delta }}=({\sigma }_{\delta \epsilon }^2, {\sigma }_{\delta v}^2)$, will represent the vector of variance parameters. The variance covariance matrix of ${\mathsf{\delta }}_{\xi }$ is:

$${\mathit{V}}_{\delta \xi }({\mathit{\theta }}_{\mathsf{\delta }})={\mathit{Z}}_{\delta \xi }{\mathit{G}}_{\delta }({\mathit{\theta }}_{\mathsf{\delta }}){\mathit{Z}}_{\delta \xi }^t+{\mathit{R}}_{\delta \xi }({\mathit{\theta }}_{\mathsf{\delta }}).$$

(4)

Model (3) is a linear mixed effect model and a special case of the basic unit-level described in [28] (pp. 174).

2.5.2. Target Parameter $\delta$-modeling Method

Under model (2) the target parameter (1) for a generic AOI $\alpha$ can be expressed as:

$${\Delta }_{\alpha }=\frac{1}{N_{\alpha }}{{\displaystyle \sum }}_{i=1}^{N_{\alpha }}{\delta }_{i\alpha }=\frac{1}{N_{\alpha }} {\mathbf{1}}^{\mathit{t}}{\mathsf{\delta }}_{\alpha }=\frac{1}{N_{\alpha }} {\mathbf{1}}^{\mathit{t}}({\mathit{X}}_{\delta \alpha }{\mathit{\beta }}_{\delta }+{\mathit{Z}}_{\delta \alpha }^t{\mathit{v}}_{\delta }+{\mathit{\epsilon }}_{\delta \alpha })={\mathit{l}}_{\delta \alpha }^t{\mathit{\beta }}_{\delta }+{\mathit{m}}_{\delta \alpha }^t{\mathit{v}}_{\delta }+{\mathit{q}}_{\delta \alpha }^t{\mathit{\epsilon }}_{\delta \alpha },$$

(5)

Thus, the target parameter is a linear target parameter similar to the one considered in [4] where ${\mathbf{1}}^{\mathit{t}}$ is a vector of ones and ${\mathit{l}}_{\delta \alpha }^t=\frac{1}{N_{\alpha }} {\mathbf{1}}^{\mathit{t}}{\mathit{X}}_{\delta \alpha }$, ${\mathit{m}}_{\delta \alpha }^t=\frac{1}{N_{\alpha }} {\mathbf{1}}^{\mathit{t}}{\mathit{Z}}_{\delta \alpha }^t$ and ${\mathit{q}}_{\delta \alpha }^t=\frac{1}{N_{\alpha }}{\mathbf{1}}^{\mathit{t}}$ are vectors of known constants for the target AOI $\alpha$.

2.5.3. Model Selection and Estimator $\delta$-modeling Method

The target parameter ${\Delta }_{\alpha }$ was estimated for all considered AOIs using ${\widehat{\Delta }}_{\alpha }$ the empirical best linear unbiased predictor (EBLUP) described in [29]. For each variable of interest, auxiliary variables included in ${\mathit{X}}_{\delta \alpha }$ were preselected using the best subset selection procedure described in [29] (pp. 179–180). When models with similar values of model root mean square error or coefficients of determination were compared, the preferred option was to select the model with smallest values of ${\sigma }_{\delta v}^2$. This criterion is appropriate to minimize the leading term of the AOI specific mean square errors [29] (pp. 176). Pre-selected models considered constant model error variances and no spatial correlation of model errors were fitted using maximum likelihood (ML). In a subsequent stage, models were re-fitted using ML including: (1) an exponential spatial correlation model for the model errors and (2) a non-constant error variance where ${\epsilon }_{\delta ij}\sim N(0,{\sigma }_{\delta \epsilon }^2{k}_{ij}^{2w_{\delta }})$. The term $k_{ij}$ is the value of the predictor included in the model most correlated to $\delta$ and $w_{\delta }$ is an additional parameter to account for heteroscedasticity. For all variables, no clear patterns of spatial correlation or non-constant variances were observed, which supports the model form described in Section 2.5.1.

Final estimates ${\widehat{\mathit{\theta }}}_{\mathsf{\delta }}$ of the variance parameters ${\mathit{\theta }}_{\mathsf{\delta }}$ were obtained using restricted maximum likelihood (REML) with the R [30] package nlme [31]. REML estimates ${\widehat{\mathit{\beta }}}_{\mathit{\delta }}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}$) of ${\mathit{\beta }}_{\delta }$ were functions of the estimated variance parameters (6):

$${\widehat{\mathit{\beta }}}_{\mathit{\delta }}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})={\{{\mathit{X}}_{\delta s}^t{\widehat{\mathit{V}}}_{\delta s}{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-1}{\mathit{X}}_{\delta s}\}}^{-1}{\mathit{X}}_{\delta s}^t{\widehat{\mathit{V}}}_{\delta s}{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-1}{\mathsf{\delta }}_s.$$

(6)

Matrices ${\widehat{\mathit{V}}}_{\delta s}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$, ${\widehat{\mathit{G}}}_{\delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ and ${\widehat{\mathit{R}}}_{\delta s}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ are obtained replacing the estimated variance parameters ${\widehat{\mathit{\theta }}}_{\mathsf{\delta }}$ in ${\mathit{V}}_{\delta s}({\mathit{\theta }}_{\mathsf{\delta }})$ ${\mathit{G}}_{\delta }({\mathit{\theta }}_{\mathsf{\delta }})$ and ${\mathit{R}}_{\delta s}({\mathit{\theta }}_{\mathsf{\delta }})$, by their REML estimates $\widehat{\mathit{\theta }}$. EBLUPs ${\widehat{\Delta }}_{\alpha }$ are also functions of $\widehat{\mathit{\gamma }}$ and are obtained using Equation (7):

$${\widehat{\Delta }}_{\alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})={\mathit{l}}_{\delta \alpha }^t{\widehat{\mathit{\beta }}}_{\mathit{\delta }}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})+{\mathit{m}}_{\delta \alpha }^t{\widehat{\mathit{v}}}_{\delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}),$$

(7)

where ${\widehat{\mathit{v}}}_{\delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ equals:

$${\widehat{\mathit{v}}}_{\delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})={\widehat{\mathit{G}}}_{\delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}){\mathit{Z}}_{\delta s}^t{V}_{\delta s}{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-1}\{{\mathsf{\delta }}_s-{\mathit{X}}_{\delta \alpha }{\widehat{\mathit{\beta }}}_{\mathit{\delta }}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})\}.$$

(8)

It is important to note that for AOIs in unsampled stands (i.e., pixels in unsampled compartments or the unsampled stands themselves), estimation will be made assuming that the model fit for the sampled stands also holds for the unsampled stands. Under that assumption, ${\mathit{m}}_{\delta \alpha }^t{\widehat{\mathit{v}}}_{\delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})=0$ and ${\widehat{\Delta }}_{\alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})={\mathit{l}}_{\delta \alpha }^t{\widehat{\mathit{\beta }}}_{\mathit{\delta }}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ is a synthetic predictor.

2.5.4. MSE Estimators for the $\delta$-modeling Method

For all AOIs, the mean squared error of the EBLUP was estimated using the estimator provided by [32] and extended in [4] to account for the fact that AOIs can contain a small number of population units. This estimator is the sum of three components where the last one, $2g_{3,\alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$, is a bias correction factor:

$$\widehat{MSE}\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})\}=g_{1\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})+g_{2\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})+2g_{3,\alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}),$$

(9)

The first term of (9) equals:

$$g_{1\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})={\mathit{m}}_{\delta \alpha }^t\{\widehat{\mathit{G}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})-\widehat{\mathit{G}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}){\mathit{Z}}_{\delta s}^t{\widehat{\mathit{V}}}_{\delta \alpha }{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-\mathbf{1}}{\mathit{Z}}_{\delta s}\widehat{\mathit{G}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})\}{\mathit{m}}_{\delta \alpha }+{\mathit{q}}_{\delta \alpha }^t{\mathit{R}}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}){\mathit{q}}_{\delta \alpha }.$$

(10)

The second term of (9) is:

$$g_{2\alpha \delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})={\mathit{d}}_{\delta \alpha }^t{\{{\mathit{X}}_{\delta s}^t{\mathit{V}}_{\delta s}{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-1}{\mathit{X}}_{\delta s}\}}^{-1}{\mathit{d}}_{\delta \alpha }$$

(11)

with ${\mathit{d}}_{\delta \alpha }^t={\mathit{l}}_{\delta \alpha }^t-{\mathit{m}}_{\delta \alpha }^t\widehat{\mathit{G}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}){\mathit{Z}}_{\delta s}^t{\widehat{\mathit{V}}}_{\delta s}{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-1}{\mathit{X}}_{\delta s}$. The term $g_{1\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ of $\widehat{MSE}\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})\}$ accounts for the uncertainty due to the estimation of the random effects while $g_{2\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ accounts for the uncertainty due to estimating ${\mathit{\beta }}_{\delta }$.

For model (2), it is possible to compute a bias correction factor for the mean square error estimator, that accounts for the uncertainty due to estimating ${\mathit{\theta }}_{\mathsf{\delta }}$. This correction factor equals:

$$g_{3\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})=tr\left\{\left({\frac{\partial {\mathit{b}}_{\delta \alpha }^t}{\partial {\mathit{\theta }}_{\mathit{\delta }}}|}_{{\widehat{\mathit{\theta }}}_{\mathit{y}}}\right){\mathit{V}}_{\delta s}{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-\mathbf{1}}{\left({\frac{\partial {\mathit{b}}_{\delta \alpha }^t}{\partial {\mathit{\theta }}_{\mathit{\delta }}}|}_{{\widehat{\mathit{\theta }}}_{\mathit{\delta }}}\right)}^t{\overline{\mathit{V}}}_{\delta s}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})\right\},$$

(12)

where, ${\mathit{b}}_{\delta \alpha }^t={\mathit{m}}_{\delta \alpha }^t\mathit{G}({\mathit{\theta }}_{\mathsf{\delta }}){\mathit{Z}}_{\delta s}^t{\mathit{V}}_{\delta s}{({\mathit{\theta }}_{\mathsf{\delta }})}^{-1}$ and ${\overline{H}}_{\delta s}(\widehat{\mathit{\theta }})$ is the Fisher information matrix for the fitted model. Explicit formulas for $g_{3\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ are provided [29] (pp. 179–180). This bias correction factor was used as a reference in comparisons with the $y$-modeling method.

All estimators of the mean square errors for AOIs in unsampled stands were made assuming that the model fitted for the sampled stands holds in the unsampled stands. Under this assumption, the leading term $g_{1\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$, of $\widehat{MSE}\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})\}$ will be larger than if the stand containing the AOI was sampled. This occurs because the negative term ${\mathit{m}}_{\delta \alpha }^t\widehat{\mathit{G}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}){\mathit{Z}}_{\delta s}^t{\widehat{\mathit{V}}}_{\delta s}{({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})}^{-\mathbf{1}}{\mathit{Z}}_{\delta \alpha }\widehat{\mathit{G}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }}){\mathit{m}}_{\delta \alpha }$ makes the term $g_{1\delta \alpha }(\widehat{\mathit{\theta }})$ smaller as the stand sample size increases.

2.6. $y$-Modeling Method

2.6.1. Model $y$-modeling Method

Models in the $y$-modeling method relate the forest structural variables in a population unit at different points in time with the auxiliary variables for that population unit. To indicate that these models directly consider the variables of interest, model parameters, stand-level random effects and model errors will include the subscript $y$. Auxiliary variables considered in the $y$-modeling method include the LiDAR auxiliary variables for 2009 and 2015, (i.e., Set 1 and Set 2, respectively) plus the incoming solar radiation and the factors structural diversity and presence of prescribed fires, Set 4 for both 2009 and 2015. The modeling for the method started obtaining models for the variable of interest for 2009 and models for the variable of interest in 2015.

For given time $t$, the variable of interest in the $j^{th}$ population unit in the $i^{th}$ stand is expressed as:

$$y_{ijt}={\mathit{x}}_{yijt}^t{\mathit{\beta }}_{yt}+u_{yit}+e_{yijt},$$

(13)

where ${\mathit{x}}_{yijt}^t$ is a vector of auxiliary variables, specific for time $t$, in which the first element takes the value 1. The term ${\mathit{\beta }}_{yt}$ is a vector of time-specific coefficients with the first element representing the model intercept. The random components of model (13) are the stand-level random effects $u_{yit}$ and the model errors $e_{yijt}$. To account for heteroscedasticity, model errors $e_{yijt}$ were of the form $e_{yijt}={\epsilon }_{yijt}{k}_{ijt}^{{\omega }_{yt}}$ with ${\epsilon }_{yijt}\sim N(0,{\sigma }_{y\epsilon t}^2)$; the term $k_{ijt}$, the predictor included in the model for time $t$, is most correlated to $y_t$, and ${\omega }_{yt}$ is a parameter to model the change in the error variance. The stand-level random effects $u_{yit}$ were assumed to be independently and identically distributed (i.i.d.) normal random variables $u_{yit}\sim N(0,{\sigma }_{yut}^2)$ for all $i=1,\dots ,D$, where $D$ is the total number of stands in the study area. Model errors were assumed independent of the stand-level random effects (i.e., $Cov(e_{yijt},u_{ykt})=0$, for all $i,j$ and $k$). Finally, model errors were consider independent with $Cov({\epsilon }_{yijt},{\epsilon }_{yklt})=0$ if $i\ne k$ or $j\ne l$ for both $t=2009$ and $t=2015$. Models with spatially correlated errors were initially considered, but discarded for both years in the selection stage as no clear spatial correlation patterns were observed in the residuals. Auxiliary variables included in the model were selected following the same procedure used in the $\delta$-modeling method, using the best subsets selection procedure described in [27].

To account for expected correlations, models for 2009 and 2015 were combined into a single model where stand-level random effects and model errors for 2009 and 2015 were allowed to be time correlated. Then for the $j^{th}$ population unit in the $i^{th}$ stand the two-dimensional vector ${\mathit{y}}_{ij}={(y_{ij09},y_{ij15})}^t$ of variables of interest was related to the auxiliary variables through the following model:

$${\mathit{y}}_{ij}={\mathit{X}}_{ij}{\mathit{\beta }}_y+{\mathit{B}}_{ij}{\mathit{v}}_{yi}+{\mathit{e}}_{yij}$$

(14)

with:

$${\mathit{X}}_{ij}=\left(\begin{array}{cc}{\mathit{x}}_{yij09}^t& {\mathbf{0}}_{p15}^t\\ {\mathbf{0}}_{p09}^t& {\mathit{x}}_{yijt15}^t\end{array}\right),{\mathit{\beta }}_y=\left(\begin{array}{c}{\mathit{\beta }}_{y09}\\ {\mathit{\beta }}_{y15}\end{array}\right), {\mathit{B}}_{ij}=\left(\begin{array}{c}\begin{array}{ccc}1& 1& 0\end{array}\\ \begin{array}{ccc}1& 0& 1\end{array}\end{array}\right), {\mathit{v}}_{yi}=\left(\begin{array}{c}{v}_{yi}\\ v_{yi09}\\ v_{yi15}\end{array}\right),{\mathit{e}}_{yij}=\left(\begin{array}{c}{e}_{yij09}\\ e_{yij15}\end{array}\right),$$

(15)

where ${\mathbf{0}}_{p2009}^t$ and ${\mathbf{0}}_{p2015}^t$ are, respectively, row vectors of zeroes of dimensions equal to ${\mathit{x}}_{yij09}^t$ and ${\mathit{x}}_{yij2015}^t$.

As with the time-specific models, to account for heteroscedasticity in the combined model, model errors $e_{yijt}$ were of the form $e_{yijt}={\epsilon }_{yijt}{k}_{ijt}^{{\omega }_{yt}}$ with ${\epsilon }_{yijt}\sim N(0,{\sigma }_{y\epsilon t}^2)$. The parameters ${\omega }_{yt}$ were updated when fitting the combined model. Spatial correlation of model errors was not found to be significant when considering each year separately, therefore, no spatial correlation patterns were considered in the combined model. The only source of correlation of model errors present in the combined model was temporal correlation. For a given location, the variables ${\epsilon }_{yijt}{}_{yij09}\sim N(0,{\sigma }_{y\epsilon 09}^2)$ and ${\epsilon }_{yij15}\sim N(0,{\sigma }_{y\epsilon 15}^2)$ were allowed to be correlated random variables. The correlation between ${\epsilon }_{yij2009}$ and ${\epsilon }_{yij2015}$ is ${\rho }_{\epsilon }$ and the variance-covariance matrix of ${\mathit{e}}_{yij}$ is:

$$Cov\left(\begin{array}{c}{e}_{yij09}\\ e_{yij15}\end{array}\right)={\mathit{R}}_{yij}=\left(\begin{array}{cc}{\sigma }_{y\epsilon 09}^2{k}_{ij09}^{2{\omega }_{y09}}& {\rho }_{\epsilon }{\sigma }_{09}{k}_{ij09}^{{\omega }_{y09}}{\sigma }_{09}{k}_{ij15}^{{\omega }_{y15}}\\ {\rho }_{\epsilon }{\sigma }_{09}{k}_{ij09}^{{\omega }_{y09}}{\sigma }_{09}{k}_{ij15}^{{\omega }_{y15}}& {\sigma }_{yv}^2{k}_{ij15}^{2{\omega }_{y15}}\end{array}\right).$$

(16)

To model correlation between stand-level random effects, three random components $v_{yi}$, $v_{yi2009}$ and $v_{yi2015}$, independent of each other, were considered. These random components had distributions $v_{yi}~N(0,{\sigma }_{yv}^2$), $v_{yi09}~N(0,{\sigma }_{yv09}^2)$ and $v_{yi15}~N(0,{\sigma }_{yv15}^2)$. Stand-level random effect for a given point, at time $t$, $u_{yit}$, are the sum of a pure stand effect, independent of time $t$, $v_{yi}$, and a time-specific stand random effect $v_{yi09}$ or $v_{yi15}$. The term ${\mathit{B}}_{ij}{\mathit{v}}_{yi}={\mathit{u}}_{yi}=\left(\begin{array}{c}{v}_{yi}+v_{yi09}\\ v_{yi}+v_{yi15}\end{array}\right)=\left(\begin{array}{c}{u}_{yi05}\\ u_{yi09}\end{array}\right)$ represents these sums. The variance covariance matrix of ${\mathit{v}}_{yi}$ is diagonal, therefore, the variance covariance matrix of ${\mathit{u}}_{yi}$ is:

$$Cov\left(\begin{array}{c}{u}_{yi09}\\ u_{yi15}\end{array}\right)={\mathit{G}}_{yi}=\left(\begin{array}{cc}{\sigma }_{yv}^2+{\sigma }_{yv15}^2& {\sigma }_{yv}^2\\ {\sigma }_{yv}^2& {\sigma }_{yv}^2+{\sigma }_{yv15}^2\end{array}\right)$$

(17)

The fact that the random effect $v_{yi}$ is present for both 2009 and 2015 results in a positive correlation of the terms of ${\mathit{u}}_{yi}$, with a correlation coefficient ${\rho }_u=\frac{{\sigma }_{yv}^2}{({\sigma }_{yv}^2+{\sigma }_{yv09}^2)({\sigma }_{yv}^2+{\sigma }_{yv15}^2)}$. In a last step, models with a simpler structure of random effects were fitted and compared to the original models using a likelihood ratio test. Simplified models contained only random effects $v_{yi}$ that did not depend on time (i.e., models did not contain time-specific random effects $v_{yi09}$ and $v_{yi15}$). For simplified models $u_{yi09}=u_{yi15}=v_{yi}$.

For a generic set of population units $\xi$, the combined model can be expressed in matrix notation as:

$${\mathit{y}}_{\xi }={\mathit{X}}_{y\xi }{\mathit{\beta }}_y+{\mathit{Z}}_{y\xi }{\mathit{v}}_y+{\mathit{e}}_{y\xi },$$

(18)

where ${\mathit{y}}_{\xi }$, ${\mathit{e}}_{y\xi },$ and ${\mathit{X}}_{y\xi }$, are obtained stacking the vectors ${\mathit{y}}_{ij}$, ${\mathit{e}}_{yij},$ or the matrices ${\mathit{X}}_{ij}$ of all units in $\xi$. As no spatial correlation patterns were found, the variance covariance matrix of ${\mathit{e}}_{y\xi }$ is, ${\mathit{R}}_{y\xi }=dia{g}_{i,j\in \xi }({\mathit{R}}_{yij})$, a block diagonal matrix of dimension $2N_{\xi }x2N_{\xi }$ with $2x2$ blocks equal to ${\mathit{R}}_{yij}$. The vector of stand-level random effects $\mathit{v}={({\mathit{v}}_{\mathit{i}1}^t,{\mathit{v}}_{\mathit{i}2}^t, \dots ,{\mathit{v}}_{\mathit{i}D}^t)}^{\mathit{t}}$ and the matrix ${\mathit{Z}}_{y\xi }$ is an incidence matrix of dimension $2N_{\xi }xD$ for the simplified models and $2N_{\xi }x3D$ for the models with time-specific random effects. The variance covariance matrix of ${\mathit{y}}_{\xi }$ can be expressed as:

$${\mathit{V}}_{y\xi }({\mathit{\theta }}_{\mathit{y}})={\mathit{Z}}_{y\xi }{\mathit{G}}_y({\mathit{\theta }}_y){\mathit{Z}}_{y\xi }^{\mathit{t}}+{\mathit{R}}_{y\xi }({\mathit{\theta }}_y).$$

(19)

In Equation (19), it is explicitly indicated that matrices ${\mathit{V}}_{y\xi }({\mathit{\theta }}_y)$, ${\mathit{G}}_y({\mathit{\theta }}_y)$ and ${\mathit{R}}_{y\xi }({\mathit{\theta }}_y)$ depend on the vector of variance-covariance parameters ${\mathit{\theta }}_y={({\sigma }_{yv}^2,{\sigma }_{yv09}^2,{\sigma }_{yv15}^2,{\sigma }_{y\epsilon 09}^2,{\sigma }_{y\epsilon 15}^2,{\rho }_{\epsilon })}^t$. For the models with simplified random effects the vector of variance covariance parameters reduces to ${\mathit{\theta }}_y={({\sigma }_{yv}^2,{\sigma }_{y\epsilon 09}^2,{\sigma }_{y\epsilon 15}^2,{\rho }_{\epsilon })}^t$. Model (18) is a special case of linear mixed effect model with block diagonal covariance structure.

2.6.2. Target Parameter $y$-modeling Method

Under model (18) the target parameter (1) for a generic AOI, $\alpha$ is a linear combination of the form:

$${\Delta }_{\alpha }=\frac{1}{6N_{\alpha }}{{\displaystyle \sum }}_{i=1}^{N_{\alpha }}(y_{i\alpha 15}-y_{i\alpha 09})={\mathit{l}}_{y\alpha }^t{\mathit{\beta }}_{\delta }+{\mathit{m}}_{y\alpha }^t{\mathit{u}}_y+{\mathit{q}}_{y\alpha }^t{\mathit{e}}_{y\alpha },$$

(20)

where ${\mathit{l}}_{y\alpha }^t={\mathit{q}}_{y\alpha }^t{\mathit{X}}_{y\alpha }$, ${\mathit{m}}_{y\alpha }^t={\mathit{q}}_{y\alpha }^t{\mathit{Z}}_{y\alpha }$ and ${\mathit{q}}_{y\alpha }^t$ are vectors of known constants for the target AOI $\alpha$, with ${\mathit{q}}_{y\alpha }^t$ a vector of dimension $2N_{\alpha }$ where the $k^{th}$ element equals $\frac{{(-1)}^k}{6N_{\alpha }}$. It is important to remark that for models with a simplified structure of stand random effects, the target parameters do not depend on ${\mathit{u}}_y$. For these models, $y_{ij15}-y_{ij15}=({\mathit{x}}_{yij15}^t{\mathit{\beta }}_{y15}-{\mathit{x}}_{yij09}^t{\mathit{\beta }}_{y09})+(e_{yij15}-e_{yij09})$, and $u_{yi09}-u_{yi15}=v_{yi}-v_{yi}=0$. For these type of models, one can expect significant gains in accuracy because it is not necessary to estimate random effects.

2.6.3. Estimator $y$-modeling Method, and Estimator of the MSE

Model (18) is a linear mixed effects model with block diagonal structure and ${\Delta }_{\alpha }$ a linear model parameter; thus, after [29] (pp. 108–110), the EBLUP ${\widehat{\Delta }}_y({\widehat{\mathit{\theta }}}_y)$ of ${\Delta }_{\alpha }$ is:

$${\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)={\mathit{l}}_{y\alpha }^t{\widehat{\mathit{\beta }}}_{\mathit{y}}({\widehat{\mathit{\theta }}}_y)+{\mathit{m}}_{y\alpha }^t{\widehat{\mathit{v}}}_y({\widehat{\mathit{\theta }}}_y),$$

(21)

where ${\widehat{\mathit{\beta }}}_{\mathit{y}}({\widehat{\mathit{\theta }}}_y$) equals:

$${\widehat{\mathit{\beta }}}_{\mathit{y}}({\widehat{\mathit{\theta }}}_y)={\{{\mathit{X}}_{ys}^t{\widehat{\mathit{V}}}_{ys}{({\widehat{\mathit{\theta }}}_y)}^{-1}{\mathit{X}}_{ys}\}}^{-1}{\mathit{X}}_{ys}^t{\widehat{\mathit{V}}}_{ys}{({\widehat{\mathit{\theta }}}_y)}^{-1}{\mathit{y}}_s.$$

(22)

Matrices ${\widehat{\mathit{V}}}_{ys}({\widehat{\mathit{\theta }}}_y)$, ${\widehat{\mathit{G}}}_y({\widehat{\mathit{\theta }}}_y)$ and ${\widehat{\mathit{R}}}_{ys}({\widehat{\mathit{\theta }}}_y)$ are obtained by replacing the estimated variance parameters ${\mathit{\theta }}_y$ in ${\mathit{V}}_{y\xi }({\mathit{\theta }}_y)$, ${\mathit{G}}_y({\mathit{\theta }}_y)$ and ${\mathit{R}}_{y\xi }({\mathit{\theta }}_y)$, by their REML estimates ${\widehat{\mathit{\theta }}}_y$. EBLUPs ${\widehat{\Delta }}_{\alpha }$ are also functions of $\widehat{\mathit{\gamma }}$ and are obtained using formula (7), where ${\widehat{\mathit{v}}}_{\delta }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ equals:

$${\widehat{\mathit{v}}}_y({\widehat{\mathit{\theta }}}_y)={\widehat{\mathit{G}}}_y({\widehat{\mathit{\theta }}}_y){\mathit{Z}}_{ys}^t{V}_{ys}{({\widehat{\mathit{\theta }}}_y)}^{-1}\{{\mathit{y}}_s-{\mathit{X}}_{y\alpha }{\widehat{\mathit{\beta }}}_y({\widehat{\mathit{\theta }}}_y)\}.$$

(23)

As with the $\delta$-modeling method, estimates for AOIs in unsampled stands were made assuming that the model fit for the sampled stands also applied in the unsampled stands, which leads to ${\mathit{m}}_{y\alpha }^t{\widehat{\mathit{v}}}_y({\widehat{\mathit{\theta }}}_y)=0$ and ${\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)={\mathit{l}}_{y\alpha }^t{\widehat{\mathit{\beta }}}_{\mathit{y}}({\widehat{\mathit{\theta }}}_y)$ is a synthetic predictor.

For all AOIs, the estimator of the mean square error of the EBLUP under the $y$-modeling method, $\widehat{MSE}\{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)\}$, is:

$$\widehat{MSE}\{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)\}=g_{1y\alpha }({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})+g_{2y\alpha }({\widehat{\mathit{\theta }}}_y).$$

(24)

The terms $g_{1y\alpha }({\widehat{\mathit{\theta }}}_y)$ and $g_{2y\alpha }({\widehat{\mathit{\theta }}}_y)$ in (24) are analogous to those in (10) and (11) and have similar interpretation. To compute $g_{1y\alpha }({\widehat{\mathit{\theta }}}_y)$ and $g_{2y\alpha }({\widehat{\mathit{\theta }}}_y)$, matrices ${\widehat{\mathit{G}}}_{\mathit{\delta }\mathit{s}}({\widehat{\mathit{\theta }}}_{\delta })$, ${\widehat{\mathit{R}}}_{\mathit{\delta }\mathit{s}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$, ${\widehat{\mathit{V}}}_{\mathit{\delta }\mathit{s}}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ and ${\widehat{\mathit{R}}}_{\mathit{\delta }\mathit{\alpha }}({\widehat{\mathit{\theta }}}_{\mathsf{\delta }})$ must be replaced by ${\widehat{\mathit{G}}}_{ys}({\widehat{\mathit{\theta }}}_y)$, ${\widehat{\mathit{R}}}_{ys}({\widehat{\mathit{\theta }}}_y)$, ${\widehat{\mathit{V}}}_{\mathit{y}\mathit{s}}({\widehat{\mathit{\theta }}}_y)$ and ${\widehat{\mathit{R}}}_{\mathit{y}\mathit{\alpha }}({\widehat{\mathit{\theta }}}_y)$. For the $y$-modeling method we did not compute the second-order correction factors.

2.7. Comparison of Methods

Methods were compared using three different criteria. First, we used general measures of accuracy providing the average error or uncertainty of prediction at the pixel-level (2.7.1); then, we compared methods using AOI-specific estimates and measures of uncertainty (2.7.2). Finally, we assessed the risk of generating biased predictions when using the $\delta$-modeling method and $y$-modeling method in unsampled stands (Section 2.7.3).

2.7.1. General Accuracy Assessment

To compare the $\delta$-modeling method and $y$-modeling method, a fist assessment was made using the cross-validated model mean squared error, $mRMSE$, and the model bias $mBias$:

$$mRMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i,j\in s}^n{({\delta }_{ij}-{\widehat{\delta }}_{ij})}^2}{n}},$$

(25)

$$mBias=\frac{{{\displaystyle \sum }}_{i,j\in s}({\delta }_{ij}-{\widehat{\delta }}_{ij})}{n},$$

(26)

where ${\delta }_{ij}$ is the observed value of change for the $j^{th}$ plot included in the $i^{th}$ sampled stand and ${\widehat{\delta }}_{ij}$ is the predicted value for that plot when model coefficients are obtained removing that plot from the training dataset. For the $y$-modeling method ${\widehat{\delta }}_{ij}$ is obtained using the observed and fitted values of the variable of interest, as ${\widehat{\delta }}_{ij}=\frac{1}{6}({\widehat{y}}_{ij15}-{\widehat{y}}_{ij09})$ where ${\widehat{y}}_{ij09}$ and ${\widehat{y}}_{ij15}$ are the predictions of $y_{ij09}$ and $y_{ij15}$ are obtained fitting the corresponding $y$-model without the observations for plot $ij$. In addition, we computed $mRMSE$ and $mBias$ in terms relative to the average changes observed in the sampled plots. These quantities are denoted as $mRRMSE$ = $mRMSE/{\widehat{\Delta }}_f$ and $mRBias=mBias/{\widehat{\Delta }}_f$ where ${\widehat{\Delta }}_f$ is the mean of the changes observed in the field plots.

2.7.2. AOI-specific Comparisons.

For each of the considered areas of interest an estimate by each method (i.e., EBLUPs using either ${\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_y)$ or ${\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)$) and their corresponding mean square error estimators (i.e., $\widehat{MSE}\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_y)\}$ or $\widehat{MSE}\{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)\}$) were available. First, for each AOI and method, we directly compared ${\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_y)$ and ${\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)$, and the square roots of $\widehat{MSE}\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\delta })$ and $\widehat{MSE}\{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)$. To simplify the notation, we will omit the subscript indicating the target AOI unless it is necessary and refer to ${\widehat{\Delta }}_{\delta }$ as ${\widehat{\Delta }}_y$. Similarly, after omitting the subscript $\alpha$, the AOI specific root mean square errors will be denoted as:

$$RMS{E}_{\delta }=\sqrt{\widehat{MSE}\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\delta })\}},$$

(27)

$$RMS{E}_y=\sqrt{\widehat{MSE}\{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)\}},$$

(28)

To perform an assessment relative to the predicted values the following coefficient of variation:

$$CV\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\delta })\}=\frac{\sqrt{\widehat{MSE}\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\delta })\}}}{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_{\delta })},$$

(29)

$$CV\{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)\}=\frac{\sqrt{\widehat{MSE}\{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)\}}}{{\widehat{\Delta }}_{y\alpha }({\widehat{\mathit{\theta }}}_y)},$$

(30)

was computed for each AOI and method. Finally, for each AOI we compared $CV\{{\widehat{\Delta }}_{\delta \alpha }({\widehat{\mathit{\theta }}}_y)\}$ to $CV\{{\widehat{\Delta }}_{y\alpha }({\widehat{\theta }}_y)\}$. To simplify the notation we will refer to these coefficients of variation as $C{V}_{\delta }$ and $C{V}_y$. Finally, for each sampled AOI we computed, using only the field information, the sample mean ${\widehat{\Delta }}_{f\alpha }$ and its standard error:

$$S{E}_{f\alpha }=\sqrt{\frac{{{\displaystyle \sum }}_{i,j\in s}^{n_{\alpha }}{({\delta }_{ij}-{\widehat{\Delta }}_{f,\alpha })}^2}{(n_{\alpha }-1)n_{\alpha }}},$$

(31)

and its coefficient of variation $C{V}_{f\alpha }$. In Equation (31), $n_{\alpha }$ is the number of field plots in the considered AOI and the sub-index $f$ is used to indicate that these quantities are calculated using only field data. Again, to simplify the notation we removed the sub-indexes $\alpha$ unless they were necessary. Finally, the sample mean and the coefficient of variation were then compared to their counterparts (7) and (29) and (21) and (30) obtained by the $\delta$-modeling method and $y$-modeling method, respectively.

2.7.3. Extrapolation to Thinned Stands

The fact that thinned stands were not represented in the sample of field plots raises the question of how applicable the models obtained are using either the $\delta$-modeling method or the $y$-modeling method to these stands. Applying the models to these stands involves a degree of extrapolation to a different population and a high risk of producing biased predictions. We assessed this risk by comparing the distributions of the LiDAR predictors included in the models for the $\delta$-modeling method and the $y$-modeling method for the sample of field plots to the distributions of the predictors in: (1) the sampled stands, (2) the unsampled and not thinned stands and (3) unsampled and thinned stands. Within each group (i.e., field plots, FP; sampled stands not thinned, SS; unsampled stands not thinned, UN; and unsampled stands subject to thinning, UT), we estimated density functions for each LiDAR predictor using a Gaussian kernel and a bandwidth determined using Silverman’s rule [33]. Note that we considered two AOIs for the largest level of aggregation, the first one is SS and the other one, SA, is the union of SS and UN. We first considered each predictor separately and graphically compared their density functions. Predictors for 2009 and 2015 in the $y$-modeling method were considered separately. For each predictor and group, we computed the area of overlap, $AO$, with the density function for the field plots which takes value 0 if there is no overlap and value 1 if the distribution of the predictor in the considered group equals the distribution for the sample.

In addition to the area of overlap and aiming to consider all predictors in a given model at once, we calculated $\overline{NT2}$, the average of Mesgaran’s novelty index $NT2$ for each model and group [34]. This quantity provides the average Mahalanobis distance from the pixels of the group to the mean of the sample of field plots, and it is expressed in terms relative to the maximum Mahalanobis distance observed in the sample. Values of $\overline{NT2}$ above one indicate that on average pixels in a group are at a distance to the mean of the field plots larger than the distance from the extreme field observation to the mean of the field plots. We also calculated ${\overline{NT2}}_{mean}$, the average of $NT2$, but using the mean Mahalanobis distance as normalizing constant instead of the maximum. The reference value of one for ${\overline{NT2}}_{mean}$ indicates that the average Mahalanobis distance from pixels, to the mean of the field plots, is the same as the average of the Mahalanobis distances observed in the sample. Means and variance covariance matrices for computation of Mahalanobis distances are always estimated using the sample of field plots.

3. Results

3.1. Selected Models $\delta$-modeling Method and $y$-modeling Method

Selected models for the $\delta$-modeling method included auxiliary variables from Set 1 and Set 3. It was possible to find alternative models including fixed effects for the diversity treatments (i.e., predictors from Set 4) with similar values of $mRMSE$ and $mRBias$; however, those models did not improve the model fit. From a practical point of view, models that only depend on the LiDAR variables but do not depend on the structural diversity treatments or the presence/absence of prescribed fires make the models more portable and applicable to stands without needing to know exactly which one of these treatments was applied. Considering that models using the structural diversity and presence of prescribed fires as predictors did not result in important gains in accuracy, we selected models that were not dependent on these treatments (

Table 2. Summary models for the $\delta$-modeling method. Model coefficients, standard errors of the model coefficients, variance parameters and general metrics for accuracy assessment are provided. Predictor acronyms are explained in Table A1. Coef is the value of the coefficient and Std.Error its corresponding standard error. V indicates volume, BA indicates basal area and B indicates biomass.

Model	Predictor	Coef	Std. Error	${\widehat{\mathit{\sigma }}}_{\mathit{\delta }\mathit{v}}^2$	${\widehat{\mathit{\sigma }}}_{\mathit{\delta }\mathit{\epsilon }}^2$	$\mathit{m}\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{m}\mathit{R}\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{m}\mathit{B}\mathit{i}\mathit{a}\mathit{s}$	$\mathit{m}\mathit{R}\mathit{B}\mathit{i}\mathit{a}\mathit{s}$
V(m3 ha−1 year−1)	Intercept	1.16	0.31	0.50	10.53	3.47	241.99%	−1.83 × 10−4	−0.01%
	$\delta$Elev_P5015-09	1.33	0.27
	$\delta$PcFstAbv215-09	0.23	0.07
BA(m2 ha−1 year−1)	Intercept	0.31	0.12	0.01	0.14	0.39	116.30%	−8.2 × 10−4	−0.25%
	$\delta$PcAllAbv215-09	0.05	0.01
	Elev_P7509	−0.03	0.01
	PcAllAbv215-09	0.02	<0.01
B(Mg ha−1 year−1)	Intercept	1.03	0.17	0.19	2.52	1.72	180.20%	−1.09 × 10−3	−0.11%
	$\delta$Elev_var15-09	0.05	0.02
	$\delta$Elev_P5015-09	1.03	0.20
	$\delta$CRR15-09	−16.67	6.58

).

For the models in the $\mathsf{\delta }$-modeling method, the variance of random the effects, ${\widehat{\sigma }}_{\delta v}^2$, was very small compared to the variance of the model errors, ${\widehat{\sigma }}_{\delta \epsilon }^2$, (Table 2). This indicated that, in this forest and for these variables, the use of synthetic estimators that do not account for the variability between stands should not cause a strong bias problem.

Models for the $y$-modeling method showed a pattern similar to that observed for the $\delta$-modeling and only included predictors from Set 1 and Set 2 (

Table 3. Summary models for the $y$-modeling method. Model coefficients, standard errors of the model coefficients, variance-covariance parameters and general metrics for accuracy assessment are provided. Covariate acronyms are explained in Table A1. Coef is the value of the coefficient and Std.Error its corresponding standard error. Coef is the value of the coefficient and Std.Error its corresponding standard error. V indicates volume, BA indicates basal area and B indicates biomass.

Model	Year	Covariate	Coef	Std.Error	${\widehat{\mathit{\sigma }}}_{\mathit{y}\mathit{u}}^2$	Kijt	${\mathit{\omega }}_{\mathit{y}\mathit{t}}$	${\widehat{\mathit{\sigma }}}_{\mathit{y}\mathit{t}\mathit{\epsilon }}^2$	${\mathit{\rho }}_{\mathit{e}}$	General Accuracy Metrics for Change Per Hectare and Year
										$\mathit{m}\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{m}\mathit{R}\mathit{R}\mathit{M}\mathit{S}\mathit{E}$	$\mathit{m}\mathit{B}\mathit{I}\mathit{A}\mathit{S}$	$\mathit{m}\mathit{R}\mathit{B}\mathit{i}\mathit{a}\mathit{s}$
V(m3 ha−1)	2009	Intercept	−19.09	10.36	640.29	Elev_mean209	0.64	3.00	0.85	3.76	262.62%	0.13	9.24%
		Elev_mean209	2.52	0.23
		PcFstAbv209	0.63	0.05
	2015	Intercept	2.69	0.23		Elev_mean215	0.61	4.17
		Elev_mean215	0.69	0.05
		PcFstAbv215	−26.30	11.10
BA(m2 ha−1)	2009	Intercept	−0.22	1.57	7.42	PcFstAbv209	0.48	0.81	0.85	0.47	138.06%	0.01	1.53%
		Elev_P1009	−1.37	0.34
		Elev_P3009	1.58	0.24
		PcFstAbv209	0.51	0.03
	2015	Intercept	−2.16	0.61		PcFstAbv215	0.45	1.12
		Elev_P1015	2.56	0.51
		Elev_P2015	0.57	0.03
		PcFstAbv215	−0.97	1.72
B(Mg ha−1)	2009	Intercept	−11.86	5.19	165.69	Elev_mean209	0.71	0.39	0.85	1.94	203.69%	0.08	8.60%
		Elev_mean209	1.19	0.11
		PcFstAbv209	0.34	0.02
	2015	Intercept	1.31	0.12		Elev_mean215	0.58	1.47
		Elev_mean215	0.37	0.03
		PcFstAbv215	−14.15	5.76

). Errors showed non-constant variance patterns for all variables. The predictor most correlated with the variable of interest (i.e., the predictor used to model the error variance) was the same for 2009 and 2015. For V and B, the variance of model errors was a function of the square of the mean LiDAR elevation (Elev_mean²), and the exponents of the error variance function were very close to those obtained in [4,17,27] for V, and in [27] for B. For BA, variance of model errors was a function of the percentage of first returns above two meters (PcFstAbv2). Based on the results of the likelihood ratio tests, that for all variables resulted in p-values larger than 0.87, simplified models were selected and used for prediction.

3.2. General Accuracy Assessment and Comparison of Methods

For all variables and modeling alternatives, values of $mBias$ and $mRBias$ were orders of magnitude smaller than $mRMSE$ and $mRRMSE$ (Table 2 and Table 3). For all variables and methods, the percentages of explained variance for the change in V, BA and B were low. For the $\delta$-modeling method, models explained 34.38%, 31.37% and 39.04% of the variance of the change in V, BA and B, respectively. For the $y$-modeling method, models explained only 10.65% and 5.37% of V and B, respectively, while for BA the prediction using the $y$-modeling method was not better than the sample mean. In addition, $\delta$-models had values of $mRBias$ lower than those obtained for the $y$-models. When instead of the change we considered the forest structural attributes with the $y$-modeling method, percentages of explained variance were 82.16% for V, 82.53% for BA and 82.93% for B. Considering only 2009, the percentage of explained variance for V, BA and B was 81.60%, 83.45%, and 82.84%, respectively. Considering only 2015, the percentage of explained variance for V, BA and B was 82.72%, 81.42% and 82.98%, respectively.

3.3. AOI-Specific Estimates

3.3.1. Entire Study Area

Estimates for the sampled stands and for the whole study area using either the $\delta$-modeling method or the $y$-modeling method were consistent with the estimates obtained using only the field information except for BA and B in SA. For the entire study area values of $RMS{E}_{\delta }$ tended to be smaller than $RMS{E}_y$. When considering the sampled stands, SS, approximate confidence intervals computed as ${\widehat{\Delta }}_f\pm 2S{E}_f$ for the field estimates, and as ${\widehat{\Delta }}_{\delta }\pm 2RMS{E}_{\delta }$ and ${\widehat{\Delta }}_y\pm 2RMS{E}_y$ for each one of the LiDAR based methods, overlapped for all variables (

Table 4. Average increments of volume V, basal area BA and biomass B in the entire study area excluding the thinned stands (SA) and for the union of the sampled stands (SS). Estimates $(\widehat{\Delta })$, root mean square errors $(RMSE)$, coefficients of variation $(CV)$ and confidence intervals ($CI$) obtained using the $\delta$-modeling method and the $y$-modeling method are compared to estimates (${\widehat{\Delta }}_f$), standard errors ($S{E}_f$) coefficients of variation ($C{V}_f$), and confidence intervals ($C{I}_f$) using only the field information.

Variable	Area	$\mathit{\delta }$-modeling Method					$\mathit{y}$-modeling Method					Field Only Estimates
		${\widehat{\Delta }}_{\mathit{\delta }}$	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{\delta }}$	$\mathit{C}{\mathit{V}}_{\mathit{\delta }}$	$\mathit{C}{\mathit{I}}_{\mathit{\delta }}$		${\widehat{\Delta }}_{\mathit{y}}$	$\mathit{R}\mathit{M}\mathit{S}{\mathit{E}}_{\mathit{y}}$	$\mathit{C}{\mathit{V}}_{\mathit{y}}$	$\mathit{C}{\mathit{I}}_{\mathit{y}}$		${\widehat{\Delta }}_{\mathit{f}}$	$\mathit{S}{\mathit{E}}_{\mathit{f}}$	$\mathit{C}{\mathit{V}}_{\mathit{f}}$	$\mathit{C}{\mathit{I}}_{\mathit{f}}$
V(m3 ha−1 year−1)	SS	1.66	0.27	16.29%	1.12	2.20	1.95	0.32	16.48%	1.31	2.60	1.43	0.32	22.21%	0.80	2.07
	SA	1.67	0.30	17.98%	1.07	2.27	1.98	0.29	14.67%	1.40	2.56
BA(m2 ha−1 year−1)	SS	0.36	0.03	8.68%	0.30	0.42	0.37	0.04	9.93%	0.30	0.45	0.34	0.04	10.87%	0.26	0.41
	SA	0.42	0.04	8.41%	0.35	0.49	0.44	0.04	9.61%	0.35	0.52
B(Mg ha−1 year−1)	SS	1.07	0.13	12.35%	0.81	1.34	1.24	0.17	13.61%	0.90	1.57	0.95	0.16	16.89%	0.63	1.28
	SA	1.15	0.16	13.66%	0.83	1.46	1.29	0.15	11.83%	0.98	1.59

) and contained estimates derived from other methods. Differences between the uncertainty of estimates obtained from LiDAR-based methods and the uncertainty of estimates obtained from field-based methods tend to be of small magnitude.

3.3.2. Stands

Estimated change of V, BA and B in the sampled stands by both the $\delta$-modeling method and the $y$-modeling method agreed with their field-based counterparts in most stands (Figure 2). However, the width of the confidence intervals obtained using the $\delta$-modeling method tended to be larger than the confidence intervals of the estimates derived using the $y$-modeling method (Figure 2).

For unsampled stands, estimates and confidence intervals had larger variability in stands where the forested area was small (Figure 3). This variability cannot be avoided, and indicates that certain sources of errors cannot be compensated if the number of pixels that are aggregated is low. Finally, for both methods, values of $RMS{E}_{\delta }$ and $RMS{E}_y$ were in the range of 0.25 to 1 m³ ha⁻¹year⁻¹ for V, of 0.02 to 0.15 m² ha⁻¹year⁻¹ for BA and of 0.10 to 0.80 Mg ha⁻¹year⁻¹ for B. However, for B and V, the $RMS{E}_y$ tended to be smaller than $RMS{E}_{\delta }$ while negligible differences between methods were observed for BA (Figure 3 and Figure 4).

For the thinned (and unsampled) stands, differences between $\delta$-modeling method and the $y$-modeling method for BA were large and their confidence intervals did not overlap (Figure 3). For these stands, the estimates for BA using the $\delta$-modeling method tended to indicate almost no changes in BA. Estimates for the thinned stands using the $\delta$-modeling method provided inconsistent results indicating gains in B, and changes close to zero for V and BA. Certain inconsistencies were also observed for stands subject to thinning when using the $y$-modeling method where predictions of the change in V and B were positive for three and five stands respectively. These inconsistencies seem to derive from the fact that the distribution of predictors in Set 3 (i.e., changes in LiDAR predictors) in the thinned stands was rather different to the distribution of these predictors in the sample of field plots, in the sampled stands and in the unsampled and not thinned stands. For predictors of the $y$-modeling method, modeled differences between thinned stands and the remaining groups were of much smaller magnitude. Results for the analysis of the extrapolation risks are presented in detail in Section 3.4.

3.3.3. Pixel-level

For both methods, inconsistencies observed at the stand level were observed at the pixel-level, especially the positive predictions of change obtained with the $\delta$-modeling method in the thinned stands (Figure A1). In addition, due to the low correlations of LiDAR predictors with the change in V, BA, and B, predictions at this level have large uncertainties. Mean and median values of $RMS{E}_{\delta }$ were 2.30 m³ ha⁻¹ year⁻¹ and 3.30 m³ ha⁻¹ year⁻¹ for V, 0.39 m² ha⁻¹ year⁻¹ and 0.38 m² ha⁻¹ year⁻¹ for BA, and 1.67 Mg ha⁻¹ year⁻¹and 1.65 Mg ha⁻¹ year⁻¹ for B. Mean and median values of $RMS{E}_y$ were 2.49 m³ ha⁻¹ year⁻¹ and 2.20 m³ ha⁻¹ year⁻¹ for V, 0.48 m² ha⁻¹ year⁻¹ and 0.48 m² ha⁻¹ year⁻¹ for BA and 1.89 Mg ha⁻¹ year⁻¹ and 1.76 Mg ha⁻¹ year⁻¹ for B (

Table 5. Minimum (Min), 5th percentile (p05), mean, median, 95th percentile (p95) and maximum (Max) of $RMS{E}_{\delta }$ (27) and $RMS{E}_y$ (28) for the pixels of the study area.

Variable	Method	Min	p05	Mean	Median	p95	Max
V(m3 ha−1 year−1)	$\delta$-modeling method	0.42	0.42	2.30	3.30	3.59	9.41
	$y$-modeling method	0.08	0.37	2.49	2.20	6.01	32.69
BA(m2 ha−1 year−1)	$\delta$-modeling method	0.38	0.38	0.39	0.38	0.40	0.59
	$y$-modeling method	0.11	0.30	0.48	0.48	0.64	1.47
B(Mg ha−1 year−1)	$\delta$-modeling method	1.62	1.63	1.67	1.65	1.76	4.57
	$y$-modeling method	0.47	1.10	1.89	1.76	3.09	10.45

and Figure A1). Predictions from the $\delta$-modeling method tend to be smoother than predictions from the $y$-modeling method. For all variables, the proportion of pixel-level predictions using the $\delta$-modeling method within the range of values observed for the field plots, was always 99.84% or larger (Figure A2). Considering that these results were obtained in the presence of thinned stands and the relatively small fraction of the forest that was sampled, obtaining less than 0.16% of the predictions outside of the measurement range seems to be a clear sign of over smoothing (see Appendix B).

3.4. Extrapolation to Thinned Stands

Estimates of change in B for the thinned stands by both methods were clearly subject to bias problems. The predicted change in B for the total area subject to thinning for the period 2009–2015, using the $\delta$-modeling method was an increase in biomass of 40,469.22 Mg. The predicted change using the $y$-modeling method was a removal of B. However, the predicted removal for the period 2009–2015 was only 1750.29 Mg while the weighted extractions for the thinned stands were orders of magnitude larger. For BA both methods estimated extractions in BA, which is consistent with the fact that these stands were thinned. Estimated changes in BA using the $\delta$-modeling method for the thinned stands ranged from −0.05 m² ha⁻¹ to −1.88 m² ha⁻¹, which seems to be a very small change in basal area. Estimated changes in BA using the $y$-modeling method ranged from −3.58 m² ha⁻¹ to −7.01 m² ha⁻¹. An advantage of the $y$-modeling method is that it allows obtaining the values of the structural attributes at a given point in time. Using the $y$-model we estimated BA for the thinned stands for 2015. For those stands where thinning prescriptions dictated leaving a residual BA of 17.22 m² ha⁻¹ to 25.25 m² ha⁻¹, estimated BA for 2015 ranged from 19.87 m² ha⁻¹ to 26.22 m² ha⁻¹, which is in accordance with the thinning prescriptions. For the remaining area subject to thinning the estimated BA for 2015 was 17.64 m2 ha⁻¹, while the prescriptions dictated leaving a residual BA ranging from 6.89 m² ha⁻¹ to 13.77 m² ha⁻¹ in 75% of the area and leaving the remaining area untouched. In general, the estimated BA for 2015 are consistent with the prescriptions, which indicates that the $y$-modeling method produces reasonable estimates of BA when extrapolating to the thinned stands. In summary, for the estimation of changes, biases derived from extrapolation seemed to be of larger magnitude for the $\delta$-modeling method although they were also present for the $y$-modeling method.

The extrapolation indexes $\overline{NT2}$ and ${\overline{NT2}}_{mean}$ showed that predictions in thinned stands, involved a large amount of extrapolation when using the $\delta$-modeling method. For the $y$-modeling method, differences between thinned stands and stands not subject to thinning were of smaller magnitude (Figure 5 and Figure A3). The inspection of the distribution of the LiDAR predictors in the field plots, the sampled and not thinned stands, unsampled and not thinned stands and the unsampled and thinned stands showed similar results for all variables, being the distributions of predictors from Set 3 (i.e., changes in LiDAR predictors for the period 2009–2015) very sensitive to the thinning operations (Figure 6, Figure A4 and Figure A5).

4. Discusion

4.1. General Accuracy Assessment and Comparison of Methods.

The smallest values of $mRMSE$ were obtained using the $\delta$-method, which is consistent with previous results reported by Poudel et al. [8] for V and B in coastal coniferous forest of Western Oregon and by Temesgen et al. [9] for B in spruce-dominated forest of Alaska (Table 2 and Table 3). We observed, however, smaller differences between methods. Additionally, as observed in previous studies, (e.g., [8,21,35,36]) where LiDAR auxiliary variables showed a much stronger correlation with structural attributes at a given point in time than with their change.

Values of $mRMSE$ for V were 3.47 m³ ha⁻¹year⁻¹ when using the $\delta$-method and 3.76 m³ ha⁻¹ year⁻¹ when using the $y$-method. These values are slightly smaller than the $mRMSE$ obtained by Poudel et al. [8] using the $\delta$-method (4.74 m³ ha⁻¹ year⁻¹) and two lidar acquisitions separated in time by five years. For B, $mRMSE$ using the $\delta$-method and the $y$-method were, respectively, 1.72 Mg ha⁻¹year⁻¹ and 1.94 Mg ha⁻¹ year⁻¹. These values were very close to those reported by Poudel et al. [8] using the $\delta$-method (1.88 Mg ha^-1 year^-1) and worse than those reported by Temesgen et al. (1.25 Mg ha⁻¹year⁻¹ and 1.63 Mg ha⁻¹ year⁻¹), also using two LiDAR acquisitions separated in time by five years. Values of $RMSE$ for BA were similar to those obtained by Næsset and Gobakken [20] in coniferous forest in Norway, using the $y$-method with log-transformed models and two LiDAR acquisitions that were two years apart from each other. In relative terms, for V and B, the values that we obtained for $mRRMSE$ were considerably larger than those obtained by Poudel et al. [8]. These differences are due to the fact that observed growth rates in Poudel et al [8] are much higher than we observed at BMEF.

4.2. AOI-Specific Estimates

4.2.1. Entire Study Area

Most studies on estimation of change of structural variables using repeated LiDAR measurements have focused on analyzing indexes of model fit and reported only global measures of accuracy developed at a plot level. There is an important difference between the values of $RMS{E}_{\delta }$ and $RMS{E}_y$ and $mRMS{E}_{\delta }$ and $mRMS{E}_y$ being $mRMS{E}_{\delta }$ and $mRMS{E}_y$ an order of magnitude larger than $RMS{E}_{\delta }$ and $RMS{E}_y$. Model root mean square errors $mRMS{E}_{\delta }$ and $mRMS{E}_y$ provide an average measure of the errors that can occur when predicting a single pixel. For large areas, there will be some level of compensation of overpredicted and underpredicted pixels. Knowing how important that compensation is requires calculating AOI-specific root mean square errors. These AOI-specific measures cannot be directly derived from $mRMSE$ because $RMS{E}_{\delta }$ and $RMS{E}_y$ consider factors such as the uncertainty in the estimation of the fixed and random effects that are not accounted for in $mRMS{E}_{\delta }$ or $mRMS{E}_y$. The effect of these factors in $RMS{E}_{\delta }$ and $RMS{E}_y$ can cause that the way two models with similar values of $mRMSE$ rank based on this metric could change when attending to $RMSE$. But the most important consequence of the disconnect between $mRMS{E}_{\delta }$ and $mRMS{E}_y$ and AOI-specific measures of uncertainty, is that the former cannot be used as quality controls in LiDAR based inventories.

While numerous studies on estimation of changes using LiDAR rely on global measures of accuracy such as $mRMSE$, exceptions to this trend can be found in the literature [10,12,18,19,20]. The last four studies used model assisted techniques to derive either landscape or stratum level changes. Reported errors in those studies changed depending on the modeling techniques and study areas, but they all were of similar magnitude for changes in biomass per hectare and year (Table 4). Errors for the methods tested in this study were smaller than those reported by [10], where changes in live carbon stocks in Norway were estimated using generalized regression estimators (GREG). Differences with the errors reported in [10] for carbon, using the same 0.5 biomass to carbon conversion factor, were in the range of 0.12–0.09 Mg ha⁻¹ year⁻¹. These differences seem to be due to multiple factors such as differences between study areas, changes in live biomass versus changes in standing biomass, time between LiDAR acquisitions and field plot sizes etc. Further investigation is needed to test if the model-based estimators studied here and the GREG estimators in [10] have a similar performance when used under the same conditions.

The study from Magnussen et al. [18] also included model based estimators using the $y$-modeling method. Reported errors were slightly larger than the ones observed here but at the same time smaller in terms relative to the observed mean change. An important result from the comparisons of [18] was the drastic improvement in model accuracy when developing stratum specific models (i.e., a set of model coefficients per stratum) as opposed to a global model for the whole study area. The mixed effect models used in this study can be used in combination with stratification if sample sizes are large. The introduction of stand level random effects allows for certain variability between AOIs that can be applied in situations where AOI sample sizes are limited.

4.2.2. Stands

One of the novelties of this study was the analysis of estimates for AOIs with small sample sizes to develop AOI specific models (i.e., stands). While at large scales both LiDAR-based and field-based estimates were very similar and had equivalent accuracies, at the stand-level, LiDAR based estimates, clearly had smaller errors than their field-based counterparts do. Qualitatively, this result for the change in V, BA and B is similar to the results obtained in [15,17] for the structural variables themselves and shows that the LiDAR auxiliary information allows for gains in efficiency when estimating changes in AOIs with small sample sizes. However, due to the low correlation of LiDAR and structural changes, values of $C{V}_{\delta }$ and $C{V}_y$ were oftentimes larger than 50%. These values of $CV$ are larger than those observed for structural variables in similar AOIs in previous studies [4,14,17]. While differences were not of large magnitude $RMS{E}_y$, tended to be smaller than $RMS{E}_{\delta }$. In addition, $RMS{E}_y$ had a larger variability because errors did not have constant variance. Finally, stand level estimates using the $\delta$-modeling method in the thinned and unsampled stands involved an important degree of extrapolation that can cause inconsistent estimates and severe biases, which indicates that the $\delta$-modeling method is more sensitive to extrapolating than the $y$-modeling method.

4.2.3. Pixel-level

For the most detailed level of disaggregation, the magnitude of the errors was very large. This is due to the low correlation between LiDAR auxiliary variables and the change in structural attributes. First-order and second-order texture indexes [37] are auxiliary variables with a promising potential for future research aiming to improve the prediction of structural changes. While for structural variables, maps at the pixel-level can provide a reliable reference about the forest structure; for growth and changes, pixel-level maps like the one in Figure A1 should be taken as mere approximations. They could be used to infer certain trends and patterns, but the high values of $RMS{E}_{\delta }$ and $RMS{E}_y$ show that estimates for a particular location made at the pixel scale can differ significantly with reality. These results clearly indicate that, predictions at such a fine scale are highly unreliable, and it is necessary either to perform some level of spatial aggregation or to increase the lapse between LiDAR acquisitions.

4.3. Advantages of Modeling Alternatives

In general, the $\delta$-modeling method was found to be a better alternative to estimate changes for the entire study area than the $y$-modeling method; however, the $y$-modeling method produced better results at the stand-level and also seemed to be advantageous to prevent problems related to extrapolation to values of the covariates outside of those included in model development.

The $\delta$-modeling method offers a faster model developments and fitting, and is significantly simpler than the modeling with the $y$-modeling method, as it is not necessary to consider differences between years and time correlations. The main disadvantage of this method is that it seems to be more prone to extrapolation errors. Predictors from Set 3 are sensitive to intense changes in the forest structure caused for example by thinning (see Figure 5, Figure A4, Figure A5 and Figure 6). The inspection of predictors of alternative models for V and B using this method revealed that inconsistencies of predictions in unthinned stands could be attenuated including more predictors from Set 1. The sensitivity to changes of predictors from Set 3 can be an advantage if all possible changes are correctly represented in the field sample. However, for relatively short periods of time between acquisitions and a low amount of forest operations, changes that are not very frequent in the landscape can be misrepresented or even not included in the sample. Thus, results for areas subject to those changes can be severely biased and inconsistent.

The more complex model development for the $y$-modeling method may be compensated by its ability to produce a richer set of outputs, by its apparently smaller risk of extrapolation and by its more accurate estimates for AOIs with small sample sizes (i.e., stands). In this study we analyzed the performance of the $y$-modeling method when estimating change, but estimates of V, BA and B for all the AOIs in 2009 and 2015 could have been readily obtained using this method. Results from our study also support the idea that the $y$-modeling method has advantages over the $\delta$-modeling method in terms of protection against inconsistent extrapolations. The distributions of predictors from Set 1 and Set 2 in thinned stands were relatively similar to the distributions observed for the sample while the distributions of predictors in Set 3 used in the $\delta$-modeling method, these distributions were rather different (see Figure 5, Figure A4, Figure A5 and Figure 6). The greater similarity between thinned stands and the sample of field plots, for predictors from Set 1 and Set 2, indicates that the effect of thinning, in terms of auxiliary information, can be seen as transition from one situation in 2009 to another in 2015, and both seem to be represented in the field sample (e.g., Figure 5). If structures before and after the thinning (or other changes) are represented in the sample, the need for extrapolation will be limited. Within certain limits, if the sampling design covers all structures present at both points in time, even if there is a particular change from one structure to another that is not represented in the sample, predictions from the $y$-modeling method will not involve large extrapolations.

5. Conclusions

The four main conclusions obtained from this study include:

• The change of structural attributes and LiDAR auxiliary information are weakly correlated. This weak correlation seems to more evident in BMEF than in previous studies because of the slower growth in the study area and the relatively short lapse of time between LiDAR acquisitions, which indicates that for future studies in similar areas it might be necessary to increase the time lags between LiDAR flights.
• In general, the $\delta$-modeling method was found to be a slightly more accurate alternative to obtain estimates of change for the whole study area; however, the $y$-modeling method was able to produce better estimates at the stand level. In addition, the $y$-modeling method method also seemed to be less prone to extrapolation problems. This indicates that field campaigns for the $\delta$-modeling method have to be carefully designed while the $y$-modeling method might be less sensitive to certain bias problems.
• Despite the weak correlations with the changes in structural attributes, LiDAR auxiliary information allows obtaining estimates of growth for stands that improve over those derived using only field information.
• The large uncertainty observed for pixel-level predictions indicated that high-resolution maps of growth, generated using LiDAR auxiliary information in similar conditions, should be taken as approximated products.