Integration of data from probability surveys and big found data for finite population inference using mass imputation
Section 4. Other techniques for mass imputation

• Previous
• Next

4.1   $k$ -nearest neighbor imputation

Instead of using a single imputed value, we now consider fractional imputation with $k$ imputed values for each missing outcome. Fractional imputation is designed to reduce the variance of the final estimator due to imputation (Kalton and Kish, 1984; Kim and Fuller, 2004).

Assume no matching ties, let $J_k\left(\text{}i\text{}\right)$ be the set of $k$ nearest neighbors for unit $i$

$$J_k\left(\text{}i\text{}\right)\mathrm{<mo>=</mo>}\left\{l\in B\mathrm{<mo>:</mo>}{\displaystyle \sum _{j\in B}{1}_{\left\{d\left(X_j\mathrm{,}{X}_i\right)\le d\left(X_l\mathrm{,}{X}_i\right)\right\}}}\le k\right\}\mathrm{<mo>=</mo>}\left\{i\left(\text{}1\text{}\right)\mathrm{,}\dots \mathrm{,}i\left(\text{}k\text{}\right)\right\}\mathrm{.}$$

The $k$ nearest neighbor approach to mass imputation can be described in the following steps: 

$${\widehat{\mu }}_{g\mathrm{,}\text{knn}}\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\in A}{\pi }_i^{-1}{\widehat{\mu }}_g\left(X_i\right)}\mathrm{.}(4.1)$$

In the non-parametric estimation literature, researchers have investigated the asymptotic properties of the $k$ nearest neighbor imputation estimators extensively. See, e.g., Mack and Rosenblatt (1979) and Mack (1981) for early references. Cheng (1994) establishes root- $n$ consistency of the $k$ nearest neighbor imputation estimator of the outcome mean when the outcome is subject to missingness. We derive the asymptotic theory for ${\widehat{\mu }}_{g\mathrm{,}\text{knn}}$ in the context of mass imputation combining a probability sample and a big data sample in the following theorem and defer its proof to the Supplementary Material.

Theorem 2. Under Assumptions 1 $–$ 4, $n{\left(k/N\right)}^{4/p}\to 0,$   $k/n\to 0,$  and $k^2/n\to \infty ,$

$$n^{1/2}\left({\widehat{\mu }}_{g\mathrm{,}\text{knn}}-{\mu }_g\right)\to N\left(\mathrm{0,}{V}_{\text{knn}}\right)\mathrm{,}(4.2)$$

where

$$V_{\text{knn}}\mathrm{<mo>=</mo>}\underset{n\to \infty }{{\displaystyle \lim }}\frac{n}{N^2}\left(E\left[{\text{var}}_p\left\{{\displaystyle \sum _{i\in A}{\pi }_i^{-1}{\mu }_g(X_i)}\right\}\right]+E\left\{\frac{1-{\pi }_B(X)}{{\pi }_B(X)}{\sigma }_g^2(X)\right\}\right)\mathrm{,}$$

and ${\pi }_B\left(\text{}X\text{}\right)\mathrm{<mo>=</mo>}P\left({\delta }_B\mathrm{<mo>=</mo>}1|X\right)$  and ${\sigma }_g^2\left(\text{}X\text{}\right)\mathrm{<mo>=</mo>}\text{var}\left\{g\left(\text{}Y\text{}\right)|X\right\}.$

If ${\pi }_B\left(\text{}X\text{}\right)$ goes to 1, $V_{\text{knn}}$ reduces to ${{\displaystyle \lim }}_{n\to \infty }\left(n/N^2\right)E\left[{\text{var}}_p\left\{{\displaystyle {\sum }_{i\in A}{\pi }_i^{-1}{\mu }_g\left(X_i\right)}\right\}\right].$ It suggests that if the big sample is a large fraction of the target population, $V_{\text{knn}}$ can be smaller than $V_{\text{nni}},$ suggesting that ${\widehat{\mu }}_{g\mathrm{,}\text{knn}}$ gains efficiency over ${\widehat{\mu }}_{g\mathrm{,}\text{nni}}.$ In finite samples, Beretta and Santaniello (2016) conduct a simulation study to compare nearest neighbor imputation and $k$ nearest neighbor imputation in the setting with independent and identically distributed data. They found that $k$ nearest neighbor imputation with a small $k$ outperforms nearest neighbor imputation in terms of mean squared error. On the one hand, a larger $k$ can use more information in the big data sample and leads to more efficiency gain; on the other hand, $k$ cannot be too large, in order to control the bias of our estimator. In practice, we suggest using data-driven methods, such as cross-validation, to choose a reasonable $k,$ and conducting sensitivity analysis varying the choice of $k.$

4.2   Generalized additive models

Nearest neighbor imputation methods are non-parametric. On the other hand, parametric models especially linear models are sensitive to model misspecification. We now consider semiparametric methods for mass imputation. Among semiparametric methods, generalized additive models (Hastie and Tibshirani, 1990) are flexible regarding model specification of the dependence of $Y$ on $X$ by specifying the model only through smooth functions rather than assuming a parametric relationship. As other non-parametric methods, the performance of generalized additive models will deteriorate as the dimension of $X$ becomes large. For $X$ with a moderate dimension, we apply generalized additive models to leverage the predictive power of the big data sample to produce a predictive model for $Y$ given $X,$ so as to facilitate mass imputation for the probability sample.

We assume that $g\left(Y_i\right)$ given $X_i$ follows some exponential family distribution, and

$$h^{-1}\left\{{\mu }_g\left(X_i\right)\right\}\mathrm{<mo>=</mo>}{f}_1\left(X_i^1\right)+f_2\left(X_i^2\right)+\dots f_p\left(X_i^p\right)\mathrm{,}(4.3)$$

where $h\left(\text{}\cdot \text{}\right)$ is an inverse link function, and each $f_k\left(\text{}\cdot \text{}\right)$ is a smooth function of $X^k,$ for $k\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}p.$ Model (4.3) allows for rather flexible specification of the dependence of $Y$ on $X.$ The estimated function $f_k\left(\text{}{X}^k\text{}\right)$ can reveal possible nonlinearities of the relationship of $Y$ and $X^k.$

There are several challenges in fitting model (4.3). First, $f_k\left(\text{}x\text{}\right)$ is an infinite-dimensional parameter, estimation of which often relies on some approximation. Second, we need to decide how smooth the $f_k\left(\text{}x\text{}\right)$ should be to balance the trade-off between model complexity and overfitting to the data at hand.

To solve the first issue, a common way to approximate $f_k\left(\text{}x\text{}\right)$ using splines. Let $B_m\left(\text{}x\text{}\right)$ be the basis spline functions for $m\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}M$ (Ruppert, Wand and Carroll, 2009). We approximate $f_k\left(\text{}x\text{}\right)$ by $f_k\left(\text{}x\text{}\right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{m\mathrm{<mo>=</mo>\; 1}}^M{\gamma }_m^k{B}_m\left(\text{}x\text{}\right)}$ with spline coefficients ${\gamma }_m^k.$ This leads to an approximation of model (4.3):

$$h^{-1}\left[\widehat{E}\left\{g\left(Y_i\right)|X_i\right\}\right]\mathrm{<mo>=</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^p{\displaystyle \sum _{m\mathrm{<mo>=</mo>\; 1}}^M{\gamma }_m^k{B}_m\left(X_i^k\right)}}\mathrm{.}(4.4)$$

In (4.4), a large $M$ allows for increased model complexity and also an increased chance of overfitting; while a small $M$ may result in an inadequate model. This trade-off is balanced by choosing a relatively large $M$ and then penalizing the model complexity in the estimation stage (Eilers and Marx, 1996). Let the vector of spline coefficients be ${\gamma }_k^{\text{T}}\mathrm{<mo>=</mo>}\left({\gamma }_1^k\mathrm{,}\dots \mathrm{,}{\gamma }_m^k\right)$ and ${\gamma }^{\text{T}}\mathrm{<mo>=</mo>}\left({\gamma }_1^{\text{T}}\mathrm{,}\dots \mathrm{,}{\gamma }_p^{\text{T}}\right).$ The estimate $\widehat{\gamma }$ is obtained by maximizing the penalized likelihood:

$$-2l\left(\gamma \right)+{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^p{\lambda }_k{\gamma }_k^{\text{T}}{S}_k{\gamma }_k}(4.5)$$

where $l\left(\gamma \right)$ is the log likelihood function of $\gamma ,$ $S_k$ is a matrix with the ${\left(m\mathrm{,}l\right)}^{\text{th}}$ component ${\displaystyle \int B_m^{{}^{\prime \text{}\prime }}\left(\text{}x\text{}\right)B_l^{{}^{\prime \text{}\prime }}\left(\text{}x\text{}\right)\text{d}x},$ ${\gamma }_k^{\text{T}}{S}_k{\gamma }_k$ regularizes $f_k$ to be smooth for which the degree of smoothness is controlled by ${\lambda }_k.$ Given the smoothing parameter ${\lambda }^{\text{T}}\mathrm{<mo>=</mo>}\left({\lambda }_1\mathrm{,}\dots \mathrm{,}{\lambda }_p\right),$ the penalized likelihood function in (4.5) is optimized by a penalized version of the iteratively reweighted least squares algorithm (Nelder and Baker, 1972; McCullagh, 1984) to obtain $\widehat{\gamma }.$ Regarding the choice of $\lambda ,$ we note that $\lambda$ controls the trade-off between model complexity and overfitting, which can be estimated separately from other model coefficients using generalized cross-validation or estimated simultaneously using restricted maximum likelihood estimation (Wood, 2006). In practice, the model performance is not sensitive to the choice of the number of basis functions as long as the number of basis functions is large relative to the sample size in the specification, but rather estimation of the smoothing parameter is critical to control the model complexity.

Once fitting the model, we can create an imputed value for each element $i$ in Sample A as

$${\widehat{\mu }}_{g\mathrm{,}\text{GAM}}\left(X_i\right)\mathrm{<mo>=</mo>}h\left\{{\widehat{f}}_1\left(X_i^1\right)+{\widehat{f}}_2\left(X_i^2\right)+\dots {\widehat{f}}_p\left(X_i^p\right)\right\}\mathrm{,}$$

where ${\widehat{f}}_k\left(\text{}x\text{}\right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{m\mathrm{<mo>=</mo>\; 1}}^M{\widehat{\gamma }}_m^k{B}_m\left(\text{}x\text{}\right)}$ for $k\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}p.$ The mass imputation estimator based on the generalized additive model is

$${\widehat{\mu }}_{g\mathrm{,}\text{GAM}}\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\in A}{\pi }_i^{-1}{\widehat{\mu }}_{g\mathrm{,}\text{GAM}}\left(X_i\right)}\mathrm{.}$$

Because in our context, the sample size of Sample B is much larger than that of Sample A, the estimation error in the imputation model can be negligible compared to the sampling variability of ${\widehat{\mu }}_{g\mathrm{,}\text{GAM}}.$

To close this subsection, it is worth commenting on the assumption of additive effects of $X$ in model (4.3). This assumption may be fairly strong one. To relax the additivity assumption, we can extend model (4.3) to include interactions through using the tensor product basis. For example, we can include a bivariate interaction surface $f_{12}\left(X^1\mathrm{,}{X}^2\right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{m\mathrm{<mo>=</mo>\; 1}}^M{\displaystyle {\sum }_{l\mathrm{<mo>=</mo>\; 1}}^L{\gamma }_{ml}{B}_m\left(\text{}{X}^1\text{}\right)B_l\left(\text{}{X}^2\text{}\right)}}.$ When using the tensor product basis, care should be taken with respect to the penalty function in order to result in appropriate effective degrees of freedom for the smoother. This topic has been investigated extensively in the literature; see, e.g., Wood (2006).

• Previous
• Article main page
• Next