6 Simulation study

Jae Kwang Kim and Changbao Wu

Previous | Next

In this section we report results from a simulation study. We consider a synthetic finite population of size $N=$ 2,248 families using a real data set of Statistics Canada's 2000 Family Expenditure Survey for the province of Ontario. For the $i^{\text{th}}$ selected family, the data set contains observations on several variables, including $x_{i1}:$ the number of persons in the family; $x_{i2}:$ the number of children (age $<15);x_{i3}:$ the number of youths (age 15 - 24); $x_{i4}:$ the total annual income after taxes; $y_{i1}:$ the total annual expenditure; $y_{i2}:$ the annual expenditure on clothing; $y_{i3}:$ the annual expenditure on furnishings and equipment.

We consider three population parameters for comparing different versions of replication variance estimators. The first is the population total of overall annual expenditures, i.e., ${\theta }_1=t_{y1}={\displaystyle {\sum }_{i=1}^N{y}_{i1}}.$ The second is the ratio of population totals of expenditures on clothing and on furnishings and equipment, i.e., ${\theta }_2=t_{y2}/t_{y3}=\left({\displaystyle {\sum }_{i=1}^N{y}_{i2}}\right)/\left({\displaystyle {\sum }_{i=1}^N{y}_{i3}}\right).$ Note that ${\theta }_2={\mu }_{y2}/{\mu }_{y3}.$ The third is the population correlation coefficient ${\theta }_3=\rho (y_1,y_2)$ between the overall annual expenditure $(y_1)$ and the annual expenditure on clothing $(y_2).$ For each parameter, several replication variance estimators were evaluated through simulation.

We investigate two approaches of replication variance estimation. For the first one, the initial $L$ sets of the replication weights are constructed using the general method described in Section 2. For the second one, the $L=n$ sets of standard delete-1 jackknife replication weights are used to produce fully efficient variance estimators.

Case I. Unequal probability samples are selected by the Rao-Sampford PPS sampling method (Rao 1965; Sampford 1967), with inclusion probabilities ${\pi }_i$ proportional to the total annual income $x_{i4}.$ One of the attractive features of the Rao-Sampford method is that the second order inclusion probabilities ${\pi }_{ij}$ can be computed exactly. The general procedure described in Section 2 is used to create $L=n$ sets of fully efficient replication weights, and the corresponding variance estimator is denoted as $v_R^{}\text{}.$ Those weights are used as the basis to compute and compare different versions of variance estimators $v_R^{(l)},l=1,2,3,4$  described in Section 3.1, based on a smaller number $L_0$ sets of weights. We restrict $L_0$ to be 25 or 50.

The calibrated replication variance estimator described in Section 3.2 is denoted as $v_C\text{}.$ The initial $L_0$ sets of weights are selected from the original $L$ sets of weights by simple random sampling, and $z_i=(x_{i1},x_{i2},x_{i3},x_{i4},y_{i2},y_{i3}{)}^{\prime }$ is used as calibration variables. Under this setting, the first parameter ${\theta }_1$ is not directly related to $z$ but the second parameter ${\theta }_2$ is defined as a nonlinear but smooth function of $t_z.$ The third parameter $\rho (y_1,y_2)$ is more complex and involves population quantities not included in $t_z.$

Case II. The population of $N=$ 2,248 units is first duplicated 10 times, to create a larger population with 22,480 units. Simple random samples of $n=$ 100, 200 or 400 are selected from the population. The sampling fractions are less than 2%. Under such scenarios the standard $n$ sets of delete-1 jackknife weights provide fully efficient variance estimator $v_J\text{}.$ Let $v_J^{(1)}$ be the variance estimator using $L_0$ sets of weights, randomly selected from the $n$ sets of jackknife weights. Let $v_J^{(C)}$ be the variance estimator using the $L_0$ sets of weights plus calibration over the $z$ variables.

For each simulated sample of size $n$ and a particular population parameter $\theta ,$ we compute design-based estimator $\widehat{\theta }$ and different versions of variance estimators. The process is repeated $B$ times, independently, with $B=$ 5,000 for Case I and $B=$ 10,000 for Case II. The true variance $V=V(\widehat{\theta })$ is approximated by $V\doteq B^{-1}{\displaystyle {\sum }_{b=1}^B\text{}}{({\widehat{\theta }}_b-\theta )}^2,$ where ${\widehat{\theta }}_b$ is calculated from the $b^{\text{th}}$ simulated sample, using another independent $B$ simulated samples. Simulation results show that the bias of $\widehat{\theta }$ is negligible for all three parameters. Performances of a variance estimator $v$ are measured by the simulated coverage probability of the 95% normal theory confidence interval, computed as $\text{CP}=B^{-1}{\displaystyle {\sum }_{b=1}^{B}I[{\widehat{\theta }}_b-1.96{(v_b)}^{1/2}\le \theta \le {\widehat{\theta }}_b+1.96{(v_b)}^{1/2}]},$  the average length of the interval $\text{AL}=B^{-1}{\displaystyle {\sum }_{b=1}^{B}2\times 1.96}{(v_b)}^{1/2},$  and the Relative Root Mean Square Error (RRMSE), computed as $\text{RRMSE}={\{\text{MSE}(v)\}}^{1/2}/V,$ where $v_b$ is the variance estimator $v$ computed from the $b^{\text{th}}$ simulated sample, and $\text{MSE}(v)=B^{-1}{\displaystyle {\sum }_{b=1}^B{(v_b-V)}^2.}$

The simulated coverage probabilities are reported in Tables 6.1 and 6.2. The fully efficient variance estimator $v_R$ and $v_J$ provides good coverage for all scenarios except for $\rho (y_1,y_2)$ with Case I where the coverage is a bit low. The variance estimators $v_R^{(l)},l=1,2,3,4$  and $v_J^{(1)}$ based on $L_0$ sets of weights seem to work for ${\theta }_1,$ to certain degree for ${\theta }_2$ as well, but none is working for ${\theta }_3=\rho (y_1,y_2).$ The calibrated estimator $v_C$ provides satisfactory coverage for all scenarios for Case I. As for the calibrated estimator $v_J^{(C)}$ with Case II, it works very well for ${\theta }_1$ and ${\theta }_2,$ but none are working well for ${\theta }_3=\rho (y_1,y_2).$

 

Table 6.1
Coverage probabilities of 95% confidence intervals (Case I)

$\theta$	$L_0$	$n$	$v_R$	$v_R^{(1)}$	$v_R^{(2)}$	$v_R^{(3)}$	$v_R^{(4)}$	$v_C(\text{AL})$
$t_{y1}$	25	50	93.9	92.9	93.1	92.4	93.1	94.3 (1.03)
		100	94.4	92.0	92.4	91.9	93.0	93.4 (1.01)
		150	95.1	91.5	91.9	92.1	93.2	93.7 (0.99)
	50	100	94.5	93.2	93.2	93.4	93.6	94.1 (1.01)
		150	95.1	93.0	93.3	93.5	93.8	94.5 (0.99)
${\mu }_{y2}/{\mu }_{y3}$	25	50	92.6	91.0	91.2	90.6	91.0	92.9 (1.02)
		100	93.7	91.1	91.2	89.6	90.8	93.7 (1.01)
		150	94.3	91.1	90.7	89.5	90.8	94.3 (1.00)
	50	100	93.6	92.6	92.5	91.9	92.5	93.8 (1.01)
		150	94.2	92.7	92.6	91.9	92.9	94.3 (1.00)
$\rho (y_1,y_2)$	25	50	89.0	85.7	85.6	79.3	81.9	91.3 (1.14)
		100	90.5	85.4	85.3	78.6	81.5	92.1 (1.17)
		150	90.7	84.6	84.5	76.9	81.6	91.9 (1.17)
	50	100	90.4	88.2	88.2	83.2	85.8	92.9 (1.18)
		150	90.7	87.5	87.6	81.7	84.8	93.4 (1.18)
$v_R:$ The fully efficient replication variance estimator (Section 2); $v_R^{(l)},l=$ 1, 2, 3, 4: replication variance estimators based on $L_0$ sets of weights (Section 3.1); $v_C:$ replication variance estimator based on $L_0$ sets of calibrated weights (Section 3.2); $AL:$ average length of the confidence interval relative to the one using $v_R.$

 

It should be noted that the definition of $\rho (y_1,y_2)$ involves population means over three derived variables $y_1^2,y_2^2$ and $y_1{y}_2.$ When those three variables are also included at the calibration stage, in addition to $z,$ the resulting variance estimator is denoted as $v_J^{(C+)}$ for Case II. It turns out that $v_J^{(C+)}$ provides much better results for ${\theta }_3$ and also improved results for ${\theta }_1$ and ${\theta }_2.$

Also included in Tables 6.1 and 6.2 are the average length of the confidence intervals using $v_C,v_J^{(C)}$ and $v_J^{(C+)}.$ The results (AL, in parentheses) are relative to the length of the interval using $v_R$ (Table 6.1) or $v_J$ (Table 6.2), with a value (say) 1.05 indicating 5% increase in length. It can be seen that the calibrated variance estimators produce confidence intervals which are either comparable in length to the intervals using $v_R(v_J)$ or slightly wider, depending on the parameter and/or sample sizes.

 

Table 6.2
Coverage probabilities of 95% confidence intervals (Case II)

$\theta$	$L_0$	$n$	$v_J$	$v_J^{(1)}$	$v_J^{(C)}(\text{AL})$	$v_J^{(C+)}(\text{AL})$
$t_{y1}$	25	100	94.4	92.0	94.4 (1.02)	94.9 (1.07)
		200	95.0	92.4	95.0 (1.01)	95.2 (1.03)
		400	95.3	92.5	95.1 (0.99)	95.3 (1.01)
	50	100	94.1	93.1	94.2 (1.02)	94.8 (1.07)
		200	94.7	93.3	94.8 (1.01)	95.0 (1.04)
		400	94.7	93.4	94.5 (0.99)	94.8 (1.02)
${\mu }_{y2}/{\mu }_{y3}$	25	100	92.6	87.3	92.6 (1.05)	93.3 (1.11)
		200	93.6	86.8	93.3 (1.02)	93.7 (1.07)
		400	94.1	86.8	93.8 (0.99)	94.1 (1.04)
	50	100	92.8	90.3	92.9 (1.06)	94.2 (1.11)
		200	93.9	89.8	93.8 (1.03)	94.3 (1.08)
		400	94.1	89.6	93.8 (1.00)	94.1 (1.05)
$\rho (y_1,y_2)$	25	100	92.5	78.0	89.4 (1.06)	91.7 (1.09)
		200	92.7	72.3	86.3 (1.00)	91.4 (1.05)
		400	93.2	71.2	84.5 (0.95)	92.1 (1.04)
	50	100	92.2	84.5	92.2 (1.09)	92.5 (1.11)
		200	92.8	80.5	90.3 (1.05)	92.2 (1.08)
		400	93.1	77.4	88.1 (1.00)	92.6 (1.06)
$v_J:$ The delete-1 jackknife variance estimator; $v_J^{(1)}:$ replication variance estimator based on $L_0$ sets of jackknife weights; $v_J^{(C)}:$ replication variance estimator based on $L_0$ sets of calibrated jackknife weights; $v_J^{(C+)}:$ replication variance estimator based on $L_0$ sets of calibrated jackknife weights, with added variables for weight-calibration; $AL:$ average length of the confidence interval relative to the one using $v_J.$

 

The relative root mean square errors (RRMSE) of variance estimators are presented in Tables 6.3 and 6.4. The results seem to depend not only on the parameter and its estimator but also the sampling design and the replication method. For Case I, the variance estimator $v_C,$ which is of primary interest, is more stable than $v_R$ for ${\theta }_1,$ almost the same for ${\theta }_2,$ and is less stable for ${\theta }_3.$ Because $y_{i1}$ is well explained by $z_i,v_C$ is quite efficient for estimating the variance of ${\widehat{\theta }}_1={\widehat{t}}_{y1}.$ For Case II, $v_J^{(C)}$ and $v_J^{(C+)}$ are similar to each other but both are less stable than $v_J.$

 

Table 6.3
Relative Root Mean Square Errors (RRMSE, Case I)

$\theta$	$L_0$	$n$	$v_R$	$v_R^{(1)}$	$v_R^{(2)}$	$v_R^{(3)}$	$v_R^{(4)}$	$v_C$
$t_{y1}$	25	50	1.84	2.76	2.24	1.99	1.86	1.43
		100	1.32	2.34	1.67	1.89	1.40	0.83
		150	1.19	1.99	1.34	1.37	1.46	0.87
	50	100	1.32	1.91	1.69	1.63	1.35	0.92
		150	1.19	1.81	1.50	1.62	1.24	0.73
${\mu }_{y2}/{\mu }_{y3}$	25	50	0.72	0.89	0.88	1.07	0.89	0.74
		100	0.45	0.78	0.77	0.99	0.72	0.46
		150	0.41	0.93	0.87	1.01	0.74	0.41
	50	100	0.46	0.60	0.60	0.77	0.56	0.46
		150	0.41	0.70	0.68	0.70	0.54	0.41
$\rho (y_1,y_2)$	25	50	0.65	0.79	0.83	1.45	1.26	0.96
		100	0.65	1.12	1.16	2.20	1.37	1.24
		150	0.59	1.29	1.34	2.27	1.43	1.50
	50	100	0.65	0.84	0.88	1.63	0.95	1.03
		150	0.59	0.88	0.94	1.48	1.05	1.12
$v_R:$ The fully efficient replication variance estimator (Section 2); $v_R^{(l)},l=$ 1, 2, 3, 4: replication variance estimators based on $L_0$ sets of weights (Section 3.1); $v_C:$ replication variance estimator based on $L_0$ sets of calibrated weights (Section 3.2).

 

Table 6.4
Relative Root Mean Square Errors (RRMSE, Case II)

$\theta$	$L_0$	$n$	$v_J$	$v_J^{(1)}$	$v_J^{(C)}$	$v_J^{(C+)}$
$t_{y1}$	25	100	0.29	0.56	0.56	0.66
		200	0.21	0.57	0.47	0.53
		400	0.15	0.56	0.20	0.41
	50	100	0.29	0.41	0.50	0.57
		200	0.21	0.41	0.39	0.44
		400	0.15	0.41	0.17	0.35
${\mu }_{y2}/{\mu }_{y3}$	25	100	0.78	1.58	1.90	1.98
		200	0.56	1.44	1.41	1.57
		400	0.39	1.54	0.87	1.40
	50	100	0.81	1.12	1.61	1.67
		200	0.57	1.10	1.22	1.32
		400	0.38	1.04	0.72	1.00
$\rho (y_1,y_2)$	25	100	1.02	2.43	2.71	2.72
		200	0.74	2.44	2.64	2.57
		400	0.47	2.51	2.52	2.42
	50	100	1.04	1.64	1.97	2.12
		200	0.71	1.75	2.01	1.96
		400	0.48	1.76	1.83	1.73
$v_J:$ The delete-1 jackknife variance estimator; $v_J^{(1)}:$ replication variance estimator based on $L_0$ sets of jackknife weights; $v_J^{(C)}:$ replication variance estimator based on $L_0$ sets of calibrated jackknife weights; $v_J^{(C+)}:$ replication variance estimator based on $L_0$ sets of calibrated jackknife weights, with added variables for weight-calibration.

 

Previous | Next