Robust variance estimators for generalized regression estimators in cluster samples
Section 3. Simulation

• Previous
• Next

We performed a series of simulation studies to test the performance of the new variance estimators in different populations. In each simulated sample, we computed the quantities listed in Table 3.1. To evaluate the variance estimators, we calculated the average of the variance estimates, compared those averages to the empirical mean square error, and computed coverage probabilities of confidence intervals based on the different variance estimates. Table 3.2 summarizes the sample designs for the 18 simulation studies. The column called Label gives the headings used in later tables. The sample designs are used in three populations described below.

3.1  Data

We conducted simulations on three populations to assess the design-based performance of the variance estimators under a variety of situations. In the first population, we investigated the performance of the variance estimators when the first-stage sampling fraction was large and the sample size was moderate. The focus of the second simulation study was on the performance of the variance estimators under a relatively messy dataset and a small first-stage sample size. The final simulation study shows the performance of the variance estimators in large samples.

3.1.1  Third grade population

The first simulation study used the Third Grade population from Appendix B.6 of Valliant et al. (2000). This dataset contained the mathematics achievement scores for 2,427 third graders in 135 schools. The relatively small number of schools in this population and the fairly constant number of students in each school made it ideal for studying samples with large sampling fractions.

We used GREG to estimate the average mathematics achievement score for third graders. Altogether, we selected 1,000 samples in each of six sample designs listed in Table 3.2. In the first sample design, we selected 1,000 simple random samples without replacement (srswor) of 25 schools. Within each sampled school, we selected exactly five students via srswor. Because the number of students in each school varied from school to school, this sample design resulted in different unconditional probabilities of selection, but a fixed sample size of 125 students. The second sample design was similar to the first, except we selected 50 schools. Selecting 50 of the 135 schools resulted in a large first-stage sampling fraction of 0.37, necessitating a finite population correction factor. Both the samples of $m\mathrm{<mo>=</mo>}$25 and 50 might be considered to be of “moderate” size.

In the third sample design, we selected 1,000 simple random samples of 25 schools without replacement. Within each sampled school, we selected students at a constant rate of $\frac{675}{\text{2,427}},$ yielding 1,000 samples with random sizes centered around 125 students. The result of this design was that each student had the same unconditional probability of selection. The fourth sample design was similar to the third, except we selected 50 schools. The sample sizes were also random under this design, with an average of 250 students. Since the third and fourth sample designs resulted in every unit getting the same chance of selection, these sample designs are labeled srs epsem (equal probability selection mechanism) in subsequent tables.

In the fifth design, we selected 1,000 samples of 25 schools with probabilities proportional to the number of students in each school. Within each sampled school, we selected exactly five students, yielding 1,000 samples with exactly 125 students each. The sixth sample design was similar to the fifth, except we selected 50 schools. We selected 1,000 samples of size 250 students using this design. The fifth and sixth designs are epsem. Like the second and fourth sample designs, this sample design also had a large sampling fraction and warranted the need for a finite population correction factor to adjust the variance estimators.

From each sample, we estimated the average achievement scores for the finite population using a GREG estimator and assuming that the number of students in the population was known. The assisting model was meant to replicate the clustered linear regression model in Section 9.6 of Valliant et al. (2000). The eleven explanatory variables used to model each student’s math achievement score were: an intercept, sex (male or female), ethnicity (White/Asian, Black, Native American/Other, or Hispanic), language spoken at home is the same as the test (Always, Sometimes/Never), type of community (Outskirts of a town or city, Village/City), and school enrollment. The total mathematics achievement estimated with the GREG estimator was divided by the number of students in the population, 2,427, to get the average achievement score. The average achievement score for the population was 477.7. For the full population, the R-squared for the student-level linear model was 0.9735, indicating a very strong linear relationship.

3.1.2  American Community Survey population

The second simulation study used Census 2000 Summary File 3 data and American Community Survey (ACS) 2005 - 2009 Summary File data. The goal was to estimate the total number of housing units in the U.S. state of Alabama as reported in the ACS Summary File. Block group counts from Census 2000 were used as covariates in the assisting model.

To create the population, first all block group data were extracted from the ACS Summary File and the Census 2000 Summary File 3. Then, the two files were merged at the block group level. Block groups with 1,000 or more housing units in Census 2000 were removed because such large block groups had different characteristics than the majority of blocks. In many sampling designs such large units would be placed in a separate, certainty stratum and not contribute to the variance of estimates. Also, block groups with extreme growth in the total number of housing units were also removed. Specifically block groups that had gained more than 10 units over twice the 2000 census count were removed.

Clusters were defined as counties and block groups were treated as units. Treating block group as a unit is motivated by the common task of selecting a sample of blocks, listing them, and then using the listings to estimate the total number of housing units in the finite population.

Clusters with fewer than 10 block groups or more than 120 block groups in them were removed from the frame of clusters. Overall, there were 61 clusters (counties) containing a total of 2,051 block groups and 1,109,499 housing units in the edited dataset. Altogether, six counties and 1,278 block groups containing 1,030,471 housing units were removed from the Alabama file.

Figure 3.1 shows two scatterplots. The first plot shows the total number of housing units in the block group as reported on the ACS summary file as a function of the 2000 Census housing unit count. Each point represents one of the 2,051 block groups in the finite population. The diagonal line is a nonparametric smoother, indicating a strong relationship between the two variables. The plot also shows some evidence of heteroscedasticity because the points appear to fan out as the 2000 census count increases. The second plot shows the residuals obtained by regressing the 2000 census housing unit count on the ACS housing unit count using ordinary least squares (OLS) plotted versus the ACS housing unit count. As the number of housing units reported on the ACS file increases, the model predictions appear to seriously underestimate the true number of housing units. This suggests some degree of nonlinearity in the mean function. In addition, there is noticeable heteroscedasticity in variance.

As in the first simulation study, we tested six different sample designs. We selected 5,000 samples in each of six different selection mechanisms listed in Table 3.2. In the first sample design, we selected 5,000 simple random samples of 3 clusters without replacement. In large national surveys, it is not uncommon to select a small number of primary sampling units in each stratum. In this case, we treat Alabama as if it were a single design stratum and its 61 counties as clusters. Three counties within that stratum were sampled. Within each cluster, we selected nine block groups using srswor. The second design was similar with 15 clusters and 9 block groups per cluster. The first two sample designs resulted in highly variable weights. The other designs (rows 9-12) were parallel to those in rows 3-6 for the Third Grade population. The sample sizes of $m=$ 3 and 15 are small so that theoretical, large sample properties are less likely to hold.

From each sample, we estimated the total number of housing units in the finite population using a GREG estimator. The assisting model included an intercept and the Census 2000 count of housing units; the heteroscedasticity noted above was not accounted for in the GREG. For the full population, the R-squared was 0.819, again indicating a strong linear relationship.

3.1.3  Simulated population

A population was created with a large number of clusters to assess the asymptotic characteristics of the variance estimators. Generated using a classic linear model, a total of 30,000 clusters were created, each with a random number of units. The number of units in each cluster was determined by adding three to a uniform random integer between 0 and 7. This created clusters ranging in size from 3 to 10 units. Altogether, the population contained 195,164 units within 30,000 clusters. For each unit, a positive covariate was created as $x_k\sim \text{1,000}\exp N\left(\mathrm{0,}1\right)$ where $N\left(\mathrm{0,}1\right)$ is a normal random variate with mean of 0 and standard deviation of 1. A random response was created such that $y_k\sim N\left(\text{1,000}+2x_k\mathrm{,}\frac{x_k}{2}\right)$ . Figure 3.2 shows scatter plots of the relationship between $x_k$ and $y_k$ for the finite population.

We selected samples using the six different probability selection mechanisms listed in rows 13-18 of Table 3.2. The types of sample designs are parallel to those used for the Third Grade and ACS populations. In designs 14, 16, and 18, we selected 100 simple random samples of 1,500 clusters without replacement. We only selected 100 samples due to the excessive amount of computer time it took to select and process each sample. The sample sizes of $m=$ 300 and 1,500 are large so that theoretical, large sample properties should hold.

From each sample, we estimated the total of the response using a GREG estimator. The true finite population total was 839,149,969. The assisting model included an intercept and $x$ with $Q=I.$ For the full population, the R-squared was 0.953, indicating a very strong linear relationship. Figure 3.2 shows a scatter plot of the population as well as a residual plot based on an OLS regression of $x_k$ on $y_k$ for the full population. There is clear evidence of heteroscedasticity of errors.

3.2  Results

We explored the bias, variability, and confidence interval coverage of the new and existing variance estimators. We only show tables for some of the simulations to conserve space. Table 3.3 shows the means of the $\pi$ -estimator and the GREG estimator as well as the ratios of the average values of the variance estimators to the empirical mse’s for all populations and sample size combinations across all simulations. Both the $\pi$ -estimator and the GREG estimator are approximately unbiased; however, the GREG estimator is much more efficient.

The performance of the variance estimators depends on the sample design and the population. Some of the estimates in Table 3.3 from the ACS population with the simple random sample of 3 clusters and 9 units in each cluster stand out as being extremely poor. The inverses of the probabilities of selection vary quite a bit for this sample design. The variability of these weights, coupled with some extreme observations in the population, causes instability for some of the variance estimators. Namely, ${\upsilon }_{J2},$ ${\upsilon }_{\text{Jack}},$ ${\upsilon }_{J1},$ ${\upsilon }_{J2}^{*},$ ${\upsilon }_{\text{Jack}}^{*},$ ${\upsilon }_{J1}^{*}$ are extreme overestimates on average. All six of these estimators contain explicit or implicit hat matrix adjustments which can be quite large and seriously inflate the variance estimators when coupled with large sampling weights. On the other hand, ${\upsilon }_D,$ which also has a hat matrix adjustment, performs reasonably well for all populations and sample sizes. Noteworthy is the result that ${\upsilon }_D$ is much less of an overestimate for the mse in the combination (ACS, srs fixed, $m\mathrm{<mo>=</mo>}$3, $n_i\mathrm{<mo>=</mo>}$9) whereas other hat-matrix adjusted estimators were extreme overestimates. The estimators, ${\upsilon }_g,$ ${\upsilon }_{wr},$ and to a lesser extent, ${\upsilon }_R$ and ${\upsilon }_{JL},$ tend to be underestimates at the smaller sample sizes in the Third Grade and ACS populations and for all sample designs in those populations, but the problem diminishes for the larger sample sizes.

The boxplots in Figure 3.3 show the variability of the estimators more clearly for srs’s of size $m\mathrm{<mo>=</mo>}$25 and 50 from the Third Grade population. The boxplots depict the estimated standard errors (SEs) as a fraction of the empirical SE for the samples in each simulation. A ratio of 1 means that the estimated variance was equal to the empirical variance. Some samples yield large SE estimates, even though the majority of samples are much closer to the empirical variance. The degree of overestimation and the incidence of extreme values decreases substantially with the larger sample size as is evident by comparing the figures. The hat-matrix adjusted estimators also tend to somewhat overestimate the true variance, as evinced by the boxes that are shifted above the reference lines drawn at 1. This can be an advantage for confidence interval coverage.

Table 3.4 shows the six-number summaries of the ratios of the SE estimates, $\sqrt{v},$  to the square root of the empirical variance, $\sqrt{v_E},$  for the Third Grade population for four of the sample designs. As indicated by the median value of the ratios for ${\upsilon }_{J2},$ ${\upsilon }_{\text{Jack}},$ ${\upsilon }_{J1},$ ${\upsilon }_{J2}^{*},$ ${\upsilon }_{\text{Jack}}^{*},$  and ${\upsilon }_{J1}^{*},$  they are generally centered near the empirical SEs, but can have extremely large values in some samples that affect their averages. (The problem of outlying values is even more severe in the ACS population; details are not shown here.) The estimators that are least affected by extremes are ${\upsilon }_g,$ ${\upsilon }_{wr},$ ${\upsilon }_{JL},$ ${\upsilon }_R,$ ${\upsilon }_D,$ ${\upsilon }_R^{*},$  and ${\upsilon }_D^{*}.$  However, the estimators that incorporate fpc’s often are underestimates except in the case of srs and $m=$ 25.

Lastly, Table 3.5 shows the 95% confidence interval coverage for all of the estimators based on t-distributions. That is, we computed, $\left[{\widehat{t}}_y^{gr}-t_{\text{0}\text{.975}\mathrm{,}m-1}\sqrt{\upsilon }\mathrm{,}{\widehat{t}}_y^{gr}+t_{\text{0}\text{.975}\mathrm{,}m-1}\sqrt{\upsilon }\right]$ where $t_{\text{0}\text{.975}\mathrm{,}m-1}$ is the 97.5^th percentile from a t-distribution with $m-1$ degrees of freedom. We then noted how often the true value fell below, above, and inside this range. In addition to the new and old estimators, Table 3.5 also shows the confidence interval coverage attained when the empirical variance, ${\upsilon }_E,$ was used to form the confidence intervals. Ideally, the population total should be within the estimated 95% confidence interval for 95% of the samples. The true total should be below the 95% confidence bounds for 2.5% of the samples and above the confidence bounds for the same percentage of samples.

The jackknife-based estimators, ${\upsilon }_D^{*},$ ${\upsilon }_{\text{Jack}}^{*},$ and ${\upsilon }_{J2},$ cover at higher rates than the other variance estimators because they are larger. In small samples, jackknife-based estimators cover above the nominal level. The traditional variance estimators, ${\upsilon }_g,$ ${\upsilon }_{wr},$ and ${\upsilon }_{JL}$ under-covered in a number of cases, although their coverage was almost always higher than 90%. Note that ${\upsilon }_D$ is generally an improvement over ${\upsilon }_R$ due to the hat-matrix adjustment that makes ${\upsilon }_D$ larger.

The variance estimators that incorporate hat matrix adjustments $({\upsilon }_D,$ ${\upsilon }_{J2},$ ${\upsilon }_{\text{Jack}},$ and ${\upsilon }_R^{*})$ generally increase CI coverage rates compared to the other choices. This advantage was especially noticeable for the ACS population where, for example, ${\upsilon }_{wr}$ covers in less than 90% of samples in the combinations, $({\upsilon }_{\text{Jack}}^{*},$ $m\mathrm{<mo>=</mo>}$3), (srs epsem, $m\mathrm{<mo>=</mo>}$3), and (srs epsem, $m\mathrm{<mo>=</mo>}$15). Although, in principal, an fpc would seem useful in some of the population and sample size combinations, CIs based on the variance estimators with fpc’s cover at lower rates than their counterparts without the fpc’s. For example, in ACS (srs epsem, $m\mathrm{<mo>=</mo>}$15) the coverage rates for ${\upsilon }_R^{*},$ ${\upsilon }_D^{*},$ ${\upsilon }_{J2}^{*},$ ${\upsilon }_{\text{Jack}}^{*},$ and ${\upsilon }_{J1}^{*}$ range from 86.1 to 90.6% while the rates for the versions without fpc’s range from 90.2 to 93.4%.

One feature of ${\upsilon }_D$ and ${\upsilon }_D^{*}$ is that both the cluster-specific contributions, ${\upsilon }_{D\text{}\mathrm{,}i}$ and ${\upsilon }_{D,i}^{*},$ as well as the overall variance estimates can be negative. In the simulations, the adjustment described after (2.11) was used to avoid negative contributions. Negative estimates were more common when the second stage sample sizes were small and the weights were quite variable. For example, for the ACS population, almost 28% of the simple random samples of 3 clusters and $m_i\mathrm{<mo>=</mo>}$9 resulted in at least one negative variance contribution for a cluster. More commonly, about 10% of the samples contained at least one negative variance estimate for a cluster. In the Third Grade population, 16% to 27% of the samples had at least one negative value of ${\upsilon }_{Di}.$ In the simulated population with large sample sizes, ${\upsilon }_{Di}$ was negative in less than 5% of the samples. With the ad hoc correction of setting $I_i-H_{ii}$ to $I_i,$ ${\upsilon }_D$ is one of the most attractive variance estimators because it tends to slightly overestimate the empirical variance, has some of the best confidence interval coverage, and has reasonable variability compared to other variance estimators.

• Previous
• Article main page
• Next