TY - JOUR
AU - Gobakken, Terje
AB - Design-based estimators for two-stage simple random sampling with regression can have a lack of precision (efficiency) when the primary sampling units (PSUs) vary in size, and PSU totals are approximately proportional to the size of a PSU. Precision and efficiency may deteriorate further for domain-specific estimators when PSUs contain elements from different domains. Design model-unbiased ratio-to-size estimators have been proposed as more efficient. This study introduces a variance estimator for a design model-unbiased ratio estimator. The estimator of variance is derived from a single-stage estimator of a variance of a ratio under two assumptions: the target variable (y) is equal to the sum of a model prediction and two error terms capturing residual errors and model estimation errors; and the existence of an unbiased estimator of model parameters. Extensive simulations confirmed the negative effects of unequal PSU sizes on the precision of the model-assisted estimators of variance and a superior performance of the proposed estimator of variance. These results were confirmed in the analysis of a regional two-stage survey of forest biomass. The proposed variance estimator was generally more efficient than an existing alternative and more stable across a large suite of design settings. sampling design, ratio-to-size estimator, forest inventory, large-scale survey, auxiliary variables The use of remotely sensed auxiliary variables has a long history in forest resource surveys (Spurr 1952, p. 424, Loetsch and Haller 1973, p. 299). Today, auxiliary variables are used on a routine basis in regional and national surveys (Tomppo 1991, Bauer et al. 1994, McRoberts et al. 2010b). Auxiliary variables obtained from airborne laser scanning (ALS) are valued for their relationship to stem volume, biomass, basal area, canopy height, canopy density, and vertical canopy structure (Næsset 2004, Hyyppä et al. 2008, Hudak et al. 2009). Numerous studies have developed robust and reliable methods for the extraction of inventory target variables (Y) and powerful predictors of Y from ALS data (e.g., Breidenbach et al. 2010, Dalponte et al. 2012, Magnussen et al. 2012). Model-assisted (MA), and model-dependent (MD) estimators of totals (means) and variances have been worked out for many types of multiphase forest survey designs (Schreuder et al. 1995, Köhl et al. 2006, Gregoire and Valentine 2008). In contrast with those of an MD estimator, the bias properties of a MA estimator do not depend on the correctness of an estimated model (Särndal et al. 1992, Chapter 13.6). Variance reductions attributed to the auxiliary variables can be substantial (Coomes et al. 2002, Opsomer et al. 2007, McRoberts et al. 2010a). Technological advances in remote sensing and progress in application research indicate a continuation of progress (Köhl et al. 2006, p. 196). In multiphase (Särndal et al. 1992, p. 344) and multistage designs (Särndal et al. 1992, p. 133), models are typically used to link auxiliary variables (X) obtained at different stages (phases) of a design to one or more target variable(s) (Y). The models are intended as a means to an end: the estimation and inference about a finite population parameter (total, mean) that lead the analyst to a design-based inference (Särndal et al. 1992, p. 513). Four issues in large regional natural resource surveys, however, can set the stage for an imprecise or a questionable choice of design-based estimators: (1) logistical and technical considerations combined with the spatial extent of the sample frame often lead to equal-probability designs with unequal-sized primary sampling units (PSUs); (2) PSUs may contain secondary sampling elements from different strata, administrative units, and categorical classes; (3) the combination of high costs associated with direct observations (measurements) of Y and relatively low unit-level costs of X favors designs in which the sample size of Y is a very small fraction of the sample size of X; and (4) a desire to obtain a spatially balanced sample (Stevens and Olsen 2004) encourages designs with a systematic selection of units and elements. The recent Hedmark County survey (Næsset et al. 2013) is a good example of how the four issues combine to pose a complex challenge for the analyst. We use this survey as an illustrative case study. This survey targeted a forested area of approximately 27,390 km2, divided into 636 parallel strips to facilitate airborne data collection. Strips were of equal width but different length, and each contained parcels of varying forest cover and administrative status. A systematic sample of strips was flown, with limited ground-based measurements for forest inventory attributes taken at regular intervals within these strips. Design-based estimators for equal-probability unequal-sized PSU designs (issue 1) can be imprecise when PSU sizes vary greatly but means per element vary little from PSU to PSU (Cochran 1977, p. 249). This follows because the expected value of the squared PSU size is greater than the squared mean PSU size (Cochran 1977, 11.24 and 11.30). Precision deteriorates as the variability in PSU size increases. Ratio-to-size estimators will often be more efficient in these situations (Cochran 1977, p. 303). Neither MD (e.g., Ståhl et al. 2011b) nor least squares prediction estimators (Royall 1976, Valliant et al. 2000, p. 266) offer a solution to this problem. For issue 2, in textbooks on sample designs (e.g., see Thompson 1992, Fuller 2009), a sample unit belongs to one and only one subpopulation. A PSU defined solely by logistic modalities surrounding the capture of auxiliary variables, without regard to the structure of the population to be surveyed, may contain elements from different strata, classes, and domains of interest. This fact compounds the first issue of unequal-sized PSUs and also confounds the use of a nominal second-stage sample inclusion probability for elements in a PSU in a two-stage design-based estimator (Särndal et al. 1992, Chapter 8). Very low sample fractions in stage two of a two-stage design (issue 3) means that second-stage π estimators (Fuller 2009, p. 8) become sensitive to outliers and random variation in the number of stage-two elements in a certain stratum, class, or domain (issue 2). Finally (issue 4), a systematic or semisystematic selection of sample units, is common in large-scale natural resource surveys (Spencer and Czaplewski 1997, Opsomer et al. 2007, Gasparini et al. 2010, Ståhl et al. 2011a) and effectively precludes the validity of a design-based estimator of variance (Gregoire and Valentine 2008, p. 55); in practice, however, the analyst frequently resorts to estimators for simple random sampling (SRS) (McRoberts et al. 2010b). An analyst faced with the challenge of analyzing results from a two-stage with regression survey of natural resources, in which one or a combination of the issues discussed above is deemed important, may look for a suitable “robust” alternative to design-based estimation (Ghosh 2008). For the specific case of a two-stage with regression sampling design, Ståhl et al. (2011b) suggested a ratio-regression estimator for the mean per element and proposed a new MD variance estimator. A direct application in Norway demonstrated considerable gains but only after poststratification (Ståhl et al. 2011b, Ene et al. 2012, Gobakken et al. 2012). In an attempt to improve on this situation, this study proposes an adaptation of the ratio-to-size estimators for a single-stage sampling design for SRS (Cochran 1977, Chapter 6.2–6.4) to a two-stage with regression SRS design. Ratio estimators can effectively address the earlier discussed issues 1–3 that can lead to imprecise design-based estimates. MD estimators are typically derived from model predictions and their properties, but here we begin with a design model-unbiased ratio estimator. To succeed with an adaptation from this starting point, we regard the true element-level value of Y as being approximately equal to the sum of a model-based prediction and two independent zero-mean random errors: one due to residual variance and one due to model estimation error. With this approximation, the ratio-to-size estimator and its variance are obtained by replacing Y with its prediction plus the two unknown random errors. Error variances and covariance are estimated from statistics of the sample and the fitted model (Draper and Smith 1981, p. 19 and 83). The proposed ratio estimator and the estimator of variance are evaluated and compared with conventional two-stage with regression MA design-based estimators (Särndal et al. 1992, Chapter 8.9) and with the MD estimator of variance proposed by Ståhl et al. (2011b). Note that our proposed ratio estimator is practically identical to the one proposed by Ståhl et al. (2011b), whereas our proposed variance estimator is different. Comparisons are made in the context of simulated two-stage sampling from two artificial populations and results from an actual survey. Materials and Methods Population and Sampling Designs We consider a finite area population (U) composed of a known number (N) of elements (U = {1, …, k,…, N}) of equal size (area = a ha) and shape. Each element belongs to one of nAU administrative units (AUs) and to one of nH classes (H) defined by intrinsic characteristics (e.g., land cover) of an element. The number of elements and areas (A) in each AU and H is known. Associated with each element is one or more attribute(s) of interest (Y). It is desirable to estimate the total of Y (ty), or the per hectare average of Y (μy), for various domains: the entire population (U), AUs, and individual classes. Estimation methods for domains follow the direct estimation principles outlined in Särndal et al. (1992, Section 10.3). Sample sizes for an estimation for subdomains, e.g., classes within a given AU, are considered too small for a direct estimation (Gregoire et al. 2010). Estimation for these subdomains will have to adopt approaches tailored specifically to small-area estimation problems (Fuller 2009, p. 311). Because the population is too large to afford a census, the desired estimates are obtained from a probability sample of elements (Hansen et al. 1983). We decided to use probability sampling in two stages. In the first sample (S1), a number (P) of low-cost auxiliary variables (X) assumed to be associated with Y are observed in a large sample of elements (s1). The second-stage sample (s2) is much smaller and forms a subset of s1. Costly observations of Y are made on each element in s2. Logistics and the nature of X call for large compact PSUs in S1 followed by an independent equal probability subsampling of elements within each PSU (Cochran 1977, p. 300). Specifically, the population is divided into L PSUs of which l are selected by simple random sampling without replacement (SRSwor). The technique used to collect X largely determines the size, shape, and orientation of individual PSUs. Neither AU nor H is considered during the partitioning of U to PSUs. Resulting PSUs will, therefore, typically contain a different number (m) of elements. In addition, the elements in a PSU may belong to different AUs and different Hs. To avoid repetitive presentations of similar estimators for different domains in U, we provide elements, estimators, and estimates with the symbol ϕ to identify the set of elements that forms the domain for which an estimate is desired. Predicting Y from X The S2 sample (s2) serves to estimate a linear function (f) for the relationship between Y (or some transform of Y) (Lahiri and Rao 1995) and the auxiliary variable(s) X. For the jth element in the ith PSU, the linear model elements in domain ϕ is assumed to be   where β0ϕ is a domain-specific intercept and βpϕ is a regression coefficient of the pth auxiliary variable in domain ϕ and εijϕ is a domain-specific residual with an expected value of 0 and a variance of σε2 (ϕ). In our study with equal probability sampling, the ordinary least squares techniques can be applied to obtain the regression coefficients in Equation 1. Otherwise we refer readers to 6.4.2 on p. 225 in Särndal et al. (1992). In the cases presented here, we used a single global regression model that includes domain-specific intercepts and domain-specific regression coefficients. Least squares estimates of the parameters in Equation 1 are denoted as , and a model-based prediction of yij is indicated as . Sample-based estimates of a parameter or a statistic are modified with a caret () to distinguish them from their actual value. With sufficient S2 sampling in each S1 PSU, the analyst might also entertain an expansion of Equation 1 with random PSU effects (Stanek and Singer 2004). When Y is transformed to achieve a linear model, the analyst must take steps to correct for any bias in final back-transformed estimates on the natural scale of Y (Li and Lahiri 2007). For an example on the back transformation of a square root transform of Y, see Gregoire et al. (2008). Regression Estimators for Two-Stage SRSwor Sampling For a two-stage equal probability design, the direct MA estimator for the total of Y in elements in domain ϕ becomes (for a generic formulation, see Särndal et al. 1992, p. 323, 8.9.6)   where mi is the number of elements in the ith S1-selected PSU, ni is the S2 sample size in the ith PSU, ij is the residual yij − and δC is an indicator value that takes the value of 1 when the event C is true and 0 otherwise. The equal probability direct variance estimator for is (for a generic formulation, see Särndal et al. 1992, p. 326, 8.9.27)   where (zk) is the estimated variance of z over indices of k (Cochran 1977, p. 26). Note that these variances are computed for elements in s2 and domain ϕ. As in Gregoire et al. (2010), they are simple (unweighted) moment estimators of variance; we recognize that a potentially more efficient g-weighting would be possible. Corresponding per unit area (hectare) estimators for domain (ϕ) become . Negative values of are a possibility because of the subtraction of the scaled estimate of the within-PSU variance (W). In addition, small domain-specific s2 subsamples can produce erratic estimates in . An alternative and more robust direct MA variance estimator, based on the variance of the regression estimates of PSU totals, has been suggested (Särndal et al. 1992, p. 154, 4.6.2)   A forestry application with this estimator is provided in Ene et al. (2012, Equation 9). Ratio-Regression Estimators for Two-Stage SRSwor Sampling When either PSU size or the number of domain-specific elements in a PSU vary greatly, the estimator in Equation 2 and, in particular, the variance estimators in Equations 3 and 4 can be of poor precision (Cochran 1977, p. 249, Rao 1988, Fuller 2009, p. 97). This occurs when domain means per element (of Y) vary little among PSUs, and domain-specific PSU totals are approximately proportional to the number of elements contributing to the total. Low S2 sample fractions, exemplified by expensive natural resource surveys, amplify the variability in estimates obtained from Equations 2, 3, and 4. Illustrative examples of these problems can be found in Ene et al. (2012) and in Næsset et al. (2013). Särndal et al. (1992, p, 327, 8.9.30) suggested a more robust design-based ratio-regression estimator of μy that addressed the issue of unequal PSU sizes but not issues related to possible side effects of excessive scaling of residuals in designs with very low S2 sampling intensities. Provided that the expected value under the sampling design of f (xij;) is unbiased (or nearly unbiased), a more suitable design model-unbiased ratio estimator comes in the form of a conventional one-stage ratio-regression estimator for μy (Cochran 1977, p. 31). The domain (ϕ)-specific direct per element ratio estimator given in Equation 5 addresses both the issue of unequal PSU sizes and low S2 sample inclusion probabilities   Note that because we have an SRS design in both stages and our regression model includes a domain-specific intercept (i.e., the sum of domain-specific regression residuals is 0), the estimator in Equation 5 is identical to the ratio regression estimator given in 7.13.7 in Särndal et al. (1992). The corresponding per unit area estimator is where a, as defined earlier, is the area of a population element. The associated estimator of a population total is , where Lϕ is the number of PSUs in the population with at least one element in domain ϕ and ϕ is the S1 sample estimate of the mean number of domain ϕ elements in PSUs with at least one element in domain ϕ (i.e., ϕ = lϕ−1Σi=1lmiϕ). Note that membership of elements to a domain is inferred from available maps. A comparison of and reveals differences in scaling to a total (i.e., L × l−1 versus Lϕ × lϕ−1) and scaling of S2 sums of residuals (i.e., mi × n−1i versus 1 for all domains). Ratio estimators, like the one in Equation 5 are known to be nearly unbiased for large sample sizes (i.e., Σi=1lΣjmi δij∈ϕ = Σi=1lmiϕ) and when PSU totals of Y are (approximately) proportional to the size of a PSU (Fuller 2009, 1.3.50). Ståhl et al. (2011b) adopted, for a two-stage with regression design, the average of model-based predictions as the estimator of the per element mean of Y (i.e., using for s2 and non-s2 elements in contrast with Equation 5). The two estimators are identical only when they include domain-specific regression variables, viz., domain indicators, which guarantees that residuals sum to 0 within a domain. The proposed variance estimator for RTϕ is a modification of a variance estimator for a ratio estimator with SRSwor of unequal size PSUs (Cochran 1977, p. 155, 6.10 or 6.13). The modification follows by replacing the observed values of the target variable (y) by a model-based prediction plus the sum of two independent error terms. The resulting direct estimator of variance is given in Equation 6, and details of derivation and notation appear in the Appendix   where φ is the variance effective number of second-stage PSUs (see A.3 in the Appendix), ϕ is the average number of ϕ elements in the S1 sample of PSUs, and fpcs is a finite population correction factor for stage S sampling (Cochran 1977, p. 303). With SRSwor S1 sampling of l PSUs included out of a total of L PSUs, the domain-specific correction factor is fpc1ϕ = (1 − lϕ/Lϕ). For SRSwor S2 sampling with n observations of yij and a model (f) with np parameters, the domain-specific correction factor was fixed at fpc2ϕ = min[1, 1 − (nϕ − np)/(Lϕ)], where nφ is the number of elements from domain ϕ in the S2 sample. Finally, is the estimate of the covariance matrix of the domain-specific regression model and is the within-PSU covariance matrix of the explanatory variables in the domain-specific regression model. The first right-hand side term of Equation 6 captures the variance due to sampling (of ) and the second term captures the covariance among PSUs of model-based estimates of PSU totals in domain ϕ. Finally, the ratio-regression variance estimator of per unit area estimate y, RTϕ becomes RT(y, RTϕ) = a−2RT(RTϕ). An Alternative Model-Dependent Variance Estimator For a ratio estimator practically identical to the estimator in Equation 5, Ståhl et al. (2011b, Equation 15) proposed the following direct variance estimator   As before, the conversion of GS(RTϕ) to a variance estimator of per unit area estimate y, RTϕ follows from GS(y, RTϕ) = a−2GS(RTϕ). A comparison between the proposed estimator in Equation 6 and the estimator in Equation 7 reveals an absence of finite population correction factors in Equation 7, the use of ϕ in lieu of lϕ(ϕ ≤ lϕ)in the denominator in Equation 6, and the absence of the term in the numerator of Equation 7. The combined effect of these differences can be either positive or negative. Assessment of Estimators The two-stage with regression estimators in Equations 2, 3, 4, and 7 and the proposed ratio estimators in Equations 5 and 6 are assessed in simulated two-stage SRSwor sampling from two artificial populations. The estimators y, R2ϕ and y, RTϕ are assessed for bias (estimated minus actual value) and root mean squared error (RMSE = square root of sampling variance plus squared bias), whereas the sampling variance is estimated as the variance of per hectare estimates in repeated (simulated) sampling. The square root of the average replicate value of the analytical estimators , , (y, RTϕ), and GS(y, RTϕ) is compared with the square root of the corresponding empirical replicate variance of and y, RTϕ [denoted as rep() and rep(y, RTϕ), respectively]. If two nearly unbiased estimators have similar RMSEs and approximately unbiased estimators of variance, then a variance estimator with a lower variability should be preferred, because it provides a more consistent estimate of variance than a more variable estimator. Thus, we also consider the sampling distribution of variance estimates. Confidence intervals provide a good summary of the uncertainty in an estimate. We therefore computed normal quantile-based 95% confidence intervals (CI95) from per hectare estimates of Y and their associated estimates of variance. We report on the coverage of these intervals, i.e., the proportion of intervals that included the true value of Y (Casella and Berger 2002, p. 454). Coverage for the two-stage with regression estimators is denoted as CI95R2 when based on Equations 2 and 3, CI95R1 when derived from Equations 2 and 4, CI95RT when computed from Equations 5 and 6, and CI95GS when computed from Equations 5 and 7. Test Population POP1 Test population POP1 has N elements split into L compact PSUs of equal size (m = N × L−1) in some design settings and of unequal sizes, distributed uniformly at random, in others. There is one target variable (Y) and one correlated auxiliary variable (X) associated with each element. Element values of Y and X were generated from a multivariate t distribution (MVT) with λ degrees of freedom. Thus, the relationship between X and Y is linear, yij = β0 + β1xij + εy|x, ij. Specifically,   with the per unit area means of Y and X fixed at μy = 100, μx = 10, and   within-PSU variances of Y and X fixed at σy2 = 140 and σx2 = 6. Conditional on the PSU effects, the correlation between Y and X is ρyx and factor γ controls the magnitude of the per element among PSU variance in Y and X relative to the per element within variance in Y and X. To examine the effects of ρxy, γ, and λ on the assessed estimators, a total of 12 different distributions of X and Y were generated by introducing four levels in each of the parameters ρxy, γ, and λ. Four distributions were generated with ρxy = 0.3, 0.5, 0.7, or 0.9, while λ and γ were held constant at 15 and 30, respectively. A second set of four distributions had γ = 15, 30, 45, or 60, with ρxy = 0.7 and λ = 15. The third set of four had λ = 5, 15, 25, or 35, with ρxy = 0.7 and γ = 30. Four population sizes, with fixed and uniform random variation in PSU sizes, were then generated from each of the above 12 joint distributions of Y and X to give a total of 96 test populations. Population sizes (N) were 30 × 106 (L = 600, = 5 × 104), 4 × 106(L = 400, = 104), 4 × 105(L = 200, = 2 × 104), and 5 × 104(L = 100, = 500). In populations with uniform random PSU sizes, the minimum size was 0.7 × and the maximum size was 1.3 × (i.e., a coefficient of variation in size of approximately 22%). Eight two-stage SRSwor sample designs were simulated, four with S1 inclusion probabilities of 0.06, 0.08, 0.10, and 0.12 with a fixed S2 sample size of 10 in each S1-selected PSU, and four with a fixed S1 inclusion probability of 0.08 and S2 sample sizes of 5, 10, 15, and 20. For each sample size, a total of 40,000 samples were taken from each of the 96 fixed populations to keep Monte Carlo errors (Koehler et al. 2009) of an estimated variance at <0.5%. Test Population POP2 Test population POP2 was intended as a realistic example of a two-stage forest inventory in a 1,608,575-ha forested region with S1 sampling of two auxiliary variables (X1, X2) and S2 sampling of a target variable Y representing stem wood volume (vol, m3 ha−1). The region has four (spatially compact) administrative units (AUs) and eight land cover classes (H) of which the first five represent elements with different forest characteristics. The five forest classes constitute the sample frame (approximately 70% of POP2). An element is a 1-ha square area (a = 1 ha). Elements were assigned an AU according to location. POP2 is divided into 3,144 100-m wide east-west-oriented PSUs (strips) to accommodate S1 sampling of X. The number of elements in a PSU varies from 24 to 782 with a mean of 512 and a SD of 177. Estimates of μy are desired for the four AUs (AU1, …, AU4) and for the five forest classes (H1, …, H5). Figure 1 shows the outline of POP2, the four AUs, and an example of a systematic sample of 66 PSUs (flight lines). Figure 1. View largeDownload slide Outline of POP2 with four administrative units (AU1, …, AU4). A first-stage sample of 66 equally spaced PSUs (flight lines) is shown. Note that the actual S1 sample size is 131. Figure 1. View largeDownload slide Outline of POP2 with four administrative units (AU1, …, AU4). A first-stage sample of 66 equally spaced PSUs (flight lines) is shown. Note that the actual S1 sample size is 131. The sampling design is, for logistical reasons, a systematic selection of 131 equally spaced PSUs (i.e., 4%) in S1 followed by a systematic S2 sampling within each S1-selected PSU (S2 sample fraction, 1:60). A minimum of one element was selected in five PSUs with <60 elements. Simulated two-stage sampling was done for all 24 possible S1 samples and 30 replications of the S2 sampling for each S1 sample. The limitation to 30 replications was dictated by the time to complete computations. The average number of S2 elements selected was 8.5 per PSU (i.e., 1,114 elements per sample). For a S1 PSU with a single S2 element, the within-PSU estimate of the sampling variance of Y (sϕi2) was computed as the squared difference between the single observed Y value and the S2 average of yij∋ϕ. POP2 is simulated as a mosaic of 178,807 polygons of variable size. Each polygon, with the elements in it, was assigned to one of the eight classes by a random draw from a multinomial distribution. Polygons vary in size from 0.3 to 16 ha with a mean of 3.4 ha (median = 3 ha, SD = 1.8 ha). The average per hectare value of vol in a forested polygon was determined by its class, a randomly assigned age and species composition (pine, spruce, and mixed wood), and a randomly assigned canopy closure (CL, 0.63≤ CL ≤ 1.0). Species and class-specific logistic functions with vol as the dependent variable and age as the independent variable were used to assign a volume to each forested polygon. Assigned polygon volumes were then multiplied by CL. A within-polygon random variation of approximately 22% in vol among polygon elements was simulated by adding a random gamma-distributed deviation to each element in a forest polygon. The simulated element-level volume values varied from 5.0 to 1,058 m3 ha−1 with a mean of 476 m3 ha−1 and an SD of 263 m3 ha−1. There are no global trends superimposed on this assignment scheme. Additional details are available from the corresponding author on request. The two auxiliary element-level variables (X1, X2) were simulated as examples of two proxies of canopy height (mean and dominant) recorded with an airborne laser scanner and correlated with vol (Næsset 2002). In a first step, the mean (Z1) and dominant height (Z2) of trees in a polygon were determined from species- and class- specific logistic functions with age as a predictor. Then (X1, X2) was obtained by adding to (Z1, Z2) a random bias term drawn from a bivariate PERT distribution (a scaled beta distribution; Fazar 1959). Element-level values of X1 varied from 0.2 to 33.2 m with a mean and SD of 12.3 and 6.3 m, respectively. Corresponding results for X2 were 5 m (minimum), 36.0 m (maximum), 14.3 m (mean), and 7.3 m (SD). Across all classes the correlation between Y and each column of X was approximately 0.95. Within a class, it dropped to approximately 0.8. A summary of the main features of the simulations leading to POP2 is provided in the form of a diagram in Figure 2. Examples of the spatial distribution of cover classes, age, and stem volume are provided in Figure 3 for a block of 120 PSUs. Additional details and POP2 data are available from the corresponding author on request. Figure 2. View largeDownload slide A diagram of programming steps to create POP2. To avoid clutter, the addition of stochastic variation is not shown. Centers are the location centroid of stands. See text for details on symbols. Figure 2. View largeDownload slide A diagram of programming steps to create POP2. To avoid clutter, the addition of stochastic variation is not shown. Centers are the location centroid of stands. See text for details on symbols. Figure 3. View large Download slide A block of 120 PSUs indicating element values of cover class (top), age (center), and volume (below). Cover class: , pine; , spruce; , mixed; , nonforest. Values of age and volume are displayed as a gray scale (white, minimum; black, maximum). The block is located in AU4 (cf. Figure 1). Figure 3. View large Download slide A block of 120 PSUs indicating element values of cover class (top), age (center), and volume (below). Cover class: , pine; , spruce; , mixed; , nonforest. Values of age and volume are displayed as a gray scale (white, minimum; black, maximum). The block is located in AU4 (cf. Figure 1). Preliminary regression analyses with various transforms of Y and X (log, square root, and Box-Cox) and a sample size of 1,114 suggested the following linear relationship between and X  The model in Equation 10 only lists terms that were statistically significant at the 5% level in ≥50 of 1,000 fitted models with forward selection of variables (Draper and Smith 1981, p. 308). Across all samples, the model in Equation 10 explained approximately 93% of the variation in . Hedmark County Survey A two-stage survey of Hedmark County, Norway, was conducted in 2005–2007 for the purpose of linking airborne laser measurements of vegetation with field data to estimate aboveground forest tree biomass (AGB) (Mg ha−1). A two-way cross-classification by administrative unit (AU) and cover type (H) of the county was conducted with the aim of producing estimates of AGB for individual classes as well as the entire county. Hedmark County is located in southeastern Norway bordering Sweden. Its land area is approximately 27,390 km2 and comprises of four AUs (AU1, …, AU4) and eight cover types. This study uses data from six forest cover types (H1, …, H6) accounting for 67.4% of the county area. Additional details can be found in Gobakken et al. (2012). The county was partitioned into 636 parallel east-west-oriented PSUs centered on a 3 km × 3 km national forest inventory grid (Gobakken et al. 2012). Each PSU was 500 m wide and tessellated into 250 m2 square elements. In S1, a random-start systematic selection of 1 of every 12 consecutive PSUs (8.3%) was performed, and a suite of auxiliary variables (X) was obtained from an airborne laser scanner, topographic maps, and Landsat Thematic Mapper data (Gobakken et al. 2012). PSU sizes varied from 4,947 elements (123.7 ha) to 221,823 elements (5,545.6 ha) with a median of 104,405 (2,610.1 ha) and an SD of 49,997 (1,249.9 ha). The S2 sample consisted of 577 elements located in the intersection of a national 3 km × 3 km grid and the S1-selected PSUs (Gregoire et al. 2010). Each element in S2 belonged to one of the six forest cover types. An estimate of AGB was available for each S2 element, 533 from the National Forest Inventory of Norway and 44 from the Norwegian University of Life Sciences (Gobakken et al. 2012). The number of S2 units per PSU varied from 1 (twice) to 22 (once) with a median of 11 (mean = 11.1) and an SD of 5.0. As was done in POP1, with just a single S2 element in an S1 PSU, the estimate of the within-PSU sampling variance of Y (sϕi2) was estimated as the squared difference between the single observed Y value and the second-stage average of yij∋ϕ. The following regression model was used to estimate AGB from the auxiliary variables   where TF1 is a first-return (ALS) metric (proportion) of canopy density at the 10th percentile height of first-return pulses, MAXs is the maximum canopy height in last-return pulses, TS5 is canopy density at the 50th percentile of the last-pulse returns, is the mean canopy height of first returns, PFq is the q percentile canopy height in first returns, ELEV is elevation above sea level (meters), and TM5 is the top-of-atmosphere corrected (Zhu and Woodcock 2012) and normalized (Canty et al. 2004) reflectance values in Landsat Thematic Mapper band 5 in a June 3 and 10, 2007, mosaic of the county. The relative RMSE of the estimated model was estimated to be 11.6% of the average . The coefficient of determination (2) was 0.92. Note the absence of AU-specific terms in 11. The sum of second stage residuals in each AU was approximately 0 (<0.1) and nonsignificant. In addition to county, AU-, and class-specific estimates of AGB, the analyst also wants to know what the sampling variance would have been without any auxiliary variables. We therefore obtained two-stage SRSwor estimates of AGB and sampling variance. The per unit area (hectare) estimator of μy was   as per Cochran (1977, p. 303, 11.21), for example. Two variance estimators for y, S2ϕ are considered: one that corresponds to the two-stage with regression estimator in Equation 3 and one that corresponds to the proposed ratio estimator RT(RTϕ) in Equation 6. The former is (Cochran 1977, p. 303, 11.24)   and the latter is (Cochran 1977, p. 305, 11.30)   where iϕ is the S2 domain average of yij∋ϕ in the ith PSU, = (Σi=1lΣmij=1yijδij∈ϕ) × (Σi=1lΣmij=1δij∈ϕ)−1 is the PSU size weighted per element average of yij in domain ϕ, and vari|j(zij) is the variance of zij within the ith PSU. Results POP1, Equal PSU Sizes In two-stage SRSwor designs with equal PSU sizes, the estimates of y, R2 and y, RT were identical and nearly unbiased. Estimates of (absolute) bias were unimportant in all settings (maximum 0.3%, mean 0.1%). Given the low level of bias, estimates of variance and RMSEs were very similar. Variances of y, R2 and y, RT were nearly identical and strongly correlated (0.985). They varied from 0.37 to 1.57 with an average over all populations and settings of 0.77. The average values over all designs and population configurations of R2(y, R2)0.5, R1(y, R2)0.5, RT(y, RT)0.5, and GS (y, RT)0.5 were 0.83, 0.74, 0.84, and 0.68, respectively. That is an overestimation of 9 and 11% for R2(y, R2)0.5 and RT(y, RT)0.5 and an underestimation of 3 and 11% for R1(y, R2)0.5 and GS (y, RT)0.5. Across design settings, the four analytical estimators and their empirical counterparts were strongly correlated (>0.98). Yet both R2(y, R2)0.5 and RT(y, RT) became increasingly conservative as the empirical sampling variance of y, R2 and y, RT increased. Conversely, the underestimation in GS(y, RT)0.5 became less pronounced in design settings with a larger empirical variance. Replication variances of estimated analytical variances conformed to the expectation (Casella and Berger 2002, p. 60): var(σ2) ∝ σ4 (correlation ≥0.97); yet the coefficient of variation in the two MA variance estimators was approximately three to four times larger than that in the two MD variance estimators. Coverage of nominal 95% confidence intervals was best for the proposed ratio-to-size estimator of variance (minimum = 0.94, maximum = 0.97, mean = 0.96), second best for R1(y, R2)0.5 (minimum = 0.88, maximum = 0.94, mean = 0.92), and worst for R2(y, R2)0.5 (minimum = 0.70, maximum = 0.96, mean = 0.91) and GS(y, RT)0.5 (minimum = 0.78, maximum = 0.95, mean = 0.91). As expected, coverage improved with increasing S1 sample size. The marked difference in coverage between RT(y, RT)0.5 and R2(y, R2)0.5 can be explained by the larger coefficient of variation in the latter. POP1, Unequal PSU Sizes Estimates of (absolute) bias for designs with unequal PSU sizes were practically identical to those obtained with equal PSU sizes. Unequal PSU sizes inflated the replicate variance of y, R2 and associated estimators of variance R2(y, R2)0.5 and R1(y, R2)0.5 by an average factor of 5.7 (minimum = 4.3, maximum = 7.5). The inflation of R2(y, R2)0.5 can be traced to the unequal scaling of within-PSU regression residuals in Equation 3. The inflation was strongest in designs with the smallest S2 sample size and weakest in designs with the largest S2 sample size. Figure 4 illustrates the marked effect of unequal PSU sizes on estimates of R2(y, R2)0.5. Figure 4. View large Download slide Scatter plots of the mean of analytical estimates of the SEs of y,R2 and y,RT plotted against their empirical (Monte Carlo) estimates. , R2(y,R2)0.5; #, GS(y,RT)0.5; *, RT (y,RT)0.5. Figure 4. View large Download slide Scatter plots of the mean of analytical estimates of the SEs of y,R2 and y,RT plotted against their empirical (Monte Carlo) estimates. , R2(y,R2)0.5; #, GS(y,RT)0.5; *, RT (y,RT)0.5. In contrast, the replicate variance of y, RT was barely affected by a variation in PSU sizes. Replication averages of RT(y, R2)0.5 were similar (minimum = 0.37, mean = 0.85, maximum = 1.56) to those reported for equal PSU sizes. GS(y, RT)0.5 was, with an overall average of 0.99, about 1.5 times larger than when PSU sizes were constant. Replicate averages of the analytical variance estimators matched their empirical counterparts to within 1%. Thus, the inflation in var(y, R2) due to unequal PSU sizes was correctly captured by the analytical estimator. For RT(y, R2)0.5, the coefficient of variation remained almost unaffected by unequal PSU sizes; it was approximately 40% lower than that for the three alternative estimators. Coverage of nominal 95% confidence intervals computed with the proposed estimator of variance (Equation 6) was virtually equal to the coverage in scenarios with equal PSU sizes. Coverage achieved with R2(y, R2) was now slightly better than that in scenarios with equal PSU sizes and matched the coverage achieved with R1 (y, R2). The inflation of R2(y, R2) in response to unequal PSU sizes explains the improvement. Coverage computed with GS (y, R2) was, on average, 0.98 (minimum = 0.96, maximum = 0.99), which marks a switch from a 0.04 undercoverage in the case of an equal-sized PSU. POP2 Estimators y, R2 and y, RT of the average stem volume (vol, m3 ha−1) in the five forested classes appear to be nearly unbiased (bias ≈ ±0.2%) (Table 1). The bias of y, RT of single domains was approximately twice as large as that for all forested classes combined, but only approximately half as large as the apparent bias of y, R2. The bias of y, RT was larger for smaller domains. However, the bias of y, R2 was seemingly not related to the size of a domain. The larger domain-specific bias of y, R2 is mostly due to the π expansion of the sum of within-PSU regression residuals by the ratio of the nominal PSU size to the nominal size of the within-PSU S2 sample (Equation 2). A scaling based on the ratio of the actual number of S1 to S2 domain-specific elements would have reduced the bias to the level in y, RT (not shown). Parameters and estimation results for POP2 domains. Table 1. Parameters and estimation results for POP2 domains. See text for definitions of statistics. *Equation(s) used to derive a table entry. †μy = mean per hectare stem volume (m3 × ha−1), V(y)0.5 = SD of volume. View Large Table 1. Parameters and estimation results for POP2 domains. See text for definitions of statistics. *Equation(s) used to derive a table entry. †μy = mean per hectare stem volume (m3 × ha−1), V(y)0.5 = SD of volume. View Large For the combined classes (H1, …, H5), the RMSE of y, R2 was, on average, 1.4 times larger than the RMSE of y, RT (Table 1). At the level of AUs, the RMSEs of y, R2 were, on average, 2.6 times larger than the RMSEs of y, RT. In forest classes, the RMSEs of y, R2 were, on average, 5.8 times larger than that for y, RT. For the proposed variance estimator (Equation 6), the average estimate of the SEM vol in the five forested classes [(y, RT)0.5] was only 3% above the corresponding empirical estimate (Table 1). In addition, nominal 95% confidence intervals achieved coverage of 0.94 (Table 1, row 16). Results with the alternative MD variance estimator (Equation 7) (Table 1, row 13) suggest an overestimation of approximately 18% in estimates of SEs and a corresponding positive bias of 0.02 in coverage (Table 1, row 17). In contrast, analytical SEs of y, R1 and (Table 1, POP2 column in rows 10 and 11) were approximately 5.3 times larger than their empirical counterparts. The poor performance of R2(y, R2)0.5 and R1(y, R2)0.5 is, for the most part, a consequence of a sensitivity to the π scaling of regression residuals. AU-specific ratio-to-size estimates of SE RT(y, RT)0.5 were, on average, 38% above the empirical estimates, and coverage of 95% confidence intervals was poor (Table 1, rows 9, 12, and 16). SEs estimated by GS(y, R2)0.5 were on average, 153% larger than the empirical estimates; computed confidence intervals were correspondingly conservative (coverage = 1). Results with the two-stage MA estimators (Table 1 rows 10 and 11) were considerably worse. Class-specific analytical estimates of SE with RT(y, RT)0.5 were, on average, 4.7% below corresponding empirical estimates. Results with GS(y, RT)0.5 were similar. In contrast, R2(y, R2)0.5 and R1 (y, R2)0.5 overestimated the empirical SE by 103 and 79%, respectively (Table 1, rows 10 and 11); coverage of 95% confidence intervals was correspondingly poor (Table 1, rows 14–17). Hedmark County Survey A summary of the Hedmark County results is presented in Table 2. With the admonition that interpretation of nonreplicated results requires caution, a few results deserve attention, whereas others appear to confirm trends seen in the simulations. Parameters and estimation results from the Hedmark County survey. Table 2. Parameters and estimation results from the Hedmark County survey. See text for definitions of statistics. *Equation(s) used to compute, in whole or partially, a table entry. View Large Table 2. Parameters and estimation results from the Hedmark County survey. See text for definitions of statistics. *Equation(s) used to compute, in whole or partially, a table entry. View Large The ratio-to-size estimates of per hectare of AGB (y, RT) were, when compared with the MA estimates (y, R2), considerably closer to estimates that would have been obtained from a two-stage design without auxiliary information (y, S2). The importance of choosing an appropriate variance estimator for two-stage equal probability designs with unequal PSU sizes was reiterated by the fact thatS2(y, S2)0.5 was, on average, approximately 69% larger than RT2(y, S2)0.5 (Table 2). Estimates of analytical MD SEs [RT(y, RT)0.5 and GS (y, RT)0.5] were consistently smaller than corresponding estimates from a two-stage design without auxiliary information [RT2 (y, S2)0.5]. Overall, it appears that the auxiliary information has lowered the SE by approximately 50%. No important difference emerged between results obtained with the proposed ratio-to-size variance estimator [RT(y, RT)] and those obtained with the estimator proposed by Ståhl et al. (2011), i.e., GS(y, RT). Two-stage with regression estimates of (analytical) SEs (R20.5) were often larger than estimates S20.5 obtained without the use of auxiliary information and, consequently, were well above MD estimates of error. Although the results with the MA estimators are counterintuitive, they vividly confirm a lack of precision in this particular application. Discussion The prevalence of a multitude of easy, accessible, remotely sensed auxiliary information sources has transformed large-scale natural resource surveys (Bauer et al. 1994, McRoberts et al. 2006). The simplest way to use the auxiliary information is for prestratification (Stehman et al. 2003, Tomppo et al. 2008) and poststratification (McRoberts et al. 2012). When the auxiliary variable is correlated with the variable(s) of interest, it can be used in single-stage regression, ratio, and calibrated regression estimators (Corona and Fattorini 2008, Stehman 2009, Næsset et al. 2011). Alternatively, the auxiliary variable(s) may be an integral part of the sampling design as in double sampling (Parker and Evans 2004, Opsomer et al. 2007, Mandallaz 2008, p. 173, remark 5, Saborowski et al. 2010) and multistage (phase) designs (Schreuder et al. 1995, Baffetta et al. 2007, Fattorini et al. 2009). Textbooks on sampling statistics contain appropriate MA estimators of totals and variance (e.g., Särndal et al. 1992, Köhl et al. 2006, Gregoire and Valentine 2008, Fuller 2009). In large, regional natural resource surveys wall-to-wall coverage with auxiliary variables may not be economically or technically feasible nor may it be needed to achieve a target precision. In those cases, a sample design for the auxiliary variables is needed, whereby the acquisition mode and spatial resolution of the auxiliary variable(s) may define both sampling units and their location. A combination of logistic, technical, and economic constraints paired with the geographic outline of the target population, may lead to a design with PSUs of unequal size, yet equal probability of sample inclusion. In this case, two-stage MA estimators can be imprecise (Ene et al. 2012, Næsset et al. 2013). Two-phase MA estimators appear to be more promising (Ene et al. 2013, Mandallaz 2013). The issue of unequal sized PSUs is generally resolved by switching to a ratio-to-size type of estimator (Cochran 1977, p. 249 and 303). With large sample sizes, as is generally the case in regional surveys, the bias in a ratio estimator is unimportant (Cochran 1977, p. 160). To our knowledge, there does not seem to be a clearcut choice of a MA estimator for two-stage designs suited to handle the issue of equal-probability PSUs with unequal sizes. If one is willing to drop the invariance requirement (Särndal et al. 1992, p. 134), one may revert to a MA two-phase estimator (Estevao and Särndal 2009). Results made public after the submission of this study (Ene et al. 2013, Mandallaz 2013) suggest that the simpler MA two-phase estimator can be efficient. In the MD ratio-to-size variance estimators, the second-stage design is ignored. Dispensing with the details of the S2 sampling design makes MD estimators more flexible (Little 2004). An MD estimator of variance should therefore be expected to perform relatively well when the efficiency of a design-based estimator is compromised (Särndal et al. 1992, p. 516). The proposed ratio estimator is flexible and applies to a suite of regional or national surveys. The proposed ratio-to-size variance estimator appears to be an improvement over the estimator proposed by Ståhl et al. (2011b). Its performance in simulated sampling was better and more consistent across designs with equal and unequal PSU sizes. The difference between the proposed ratio-to-size estimator of variance and that of Ståhl et al. (2011) arises from the following three terms included in the former but not the latter: finite population correction factors; a variance effective number of PSUs; and the numerator term . The first two terms reflect our choices, and the last term captures the effects of sampling x. An analyst has good arguments for preferring design-based estimators over MD estimators (Särndal et al. 1992, p. 516, Gregoire 1998, Little 2004). However, reliance on attractive theoretical properties of design-based estimators, when the analyst is faced with the prospect of imprecise, inefficient, and even counterintuitive two-stage MA estimates, for reasons beyond his or her control, is a risky proposition (Valliant et al. 2000, p. 14). In these situations, we recommend ratio-to-size estimator(s). In large, two-stage natural resource surveys with a sufficiently large probability sample of Y, the properties of models used in MD and MA estimators should, however, be similar and qualify as approximately unbiased (Hansen et al. 1983). A drawback of ratio estimators is that they are not additive (Särndal et al. 1992, p. 403). When, as here, a PSU can straddle two or more domains, neither ratio- nor design-based estimators are additive. The latter lack the additive property because the nominal design-based π scaling applied to individual domains is biased. The reason is that expected realized sample sizes in stage one and two do not match the nominal sample sizes. In their analysis of the Hedmark County survey, Ståhl et al. (2011) obtained ratio-to-size estimates for separate domains (administrative units and cover types) and then combined these estimates and their variances for an estimate for the county (Ståhl et al. 2011b). Our direct approach to estimating a population level per element average leads to more efficient estimators, albeit at the loss of additivity of individual domain estimates. In POP2 and in the Hedmark County survey, the selection of stage one units and stage two elements was systematic. The application of SRSwor estimators is pragmatic. Results from POP2 suggest that ratio-to-size estimators performed better than a two-stage MA estimator. Still, alternative estimators for systematic sampling (Sherman 1996, D'Orazio 2003, Ekström and Sjöstedt-de Luna 2004) or resampling estimators (Sitter 1992, Shao 1996, 2010, Lahiri 2003a, 2003b) deserve exploration. Conclusions In large-scale natural resource probability sampling, the use of remotely sensed auxiliary information has significantly improved precision and efficiency. At times, however, various constraints (e.g., costs, time, logistics, and technical issues), as well as the physical structure and composition of a population, can lead to designs for which design-based estimators may be imprecise. Faced with this prospect, alternative estimators, such as the two-phase with regression (Mandallaz 2013) or first-difference estimator (Nelson et al. 2008, Fuller 2009, p. 306), should be explored. For two-stage SRSwor designs with regression and PSUs of unequal size, ratio-to-size estimators appear attractive. Acknowledgments: We thank the associate editor and three anonymous reviewers for many helpful suggestions and improvements to earlier drafts of this article. We also acknowledge Wayne Fuller for a helpful discussion of the proposed variance estimator. Literature Cited Baffetta F. Bacaro G. Fattorini L. Rocchini D. Chiarucci A. . 2007. Multi-stage cluster sampling for estimating average species richness at different spatial grains. Community Ecol . 8( 1): 119– 127. Google Scholar CrossRef Search ADS   Bauer M.E. Burk T.E. Ek A.R. Coppin R.R.S.D. Lime S.D. Walsh T.A. Walters D.K. Befort W. Heinzen D.F. . 1994. Satellite inventory of Minnesota forest resources. Photogr. Eng. Remote Sens . 60( 3): 287– 298. Breidenbach J. Næsset E. Lien V. Gobakken T. Solberg S. . 2010. Prediction of species specific forest inventory attributes using a nonparametric semi-individual tree crown approach based on fused airborne laser scanning and multispectral data. Remote Sens. Environ . 114( 4): 911– 924. Google Scholar CrossRef Search ADS   Canty M.J. Nielsen A.A. Schmidt M. . 2004. Automatic radiometric normalization of multitemporal satellite imagery. Remote Sens. Environ . 91( 3–4): 441– 451. Google Scholar CrossRef Search ADS   Casella G. Berger R.L. . 2002. Statistical inference . Duxbury Press, Pacific Grove, CA. 660 p. Cochran W.G. 1977. Sampling techniques . John Wiley & Sons, New York. 380 p. Coomes D.A. Allen R.B. Scott N.A. Goulding C. Beets P. . 2002. Designing systems to monitor carbon stocks in forests and shrublands. For. Ecol. Manage . 164( 1–3): 89– 108. Google Scholar CrossRef Search ADS   Corona P. Fattorini L. . 2008. Area-based lidar-assisted estimation of forest standing volume. Can. J. For. Res . 38( 11): 2911– 2916. Google Scholar CrossRef Search ADS   D'Orazio M. 2003. Estimating the variance of the sample mean in two-dimensional systematic sampling. J. Agr. Biol. Environ. Stat . 8( 3): 280– 295. Google Scholar CrossRef Search ADS   Dalponte M. Bruzzone L. Gianelle D. . 2012. Tree species classification in the Southern Alps based on the fusion of very high geometrical resolution multispectral/hyperspectral images and LiDAR data. Remote Sens. Environ . 123: 258– 270. Google Scholar CrossRef Search ADS   Draper N.R. Smith H. . 1981. Applied regression analysis . John Wiley & Sons, New York. 709 p. Ekström M. Sjöstedt-de Luna S. . 2004. Subsampling methods to estimate the variance of sample means based on non-stationary spatial data with varying expected values. J. Am. Stat. Assoc . 99( 465): 82– 95. Google Scholar CrossRef Search ADS   Ene L.T. Næsset E. Gobakken T. Gregoire T.G. Ståhl G. Holm S. . 2013. A simulation approach for accuracy assessment of two-phase post-stratified estimation in large-area LiDAR biomass surveys. Remote Sens. Environ . 133: 210– 224. Google Scholar CrossRef Search ADS   Ene L.T. Næsset E. Gobakken T. Gregoire T.G. Ståhl G. Nelson R. . 2012. Assessing the accuracy of regional LiDAR-based biomass estimation using a simulation approach. Remote Sens. Environ . 123: 579– 592. Google Scholar CrossRef Search ADS   Estevao V.M. Särndal C.E. . 2009. A new face on two-phase sampling with calibration estimators. Surv. Methods  35( 1): 3. Faes C. Molenberghs G. Aerts M. Verbeke G. Kenward M.G. . 2009. The effective sample size and an alternative small-sample degrees-of-freedom method. Am. Stat . 63( 4): 389– 399. Google Scholar CrossRef Search ADS   Fattorini L. Franceschi S. Pisani C. . 2009. A two-phase sampling strategy for large-scale forest carbon budgets. J. Stat. Plan. Infer . 139: 1045– 1055. Google Scholar CrossRef Search ADS   Fazar W. 1959. Program evaluation and review technique. Am. Stat . 13( 2): 10– 16. Fuller W.A. 2009. Sampling statistics . John Wiley & Sons, New York. 454 p. Google Scholar CrossRef Search ADS   Gasparini P. Gregori E. Pompei E. Rodeghiero M. . 2010. The Italian National forest inventory: Survey methods for carbon pools assessment. Sherwood For. Alberi Oggi  ( 168): 13– 18. Ghosh M. 2008. Robust estimation in finite populations. P. 116– 122 in Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen . International Mathematics Society, Beachwood, OH. Gobakken T. Næsset E. Nelson R. Bollandsås O.M. Gregoire T.G. Ståhl G. Holm S. Ørka H.O. Astrup R. . 2012. Estimating biomass in Hedmark County, Norway using national forest inventory field plots and airborne laser scanning. Remote Sens. Environ . 123: 443– 456. Google Scholar CrossRef Search ADS   Gregoire T.G. 1998. Design-based and model-based inference in survey sampling: Appreciating the difference. Can. J. For. Res . 28: 1429– 1447. Google Scholar CrossRef Search ADS   Gregoire T.G. Lin Q.F. Boudreau J. Nelson R. . 2008. Regression estimation following the square-root transformation of the response. For. Sci . 54( 6): 597– 606. Gregoire T.G. Ståhl G. Næsset E. Gobakken T. Nelson R. Holm S. . 2010. Model-assisted estimation of biomass in a LiDAR sample survey in Hedmark County, Norway. Can. J. For. Res . 41( 1): 83– 95. Google Scholar CrossRef Search ADS   Gregoire T.G. Valentine H.T. . 2008. Sampling strategies for natural resources and the environment . Chapman & Hall/CRC, Boca Raton, FL. 465 p. Hansen M.H. Madow W.G. Tepping B.J. . 1983. An evaluation of model-dependent and probability-sampling inferences in sample surveys. J. Am. Stat. Assoc . 78( 384): 776– 793. Google Scholar CrossRef Search ADS   Hudak A.T. Evans J.S. Stuart Smith A.M. . 2009. LiDAR utility for natural resource managers. Remote Sens . 1( 4): 934– 951. Google Scholar CrossRef Search ADS   Hyyppä J. Hyyppä H. Leckie D. Gougeon F. Yu X. Maltamo M. . 2008. Review of methods of small-footprint airborne laser scanning for extracting forest inventory data in boreal forests. Int. J. Remote Sens . 29( 5): 1339– 1366. Google Scholar CrossRef Search ADS   Koehler E. Brown E. Haneuse J.-P.A. . 2009. On the assessment of Monte Carlo error in simulation-based statistical analyses. Am. Stat . 63( 2): 155– 162. Google Scholar CrossRef Search ADS PubMed  Köhl M. Magnussen S. Marchetti M. . 2006. Sampling methods, remote sensing and GIS multiresource forest inventory . Springer, Berlin, Germany. 374 p. Google Scholar CrossRef Search ADS   Lahiri P. Rao J.N.K. . 1995. Robust estimation of mean squared error of small area estimators. J. Am. Stat. Assoc . 90( 430): 758– 766. Google Scholar CrossRef Search ADS   Lahiri R. 2003a. On the impact of bootstrap in survey sampling and small-area estimation. Stat. Sci . 18( 2): 199– 210. Google Scholar CrossRef Search ADS   Lahiri S.N. 2003b. Resampling methods for dependent data . Springer, New York. 374 p. Google Scholar CrossRef Search ADS   Li Y. Lahiri P. . 2007. Robust model-based and model-assisted predictors of the finite population total. J. Am. Stat. Assoc . 102( 478): 664– 673. Google Scholar CrossRef Search ADS   Little R.J.A. 2004. To model or not to model? Competing modes of inference for finite population sampling. J. Am. Stat. Assoc . 99( 466): 546– 556. Google Scholar CrossRef Search ADS   Loetsch F. Haller K.E. . 1973. Forest inventory . BLV-Verlagsgesellschaft, Munich, Germany. 436 p. Magnussen S. Næsset E. Gobakken T. Frazer G. . 2012. A fine-scale model for area-based predictions of tree-size-related attributes derived from LiDAR canopy heights. Scand. J. For. Res . 27( 3): 1– 11. Google Scholar CrossRef Search ADS   Mandallaz D. 2008. Sampling techniques for forest inventories . Chapman & Hall, Boca Raton, FL. 251 p. Mandallaz D. 2013. Design-based properties of some small-area estimators in forest inventory with two-phase sampling. Can. J. For. Res . doi:10.1139/cjfr-2012-0381. McRoberts R.E. Cohen W.B. Næsset E. Stehman S.V. Tomppo E.O. . 2010a. Using remotely sensed data to construct and assess forest attribute maps and related spatial products. Scand. J. For. Res . 25( 4): 340– 367. Google Scholar CrossRef Search ADS   McRoberts R.E. Gobakken T. Næsset E. . 2012. Post-stratified estimation of forest area and growing stock volume using lidar-based stratifications. Remote Sens. Environ . 125: 157– 166. Google Scholar CrossRef Search ADS   McRoberts R.E. Holden G.R. Nelson M.D. Liknes G.C. Gormanson D.D. . 2006. Using satellite imagery as ancillary data for increasing the precision of estimates for the Forest Inventory and Analysis program of USDA Forest Service. Can. J. For. Res . 36: 2968– 2980. McRoberts R.E. Tomppo E.O. Næsset E. . 2010b. Advances and emerging issues in national forest inventories. Scand. J. For. Res . 25( 4): 368– 381. Google Scholar CrossRef Search ADS   Næsset E. 2002. Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data. Remote Sens. Environ . 80( 1): 88– 99. Google Scholar CrossRef Search ADS   Næsset E. 2004. Accuracy of forest inventory using airborne laser scanning: Evaluating the first Nordic full-scale operational project. Scand. J. For. Res . 19( 6): 554– 557. Google Scholar CrossRef Search ADS   Næsset E. Gobakken T. Bollandsås O.M. Gregoire T.G. Nelson R. St.åhl G. . 2013. Comparison of precision of biomass estimates in regional field sample surveys and airborne LiDAR-assisted surveys in Hedmark County, Norway. Remote Sens. Environ . 130: 108– 120. Google Scholar CrossRef Search ADS   Næsset E. Gobakken T. Solberg S. Gregoire T.G. Nelson R. Ståhl G. Weydahl D. . 2011. Model-assisted regional forest biomass estimation using LiDAR and InSAR as auxiliary data: A case study from a boreal forest area. Remote Sens. Environ . 115( 12): 3599– 3614. Google Scholar CrossRef Search ADS   Nelson R. Krabill W. Tonelli J. . 1988. Estimating forest biomass and volume using airborne laser data. Remote Sens. Environ . 24: 247– 267. Google Scholar CrossRef Search ADS   Opsomer J.D. Breidt F.J. Moisen G.G. Kauermann G. . 2007. Model-assisted estimation of forest resources with generalized additive models. J. Am. Stat. Assoc . 102( 478): 400– 409. Google Scholar CrossRef Search ADS   Parker R.C. Evans D.L. . 2004. An application of LiDAR in a double-sample forest inventory. W. J. Appl. For . 19( 2): 95– 101. Rao P.S.R.S. 1988. Ratio and regression estimators. P. 449– 468 in Handbook of Statistics, vol. 6: Sampling , Krishnaiah P.R. Rao C.R. (eds.). Elsevier, Amsterdam, The Netherlands. Google Scholar CrossRef Search ADS   Royall R.M. 1976. The linear least-squares prediction approach to two-stage sampling. J. Am. Stat. Assoc . 71( 355): 657– 664. Google Scholar CrossRef Search ADS   Saborowski J. Marx A. Nagel J. Böckmann T. . 2010. Double sampling for stratification in periodic inventories—Infinite population approach. For. Ecol. Manage . 260( 10): 1886– 1895. Google Scholar CrossRef Search ADS   Särndal C.E. Swensson B. Wretman J. . 1992. Model assisted survey sampling . Springer, New York. 694 p. Google Scholar CrossRef Search ADS   Satterthwaite F.E. 1946. An approximate distribution of estimates of variance components. Biometrics Bull . 2: 110– 114. Google Scholar CrossRef Search ADS   Schreuder H.T. LaBau V.J. Hazard J.W. . 1995. The Alaska four-phase forest inventory sampling design using remote sensing and ground sampling. Photogram. Eng. Remote Sens . 61( 3): 291– 297. Shao J. 1996. Resampling methods in sample surveys. Statistics  27: 203– 254. Google Scholar CrossRef Search ADS   Shao X. 2010. The dependent wild bootstrap. J. Am. Stat. Assoc . 105( 489): 218– 235. Google Scholar CrossRef Search ADS   Sherman M. 1996. Variance estimation for statistics computed from spatial lattice data. J. R. Stat. Soc. Ser. B . 58( 3): 509– 523. Sitter R.R. 1992. A resampling procedure for complex survey data. J. Am. Stat. Assoc . 87( 419): 10. Google Scholar CrossRef Search ADS   Spencer R.D. Czaplewski R.L. . 1997. National forest inventory in the USA: An outline of the procedure. Aust. For . 60( 1): 56– 66. Google Scholar CrossRef Search ADS   Spurr S.H. 1952. Forest inventory . Ronald Press, New York. 476 p. Ståhl G. Allard A. Esseen P.A. Glimskär A. Ringvall A. Svensson J. Sundquist S. et al.  . 2011a. National Inventory of Landscapes in Sweden (NILS)-scope, design, and experiences from establishing a multiscale biodiversity monitoring system. Environ. Monit. Assess . 173( 1–4): 579– 595. Google Scholar CrossRef Search ADS   Ståhl G. Holm S. Gregoire T.G. Gobakken T. Næsset E. Nelson R. . 2011b. Model-based inference for biomass estimation in a LiDAR sample survey in Hedmark County, Norway. Can. J. For. Res . 41( 1): 96– 107. Google Scholar CrossRef Search ADS   Stanek E.J. Singer J.M. . 2004. Predicting random effects from finite population clustered samples with response error. J. Am. Stat. Assoc . 99( 468): 1119– 1130. Google Scholar CrossRef Search ADS   Stehman S.V. 2009. Model-assisted estimation as a unifying framework for estimating the area of land cover and land-cover change from remote sensing. Remote Sens. Environ . 113( 11): 2455– 2462. Google Scholar CrossRef Search ADS   Stehman S.V. Sohl T.L. Loveland T.R. . 2003. Statistical sampling to characterize recent United States land-cover change. Remote Sens. Environ . 86: 517– 529. Google Scholar CrossRef Search ADS   Stevens D.L. Olsen A.R. . 2004. Spatially balanced sampling of natural resources. J. Am. Stat. Assoc . 99( 465): 262– 278. Google Scholar CrossRef Search ADS   Thompson S.K. 1992. Sampling . John Wiley & Sons, New York. 343 p. Tomppo E. 1991. Satellite image-based national forest inventory of Finland. P. 419– 424 in Proc. of the Symposium on global and environmental monitoring, techniques and impacts. International Society for Photogrammetry and Remote Sensing, Victoria, BC, Canada. Tomppo E. Olsson H. Stahl G. Nilsson M. Hagner O. Katila M. . 2008. Combining national forest inventory field plots and remote sensing data for forest databases. Remote Sens. Environ . 112( 5): 1982– 1999. Google Scholar CrossRef Search ADS   Valliant R. Dorfman A.H. Royall R.M. . 2000. Finite population sampling and inference. A prediction approach . John Wiley & Sons, New York. 504 p. Zhu Z. Woodcock C.E. . 2012. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ . 18( 1): 83– 94. Google Scholar CrossRef Search ADS   Appendix Had we observed yij in all S1 elements, the variance estimator for RTϕ in Equation 5, given the SRSwor two-stage design, would be (Cochran 1977, p, 155, 6.10 or 6.13)   Derivations of the proposed estimator of variance in Equation 6 are given next. To simplify notation, we have dropped the reference (ϕ) to a specific domain. The first step in the derivations is an MD assumption that yij is equal to plus two independent random error terms (ε1ij, ε2ij) due to residual error (ε1ij) and model estimation error (ε2ij). It is furthermore assumed that the model is approximately unbiased for all s1 and s2 elements in each domain, i.e., the sum of both ε1ij and ε2ij over all s1 and s2 elements is approximately 0. Residual errors are defined as ε1ij ≡ yij − f (xij;β), and ε1ij is assumed to be independently and identically distributed with a mean of zero and a variance of σε12. Model estimation errors are defined as   where f (xij;β) is a least squares estimate of yij, ε is a p × 1 row vector of the differences (β − ) between the (true) finite population regression coefficients and the S2-based estimates of β, and εt is the transpose of ε. From A.1 it is clear that ε2ij depends on xij and ε, which is an s2-specific (unknown) model error. Note that the expression in A.1 does not extend without modifications to a nonlinear model. Inserting yij = + ε1ij = + ε1ij + ε2ij into A.0 and making a few additional modifications (see below) one obtains   where the expectation (E[ ]) is over repeated two-stage SRSwor samples, mi is the number of elements in the ith PSU, is the average number of elements per PSU, l is the number of S1 PSUs with at least one element, and is the variance effective equivalent of l (Satterthwaite 1946, Cochran 1977, p 96, 5.16, Faes et al. 2009). In simulations (not shown), we found that replacing l with in the denominator of A.2 lowered the chance of underestimating the variance of RT. The formula used here for was   where is the variance of the residuals in the fitted regression model and is the estimated covariance matrix of the regression coefficients in the model. The expectation operator in A.2 is necessary because we have replaced a fixed quantity yij by a sum of three random terms. Over repeated samples, the regression coefficients will vary, which makes random. To compute the expected variance in A.2, we appeal to the assumption that f(xij;β) is an unbiased estimator of f(xij;β), and we extend to + ε1ij + ε2ij the property that if yij had been observed in all sampled elements the covariance among PSU totals would be 0 by virtue of the SRS sampling in stage one. All three terms in the numerator of A.2 contain random variables. We shall now derive the expectations for all three terms. The first term in the numerator of A.2 expands to   To simplify the expression in A.4, we appeal to the assumptions of the independence of the model estimation and the residual error and the approximately unbiased property of the model to get   A parallel to A.5 is the zero covariance between model estimation errors (ε2ij) and ordinary least squares residuals (ε1ij + ε2ij) (Draper and Smith 1981, p 148). Inserting the approximation in A.5 into A.4 gives   A second appeal to the approximate unbiased property of the model allows a simplification of A.6. Specifically,   with an equation sign only if applied to ordinary least squares residuals (ε1ij + ε2ij) (Draper and Smith, 1981, p. 148, #816). With the approximations in A.5 and A.7 inserted into A.4, we get   In A.8, the difference ε1ijε1ik − ε2ijε2ik, when summed over all s1 and s2 elements, is approximately 0   The sequence of expansions and approximations from A.4 to A.9 leaves us with the term ijik summed over j and k. However, because ijik is the product of two random variables (due to errors in ), we have to account for the uncertainty in and how it adds to the variance of RTϕ. If X was fixed, the expected value of ijik is E[] E[ik] + xijtΣxij, where Σ is the estimated covariance matrix for the model parameters β. Conversely, if we fix β and let X vary from sample to sample, the expected value of ijik over all possible samples would be where Σxi is the sample covariance matrix of xij. Because we do not know the expectations, we simply replace them with our sample-based estimates, which gives   Note that a robust design-consistent and asymptotically unbiased estimator of Σβ could have been used instead of the least squares estimate (Mandallaz 2013). To derive expectations for the two remaining terms in the numerator of A.2 [i.e., (RT)2mi2 and 2miRT Σj=1mi/( + ε1ij + ε2ij)], we reapply the derivations leading from A.4 to A.10, specifically:   and   Now, if we expand RTϕ to a ratio of two sums, as in Equation 5, we have two additional terms in the numerator of A.2 that contains sums of products of two model predictions (i.e., ijik). After replacing the three random terms in the numerator of A.2 with their expectations in (cf. A.10, A.11, and A12), reinstating domain indicator(s) and subscripts and adding the final population correction factors, we arrive (after some simplification) at the proposed variance estimator in Equation 6. Copyright © 2014 Society of American Foresters
TI - An Estimator of Variance for Two-Stage Ratio Regression Estimators
JF - Forest Science
DO - 10.5849/forsci.12-163
DA - 2014-08-01
UR - https://www.deepdyve.com/lp/springer-journals/an-estimator-of-variance-for-two-stage-ratio-regression-estimators-TAdChWS0M0
SP - 663
EP - 676
VL - 60
IS - 4
DP - DeepDyve
ER -