Estimation and inference of domain means subject to qualitative constraints
Section 4. Performance of constrained estimator

• Previous
• Next

4.1  Simulations

We run simulation experiments to measure the performance of the proposed methodology to carry out estimation and inference of population domain means. Given a pair of natural numbers $D_1$ and $D_2,$ we generate the limiting domain means ${\mu }_d$ from the monotone bivariate function $\mu \left(x_1\mathrm{,}{x}_2\right)$ given by

$$\mu \left(x_1\mathrm{,}{x}_2\right)\mathrm{<mo>=</mo>}\sqrt{1+4x_1/D_1}+\frac{4\exp \left(0.5+2x_2/D_2\right)}{1+\exp \left(0.5+2x_2/D_2\right)}\mathrm{.}$$

The ${\mu }_d$ are created by evaluating $\mu \left(x_1\mathrm{,}{x}_2\right)$ at every combination of $x_1\mathrm{<mo>=</mo>}\mathrm{1,}\mathrm{2,}\dots \mathrm{,}{D}_1$ and $x_2\mathrm{<mo>=</mo>}\mathrm{1,}\mathrm{2,}\dots \mathrm{,}{D}_2,$ producing a total number of domains equal to $D\mathrm{<mo>=</mo>}{D}_1{D}_2.$ We set $D_1\mathrm{<mo>=</mo>}6$ and $D_2\mathrm{<mo>=</mo>}4.$ Note that although the function $\mu \left(x_1\mathrm{,}{x}_2\right)$ produces a matrix rather than a vector of domain means, it can be vectorized in order to represent the limiting domain means as the vector $\mu .$ For each domain $d,$ we generate its $N_d\mathrm{<mo>=</mo>}N/D\mathrm{<mo>=</mo>}$ 400 elements by adding independent and normally distributed noise with mean 0 and variance ${\sigma }^2$ to the ${\mu }_d.$ Once the elements of the population have been simulated, then the population domain means ${\overline{y}}_U$ are computed. The population domain means used for simulations when $\sigma \mathrm{<mo>=</mo>}1$ are displayed in Figure 4.1. Observe that these domain means are reasonably (not strictly) monotone with respect to $x_1$ and $x_2.$

Samples are drawn from a stratified sampling design without replacement, with 4 strata that cut across the $D$ domains. Strata are constructed using an auxiliary variable $\nu$ that is correlated with the variable of interest $y.$ The vector $\nu$ is created by adding independent standard normally distributed noise to $\sigma d/D,$ for each element in domain $d.$ Then, stratum membership is assigned by sorting the vector $\nu$ and creating 4 blocks of $N/4\mathrm{<mo>=</mo>}$ 2,400 elements each based on the sorted $\nu .$ To make the design informative, we sample $n\mathrm{<mo>=</mo>}$ 480 elements divided across strata in (60, 120, 120, 180). This probability sampling design is similar to the one described in Wu et al. (2016).

We consider 4 different scenarios obtained from the combination of two possible types of shape constraints and $\sigma \mathrm{<mo>=</mo>}1$ or 2. The first type of constraints assumes the population domain means are monotone increasing with respect to both $x_1$ and $x_2$ (double monotone), while the second type of constraints assumes monotonicity only with respect to $x_1$ (only $x_1$ monotone). For a fixed $\sigma ,$ the exact same population is considered for the two possible types of constraints. For each scenario, the unconstrained ${\tilde{y}}_s$ and constrained ${\tilde{\theta }}_s$ estimates are computed along with their linearization-based variance estimates (see (2.11)). Constrained estimates are computed using the CPA, and their variance estimates are computed by relying on the sample-selected set $J\in G_s.$ In addition, 95% Wald confidence intervals based on the normal distribution are constructed for both estimators.

To measure the precision of ${\tilde{y}}_s$ and ${\tilde{\theta }}_s$ as estimators of the population domain means ${\overline{y}}_U,$ we consider the Weighted Mean Squared Error (WMSE) given by

$$\text{WMSE}\left({\tilde{\phi }}_s\right)\mathrm{<mo>=</mo>}E\left[{\left({\tilde{\phi }}_s-{\overline{y}}_U\right)}^T{W}_U\left({\tilde{\phi }}_s-{\overline{y}}_U\right)\right]\mathrm{,}$$

where ${\tilde{\phi }}_s$ could be either the unconstrained or constrained estimator and $W_U$ is the diagonal matrix with elements $N_d/N,$ $d\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}D.$ The WMSE values are approximated by simulations as

$$\frac{1}{B}{\displaystyle \sum _{b\mathrm{<mo>=</mo>\; 1}}^B}{\left({\tilde{\phi }}_s^{\left(b\right)}-{\overline{y}}_U\right)}^T{W}_U\left({\tilde{\phi }}_s^{\left(b\right)}-{\overline{y}}_U\right)\mathrm{,}$$

where $B$ is the number of simulations, and ${\tilde{\phi }}_s^{\left(b\right)}$ is the estimator for the $b^{\text{th}}$ sample.

Simulation results are summarized in Figures 4.2 - 4.5, and are based on $R\mathrm{<mo>=</mo>}$ 10,000 replications. These display the 24 domains divided in groups of 6, where each group is assumed to be monotone. For the double monotone scenario, similar plots with groups of 4 monotone domains each can be also pictured. As illustrated in the fits of a single sample in these figures, it can be seen that the constrained estimates can be exactly equal to the unconstrained estimates for some domains. In those cases, their variance estimates are also equal. Overall, confidence intervals for the constrained estimator tend to be tighter in comparison with those for the unconstrained estimator. On average, the constrained estimator behaves slightly differently than the population domain means, due to the latter’s non-strict monotonicity. As an advantage, the percentiles for the constrained estimator are narrower, demonstrating that the distribution of the proposed estimator is tighter than the distribution of the unconstrained estimator. For small values of $\sigma ,$ the unconstrained estimates are more likely to satisfy the assumed restrictions, which leads to small improvements on the constrained estimator over the unconstrained. In contrast, shape assumptions tend to be more severely violated in unconstrained estimates for larger values of $\sigma ,$ allowing the proposed estimator to gain much more efficiency on these cases. This latter property can be noted by observing that the constrained estimator percentile band gets farther away from the unconstrained estimator band as $\sigma$ increases.

In terms of variability, the constrained estimator has the smaller variance of the two estimators. Interestingly, it gets overestimated by its corresponding linearization-based variance estimate. In contrast, the variance estimate of the unconstrained estimator underestimates the true variance, which is a known and often observed drawback of linearization variances. Despite this difference, confidence intervals for both estimators demonstrate a similar good coverage rate when $\sigma \mathrm{<mo>=</mo>}1,$ meanwhile such coverage gets slightly improved by the constrained estimator when $\sigma \mathrm{<mo>=</mo>}2.$

Table 4.1 shows that the constrained estimator is more precise on average than the unconstrained estimator. The precision of the constrained estimator improves when the monotonicity with respect to the two variables is assumed, instead of only with respect to $x_1.$ This is expected here, because the underlying surface is indeed doubly monotone, so that the estimator benefits from imposing the stronger constraint.

4.2  Replication methods for variance estimation

In practice, it is common for large-scale surveys to use replication-based methods for variance estimation. Examples of such surveys are the last editions of the NHANES and the National Survey of College Graduates (NSCG). To study the performance of replication-based variance estimators under the proposed constrained methodology, we carry out simulation studies based on the delete-a-group Jackknife (DAGJK) variance estimator proposed by Kott (2001).

We perform replication-based simulation experiments using the setting described in Section 4.1. To compute the DAGJK variance estimator, we first randomly create $G$ equal-sized groups within each of the 4 strata. Then, for each replicate $g\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}G,$ we delete the $g^{\text{th}}$ group in each of the strata, adjust the remaining weights by $w_k^{\left(g\right)}\mathrm{<mo>=</mo>}\left(\frac{G}{G-1}\right)w_k,$ where $w_k\mathrm{<mo>=</mo>}{\pi }_k^{-1};$ and compute the replicate constrained estimate ${\tilde{\theta }}_s^{\left(g\right)}$ using the adjusted weights. The DAGJK variance estimate of ${\tilde{\theta }}_{s_d},$ ${\widehat{V}}_{\text{JK}}\left({\tilde{\theta }}_{s_d}\right),$ is obtained by calculating

$${\widehat{V}}_{\text{JK}}\left({\tilde{\theta }}_{s_d}\right)\mathrm{<mo>=</mo>}\frac{G-1}{G}{\displaystyle \sum _{g\mathrm{<mo>=</mo>\; 1}}^G}{\left({\tilde{\theta }}_{s_d}^{\left(g\right)}-{\tilde{\theta }}_{s_d}\right)}^2\text{}\mathrm{.}$$

A replication-based variance estimator of ${\tilde{y}}_{s_d}$ is obtained by substituting ${\tilde{\theta }}_s$ by ${\tilde{y}}_s.$

Our simulations consider only the double monotone scenario, with $\sigma \mathrm{<mo>=</mo>}$ 1 or 2, and $G\mathrm{<mo>=</mo>}$ 10, 20 or 30. The sample size is set to either $n\mathrm{<mo>=</mo>}$ 480 or $n\mathrm{<mo>=</mo>}$ 960, where the latter case is obtained by doubling the original sample size in each strata. Figures 4.6 - 4.9 contain simulation results based on 10,000 replications. In contrast to the behavior of the linearization-based variance estimates, it can be seen that the DAGJK estimates tend to overestimate the variance of the unconstrained estimator, as is often observed in practice. Both replication-based and linearization-based variance estimates of the constrained estimator overestimate the true variance, so that the results are more consistent across variance estimation methods. As the number of groups $G$ increases, DAGJK estimates tend to be greater, especially for small values of $\sigma .$ Such increments on DAGJK estimates have the direct consequence of increasing the coverage rate as $G$ gets larger. In addition, the coverage rate for both estimators is improved (closer to 0.95) when the sample size is increased. Overall, it appears that replication variance estimation is a practical alternative to linearization.

• Previous
• Article main page
• Next