Multiple imputation of missing values in household data with structural zeros
Section 6. Discussion

• Previous
• Next

The empirical study suggests that the NDPMPM can provide high quality imputations for categorical data nested within households. To our knowledge, this is the first parametric imputation engine for nested multivariate categorical data. The study also illustrates that, with modest sample sizes, agencies should not expect the NDPMPM to preserve all features of the joint distribution. Of course, this is the case with any imputation engine. For the NDPMPM, agencies may be able to improve accuracy for targeted quantities by recoding the data used to fit the model. For example, one can create a new household-level variable that equals one when everyone has the same race and equals zero otherwise, and replace the individual race variable with a new variable that has levels “1 = race is the same as race of household head”, “2 = race is white and differs from race of household head”, “3 = race is black and differs from race of household head”, and so on. The NDPMPM would be estimated with the household-level same race variable and the new individual-level race variable. This would encourage the NDPMPM to estimate the percentages with the same race very accurately, as it would be just another household-level variable like home ownership. It also would add structural zeros involving race to the computation. Evaluating the trade offs in accuracy and computational costs of such recodings is a topic for future research.

The NDPMPM can be computationally expensive, even with the speed-ups presented in this article. The expensive parts of the algorithm are the rejection sampling steps. Fortunately, these can be done easily by parallel processing. For example, we can require each processor to generate a fraction of the impossible cases in Section 2.2. We also can spread the rejection steps for the imputations over many processors. These steps should cut run time by a factor roughly equal to the number of processors available.

The empirical study used households up to size four. We have run the model on data with households up to size seven in reasonable time (a few hours on a standard laptop). Accuracy results are similar qualitatively. As the household sizes get large, the model can generate hundreds or even thousands times as many impossible households as there are feasible ones, slowing the algorithm. In such cases, the cap-and-weight approach is essential for practical applications.

Acknowledgements

This research was supported by grants from the National Science Foundation (NSF SES 1131897) and the Alfred P. Sloan Foundation (G-2-15-20166003).

Appendix

This is an Appendix to the paper. It contains proof that the rejection sampling step S9' in Section 3 generates samples from the correct posterior distribution. It also contains the modified Gibbs sampler for the cap-and-weight approach and a list of the structural zero rules used in fitting the NDPMPM model. Finally, we include empirical results for the speedup approaches mentioned in the paper, using synthetic data, and additional results for handling missing data using the NDPMPM under a missing completely at random scenario.

A.1 Proof that the rejection sampling step S9' in Section 3 generates samples from the correct posterior distribution

The $X_{ik}^1$ and $X_{ijk}^1$ values generated using the rejection sampler in Step S9' are generated from the full conditionals, resulting in a valid Gibbs sampler. The proof follows from the properties of rejection sampling (or simple accept reject). The target distribution is the full conditional for $X_i^{\text{mis}}\text{}.$ It can be re-expressed as

$$p\left(X_i^{\text{mis}}\right)\mathrm{<mo>=<mtext>\hspace{0.17em}</mtext></mo>}\frac{1\left\{X_i^1\notin {\mathcal{S}}_h\right\}}{\mathrm{Pr}\left(X_i\notin {\mathcal{S}}_h|\theta \right)}g\left(X_i^{\text{mis}}\right)$$

where

$$g\left(X_i^{\text{mis}}\right)={\pi }_{G_i^1}{\displaystyle \prod _{k|a_{ik}\mathrm{<mo>=</mo>\; 1}}^{p+q}}{\lambda }_{G_i^1{X}_{ik}^1}^{\left(k\right)}\left({\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{n_i}}{\omega }_{G_i^1{M}_{ij}^1}{\displaystyle \prod _{k|b_{ijk}\mathrm{<mo>=</mo>\; 1}}^p}{\varphi }_{G_i^1{M}_{ij}^1{X}_{ijk}^1}^{\left(k\right)}\right)\mathrm{.}$$

Our rejection scheme uses $g\left(X_i^{\text{mis}}\right)$ as a proposal for $p\left(X_i^{\text{mis}}\right).$ To show that the draws are indeed from $p\left(X_i^{\text{mis}}\right),$ we need to verify that $w\left(X_i^{\text{mis}}\right)\mathrm{<mo>=</mo>}p\left(X_i^{\text{mis}}\right)/g\left(X_i^{\text{mis}}\right)\mathrm{<mo><</mo>}M,$ where $\mathrm{1\; <mo><</mo>}M\mathrm{<mo><</mo>}\infty ,$ and that we are accepting each sample with probability $w\left(X_i^{\text{mis}}\right)/M\text{}.$ In our case,

1. $w\left(X_i^{\text{mis}}\right)\mathrm{<mo>=</mo>}p\left(X_i^{\text{mis}}\right)/g\left(X_i^{\text{mis}}\right)\mathrm{<mo>=</mo>}1\left\{X_i^1\notin {\mathcal{S}}_h\right\}/\mathrm{Pr}\left(X_i\notin {\mathcal{S}}_h|\theta \right)\le 1/\mathrm{Pr}\left(X_i\notin {\mathcal{S}}_h|\theta \right),$ and $\mathrm{0\; <mo><</mo>}\mathrm{Pr}\left(X_i\notin {\mathcal{S}}_h|\theta \right)\mathrm{<mo><</mo>\; 1}\Rightarrow \mathrm{1\; <mo><</mo>}1/\mathrm{Pr}\left(X_i\notin {\mathcal{S}}_h|\theta \right)\mathrm{<mo><</mo>}\infty$ necessarily.
2. By sampling until we obtain a valid sample that satisfies $X_i^1\notin {\mathcal{S}}_h,$ we are indeed sampling with probability $w\left(X_i^{\text{mis}}\right)/M\mathrm{<mo>=</mo>}1\left\{X_i^1\notin {\mathcal{S}}_h\right\}.$

A.2 Modified Gibbs sampler for the cap-and-weight approach

The modified Gibbs sampler for the cap-and-weight approach replaces steps S1, S3, S4, S5 and S6 of the Gibbs sampler in the main text as follows.

• S1*. For each $h\in \mathcal{H},$ repeat steps S1(a) to S1(e) as before but modify step S1(f) to: if $t_1\mathrm{<mo><</mo>}\lceil n_{1h}\times {\psi }_h\rceil ,$ return to step (b). Otherwise, set $n_{0h}\mathrm{<mo>=</mo>}{t}_0.$

• S3*. Set $u_F\mathrm{<mo>=</mo>\; 1.}$ Sample

$$u_g|-\sim \text{Beta}\left(1+U_g\text{}\mathrm{,}\alpha +{\displaystyle \sum _{f\mathrm{<mo>=</mo>}g+1}^F}{U}_f\right)\mathrm{,}{\pi }_g\mathrm{<mo>=</mo>}{u}_g{\displaystyle \prod _{f\mathrm{<mo><</mo>}g}}\left(1-u_f\right)$$

$$U_g\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^n}1\left(G_i^1\mathrm{<mo>=</mo>}g\right)+{\displaystyle \sum _{h\in \mathcal{H}}}\frac{1}{{\psi }_h}{\displaystyle \sum _{i|n_i^0\mathrm{<mo>=</mo>}h}}1\left(G_i^0\mathrm{<mo>=</mo>}g\right)$$

• S4*. Set $v_{gM}\mathrm{<mo>=</mo>\; 1}$ for $g\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}F.$ Sample

$$v_{gm}|-\sim \text{Beta}\left(1+V_{gm}\mathrm{,}\beta +{\displaystyle \sum _{s\mathrm{<mo>=</mo>}m+1}^S}{V}_{gs}\right)\mathrm{,}{\omega }_{gm}\mathrm{<mo>=</mo>}{v}_{gm}{\displaystyle \prod _{s\mathrm{<mo><</mo>}m}}\left(1-v_{gs}\right)$$

$$V_{gm}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^n}1\left(M_{ij}^1\mathrm{<mo>=</mo>}m\mathrm{,}{G}_i^1\mathrm{<mo>=</mo>}g\right)+{\displaystyle \sum _{h\in \mathcal{H}}}\frac{1}{{\psi }_h}{\displaystyle \sum _{i|n_i^0\mathrm{<mo>=</mo>}h}}1\left(M_{ij}^0\mathrm{<mo>=</mo>}m\mathrm{,}{G}_i^0\mathrm{<mo>=</mo>}g\right)$$

• S5*. Sample

$${\lambda }_g^{\left(k\right)}|-\sim \text{Dirichlet}\left(1+{\eta }_{g1}^{\left(k\right)}\text{}\mathrm{,}\dots \mathrm{,}1+{\eta }_{g{d}_k}^{\left(k\right)}\right)$$

$${\eta }_{gc}^{\left(k\right)}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i|G_i^1\mathrm{<mo>=</mo>}g}^n}1\left(X_{ik}^1=c\right)+{\displaystyle \sum _{h\in \mathcal{H}}}\frac{1}{{\psi }_h}{\displaystyle \sum _{i|\begin{array}{c}{n}_i^0\mathrm{<mo>=</mo>}h\mathrm{,}{G}_i^0\mathrm{<mo>=</mo>}g\end{array}}}1\left(X_{ik}^0\mathrm{<mo>=</mo>}c\right)$$

• S6*. Sample

$${\varphi }_{gm}^{\left(k\right)}|-\sim \text{Dirichlet}\left(1+{\nu }_{gm1}^{\left(k\right)}\mathrm{,}\dots \mathrm{,}1+{\nu }_{gm{d}_k}^{\left(k\right)}\right)$$

$${\nu }_{gmc}^{\left(k\right)}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i|\begin{array}{c}{G}_i^1\mathrm{<mo>=</mo>}g\mathrm{,}{M}_{ij}^1\mathrm{<mo>=</mo>}m\end{array}}^n}1\left(X_{ijk}^1\mathrm{<mo>=</mo>}c\right)+{\displaystyle \sum _{h\in \mathcal{H}}}\frac{1}{{\psi }_h}{\displaystyle \sum _{i|\begin{array}{c}{n}_i^0\mathrm{<mo>=</mo>}h\mathrm{,}{G}_i^0\mathrm{<mo>=</mo>}g\mathrm{,}{M}_{ij}^0\mathrm{<mo>=</mo>}m\end{array}}}1\left(X_{ijk}^0\mathrm{<mo>=</mo>}c\right)$$

A.3 List of structural zeros

We fit the NDPMPM model using structural zeros which involve ages and relationships of individuals in the same house. The full list of the rules used is presented in Table A.1. These rules were derived from the 2012 ACS by identifying combinations involving the relationship variable that do not appear in the constructed population. This list should not be interpreted as a “true” list of impossible combinations in census data.

A.4 Empirical study of the speedup approaches

We evaluate the performance of the two speedup approaches mentioned in the main text using synthetic data. We use data from the public use microdata files from the 2012 ACS, available for download from the United States Census Bureau (http://www2.census.gov/acs2012_1yr/pums/) to construct a population of 857,018 households of sizes $\mathcal{H}\mathrm{<mo>=</mo>}\left\{\mathrm{2,}\mathrm{3,}\mathrm{4,}\mathrm{5,}6\right\},$ from which we sample $n\mathrm{<mo>=</mo>}\text{10,000}$ households comprising $N\mathrm{<mo>=</mo>}\text{29,117}$ individuals. We work with the variables described in Table A.2. We evaluate the approaches using probabilities that depend on within household relationships and the household head.

We consider the NDPMPM using two approaches, both moving the values of the household head to the household level as in Section 4.1 of the main text and also using the cap-and-weight approach in Section 4.2 of the main text. The first approach considers ${\psi }_2\mathrm{<mo>=</mo>}{\psi }_3\mathrm{<mo>=</mo>}{\psi }_4\mathrm{<mo>=</mo>}{\psi }_5\mathrm{<mo>=</mo>}{\psi }_6\mathrm{<mo>=</mo>\; 1}$ while the second approach considers ${\psi }_2\mathrm{<mo>=</mo>}{\psi }_3\mathrm{<mo>=</mo>}1/2$ and ${\psi }_4\mathrm{<mo>=</mo>}{\psi }_5\mathrm{<mo>=</mo>}{\psi }_6\mathrm{<mo>=</mo>}1/3.$ We compare these approaches to the NDPMPM as presented in Hu et al., 2018. For each approach, we create $L\mathrm{<mo>=</mo>\; 50}$ synthetic datasets, $Z\mathrm{<mo>=</mo>}\left(Z^{\left(1\right)}\text{}\mathrm{,}\dots \mathrm{,}{Z}^{\left(50\right)}\right).$ We generate the synthetic datasets so that the number of households of size $h\in \mathcal{H}$ in each $Z^{\left(l\right)}$ exactly matches $n_h$ from the observed data. Thus, $Z$ comprises partially synthetic data (Little, 1993; Reiter, 2003), even though every released $Z_{ijk}$ is a simulated value. We combine the estimates using using the approach in Reiter (2003). As a brief review, let $q$ be the point estimator of some estimand $Q,$ and let $u$ be the estimator of variance associated with $q.$ For $l\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}L,$ let $q_l$ and $u_l$ be the values of $q$ and $u$ in synthetic dataset $Z^{\left(l\right)}\text{}.$ We use $\overline{q}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{l\mathrm{<mo>=</mo>\; 1}}^L}{q}_l/L$ as the point estimate of $Q$ and $T\mathrm{<mo>=</mo>}\overline{u}+b/L$ as the estimated variance of $\overline{q},$ where $b\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{l\mathrm{<mo>=</mo>\; 1}}^L}{\left(q_l-\overline{q}\right)}^2/\left(L-1\right)$ and $\overline{u}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{l\mathrm{<mo>=</mo>\; 1}}^L}{u}_l/L.$ We make inference about $Q$ using $\left(\overline{q}-Q\right)\sim t_v\left(\mathrm{0,}T\right),$ where $t_v$ is a $t-$ distribution with $v\mathrm{<mo>=</mo>}\left(L-1\right){\left(1+L\overline{u}/b\right)}^2$ degrees of freedom.

For each approach, we run the MCMC sampler for 20,000 iterations, discarding the first 10,000 as burn-in and thinning the remaining samples every five iterations, resulting in 2,000 MCMC post burn-in iterates. We create the $L\mathrm{<mo>=</mo>\; 50}$ synthetic datasets by randomly sampling from the 2,000 iterates. We set $F\mathrm{<mo>=</mo>\; 40}$ and $S\mathrm{<mo>=</mo>\; 15}$ for each approach based on initial tuning runs. For convergence, we examined trace plots of $\alpha$ , $\beta$ and weighted averages of a random sample of the multinomial probabilities in the NDPMPM likelihood. Across the approaches, the effective number of occupied household-level clusters usually ranges from 20 to 33 with a maximum of 38, while the effective number of occupied individual-level clusters across all household-level clusters ranges from 5 to 9 with a maximum of 12.

Based on MCMC runs on a standard laptop, moving household heads’ data values to the household level alone results in a speedup of about 63% on the default rejection sampler while the cap-and-weight approach alone results in a speedup of about 40%.

Table A.3 shows the 95% confidence intervals for each approach. Essentially, all three approaches result in similar confidence intervals, suggesting not much loss in accuracy from the speedups. Most intervals also are reasonably similar to confidence intervals based on the original data, except for the percentage of same age couples. The last row is a rigorous test of how well each method can estimate a probability that can be fairly difficult to estimate accurately. In this case, the probability that a household head and spouse are the same age can be difficult to estimate since each individual’s age can take 96 different values. All three approaches are thus off from the estimate from the original data in this case. These results suggest that we can significantly speedup the sampler with minimal loss in accuracy of estimates and confidence intervals of population estimands.

A.5 Empirical study of missing data imputation under MCAR

We also evaluate the performance of the NDPMPM as an imputation method under a missing completely at random (MCAR) scenario. We use the same data as in Section 5 of the main text. As a reminder, the data contains $n\mathrm{<mo>=</mo>}\text{5,000}$ households of sizes $\mathcal{H}\mathrm{<mo>=</mo>}\left\{\mathrm{2,}\mathrm{3,}4\right\},$ comprising $N\mathrm{<mo>=</mo>}\text{13,181}$ individuals. We introduce missing values using a MCAR scenario. We randomly select 80% households to be complete cases for all variables. For the remaining 20%, we let the variable “household size” be fully observed and randomly $–$ and independently $–$ blank 50% of each variable for the remaining household-level and individual-level variables. We use these low rates to mimic the actual rates of item nonresponse in census data.

Similar to the main text, we estimate the NDPMPM using two approaches, both combining the rejection step in Section 4.1 of the main text with the cap-and-weight approach in Section 4.2 of the main text. The first approach considers ${\psi }_2\mathrm{<mo>=</mo>}{\psi }_3\mathrm{<mo>=</mo>}{\psi }_4\mathrm{<mo>=</mo>\; 1}$ while the second approach considers ${\psi }_2\mathrm{<mo>=</mo>}{\psi }_3\mathrm{<mo>=</mo>}1/2$ and ${\psi }_4\mathrm{<mo>=</mo>}1/3.$ For each approach, we run the MCMC sampler for 10,000 iterations, discarding the first 5,000 as burn-in and thinning the remaining samples every five iterations, resulting in 1,000 MCMC post burn-in iterates. We set $F\mathrm{<mo>=</mo>\; 30}$ and $S\mathrm{<mo>=</mo>\; 15}$ for each approach based on initial tuning runs. We monitor convergence as in the main text. For both methods, we generate $L\mathrm{<mo>=</mo>\; 50}$ completed datasets, $Z\mathrm{<mo>=</mo>}\left(Z^{\left(1\right)}\text{}\mathrm{,}\dots \mathrm{,}{Z}^{\left(50\right)}\right),$ using the posterior predictive distribution of the NDPMPM, from which we estimate the same probabilities as in the main text.

Figures A.1 and A.2 display each estimated marginal, bivariate and trivariate probability ${\overline{q}}_{50}$ plotted against its corresponding estimate from the original data, without missing values. Figure A.1 shows the results for the NDPMPM with the rejection sampler, and Figure A.2 shows the results for the NDPMPM using the cap-and-weight approach. For both approaches, the NDPMPM does a good job of capturing important features of the joint distribution of the variables as the point estimates are very close to those from the data before introducing missing values. In short, the results are very similar to those in the main text, though more accurate.

Table A.4 displays 95% confidence intervals for selected probabilities involving within-household relationships, as well as the value in the full population of 764,580 households. The intervals include the two based on the NDPMPM imputation engines and the interval from the data before introducing missingness. The intervals are generally more accurate than those presented in the main text. This is expected since we use lower rates of missingness in the MCAR scenario. For the most part, the intervals from the NDPMPM with the two approaches tend to include the true population quantity. Again, the NDPMPM imputation engine results in downward bias for the percentages of households where everyone is the same race. As mentioned in the main text, this is a challenging estimand to estimate accurately via imputation, particularly for larger households.

References

Andridge, R.R., and Little, R.J.A. (2010). A review of hot deck imputation for survey non-response. International Statistical Review, 78(1), 40-64.

Bennink, M., Croon, M.A., Kroon, B. and Vermunt, J.K. (2016). Micro-macro multilevel latent class models with multiple discrete individual-level variables. Advances in Data Analysis and Classification.

Chambers, R., and Skinner, C. (2003). Analysis of Survey Data, Wiley Series in Survey Methodology, Wiley.

Dunson, D.B., and Xing, C. (2009). Nonparametric Bayes modeling of multivariate categorical data. Journal of the American Statistical Association, 104, 1042-1051.

Hu, J., Reiter, J.P. and Wang, Q. (2018). Dirichlet process mixture models for modeling and generating synthetic versions of nested categorical data. Bayesian Analysis, 13, 183-200.

Ishwaran, H., and James, L.F. (2001). Gibbs sampling methods for stick-breaking priors. Journal of the American Statistical Association, 161-173.

Kalton, G., and Kasprzyk, D. (1986). The treatment of missing survey data. Survey Methodology, 12, 1, 1-16. Paper available at https://www150.statcan.gc.ca/n1/en/pub/12-001-x/1986001/article/14404-eng.pdf.

Little, R.J.A. (1993). Statistical analysis of masked data. Journal of Official Statistics, 9, 407-426.

Manrique-Vallier, D., and Reiter, J.P. (2014). Bayesian estimation of discrete multivariate latent structure models with structural zeros. Journal of Computational and Graphical Statistics, 23, 1061-1079.

Murray, J.S., and Reiter, J.P. (2016). Multiple imputation of missing categorical and continuous values via Bayesian mixture models with local dependence (forthcoming). Journal of the American Statistical Association.

Raghunathan, T.E., and Rubin, D.B. (2001). Multiple imputation for statistical disclosure limitation. Technical Report.

Reiter, J.P. (2003). Inference for partially synthetic, public use microdata sets. Survey Methodology, 29, 2, 181-188. Paper available at https://www150.statcan.gc.ca/n1/pub/12-001-x/2003002/article/6785-eng.pdf.

Reiter, J.P., and Raghunathan, T.E. (2007). The multiple adaptations of multiple imputation. Journal of the American Statistical Association, 102, 1462-1471.

Rubin, D.B. (1976). Inference and missing data (with discussion). Biometrika, 63, 581-592.

Rubin, D.B. (1987). Multiple imputation for nonresponse in surveys, New York: John Wiley & Sons, Inc.

Rubin, D.B. (1993). Discussion: Statistical disclosure limitation. Journal of Official Statistics, 9, 462-468.

Savitsky, T.D., and Toth, D. (2016). Bayesian estimation under informative sampling. Electronic Journal of Statistics, 10.1, 1677-1708.

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica, 4, 639-650.

Si, Y., and Reiter, J.P. (2013). Nonparametric Bayesian multiple imputation for incomplete categorical variables in large-scale assessment surveys. Journal of Educational and Behavioral Statistics, 38.5, 199-521.

Vermunt, J.K. (2003). Multilevel latent class models. Sociological Methodology, 213-239.

Vermunt, J.K. (2008). Latent class and finite mixture models for multilevel data sets. Statistical Methods in Medical Research, 33-51.

Walker, S.G. (2007). Sampling the Dirichlet mixture model with slices. Communications in Statistics - Simulation and Computation, 1, 45-54.

Wang, Q., Akande, O., Hu, J., Reiter, J. and Barrientos, A. (2016). NestedCategBayesImpute: Modeling and Generating Synthetic Versions of Nested Categorical Data in the Presence of Impossible Combinations. The Comprehensive R Archive Network.

• Previous
• Article main page
• Next