11-Twophase.Rmd

# Two-phase random sampling {#Twophase}

The regression and ratio estimators of Chapter \@ref(Modelassisted) require that the means of the ancillary variables are known. If these are unknown, but the ancillary variables can be measured cheaply, one may decide to estimate the population means of the ancillary variables from a large sample. The study variable is measured in a random subsample of this large sample only. This technique is known in the sampling literature as two-phase random sampling\index{Two-phase random sampling} or double sampling\index{Double sampling}. Another application of two-phase sampling is two-phase sampling for stratification. Stratified random sampling (Chapter \@ref(STSI)) requires a map with the strata. The poststratified estimator of Subsection \@ref(PoststratifiedEstimator) requires that the sizes of the strata are known. With two-phase sampling for stratification, neither a map of the strata nor knowledge of the stratum sizes is required. Note that the term `phase' does not refer to a period of time; all data can be collected in one sampling campaign. Let me also explain the difference with two-stage cluster sampling (Chapter \@ref(Twostage)). In two-stage cluster random sampling, we have two types of sampling units, clusters of population units and individual population units. In two-phase sampling, we have one type of sampling unit only, the objects of a discrete population or the elementary sampling units of a continuous population (Section \@ref(BasicConcepts)). 

In two-phase sampling for regression and two-phase sampling for stratification, the two phases have the same aim, i.e., to estimate the population mean of the study variable. The observations of the covariate(s) and/or strata in the first phase are merely done to increase the precision of the estimated mean of the study variable. Another application of two-phase sampling is subsampling an existing probability sample designed for a different aim. So, in this case the study variable observed in the second-phase sample may not be related to the variables observed in the first-phase sample. 

An example is LUCAS-Topsoil [@Ballabio2019] which is a subsample of approximately 22,000 units sampled from a much larger sample, the LUCAS sample, designed for estimating totals of land use and land cover classes across the European Union. It was not feasible to observe the soil properties at all sites of the LUCAS sample, and for this reason a subsample was selected. Regrettably, this subsample is not a probability sample from the LUCAS sample: the inclusion probabilities are either zero or unknown. Design-based or model-assisted estimation of means of soil properties for domains of interest is not feasible. The only option is model-based prediction.

In case the subsample is a probability subsample from the first-phase sample and no variable observed in the first-phase sample\index{First-phase sample} is of use for estimating the total or mean of the study variable observed in the subsample, the population total can be estimated by the $\pi$ estimator:

\begin{equation}
\hat{t}(z) =\sum_{k \in \mathcal{S}_2}\frac{z_k}{\pi_{1k}\pi_{k|\mathcal{S}_1}} =\sum_{k \in \mathcal{S}_2}\frac{z_k}{\pi^*_{k}}
\;,
(\#eq:HTestimatorTotalDoubleSampling)
\end{equation}

with $\pi_{1k}$ the probability that unit $k$ is selected in the first phase, and $\pi_{k|\mathcal{S}_1}$ the probability that unit $k$ is selected in the second phase\index{Second-phase sample}, given the first-phase sample $\mathcal{S}_1$. This general $\pi$ estimator for two-phase sampling, referred to as the $\pi^*$ estimator\index{$\pi^*$ estimator} by @sar92, can be used for any combination of probability sampling designs in the first and second phase.

To derive the variance, it is convenient to write the total estimation error as the sum of two errors:

\begin{equation}
\begin{split}
\hat{t}(z)-t(z) &=\left(\sum_{k \in \mathcal{S}_1}\frac{z_k}{\pi_{1k}}-t(z)\right)+\left(\sum_{k \in \mathcal{S}_2}\frac{z_k}{\pi^*_{k}}-\sum_{k \in \mathcal{S}_1}\frac{z_k}{\pi_{1k}}\right)\\
& =e_1+e_2 \;.
\end{split}
(\#eq:DecomposeErrorsDoubleSampling)
\end{equation}

The first error $e_1$ is the error in the estimated population total, as estimated by the usual $\pi$ estimator using the study variable values for the units in the first-phase sample. This estimator cannot be computed in practice, as the study variable values are only known for a subset of the units in the first-phase sample. The second error $e_2$ is the difference between the $\pi^*$ estimator using the study variable values for the units in the subsample only, and the $\pi$ estimator using the study variable values for all units in the first-phase sample.

The variance of the $\pi^*$ estimator can be decomposed into the variance of these two errors as follows:

\begin{equation}
V_{p_1,p_2}(\hat{t})=V_{p_1}E_{p_2}(\hat{t}|\mathcal{S}_1)+E_{p_1}V_{p_2}(\hat{t}|\mathcal{S}_1)=V_{p_1}(e_1)+E_{p_1}V_{p_2}(e_2|\mathcal{S}_1)\;,
(\#eq:VarHTestimatorDoubleSampling)
\end{equation}

with $V_{p_1}$ and $E_{p_1}$ the variance and expectation of the estimator for the population total over repeated sampling with the design of the first phase, respectively, and $V_{p_2}$ and $E_{p_2}$ the variance and expectation of the estimator for the population total over repeated sampling with the design of the second phase, respectively. The population mean can be estimated by the estimated total divided by the population size $N$. 

## Two-phase random sampling for stratification {#TwophaseStratification}

In two-phase sampling for stratification\index{Two-phase sampling!for stratification}, in the first phase a large sample is taken and the selected sampling units are all classified. The classes thus formed are then used as strata in the second sampling phase. A stratified subsample is selected, and the study variable is observed on the units in the subsample only.

```{block2, type='rmdnote'}
This sampling design is applied, for instance, to monitor land use and land cover in the European Union by the LUCAS monitoring network mentioned above. In the first phase, a systematic random sample is selected, consisting of the nodes of a square sampling grid with a spacing of 2 km. Land use and land cover (LULC) are then determined at the selected grid nodes, using orthophotographs, satellite imagery, and fieldwork. The idea is that this procedure results in a more accurate classification of LULC at the selected units than by overlaying the grid nodes with an existing LULC map such as the Corine Land Cover map. The site-specific determinations of LULC classes are then used to select a stratified random subsample (second-phase sample). In 2018 the monitoring network was redesigned [@Eurostat2018].
```

Two-phase sampling for stratification is now illustrated with study area Voorst. A map with five combinations of soil type and land use is available of this study area. These combinations were used as strata in Chapter \@ref(STSI), and the stratum sizes were used in the poststratified estimator of Subsection \@ref(PoststratifiedEstimator). Here, we consider the situation that we do not have this map and that we do not know the sizes of these strata either. In the first phase, a simple random sample of size 100 is selected. In the field, the soil-land use combination is determined for the selected points, see Figure \@ref(fig:DoubleSampleVoorst). This time we assume that the field determinations are equal to the classes as shown on the map.

```{r}
n1 <- 100
set.seed(123)
N <- nrow(grdVoorst)
units <- sample(N, size = n1, replace = FALSE)
mysample <- grdVoorst[units, ]
```

The simple random sample is subsampled by stratified simple random sampling, using the soil-land use classes as strata. The total sample size of the second phase is set to 40. The number of points in the simple random sample per stratum is determined. Then the subsample size per stratum is computed for proportional allocation. Finally, function `strata` of package **sampling** [@Tille2016] is used to select a stratified simple random sample without replacement, see Chapter \@ref(STSI) for details. At the 40 points of the second-phase, the soil organic matter (SOM) concentration is measured.

```{r}
library(sampling)
n2 <- 40
n1_h <- tapply(mysample$z, INDEX = mysample$stratum, FUN = length)
n2_h <- round(n1_h / n1 * n2, 0)
units <- sampling::strata(mysample, stratanames = "stratum",
  size = n2_h[unique(mysample$stratum)], method = "srswor")
mysubsample <- getdata(mysample, units)
table(mysubsample$stratum)
```

```{r DoubleSampleVoorst, echo = FALSE, out.width = "100%", fig.cap = "Two-phase random sample for stratification from Voorst. Coloured dots: first-phase sample of 100 points selected by simple random sampling, with observations of the soil-land use combination.  Triangles: second-phase sample of 40 points selected by stratified simple random subsampling of the first-phase sample, using the soil-land use combinations as strata, with measurements of SOM."}
ggplot() +
  geom_raster(data = grdVoorst, mapping = aes(x = s1 / 1000, y = s2 / 1000), fill = "grey") +
  geom_point(data = mysample, mapping = aes(x = s1 / 1000, y = s2 / 1000, colour = stratum), size = 1.5) +
  scale_colour_manual(name = "Stratum", values = c(BA = "darkgreen", EA = "brown", PA = "orange", RA = "green", XF = "grey")) +
  geom_point(data = mysubsample, mapping = aes(x = s1 / 1000, y = s2 / 1000), size = 2.5, shape = 2) +
  scale_x_continuous(name = "Easting (km)") +
  scale_y_continuous(name = "Northing (km)") +
  coord_fixed() +
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.background = element_blank())
```

With simple random sampling in the first phase and stratified simple random sampling in the second phase, the population mean can be estimated by

\begin{equation}
\hat{\bar{z}}=
\sum_{h=1}^{H_{\mathcal{S}_1}}\frac{n_{1h}}{n_{1}}\,\bar{z}_{\mathcal{S}_{2h}} 
\;,
(\#eq:EstimatedMeanDoubleStratification)
\end{equation}

where $H_{\mathcal{S}_1}$ is the number of strata used for stratification of the first-phase sample, $n_{1h}$ is the number of units in the first-phase sample that form stratum $h$ in the second phase, $n_{1}$ is the total number of units of the first-phase sample, and $\bar{z}_{\mathcal{S}_{2h}}$ is the mean of the subsample from stratum $h$.

```{r}
mz_h_subsam <- tapply(mysubsample$z, INDEX = mysubsample$stratum, FUN = mean)
mz <- sum(n1_h / n1 * mz_h_subsam, na.rm = TRUE)
```

The estimated population mean equals `r formatC(mz, 1, format = "f")` g kg^-1^. The sampling variance over repeated sampling with both designs can be approximated^[In the approximation, it is assumed that $N$ is much larger than $n_1$, and $(n_{1h}-1)/(n_1-1)$ is replaced by $n_{1h}/n_1$.] by (@sar92, equation at bottom of p. 353)

\begin{equation}
\widehat{V}\!(\hat{\bar{z}}) = \sum_{h=1}^{H_{\mathcal{S}_1}}\left( \frac{n_{1h}}{n_1}\right)^2
\frac{\widehat{S^2}_{\mathcal{S}_{2h}}}{n_{2h}} + \frac{1}{n_1}\sum_{h=1}^{H_{\mathcal{S}_1}} \frac{n_{1h}}{n_1}\left( \bar{z}_{\mathcal{S}_{2h}}-\hat{\bar{z}}\right)^2 
\;,
(\#eq:VarMeanDoubleStratification)
\end{equation}

with $\widehat{S^2}_{\mathcal{S}_{2h}}$ the variance of $z$ in the subsample from stratum $h$.


```{r}
S2z_h_subsam <- tapply(mysubsample$z, INDEX = mysubsample$stratum, FUN = var)
w1_h <- n1_h / n1
v_mz_1 <- sum(w1_h^2 * S2z_h_subsam / n2_h)
v_mz_2 <- 1 / n1 * sum(w1_h * (mz_h_subsam - mz)^2)
se_mz <- sqrt(v_mz_1 + v_mz_2)
```

The estimated standard error equals `r formatC(se_mz, 1, format = "f")` g kg^-1^.

The mean and its standard error can be estimated with functions `twophase` and `svymean` of package **survey** [@Lumley2020]. A data frame with the first-phase sample is passed to function `twophase` using argument `data`. A variable in this data frame, passed to function `twophase` with argument `subset`, is an indicator with value `TRUE` if this unit is selected in the second phase, and `FALSE` otherwise. 

```{r}
library(survey)
lut <- data.frame(stratum = sort(unique(mysample$stratum)), fpc2 = n1_h)
mysample <- mysample %>%
  mutate(ind = FALSE,
         fpc1 = N) %>%
  left_join(lut, by = "stratum")
mysample$ind[units$ID_unit] <- TRUE
design_2phase <- survey::twophase(
  id = list(~ 1, ~ 1), strata = list(NULL, ~ stratum),
  data = mysample, subset = ~ ind, fpc = list(~ fpc1, ~ fpc2))
svymean(~ z, design_2phase)
```


As shown in the next code chunk, the standard error is computed with the original variance estimator, without approximation (@sar92, equation (9.4.14)).

```{r}
v_mz_1 <- 1 / N^2 * N * (N - 1) *
  sum((((n1_h - 1) / (n1 - 1)) - ((n2_h - 1) / (N - 1)))
      * w1_h * S2z_h_subsam / n2_h)
v_mz_2 <- 1 / N^2 * (N * (N - n1)) / (n1 - 1) *
  sum(w1_h * (mz_h_subsam - mz)^2)
sqrt(v_mz_1 + v_mz_2)
```

## Two-phase random sampling for regression  {#TwophaseRegression}

The simple regression estimator of Equation \@ref(eq:SimpleRegressionEstimatorSI) requires that the population mean of the ancillary variable $x$ is known. This section, however, is about applying the regression estimator in situations where the mean of $x$ is unknown\index{Two-phase sampling!for regression}. A possible application is estimating the soil organic carbon (SOC) stock in an area. To estimate this carbon stock, soil samples are collected and analysed in a laboratory. The laboratory measurements can be very accurate, but also expensive. Proximal sensors can be used to derive soil carbon concentrations from the spectra. Compared to laboratory measurements of soil the proximal sensor determinations are much cheaper, but also less accurate. If there is a relation between the laboratory and the proximal sensing determinations of SOC, then we expect that the regression estimator of the carbon stock will be more accurate than the $\pi$ estimator which does not exploit the proximal sensing measurements. However, the population mean of the proximal sensing determinations is unknown. What we can do is estimate this mean from a large sample. Additionally, for a subsample of this large sample, SOC concentration is also measured in the laboratory. This is another example of two-phase sampling.

Intuitively, we understand that with two-phase sampling the variance of the regression estimator of the total carbon stock is larger than when the population mean of the proximal sensing determinations is known. There is a sampling error in the estimated population mean of the proximal sensing determinations, estimated from the large first-phase sample, and this error propagates to the error in the estimated total carbon stock. 

Two-phase sampling for regression is now illustrated with Eastern Amazonia (Subsection \@ref(Amazonia)). The study variable is the aboveground biomass (AGB), and lnSWIR2 is used here as a covariate. We do have a full coverage map of lnSWIR2, so two-phase sampling with a large first-phase sample to estimate the population mean of lnSWIR2 is not needed. Nevertheless, hereafter a two-phase sample is selected, and the population mean of lnSWIR2 is estimated from the first-phase sample. In doing so, the effect of ignorance of the population mean of the covariate on the variance of the regression estimator becomes apparent.

In the next code chunk, a first-phase sample of 250 units, the dots in the plot, is selected by simple random sampling without replacement. In the second phase a subsample of 100 units, the triangles in the plot, is selected from the 250 units by simple random sampling without replacement. At all 250 units of the first-phase sample the covariate lnSWIR2 is measured, whereas AGB is measured at the 100 subsample units only. 

```{r twophasesample}
grdAmazonia <- grdAmazonia %>%
  mutate(lnSWIR2 = log(SWIR2))
n1 <- 250; n2 <- 100
set.seed(314)
units_1 <- sample(nrow(grdAmazonia), size = n1, replace = FALSE)
mysample <- grdAmazonia[units_1, ]
units_2 <- sample(n1, size = n2, replace = FALSE)
mysubsample <- mysample[units_2, ]
```

Figure \@ref(fig:twophaseAmazonia) shows the selected two-phase sample. 

```{r twophaseAmazonia, echo = FALSE, out.width = "100%", fig.cap = "Two-phase random sample for the regression estimator of the mean AGB in Eastern Amazonia. Coloured dots: simple random sample without replacement of 250 units with measurements of covariate lnSWIR2 (first-phase sample). Triangles: simple random subsample without replacement of 100 units with measurements of AGB (second-phase sample)."}
ggplot(data = grdAmazonia) +
  geom_raster(mapping = aes(x = x1 / 1000, y = x2 / 1000), fill = "grey") +
  geom_point(data = mysample, mapping = aes(x = x1 / 1000, y = x2 / 1000, colour = lnSWIR2), size = 1.5) +
  scale_colour_viridis_c(name = "lnSWIR2") +
  geom_point(data = mysubsample, mapping = aes(x = x1 / 1000, y = x2 / 1000), shape = 2, size = 3) +
  scale_x_continuous(name = "Easting (km)") +
  scale_y_continuous(name = "Northing (km)") +
  coord_fixed()
```

Estimation of the population mean or total by the regression estimator from a two-phase sample is very similar to estimation when the covariate mean is known, as described in Subsection \@ref(RegressionEstimator) (Equation \@ref(eq:SimpleRegressionEstimatorSI))\index{Regression estimator!for two-phase sampling}. The observations of the *subsample* can be used to estimate the regression coefficient $b$. The true population mean of the ancillary variable, $\bar{x}$ in Equation \@ref(eq:SimpleRegressionEstimatorSI), is unknown now. This true mean is replaced by the mean as estimated from the relatively large first-phase sample, $\bar{x}_{\mathcal{S}_1}$. The estimated mean of the covariate, $\bar{x}_{\mathcal{S}}$ in Equation \@ref(eq:SimpleRegressionEstimatorSI), is estimated from the subsample, $\bar{x}_{\mathcal{S}_2}$. This leads to the following estimator:

\begin{equation}
\hat{\bar{z}}= \bar{z}_{\mathcal{S}_2}+\hat{b}\left( \bar{x}_{\mathcal{S}_1}-\bar{x}_{\mathcal{S}_2}\right) \;,
(\#eq:RegressionEstimatorTwoPhase)
\end{equation}

where $\bar{z}_{\mathcal{S}_2}$ is the subsample mean of the study variable, and $\bar{x}_{\mathcal{S}_1}$ and $\bar{x}_{\mathcal{S}_2}$ are the means of the covariate in the first-phase sample and the subsample (i.e., the second-phase sample), respectively.

The sampling variance is larger than that of the regression estimator with known mean of $x$. The variance can be decomposed into two components. The first component is equal to the sampling variance of the $\pi$ estimator of the mean of $z$ with the sampling design of the first phase (in this case, simple random sampling without replacement), supposing that the study variable is observed on all units of the first-phase sample. The second component is equal to the sampling variance of the regression estimator of the mean of $z$ in the first-phase sample, with the design of the second-phase sample (again simple random sampling without replacement in this case):

\begin{equation}
\widehat{V}\!\left(\hat{\bar{z}}\right)=(1-\frac{n_1}{N})\frac{\widehat{S^{2}}(z)}{n_1} + (1-\frac{n_2}{n_1}) \frac{\widehat{S^{2}}(e)}{n_2} \;,
(\#eq:VarianceRegressionEstimatorTwoPhase)
\end{equation}

with $\widehat{S^{2}}(e)$ the variance of the regression residuals as estimated from the subsample:

\begin{equation}
\widehat{S^{2}}(e)=\frac{1}{(n_2-1)}\sum_{k \in \mathcal{S}_2}e_{k}^2 \;.
(\#eq:VarianceResidualsTwoPhase)
\end{equation}

The ratios $(1-n_1/N)$ and $(1-n_2/n_1)$ in Equation \@ref(eq:VarianceRegressionEstimatorTwoPhase) are finite population corrections (fpcs). These fpcs account for the reduced variance due to sampling the finite population and subsampling the first-phase sample without replacement.

```{r}
lm_subsample <- lm(AGB ~ lnSWIR2, data = mysubsample)
ab <- coef(lm_subsample)
mx_sam <- mean(mysample$lnSWIR2)
mx_subsam <- mean(mysubsample$lnSWIR2)
mz_subsam <- mean(mysubsample$AGB)
mz_reg2ph <- mz_subsam + ab[2] * (mx_sam - mx_subsam)
```

The estimated population mean equals `r formatC(mz_reg2ph, 1, format = "f")` 10^9^ kg ha^-1^. The standard error can be approximated as follows.

```{r}
e <- residuals(lm_subsample)
S2e <- sum(e^2) / (n2 - 1)
S2z <- var(mysubsample$AGB)
N <- nrow(grdAmazonia)
se_mz_reg2ph <- sqrt((1 - n1 / N) * S2z / n1 + (1 - n2 / n1) * S2e / n2)
```

The estimated standard error equals `r formatC(se_mz_reg2ph, 2, format = "f")` 10^9^ kg ha^-1^.

The regression estimator for two-phase sampling and its standard error can also be computed with package **survey**, as shown below. The standard error differs from the standard error computed above because it is computed with the g-weights, see Subsection \@ref(RegressionEstimator). Note argument `fpc = list(~ N, NULL)`. There is no need to add the first-phase sample size as a second element of the list, because this sample size is simply the number of rows of the data frame. Setting the second element of the list to `NULL` does not mean that the standard error is computed for sampling with replacement in the second phase. Function `twophase` assumes that the second-phase units are always selected without replacement.

```{r}
mysample <- mysample %>%
  mutate(id = row_number(),
         N = N,
         ind = id %in% units_2)
design_2phase <- survey::twophase(
  id = list(~ 1, ~ 1), data = mysample, subset = ~ ind, fpc = list(~ N, NULL))
mysample_cal <- calibrate(
  design_2phase, formula = ~ lnSWIR2, calfun = "linear", phase = 2)
svymean(~ AGB, mysample_cal)
```

#### Exercises {-}

1. Write an **R** script to select a simple random sample without replacement of 250 units from Eastern Amazonia and a subsample of 100 units by simple random sampling without replacement. Repeat this 1,000 times in a for-loop. 
    +  Use each sample selected in the first phase (sample of 250 units) to estimate the population mean of AGB by the regression estimator (Equation \@ref(eq:SimpleRegressionEstimatorSI) in Chapter \@ref(Modelassisted)). Assume that the AGB data are known for all units selected in the first phase. Use lnSWIR2 as a covariate.
    +  Compute the variance of the 10,000 regression estimates of the population mean of AGB.
    +  Use each two-phase sample to compute the regression estimator of AGB for two-phase sampling (Equation \@ref(eq:RegressionEstimatorTwoPhase)). Now only use the AGB data of the subsample. Estimate the population mean of lnSWIR2 from the first-phase sample of 250 units. Approximate the variance of the regression estimator for two-phase sampling (Equation \@ref(eq:VarianceRegressionEstimatorTwoPhase)).
    +  Compute the variance of the 10,000 regression estimates of the population mean of AGB for the two-phase sampling design.
    +  Compare the two variances and explain the difference.
    +  Compute the average of the 10,000 approximate variances and compare the result with the variance of the 10,000 estimated means, as estimated by the regression estimator for two-phase sampling. 
  
```{r, echo = FALSE}
rm(list = ls())
```