Model-assisted sample design is minimax for model-based prediction
Section 3. Simulation study

• Previous
• Next

A simulation study was conducted to compare the AV of the BLUP and its upper bound in situations where $Y$ has a nonlinear relationship with a continuous auxiliary variable $x_1.$ A second auxiliary variable, $x_2,$ is binary and independent of $x_1$ and $Y.$ Probabilities of selection depend on $x_1$ and $x_2$ in various ways. For each scenario, 5,000 populations and samples were generated, with a population size of 6,000 and sample sizes of 500 and 1,500, and with a population size of 100,000 and a sample size of 25,000. All code is available at www.github.com/rgcstats/AVLB.

3.1  Simulation of populations

The population values of $x_1$ were the $\left\{j/\left(101\right):j\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}100\right\}$ quantiles of a lognormal distribution with mean $-1/32$ and standard deviation 0.25, with equal frequencies of each of these 100 values. This means that $x_1$ are positive and right-skewed with a mean of 1 and a range of 0.54 to 1.74. The values of $x_1$ were non-stochastic and discretized in order to speed up computation, to simplify the generation of smooth models (see below), and to facilitate comparison of AVs to the GJLB by making the GJLB constant across simulations. A second binary variable $x_2$ took on the values 0 and 1 with equal frequency within each value of $x_1,$ so that the two covariates were orthogonal.

Conditional on $x_1$ and $x_2,$ the population values of $Y$ were generated independently as

$$E_M{Y}_i\mathrm{<mo>=</mo>}\mu \left(x_{1i}\right)+{ϵ}_i(3.1)$$

where

$${ϵ}_i\sim N\left(\mathrm{0,}\text{0}\text{.25}{x}_{1i}\right)\mathrm{.}(3.2)$$

The mean function $\mu \left(\mathrm{.}\right)$ was a smooth but nonlinear function,

$$\mu \left(x_1\right)\mathrm{<mo>=</mo>}4x_1+\sin \left(x_{1}2\pi /h\right)\mathrm{,}(3.3)$$

consisting of a linear term and a sinusoidal term with period $h.$ When $h$ is large, $\mu \left(\right)$ is close to linear over the range of $x_1$ in the population, while for $h$ small there are frequent cycles in the function. Figure 3.1 shows the mean function for the periods used in the simulation (0.5, 1, 2, 5).

3.2  Simulated sampling

The probabilities of selection, ${\pi }_i,$ were set to

$${\pi }_i\propto x_{1i}^b\left(1+c{x}_{2i}\right)(3.4)$$

where $b$ was 0.5, 1 or 2, reflecting light, medium or high dependence on $x_{1i}.$ The values of $c$ were 0, 0.5 or 1.5, reflecting no, medium or high dependence on $x_{2i}.$ (Other values of $c$ were also used for the purposes of Figure 3.2 only.)

The second auxiliary variable, $x_2,$ is unrelated to $Y$ but may affect the selection probabilities. One might expect the BLUP to do better compared to the GJLB when probabilities of selection depend on $x_2,$ since the BLUP may be based on a model omitting $x_2,$ potentially leading to lower variance. Of course, there is also the possibility of the working model omitting $x_1,$ leading to the BLUP being biased, however this never occurred in any of the simulations. To explore robustness to incorrect omission of $x_1,$ it would be of interest to consider relationships weaker than those shown in Figure 3.1, but this was beyond the scope of the paper.

Inclusion probabilities are forced to obey the proportionality in (3.4) but are truncated above at 1 and below at 1/40 and scaled such that they add to the required sample size after truncation. Samples are selected by unequal probability systematic sampling with random ordering using the sampling package in R (Tillé and Matei, 2016).

3.3  Estimation of the population total of $Y$

A linear model in $x_1$ and spline models in $x_1$ with between 1 and 10 interior knots are fitted to each sample. (See for example Breidt, Claeskens and Opsomer, 2005 for the use of splines in model-assisted survey estimation.) Another 11 models are defined by also including $x_2$ as an additive covariate. The model with the lowest BIC is then used to calculate a BLUP of $t_y.$ This model selection step would be expected to increase the variability of this predictor. The BLUP based on the simple linear model in $x_1$ is also calculated. The process is repeated with working models including the correct variance specification ${\text{var}}_M\left(y_i\right)\propto x_{1i}$ and a mis-specification ${\text{var}}_M\left(y_i\right)\propto x_{1i}^2.$

3.4  Simulation results

Tables 3.1-3.4 show the ratios of the prediction MSEs of various BLUP estimators to the GJLB (the right hand side of equation 2.9) for various sample designs and choices of $\mu \left(x_1\right).$ The prediction MSEs are the means over all simulations of ${\left({\widehat{t}}_y-t_y\right)}^2$ so they are with respect to both model and design, as are the results in Theorems 1 and 2.

Table 3.1a evaluates the BLUP corresponding to the lowest-BIC model for a sample size of 500 with correctly specified variances. Nine sample designs are shown corresponding to three choices for $b$ (low/medium/high dependency of the selection probabilities on $x_1)$ and $c$ (no/some/high dependency of selection probabilities on $x_2).$ The period $h$ of the sinusoidal component of $\mu \left(x_1\right)$ is also shown (see equation 3.3). The table shows that:

• The ratio is always either less than or equal to 1 or slightly above 1, consistent with Theorem 2. Its range is 0.815 to 1.086.
• The ratio decreases slightly as $h$ increases. So, when the true $E_M\left(y|x_1\right)$ is close to linear, the model-based BLUP has lower MSE relative to the GJLB, while for more nonlinear models the ratio is closer to 1.
• The ratio depends on $b$ to a degree, although the pattern depends on the other parameters.
• The ratio decreases dramatically as $c$ increases, with reductions of up to 20% from $c\mathrm{<mo>=</mo>}0.$ This shows that the BLUP does much better relative to the GJLB when there is a covariate $x_2$ which is relevant in the design but not relevant to $Y$ so that it can be omitted from the estimation model.

Table 3.1b shows results for a much larger sample size of 25,000 from a population of 100,000. These results were included to see whether the ratios are less than or equal to 1 for large $n$ as predicted by the theory in Section 2. A larger number of simulations were also used for this panel (15,000 rather than 5,000). Results are only shown for $c\mathrm{<mo>=</mo>}0$ since these were the designs with the highest ratios in Table 3.1a. The ratios in 3.1b range from 0.984 to 1.011. The values slightly above 1 may reflect that the working spline model does not perfectly capture the sinusoidal functions used to generate the data.

Table 3.2 shows the same scenarios as Table 3.1a except that the BLUP is based on a mis-specified variance model with ${\sigma }_i^2\propto x_i^2$ (the generating model has ${\sigma }_i^2\propto x_i).$ The ratios of the MSE of the lowest-BIC-model BLUP to the GJLB are on the whole slightly higher than Table 3.1 (generally by less than one percentage point). The ratio is still almost always less than or equal to 1, with a maximum value of 1.089.

Table 3.3 is similar to Table 3.1a except that the sample size is 1,500 rather than 500. The ratios are almost always lower than in Table 3.1a. The maximum ratio is 1.030.

Table 3.4 shows results for the BLUP based on the simple linear model containing only $x_1$ with mis-specified variance. The sample size is 500. When the period is 5, so that the true model is virtually linear in $x_1,$ this BLUP does very well. The ratios are then always below 1.1 and can be as low as 0.790. When the period is 2, there is visible curvature in $E\left(y|x_1\right)$ (as shown in Figure 3.1), but the simple BLUP still does well, with all ratios less than 1.2. However, for periods 0.5 and 1, the ratios are well above 1, with a maximum value of 3.4. This shows the substantial bias of the BLUP when the model is badly mis-specified.

The extent to which the selection probabilities depend on $x_2$ is the major factor determining the ratio of the MSE to the GJLB, as shown by Tables 3.1-3.3. Figure 3.2 shows this phenomenon in more detail for correctly specified variance. The sample size is 1,500 with PP $\sigma$ sampling so that ${\pi }_i\propto x_{1i}^{\text{0}\text{.5}}.$ Values of $c$ (0, 0.25, ..., 3) are on the x-axis and results are shown for different periods $h.$ The figure shows that the ratio is slightly above 1 for $c\mathrm{<mo>=</mo>}0$ and decreases smoothly with $c,$ to about 0.6 when $c\mathrm{<mo>=</mo>}3.$ Higher periods $h$ (reflecting a smoother relationship between $Y$ and $x_1)$ are also associated with lower ratios, but the differences are so small as to be almost indiscernible in Figure 3.2.

• Previous
• Article main page
• Next