4 Simulation results for $H=1$ 

Jun Shao, Eric Slud, Yang Cheng, Sheng Wang, and Carma Hogue

Previous | Next

Large sample theory as presented above is not adequate to tell whether the asymptotic results adequately describe the behavior of the estimators ${\widehat{Y}}_{\text{reg}\mathrm{,1}},$ ${\widehat{Y}}_{\text{reg}\mathrm{,2}}$ and ${\widehat{Y}}_{\text{dec}}$ and their variance estimators in moderate samples, or whether ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ and ${\widehat{Y}}_{\text{dec}}$ever provide useful Mean-Squared-Error improvements in moderate sized samples. We present some simulation results to study these questions, as well as the small-sample issues arising in applying these methods in the context of the ASPEP survey.

In the simulations, values in the frame population $U$ are either generated under some model or are taken from the 2002 and 2007 Government censuses with 2007 ASPEP sample weights. The first set of simulations (reported in Tables 4.1-4.6) summarizes average behavior over many model-generated frame populations. In the second set of artificial-data simulations, summarized in Table 4.8, the frame population remains fixed throughout the simulation. All frame populations consist of a single stratum $\left(H=1\right)$ broken into two substrata $\left(j=1,2\right)$ according as a size variable falls below or above a specific quantile, usually the 0.8 quantile. Sampling from the frame populations is done PPS with-replacement in all simulations in this section.

4.1 Small sample considerations

Before proceeding to describe the simulations, we discuss some special features of PPS with-replacement (PPSWR) sampling which, when done in settings with small samples and unbalanced size variables, requires special computational handling. Numerically erratic results can arise when the small drawn samples are used stratumwise and then bootstrapped to estimate variances.

The weights ${\pi }_i^{-1}=1/\left(n_j{p}_{ij}\right)$ in PPSWR are all greater than 1 only when the single-draw probabilities $p_{ij}=z_i/{\displaystyle {\sum }_{i^{\prime }\in U_j}}{z}_{i^{\prime }}$ are all below $1/n_j.$ To avoid anomalous small-sample results, and to maintain the relevance of PPSWR designs in imitating PPS without-replacement designs, any units $i\in U_j$ with $n_j{p}_{ij}\ge 1$ are made self-representing (SR), i.e., are sampled with certainty but only once, and if there are $k_j$ such units, then the probabilities $\left\{p_{ij}:i\in U_j\mathrm{,}{n}_j{p}_{ij}\text{<}1\right\}$ are renormalized to draw a size $n_j-k_j$ PPSWR sample. If any of the remaining renormalized probabilities are $\ge 1/(n_j-k_j),$ then their units also become self-representing and a new renormalization is done. This is repeated as often as necessary. Thus, small samples with very unequal size-variable distributions may not be compatible with PPSWR sampling, a condition arising in some of the real-data ASPEP cases considered below.

Although a different choice could have been made, we conform with ASPEP practice in including all SR units in the fitting of the survey-weighted regression estimators ${\widehat{\beta }}_2$ and $\widehat{\beta }.$ However, with this choice, PPSWR sampling followed by bootstrap resampling of small samples can lead to extremely erratic behavior, which must be recognized in summarizing the behavior of bootstrap variance estimators. The problem is that when a small number $m$ of non-self-representing items are sampled PPSWR, in addition to a set of SR items, and then bootstrapped, the probability can be surprisingly large that there is only one unique non-SR item in the bootstrap sample, leading to very high bootstrap variability. This phenomenon was observed in the simulations reported below, with large-size substratum containing 20 or fewer elements and very skewed size-variables, either in the cases with lognormal or ASPEP $x_i$variables.

4.2 Artificial model-generated data

All of the artificial frame populations were generated with $N=$ 2,000 iid triples $(x_i\mathrm{,}{y}_i\mathrm{,}{z}_i)$ satisfying (2.7), for $U_1$ consisting of the $N_1=$ 1,600 for which $x_i$ fell below their empirical 80'th percentile $c=\left(x_{\left(1,600\right)}+x_{\left(1,601\right)}\right)/2,$ and $U_2$ consisting of the other 400 indices. In most cases, $z_i$ were generated as $N\left(30+x_i\mathrm{,100}\right)$ variates conditioned to be positive (which required occasional re-simulation in the lognormal- $x_i$ models below) and were conditionally independent of $y_i$ given $x_i.$ (However, in some cases, unweighted samples were drawn by taking $z_i$ identically equal.) PPS with-replacement stratified samples of sizes $\left(n_1\mathrm{,}{n}_2\right)=\left(\mathrm{100,50}\right),\left(\mathrm{100,20}\right),$ or $\left(\mathrm{50,20}\right)$ were drawn in successive simulation runs, with size-variables $z_i,$from the same frame.

The models generating $(x_i\mathrm{,}{y}_i)$ are indexed as follows. In those with prefix $M1,$ the predictors $x_i$ are Gamma $(4,0.1)$ distributed, with 0.8 quantile 55.2, while in the models $M2,$ the $x_i$ variables are Lognormal(1,6.25), with 0.8 quantile 22.3. The $M1$ populations, and the $M2$ models with suffix $E,$ have conditional variance 100 for $y_i$ given $x_i,$ while the $M2$ models without suffix $E$ have conditional variance $20x_i.$ Conditional means $E\left(y_i|x_i\right)$ are all linear, equal to $20+1.5x_i$ in models indexed $H0$ and to $20+x_i+0.5\left(x_i-c\right)I_{\left[j=2\right]}$ within the substratum $U_j$ in models $H1.$ The intercepts of the regression models are so chosen that whether or not the slopes are the same, the lines intersect at $x=c$ (see the discussion in Section 1). Table 4.1 exhibits the average and standard deviation for the totals $Y$ generated from the frame-population attributes ${\left\{y_i\right\}}_{i=1}^{\text{2,000}}$ under the various models. The variates $x_i$ as well as the totals $Y$are much longer-tailed under the Lognormal models.

Table 4.1
Means and standard deviations for totals Y under simulation models.
Table summary
This table displays the results of Means and standard deviations for totals Y under simulation models.. The information is grouped by Model (appearing as row headers), Gamma and Lognormal (appearing as column headers).

Model	Gamma		Lognormal
	M1.H0	M1.H1	M2.H0	M2.H0E	M2.H1	M2.H1E
E(Y)	160,000	123,177	225,603	225,603	173,485	173,485
SD(Y)	1,414.20	653.5	94,380	94,368	62,362	62,344

Simulated population models

$M1.H0$ :$x_i\sim \text{Gamma}\left(\mathrm{4,0.1}\right)$ (shape parameter 4, scale 10),
$y_i\sim N\left(20+1.5x_i\mathrm{,100}\right)$ (variance 100), all $i\in U.$

$M1.H1$ :$x_i\sim \text{Gamma}\left(\mathrm{4,0.1}\right),$ $y_i\sim N\left(20+x_i+0.5\left(x_i-c\right)I_{\left[x_i\ge c\right]}\mathrm{,100}\right),$ all $i.$

$M2.H0$ :$\log \left(x_i\right)\sim N\left(\mathrm{1,6.25}\right)\mathrm{,}{y}_i\sim N\left(20+1.5x_i\mathrm{,20}{x}_i\right)\mathrm{,}\text{all}i.$

$M2.H0E$ :$\log \left(x_i\right)\sim N\left(\mathrm{1,6.25}\right)\mathrm{,}{y}_i\sim N\left(20+1.5x_i\mathrm{,100}\right)\mathrm{,}\text{all}i.$

$M2.H1$ :$\log \left(x_i\right)\sim N\left(\mathrm{1,6.25}\right)\mathrm{,}{y}_i\sim N\left(20+x_i+0.5\left(x_i-c\right)I_{\left[x_i\ge c\right]}\mathrm{,20}{x}_i\right)\mathrm{,}\text{all}i.$

$M2.H1E$ :$\log \left(x_i\right)\sim N\left(\mathrm{1,6.25}\right)\mathrm{,}{y}_i\sim N\left(20+x_i+0.5\left(x_i-c\right)I_{\left[x_i\ge c\right]}\mathrm{,100}\right)\mathrm{,}\text{all}i.$

The simulation and bootstrap results in Tables 4.2-4.5 were generated by the following design and reporting scheme. For each population type, 60 distinct frame populations were generated, and 50 independent sampling experiments were conducted with each of those. In those cases where results of weighted and unweighted sampling were compared, these samples were drawn independently from the same set of 60 frame populations. Thus there were 3,000 independent replications for Monte Carlo averaging of statistical results, done for each of three different stratified sample sizes, and 400 bootstrap iterations were performed for each such generated sample.

Table 4.2
Empirical and estimated SD’s and CI coverage, from model M1 simulations
Table summary
This table displays the results of Empirical and estimated SD’s and CI coverage. The information is grouped by Sizes (appearing as row headers), Stat, M1.H0 and M1.H1 (appearing as column headers).

Sizes	Stat	M1.H0			M1.H1
		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	${\widehat{Y}}_{\text{dec}}$	${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	${\widehat{Y}}_{\text{dec}}$
100,50	${\text{SD}}_{MC}$	1,785.5	1,794.3	1,788.0	1,817.6	1,773.5	1,774.4
	${\widehat{\text{SD}}}_S$	1,757.1	1,751.5	1,755.6	1,794.6	1,735.2	1,735.8
	${\widehat{\text{SD}}}_B$	1,752.4	1,762.0	1,758.4	1,788.1	1,742.9	1,747.0
	${\text{CP}}_S$	94.47	94.37	94.50	93.93	93.73	93.77
	${\text{CP}}_B$	94.60	94.53	94.67	93.93	94.03	94.07
100,20	${\text{SD}}_{MC}$	1,930.0	1,944.8	1,934.0	2,008.4	1,944.4	1,960.4
	${\widehat{\text{SD}}}_S$	1,888.3	1,876.6	1,884.1	1,944.4	1,861.0	1,866.5
	${\widehat{\text{SD}}}_B$	1,878.8	1,901.4	1,895.8	1,936.1	1,885.6	1,897.9
	${\text{CP}}_S$	94.20	93.83	94.13	93.53	93.20	93.07
	${\text{CP}}_B$	93.80	94.00	93.97	93.60	93.83	93.97
50,20	${\text{SD}}_{MC}$	2,583.5	2,610.7	2,593.5	2,591.3	2,522.8	2,535.4
	${\widehat{\text{SD}}}_S$	2,509.2	2,490.8	2,505.1	2,562.2	2,465.0	2,474.5
	${\widehat{\text{SD}}}_B$	2,498.5	2,538.0	2,522.9	2,550.3	2,508.5	2,525.6
	${\text{CP}}_S$	93.70	93.13	93.57	93.97	93.63	93.43
	${\text{CP}}_B$	93.63	93.73	93.87	93.83	93.77	94.10

 

Table 4.3
Empirical and estimated SD’s and CI coverage, from model M2 simulations
Table summary
This table displays the results of Empirical and estimated SD’s and CI coverage. The information is grouped by Sizes (appearing as row headers), Stat, M2.H0 and M2.H1 (appearing as column headers).

Sizes	Stat	M2.H0			M2.H1
		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	${\widehat{Y}}_{\text{dec}}$	${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	${\widehat{Y}}_{\text{dec}}$
100,50	${\text{SD}}_{MC}$	3,400.1	3,475.4	3,406.8	3,481.9	3,483.8	3,482.2
	${\widehat{\text{SD}}}_S$	3,420.6	3,400.0	3,417.0	3,537.8	3,405.0	3,463.7
	${\widehat{\text{SD}}}_B$	3,590.0	3,715.2	3,623.4	3,852.0	3,921.9	3,898.4
	${\text{CP}}_S$	95.10	93.43	94.83	95.03	93.40	94.13
	${\text{CP}}_B$	95.67	95.77	95.77	95.63	95.77	95.70
100,20	${\text{SD}}_{MC}$	5,655.2	6,184.0	5,698.6	5,853.0	6,181.1	5,955.6
	${\widehat{\text{SD}}}_S$	5,644.9	5,575.7	5,640.9	5,798.3	5,587.3	5,697.3
	${\widehat{\text{SD}}}_B$	5,565.1	6,687.3	5,857.8	5,907.8	6,838.0	6,466.6
	${\text{CP}}_S$	93.83	88.47	93.40	92.77	88.30	90.70
	${\text{CP}}_B$	92.33	93.67	93.37	92.63	94.33	94.17
50,20	${\text{SD}}_{MC}$	5,773.2	6,319.2	5,833.9	5,934.2	6,230.6	6,009.8
	${\widehat{\text{SD}}}_S$	5,800.2	5,677.2	5,785.8	6,012.6	5,755.4	5,919.2
	${\widehat{\text{SD}}}_B$	5,728.5	6,825.2	6,086.0	6,102.2	6,978.1	6,522.1
	${\text{CP}}_S$	94.60	88.67	93.97	94.07	89.37	92.27
	${\text{CP}}_B$	93.40	94.23	94.27	93.47	95.03	94.80

 

Table 4.4
SD’s for ${\widehat{Y}}_{\text{HT}}$ vs. ${\widehat{Y}}_{\text{dec}},$ and coverage for Bootstrap Percentile Confidence Intervals for ${\widehat{Y}}_{\text{dec}},$ for $\tau =0.05$ vs. $0.20$, for models M1 and M2, H0 and H1
Table summary
This table displays the results of SD’s for XXX vs. XXX and coverage for Bootstrap Percentile Confidence Intervals for XXX for XXX vs. 0.20. The information is grouped by Model (appearing as row headers), Samples and xxx (appearing as column headers).

Model	Samples	${\widehat{Y}}_{\text{dec}}\mathrm{,}\tau =0.05$		${\widehat{Y}}_{\text{HT}}$	${\widehat{Y}}_{\text{dec}}\mathrm{,}\tau =0.20$
		${\text{SD}}_{MC}$	${\text{CP}}_{BP}$	${\text{SD}}_{\text{HT}}$	${\text{SD}}_{MC}$	${\text{CP}}_{BP}$
M1.H0	100,50	1,788.0	94.23	2,774.0	1,745.5	94.60
	100,20	1,934.0	93.50	3,032.6	1,915.9	94.10
	50,20	2,593.5	93.17	3,000.7	2,500.1	94.43
M1.H1	100,50	1,774.4	93.70	2,387.3	1,737.3	94.43
	100,20	1,960.4	93.27	2,678.9	1,948.0	93.23
	50,20	2,535.4	93.90	3,035.0	2,509.8	94.23
M2.H0	100,50	3,406.8	95.20	4,160.0	3,398.8	94.83
	100,20	5,698.6	91.13	6,720.2	5,705.7	92.57
	50,20	5,833.9	92.60	7,080.0	5,979.8	92.17
M2.H1	100,50	3,482.2	95.13	4,393.6	3,423.9	94.03
	100,20	5,955.6	92.07	7,413.1	5,917.3	92.40
	50,20	6,009.8	92.33	7,840.4	6,105.6	92.17

 

Table 4.5
Comparisons of SD estimates and CI coverage for H0 and H1 for three lognormal settings, weighted (W) and unweighted (U) within M2, and weighted (E) within M2.E. CI % coverages are given for both the Bootstrap SD and Percentile Intervals.
Table summary
This table displays the results of Comparisons of SD estimates and CI coverage for H0 and H1 for three lognormal settings. The information is grouped by Model (appearing as row headers), Size, Stat, SD and xxx (appearing as column headers).

Model	Size	Stat	SD	${\widehat{\text{SD}}}_S$	${\widehat{\text{SD}}}_B$	${\text{CP}}_S$	${\text{CP}}_B$	${\text{CP}}_{BP}$
H0.W	100,50	${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	3,400.1	3,420.6	3,590.0	95.10	95.67	94.93
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	3,475.4	3,400.0	3,715.2	93.43	95.17	95.33
		${\widehat{Y}}_{\text{dec}}$	3,406.8	3,417.0	3,623.4	94.83	95.77	95.20
H0.U		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	5,481.6	3,674.8	5,571.9	81.43	93.50	92.07
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	5,782.8	3,646.6	6,076.3	80.13	93.67	91.90
		${\widehat{Y}}_{\text{dec}}$	5,525.5	3,669.0	5,726.8	81.07	93.83	92.20
H0.E		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	1,888.8	1,930.1	1,904.7	94.73	94.53	94.23
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	1,888.6	1,911.1	1,893.2	94.43	94.30	94.20
		${\widehat{Y}}_{\text{dec}}$	1,892.9	1,926.5	1,905.0	94.67	94.57	94.20
H0.W	50,20	${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	5,773.2	5,800.2	5,728.5	94.60	93.40	92.00
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	6,319.2	5,677.2	6,825.2	88.67	94.23	92.60
		${\widehat{Y}}_{\text{dec}}$	5,833.9	5,785.8	6,086.0	93.97	94.27	92.60
H0.U		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	10,000.3	5,136.5	9,905.6	71.10	90.73	89.80
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	11,192.8	5,085.0	12,806.8	68.70	92.90	89.37
		${\widehat{Y}}_{\text{dec}}$	10,134.1	5,120.7	11,245.9	70.73	92.37	90.27
H0.E		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	2,811.4	2,831.6	2,769.5	94.13	94.00	93.93
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	2,811.9	2,753.8	2,741.1	93.47	93.77	93.30
		${\widehat{Y}}_{\text{dec}}$	2,817.4	2,821.8	2,777.0	93.83	93.90	93.77
H1.W	100,50	${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	3,481.9	3,537.8	3,852.0	95.03	95.63	95.27
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	3,483.8	3,405.0	3,921.9	93.40	95.77	95.10
		${\widehat{Y}}_{\text{dec}}$	3,482.2	3,463.7	3,898.4	94.13	95.70	95.13
H1.U		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	5,631.4	3,774.8	5,614.6	80.90	92.33	91.07
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	5,838.3	3,699.6	6,010.5	79.13	92.73	91.37
		${\widehat{Y}}_{\text{dec}}$	5,727.0	3,732.8	5,870.5	80.40	92.93	91.63
H1.E		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	2,005.5	2,094.2	2,019.1	95.60	94.97	94.60
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	1,909.9	1,908.2	1,892.5	94.83	94.77	94.17
		${\widehat{Y}}_{\text{dec}}$	1,931.9	1,941.7	1,934.6	94.97	95.20	94.83
H1.W	50,20	${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	5,934.2	6,012.6	6,102.2	94.07	93.47	91.97
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	6,230.6	5,755.4	6,978.1	89.37	95.03	92.23
		${\widehat{Y}}_{\text{dec}}$	6,009.8	5,919.2	6,522.1	92.27	94.80	92.33
H1.U		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	9,315.8	5,350.9	10,040.0	74.17	93.10	90.57
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	10,583.8	5,229.6	12,476.8	71.23	94.57	90.87
		${\widehat{Y}}_{\text{dec}}$	9,989.6	5,295.4	11,479.5	72.53	94.33	91.47
H1.E		${\widehat{Y}}_{\text{reg}\mathrm{,1}}$	3,096.1	3,137.7	2,795.6	94.63	93.43	93.37
		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$	2,880.6	2,766.8	2,745.7	93.10	93.40	93.47
		${\widehat{Y}}_{\text{dec}}$	2,977.3	2,929.2	2,882.0	93.77	93.77	93.77

We calculated the following quantities for each combination of model, weighting, and sample size: the percentage biases of ${\widehat{Y}}_{\text{reg}\mathrm{,1}},$ ${\widehat{Y}}_{\text{reg}\mathrm{,2}},$ ${\widehat{Y}}_{\text{dec}}$ (with $\tau =$ 0.05 in all tables except Table 4.4, and $\tau =$ 0.05 or 0.20 in Table 4.4) as estimators of $Y;$ the Monte Carlo standard deviations (SD), ${\text{SD}}_{MC},$ of these three estimators; the estimated SD's of the estimators, respectively using the substitution $\left({\widehat{\text{SD}}}_S\right)$ and bootstrap $\left({\widehat{\text{SD}}}_B\right)$ SD estimators described in Section 3; the coverage probability, $\text{C}{P}_u,$ of the nominal 95% confidence intervals for $Y:$ $\widehat{Y}\pm 1.960\cdot {\widehat{\text{SD}}}_u,$ where $\widehat{Y}$ is one of the three estimators of $Y,$ and $u=S$ or $B;$ and the bootstrap percentile confidence intervals (and their coverage percentages $\text{C}{P}_{BP}$ ) obtained from the empirical 0.025 and 0.975 quantiles of the (400) bootstrapped values of each of the three estimators $\widehat{Y}$ of $Y.$ In addition, we calculated empirical biases of the Horvitz-Thompson estimates ${\widehat{Y}}_{\text{HT}}$ in (1.1) and their empirical standard deviations ${\text{SD}}_{\text{HT}}$ . (Of these calculated quantities, only the biases are not shown, since all of the biases were well below 0.5% except in the model $M2\text{.}H1\text{.}U,$ and even there the largest magnitude of bias was about 1%.) Two further statistics, computed and displayed in Table 4.6 for each of the estimators $\widehat{Y}$ of $Y,$ are the standard errors across randomly generated frame populations of the Monte Carlo and Bootstrap within-population estimated SD's of estimators $\widehat{Y}.$

Table 4.6
Cross-population Standard errors of Empirical and Bootstrap SD’s estimated for the estimators ${\widehat{Y}}_{\text{reg}\mathrm{,1}}\mathrm{,}{\widehat{Y}}_{\text{reg}\mathrm{,2}},$ and ${\widehat{Y}}_{\text{dec}},$ for selected models and weighting.
Table summary
This table displays the results of Cross-population Standard errors of Empirical and Bootstrap SD’s estimated for the estimators XXX and XXX for selected models and weighting.. The information is grouped by Model (appearing as row headers), Sizes and xxx (appearing as column headers).

Model	Sizes	${\widehat{Y}}_{\text{reg}\mathrm{,1}}$		${\widehat{Y}}_{\text{reg}\mathrm{,2}}$		${\widehat{Y}}_{\text{dec}}$
		SD	${\widehat{\text{SD}}}_B$	SD	${\widehat{\text{SD}}}_B$	SD	${\widehat{\text{SD}}}_B$
M1.H0	100,50	198	35	196	35	197	35
	50,20	210	52	208	51	210	51
M1.H1	100,50	204	39	183	40	184	41
	50,20	319	57	298	62	302	62
M2.H0	100,50	404	345	450	383	405	351
	50,20	825	518	1,075	916	889	631
M2.H0.E	100,50	187	49	185	45	184	47
	50,20	294	85	293	71	298	82
M2.H1	100,50	409	409	410	421	408	414
	50,20	767	624	946	929	841	730
M2.H1.E	100,50	208	59	196	46	204	50
	50,20	258	141	261	82	239	102
M2.H1.U	100,50	1,676	1,351	1,773	1,539	1,726	1,467
	50,20	2,397	2,543	3,425	3,454	3,102	3,159

4.3 Real government-census data

Our simulations based on repeated sampling from real-data frames rely on a national state-wise dataset assembled by Yang Cheng. For the ASPEP survey of governments for sample year 2007, which was also a census year, the ASPEP frame is the same as the 2007 Census of Governments file. Our dataset consists of the 2002 and 2007 ASPEP variable values (full- and part-time employees, payroll and hours) derived from the censuses in those years, plus the 2007 sample weights and in-sample indicators for ASPEP. Weights equal to 1 imply that governmental units were self-representing (SR), in the sense that they were chosen for inclusion with certainty in ASPEP. The size-variable $z_i$ for PPS sampling within ASPEP is the sum of full- and part-time payroll from the most recent census, so we restrict attention to the 53,402 governmental units in the file for which this variable was positive. Table 4.7 gives the subcounty and special-district governmental types (the only ones that are subdivided into Small and Large unit substrata) in nine selected states, giving also the SR counts and numbers sampled in 2007. As mentioned in subsection 4.1, the final SR count for PPS with-replacement sampling can exceed the number of units initially chosen for certain inclusion, and the larger numbers, corresponding to the sample size actually drawn in 2007, are shown in the SR columns of Table 4.7. Inspection of this Table shows that several of the state by type combinations either have no population in a substratum or have too few non-SR units to be useful in simulating repeated samples. We take 15 as a rule-of-thumb minimum for the number of non-SR units, and suggest that substratum pairs with fewer non-SR units in the large-unit stratum should be collapsed without recourse to the decision-based strategy studied in this paper.

Table 4.7
Census population, ASPEP sample sizes and SR counts of Subcounty and Special-District governmental units by substratum in 2007, for 9 selected states.
Table summary
This table displays the results of Census population. The information is grouped by (appearing as row headers), Subcounty , Special District , Small and Large (appearing as column headers).

	Subcounty					Special District
	Small		Large			Small		Large
	Pop	Samp	Pop	Samp	SR	Pop	Samp	Pop	Samp	SR
AL	378	15	55	45	26	0	0	400	102	64
CA	0	0	475	104	86	1,595	39	107	107	107
CO	0	0	265	34	18	627	16	65	55	33
FL	317	16	81	54	36	0	0	330	48	24
GA	461	17	49	36	20	0	0	293	70	32
MO	980	25	101	101	101	799	27	106	66	42
NY	1,473	25	69	69	69	606	16	33	23	4
PA	2,409	55	123	81	31	921	21	37	37	37
WI	1,702	36	129	71	44	281	16	61	40	20

For nine government-by-type combinations with 15 or more non-SR units and at least 17 non-sampled non-SR large-substratum units (except for AL, CO, and GA for which there were respectively 9, 10, and 11 non-sampled non-SR units), Table 4.8 displays results for the decision-based estimators and variance estimates in substratum pairs. In each of the state-type combinations, 3,000 stratified PPSWR samples of the indicated sizes were drawn from the ASPEP and government census frame described above, with $x_i$ and $y_i$ respectively the full-time payroll amount for the governmental unit as recorded in the 2002 and 2007 governmental censuses, and $z_i$ the total (full-time plus part-time) payroll in 2002. Within each simulated sample, the estimators ${\widehat{Y}}_{\text{reg}\mathrm{,1}}\mathrm{,}{\widehat{Y}}_{\text{reg}\mathrm{,2}}\mathrm{,}{\widehat{Y}}_{\text{dec}}$ were calculated, and the empirical variances estimated. The variance of ${\widehat{Y}}_{\text{dec}}$ was also estimated by the substitution formula and bootstrap methods as in the artificial-data simulations. (But note that, as described above, the bootstrap samples were drawn only from the non-SR units in each substratum sample.) The results are shown in Table 4.8. The relative efficiencies between the combined and separate stratified regression estimators can be gleaned from the corresponding ratio of SD's given in column 5 of the table. The remaining SD's shown are the empirical, substitution, and bootstrap SD estimators of ${\widehat{Y}}_{\text{dec}}.$

Table 4.8
Summary of repeated-sampling simulations from ASPEP 2007 frame. Total full-time pay (Y) given in multiple of $100 million, and estimated SD’s of ${\widehat{Y}}_{\text{dec}}$ given in columns 6-8 in units of $1 million. ${\text{SD}}_1/{\text{SD}}_2$ in column 5 is ratio of empirical SD of ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ over that of ${\widehat{Y}}_{\text{reg}\mathrm{,2}}$
Table summary
This table displays the results of Summary of repeated-sampling simulations from ASPEP 2007 frame. Total full-time pay (Y) given in multiple of $100 million. The information is grouped by State (appearing as row headers), Stratum, Y, Size and xxx (appearing as column headers).

State	Stratum	Y	Size	${\text{SD}}_1/{\text{SD}}_2$	$\widehat{\text{SD}}$	${\widehat{\text{SD}}}_S$	${\widehat{\text{SD}}}_B$
AL	SubCty	1.2	25,46	2.14	4.90	3.67	5.71
CA	SpcDst	4.3	30,90	0.98	29.4	21.2	26.8
CO	SpcDst	0.6	25,55	1.14	3.77	2.58	3.00
FL	SubCty	4.3	25,54	1.16	11.9	9.4	12.2
GA	SubCty	1.5	25,38	1.15	4.38	3.26	4.88
MO	SpcDst	0.6	40,70	2.13	2.99	2.20	2.99
NY	SubCty	23.6	35,52	1.53	13.6	12.0	14.1
PA	SubCty	3.0	40,70	1.12	7.28	5.79	7.60
WI	SubCty	1.4	40,70	2.06	5.00	4.45	5.17

4.4 Discussion of simulation results

The following is a summary and interpretation of the results in the Tables, as well as of other results not shown.

(I) Many of the artificial-data simulations serve to confirm the large-sample theoretical results of the Theorems. It has already been mentioned that in Tables 4.2 and 4.3 the biases for all three $Y­$ estimators $\left({\widehat{Y}}_{\text{reg}\mathrm{,1}}\mathrm{,}{\widehat{Y}}_{\text{reg}\mathrm{,2}}\mathrm{,}{\widehat{Y}}_{\text{dec}}\right)$ are generally small. Within Table 4.2, referring to models with predictors and weights related to the Gamma distribution in models $M1,$ the substitution and bootstrap variance estimators for each $Y­$ estimator are quite accurate and close to one another, and the confidence intervals all have close to nominal coverage. Under both $M1.H0$ and $M1.H1,$ there is a tendency with smaller $n_2$ sample size for the ${\widehat{\text{SD}}}_S$ and ${\widehat{\text{SD}}}_B$ estimators to be slight underestimates of the actual or empirical SD's, but ${\widehat{\text{SD}}}_B$ seems to track SD more closely than ${\widehat{\text{SD}}}_S$ for ${\widehat{Y}}_{\text{reg}\mathrm{,2}}$ and ${\widehat{Y}}_{\text{dec}}.$

(II) The lognormal $x_i$ values in models $M2$ are much more dispersed and skewed than the values in $M1,$ but the simulation results still support the asymptotic theory when $n_2=50,$ although not when $n_2=20.$ The substitution-estimator based confidence intervals for $Y$ in terms of ${\widehat{Y}}_{\text{reg}\mathrm{,2}}$ have coverage probability far too small when the substitution variance estimator is used. In Table 4.3, for each type of $Y­$estimator there is a pronounced tendency for the substitution variance estimator to underestimate the true (empirical) variance, and for the bootstrap estimator to overestimate.

Table 4.5 clarifies that the extreme behavior of variance estimators under models $M2$ occurs partly because the predictors and $y_i$ are dispersed and skewed, and partly because the size-variable used in PPS weighting shares these properties. The cases with suffix $W$ in this Table are the same as in Table 4.3. The cases with suffix $E$ have $\left(x_i\mathrm{,}{z}_i\right)$ the same as in Table 4.3, but the conditional variances of $y_i$ given $x_i$ have the constant value of 100; and with this change, the erratic behavior of SD estimators disappears. However, when the conditional $y_i$ variances are as in the basic model $M2$ but the PPSWR sampling is done unweighted, i.e., with all $z_i$ replaced by 1, the empirical and bootstrap SD estimators track each other and are very large, while the substitution variance estimator is too low by dramatic factors of $1/2$ to $3/4.$ This weird phenomenon applies equally to all three $Y$ estimators. (However, an unweighted-sampling variant in model $M1$does not materially change the results from those shown in Table 4.2.)

(III) One objective of the simulations was to learn whether there can ever be any Mean-squared Error (MSE) benefit in using ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ rather than ${\widehat{Y}}_{\text{reg}\mathrm{,2}},$ without which there would be little motivation for ${\widehat{Y}}_{\text{dec}}$ . Indeed, the large-sample Theorems say that the main large-sample variance term is always optimal for ${\widehat{Y}}_{\text{reg}\mathrm{,2}}$ (whether because it is the same as for ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ under the null hypothesis or because it is strictly better under model (2.7) with distinct slopes). However, we indicated following Theorem 3, in the bound (2.9), that ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ can have smaller second-order MSE than ${\widehat{Y}}_{\text{reg}\mathrm{,2}},$ and the $H0$ columns of Tables 4.2 and 4.3 do show a small but consistent SD advantage for ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ versus ${\widehat{Y}}_{\text{reg}\mathrm{,2}},$ an advantage which is more pronounced in $M2.$ This advantage disappears under the fixed alternative $M1.H1$ but interestingly, not under $M2.H1.$ The slight but real conditional MSE advantage for ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$when the substratum slopes are very close to equality is discussed further by Slud (2012).

The estimators ${\widehat{Y}}_{\text{reg}\mathrm{,1}}\mathrm{,}{\widehat{Y}}_{\text{reg}\mathrm{,2}}\mathrm{,}{\widehat{Y}}_{\text{dec}}$ considered here are of regression type, and it may be of interest to compare their MSE behavior in the simulated populations with that of the simpler Horvitz-Thompson estimator ${\widehat{Y}}_{\text{HT}}$ in (1.1). All of these estimators are nearly unbiased, so that MSE's are essentially the same as variances, and a comparison of the third and fifth columns of Table 4.4 shows that the ${\widehat{Y}}_{\text{HT}}$ variances are considerably larger than those of ${\widehat{Y}}_{\text{dec}}.$ The difference is least pronounced with the larger sample sizes, but even there is 30-55%. The advantage of ${\widehat{Y}}_{\text{dec}}$ is still very pronounced in model $M2,$ where model variances and distributional skewness are larger, but less so than in model $M1.$

(IV) The definition of ${\widehat{Y}}_{\text{dec}}$ contains the arbitrary nominal significance level $\tau ,$ which in all tables other than Table 4.4 was taken to be 0.05. As the large-sample theory suggests, the properties of the decision-based estimator fall between those of ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ and ${\widehat{Y}}_{\text{reg}\mathrm{,2}},$ and larger values of $\tau$ make ${\widehat{Y}}_{\text{dec}}$ more often equal to ${\widehat{Y}}_{\text{reg}\mathrm{,1}}.$ As can be seen from comparison of columns 6 and 7 of Table 4.4, the choice $\tau =$ 0.20 seems in the simulated models to lead to very slightly smaller SD of ${\widehat{Y}}_{\text{dec}}$ under model $M1,$ but in model $M2$ the SD is if anything larger at the smaller sample sizes. The conclusion is weak because the differences are quite small compared to the differences between SD's from one frame population to another. Our preference is to let smaller $\tau$ dictate the frequent pooling of substrata except when there are pronounced differences in estimated slope between the substrata. This finding that larger significance levels $\tau$ do not improve performance of ${\widehat{Y}}_{\text{dec}}$differs from the finding of Saleh (2006) that larger significance levels are highly beneficial in other preliminary-testing contexts.

(V) Table 4.6 gives information about the variability across frame populations of SD estimators for the $Y$ estimators. The bootstrap variance estimators appear less susceptible to variation across frame populations, because the bootstrap averaging stabilizes them. The key finding in this table seems to be that the variability across frame populations is moderate except in the unweighted $M2$ setting, where it is remarkably large. This seems to account for the extreme inflation of variances under $M2.U$ seen in Table 4.5.

(VI) In many bootstrap applications with approximately normally distributed statistics, failure of coverage of normal-theory-based confidence intervals due to nonnormality of the bootstrapped statistic can be mitigated by using the bootstrap percentile (BP) intervals (Shao and Tu 1995, Section 4.1). In the present simulations, Table 4.4 (columns 4 and 6) gives the coverage percentages of BP intervals for ${\widehat{Y}}_{\text{dec}}$ in settings where Tables 4.2 and 4.3 give the coverages of normal-theory CI's based on the bootstrap-estimated SD. For whatever reason, the tables show that the normal-theory coverage ${\text{CP}}_B$ tends systematically to be slightly below nominal and yet slightly larger than the BP interval coverage ${\text{CP}}_{BP}.$ Thus, our simulations indicate the preference in this setting for the simpler interval ${\widehat{Y}}_{\text{dec}}\pm 1.96\cdot {\widehat{\text{SD}}}_B.$

(VII) It remains to draw lessons from the simulations with real government-census data in Section 4.3. The first necessary comment is that the spread and skewness of the full-time payroll predictors $x_i$ and the total-payroll size-variable $z_i$ are very large, much more like the lognormal models $M2$ than the gamma models $M1.$ Table 4.8 indicates (in column 5) a consistent MSE advantage for ${\widehat{Y}}_{\text{reg}\mathrm{,2}}$ over ${\widehat{Y}}_{\text{reg}\mathrm{,1}}$ except in the CA Special-district case, although the difference is small in the CO Special-district and the FL, GA and PA Subcounty cases. It is notable in almost all of these examples that the bootstrap SD estimator for ${\widehat{Y}}_{\text{dec}}$ is more accurate than the substitution-formula estimator, despite the rather small numbers of sampled and unsampled non-SR units and (in several cases, as shown in Table 4.7) relatively large numbers of SR units. The substitution SD estimates are consistently too small while the bootstrap estimates are usually slightly high (i.e., generally ${\widehat{\text{SD}}}_S\text{<}\widehat{\text{SD}}\text{<}{\widehat{\text{SD}}}_B$ ). The relative error of ${\widehat{\text{SD}}}_B$ versus $\widehat{\text{SD}}$ is no more than about 5% in these examples, except in the cases (AL, CO, GA) where there are particularly few non-sampled non-SR units in the large-unit substratum.

The large-unit substrata in ASPEP usually have small total frame population and often have relatively large numbers of SR units. While we have seen in these simulations that this does not quite invalidate inferences drawn with ${\widehat{Y}}_{\text{reg}\mathrm{,1}}\mathrm{,}{\widehat{Y}}_{\text{reg}\mathrm{,2}}$ or ${\widehat{Y}}_{\text{dec}},$these statistics have distributions rather different from those of large-sample theory, and perhaps future substratum splits should allow slightly larger large-unit substrata for well-behaved statistical inferences.

More broadly, the simulation results indicate that the decision-based estimator with interval estimator defined from bootstrap variances is well-behaved and can be recommended except in extremely dispersed and skewed populations or in populations with large-unit sample sizes less than 20-25.

Previous | Next