Semiparametric quantile regression imputation for a complex survey with application to the Conservation Effects Assessment Project
Section 5. Simulations

• Previous
• Next

We construct a simulation study to represent properties of the CEAP data and design. An extended set of simulations using the simulation models of Chen and Yu (2016) yields similar results and is not presented here for brevity. The objectives of the simulations are to evaluate the variance estimator and to compare QRI to nonparametric and fully parametric alternatives.

The fully parametric imputation procedure is parametric fractional imputation (Kim, 2011). The imputation model specified for parametric fractional imputation (PFI) is $y_i\mathrm{<mo>=</mo>}{\gamma }_0+{\gamma }_1{x}_i+{ϵ}_i,$ where ${ϵ}_i\sim N\left(\mathrm{0,}{\sigma }_{ϵ}^2\right).$ The imputed values for PFI are generated as, $y_{ij}^{*}\sim N\left({\widehat{\gamma }}_0+{\widehat{\gamma }}_1{x}_i\mathrm{,}{\widehat{\sigma }}_{ϵ}^2\right),$ where $\widehat{\gamma }={\left({\widehat{\gamma }}_0\mathrm{,}{\widehat{\gamma }}_1\mathrm{,}{\widehat{\sigma }}_{ϵ}^2\right)}^{\prime }$ satisfies $S_w\left(\widehat{\gamma }\right)\mathrm{<mo>=</mo>}0,$

$$S_w\left(\gamma \right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^n}{\pi }_i^{-1}{\delta }_i{d}_i\mathrm{,}(5.1)$$

and $d_i={\left(y_i-{\gamma }_0-{\gamma }_1{x}_i\mathrm{,}\left(y_i-{\gamma }_0-{\gamma }_1{x}_i\right)x_i\mathrm{,}{\left(y_i-{\gamma }_0-{\gamma }_1{x}_i\right)}^2/{\sigma }_{ϵ}^2-1\right)}^{\prime }\text{}.$ By incorporating ${\pi }_i^{-1}$ in the score function (5.1), the estimator is consistent if the population model is a linear model with $iid$ normally distributed errors and either the MAR assumption in (2.7) or (2.6) holds.

The non-parametric imputation (NPI) procedure is based on Wang and Chen (2009). For NPI, the $j^{\text{th}}$ imputed value for nonrespondent $i,$ $y_{ij}^{*},$ is generated from a multinomial distribution with sample space $\left\{y_s\text{}\mathrm{<mo>:</mo>}{I}_s\mathrm{<mo>=</mo>}{\delta }_s\mathrm{<mo>=</mo>\; 1}\right\}.$ Specifically,

$$P\left(y_{ij}^{*}\mathrm{<mo>=</mo>}{y}_s\right)\mathrm{<mo>=</mo>}\frac{{\pi }_i^{-1}K\left\{\mathrm{<mo\; stretchy="false">(</mo>}{x}_i-x_s\mathrm{<mo\; stretchy="false">)</mo>}/h\right\}}{{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^N{I}_j{\delta }_j{\pi }_j^{-1}}K\left\{\mathrm{<mo\; stretchy="false">(</mo>}{x}_i-x_j\mathrm{<mo\; stretchy="false">)</mo>}/h\right\}}\mathrm{,}(5.2)$$

where $K(\cdot )$ is a normal kernel with bandwidth $h$ selected by applying the method of Sheather and Jones (1991), as implemented in the R function $dpik,$ to $\left\{x_i\text{}\mathrm{<mo>:</mo>}{I}_s\mathrm{<mo>=</mo>}{\delta }_s\mathrm{<mo>=</mo>\; 1}\right\}.$

The QRI procedure is implemented as described in Sections 2-3. To define the penalized B-spline, we set $p\mathrm{<mo>=</mo>\; 3},$ $m\mathrm{<mo>=</mo>\; 2},$ $K_n\mathrm{<mo>=</mo>\; 16},$ and $\lambda \mathrm{<mo>=</mo>}\text{0}\text{.004}.$ The value of $\lambda \mathrm{<mo>=</mo>}\text{0}\text{.004}$ is the median of the values selected using the R function “cobbs” across 1,000 samples of a preliminary simulation. To select $\lambda$ using “cobbs”, we first use the R function “cobbs” to obtain ${\lambda }_{{\tau }_j}$ for ${\tau }_1\mathrm{,}\dots \mathrm{,}{\tau }_J.$ The selected $\lambda$ is the minimum of the $\left\{{\lambda }_{{\tau }_j}\text{}\mathrm{<mo>:</mo>}j\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}J\right\},$ which introduces the least amount of smoothing from among the selected ${\lambda }_{{\tau }_j}\text{}.$

In simulations not presented here, we also consider multiple imputation. Modifications to standard multiple imputation procedures are needed to produce unbiased estimators for a situation in the sample missing at random assumption (2.6) does not hold (Berg et al., 2016; Reiter, Raghunathan and Kinney, 2006). Because an exploration of the modifications to multiple imputation needed to ensure consistent estimation is beyond the scope of this study, we restrict attention to PFI, NPI, and QRI.

For all three imputation procedures, GMM based on the imputed values is used to estimate the parameters. Note that this differs from Wang and Chen (2009), which uses empirical likelihood instead of GMM. The number of imputations for the simulation is $J\mathrm{<mo>=</mo>\; 50.}$ The Monte Carlo (MC) sample size is 1,000.

We consider estimation of several parameters: ${\theta }_1\mathrm{<mo>=</mo>}E\left[y_i\right],$ ${\theta }_2\mathrm{<mo>=</mo>}V\left\{y_i\right\},$ ${\theta }_3\mathrm{<mo>=</mo>}\text{Cor}\left\{y_i\mathrm{,}{x}_i\right\},$ ${\theta }_4\mathrm{<mo>=</mo>}E\left[E\left[y_i|x_i\le \text{0}\text{.65}\right]\right],$ and ${\theta }_5\mathrm{<mo>=</mo>}P\left(y_i\le 8\right).$ With the exception of ${\theta }_5,$ GMM estimators of these parameters satisfy the assumptions required for the theory of Section 3. In particular, the function $g_i\left(\cdot \mathrm{<mo>;</mo>}\theta \right)$ defining the estimator of $\left({\theta }_1\mathrm{,}{\theta }_2\mathrm{,}{\theta }_3\mathrm{,}{\theta }_4\right)$ has two continuous derivatives. The estimator of ${\theta }_5$ does not fall in the framework of Section 3 because $I\left[a\le 8\right]$ is a non-smooth function of $a;$ however, we evaluate the empirical properties of ${\widehat{\theta }}_5$ defined as

$${\widehat{\theta }}_5={\left({\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^n}{\pi }_i^{-1}\right)}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^n}{\pi }_i^{-1}\left\{{\delta }_{i}I\left[y_i\le 8\right]+\left(1-{\delta }_i\right)J^{-1}{\displaystyle \sum _{j\mathrm{<mo>=</mo>\; 1}}^J}I\left[y_{ij}^{*}\le 8\right]\right\}\mathrm{.}(5.3)$$

For details on the function $g_i\left(\cdot ;\theta \right)$ defining the estimators for the simulation, see Section D of the online supplement https://github.com/emilyjb/Semiparametric-QRI-Supplement/blob/master/SupplementToQRI.pdf, (Berg and Yu, 2016).

5.1  Superpopulation model and design for simulations

The superpopulation model represents four aspects of the CEAP data and survey: (1) the shape of the expectation function, (2) the inclusion of a mean-variance relationship, (3) the use of probability proportional to size (PPS) with-replacement sampling, and (4) the sample sizes and response rates. The specific model for the simulation is $y_i\mathrm{<mo>=</mo>}m\left(x_i\right)+e_i,$ where $e_i\sim N\left(\mathrm{0,}{\sigma }_e^{2}m{\left(x_i\right)}^2\right),$ $m\left(x_i\right)\mathrm{<mo>=</mo>\; 2}+10{\left(1+8\text{exp}\left(-5x_i\right)\right)}^{-\frac{5}{4}}\text{},$ and $x_i\sim \text{Trunc}\text{.}\text{Norm}\text{.}\left(\text{0}\text{.5}\mathrm{,}\text{0}\text{.3}\right).$ The sample design is PPS with replacement, where the probability of selecting unit $i$ on a single draw is ${\left({\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^N}{\tilde{\psi }}_i\right)}^{-1}{\tilde{\psi }}_i\text{},$ $\text{logit}\left({\tilde{\psi }}_i\right)\mathrm{<mo>=</mo>}-3-\text{0}\text{.33}{z}_i+\text{0}\text{.1}{y}_i\text{},$ $z_i\sim \text{Trunc}\text{.}\text{Norm}\text{.}\left(\text{0}\text{.5}\mathrm{,}\text{0}\text{.3}\right),$ and $N\mathrm{<mo>=</mo>}\text{50,000}.$ The number of draws is $n\mathrm{<mo>=</mo>}\text{1,500},$ leading to a median sample size of 1,477, where the sample size is the number of unique units in the sample. The first and second order selection probabilities corresponding to ${\tilde{\psi }}_i$ are, ${\pi }_i\mathrm{<mo>=</mo>\; 1}-{\left(1-{\tilde{\psi }}_i\right)}^n\text{},$ and ${\pi }_{ij}\mathrm{<mo>=</mo>\; 1}-{\left(1-{\tilde{\psi }}_i\right)}^n-{\left(1-{\tilde{\psi }}_j\right)}^n+{\left(1-{\tilde{\psi }}_j-{\tilde{\psi }}_i\right)}^n\text{}\mathrm{.}$ The response indicator ${\delta }_i\sim \text{Bernoulli}\left(p_i\right),$ where $\text{logit}\left(p_i\right)\mathrm{<mo>=</mo>}\text{0}\text{.5}{x}_i+\text{1}\text{.5}{z}_i,$ which yields a median response rate of 0.631.

By the model for $y_i$ given $x_i,$ the assumption of population missing at random (2.7) holds for this simulation. Incorporating $z_i$ in the models for $p_i$ and ${\pi }_i$ is the approach used in Berg et al. (2016) that causes the sample missing at random assumption (2.6) to fail. The variable $z_i$ can be interpreted a design variable that is omitted from the imputation model.

5.2  Results

Table 5.1 contains three measures for comparing the QRI estimator to the PFI and NPI estimators. The percent relative MC MSE for estimator $k\left(k=\text{PFI},\text{NPI}\right)$ is defined,

$$\text{Pct}\text{.}\text{Rel}\text{.}\text{MSE}\left(k\right)\mathrm{<mo>=</mo>\; 100}\frac{{\text{MSE}}_{\text{MC}}\left(\widehat{\theta }\left(k\right)\right)-{\text{MSE}}_{\text{MC}}\left(\widehat{\theta }\left(\text{QRI}\right)\right)}{{\text{MSE}}_{\text{MC}}\left(\widehat{\theta }\left(\text{QRI}\right)\right)}\mathrm{,}(5.4)$$

where $\widehat{\theta }\left(k\right)$ is the estimator based on imputation procedure $k.$ The percent relative variance for estimator $k$ is defined,

$$\text{Pct}\text{.}\text{Rel}\text{.}\text{Var}\left(k\right)\mathrm{<mo>=</mo>\; 100}\frac{{\text{Var}}_{\text{MC}}\left(\widehat{\theta }\left(k\right)\right)-{\text{Var}}_{\text{MC}}\left(\widehat{\theta }\left(\text{QRI}\right)\right)}{{\text{Var}}_{\text{MC}}\left(\widehat{\theta }\left(\text{QRI}\right)\right)}\mathrm{,}(5.5)$$

for $k\mathrm{<mo>=</mo>}\text{NPI},\text{PFI}\text{.}$ The percent of mean squared error due to squared bias is defined by

$$\text{Pct}\text{.}\text{Bias}\left(k\right)\mathrm{<mo>=</mo>\; 100}\frac{{\left(E_{\text{MC}}\left(\widehat{\theta }\left(k\right)\right)-\theta \right)}^2}{{\text{MSE}}_{\text{MC}}\left(\widehat{\theta }\left(k\right)\right)}\mathrm{,}(5.6)$$

where $k\mathrm{<mo>=</mo>}\text{NPI},\text{PFI,}\text{QRI}\text{.}$ The MSE of the QRI estimator is smaller than the MSE of the NPI and PFI estimators for all parameters. The PFI estimator is biased because the model underlying the PFI procedure does not account for the nonlinearity in the quantile curves or the nonconstant variances. The NPI procedure has a relatively large variance for sample sizes such as those obtained in the CEAP survey. The squared MC bias of the QRI procedure is less than 0.5% of MC MSE for all parameters.

The last two columns of Table 5.1 contain the relative bias of the variance estimator and the empirical coverage of normal theory 95% confidence intervals. The relative bias of the variance estimator defined as

$$\text{Rel}\text{.}\text{Bias}\mathrm{<mo>=</mo>}\frac{E_{\text{MC}}\left[\widehat{V}\left(\widehat{\theta }\right)\right]-V_{\text{MC}}\left(\widehat{\theta }\right)}{V_{\text{MC}}\mathrm{,}\left(\widehat{\theta }\right)}\mathrm{,}(5.7)$$

where $E_{\text{MC}}\left[\widehat{V}\left(\widehat{\theta }\right)\right]$ is the MC mean of the variance estimators and $V_{\text{MC}}\left(\widehat{\theta }\right)$ is the MC variance of the QRI estimator. The MC relative bias of the variance estimator for the QRI estimator is between -6% and -1%. Empirical coverages of normal theory confidence intervals are within 1% of the nominal 95% level.

• Previous
• Article main page
• Next