Response Model Selection in Small Area Estimation Under not Missing at Random Nonresponse

Abstract

Sverchkov and Pfeffermann^[16] consider Small Area Estimation (SAE) under informative probability sampling of areas and within the sampled areas, and not missing at random (NMAR) nonresponse. To account for the nonresponse, the authors assume a given response model, which contains the outcome values as one of the covariates and estimate the corresponding response probabilities by application of the Missing Information Principle, which consists of defining the likelihood as if there was complete response and then integrating out the unobserved outcomes from the likelihood by employing the relationship between the distributions of the observed and the missing data. A key condition for the success of this approach is the ‘correct’ specification of the response model. In this article, we consider the likelihood ratio test and information criteria based on the appropriate likelihood and show how they can be used for the selection of the response model. We illustrate the approach by a small simulation study.

AMS Subject Classification: 62D05, 62D10, 62F10

Keywords

AIC BIC information criteria likelihood ratio test population distribution respondents’ model sample distribution

1. Introduction

There exists almost no survey without nonresponse, but in practice most methods that deal with this problem assume either explicitly or implicitly that the missing data are ‘missing at random’ (MAR). However, in many practical situations, this assumption is not valid, since the probability to respond often depends on the outcome value, even after conditioning on available covariate information. In such cases, the use of methods that assume that the nonresponse is MAR can lead to large bias of parameter estimators and distort subsequent inference.

The case where the missing data are not MAR (NMAR) can be treated by postulating a parametric model for the distribution of the outcomes before non-response and a model for the response mechanism. These two models define a parametric model for the observed outcomes, so that the parameters of these models can be estimated from the observed data. See, for example, Pfeffermann and Sverchkov^[7] for details, with overview of related literature.

Modeling the distribution of the outcomes before non-response can be problematic since only the observed data are available. Sverchkov^[12] proposes an alternative approach, which allows to estimate the parameters of the response model without postulating a parametric model for the distribution of the outcomes before nonresponse. To account for the nonresponse, Sverchkov^[12] assumes a given response model and estimates the corresponding response probabilities by application of the missing information principle (MIP), which consists of defining the likelihood as if there was complete response, and then integrating out the unobserved outcomes from the likelihood, employing the relationship between the distributions of the observed and unobserved data. Sverchkov and Pfeffermann^[16] apply this approach for small area estimation (SAE) under informative probability sampling of areas and within the sampled areas, and NMAR nonresponse. We describe the main steps of this approach in Sections 2 and 3.

A key condition for the success of this approach is the ‘correct’ specification of the response model. In section 4 we consider the likelihood ratio test and information criteria based on the appropriate likelihood and show how they can be used for the selection of the response model. Section 5 illustrates the application of the approach by a small simulation study.

2. Notation and Models

Let ${y_{i j}, x_{i j}; i = 1, ..., M, j = 1, ..., N_{i}}$ represent the data in a finite population of N units, comprised of M areas with N_i units in area i, $\sum_{i = 1}^{M} N_{i} = N$ , where Y_ij is the value of the outcome variable for unit j in area i and ${x^{'}}_{i j} = (x_{i j, l}, ..., x_{i j, K})$ is a vector of corresponding K covariates. We assume that the covariates are known for every unit in the population. Suppose that the population outcome values follow the generic two-level model:

\begin{array}{l} y_{i j} | x_{i j}, u_{i}^{U} ~ f (y_{i j} | x_{i j}, u_{i}^{U}), i = 1, ..., M, j = 1, ..., N_{i} \\ u_{i}^{U} ~ f (u_{i}^{U}); E (u_{i}^{U}) = 0, V (u_{i}^{U}) = σ_{u^{U}}^{2}, \end{array}

(2.1)

where $u_{i}^{U}$ is the i^th area level random effect. The target is to estimate the area means ${\bar{Y}}_{i} = N_{i}^{- 1} \sum_{j = 1}^{N_{i}} y_{i j}, i = 1, ..., M$ , based on a sample obtained by the following two-stage sampling scheme: i) select a sample s of m out of the M population areas with inclusion probabilities $π_{i} = P r (i \in s)$ ; ii) select a sample s_i of $n_{i} > 0$ units from selected area i with probabilities $π_{j | i} = \Pr (j \in s_{i} | i \in s)$ . Denote by I_i, I_ij the sample indicators; $I_{i} = 1$ if area i is selected in the first stage and 0 otherwise, $I_{i j} = 1$ if unit j of selected area i is sampled in the second stage and $I_{i j} = 0$ otherwise. Let $w_{i} = 1 / π_{i}, w_{j | i} = 1 / π_{j | i}$ denote the first- and second-stage base sampling weights.

In practice, not every unit in the sample responds. Define the response indicator; $R_{i j} = 1$ if unit $j \in s_{i}$ responds and $R_{i j} = 0$ otherwise. The sample of respondents is thus $R = {(i, j) : I_{i} = 1, I_{i j} = 1, R_{i j} = 1}$ and the sample of nonrespondents among the sampled units is $R^{c} = {(i, k) : I_{i} = 1, I_{i k} = 1, R_{i k} = 0}$ . The response process is assumed to occur stochastically, independently between units. We assume $\sum_{j = 1}^{n_{i}} R_{i j} > 0$ in all the sampled areas. The sample of respondents defines therefore a third, self-selected stage of the sampling process with unknown response probabilities. (Särndal and Swensson^[11]).

Define, $u_{i} = u_{i}^{U} - E (u_{i}^{U} | i \in s)$ . Then, under the population model (2.1), the observed data follow the two-level ‘respondents’ model:

\begin{array}{l} f_{R} (y_{i j} | x_{i j}, u_{i}) = f (y_{i j} | x_{i j}, u_{i}, (i, j) \in R); \\ u_{i} ~ f (u_{i} | i \in s), E (u_{i} | i \in s) = 0. \end{array}

(2.2)

The model (2.2) is again general and all that we state at this stage is that under informative sampling and/or NMAR nonresponse, the population and the respondents’ models differ;

f_{R} (y_{i j} | x_{i j}, u_{i}) \neq f (y_{i j} | x_{i j}, u_{i}^{U})

Remark 1. The respondents’ model refers to the observed data and hence can be estimated and tested by standard SAE methods. See Pfeffermann^[3] and Rao and Molina^[9] for estimation and testing procedures in SAE, with references.

Let $p_{r} (y_{i j}, x_{i j}) = \Pr [R_{i j} = 1 | y_{i j}, x_{i j}, i \in s, j \in s_{i}] \Pr [R_{i j} = 1 | y_{i j}, x_{i j}, i \in s, j \in s_{i}]$ . If the probabilities $p_{r} (y_{i j}, x_{i j})$ were known, the sample of respondents could be considered as a two stage sample from the finite population with known sampling probabilities ύi and ${\tilde{π}}_{j | i} = π_{j | i} p_{r} (y_{i j}, x_{i j})$ . In this case, the area means ${\bar{Y}}_{i}$ can be estimated as in Pfeffermann and Sverchkov.^[6] Also, if known, the response probabilities could be used for imputation of the missing data within the selected areas, by application of the relationship between the sample and sample-complement distributions, (Sverchkov and Pfeffermann^[15]);

f (y_{i j} | x_{i j}, u_{i}, (i, j) \in R^{c}) = \frac{[p_{r}^{- 1} (y_{i j}, x_{i j}) - 1] f (y_{i j} | x_{i j}, u_{i}, (i, j) \in R)}{E {[p_{r}^{- 1} (y_{i j}, x_{i j}) - 1] | x_{i j}, u_{i}, (i, j) \in R}}

(2.3)

See Sverchkov and Pfeffermann,^[16] and Pfeffermann and Sverchkov^[8] for details.

3. Estimation of Response Probabilities

Unlike the sampling probabilities, the response probabilities are generally unknown. We assume therefore a parametric model, which is allowed to depend on the outcome and the covariate values; $\Pr [R_{i j} = 1 | y_{i j}, x_{i j}, i \in s, j \in s_{i}; γ] = p_{r} (y_{i j}, x_{i j}; γ)$ , where γ is a vector of unknown coefficients. We assume that $p_{r} (y_{i j}, x_{i j}; γ)$ is differentiable with respect to γ and satisfies the same mild regularity conditions as in Sverchkov and Pfeffermann.^[16]

Under these assumptions, if the missing outcome values were observed, γ could be estimated by solving the likelihood equations:

\sum_{(i, j) \in R} \frac{\partial \log p_{r} (y_{i j}, x_{i j}; γ)}{\partial γ} + \sum_{(i, k) \in R^{c}} \frac{\partial \log [1 - p_{r} (y_{i k}, x_{i k}; γ)]}{\partial γ} = 0

(3.1)

In practice, the missing data are unobserved and hence the likelihood equations (3.1) are not operational. However, one may apply in this case the missing information principle:

Missing Information Principle (MIP, Cepillini et al.^[1], Orchard and Woodbury^[2]): Let $O = {y_{i j}, n_{i}, (i, j) \in R; x_{h t}, h = 1, ..., M, t = 1, ..., N_{i}}$ represent the known observed data used below. Since no observations are available for (i, k)∈ R^c, solve instead,

E_{U} \overset{b y (2.3)}{=} \sum_{(i, j) \in R} \frac{\partial \log p_{r} (y_{i j}, x_{i j}; γ)}{\partial γ} \{[\sum_{(i, j) \in R} \frac{\partial \log p_{r} (y_{i j}, x_{i j}; γ)}{\partial γ} + \sum_{(i, k) \in R^{c}} \frac{\partial \log [1 - p_{r} (y_{i k}, x_{i k}; γ)]}{\partial γ}]| O\} + \sum_{(i, k) \in R^{c}} E_{s} (\frac{E_{r e} \{[p_{r}^{- 1} (y_{i k}, x_{i k}; γ) - 1] \frac{\partial \log [1 - p_{r} (y_{i k}, x_{i k}; γ)]}{\partial γ}| | x_{i k}, u_{i}, (i, k) \in R\}}{E_{r e} {[p_{r}^{- 1} (y_{i k}, x_{i k}; γ) - 1] | x_{i k}, u_{i}, (i, k) \in R}}| O) = 0.

(3.2)

See Sverchkov^[12] and Sverchkov and Pfeffermann^[16] for derivation of (3.2). In these equations, E_U,E_s,E_re define respectively expectations with respect to the population distribution, the sample distribution and the respondents’ distribution. Notice that the internal expectations in the last expression are with respect to the model holding for the observed data for the respondents.

Remark 2. When the response probabilities $p_{r} (y_{i j}, x_{i j}; γ)$ depend on only x_ij, they are referred to as propensity scores, and the missing data are missing at random. This kind of response mechanism may hold in establishment surveys, for example, when the response probability is related to the known size of the establishment. The estimating equations in (3.2) reduce in this case to the common log-likelihood equations,

\sum_{(i, j) \in R} \frac{\partial \log p_{r} (x_{i j}; γ)}{\partial γ} + \sum_{(i, k) \in R^{c}} \frac{\partial \log [1 - p_{r} (x_{i k}; γ)]}{\partial γ} = 0,

(3.3)

where $p_{r} (x_{i j}; γ) = P r (R_{i j} = 1 | x_{i j}; γ)$ .

Sverchkov and Pfeffermann^[16] propose to solve the equations (3.2) by maximizing the log-likelihood leading to them, that is, maximizing,

l (γ) = \sum_{(i, j) \in R} \log p_{r} (y_{i j}, x_{i j}; γ) + \sum_{(i, k) \in R^{c}} E_{s} (\frac{E_{r e} {[p_{r}^{- 1} (y_{i k}, x_{i k}; γ^{*}) - 1] l o g [1 - p_{r} (y_{i k}, x_{i k}; γ)] | x_{i k}, u_{i}, (i, k) \in R}}{E_{r e} {[p_{r}^{- 1} (y_{i k}, x_{i k}; γ^{*}) - 1] | x_{i k}, u_{i}, (i, k) \in R}}| O)

(3.4)

We distinguish between γ^* and γ because by (3.2), the derivatives should only be taken with respect to γ.

We maximize the likelihood (3.4) by replacing u_i by û_i, obtained by fitting a model of the form (2.2), and dropping the external expectation-E_s. The maximization is carried out iteratively, by maximizing in the (q+1) iteration the expression,

\sum_{(i, j) \in R} \log p_{r} (y_{i j}, x_{i j}; γ^{(q + 1)}) + \sum_{(i, k) \in R^{c}} \frac{E_{r e} {[p_{r}^{- 1} (y_{i k}, x_{i k}; γ^{(q)}) - 1] \log [1 - p_{r} (y_{i k}, x_{i k}; γ^{(q + 1)})] | x_{i k}, {\hat{u}}_{i}, (i, k) \in R}}{E_{r e} {[p_{r}^{- 1} (y_{i k}, x_{i k}; γ^{(q)}) - 1] | x_{i k}, {\hat{u}}_{i}, (i, k) \in R}}

(3.5)

with respect to γ^(q+1). The maximization can be carried out, for example, by SAS Proc NLIN. See Sverchkov^[14] and the examples following Remark 3 for details.

Remark 3. A fundamental question regarding the solution of the MIP equations is the existence of a unique solution or more generally, the identifiability of the response model. Riddles et al.^[10] propose a similar approach to deal with NMAR nonresponse in the general context of survey sampling inference and establish the following fundamental condition for the response model identifiability: the covariates x can be decomposed as x = (x₁,x₂) with $d i m (x_{2}) \geq 1$ , such that $\Pr (R_{i j} = 1 | y_{i j}, x_{i j}) = \Pr (R_{i j} = 1 | y_{i j}, x_{1 i j})$ . In other words, the covariates in x₂ that appear in the outcome model do not affect the response probabilities, given the outcome and the other covariates. Variable(s) 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 effect, given the outcome value. In practice, the random effects are unobservable, but we estimate them and then solve the equations (3.5) by conditioning on the estimated effects. So, it is actually the estimated random effects that play the role of the covariates $x_{2}$ . (Other covariates that are predictive of the outcome but not of the response might exist as well).

Clearly, the larger are the absolute values of the random effects, the more they affect the values of the outcome values and hence also the values of the response probabilities. In the simulation study of Sverchkov and Pfeffermann,^[16] the authors study the effect of the magnitude of the variance of the random effects on the prediction of the area means. The conclusions from that study is that although the estimators of response model parameters become biased as the variance of the random effects increases, the biases are relatively very small and so are the standard deviations of the estimators. Furthermore, increasing the variance of the random effects has negligible effect on the estimation of the true response probabilities and the predictors of the true small area means remain virtually unbiased in each of the areas.

Riddles et al.^[10] prove asymptotic normality of the estimator $\hat{γ}$ under general regularity conditions.

Example 1. (Sverchkov and Pfeffermann^[16]): Mixed logistic model for the outcome variable.

Suppose that the model fitted to the observed data of the respondents is the mixed generalized logistic model,

p_{y} (x_{i j}, u_{i}) = \Pr (y_{i j} = 1 | x_{i j}, u_{i}, (i, j) \in R; β) = \frac{\exp (β_{0} + β_{1} x_{i j} + u_{i})}{1 + \exp (β_{0} + β_{1} x_{i j} + u_{i})}, u_{i} \overset{i.i.d.}{~} N (0, σ_{u}^{2}) .

Consider a generic response model, $p_{r} (y_{i j}, x_{i j}; γ) = \Pr [R_{i j} = 1 | y_{i j}, x_{i j}, i \in s, j \in s_{i}; γ]$ .

The components of (3.2) can be written in this case as,

\begin{array}{l} E_{r e} \{[p_{r}^{- 1} (y_{i j}, x_{i j}; γ) - 1] \frac{\partial \log [1 - p_{r} (y_{i j}, x_{i j}; γ)]}{\partial γ}| x_{i j}, u_{i}, (i, j) \in R\} = \\ p_{y} (x_{i j}, u_{i}) [p_{r}^{- 1} (1, x_{i j}; γ) - 1] \frac{\partial \log [1 - p_{r} (1, x_{i j}; γ)]}{\partial γ} + \\ [1 - p_{y} (x_{i j}, u_{i})] [p_{r}^{- 1} (0, x_{i j}; γ) - 1] \frac{\partial \log [1 - p_{r} (0, x_{i j}; γ)]}{\partial γ}; \end{array}

(3.6)

E_{r e} {[p_{r}^{- 1} (y_{i j}, x_{i j}; γ) - 1] | x_{i j}, u_{i}, (i, j) \in R} = p_{y} (x_{i j}, u_{i}) [p_{r}^{- 1} (1, x_{i j}; γ) - 1] + [1 - p_{y} (x_{i j}, u_{i})] [p_{r}^{- 1} (0, x_{i j}; γ) - 1]

(3.7)

The random effects u_i and the logistic probabilities p_y(x_ij,u) can be estimated by use of the SAS procedure PROC NLMIX.

Example 2. (Sverchkov^[14]): General continuous model.

In Example 1, the outcomes follow a discrete distribution. In this section, we consider continuous outcomes. The proposed algorithm consists of two parts:

Part 1: Fit (estimate) the model (2.2). The output of this part (input for Part 2) contains the model parameter estimates, the estimated random effects, û_i, and for each (i,j) R, estimates of $p_{y}^{(l)} (x_{i j}, {\hat{u}}_{i}) = P_{R} (a_{l} \leq y_{i j} < a_{l + 1} | x_{i j}, {\hat{u}}_{i}, (i, j) \in R)$ , $l = 0, ..., L + 2$ ; $a_{0} = - \infty$ , $a_{L + 2} = \infty$ , $a_{l} = \min (y_{i j}) + (l - 1) \frac{\max (y_{i j}) - \min (y_{i j})}{L}$ , $l = 1, ..., L + 1$ . The max and min are over all the observed values y_ij.

Part 2: Approximate the expectations in (3.2) similarly to (3.6) and (3.7):

\begin{array}{l} E_{r e} \{[p_{r}^{- 1} (y_{i j}, x_{i j}; γ) - 1] \frac{\partial \log [1 - p_{r} (y_{i j}, x_{i j}; γ)]}{\partial γ}| x_{i j}, u_{i}, (i, j) \in R\} ≅ \\ \sum_{l = 1}^{L + 1} {\hat{p}}_{y}^{(l)} (x_{i j}, {\hat{u}}_{i}) [p_{r}^{- 1} (a_{l}, x_{i j}; γ) - 1] \frac{\partial \log [1 - p_{r} (a_{l}, x_{i j}; γ)]}{\partial γ} \end{array}

(3.8)

E_{r e} \{[p_{r}^{- 1} (y_{i j}, x_{i j}; γ) - 1]| x_{i j}, u_{i}, (i, j) \in R\} ≅ \sum_{l = 1}^{L + 1} {\hat{p}}_{y}^{(l)} (x_{i j}, {\hat{u}}_{i}) [p_{r}^{- 1} (a_{l}, x_{i j}; γ) - 1]

(3.9)

where $p_{r} (a_{l}, x_{i j}; γ) = \Pr [R_{i j} = 1 | y_{i j} = a_{l}, x_{i j}, i \in s, j \in s_{i}; γ]$ .

Substitute (3.8) and (3.9) into (3.2) and estimate γ by iteratively maximizing (3.5).

4. Selection of the Response Model

There is no direct way to test the appropriateness of a chosen response model because the outcome values, which are part of the model, are unknown for the nonresponding units. If the model for the outcomes before nonresponse was known, one could derive the model holding for the observed outcomes based on this model and the model assumed for the responding units, and test the resulting model by use of standard tests that compare the cumulative hypothesized distribution of the observed data with the corresponding empirical distribution, and/or by testing moments of the assumed model. See, for example, Pfeffermann and Landsman^[4] and Pfeffermann and Sikov.^[5] However, in the approach described in Section 3, we start with a model fitted to the observed outcomes, which does not include the response model and therefore, we cannot use a similar strategy.

When following the approach proposed in Section 3, the likelihood (3.4) suggests at least two procedures for the selection of the response model in SAE under NMAR nonresponse. (a) Compare different models based on information criteria such as the Akaike information criterion, AIC = –2l() + 2dim(γ), or Schwarz information criterion, $BIC= - 2 l (γ) + \dim (γ) \log (n), n = \sum_{i \in s} n_{i}$ ; (b) test a saturated versus a nested model based on the likelihood ratio test. In Section 5 we illustrate via a simulation study how the likelihood (3.4) can be used for the application of these selection procedures.

5. Simulation Study

5.1 Simulation Set-up

We start by defining the sample model before nonresponse because as stated in Section 4, our approach for estimating the response model is based on fitting a model to the observed outcomes, which does not include the response model. For convenience, we assume noninformative sampling of areas and within the areas, such that the sample model before nonresponse is the same as the population model. Note that although the sampling design defines the observed model (2.2), once this model is estimated, the sampling design does not affect the estimation of the response probabilities in Section 3.

The simulation study consists of the following steps:

Generate auxiliary values $x_{i j}, i = 1, ..., 100, j = 1, ..., 20$ from a Uniform(0,2) distribution. Next, generate sample values from the small area model,

y_{i j} | x_{i j}, u_{i} ~ N (x_{i j} + u_{i}^{}, 1), i = 1, ..., 100, j = 1, ..., 20; u_{i} ~ N (0, 1)

(5.1)

Consider three unit response models (no selection of areas):

p_{r}^{(1)} (y_{i j}, x_{i j}) = logit (- x_{i j} / 2 + 2 y_{i j}),

p_{r}^{(2)} (y_{i j}, x_{i j}) = logit (- x_{i j} / 2 + 2 y_{i j} - 0.3 y_{i j}^{2}),

p_{r}^{(3)} (y_{i j}, x_{i j}) = logit (1.5 x_{i j}) .

Select 3 sets of respondents:

R1 uses Poisson sampling, independently between the units with response probabilities $p_{r}^{(1)} (y_{i j}, x_{i j})$ ,

R2 is the same as R1 but with response probabilities $p_{r}^{(2)} (y_{i j}, x_{i j})$ ,

R3 is the same as R1 but with response probabilities $p_{r}^{(3)} (y_{i j}, x_{i j})$ .

The 3 response probabilities yield similar response rates of 65–75 per cent.

The working model for the observed data for the responding units is,

y_{i j} | x_{i j}, u_{i} ~ N (θ_{0} + θ_{1} x_{i j} + u_{i}, θ_{2}), i = 1...100, j = 1...20; u_{i} ~ N (0, σ_{u}^{2}) .

(5.2)

Remark 4: The working model (5.2) is correct for the observed sample R3 that corresponds to MAR nonresponse, but not for R1 and R2, under which the nonresponse is NMAR.

Define three working response models:

M1: p_{r} (y_{i j}, x_{i j}; γ^{1}) = logit (γ_{0} + γ_{1} x_{i j} + γ_{2} y_{i j}),

M2: p_{r} (y_{i j}, x_{i j}; γ^{2}) = logit (γ_{0} + γ_{1} x_{i j} + γ_{2} y_{i j} + γ_{3} y_{i j}^{2}),

M3: p_{r} (y_{i j}, x_{i j}; γ^{3}) = logit (γ_{0} + γ_{1} x_{i j}),

Note that p_r(y_ij,x_ij;γ³) is nested in p_r(y_ij,x_ij;γ¹) and p_r(y_ij,x_ij;γ²), and p_r(y_ij,x_ij;γ¹) is nested in p_r(y_ij,x_ij;γ²). The response probability p_r(y_ij,x_ij;γ³) defines MAR nonresponse and hence, can be estimated by solving (3.3).

Estimate the unknown parameters θ₀,θ₁,θ₂ in (5.2) by SAS Proc NMIX, and then estimate γ by maximizing (3.4), as described in Section 3. The maximization was carried out by use of SAS Proc NLIN under the following 9 scenarios, as defined by the true response model and the assumed working response model:

S1: R1 set of respondents, M1 working response model.

S2: R1 set of respondents, M2 working response model.

S3: R1 set of respondents, M3 working response model.

S4: R2 set of respondents, M1 working response model.

S5: R2 set of respondents, M2 working response model.

S6: R2 set of respondents, M3 working response model.

S7: R3 set of respondents, M1 working response model.

S8: R3 set of respondents, M2 working response model.

S9: R3 set of respondents, M3 working response model.

Select the response model based on:

1. The Likelihood Ratio Test (LRT); test a saturated model [ $l (γ^{* *})$ ] against a nested model [ $l (γ^{*})$ ], assuming the $χ^{2}$ distribution under the null hypothesis H₀ that the nested model with a smaller number of parameters is correct. The test statistic is $λ_{L R T} = - 2 [l ({\hat{γ}}^{*}) - l ({\hat{γ}}^{* *})] ~ χ^{2}_{[\dim (γ^{* *}) - \dim (γ^{*})]}$ . Reject H₀ at the ɑ = .05 level.

2. AIC selection criterion: compare the values of the AIC as obtained for competing two models;

3. BIC selection criterion: compare the values of the BIC as obtained for competing two models.

Repeat the whole process independently 500 times.

5.2 Results

S1 Vs S2 (R1 – set of respondents, M1 – correct model, M2 – saturated model). Note that although M2 is a saturated model, it is also correct but with an additional term. The LRT selects the model M1 in 368 out of the 500 simulations. AIC selects M1 in 305 out of 500 simulations, BIC selects M1 in 324 simulations.

S1 Vs S3 (R1 – set of respondents, M1 – correct model, M3 – incorrect nested model). The LRT selects the correct model M1 in 500 out of the 500 simulations. AIC and BIC likewise select M1 in all the 500 simulations.

S4 Vs S5 (R2 – set of respondents, M2 – correct model, M1 – incorrect nested model). The LRT selects the correct model M2 in 433 out of the 500 simulations. AIC and BIC select M2 in 483 simulations.

S4 Vs S6 (R2 – set of respondents, M2 – correct model, M3 – incorrect nested model). The LRT selects the correct model M2 in 490 out of the 500 simulations. AIC and BIC select the correct model in all the simulations.

S7 Vs S8 (R3 – set of respondents, M3 – correct model, M1 – also correct but a saturated model). The LRT selects the model M3 in 241 out of 500 simulations. AIC selects M3 in 225 out of 500 simulations; BIC selects M3 in 420 simulations.

S7 Vs S9 (R3 – set of respondents, M3 – correct model, M2 – also correct but a saturated model). The LRT selects the M3 model in 241 out of 500 simulations. AIC selects M3 in 361 out of 500 simulations, BIC selects M3 in 450 simulations.

Note that when R3 is the set of respondents and M3 is the correct model, M1 and M2 also produce correct estimates of the response probabilities, although with additional estimated parameters. Thus, the fact that the LRT test and the AIC select the M3 model in about half of the simulations is not surprising. The use of the BIC criterion performs better in these cases.

The results are summarized in Table 1.

Table 1.

Percentages Out of 500 Simulations in which each of the Three Selection Procedures Selected the Correct Model, for Different Combinations of Correct (Rows) and Working (Columns) Response Probability Models.

	M1LRT AIC BIC	M2LRT AIC BIC	M3LRT AIC BIC
R1, M1 correct	– – –	73.6 61.0 64.8	100 100 100
R2, M2 correct	86.6 96.6 96.6	– – –	98 100 100
R3, M3 correct	48.2 45.0 82.0	48.2 72.2 90.0	– – –

Finally, we consider the case where a working model is incorrect but might be a good approximation of the correct model: let R1 be the set of respondents such that M1 is the correct working model. Let M4 be the following working model: $p_{r} (y_{i j}, x_{i j}; γ^{4}) = logit (γ_{0} + γ_{1} x_{i j} + γ_{2} y_{i j}^{2} + γ_{3} y_{i j}^{3})$ . Compare M4 with M1 (correct model). In this case, AIC selects the correct M1 model in 430 out of the 500 simulations and BIC selects the correct model in 431 simulations.

Sverchkov^[13] suggested testing whether the response is NMAR or MAR by testing the significance of the corresponding estimated coefficients in the saturated response model. We applied this idea by testing the significance of the estimated coefficients ${\hat{γ}}_{2}$ and ${\hat{γ}}_{3}$ under the response models $p_{r}^{(1)} (y_{i j}, x_{i j})$ and $p_{r}^{(2)} (y_{i j}, x_{i j})$ , (both assume NMAR nonresponse), when in fact the true response model is $p_{r}^{(3)} (y_{i j}, x_{i j})$ (MAR) or $p_{r}^{(1)} (y_{i j}, x_{i j})$ (NMAR), using the standard t-tests. (SAS Proc NLIN provides standard errors of the estimated coefficients.)

We considered two samples of respondents, R3 and R1. For R3, we found that when testing the working response model $p_{r}^{(1)} (y_{i j}, x_{i j})$ , in 432 out of the 500 simulations, ${\hat{γ}}_{2}$ was not significant at the 0.05 level. When testing the working response model $p_{r}^{(2)} (y_{i j}, x_{i j})$ , in 350 out of the 500 simulations, ${\hat{γ}}_{2}$ was not significant at the 0.05 level, and in 398 simulations ${\hat{γ}}_{3}$ was not significant. Recall that for the respondents’ sample R3, the working SAE model (5.2) for the observed outcomes is correct since the response is MAR.

For the respondents’ sample R1, the response model $p_{r}^{(1)} (y_{i j}, x_{i j})$ is correct and in all the 500 simulations, the estimator ${\hat{γ}}_{2}$ was found significant. Also, when testing the response model $p_{r}^{(2)} (y_{i j}, x_{i j})$ , the estimator was significant in 498 out of the 500 simulations. However, in 388 simulations the estimator ${\hat{γ}}_{3}$ was significant, even at the 0.01 level. This result might be explained by the fact that the working outcome model (5.2) is not correct when the response model is NMAR and thus, the likelihood (3.5), which conditions on the estimated random effects for the estimation of the γ coefficients is incorrect.

6. Summary

In this article we investigate the use of the likelihood function of the observed respondents’ data for selecting an appropriate response model under possible NMAR nonresponse. For estimating the hypothesized model, we applied the missing information principle. Despite of what seems to be a rather complex estimation process, we find in our simulation study that the AIC and BIC information criteria and the LRT test, when applicable, perform well for model selection. Clearly, the use of other likelihood-based tests and selection criteria should be investigated as well.

Footnotes

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The authors received no financial support for the research, authorship and/or publication of this article.

ORCID iD

Michael Sverchkov

References

Cepillini

, Siniscialco

and Smith

CAB.

The estimation of gene frequencies in a random mating population. Ann Hum Genet 1955; 20: 97–115.

Orchard

and Woodbury

MA.

A missing information principle: theory and application. Proc 6^th Berkeley Symp Math Stat Prob 1972; 1: 697–715.

Pfeffermann

New Important Developments in Small Area Estimation. Statistical Science 2013; 28: 40–68.

Pfeffermann

and Landsman

Are Private Schools Better than Public Schools? Appraisal for Ireland by Methods for Observational Studies. Ann App Stat 2011; 5: 1726–1751.

Pfeffermann

and Sikov

Imputation and estimation under nonignorable nonresponse in household surveys with missing covariate information. J Off Stat 2011; 27: 181–209.

Pfeffermann

and Sverchkov

Small-area estimation under informative probability sampling of areas and within selected areas. J Amer Stat Assoc 2007; 102: 1427–1439.

Pfeffermann

and Sverchkov

Inference under informative sampling. In: Pfeffermann

and Rao

, editors. Handbook of statistics 29B; Sample surveys: Inference and snalysis . North-Holland 2009, pp. 455–487.

Pfeffermann

and Sverchkov

Multivariate small area estimation under nonignorable nonresponse, Stat Theo Rel Fiel 2019; 3: 213–223.

Rao

JNK

and Molina

Small area estimation , 2nd ed, Wiley 2015.

10.

Riddles

, Kim

and Im

A propensity-score adjustment method for nonignorable nonresponse. J Surv Stat Meth 2016; 4: 215–245.

11.

Särndal

and Swensson

A general view of estimation for two phases of selection with applications to two-phase sampling and nonresponse. Int Stat Rev 1987; 55: 279–294.

12.

Sverchkov

A new approach to estimation of response probabilities when missing data are not missing at random. Joint Stat Meet, Proc Sect Surv Res Meth 2008; 867–874. URL: https://www.bls.gov/osmr/research-papers/2008/pdf/st080050.pdf

13.

Sverchkov

Is it MAR or NMAR? In; 2013 JSM Meetings, Proceedings of the Section on Survey Methods Research , 2013, pp. 2307–2311.

14.

Sverchkov

An algorithm for small area estimation under not missing at random non-response. Joint Stat Meet, Proc Sect Sur Res Meth 2022; 1735–1745.

15.

Sverchkov

and Pfeffermann

Prediction of finite population totals based on the sample distribution. Surv Meth 2004; 30: 79–92.

16.

Sverchkov

and Pfeffermann

Small area estimation under informative sampling and not missing at random non-response. J Royal Stat Soc, Series A 2018; 181: 981–1008.