A sample coordination method to monitor totals of environmental variables

A new sampling strategy for design‐based environmental monitoring is proposed. It has the potential to produce superior estimators for totals of environmental variables and their changes over time. In the strategy, we combine two concepts known as spatially balanced sampling and coordination of samples. Spatially balanced sampling can provide superior estimators of totals, while coordination of samples over time is often used to improve estimators of change. Compared with reference strategies, we show that the new strategy can improve the precision of the estimators of totals and their change simultaneously. A forest inventory application is used to illustrate the new strategy and the results can be summarized as (i) using auxiliary information to spread the sample can improve the estimators of totals; (ii) a positive coordination of the samples reduced the variance of the estimator of change by more than 37% compared with independent samples; and (iii) a high overlap between successive samples does not guarantee a good estimator of change. The presented strategy can be used to develop more efficient environmental monitoring programs.


INTRODUCTION
Environmental monitoring is defined as the observation and study of the environment (Awange, 2012, ch. 1). The approach for environmental monitoring is to collect and analyze a subset that represents the environment in space and time (Artiola, Pepper, & Brusseau, 2004, ch. 2). In the whole process, sampling is usually employed as a tool to select a representative portion from the population in order to do the analysis. As an important component of environmental monitoring, it provides the foundation of data required for assessments of environmental variables.
In this article, we propose a design-based sampling strategy for monitoring totals of environmental variables. Two concepts are combined in the strategy, spatially balanced sampling and coordination of samples over time. Spatially balanced sampling designs can provide representative samples, and sample coordination is a method to statistically control the overlap of successive samples. Spatially correlated Poisson sampling (SCPS) and sample coordination based on permanent random numbers (PRNs), introduced by Brewer, Early, and Joyce (1972), are used in the algorithm for the new strategy. SCPS was first presented by Grafström (2012) as a spatial sampling method for selecting well-spread samples from finite populations. We show how sample coordination can be applied within a continuous population framework by using a double (two-phase) sampling approach. Auxiliary information is derived for a very large first-phase sample, and from the finite first-phase sample we select well-spread and positively coordinated samples using SCPS with the aid of the auxiliary variables.
Two common objectives of environmental monitoring are to characterize the current state of some resource as well as the change or trend in the state over time and space (Marker & Stevens Jr., 2009). Here, we consider only states that can be expressed as totals of environmental variables. Good estimators of such states can be achieved by using representative samples of the population at different time occasions, which means that states are generally best detected by updating the samples according to the current population. A high degree of overlap level between successive samples can produce more precise estimators of change over time (see Qualité & Tillé, 2008;Sen, 1973). Permanent samples have traditionally been used to address this issue. The samples are then often systematically distributed over the landscape (Scott, 1998).
Huge amount of financial resources are spent on environmental monitoring all over the world and it is very important to apply more efficient sampling strategies that can increase the quality of estimates while potentially saving a considerable amount of costs. McDonald (2003) provided a very detailed review of different survey designs for large-scale environmental monitoring programs, and the split panel designs proposed by Kish (1983Kish ( , 1986 were recommended among many other designs. One feature of environmental populations is that they exist in a spatial context. Commonly, the responses of nearby locations tend to be more similar than the responses of locations which are farther apart. Thus, the spatial distribution of these populations can be used as important information when designing the sample. Several efficient spatial sampling methods have recently been developed for sampling from georeferenced populations. One of the first and the most widely used method is called Generalized Random Tessellation Stratified (GRTS) design proposed by Stevens and Olsen (2004). This method uses a random mapping to map the two-dimensional locations into one dimension while preserving some spatial relationship. Then the systematic ps sampling is applied and the sample is mapped back to the two-dimensional original space. Using this design the samples can be spread evenly over the geographical space. Robertson, Brown, McDonald, and Jaksons (2013) extended the idea of GRTS to a new design called balanced acceptance sampling (BAS). To select a spatially balanced sample using BAS, we need first to specify a d-dimensional hyperrectangular box that encloses the population. For example, it could be a rectangle that encloses a two-dimensional geographical study area. Then a two-dimensional random start Halton sequence with uniformly spread points in the rectangle is generated. The first n points that are observed in the study area will constitute the sample of size n. An alternative design to BAS which is called Halton iterative partitioning (HIP) was introduced by Robertson, McDonald, Price, and Brown (2018) to overcome some drawbacks of BAS for finite populations.
Another feature of environmental population is that there are often some auxiliary information (other than geographical coordinates) available. Nowadays, there is a wealth of remotely sensed information available from satellite, aerial photography, or laser scanning that can be used to efficiently distribute the sample units. Therefore, we shall take the full advantage of the available auxiliary information and the properties of the environmental populations when we derive the sampling strategies. The local pivotal method (LPM) and SCPS proposed by Grafström (2012) and Grafström, Lundström, and Schelin (2012) are two spatial sampling designs that employ auxiliary variables (often including geographical coordinates plus several other attribute variables) to spread the samples based on distances. By adding more variables that are related to the target variables when spreading the samples, we may get more representative samples and then improve the precision of the estimators. Since the samples are spread in all of the auxiliary variables, they will also provide a good basis for model-based inference. Grafström (2012) and Benedetti, Piersimoni, and Postiglione (2015, ch. 7) showed that SCPS was more efficient than GRTS when the auxiliary variables were only the geographical coordinates. Robertson et al. (2013) argued that the statistical performance of BAS was comparable with LPM and SCPS when we only spread the samples in the geographical space. Compared with GRTS and BAS (HIP), the biggest advantage of the LPM and SCPS designs is that they can be used with any type and any number of auxiliary variables. Different from BAS (HIP), the selections are based on distances. By applying LPM and SCPS, we can spread the samples in the auxiliary variables even if they do not constitute dimensions of the population (Grafström & Matei, 2018a).
For environmental monitoring programs that cover large areas, many target variables usually vary rapidly over the landscape with respect to the low sampling intensity (Dobbie, Henderson, & Stevens, 2008). This means that geographically spread samples over the landscape may not optimally capture the distribution of the target variables (Grafström, Zhao, Nylander, & Petersson, 2017). A sample which is well spread in auxiliary variables implies that the sample has the potential to be well spread also for the target variables that are related to the auxiliary variables. Moreover, a sampling design that uses auxiliary variables to spread the sample is particularly useful for multipurpose monitoring programs (Grafström & Schelin, 2014). To avoid misunderstanding, we would like to clarify that a well-spread sample in this article means that the sample is well spread in the auxiliary variables.
Various sampling strategies have been developed based on sample coordination (e.g., Ernst, 1999;Keyfitz, 1951;Kish & Scott, 1971;Patterson, 1950). There are two kinds of coordination: positive coordination and negative coordination. They intend to maximize or minimize the overlap of two or more samples, respectively. We will focus only on positive coordination in order to produce good estimators of change.
Consider sampling from a dynamic population, the coordinated samples selected at different time occasions are dependent. The degree of coordination is measured by the expected size of the overlap between samples. In essence, the positive coordination method based on PRNs consists of assigning a uniformly distributed random number on the interval [0,1] to each unit in the frame. The numbers assigned remain with the units over time, and such a number is called a PRN. These PRNs are used to decide the sampling outcomes (inclusions or exclusions) at each time point.
As the population changes over time, there is a need to update the sample to account for the changes in the auxiliary variables. If we use independent well-spread samples at different time occasions, we can have good estimators of states. However, since the variance of the estimator of change equals the summation of the variances for the two state estimators minus two times their covariance, the estimator of change may not be the best as the covariance between the state estimators will be zero. A permanent sample has a major drawback, as new auxiliary information cannot be entered into the design. A sample that matches the distribution for "today" will be unlikely to match the distribution for "tomorrow." This is because the developments of the sample and the population may be different over time. Thus, if we select a sample that is good at the time of selection, the quality of such a sample is likely to deteriorate over time, which impacts both state and change estimators. Therefore, we need to have an adaptive design that update the sample to maintain the representativeness of the changing population. When a new sample is selected, we want to have overlap (common units) between the new and the previous samples, so that we also improve estimators of change.
Compared with four reference strategies, (i) a strategy that uses permanent geographical-spread samples; (ii) a strategy that employs permanent well-spread samples; (iii) a strategy with independent well-spread samples; (iv) a strategy that applies a split-panel design), the superiority of the new strategy is illustrated by a forest inventory application. By using Monte-Carlo simulations, we show that the new strategy can outperform the reference strategies for both state and change estimators.
This article is mainly intended for those with a solid background in statistics and those who work with survey designs related to environmental monitoring. The rest of the article is organized as follows. In Section 2, a continuous sampling framework and the double sampling approach are introduced, sample coordination for a finite population is explained, and spatial balance is defined. Details about the new strategy are presented in Section 3. In Section 4, we use a simulation study to compare the new strategy against the four reference strategies. Conclusions and comments are given in Section 5.

SAMPLING FRAMEWORK
Let F denote the continuous population that we sample from, and F is assumed to be a bounded open subset of the Euclidean plane R 2 , with surface area (F). The set F is considered to be fixed over time. The response of a target variable for a point x ∈ F at time t is denoted as y t (x). The population total of the response for the target variable at time t can hence be expressed as Y (t) = ∫ F y t (x) dx. The sampling design of size n(t) on F at time t is specified by a joint distribution of n(t) random variables x 1 , … , x n(t) . The sampling intensity at time t is given by is the marginal probability density function of x i at time t. We have ∫ F t (x) dx = n(t).

Double sampling approach
To use auxiliary information in the sampling design, a double sampling approach can be employed for the continuous population. In the first-phase sample, a large number N of locations are selected independently using a constant sampling intensity, (x) = N∕ (F). We let U be the indexes of the geographical locations, U = {1, … , i, … , N}. Then U serve as a dynamic (but permanent with respect to locations and indexes) frame over time. Auxiliary information is then derived for each unit from U at different time occasions. According to the Glivenko-Cantelli theorem and its multivariate generalizations, the empirical distribution of the auxiliary variables in the first-phase sample converges uniformly almost surely to the population distribution as the size of the sample increases (see Dehardt, 1971;Wolfowitz, 1954). Because of the large sample size, the empirical distribution of any variable in the first-phase sample U will closely match the population distribution. The realization of the first-phase sample of size N can be treated as a permanent frame over all repeated surveys. Then, we select the second-phase samples S(t) at different time occasions, S(t) ⊂ U. Denote the first-phase sample of N random locations as S = {x 1 , … , x N }. Conditioned on S, the inclusion probability of unit i ∈ U to be included in Then, the total of the target variable in the first-phase sample at time t can be written as To preserve the relationship between the auxiliary variables and the target variables, we define the response of the auxiliary variables as The target responses are observed for the locations in the second-phase sample in order to estimate the state Y(t). The unbiased Horvitz-Thompson (HT) estimator is then defined aŝ As N is supposed to be much larger than the sample sizes, we allow ourself to do estimation and variance estimation conditioned on the first-phase sample. In that case, the variance of the HT-estimator (1) can be expressed as where ij (t) = Pr (i ∈ S(t), j ∈ S(t)) is the second-order inclusion probability for a pair of points (i, j) at time t. An estimator of (2) isv . ( The estimator (3) is unbiased for (2) provided that all second-order inclusion probabilities are strictly positive.
We define the change of the states between two time points as Δ Y (1, 2) = Y(2) − Y(1) and its estimator asΔ Y (1, 2) = Y (2) −Ŷ (1). SinceŶ (t) is unbiased at any time t, the estimator of change is unbiased as well. When estimating the variance of change from one time occasion to another time occasion, say Occasion 1 and 2, we are also interested in estimating the covariance cov(Ŷ (1),Ŷ (2)). The variance of the change estimator is ) .
The covariance term in (4) is given by where ij (1, 2) = Pr (i ∈ S(1), j ∈ S(2)). It is possible to construct the HT-estimator of (5) based on the two samples, that is, The estimator (6) is unbiased for (5) provided that the ij (1, 2) are strictly positive for all i, j.

F I G U R E 1
The general framework of the strategy for two time occasions. F is the continuous population frame. U stands for the first-phase sample with N = 10, 000. S(1) and S(2) represent the positively coordinated second-phase samples, n(1) = n(2) = 100. The overlapped sample points are marked with solid circles, and the points that are on different locations are marked with triangles and crosses for S(1) and S(2), respectively

Sample coordination for a finite population
The general framework of sample coordination was summarized by Grafström and Matei (2015). Here, we only consider the situation where we select samples from the same population U at different time occasions. Consider only two time occasions, that is, Occasion 1 and 2, the overall sampling design p is defined on U × U, with marginal designs p 1 and p 2 . A random sample of n(1) locations selected at time Point 1 is denoted by S(1), and a random sample of size n(2) selected at time Point 2 is denoted by S(2). The overall sampling design is said to be coordinated if the joint probability of selection of two samples is not equal to the product of the probabilities of selecting each separate sample, that is, if p (s(1), s(2)) ≠ p 1 (s(1)) p 2 (s(2)) (see Cotton & Hesse, 1992;Mach, Reiss, &Şchiopu-Kratina, 2006). Our aim of coordination is to maximize the overlap between several samples drawn successively from U. Therefore, the selection of a new sample will depend on the samples previously drawn. In order to obtain a larger or a smaller overlap of the samples, a dependence between the samples must be introduced. This dependence will determine the expected number of common units in the selected samples. Let  denote the random variable "size of the overlap,"  = ∑ i∈U I i (1)I i (2). The coordination degree between two samples is measured by the expected size of the overlap where i (1, 2) = Pr (i ∈ S(1) ∩ S (2)) is the probability for unit i to be included in both S(1) and S(2). According to Mach et al. (2006), the expected size of the overlap may also have impact on the precision of the change estimators between two occasions. Figure 1 is an illustration of how we construct the finite framework from F and how we select positively coordinated samples from U for two time occasions. In this example, percentage of overlap of the two samples is 74%.

Spatial balance
Spatial balance is a measure to check the spread of a spatial sample. It is often used when the auxiliary space is multidimensional. The measure is based on Voronoi polytopes (Stevens & Olsen, 2004). For a sample of size n(t), we need to construct n(t) polytopes. For each i ∈ s(t), the polytope i (t) includes all units in the population closer to i than to any other sample unit j ∈ s(t), j ≠ i. The distance used when we construct the polytopes is the Euclidean distance on standardized variables. Spatial balance of a sample at time occasion t can be measured by where is the total probability mass within the polytope at time occasion t. Spatial balance can be interpreted as a measure of the variance of the total probability mass within the polytopes. A small value of B(t) indicates the sample is well spread at time t. If the sample is spread perfectly, the total probability mass within i (t) equals to 1. To measure how well a design succeeds in selecting spatially balanced samples, simulation to find the expected value of B(t) under the design is needed.

SAMPLE COORDINATION FOR SPATIALLY BALANCED SAMPLES
By using well-spread samples at each point in time, we are likely to reduce the variances of the state estimators. Increasing the expected size of the overlap between the samples, by a coordinated sample selection, may lead to a higher positive covariance between the state estimators. As Duncan and Kalton (1987) said, the reason for the increased covariance is that many sample units are the same in the two samples and their values tend to be similar at the two time occasions. Thus, by also introducing coordination, we will most likely achieve a smaller variance of the change estimator. Under our framework of coordination for two (or more) samples, the first time occasion sample S(1) and the second time occasion sample S(2) have fixed sample sizes n(1) and n(2), respectively. The sample overlap S(1, 2) = S(1) ∩ S(2) contains n(1, 2) units and n(1, 2) is not fixed. We would like to achieve a high degree of overlap without loosing the spatial balance compared with independent selection of samples.

Spatially correlated Poisson sampling
With the aid of the auxiliary variables derived from the first-phase sample at time occasion t, we select a second-phase sample S(t) of size n(t) from the large first-phase sample. We would like to select a second-phase sample whose distribution of the auxiliary variables matches the distribution of the first-phase sample, and thus also match the population distribution at time t. Assuming the dependence between sampling units decrease as the distance between them increase, to minimize the sampling variance, we should select the sampling units so that we maximize the distance between them (Benedetti et al., 2015, ch. 7). In other words, the sample should be well spread over the auxiliary variables. SCPS is a list-sequential sampling method of selecting well-spread samples. The method was first derived as a spatial application of correlated Poisson sampling, proposed by Bondesson and Thorburn (2008). It is a fixed size ps design that achieves a good spread of the selected samples by using the auxiliary variables. The main idea of SCPS is motivated by a generalization of Tobler's first law of geography. According to that law, geographically nearby locations tend to have more similar properties than locations farther apart. As the distance measure applied in SCPS's algorithm is the standardized Euclidean distance in the space of the auxiliary variables, we can call it "law of auxiliary variables" instead. If the auxiliary variables have high explanatory power for the target variables, then two units with a small distance in the auxiliary space will tend to have more similar values on the target variables than two units farther apart. Generally, in SCPS, the selection of nearby units is avoided to the furthest extent possible, which creates well-spread samples. This is guaranteed by creating a strong negative correlation between the inclusion indicators of nearby units.

Algorithm of SCPS
It is assumed that we have a list U of the units to be sampled. The sampling outcome is first decided for the first unit in the list and then for the second, and so forth. After each sampling decision, the inclusion probabilities for the remaining units in the list are updated. Denote the prescribed inclusion probability of each unit i at time t as i (t), i = 1, 2, … , N, with ∑ N i=1 i (t) = n(t). Then, we have a starting vector of inclusion probabilities ( 1 (t), … , N (t)). In the end of the algorithm, we will get a vector of inclusion indicators by gradually updating the vector of inclusion probabilities in maximum of N steps. The updating can be illustrated as . (9) The first unit is included with probability (0) 1 (t) = 1 (t) at time occasion t. If the first unit was included, we set I 1 (t) = 1, otherwise I 1 (t) = 0. The sampling outcome at each step is decided by comparing the inclusion probability of the step unit with a random number associated with the unit. Denote the random number associated with the unit i ∈ U at time t as r i (t), with r 1 (t), r 2 (t), … , r N (t) i.i.d. U(0, 1). When the values for I 1 (t), … , I j−1 (t) have been decided for the first j − 1 units in the list at time occasion t, the step unit j is included in the sample S(t) at time t, that is, I j (t) = 1, if r j (t) < (j−1) j (t), and I j (t) = 0 otherwise. The inclusion probabilities for the rest of the units in the list are updated according to where i (t) depends on the sampling outcomes of the first j − 1 units. To make sure that 0 ≤ (j) i (t) ≤ 1 holds, the weights need to satisfy the following restrictions −min In the SCPS, the decision (include or not) is always made for the step unit j, that is, the outcome of I j (t) is decided in step j. From Equation (11), we can see that it is possible for the weights to be negative. However, to achieve spatial balance, the weights need to be positive. When updating the inclusion probability for unit i, the weight it receives depends on the distance between i and the step unit j. If the inclusion indicator for the step unit is 1, then based on the inclusion probability in step j − 1, the unit j needs to "steal" more probability from its nearby units. The closer the unit to the step unit, the more probability mass will be "stolen" by the step unit until the updated inclusion probability of the step unit becomes 1, that is, the nearest unit to j will receive as much weight as possible from j, then as much as possible weight will received by the second nearest unit from j, and so forth. By contrast, if the step unit has an inclusion indicator equal to 0, it will give away all its probability mass to its neighbors in a similar way. The above strategy is called maximal weight strategy. Thus, with the restriction of the maximum weight of each unit and the sum of the weights equal to 1, we will have a sample of a fixed size. The algorithm as well as an example of updating the inclusion probabilities of the SCPS can be found in Grafström (2012).

Positive coordination under SCPS
Coordination of well spread samples using the SCPS was introduced for finite populations by Grafström and Matei (2018b). A coordination method based on PRN is used in the algorithm, where the random number associated with each unit in the algorithm of SCPS remains over time, that is, r i (1) = r i (2) = r i for two time occasions. Positive coordination is achieved by using the same comparison rule for all time occasions at each step, that is, I j (t) = 1, if r j < (j−1) j (t), and I j (t) = 0 otherwise. The implementation of positive coordination using SCPS can be found in the R package "BalancedSampling" (Grafström & Lisic, 2019). The algorithm for two time occasions is illustrated in Example 1. i (t) represents the prescribed inclusion probability of unit i at time t, it does not change from Time 1 to Time 2 in this example. r i is the permanent random number associated with unit i. T A B L E 1 Auxiliary variables, prescribed inclusion probabilities, as well as the random numbers associated with the units for two time occasions in Example 1

F I G U R E 2 The distances between units of the auxiliary variable for two time occasions in Example 1
Example 1. Suppose we have a population of size N = 4, we want to select two positively coordinated samples of size n = 2 using SCPS at two time occasions. Table 1 presents the auxiliary value, prescribed inclusion probability as well as the random number associated with each unit at two time occasions. Besides the random number, the prescribed inclusion probability for each unit is also same at each time. According to Table 1, the distances between different units at both time occasions can be calculated, and they are illustrated by Figure 2. At each time occasion, the distance between units is calculated by comparing the value of the auxiliary variable for each unit. The smaller the differences for the values, the closer the units will be. The visiting order is chosen to be 1, 2, 3, 4, the decision is made for each unit according to the same order at both time occasions. The updating of inclusion probabilities for the units at both time occasions for positive coordination of SCPS is illustrated by Figure 3. The maximum weight and the updated inclusion probability for each unit at each step is calculated by using (10) and (11) at both time occasions. According to Figure 3, the overlap for the selected samples is 50% in the example.
When it comes to estimation under positively coordinated samples selected with SCPS, we may use the unbiased HT-estimator (1). However, for designs that produce well-spread samples, many of the second-order inclusion probabilities may be zero. Hence, it is not possible to have a design-based unbiased variance estimator. Based on squared local deviations, Grafström and Schelin (2014) derived an approximate variance estimator for the HT-estimator, under spatially balanced sampling. It can be expressed asv where i ′ is the nearest neighbor to i in S(t) at time t. The distance measure we apply to find i ′ is the same as we used when selecting the sample. Estimation of the covariance between successive state estimators is difficult. The estimator (6) is only unbiased for (5) provided the ij (1, 2) are strictly positive for all i, j. However, that requirement does not hold in general for positively coordinated and spatially balanced samples. The reason is that it is likely that inclusion of a unit i at Time 1 often imply inclusion of unit i at Time 2 and hence also exclusion of neighboring units at Time 2. In such a case, the ij (1, 2) cannot be guaranteed to be strictly positive for all i, j. Finding a suitable estimator for the covariance (5) under positively coordinated and spatially balanced samples remains a challenging problem for the future. F I G U R E 3 Sampling procedure for positively coordinated samples selected using spatially correlated Poisson sampling in Example 1. At each time and step, the order for updating the inclusion probabilities for the remaining units is decided by the distance between the step unit and the remaining units. The distances are measured in auxiliary variables

A FOREST INVENTORY APPLICATION
In forest inventories, we usually sample circular plots or clusters of circular plots. The field inventories are based on the samples selected. Before each survey, the field staff need to find the center of the plots (clusters), then survey each plot with a fixed radius (survey each cluster with a fixed configuration). Since we lack a list frame for the individual objects, the continuous sampling framework can be applied in forest inventories (see e.g. Mandallaz, 2007, ch. 4;Grafström, Schnell, Saarela, Hubbell, & Condit, 2017). Suppose we have a finite population U(t) of objects (e.g., trees) at time t, The number of objects can be different at different time occasions because there will be births and deaths. As the finite population cannot be partitioned into circular plots, it is impractical to sample the objects directly. Thus, we need to construct a continuous framework from the original finite framework, then do the sample selections based on the continuous framework. In this section, the new strategy is compared with four reference strategies by a Swedish National Forest Inventory (NFI) example using simulation. The new sampling strategy as well as the four reference strategies are listed here.
• Strategy 1: The new strategy, which employs well-spread and positively coordinated samples selected by SCPS over time. • Strategy 2: The first reference strategy. Use a permanent geographical-spread sample over time.
• Strategy 3: The second reference strategy. A sampling strategy that uses a permanent sample selected by SCPS, which is well spread at the first time occasion. • Strategy 4: The third reference strategy, which uses independent well-spread samples selected by SCPS over time without sample coordination.
• Strategy 5: The last reference strategy. Use split-panel designs to split the sample into two parts: a panel with a permanent geographical-spread sample and a panel with well-spread samples. This strategy is similar to the current strategy of the Swedish NFI.
A region in the middle of Sweden is selected as our study region. In this region, each cluster consists of 12 circular plots of 7-m radius. Plots in a cluster are placed along a square formation with a side-length of 1,500 m and with 500 m between plots. Denote the response of target at time t as i (t), then the population total at time occasion t can be expressed as Y (t) = ∑ i∈U(t) i (t). When sampling clusters with a given configuration and a fixed orientation over time, the inclusion zone for an object i on location x i is the collection of potential sample points, which lead to inclusions of the object. Mathematically, it can be denoted by K i ⊂ F, K i = K(x i ) = {x ∈ F ∶ x i ∈ (x)}, and (x) is the cluster centered on x. Any cluster (x) with its cluster center x within K i includes the object in one of its plots. Details about inclusion zone and the cluster configuration of the study region can be found in Grafström, Zhao, et al. (2017). Different formulations for the density function of the target variable have been discussed by Grafström, Schnell, et al. (2017). For constructing the continuous framework of the NFI application, the density function at time t can be expressed as where I ti (x) = 1 if x ∈ K i at time t, and 0 otherwise, (K i ) is the area of the inclusion zone for object i. By using the expression (13), the continuous population total is identical to the corresponding finite population total. This is because Once we have constructed the continuous framework, we can easily employ the general framework in Section 2 for our example.
First, a number of 100,000 clusters were independently selected as an initial first-phase sample, then a subset of size 10,000 clusters were selected using SCPS as the first-phase sample, since the algorithm is quite computationally intensive. (This is only needed to do the simulation, we use the initial first-phase sample as the first-phase sample in reality because we only select the second-phase sample once.) Then, from the first-phase sample, a sample of size 100 was selected as the second-phase sample for both time occasions, respectively, using different strategies.
The auxiliary variables we used were the geographical coordinates, the mean tree height, the mean basal area, and the mean elevation. Denote the five auxiliary variables at the plot level as q k (t) = ( q x k (t), q y k (t), q h k (t), q b k (t), q e k (t) ) T ∈ R 5 .
According to Grafström, Zhao, et al. (2017), we only had auxiliary information for one time occasion. Stand-level growth models were applied to generate the 5-year's growth in the plot level for the mean tree height and the mean basal area. The growth models were based on data from permanent samples of the Swedish NFI established during 1983-1987 and reinventoried three to four times between 1988 and 2010 (Fridman, Holm, Nilsson, Ringvall, & Ståhl, 2014). We also applied a clear cutting with a rate of 5% for the second time occasion. First, 20% of the plots who had the highest mean tree height values were selected as the potential plots. Then, 1/4 of the plots among them were randomly selected for clear cutting. Based on the growth and the clear cutting, we got the mean tree height and the mean basal area at the next time occasion with a 5-years' time difference. The geographical coordinates as well as the elevation remained the same at time Occasion 2. The 5-years' growth models are expressed for the mean tree height and the mean basal area in the plot level in Equations (15) and (16), respectively.
where h k ∼ N(0, h k ) and b k ∼ N(0, b k ). Before employing the auxiliary information to the sampling design, we aggregated the plot level auxiliary information q k (t) to cluster-level z i (t). Whether a plot contributes to the aggregated cluster value or not depends on if the plot center is inside of the region or not.

TA B L E 2 Example National Forest Inventory
) V ) stands for the variance of the estimator of change. P ij is the split panel design, the value of i corresponding to the percentage of the permanent sample and value of j corresponding to the percentage of temporary sample. For example, P 28 means the permanent panel is 20% and the temporary panel is 80%.
We applied the growth models since we would like to generate more realistic auxiliary variables at the second time occasion. Separate research can be done within this topic. Since we do not have any target variable in the simulation, results are only presented for auxiliary variables. As a multipurpose inventory, there are a wide range of target variables in the Swedish NFI, some examples can be the proportions of different land types, volume and number of trees per hectare, mean age, damages, amount of different berries, and so forth. Among them, most of the variables are related to the selected auxiliary variables.
The simulation results of the example is shown in Table 2 for the new strategy and the four reference strategies. According to Dobbie et al. (2008), the optimal panel design depends on the balance between needing to detect trend and report on the states, four different partition schemes are applied for Strategy 5 to test if there is an optimal way to split the two panels. The mean of spatial balance for samples at two time occasions, the percentage of overlap, the empirical variances of the estimator of the two states as well as change are listed for mean tree height and mean basal area.
We can clearly see that both state and change estimators obtained by employing the new strategy (Strategy 1) are better than what are possible by using reference strategies. There is no any optimal combination of the two panels for the current Swedish NFI strategy. The higher the proportions of the permanent panels, the worse the estimators will be. The main results can be summarized from two aspects: effect of using spatially balanced sampling designs and effect of the amount of overlap.
Strategies 1 and 4 that apply spatially balanced sampling designs produce the best estimators for states. Theoretically, the values of the estimators tend to be the same for the two strategies when the number of repetitions in the simulation is large enough. Strategy 2, which only spreads the sample in the geographical space, leads to the worst estimators for states. This confirms the importance of spreading the samples also in auxiliary variables other than the geographical space. Comparing Strategies 1 and 4 against Strategy 2, reduction of variances of the state estimators for both auxiliary variables is more than 92%, and the decrease in width of confidence intervals for both auxiliary variables is more than 70%. For Strategy 3, the quality of the state estimators reduced at the second time occasion, since the permanent sample is not as well spread anymore at that time point. For the strategy that applies split panel designs, the value of the variances for state estimators reduced as we increased the proportions of the well-spread panel.
Comparing Strategy 1 against Strategy 4, the reduction of variances for the estimators of change is more than 37%, and the decrease in width of confidence intervals is more than 20% for both auxiliary variables. This is because we get a much higher overlap (62%) between the samples selected at both time occasions using Strategy 1, which leads to higher covariances between the two state estimators compared with Strategy 4. For mean tree height, the estimate for covariance is 0.288 and −0.012 for Strategy 1 and Strategy 4, respectively. For mean basal area, the values are 0.006 and 8.21 × 10 −5 , respectively. Based on Equation (4), when fixing the state variances, the strategy that has the larger covariance will produce a smaller variance for the estimator of change. At first glance, results produced by Strategy 5 do not seem to be reasonable in terms of the estimators of change. It appears to be a contradictory statement that the variances of estimators of change increase as the percentages of overlap increase. However, by careful observations, we can find that the reason for the increased variances of change estimators is not because we do not have higher values of the covariances. For tree height, the covariance between two state estimators increased from 1.115 to 3.291 when we increased the overlap from 20% to 80%. For basal area, the covariance increased from 0.001 to 0.075. The reason of the increased variances of the change estimators are the increased values of variances for state estimators. According to Equation (4), when the increase in covariance cannot compensate the increase in the two variances of state estimators, the variance of estimator of change will still increase, even if we have a very high percentage of overlap.

CONCLUSION AND DISCUSSION
We proposed a new sampling strategy with its main focus on monitoring totals of environmental variables. In practice, it is not restricted to monitor only the totals, it can also be applied for parameters such as quantiles (Grafström & Schelin, 2014). Positive coordination is studied by using SCPS within a continuous population framework. Based on an application, with settings similar to the Swedish NFI, we illustrated that the proposed new strategy performed better than all reference strategies. When matching the sample distribution of the auxiliary variables to the population distribution and at the same time use positively coordinated samples, we improve the precision for both the state and the change estimators.
If we use a sample, which is well spread only at the first time occasion as a permanent sample, then there is a big risk that the sample evolves differently from the population. The sample may become less balanced over time and, as a result, the state estimators also become less efficient over time. Although the sample overlap is 100% for a permanent sample, the estimator of change will also gradually become worse over time as the sample can differ more and more from the population in terms of the distributions of the auxiliary variables (and hence the target variables).
For the current Swedish NFI, two different types of clusters are used: permanent clusters and temporary clusters. In fact, the Swedish NFI design that combines temporary and permanent clusters has almost become an international standard toward which NFIs in other countries' aim (Fridman et al., 2014). The permanent clusters primarily aim to increase the accuracy of change estimation, and they are resurveyed regularly, whereas temporary ones are mainly intended to capture the current state of the forest and are only surveyed once (Tomppo, Gschwantner, Lawrence, & McRoberts, 2010, ch. 35). In the NFI example, we use positively coordinated and spatially balanced samples to target change and the current states of two time occasions simultaneously. Based on the simulation, we can see that the new strategy successfully improves the precision of the estimator for both state and change for the auxiliary variables. Thus, it has potential to also improve the precision of the estimators for the target variables that are related to the auxiliary variables. For those target variables that are not related to the selected auxiliary variables, by spreading the samples in the auxiliary variables, we get similar results as independent sample selections. Therefore, there is a potential to change the current design of the Swedish NFI. Instead of using the complex and less efficient combination of the permanent and the temporary samples, employing positively coordinated and well-spread samples can achieve both goals within a single sampling strategy.
A planner has plenty of options in choosing a sampling strategy. The main properties considered often include precision, unbiasedness, cost-efficiency, and simplicity to apply. As Scott (1984) mentioned, when estimating both state and change, a combination of remeasured (matched) plots, plots not remeasured (unmatched), and replacement (new) plots is generally the most cost-effective alternative. With the suggested strategy, at each point in time, the sample is a sample of the SCPS design with the prescribed inclusion probabilities. Thus, we make no compromise on the level of spatial balance of the different samples. Yet, we achieve a quite high degree of positive coordination. As the sample size can be varied over time, the strategy is also flexible for budget changes over time. With auxiliary variables available, they can be seen as proxies for the target variables. The sampling strategy that is superior for estimating the state and change for the auxiliary variables is likely to be superior for estimating the state and change of target variables related to those auxiliary variables.
Many spatially balance sampling designs can improve the state estimators compared with traditional designs (Benedetti et al., 2015, ch. 7). We focus only on SCPS among others, since it is efficient and easy to apply when it comes to positive sample coordination. The main reason of not including the GRTS or BAS designs as reference strategies in the application is because we have additional auxiliary variables beside the geographical coordinates, which are not dimensions of the population. If we add an auxiliary variable such as elevation to the geographical coordinates, then we get a surface with zero three-dimensional volume. By enclosing the surface in a three-dimensional rectangular box and generate random points in the box, there is a zero probability to get points that lie on the surface. Thus, for example, BAS fails to use the additional information and can only spread the sample in the geographical coordinates.
The sample coordination method presented here is a probabilistic way for SCPS to define panels. This means by using this method we cannot fix the panels beforehand as we could do by using the traditional panel designs. The overlap between the successive samples at two time occasions depends on the change in the auxiliary variables between the two time occasions. Further studies on how to select well-spread samples with a prescribed percentage of overlap is of a great interest to us.