Statistical matching using fractional imputation 6. Simulation studyStatistical matching using fractional imputation 6. Simulation study

• Previous
• Next

To test our theory, we present two limited simulation studies. The first simulation study considers the setup of combining two independent surveys of partial observation to obtain joint analysis. The second simulation study considers the setup of measurement error models with external calibration.

6.1  Simulation one

To compare the proposed methods with the existing methods, we generate 5,000 Monte Carlo samples of $\left(x_i\mathrm{,}{y}_{1i}\mathrm{,}{y}_{2i}\right)$ with size $n\mathrm{<mo>=</mo>\; 400},$ where

$$\left(\begin{array}{c}{y}_{1i}\\ x_i\end{array}\right)\sim N\left(\left[\begin{array}{l}2\hfill \\ 3\hfill \end{array}\right]\mathrm{,}\left[\begin{array}{cc}1& 0.7\\ 0.7& 1\end{array}\right]\right)\mathrm{,}$$

$$y_{2i}\mathrm{<mo>=</mo>}{\beta }_0+{\beta }_1{y}_{1i}+e_i\mathrm{,}(6.1)$$

$e_i\sim N\left(\mathrm{0,}{\sigma }^2\right),$ and $\beta \mathrm{<mo>=</mo>}({\beta }_0\mathrm{,}{\beta }_1\mathrm{,}{\sigma }^2{)}^{\prime }\mathrm{<mo>=</mo>}(\mathrm{1,1,1}{)}^{\prime }.$ Note that, in this setup, we have $f\left(y_2|x\mathrm{,}{y}_1\right)\mathrm{<mo>=</mo>}f\left(y_2|y_1\right)$ and so the variable $x$ plays the role of the instrumental variable for $y_1.$

Instead of observing $\left(x_i\mathrm{,}{y}_{1i}\mathrm{,}{y}_{2i}\right)$ jointly, we assume that only $\left(y_1\mathrm{,}x\right)$ are observed in Sample A and only $\left(y_2\mathrm{,}x\right)$ are observed in Sample B, where Sample A is obtained by taking the first $n_a\mathrm{<mo>=</mo>\; 400}$ elements and Sample B is obtained by taking the remaining $n_b\mathrm{<mo>=</mo>\; 400}$ elements from the original sample. We are interested in estimating four parameters: three regression parameters ${\beta }_0\mathrm{,}{\beta }_1\mathrm{,}{\sigma }^2$ and $\pi \mathrm{<mo>=</mo>}P\left(y_1\mathrm{<mo><</mo>\; 2,}{y}_2\mathrm{<mo><</mo>\; 3}\right),$ the proportion of $y_1\mathrm{<mo><</mo>\; 2}$ and $y_2\mathrm{<mo><</mo>\; 3.}$ Four methods are considered in estimating the parameters:

1. Full sample estimation (Full): Uses the complete observation of $\left(y_{1i}\mathrm{,}{y}_{2i}\right)$ in Sample B.
2. Stochastic regression imputation (SRI): Use the regression of $y_1$ on $x$ from Sample A to obtain $\left({\widehat{\alpha }}_0\mathrm{,}{\widehat{\alpha }}_1\mathrm{,}{\widehat{\sigma }}_1^2\right),$ where the regression model is $y_1\mathrm{<mo>=</mo>}{\alpha }_0+{\alpha }_{1}x+e_1$ with $e_1\sim \left(\mathrm{0,}{\sigma }_1^2\right).$ For each $i\in B,$ $m\mathrm{<mo>=</mo>\; 10}$ imputed values are generated by $y_{1i}^{\mathrm{<mo>*</mo>}\left(j\right)}\mathrm{<mo>=</mo>}{\widehat{\alpha }}_0+{\widehat{\alpha }}_1{x}_i+e_i^{\mathrm{<mo>*</mo>}\left(j\right)}$ where $e_i^{\mathrm{<mo>*</mo>}\left(j\right)}\sim N\left(\mathrm{0,}{\widehat{\sigma }}_1^2\right).$
3. Parametric fractional imputation (PFI) with $m\mathrm{<mo>=</mo>\; 10}$ using the instrumental variable assumption.
4. Hot-deck fractional imputation (HDFI) with $m\mathrm{<mo>=</mo>\; 10}$ using the instrumental variable assumption.

Table 6.1 presents Monte Carlo means and Monte Carlo variances of the point estimators of the four parameters of interest. SRI shows large biases for all parameters considered because it is based on the conditional independence assumption. Both PFI and HDFI provide nearly unbiased estimators for all parameters. Estimators from HDFI are slightly more efficient than those from PFI because the two-step procedure in HDFI uses the Full set of respondents in the first step. The theoretical asymptotic variance of ${\widehat{\beta }}_1$ computed from PFI is

$$V\left({\widehat{\beta }}_1\right)\doteq \frac{1}{{\left(0.7\right)}^2}\frac{1}{400}2\left(1-\frac{{0.7}^2}{2}\right)+\frac{1}{{\left(0.7\right)}^2}\frac{1}{400}\left(1-{0.7}^2\right)\doteq 0.0103$$

which is consistent with the simulation result in Table 6.1. In addition to point estimation, we also compute variance estimators for PFI and HDFI methods. Variance estimators show small relative biases (less than 5% in absolute values) for all parameters. Variance estimation results are not presented here for brevity.

The proposed method is based on the instrumental variable assumption. To study the sensitivity of the proposed fractional imputation method to violations of the instrumental variable assumption, we performed an additional simulation study. Now, instead of generating $y_{2i}$ from (6.1), we use

$$y_{2i}\mathrm{<mo>=</mo>\; 0.5}+y_{1i}+\rho \left(x_i-3\right)+e_i\mathrm{,}(6.2)$$

where $e_i\sim N\left(\mathrm{0,1}\right)$ and $\rho$ can take non-zero values. We use three values of $\rho ,$ $\rho \in \left\{\mathrm{0,0.1,0.2}\right\},$ in the sensitivity analysis and apply the same PFI and HDFI procedure that is based on the assumption that $x$ is an instrumental variable for $y_1.$ Such assumption is satisfied for $\rho \mathrm{<mo>=</mo>\; 0},$ but it is weakly violated for $\rho \mathrm{<mo>=</mo>\; 0.1}$ or $\rho \mathrm{<mo>=</mo>\; 0.2.}$ Using the fractionally imputed data in sample B, we estimated three parameters, ${\theta }_1\mathrm{<mo>=</mo>}E\left(Y_1\right),$ ${\theta }_2$ is the slope for the simple regression of $y_2$ on $y_1,$ and ${\theta }_3\mathrm{<mo>=</mo>}P\left(y_1\mathrm{<mo><</mo>\; 2,}{y}_2\mathrm{<mo><</mo>\; 3}\right),$ the proportion of $y_1\mathrm{<mo><</mo>\; 2}$ and $y_2\mathrm{<mo><</mo>\; 3.}$ Table 6.2 presents Monte Carlo means and variances of the point estimators for three parameters under three different models. In Table 6.2, the absolute values of the difference between the fractionally imputed estimator and the full sample estimator increase as the value of $\rho$ increases, which is expected as the instrumental variable assumption is more severely violated for larger values of $\rho ,$ but the differences are relatively small for all cases. In particular, the estimator of ${\theta }_1$ is not affected by the departure from the instrumental variable assumption. This is because the imputation estimator under the incorrect imputation model still provides an unbiased estimator for the population mean as long as the regression imputation model contains an intercept term (Kim and Rao 2012). Thus, this limited sensitivity analysis suggests that the proposed method seems to provide comparable estimates when the instrumental variable assumption is weakly violated.

6.2  Simulation two

In the second simulation study, we consider a binary response variable $y_{2i},$ where

$$y_{2i}\sim \text{Bernoulli}\left(p_i\right)\mathrm{,}(6.3)$$

$$\text{logit}\left(p_i\right)\mathrm{<mo>=</mo>}{\gamma }_0+{\gamma }_1{y}_{1i}\mathrm{,}$$

and $y_{1i}\sim N\left({\mu }_1\mathrm{,}{\sigma }_1^2\right).$ In the main sample, denoted by $B,$ instead of observing $\left(y_{1i}\mathrm{,}{y}_{2i}\right),$ we observe $\left(x_i\mathrm{,}{y}_{2i}\right),$ where

$$x_i\mathrm{<mo>=</mo>}{\beta }_0+{\beta }_1{y}_{1i}+u_i\mathrm{,}(6.4)$$

and $u_i\sim N\left(\mathrm{0,}{\sigma }^2{\left|y_{1i}\right|}^{2\alpha }\right).$ We observe $\left(x_i\mathrm{,}{y}_{1i}\right),$ $i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{n}_A$ in a calibration sample, denoted by A. For the simulation, $n_A\mathrm{<mo>=</mo>}{n}_B\mathrm{<mo>=</mo>\; 800},$ ${\gamma }_0\mathrm{<mo>=</mo>\; 1},$ ${\gamma }_1\mathrm{<mo>=</mo>\; 1},$ ${\beta }_0\mathrm{<mo>=</mo>\; 0},$ ${\beta }_1\mathrm{<mo>=</mo>\; 1},$ ${\sigma }^2\mathrm{<mo>=</mo>\; 0.25},$ $\alpha \mathrm{<mo>=</mo>\; 0.4},$ ${\mu }_1\mathrm{<mo>=</mo>\; 0},$ and ${\sigma }_1^2\mathrm{<mo>=</mo>\; 1.}$ Primary interest is in estimation of ${\gamma }_1$ and testing the null hypothesis that ${\gamma }_1\mathrm{<mo>=</mo>\; 1.}$ The Monte Carlo (MC) sample size is 1,000.

We compare the PFI estimators of ${\gamma }_1$ to three other estimators. Because the conditional distribution of $y_{1i}$ given $x_i$ is non-standard, we use the weights of (5.3) to implement PFI, where the proposal distribution ${\widehat{h}}_a\left(y_{1i}|x_i\right)$ is an estimate of the marginal distribution of $y_{1i}$ based on the data from Sample A. We consider the following four estimators:

1. PFI: For PFI, the proposal distribution for generating $y_{1i}^{\mathrm{<mo>*</mo>}\left(j\right)}$ is a normal distribution with mean ${\widehat{\mu }}_1$ and variance ${\widehat{\sigma }}_1^2,$ where ${\widehat{\mu }}_1$ and ${\widehat{\sigma }}_1^2$ are the maximum likelihood estimates of ${\mu }_1$ and ${\sigma }_1^2,$ respectively, based on Sample A. The fractional weights defined in (5.3) has the form
2. $$w_{ij}^{*}\propto {\widehat{p}}_{ij}^{y_{2i}}{\left(1-{\widehat{p}}_{ij}\right)}^{1-y_{2i}}{\widehat{f}}_a\left(y_{1i}^{\mathrm{<mo>*</mo>}\left(j\right)}|x_i\right)\mathrm{,}(6.5)$$
3. where ${\widehat{p}}_{ij}\mathrm{<mo>=</mo>}{\left\{1+\exp \left(-{\widehat{\gamma }}_0-{\widehat{\gamma }}_1{y}_{1i}^{\mathrm{<mo>*</mo>}\left(j\right)}\right)\right\}}^{-1}$ and ${\widehat{f}}_a\left(y_{1i}|x_i\right)$ is the estimate of $f\left(y_{1i}|x_i\right)$ based on maximum likelihood estimation with the Sample A data. The imputation size $m\mathrm{<mo>=</mo>\; 800.}$

2. Naive: A naive estimator is the estimator of the slope in the logistic regression of $y_{2i}$ on $x_i$ for $i\in B.$
3. Bayes: We use the approach of Guo and Little (2011) to define a Bayes estimator. The model for this simulation differs from the model of Guo and Little (2011) in that the response of interest is binary. We implement GIBBS sampling with JAGS (Plummer 2003), specifying diffuse proper prior distributions for the parameters of the model. Letting
4. $${\theta }_1\mathrm{<mo>=</mo>}\left(\text{log}\left({\sigma }^2\right)\mathrm{,}\text{log}\left({\sigma }_1^2\right)\mathrm{,}{\mu }_1\mathrm{,}{\beta }_0\mathrm{,}{\beta }_1\mathrm{,}{\gamma }_0\mathrm{,}{\gamma }_1\right)\mathrm{,}$$
5. we assume a priori that ${\theta }_1\sim N\left({\mathrm{0,10}}^6{I}_7\right),$ where $I_7$ is a $7\times 7$ identity matrix, and the notation $N\left(\mathrm{0,}V\right)$ denotes a normal distribution with mean 0 and covariance matrix $V.$ The prior distribution for the power $\alpha$ is uniform on the interval $\left[-\mathrm{5,5}\right].$
6. To evaluate convergence, we examine trace plots and potential scale reduction factors defined in Gelman, Carlin, Stern and Rubin (2003) for 10 preliminary simulated data sets. We initiate three MCMC chains, each of length 1,500 from random starting values and discard the first 500 iterations as burn-in. The potential scale reduction factors across the 10 simulated data sets range from 1.0 to 1.1, and the trace plots indicate that the chains mix well. To reduce computing time, we use 1,000 iterations of a single chain for the main simulation, after discarding the first 500 for burn-in.

4. A Weighted Regression Calibration (WRC) estimator. The WRC estimator is a modification of the weighted regression calibration estimator defined in Guo and Little (2011) for a binary response variable. The computation for the weighted regression calibration estimator involves the following steps:
5. (i)   Using ordinary least squares (OLS), regress $y_{1i}$ on $x_i$ for the calibration sample.
6. (ii)  Regress the logarithm of the squared residuals from step (i) on the logarithm of $x_i^2$ for the calibration sample. Let $\widehat{\lambda }$ denote the estimated slope from the regression.
7. (iii) Using weighted least squares (WLS) with weight ${\left|x_i\right|}^{2\widehat{\lambda }}\text{},$ regress $y_{1i}$ on $x_i$ for the calibration sample. Let ${\widehat{\eta }}_0$ and ${\widehat{\eta }}_1$ be the estimated intercept and slope, respectively, from the WLS regression.
8. (iv) For each unit $i$ in the main sample, let ${\widehat{y}}_{1i}\mathrm{<mo>=</mo>}{\widehat{\eta }}_0+{\widehat{\eta }}_1{x}_i.$
9. (v)  The estimate of $\left({\gamma }_0\mathrm{,}{\gamma }_1\right)$ is obtained from the logistic regression of $y_{2i}$ on ${\widehat{y}}_{1i}$ in the main sample.

Table 6.3 contains the MC bias, variance, and MSE of the four estimators of ${\gamma }_1.$ The naive estimator has a negative bias because $x_i$ is a contaminated version of $y_{1i}.$ The PFI estimator is superior to the Bayes and WRC estimators.

We compute an estimate of the variance of the PFI estimators of ${\gamma }_1$ using the variance expression based on the linear approximation. We define the MC relative bias as the ratio of the difference between the MC mean of the variance estimator and the MC variance of the estimator to the MC variance of the estimator. The MC relative bias of the variance estimators for PFI is negligible (less than 2% in absolute values).

• Previous
• Next