Integration of data from probability surveys and big found data for finite population inference using mass imputation
Section 6. Empirical experiments

• Previous
• Next

In this section, we evaluate the finite sample performance of the proposed estimator using simulation studies, one based on artificial data using simple random sampling and the other based on a synthetic population file from a single month sample of the U.S. Census Bureau’s Monthly Retail Trade Survey using stratified sampling.

6.1   Kim-Wang example

We use the simulation example in Kim and Wang (2019) to compare various estimators. We generate the data according to the following mechanism. We first generate a finite population $F_N\mathrm{<mo>=</mo>}\left\{X_i\mathrm{<mo>=</mo>}\left(X_{1i}\mathrm{,}{X}_{2i}\right)\mathrm{,}{Y}_i\mathrm{<mo>=</mo>}(Y_{1i}\mathrm{,}{Y}_{2i}):i\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,N\right\}$ with size $N\mathrm{<mo>=</mo>}$ 1,000,000, where $Y_{1i}$ is a continuous outcome and $Y_{2i}$ is a binary outcome. From the finite population, we select a big data Sample B where the inclusion indicator ${\delta }_{Bi}~\text{Ber}\left(p_i\right)$ with $p_i$ the inclusion probability for unit $i$ with the sample size around 700,000. We obtain a representative Sample A of size $n\mathrm{<mo>=</mo>}$ 1,000 using simple random sampling. The parameters of interest are the population mean $N^{-1}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^N{Y}_i}$ and the conditional population mean of $Y_1$ given $Y_2\mathrm{<mo>=</mo>}1.$

For generating the finite population, we consider linear models

$$Y_{1i}\mathrm{<mo>=</mo>}1+X_{1i}+X_{2i}+{\alpha }_i+{\epsilon }_i\mathrm{,}(6.1)$$

$$P\left(Y_{2i}\mathrm{<mo>=</mo>}1|X_{1i}\mathrm{,}{X}_{2i}\mathrm{<mo>;</mo>}{\alpha }_i\right)\mathrm{<mo>=</mo>}\text{logit}\left(1+X_{1i}+X_{2i}+{\alpha }_i\right)\mathrm{,}$$

and nonlinear models

$$Y_{1i}\mathrm{<mo>=</mo>}\text{0}\text{.5}{\left(X_{1i}-\text{1}\text{.5}\right)}^2+X_{2i}^2+{\alpha }_i+{\epsilon }_i\mathrm{,}(6.2)$$

$$P\left(Y_{2i}\mathrm{<mo>=</mo>}1|X_{1i}\mathrm{,}{X}_{2i}\mathrm{<mo>;</mo>}{\alpha }_i\right)\mathrm{<mo>=</mo>}\text{logit}\left\{\text{0}\text{.5}{\left(X_{1i}-\text{1}\text{.5}\right)}^2+X_{2i}^2+{\alpha }_i\right\}\mathrm{,}$$

where $X_{1i}~N\left(\mathrm{1,}1\right),$ $X_{2i}~\text{Exp}\left(\text{}1\text{}\right),$ ${\alpha }_i~N\left(\mathrm{0,}1\right),$ ${\epsilon }_i~N\left(\mathrm{0,}1\right),$ and $X_{1i},$ $X_{2i},$ ${\alpha }_i$ and ${\epsilon }_i$ are mutually independent. The variables ${\alpha }_i$ induce the dependence of $Y_{1i}$ and $Y_{2i}$ even adjusting for $X_{1i}$ and $X_{2i}.$ For the big-data inclusion probability, we also consider a logistic linear model

$$\text{logit}\left(p_i\right)\mathrm{<mo>=</mo>}{X}_{2i}\mathrm{,}(6.3)$$

and a nonlinear logistic model

$$\text{logit}\left(p_i\right)\mathrm{<mo>=</mo>}-3+{\left(X_{1i}-\text{1}\text{.5}\right)}^2+{\left(X_{2i}-2\right)}^2\mathrm{.}(6.4)$$

We consider the following combinations: I. (6.1) and (6.3); II. (6.1) and (6.4); II. (6.2) and (6.3); and IV. (6.2) and (6.4) for data generating mechanisms. Therefore, the simulation setup is a $2\times 2$ factorial design with two levels in each factor.

Chen, Li and Wu (2020) propose the inverse propensity score weighting estimator using the estimated probability of selection into Sample B and the doubly robust estimator which further incorporates an outcome regression model. To evaluate the robustness and efficiency, we compare the following estimators: 

1. ${\widehat{\mu }}_{\text{HT}},$ the Horvitz $–$ Thompson estimator assuming the $Y_i’s$ were observed in Sample A for the purpose of benchmark comparison;
2. ${\widehat{\mu }}_{\text{ipw}},$ the inverse propensity score weighting estimator,

$${\widehat{\mu }}_{\text{ipw}}\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\in B}\frac{1}{p_i\left(\widehat{\eta }\right)}{Y}_i}\mathrm{,}$$

• where $p_i\left(\eta \right)\mathrm{<mo>=</mo>}P\left({\delta }_{Bi}\mathrm{<mo>=</mo>}1|X_{2i}\mathrm{<mo>;</mo>}\eta \right)$ is a logistic regression model with the linear predictor $X_{2i}$ with an unknown parameter $\eta ,$ and $\widehat{\eta }$ is an estimator of $\eta$ obtained by maximizing the modified likelihood function of $\eta$ (Chen et al., 2019) based on Samples A and B;

3. ${\widehat{\mu }}_{\text{dr}},$ the doubly robust estimator of Chen et al. (2019),

$${\widehat{\mu }}_{\text{dr}}\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\in B}\frac{1}{p_i\left(\widehat{\eta }\right)}\left(Y_i-X_i^{\text{T}}\widehat{\beta }\right)}+\frac{1}{n}{\displaystyle \sum _{i\in A}{X}_i^{\text{T}}\widehat{\beta }}\mathrm{,}$$

• where $\widehat{\beta }$ is the estimated regression coefficients using (6.1) as the working outcome regression model based on Sample B;

4. ${\widehat{\mu }}_{\text{nni}},$ the nearest neighbor imputation estimator;
5. ${\widehat{\mu }}_{\text{knn}},$ the $k$ nearest neighbor imputation estimator with $k\mathrm{<mo>=</mo>}5;$
6. ${\widehat{\mu }}_{\text{GAM}},$ the generalized additive model imputation estimator;
7. ${\widehat{\mu }}_{\text{RC}},$ the regression calibration estimator based on ${\widehat{\mu }}_{\text{nni}}$ with calibration variables $H\left({\delta }_B\mathrm{,}X\mathrm{,}Y\right)\mathrm{<mo>=</mo>}{\left({\delta }_B\mathrm{,}1-{\delta }_B\mathrm{,}{\delta }_{B}X\mathrm{,}{\delta }_{B}Y\right)}^{\text{T}}.$

All simulation results are based on 1,000 Monte Carlo runs. Table 6.1 summarizes the simulation results with biases, standard errors, and coverage rates of 95% confidence intervals using asymptotic normality of the point estimators. The following observations can be made from Table 6.1. ${\widehat{\mu }}_{\text{ipw}}$ has large biases when the propensity score is misspecified. ${\widehat{\mu }}_{\text{dr}}$ gains robustness over ${\widehat{\mu }}_{\text{ipw}}$ if one of the outcome regression model or the propensity score is correctly specified. However, if both models are misspecified, ${\widehat{\mu }}_{\text{dr}}$ has a larger bias. ${\widehat{\mu }}_{\text{nni}}$ has small biases across four scenarios, which shows its robustness. Importantly, the performance of ${\widehat{\mu }}_{\text{nni}}$ is close to that of ${\widehat{\mu }}_{\text{HT}}$ in terms of standard errors and coverage rates, which is consistent with our theory in Theorem 1. Moreover, as predicted by our theoretical results, ${\widehat{\mu }}_{\text{knn}}$ improves ${\widehat{\mu }}_{\text{nni}}$ in terms of efficiency. Also, ${\widehat{\mu }}_{\text{GAM}}$ shows robustness because of the flexibility of the model specification. The regression calibration estimator ${\widehat{\mu }}_{\text{RC}}$ has small biases across all scenarios and therefore shows robustness against model specifications for sampling score and outcome. Moreover, it has smaller standard errors than both ${\widehat{\mu }}_{\text{nni}}$ and ${\widehat{\mu }}_{\text{knn}}.$ The coverage rates are all close to the nominal level. 

6.2   Monthly retail trade survey

To demonstrate the practical relevance, we consider the U.S. Census Bureau’s 2014 Monthly Retail Trade Survey (Mulry, Oliver and Kaputa, 2014). The Monthly Retail Trade Survey is an economic indicator survey whose monthly estimates are inputs to the Gross Domestic Product estimates. This survey selects a sample of about 12,000 retail businesses each month with paid employees to collect data on sales and inventories. It employs an one-stage stratified sample with stratification based on major industry, further substratified by the estimated annual sales referred to as the size variable.

For simulation purpose, we use the simulated data from the 2014 Monthly Retail Trade Survey to suggest the data generating model and the true parameter values (https://ww2.amstat.org/meetings/ices/2016/contests.cfm). We generate a finite population of $N\mathrm{<mo>=</mo>}$ 812,765 retail businesses with 16 strata with a stratum identifier $h,$ sales $Y,$ inventories $X,$ and a size variable $Z$ on the log scale. Table 6.2 reports some summary statistics. We generate the inventory data from $X_{hi}~N\left({\mu }_{X\mathrm{,}h}\mathrm{,}{\sigma }_{X\mathrm{,}h}^2\right)$ for $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}{N}_h$ and $h\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}16,$ and the sales data from a linear model

$$Y_{hi}\mathrm{<mo>=</mo>}{\beta }_0+X_{hi}+{\epsilon }_{hi}\mathrm{,}(6.5)$$

and a nonlinear model

$$Y_{hi}\mathrm{<mo>=</mo>}{\beta }_0+\text{0}\text{.5}{X}_{hi}^2+{\epsilon }_{hi}\mathrm{,}(6.6)$$

where ${\epsilon }_{hi}~N\left(\mathrm{0,}\text{0}\text{.25}\right).$ In (6.5) and (6.6), we specify different values for ${\beta }_0$ so that the parameter of interest, $\mu \mathrm{<mo>=</mo>}{N}^{-1}{\displaystyle {\sum }_{h\mathrm{<mo>=</mo>\; 1}}^{16}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^{N_h}{Y}_{hi}}},$ matches with the true population mean 12.73.

We also generate a big data sample $S_B$ where the inclusion indicator ${\delta }_{hi}~\text{Ber}\left(p_{hi}\right)$ with the inclusion probability $p_{hi}$ for unit $i$ in stratum $h.$ The big data sample in practice is often available from E-commercial companies who monitor inventories and sales for retail businesses. For the big data inclusion probability, let $Z_{hi}~N\left({\mu }_{Z\mathrm{,}h}\mathrm{,}{\sigma }_{Z\mathrm{,}h}^2\right),$ for $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}{N}_h$ and $h\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}16.$ We consider a logistic linear model

$$\text{logit}\left(p_{hi}\right)\mathrm{<mo>=</mo>}{\alpha }_0+Z_{hi}\mathrm{,}(6.7)$$

and a nonlinear logistic model

$$\text{logit}\left(p_{hi}\right)\mathrm{<mo>=</mo>}{\alpha }_0+X_{hi}+Z_{hi}^2\mathrm{,}(6.8)$$

where we specify different values for ${\alpha }_0$ so that the mean inclusion probability is about 30%. Lastly, we generate a representative sample $S_A$ by stratified sampling with simple random sampling within strata without replacement; see Table 6.2 for the sample allocation.

We consider the seven estimators in Section 6.1 adopted for stratified sampling. In each mass imputed dataset, we apply the following point estimator and variance estimator: $\widehat{\mu }\mathrm{<mo>=</mo>}{N}^{-1}{\displaystyle {\sum }_{h\mathrm{<mo>=</mo>\; 1}}^H{N}_h{\overline{y}}_{n_h}}$ with ${\overline{y}}_{n_h}$ is the sample mean of $y$ in the $h^{\text{th}}$ stratum, $\widehat{V}\mathrm{<mo\; stretchy="false">(</mo>}\widehat{\mu }\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{<mo>=</mo>}{N}^{-2}{\displaystyle {\sum }_{h\mathrm{<mo>=</mo>\; 1}}^H{N}_h^2\left(1-n_h/N_h\right)s_{n_h}^2}/n_h$ with $s_{n_h}^2\mathrm{<mo>=</mo>}{\left(n_h-1\right)}^{-1}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^{n_h}{\mathrm{<mo\; stretchy="false">(</mo>}{y}_{hi}-{\overline{y}}_{n_h})}^2}.$

Table 6.3 summarizes the simulation results. A similar discussion to Section 6.1 applies. ${\widehat{\mu }}_{\text{ipw}}$ is sensitive to misspecification of the selection model; while ${\widehat{\mu }}_{\text{dr}}$ has double robustness feature, which still relies on at least one model to be correctly specified. Mass imputation based on nearest neighbor imputation, $k$ nearest neighbor imputation and generalized additive model shows good performances by leveraging the representativeness of the survey sample and the predictive power of the big data sample. In addition, if the big data membership is known throughout the survey data, the regression calibration estimator gains efficiency while maintaining the robustness against model misspecification. 

• Previous
• Article main page
• Next