Multivariate small area estimation under nonignorable nonresponse^*

* The opinions expressed in this paper are of the authors and do not necessarily represent the policies of the U.S. Bureau of Labor Statistics and the Israel Central Bureau of Statistics.

ABSTRACT

We consider multivariate small area estimation under nonignorable, not missing at random (NMAR) nonresponse. We assume a response model that accounts for the different patterns of the observed outcomes, (which values are observed and which ones are missing), and estimate the response probabilities by application of the Missing Information Principle (MIP). By this principle, we first derive the likelihood score equations for the case where the missing outcomes are actually observed, and then integrate out the unobserved outcomes from the score equations with respect to the distribution holding for the missing data. The latter distribution is defined by the distribution fitted to the observed data for the respondents and the response model. The integrated score equations are then solved with respect to the unknown parameters indexing the response model. Once the response probabilities have been estimated, we impute the missing outcomes from their appropriate distribution, yielding a complete data set with no missing values, which is used for predicting the target area means. A parametric bootstrap procedure is developed for assessing the mean squared errors (MSE) of the resulting predictors. We illustrate the approach by a small simulation study.

KEYWORDS:

• Distribution of missing data
• imputation under nonignorable nonresponse
• missing information principle
• MSE estimation
• NMAR nonresponse

1. Introduction, models and assumptions

Let $\{{\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};i=1,\dots ,M,j=1,\dots ,N_i\}$ represent the data in a finite population of $N$ units, belonging to $M$ areas, with $N_i$ units in area $i$, ${\sum }_{i=1}^M{N}_i=N$, where ${\mathbf{\text{y}}}_{ij}=(y_{ij,1}^{},\dots ,y_{ij,K}^{}{)}^{\mathrm{\prime }}$ is the vector of outcome values for unit $j$ in area $i$ and ${\mathbf{\text{x}}}_{ij}=(x_{ij,1}^{},\dots ,x_{ij,L}^{}{)}^{\mathrm{\prime }}$ is a vector of corresponding $L$ covariates. Note that the use of a single vector ${\mathbf{x}}_{ij}$ for the covariates accommodates situations where in fact different covariates, possibly of different dimension, apply to different observations. We assume that the covariates are known for every unit in the population, from a recent census or some administrative files. Suppose that the outcome values follow the generic two-level population model: (1) $\begin{array}{rl}& {\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i^U\stackrel{ind}{\sim }f({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i^U),i=1,\dots ,M,j=1,\dots ,N_i\\ & {\mathbf{\text{u}}}_i^U\stackrel{ind}{\sim }f\left({\mathbf{\text{u}}}_i^U\right);E\left({\mathbf{\text{u}}}_i^U\right)=\mathbf{\text{0}}=(0,\dots ,0{)}^{\mathrm{\prime }},V({\mathbf{\text{u}}}_i^U)={\mathrm{\Sigma }}^U,\end{array}$(1) where ${\mathbf{\text{u}}}_i^U=(u_{i,1}^U,\dots ,u_{i,K}^U{)}^{\mathrm{\prime }}$ is a K-dimensional latent random effect.

In the present article we assume that a noninformative sample has been drawn from the above population, but the observed data is incomplete because of not missing at random (NMAR) nonresponse. By noninformative sampling we mean that the sampling probabilities are not related to the outcome variable of interest after conditioning on the model covariates, such that the conditional distribution of the outcome variable in the sample, given the covariates, is the same as the corresponding distribution in the population from which the sample is taken.

In practice, the observed data in a sample are almost never complete due to non-response. The extent of the non-response may differ from unit to unit within an area, with some units providing all the requested information, while others only providing part of it, with different units answering different questions. And to make matters worse, the non-response is NMAR, that is, the probability of some target component of a unit being missing may depend, at least in part, on the missing target value, as well as the other target values for that unit, whether observed or missing. See e.g., Equation (10) for a simple example. As a consequence, approaches that ignore the non-response and just use the complete responses or those that model the non-response only as functions of the observed covariates may yield biased small area predictors. See the simulation study in Section 5.

As a practical example, consider the Household Expenditure Survey (HES) carried out by Israel's Central Bureau of Statistics. The survey collects information on socio-demographic characteristics, as well as information on income and expenditure. The sample consists of households selected with equal probabilities by a two-stage sampling design. Three important questions asked in this survey (and in other similar surveys across the world) relate to the salary in each of the three months preceding the month of the interview. Table 1 presents the distribution of the observed response patterns of the three variable in the 2017 survey, with “1” defining response and “0” nonresponse. The first position to the left defines the response regarding the salary in the month preceding the interview, the middle position defines the response regarding the salary 2 months before the interview, and the third position defines the response regarding the salary 3 months ago.

Table 1. Response patterns on 3 salary variables in Israel's HES. 2017.

Res. Pattern	000	001	010	011	100	101	110	111	Total
Count	885	23	14	308	20	9	40	9,664	10,963
Percentage	8.1	0.2	0.1	2.8	0.2	0.1	0.4	88.2	100

Pfeffermann and Sikov (2011) found that the response to salary questions is informative but they did not consider SAE and restricted to a single target variable. For further discussion and illustrations of NMAR nonresponse and related concepts, see, Rubin (1976), Little (1982), Little and Rubin (2002), Pfeffermann and Sikov (2011), and references therein.

Returning to the present article, the target is to impute the missing data and use the observed and missing data for estimating the small area means, or other summary measures of interest. It may come as a surprise, but we are not familiar with published articles considering small area estimation under NMAR nonresponse, except for Sverchkov and Pfeffermann (2018), which treats the case of univariate outcomes. The present paper extends the methodology developed in that article. See Pfeffermann and Sikov (2011) and Riddles, Kim, and Im (2016) for reviews and many references addressing the problem of NMAR nonresponse when fitting models to survey data, but with no attention to SAE applications.

Define the response indicator $R_{ij,k}=1\left(0\right)$ if $y_{ij,k}$ is observed (unobserved), and let ${\mathbf{\text{R}}}_{ij}=(R_{ij,1},\dots ,R_{ij,K}{)}^{\mathrm{\prime }}$.

Assumption 1.1:

(1a) The response occurs independently between the units,

(1b) $Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|({\mathbf{\text{y}}}_{i^{\ast }{j}^{\ast }},{\mathbf{\text{x}}}_{i^{\ast }{j}^{\ast }},{\mathbf{\text{u}}}_{i^{\ast }}^U),i^{\ast }=1,\dots .M,j^{\ast }=1,\dots ,N_i]=Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|{\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij}]$.

As noted in Sverchkov and Pfeffermann (2018), Assumption 1b is very reasonable. In particular, it states that the probability to respond to the target variable ${\mathit{y}}_{ij}$ does not depend on the corresponding random effect given ${\mathit{y}}_{ij}$, $Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|{\mathbf{\text{y}}}_{ij},u_i^U,{\mathbf{\text{x}}}_{ij}]=Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|{\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij}]$. Furthermore, it guarantees the identification of the response model. See Remark 2.2 in Section 2 for further discussion.

Note that under (1) and Assumption 1.1, (2) $\begin{array}{rl}& f[{\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i^U,{\mathbf{\text{R}}}_{ij},\{({\mathbf{\text{y}}}_{i^{\ast }{j}^{\ast }},{\mathbf{\text{x}}}_{i^{\ast }{j}^{\ast }},{\mathbf{\text{R}}}_{i^{\ast }{j}^{\ast }},{\mathbf{\text{u}}}_{i^{\ast }}^U),i^{\ast }\\ & =1\dots M,j^{\ast }=1\dots N_i;(i^{\ast },j^{\ast })\ne (i,j)\left\}\right]\\ & =f({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i^U,{\mathbf{\text{R}}}_{ij}).\end{array}$(2) We assume a parametric form for the completely observed outcomes, (3) $\begin{array}{rl}& {\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}}=(1,\dots ,1{)}^{\mathrm{\prime }}\sim f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i;{\theta }_1)\\ & =f({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}};{\theta }_1);\\ & {\mathbf{\text{u}}}_i={\mathbf{\text{u}}}_i^U-E({\mathbf{\text{u}}}_i^U|{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})\sim f_R({\mathbf{\text{u}}}_i;{\theta }_2)\\ & =f({\mathbf{\text{u}}}_i|{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}};{\theta }_2),E_R({\mathbf{\text{u}}}_i;{\theta }_2)=0.\end{array}$(3) Note that in general, ${\mathbf{\text{u}}}_i^U$ and ${\mathbf{\text{u}}}_i^{}$ are different if the nonresponse is NMAR.

Assumption 1.2:

The subset $\left\{\right(i,j):{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}}\}$ is not empty for every sampled area, such that the parameters $\theta =({\theta }_1,{\theta }_2)$ can be estimated by restricting to the fully observed data (units with no missing data), using classical small area estimation (SAE) procedures.

Remark 1.1:

Assumption 1.2 is for convenience and it is sufficient for our present approach to have fully observed data in only sufficient number of areas to allow efficient estimation of the parameters $\theta =({\theta }_1,{\theta }_2)$. Additionally, for a general response model under which the response to any given component of the multivariate target variable $\mathit{y}$ may depend on the component itself as well as the other components, with possibly different coefficients for each component, (see for example Equation (10) in Section 4), we also require sufficient number of observations for each response pattern ${\mathit{R}}_{\mathit{i}\mathit{j}}$, thus allowing efficient estimation of the response model for each component.

Denote by $\hat{\theta }=({\hat{\theta }}_1,{\hat{\theta }}_2)$ the estimate of $\theta$ obtained that way. For known $\theta$, the best predictor of the random effect ${\mathit{u}}_i$ given the completely observed data, $O_C=\left\{\right({\mathbf{\text{y}}}_{ij}:{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}}),{\mathbf{\text{x}}}_{ij},i=1,\dots ,M,j=1,\dots ,N_i\}$, is $E({\mathbf{\text{u}}}_i^{}|O_C;\theta )$. We predict, ${\widehat{\mathbf{\text{u}}}}_i^{}=E({\mathbf{\text{u}}}_i^{}|O_C;\theta =\hat{\theta })$.

Our proposed procedure to deal with the multivariate informative (NMAR) nonresponse consists of the following steps:

1. Fit a parametric model for the completely observed outcomes, (Equation (3)).

2. Fit an appropriate parametric model for the response probabilities, which may depend on the outcome and the covariates (Assumption 1b), indexed by the unknown vector parameter $\gamma$; $p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )=Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|{\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma ]$, with $p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )$ differentiable with respect to $\gamma$. See Section 2 for details.

3. Impute the missing outcomes from their appropriate distribution with the unknown parameters $({\theta }_1,{\theta }_2,\gamma )$ replaced by their sample estimates, and then use the ‘complete’ sample data (observed and imputed values), to predict the small area means or other area measures of interest. See Section 3 for the imputation equations under the model. Since we assume noninformative sampling such that if there was no nonresponse, the sample data would follow the same model as in the population, in what follows we do not distinguish between the population and sample data and consider the population data as our sample. The results of the present study can easily be generalised to the case where first a sample is selected from the finite population by some non-informative or informative sampling scheme, and then nonresponse occurs. In this case one can use the estimated distribution (3) and the estimated response model for imputation of the missing sample data as defined in the present article. Once the missing sample data are imputed, the small area means of interest can be estimated using the approach of Pfeffermann and Sverchkov (2007).

In the next section, we apply the MIP principle for estimating the response model parameters and discuss some related questions. In Section 3, we develop the imputation equations for the missing data, which, when combined with the observed data, permit simple estimation of the small area means or other area parameters of interest. In Section 4, we propose a parametric bootstrap procedure for estimating the prediction Root MSE of the resulting predictors. We illustrate our approach with a small simulation study in Section 5 and conclude with a summary of the main outcomes in Section 6.

2. Estimation of response model parameters

If the missing outcome values were actually observed, the vector parameter $\gamma$, indexing the response probabilities model, could be estimated by solving the likelihood equations: (4) $\sum _{\mathbf{\text{r}}=(0,\dots ,0{)}^{\mathrm{\prime }}}^{(1,\dots ,1{)}^{\mathrm{\prime }}}\sum _{(i,j):{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}}\frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\mathrm{\partial }\gamma }=0,$(4) where the external summation is over all the K-dimension vectors with 0,1 elements.

In practice, the missing data are unobserved for ${\mathbf{\text{R}}}_{ij}\ne \mathbf{\text{1}}$ and hence the likelihood equations (4) are not operational. However, one may apply in this case the missing information principle (MIP; Cepillini, Siniscialco, & Smith, 1955; Orchard & Woodbury, 1972). See, in particular, Sverchkov (2008), Sverchkov and Pfeffermann (2018), and Riddles et al. (2016) for recent applications of the principle to handle univariate NMAR nonresponse.

Missing Information Principle: Let $O=\left\{\right(y_{ij.k}:R_{ij,k}=1),{\mathbf{\text{x}}}_{ij},i=1,\dots ,M,j=1,\dots ,N_i\}$ denote all the observed data. Since no observations are available for elements $(ij,k):R_{ij,k}=0$, solve instead the best predictor of (4) given the observed data: (5) $\begin{array}{rl}& E\left(\left.\sum _{\mathbf{\text{r}}=(0,\dots ,0{)}^{\mathrm{\prime }}}^{(1,\dots ,1{)}^{\mathrm{\prime }}}\sum _{(i,j):{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}}\frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\mathrm{\partial }\gamma }\right|O\right)\\ & =E\left[\sum _{\mathbf{\text{r}}=(0,\dots ,0{)}^{\mathrm{\prime }}}^{(1,\dots ,1{)}^{\mathrm{\prime }}}\sum _{(i,j):{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}}\right.\\ & \times \left.\left.E\left(\left.\frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\mathrm{\partial }\gamma }\right|O,{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}\right)\right|O\right]=0.\end{array}$(5) The expectation $E\left(\left(\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )/\mathrm{\partial }\gamma \right)|O,{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}\right)$ can be approximated and solved as follows: Let $\alpha$ denote the set of indexes with observed values ${\mathbf{y}}_{ij,k}$ and $\text{β}$ denote the complement of $\alpha$, i.e., ${\mathbf{\text{y}}}_{ij,\alpha }=\{y_{ij,k};r_k=1\}$, ${\mathbf{\text{y}}}_{ij,\text{β}}=\{y_{ij,k};r_k=0\}$. Denote, ${\mathbf{\text{R}}}_{ij,\alpha }=(R_{ij,k}:k\in \alpha )$, ${\mathbf{\text{R}}}_{ij,\text{β}}=(R_{ij,k}:k\in \text{β})$ and define by ${\mathbf{\text{1}}}_{\text{β}},{\mathbf{\text{1}}}_{\alpha }$ the corresponding unit vectors of respective dimensions. By Assumption (1b), (6) $\begin{array}{rl}& E\left(\left.\frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\mathrm{\partial }\gamma }\right|O,{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}\right)\\ & =\int \frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\mathrm{\partial }\gamma }\\ & \times f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}})d{\mathbf{\text{y}}}_{ij,\text{β}}\\ & =\int \frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\mathrm{\partial }\gamma }\\ & \times \frac{\begin{array}{c}\{{[Pr({\mathbf{\text{R}}}_{ij,\text{β}}={\mathbf{\text{1}}}_{\text{β}}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij,\alpha }={\mathbf{\text{1}}}_{\alpha },{\mathbf{\text{y}}}_{ij}\left)\right]}^{-1}-1\}\\ f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}\end{array}}{\begin{array}{c}\int {[Pr({\mathbf{\text{R}}}_{ij,\text{β}}={\mathbf{\text{1}}}_{\text{β}}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij,\alpha }={\mathbf{\text{1}}}_{\alpha },{\mathbf{\text{y}}}_{ij}\left)\right]}^{-1}\\ f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}-1\end{array}};\\ & Pr({\mathbf{\text{R}}}_{ij,\text{β}}={\mathbf{\text{1}}}_{\text{β}}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij,\alpha }={\mathbf{\text{1}}}_{\alpha },{\mathbf{\text{y}}}_{ij})\\ & =\frac{p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\int p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}}.\end{array}$(6) Finally, solve (5) with respect to $\gamma$ by substituting $f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}};{\hat{\theta }}_1)=\left(f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i;{\hat{\theta }}_1)/\right.\left.\int f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i;{\hat{\theta }}_1)d{\mathbf{\text{y}}}_{ij,\text{β}}\right)$ for $f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})$, replacing ${\mathbf{\text{u}}}_i$ by ${\widehat{\mathbf{\text{u}}}}_i$ and dropping the external expectation. See Sverchkov and Pfeffermann (2018) for a similar approximation in the univariate case.

The last equality (product) in (6) extends to the multivariate case the following fundamental relationship between the sample and sample-complement distributions, derived in Sverchkov and Pfeffermann (2004) for the univariate case: (7) $\begin{array}{rl}& f(y_{ij}|x_{ij},u_i,R_{ij}=0)\\ & =\frac{\left[p_r^{-1}\right(y_{ij},x_{ij})-1]f(y_{ij}|x_{ij},u_i,R_{ij}=1)}{E\left\{\right[p_r^{-1}(y_{ij},x_{ij})-1]|x_{ij},u_i,R_{ij}=1\}}.\end{array}$(7) Equation (7) and its multivariate extension in Equation (6) form the basis for our proposed approach. It states that the distribution of an unobserved (missing) value $y_{ij}$ is defined mathematically by the distribution of $y_{ij}$ if it was observed, and the response model. Notice that under NMAR nonresponse, the distribution of $y_{ij}$ given that the unit responded is different from the distribution of $y_{ij}$ given that the unit did not respond, and also different from the population distribution of $y_{ij}$, before nonresponse takes place. The proof of the multivariate extension applied in (6) follows the same simple steps of the proof of (7) in Sverchkov and Pfeffermann (2004), utilising Bayes theorem. See also Sverchkov (2008) and Riddles et al. (2016).

In the Appendix, we illustrate the construction of Equation (6) under the mixed logistic model for the outcome variable.

Remark 2.1:

The dimension of the set of equations in (5) is equal to the dimension of $\gamma$ indexing the response model and hence it is impossible to estimate the parameters $\gamma$ and the parameters $\theta =({\theta }_1,{\theta }_2)$ of the outcome model defined by (3), by solely solving this set.

Remark 2.2:

A fundamental question regarding the use of the MIP equations (5) is the existence of a unique solution, or more generally, the identifiability of the response model. For the univariate case, Riddles et al. (2016) deal with NMAR nonresponse in the general context of sample surveys by following an approach proposed by Sverchkov (2008), which is similar to our present approach. Riddles et al. (2016) established the following fundamental condition for the response model identifiability: the covariates $\mathbf{\text{x}}$ can be decomposed as $\mathbf{\text{x}}=({\mathbf{\text{x}}}_{\mathbf{\text{1}}},{\mathbf{\text{x}}}_{\mathbf{\text{2}}}\mathbf{\text{)}}$, with $dim\text{(}{\mathbf{\text{x}}}_2)\ge 1$, such that $Pr(R_{ij}=1|y_{ij},{\mathbf{\text{x}}}_{ij})=Pr(R_{ij}=1|y_{ij},{\mathbf{\text{x}}}_{1ij})$. In other words, the covariates in ${\mathbf{\text{x}}}_2$ that appear in the outcome model do not affect the response probabilities, given the outcome and the other covariates. Covariates of this property may or may not exist in a general set up, but interesting enough, SAE models actually contain such a variable, namely, the random effects. The random effects play a fundamental role in SAE models so the outcome clearly depends on them, but it is reasonable to assume that the response probabilities do not depend on the random effects, given the outcome value, (which depends on the random effects). In practice, the random effects are unobservable but we estimate them and then solve the equations (5) by conditioning on the estimated effects. So, it is actually the estimated random effects that play the role of the covariates ${\mathbf{\text{x}}}_2$. In practice, other covariates that are predictive of the outcome but not of the response might exist as well.

3. Imputation of the missing data.

Once the parameters $\theta$ and $\gamma$ are estimated, the estimates can be substituted (together with ${\widehat{\mathbf{\text{u}}}}_i$) into the model holding for the missing data, using the relationship used in (6), yielding the following estimated distribution. Let ${\mathbf{\text{y}}}_{ij,\text{β}}=\{y_{ij,k};r_k=0\}$ define, as before, the unobserved data. (8) $\begin{array}{rl}& f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}};\hat{\gamma },\hat{\theta })\\ & =\frac{\begin{array}{c}\left[{\left(\frac{p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\hat{\gamma })}{\int p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\hat{\gamma })f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}}\right)}^{-1}-1\right]\\ \frac{f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i;{\hat{\theta }}_1)}{\int f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i;{\hat{\theta }}_1)d{\mathbf{\text{y}}}_{ij,\text{β}}}\end{array}}{\begin{array}{c}\int {\left(\frac{p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\hat{\gamma })}{\int p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\hat{\gamma })f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}}\right)}^{-1}\\ \frac{f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i;{\hat{\theta }}_1)}{\int f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i;{\hat{\theta }}_1)d{\mathbf{\text{y}}}_{ij,\text{β}}}d{\mathbf{\text{y}}}_{ij,\text{β}}-1\end{array}}.\end{array}$(8) Note again that the distribution $f_R({\mathbf{\text{y}}}_{ij}|{\mathbf{\text{x}}}_{ij},{\widehat{\mathbf{\text{u}}}}_i;{\hat{\theta }}_1)$ is of the observed data and can thus be estimated from the data using standard SAE model fitting procedures.

Imputation of the missing data can be carried out by drawing at random from the distribution (8). One may draw a single observation or multiple observations.

Once the missing observations are imputed, prediction of the true population mean of the outcome variable or other measures of interest is carried out by application of standard procedures. See the empirical study in Section 5.

Remark 2.3:

By Assumption 1.1, the response occurs independently between units.

4. Estimation of prediction MSE

As in any other statistical inference problem, one has to assess the error of the resulting predictors. In SAE applications under the frequentist paradigm, it is common to estimate the Root Prediction Mean Squared Error (RPMSE). It is quite obvious that no analytic expression of the RPMSE can be derived, given the complexity of the prediction procedure, and we therefore propose a bootstrap procedure. As before, we assume for convenience no sampling, such that the sample consists of all the population units. See Remark 4.1 below. The proposed bootstrap procedure consists of the following steps:

• B0. Impute the missing values as developed in Section 3. Consider the pseudo-population of complete responses as the ‘true’ population and calculate the corresponding true-pseudo area means.

• B1. For each unit $(i,j)$ with complete observation ${\mathit{y}}_{\mathit{i}\mathit{j}}^{\mathit{c}}$ generated in Step B0, draw observed outcomes with probabilities $p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{\mathit{i}\mathit{j}}^{\mathit{c}},{\mathbf{\text{x}}}_{ij};\hat{\gamma })$.

• B2. Apply all estimation and imputation procedures described in Sections 2 and 3 to the observed sample obtained in Step B1. Estimate all the area means.

• B3. Repeat Steps B1 and B2 independently B times (B large) and compute for each area $i$ the bootstrap RPMSE, (9) $\begin{array}{rl}RPMS{E}_{m,k}& =\frac{1}{B}{\sum }_{b=1}^B{({\hat{\overline{Y}}}_{m,k,b}-{\overline{Y}}_{m,k}^{B0})}^2;\\ m& =1,\dots ,M,b=1,\dots B,\end{array}$(9) where ${\hat{\overline{Y}}}_{m,k,b}$ is the predictor obtained from bootstrap sample $b$ for the mean of the k-th component of the outcome variable in area m and ${\overline{Y}}_{m,k}^{B0}$ is the corresponding pseudo mean in area $m$ as obtained in Step B0.

Remark 4.1:

The bootstrap procedure outlined above is partly design-based in the sense that we consider a single pseudo population and the models are used only for estimating the response probabilities and the model holding for the completely observed data. The procedure can easily be extended in two ways. First, we may generate a new pseudo population for each bootstrap sample, thus accounting also for the variability induced by the random generation of the population values. Second, we may extended the procedure to the case where a sample is selected from the population and nonresponse occurs in the sample, by first obtaining complete sample observations as in Step B0 and then generating a pseudo population using the procedure of Sverchkov and Pfeffermann (2004). Thereafter, a sample is drawn from the pseudo population with the same sampling design that was used for drawing the original sample. The other steps follow Steps B1-B3 above (with or without accounting for the generation of the pseudo population, i.e., by generating only one pseudo population or generating a new population each time).

5. Simulation study

In this section we describe the results of a simulation experiment when applying the procedures proposed in Sections 2, 3 and 4 (assuming no sampling and a single pseudo population).

The experiment consists of the following steps:

S1. Generation of population values: generate for each area $i,i=1,\dots ,300$ and for each unit $j,j=1,\dots ,50$ binary covariate values $x_{ij}$ with $Pr(x_{ij}=1)=Pr(x_{ij}=0)=0.5$, random effects ${\mathbf{\text{u}}}_i=(u_{i,1},u_{i,2}{)}^{\mathrm{\prime }}\sim N(\mathbf{\text{0}},\mathbf{\text{I}})$, $i=1,\dots ,300$, and corresponding independent outcome values from the mixed logistic model, (9) $\begin{array}{rl}& p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)=Pr(y_{ij,1}=1|x_{ij},{\mathbf{\text{u}}}_i)\\ & =\exp (-.1-x_{ij}+u_{i,1})/[1+\exp (-.1-x_{ij}+u_{i,1}\left)\right],\\ & p_{y_2}(x_{ij},{\mathbf{\text{u}}}_i)=Pr(y_{ij,2}=1|x_{ij},{\mathbf{\text{u}}}_i)\\ & =\exp (.9+u_{i,2})/[1+\exp (.9+u_{i,2}\left)\right]\end{array}$(9)

Remark 5.1:

The random effects are generated independently but they are not assumed to be independent in the estimation process.

S2. Response mechanism: compute response probabilities for unit $j$ in area $i$ as: (10) $\begin{array}{rl}& p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},x_{ij},\gamma )\\ & =\left\{\begin{array}{l}C(x_{ij},{\mathbf{\text{y}}}_{ij})\exp ({\gamma }_0+{\gamma }_{\text{1}}{x}_{ij}+{\gamma }_2{y}_{ij,1}+{\gamma }_3{y}_{ij,2}),\\ \text{if}\mathbf{\text{r}}=(1,1{)}^{\mathrm{\prime }}\\ C(x_{ij},{\mathbf{\text{y}}}_{ij})\exp ({\gamma }_4+{\gamma }_5{x}_{ij}+{\gamma }_6{y}_{ij,1}+{\gamma }_7{y}_{ij,2}),\\ \text{if}\mathbf{\text{r}}=(1,0{)}^{\mathrm{\prime }}\\ C(x_{ij},{\mathbf{\text{y}}}_{ij})\exp ({\gamma }_8+{\gamma }_9{x}_{ij}+{\gamma }_{10}{y}_{ij,1}+{\gamma }_{11}{y}_{ij,2}),\\ \text{if}\mathbf{\text{r}}=(0,1{)}^{\mathrm{\prime }}\\ C(x_{ij},{\mathbf{\text{y}}}_{ij}),\\ \text{if}\mathbf{\text{r}}=(0,0{)}^{\mathrm{\prime }}\end{array}\right.;\\ & C(x_{ij},{\mathbf{\text{y}}}_{ij})=[1+\exp ({\gamma }_0+{\gamma }_{\text{1}}{x}_{ij}+{\gamma }_2{y}_{ij,1}+{\gamma }_3{y}_{ij,2})\\ & +\exp ({\gamma }_4+{\gamma }_5{x}_{ij}+{\gamma }_6{y}_{ij,1}+{\gamma }_7{y}_{ij,2})\\ & +\exp ({\gamma }_8+{\gamma }_9{x}_{ij}+{\gamma }_{10}{y}_{ij,1}+{\gamma }_{11}{y}_{ij,2}){]}^{-1};\end{array}$(10) ${\gamma }_{\text{0}}\text{ = 0},{\gamma }_1\text{ = }-\text{.5},{\gamma }_2\text{ = 3},{\gamma }_3\text{ = }-\text{3},{\gamma }_4\text{ = 0},{\gamma }_5\text{ = }-\text{.5},{\gamma }_6\text{ = 2},{\gamma }_7\text{ = }-\text{2},{\gamma }_8\text{ = 0},{\gamma }_9\text{ = }-\text{.5},{\gamma }_{10}\text{ = 1},{\gamma }_{11}\text{ = }-\text{1}\text{.}$ Clearly, the nonresponse is NMAR since the response probabilities depend on the outcomes. Notice that the response for $y_{ij,1},y_{ij,2}$ is generated independently between units.

Remark 5.2:

We generated a single (finite) population and hence, a single set of response probabilities.

S3. Generating responses: generate responses from the (single) population generated in S1, with response probabilities defined in S2 (Equation (10)).

S4. Fitting respondents’ model: estimate ${\hat{p}}_{y_1}^{}(x_{ij},{\mathbf{\text{u}}}_i)=\hat{P}r(y_{ij,1}=1|x_{ij},{\widehat{\mathbf{\text{u}}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})$, ${\hat{p}}_{y_2}^{}(x_{ij},{\mathbf{\text{u}}}_i)=\hat{P}r(y_{ij,2}=1|x_{ij},{\widehat{\mathbf{\text{u}}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})$ by fitting the mixed logistic model (9), using PROC NLMIX in SAS with default options. Notice that the model (9) is not the true respondents’ model under the response model (10), because of the NMAR nonresponse.

S5. Estimation of response probabilities: assume the parametric response model (10), compute the expectations in (6) under the estimated models ${\hat{p}}_{y_1}^{}(x_{ij},{\widehat{\mathbf{\text{u}}}}_i)$, ${\hat{p}}_{y_2}^{}(x_{ij},{\widehat{\mathbf{\text{u}}}}_i)$ in Step S4 and estimate $\gamma$, using the procedure described in Section 2. See Sverchkov and Pfeffermann (2018) for numerical details.

S6. Imputation of missing data: impute the unobserved data from the distribution of the missing data defined in Section 3, which in the present case reduces to: $\begin{array}{rl}& f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}})\\ & =\frac{\begin{array}{c}\{{[Pr({\mathbf{\text{R}}}_{ij,\text{β}}={\mathbf{\text{1}}}_{\text{β}}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij,\alpha }={\mathbf{\text{1}}}_{\alpha },{\mathbf{\text{y}}}_{ij}\left)\right]}^{-1}-1\}\\ f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}\end{array}}{\begin{array}{c}\int {[Pr({\mathbf{\text{R}}}_{ij,\text{β}}={\mathbf{\text{1}}}_{\text{β}}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij,\alpha }={\mathbf{\text{1}}}_{\alpha },{\mathbf{\text{y}}}_{ij}\left)\right]}^{-1}\\ f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}-1\end{array}}.\end{array}$

Remark 5.3:

We imputed a single value for each missing value but one may impute several values, using a multiple imputation approach.

Repeat Steps S3–S6 independently 500 times.

Predictors considered: compute the following predictors for each area on each simulation. $\begin{array}{rl}1.{\hat{\overline{Y}}}_{i,1}^{ign}& =N_i^{-1}\left\{\sum _{j,R_{ij,1}=1}{y}_{ij,1}+\sum _{k=1,R_{ik,1}=0}^{N_i}{\hat{p}}_{y_1}(x_{ik},{\widehat{\mathbf{\text{u}}}}_i)\right\},\\ {\hat{\overline{Y}}}_{i,2}^{ign}& =N_i^{-1}\left\{\sum _{j,R_{ij,2}=1}{y}_{ij,2}+\sum _{k=1,R_{ik,2}=0}^{N_i}{\hat{p}}_{y_2}(x_{ik},{\widehat{\mathbf{\text{u}}}}_i)\right\}.\end{array}$ The predictors ${\hat{\overline{Y}}}_{i,1}^{ign},{\hat{\overline{Y}}}_{i,2}^{ign}$ ignore the response process and ‘assume’ that the population distribution holds also for the observed outcomes.

2. ${\hat{\overline{Y}}}_{i,1}^{new}=N_i^{-1}{\sum }_{j=1}^{N_i}{y}_{ij,1}^{imp}$, ${\hat{\overline{Y}}}_{i,2}^{new}=N_i^{-1}{\sum }_{j=1}^{N_i}{y}_{ij,2}^{imp}$, where $y_{ij,k}^{imp}=y_{ij,k}^{}$ if $y_{ij,k}^{}$ is observed, and $y_{ij,k}^{imp}$ is the imputed value from Step S6 if $y_{ij,k}^{}$ is missing ($k=1,2$).

The estimators ${\hat{\overline{Y}}}_{i,1}^{new},{\hat{\overline{Y}}}_{i,2}^{new}$ are our proposed estimators, accounting for the multivariate NMAR nonresponse.

Statistics considered for assessment of the of predictors’ performance

Denote by ${\overline{Y}}_{i,k,r}$ the true mean of area $i$ on the r-th simulation (for first or second coordinate, k = 1 or 2), and let ${\hat{\overline{Y}}}_{i,k,r}$ represent the first or second predictors defined above, $r=1,\dots ,500$. $\begin{array}{rl}Bia{s}_{i,k}& =\frac{{\sum }_{r=1}^{500}({\hat{\overline{Y}}}_{i,k,r}-{\overline{Y}}_{i,k,r})}{500};\\ RPMS{E}_{i,k}& =\frac{{\sum }_{r=1}^{500}{({\hat{\overline{Y}}}_{i,k,r}-{\overline{Y}}_{i,k,r})}^2}{500};\\ RelBia{s}_{i,k}& =\frac{Bia{s}_{i,k}}{\sqrt{V_{i,k}}};\\ V_{i,k}& =\frac{{\sum }_{r=1}^{500}{({\hat{\overline{Y}}}_{i,k,r}-\frac{1}{500}{\sum }_{r=1}^{500}{\hat{\overline{Y}}}_{i,k,r})}^2}{500};\\ RelRPMS{E}_{i,k}& =\frac{\sqrt{RPMS{E}_{i,k}}}{\left(\frac{1}{500}{\sum }_{r=1}^{500}{\overline{Y}}_{i,k,r}\right)}.\end{array}$ We calculated for each area the average (over the 500 simulations) of the number of complete responses and ordered the areas by these averages (the smallest mean number of complete responses is 2.3, the largest is 28.1).

S7. Estimation of the Root Prediction MSE (RPMSE): compute bootstrap estimates of RPMSE following the steps B0–B3 in Section 4.

In the following four figures we show the results for $RelBia{s}_{i,k}$ and $RelRMS{E}_{i,k}$, $k=1,2$ for each area, with the areas ordered as above, starting with the area with the smallest number of complete responses.

Figures 1 and 2 show how the proposed method reduces very significantly the bias due to NMAR nonresponse. As expected, the bias of both set of predictors decreases as the number of complete responses increases but our proposed predictors are seen to be much less biased.

Figure 1. $\text{Rel}Bia{s}_{i,1}$ of ${\hat{\overline{Y}}}_{i,1}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,1}^{new}$ (‘+’).

Figure 2. $\text{Rel}Bia{s}_{i,2}$ of ${\hat{\overline{Y}}}_{i,2}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,2}^{new}$ (‘+’).

The reduction in RelRMSE by accounting for the NMAR nonresponse in Figure 3 is not big, which is explained by the fact that the bias of the predictors that ignore the nonresponse is not very high in this case. Notice in this respect that the average number of missing values $y_{ij,1}$ over the 500 simulations is 5531.5, compared to an average number of 6014.6 missing values of $y_{ij,2}$. Nonetheless, when averaging the $RelRMS{E}_{i,1}$ over all the areas we find that, $\begin{array}{rl}& \text{Average[}RelRPMS{E}_{i,1}\left({\hat{\overline{Y}}}_{i,1}^{ign}\right)]\\ & =\frac{1}{300}\sum _{i=1}^{300}RelRPMS{E}_{i,1}\left({\hat{\overline{Y}}}_{i,1}^{ign}\right)=0.51,\\ & \text{Average[}RelRPMS{E}_{i,1}\left({\hat{\overline{Y}}}_{i,1}^{new}\right)]=\mathrm{0.44.}\end{array}$ When estimating ${\overline{Y}}_{i,2}$ in Figure 4, the reduction in the $RelRPMSE$ by use of the proposed procedure is much more drastic, particularly in the areas with small numbers of complete responses, due to the large bias when ignoring the NMAR nonresponse. $\begin{array}{rl}& \text{Average[}RelRPMS{E}_{i,2}\left({\hat{\overline{Y}}}_{i,2}^{ign}\right)]=0.93,\\ & \text{Average[}RelRPMS{E}_{i,2}\left({\hat{\overline{Y}}}_{i,2}^{new}\right)]=\mathrm{0.28.}\end{array}$ Next, we study the sensitivity of the proposed approach to correct specification of the response model. For this, we repeated the same simulation study, but by computing the response probabilities as: (11) $\begin{array}{rl}& p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},x_{ij},\gamma )\\ & =\left\{\begin{array}{l}C(x_{ij},{\mathbf{\text{y}}}_{ij})\exp [{\gamma }_0+{\gamma }_{\text{1}}{x}_{ij}({\gamma }_2{y}_{ij,1}+{\gamma }_3{y}_{ij,2}\left)\right],\\ \text{if}\mathbf{\text{r}}=(1,1{)}^{\mathrm{\prime }}\\ C(x_{ij},{\mathbf{\text{y}}}_{ij})\exp [{\gamma }_4+{\gamma }_5{x}_{ij}({\gamma }_6{y}_{ij,1}+{\gamma }_7{y}_{ij,2}\left)\right],\\ \text{if}\mathbf{\text{r}}=(1,0{)}^{\mathrm{\prime }}\\ C(x_{ij},{\mathbf{\text{y}}}_{ij})\exp \left[\right({\gamma }_8+{\gamma }_9{x}_{ij}({\gamma }_{10}{y}_{ij,1}+{\gamma }_{11}{y}_{ij,2})],\\ \text{if}\mathbf{\text{r}}=(0,1{)}^{\mathrm{\prime }}\\ C(x_{ij},{\mathbf{\text{y}}}_{ij}),\\ \text{if}\mathbf{\text{r}}=(0,0{)}^{\mathrm{\prime }}\end{array}\right.;\end{array}$(11) $\begin{array}{rl}C(x_{ij},{\mathbf{\text{y}}}_{ij})& =\{1+\exp [{\gamma }_0+{\gamma }_{\text{1}}{x}_{ij}({\gamma }_2{y}_{ij,1}+{\gamma }_3{y}_{ij,2})]\\ & +\exp [{\gamma }_4+{\gamma }_5{x}_{ij}({\gamma }_6{y}_{ij,1}+{\gamma }_7{y}_{ij,2}\left)\right]\\ & \exp [{\gamma }_8+{\gamma }_9{x}_{ij}({\gamma }_{10}{y}_{ij,1}+{\gamma }_{11}{y}_{ij,2}\left)\right]{\}}^{-1},\end{array}$ with the same coefficients as in (10).

Figure 3. $\text{RelRPMS}{\text{E}}_{i,1}$ of ${\hat{\overline{Y}}}_{i,1}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,1}^{new}$ (‘+’).

Figure 4. $\text{RelRPMS}{\text{E}}_{i,2}$ of ${\hat{\overline{Y}}}_{i,2}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,2}^{new}$ (‘+’).

With these response probabilities, the number of complete responses in an area (averaged over the 500 simulations), is in the range [9.3, 18.3].

When estimating the response model parameters in Step S5 of the simulation, we still use the model (10) as the working model, so that the model for the response is misspecified, and so is the model estimated for the missing data. (As mentioned before, the model estimated for the completely observed outcomes is also not correct).

Table 2 compares the true response probabilities (Equation (11)) with the average of the estimated response probabilities over the 500 simulations under the misspecified response model (Equation (10)). Notice that except in a few cases, the averages of the estimated response probabilities under the misspecified model are close to the true response probabilities, already illustrating lack of sensitivity of our proposed approach to correct specification of the response model, although the differences between the true response probabilities and their estimates are occasionally larger for any given sample.

Table 2. True response probabilities, $p_{r,(i,j)}$, and average of estimated response probabilities, $A{\hat{p}}_{r,(i,j)}$under misspecified response model, for different response patterns $r_{ij};(i=0,1,j=0,1)$.

$x$	$y_1$	$y_2$	$p_{r,(1,1)}$	$p_{r,(1,0)}$	$p_{r,(0,1)}$	$A{\hat{p}}_{r,(1,1)}$	$A{\hat{p}}_{r,(1,0)}$	$A{\hat{p}}_{r,(0,1)}$
0	0	0	.25	.25	.25	.25	.25	.25
0	0	1	.25	.25	.25	.31	.26	.21
0	1	0	.25	.25	.25	.19	.22	.30
0	1	1	.25	.25	.25	.24	.26	.24
1	0	0	.25	.25	.25	.30	.24	.24
1	0	1	.46	.27	.17	.36	.26	.19
1	1	0	.10	.17	.27	.23	.22	.29
1	1	1	.25	.25	.25	.28	.26	.23

As seen in Figures 5 and 6, with the misspecified response model (and the model for the completely observed data), there are no big biases even when ignoring the NMAR nonresponse. Nonetheless, even in this case, when averaging over all the areas, $\text{Average}\left[|Relbias\right({\hat{\overline{Y}}}_{i,1}^{ign}\left)|\right]=1.09$, $\text{Average}\left[|Relbias\right({\hat{\overline{Y}}}_{i,1}^{new}\left)|\right]=0.50$, $\text{Average}\left[|Relbias\right({\hat{\overline{Y}}}_{i,2}^{ign}\left)|\right]=1.17$, $\text{Average}\left[|Relbias\right({\hat{\overline{Y}}}_{i,2}^{new}\left)|\right]=0.47$.

Figure 5. $\text{Rel}Bia{s}_{i,1}$ of ${\hat{\overline{Y}}}_{i,1}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,1}^{new}$ (‘+’), response model misspecified.

Figure 6. $\text{Rel}Bia{s}_{i,2}$ of ${\hat{\overline{Y}}}_{i,2}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,2}^{new}$ (‘+’), response model misspecified.

Next we compare the RelRPMSEs of the two estimators.

Figures 7 and 8 show reduction in the RelRPMSEs when accounting for the NMAR nonresponse in the areas with small number of complete responses. When averaging over all the areas, $\begin{array}{rl}& \text{Average[}RelRPMS{E}_{i,1}\left({\hat{\overline{Y}}}_{i,1}^{ign}\right)]=0.41,\\ & \text{Average[}RelRPMS{E}_{i,1}\left({\hat{\overline{Y}}}_{i,1}^{new}\right)]=0.34;\\ & \text{Average[}RelRPMS{E}_{i,2}\left({\hat{\overline{Y}}}_{i,2}^{ign}\right)]=0.17,\\ & \text{Average[}RelRPMS{E}_{i,2}\left({\hat{\overline{Y}}}_{i,2}^{new}\right)]=\mathrm{0.14.}\end{array}$ We conclude that even under the misspecified models, our approach generally yields predictors with smaller RelRMSEs than when ignoring the NMAR nonresponse (Figures 5 and 6). Clearly, the predictors obtained under this approach have larger variances than when ignoring the NMAR nonresponse, due to all the complex computations involved, so that the large differences in the bias do not always translate into corresponding large differences in the RelRMSEs.

Figure 7. $\text{RelRPMS}{\text{E}}_{i,1}$ of ${\hat{\overline{Y}}}_{i,1}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,1}^{new}$ (‘+’), response model misspecified.

Figure 8. $\text{RelRPMS}{\text{E}}_{i,2}$ of ${\hat{\overline{Y}}}_{i,2}^{ign}$ (‘o’) and ${\hat{\overline{Y}}}_{i,2}^{new}$(‘+’), response model misspecified.

Finally, we report the results of RelRPMSE estimation. Due to time limitation, the results so far are based on only 100 parent samples and 50 bootstrap samples for each parent sample. Figures 9 and 10 compare the ‘true’ (empirical) RelRPMSEs over the 100 parent samples, with the mean of the corresponding bootstraps estimates.

Figure 9. $\text{RelRPMS}{\text{E}}_{i,1}$ of ${\hat{\overline{Y}}}_{i,1}^{new}$ (‘+’), and bootstrap estimates (‘o’).

Figure 10. $\text{RelRPMS}{\text{E}}_{i,2}$ of ${\hat{\overline{Y}}}_{i,2}^{new}$ (‘+’), and bootstrap estimates (‘o’).

The results in Figures 9 and 10 show for most areas good performance of the bootstrap estimators and we believe that with more parent samples and bootstrap samples, the results will look even better. Even with the current runs, when averaging over all the areas, $\begin{array}{rl}& \text{Average[}RelRPMS{E}_{i,1}]=0.38,\\ & \text{Average}\text{Bootstrap[}RelRPMS{E}_{i,1}]=0.35,\\ & \text{Average[}RelRPMS{E}_{i,2}]=0.41,\\ & \text{Average}\text{Bootstrap[}RelRPMS{E}_{i,2}]=0.41,\end{array}$ illustrating the unbiasedness of the bootstrap estimators when averaging over all the areas.

We compared the empirical RelRPMSE’s with the bootstrap estimates also for the case of the misspecified response model and obtained similar results. To save in space, we don’t show the corresponding figures.

6. Summary

In this paper we propose a general approach for multivariate SAE under NMAR nonresponse within the selected areas. The approach consists of fitting a model for the observed data and using this model for estimating a postulated multivariate response model by application of the missing information principle. Once the response model is estimated, we derive the model holding for the missing data, which is used for imputing the missing data, thus obtaining a complete file of sample data that is used for estimating the unknown small area parameters. A bootstrap procedure is proposed for estimating the root prediction mean squared errors of the small area predictors, which consists of generating a pseudo population with similar behaviour to the behaviour of the true underlying population, and selecting many samples from the pseudo population and many sets of responses for each sample.

A simulation study shows good performance of our approach in terms of point and RPMSE estimation. The simulation study also illustrates certain robustness to misspecification of the response model. The empirical study in this paper considers the case where the models that are fitted for the responding units and the response probabilities are logistic, but the theoretical derivations assume general models for the observed data and the response mechanism. Thus, we encourage researchers of SAE to apply the procedure to simulated and real data sets, with possibly different models assumed for the observed data and the response probabilities.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Notes on contributors

Danny Pfeffermann

Danny Pfeffermann is currently the National Statistician and General Director of the Central Bureau of Statistics of Israel. His main research areas are analytic inference from complex sample surveys and in particular, informative samples with non-ignorable nonresponse, small area estimation, seasonal adjustment and trend estimation and recently, accounting for mode effects and proxy surveys.

Michael Sverchkov

Michael Sverchkov is a Research Mathematical Statistician at the US Bureau of Labor Statistics. His main research areas are analytic inference from complex sample surveys and in particular, informative samples with non-ignorable nonresponse, small area estimation, seasonal adjustment.

Notes

* The opinions expressed in this paper are of the authors and do not necessarily represent the policies of the U.S. Bureau of Labor Statistics and the Israel Central Bureau of Statistics.

Appendix. Illustration of the use of Equation (6) for Estimation of the response probabilities:

Mixed logistic model for the outcome variables with a single covariate.

Consider bivariate variables ${\mathit{y}}_{\mathit{i}\mathit{j}}=(y_{ij1},y_{ij2})$, and suppose that the model fitted to the observed data of the respondents is the mixed generalised logistic model, (A1) $\begin{array}{rl}{p}_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)& =Pr(y_{ij,1}=1|x_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{\mathbf{\text{i}}j}=1)\\ & =\exp ({\alpha }_1+{\text{β}}_1{x}_{ij}+u_{i,1})\\ & \times [1+\exp ({\alpha }_1+{\text{β}}_1{x}_{ij}+u_{i,1}){]}^{-1}\\ p_{y_2}(x_{ij},{\mathbf{\text{u}}}_i)& =Pr(y_{ij,2}=1|x_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{\mathbf{\text{i}}j}=1)\\ & =\exp ({\alpha }_2+{\text{β}}_2{x}_{ij}+u_{i,2})\\ & \times [1+\exp ({\alpha }_2+{\text{β}}_2{x}_{ij}+u_{i,2}){]}^{-1}\\ {\mathbf{\text{u}}}_i& =(u_{i,1},u_{i,2}{)}^{\mathrm{\prime }}\sim N(\mathbf{\text{0}},\mathrm{\Sigma })\end{array}$(A1) Suppose a generic response model, $p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},x_{ij};\gamma )$=$Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|{\mathbf{\text{y}}}_{ij},x_{ij};\gamma ]$.

We assume that $y_{ij,1}$ and $y_{ij,2}$ are independent given $x_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{\mathbf{i}\mathbf{j}}=1$, and that $Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|{\mathbf{\text{y}}}_{ij},x_{ij},{\mathbf{\text{u}}}_i;\gamma ]=Pr[{\mathbf{\text{R}}}_{ij}=\mathbf{\text{r}}|{\mathbf{\text{y}}}_{ij},x_{ij};\gamma ]$.

Then, for example, for $\mathbf{\text{r}}=(0,1{)}^{\mathrm{\prime }}$, the components of (6) can be written as, $\begin{array}{rl}& \int \frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}({\mathbf{\text{y}}}_{ij},{\mathbf{\text{x}}}_{ij};\gamma )}{\mathrm{\partial }\gamma }\\ & \times \{[Pr({\mathbf{\text{R}}}_{ij,\text{β}}={\mathbf{\text{1}}}_{\text{β}}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij,\alpha }={\mathbf{\text{1}}}_{\alpha },{\mathbf{\text{y}}}_{ij}){]}^{-1}-1\}\\ & \times f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}\\ & =\frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}\left[\right(1,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma ]}{\mathrm{\partial }\gamma }\\ & \times \{[\frac{p_{\mathbf{\text{r}}}\left[\right(1,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma ]}{\begin{array}{c}{p}_{\mathbf{\text{r}}}\left[\right(1,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]p_{y_1}\right(x_{ij},{\mathbf{\text{u}}}_i)\\ +p_{\mathbf{\text{r}}}\left[\right(0,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]\right[1-p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)]\end{array}}{]}^{-1}-1\}\\ & \times p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)+\frac{\mathrm{\partial }\log p_{\mathbf{\text{r}}}\left[\right(0,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma ]}{\mathrm{\partial }\gamma }\\ & \times \{[\frac{p_{\mathbf{\text{r}}}\left[\right(0,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma ]}{\begin{array}{c}{p}_{\mathbf{\text{r}}}\left[\right(1,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]p_{y_1}\right(x_{ij},{\mathbf{\text{u}}}_i)\\ +p_{\mathbf{\text{r}}}\left[\right(0,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]\right[1-p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)]\end{array}}{]}^{-1}-1\}\\ & \times [1-p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i\left)\right],\\ & \int \{{[Pr({\mathbf{\text{R}}}_{ij,\text{β}}={\mathbf{\text{1}}}_{\text{β}}|{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij,\alpha }={\mathbf{\text{1}}}_{\alpha },{\mathbf{\text{y}}}_{ij}\left)\right]}^{-1}-1\}\\ & \times f({\mathbf{\text{y}}}_{ij,\text{β}}|{\mathbf{\text{y}}}_{ij,\alpha },{\mathbf{\text{x}}}_{ij},{\mathbf{\text{u}}}_i,{\mathbf{\text{R}}}_{ij}=\mathbf{\text{1}})d{\mathbf{\text{y}}}_{ij,\text{β}}\\ & =\{[\frac{p_{\mathbf{\text{r}}}\left[\right(1,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma ]}{\begin{array}{c}{p}_{\mathbf{\text{r}}}\left[\right(1,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]p_{y_1}\right(x_{ij},{\mathbf{\text{u}}}_i)\\ +p_{\mathbf{\text{r}}}\left[\right(0,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]\right[1-p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)]\end{array}}{]}^{-1}-1\}\\ & \times p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)\\ & +\{[\frac{p_{\mathbf{\text{r}}}\left[\right(0,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma ]}{\begin{array}{c}{p}_{\mathbf{\text{r}}}\left[\right(1,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]p_{y_1}\right(x_{ij},{\mathbf{\text{u}}}_i)\\ +p_{\mathbf{\text{r}}}\left[\right(0,y_{ij,2}{)}^{\mathrm{\prime }},{\mathbf{\text{x}}}_{ij};\gamma \left]\right[1-p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i)]\end{array}}{]}^{-1}-1\}\\ & \times [1-p_{y_1}(x_{ij},{\mathbf{\text{u}}}_i\left)\right]\end{array}$

Similar expressions are obtained for $\mathbf{\text{r}}=(1,0{)}^{\mathrm{\prime }}$ and $\mathbf{\text{r}}=(0,0{)}^{\mathrm{\prime }}$.