1

RESPONSE MODEL SELECTION IN SMALL AREA ESTIMATION UNDER 

NOT MISSING AT RANDOM NONRESPONSE 

Michael Sverchkov

Bureau of Labor Statistics, Washington DC, USA

Danny Pfeffermann

Hebrew University of Jerusalem, Israel, and University of Southampton, UK.

ABSTRACT

Sverchkov and Pfeffermann (2018) 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 unobserved 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.

Key words: AIC, BIC information criteria, likelihood ratio test, population 

distribution, respondents’ model, sample distribution, 
2

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 nonresponse 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 (2009) 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 (2008) 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 (2008) 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 

(2018) 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.
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 

{ , ; 1,..., , 1,..., } ij ij i y i M j N x = =

represent the data in a finite population of 

N

units, comprised of

M

areas with 

Ni

units in area 

i , 

1

M

i

i

N N

=  = , where 

ij y

is 

the value of the outcome variable for unit 

j

in area 

i

and

,1 , ( ,..., ) ij ij ij K x

= x x

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: 

 

2

| , ~ ( | , ), 1,..., , 1,...,

~ ( ); ( ) 0, ( ) , U

U U

ij ij i ij ij i i

U U U U

i i i i u

y u f y u i M j N

u f u E u V u 

= =

= =

x x

 (2.1)

where 

U

i

u

is the i

th area level random effect. The target is to estimate the area 

means 

1

1

, 1,...,

Ni

i i ij

j

Y N y i 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 

Pr( ) i  = i s

; ii) select a sample 

i

s

of 

0 in 

units from selected area

i

with probabilities 

| Pr( | ) j i i  =   j s i s . Denote by

i

I , 

ij I

the sample indicators; 

1 iI =

if area 

i

is selected in the first stage and 0 

otherwise, 

1 ij I =

if unit 

j

of selected area 

i

is sampled in the second stage and 

0 ij I =

otherwise. Let 

1/ wi i=  ,

| | 1/ wj i j i

= 

denote the first- and second-stage 

sampling weights. 

In practice, not every unit in the sample responds. Define the response indicator;

1 Rij=

if unit 

i

j s 

responds and

0 Rij=

otherwise. The sample of respondents is 
4

thus

{( , ) : 1, 1, 1} R i j I I R = = = =

i ij ij

and the sample of nonrespondents among the 

sampled units is 

{( , ) : 1, 1, 0} c R i k I I R = = = =

i ik ik

. The response process is 

assumed to occur stochastically, independently between units. We assume 

1

0

i

n

ij j

R

=  

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, 1987). 

Define,

( | ) U U

i i i u u E u i s = − 

. Then, under the population model (2.1), the 

observed data follow the two-level ‘respondents’ model: 

 

( | , ) ( | , ,( , ) );

~ ( | ), ( | ) 0.

R ij ij i ij ij i

i i i

f y u f y u i j R

u f u i s E u i s

= 

  =

x x

 (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;

( | , ) ( | , ) U

R ij ij i ij ij i f y u f y u x x  .

Remark 1. The respondents’ model refers to the observed data and hence can be 

estimated and tested by standard SAE methods. See Pfeffermann (2013) and Rao 

and Molina (2015) for estimation and testing procedures in SAE, with references.

Let 

( , ) r ij ij p y x = Pr[ 1| , , , ] R y i s j s ij ij ij i =   x

. If the probabilities 

( , ) r ij ij p y x

were 

known, the sample of respondents could be considered as a two-stage sample 

from the finite population with known sampling probabilities 

 i

and 

| | ( , ) j i j i r ij ij   = p y x . In this case, the area means 

Yi

can be estimated as in 

Pfeffermann and Sverchkov (2007). 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, 2004);

( | , ,( , ) )c

ij ij i f y u i j R x  =

1

1

[ ( , ) 1] ( | , ,( , ) )

{[ ( , ) 1]| , ,( , ) }

r ij ij ij ij i

r ij ij ij i

p y f y u i j R

E p y u i j R

−

−

− 

− 

x x

x x

. (2.3)
5

See Sverchkov and Pfeffermann (2018), and Pfeffermann and Sverchkov (2019)

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[ 1| , , , ; ] R y i s j s ij ij ij i =   x γ

( , ; ) r ij ij

= p y x γ

, where 

γ

is a vector of unknown coefficients. We assume that 

( , ; ) r ij ij p y x γ

is differentiable with respect to 

γ

and satisfies the same mild

regularity conditions as in Sverchkov and Pfeffermann (2018). 

Under these assumptions, if the missing outcome values were observed, γcould 

be estimated by solving the likelihood equations: 

 

( , )

log ( , ; ) r ij ij

i j R

p y





+





x γ

γ ( ,k)

log[1 ( , ; )] 0

c

r ik ik

i R

p y



 −

=





x γ

γ

. (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. 1955, Orchard and Woodbury, 

1972): Let 

{ , ,( , ) ; , 1,..., , t 1,..., } O y n i j R h M N =  = = ij i ht i x

represent the known 

observed data used below. Since no observations are available for 

( , ) ci k R  , 

solve instead,
6

( , ) ( ,k)

log ( , ; ) log[1 ( , ; )]

c

r ij ij r ik ik

U

i j R i R

p y p y E O

 

        −     +

            x γ x γ

γ γ

(2.3)

( , )

log ( , ; ) by

r ij ij

i j R

p y





=





x γ

γ

 

1

1

( , )

log[1 ( , ; )] [ ( , ; ) 1] | , ,( , )

0.

c {[ ( , ; ) 1] | , ,( , ) }

r ik ik

re r ik ik ik i

s

i k R re r ik ik ik i

p y E p y u i k R

E O

E p y u i k R

−

−



     −     − 

     + =   −       



x γ

x γ x

γ

x γ x

 

 (3.2) 

See Sverchkov (2008) and Sverchkov and Pfeffermann (2018) for derivation of 

(3.2). In these equations,

EU s re ,E ,E

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 

( , ; ) r ij ij p y x γ

depend on only 

ij x , 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,

 

( , )

log ( ; ) r ij

i j R

p





+





x γ

γ ( , )

log[1 ( ; )] 0,

c

r ik

i k R

p



 −

=





x γ

γ

 (3.3) 

where 

( ; ) Pr(R 1| ; ) r ij ij ij p x γ x γ = = .

Sverchkov and Pfeffermann (2018) propose to solve the equations (3.2) by 

maximizing the log-likelihood leading to them, i.e., maximizing,
7

( , )

( ) log ( , ; ) r ij ij

i j R

l p y



γ x γ =  

1 *

1 *

( , )

{[ ( , ; ) 1]log[1 ( , ; )]| , ,( , ) }

c {[ ( , ; ) 1]| , ,( , ) }

re r ik ik r ik ik ik i

s

i k R re r ik ik ik i

E p y p y u i k R E O

E p y u i k R

−

−



  − − 

+       − 



x γ x γ x

x γ x

. (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 

i u by 

ˆ

i u , 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,

( 1)

( , )

log ( , ; ) q

r ij ij

i j R

p y

+



 x γ

1 ( ) ( 1)

1 ( )

( , )

{[ ( , ; ) 1]log[1 ( , ; )]| , ,( , ) } ˆ

c {[ ( , ; ) 1]| , ,( , ) } ˆ

q q

re r ik ik r ik ik ik i

q

i k R re r ik ik ik i

E p y p y u i k R

E p y u i k R

− +

−



− − 

+

− 



x γ x γ x

x γ x

 

 (3.5) 

with respect to 

( 1) q+

γ . The maximization can be carried out, for example, by SAS 

Proc NLIN. See Sverchkov (2022) 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. (2016) 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 

1 2 x x x = ( , )

with 

2

dim( ) 1 x 

, such that

Pr( 1| , ) R y ij ij ij= x = = Pr( 1| , ) R y ij ij ij 1

x

. In other words, the covariates in 

2

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 
8

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 

2

x

. (Other covariates that are predictive of the outcome but not of the 

response might exist as well). 

Clearly, the larger is 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 (2018), 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. 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 are means remain virtually unbiased in each of 

the areas. 

Riddles et al. (2016) prove asymptotic normality of the estimate

γ

ˆ

under general 

regularity conditions. 

Example 1. (Sverchkov and Pfeffermann 2018): 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,

 ( , ) y ij i p x u = Pr( 1| , ,( , ) ; ) ij ij i y x u i j R =  β

0 1

0 1

exp( )

1 exp( )

ij i

ij i

x u

x u

 

 

+ +

=

+ + +

, 

i.i.d.

2

~ (0, ) i u u N  . 

Consider a generic response model, 

( , ; ) r ij ij p y x γ = Pr[ 1| , , , ; ] R y x i s j s ij ij ij i =   γ . 

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

1

log[1 ( , ; )]

[ ( , ; ) 1] , ,( , ) r ij ij

re r ij ij ij i

p y x

E p y x x u i j R −

  −      −  =

    

γ

γ

γ
9

1

log[1 (1, ; )]

( , )[ (1, ; ) 1] r ij

y ij i r ij

p x

p x u p x



−

 −

− +



γ

γ

1

log[1 (0, ; )]

[1 ( , )][ (0, ; ) 1] r ij

y ij i r ij

p x

p x u p x

−

 −

− −



γ

γ

γ

; (3.6)

1

{[ ( , ; ) 1]| , ,( , ) } E p y x x u i j R re r ij ij ij i

−

γ −  = 1( , )[ (1, ; ) 1] y ij i r ij p x u p x

−

γ − + 1

[1 ( , )][ (0, ; ) 1] y ij i r ij p x u p x

−

− − γ . (3.7)

The random effects 

ui

and the logistic probabilities 

( , u ) y ij i p x

can be estimated by 

use of the SAS procedure PROC NLMIX. 

Example 2. (Sverchkov 2022): 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, ˆ

ui

, and

for each 

( , ) i j R  , estimates of 

( )

1

( , ) ( | , ,( , ) ) ˆ ˆ

l

y ij i R l ij l ij i p u P a y a u i j R =    +

x x , 

l L = + 0,..., 2

;

0

a = − , 

L 2

a +=  , 

max( ) min( )

min( ) ( 1) ij ij

l ij

y y

a y l

L

−

= + − , 

l L = + 1,..., 1. The max and min are over all the observed values 

ij y .

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

1

log[1 ( , ; )]

[ ( , ; ) 1] , ,( , ) r ij ij

re r ij ij ij i

p y

E p y u i j R −

  −      −  

    

x γ

x γ x

γ

1

( ) 1

1

log[1 ( , ; )]

ˆ ˆ ( , )[ ( , ; ) 1]

L

l r l ij

y ij i r l ij

l

p a

p u p a

+

−

=

 −

−





x γ

x x γ

γ

, (3.8) 

 

1

[ ( , ; ) 1] , ,( , ) E p y u i j R re r ij ij ij i

−

x γ x −  

1

( ) 1

1

ˆ ˆ ( , )[ ( , ; ) 1]

L

l

y ij i r l ij

l

p u p a

+

−

=

 x x γ − , (3.9)

where 

( , ; ) Pr[ 1| , , , ; ] r l ij ij ij l ij i p a R y a i s j s x γ x γ = = =   .

Substitute (3.8) and (3.9) into (3.2) and estimate 

γ

by iteratively maximizing (3.5).
10

4. SELECTION OF A 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, e.g., Pfeffermann and 

Landsman (2011) and Pfeffermann and Sikov (2011). 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. 1- Compare different models based on information criteria

such as the Akaike information criterion, AIC= 2 ( ) 2dim( ) − + l γ γ , or Schwarz 

information criterion, 

BIC= 2 ( ) dim( )log( ), i

i s

l n n n



− + = γ γ 

; 2- 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),
11

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 

, 1,...,100, 1,...,20 ij x i j = =

from a Uniform(0,2)

distribution. Next, generate sample values from the small area model,

 | , ~ ( ,1), 1,...,100, 1,...,20 ij ij i ij i y x u N x u i j + = =;

~ (0,1) iu N . (5.1) 

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

(1) ( , ) logit( / 2 2 ) r ij ij ij ij p y x x y = − + ,

(2) 2 ( , ) logit( / 2 2 0.3 ) r ij ij ij ij ij p y x x y y = − + − ,

(3) ( , ) logit(1.5 ) r ij ij ij p y x x = .

Select 3 sets of respondents: 

R1 uses Poisson sampling, independently between the units with response 

probabilities

(1) ( , ) r ij ij p y x ,

R2 is the same as R1 but with response probabilities 

(2) ( , ) r ij ij p y x ,

R3 is the same as R1 but with response probabilities 

(3) ( , ) r ij ij p y x .

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,

 

0 1 2 | , ~ ( , ), 1...100, 1...20 ij ij i ij i y x u N x u i j    + + = =

; 

2

~ (0, ) i u u N  . (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: 

1

0 1 2 ( , ; ) logit( ) r ij ij ij ij p y x x y γ = + +    ,

M2: 

2 2

0 1 2 3 ( , ; ) logit( ) r ij ij ij ij ij p y x x y y γ = + + +     ,
12

M3: 

3

0 1 ( , ; ) logit( ) r ij ij ij p y x x γ = +   ,

Note that

3

( , ; ) r ij ij p y x γ

is nested in 

1

( , ; ) r ij ij p y x γ

and

2

( , ; ) r ij ij p y x γ , and

1

( , ; ) r ij ij p y x γ

is nested in 

2

( , ; ) r ij ij p y x γ . The response probability 

3

( , ; ) r ij ij p y x γ

defines MAR nonresponse and hence, can be estimated by solving (3.3). 

Estimate the unknown parameters

0 1 2    , ,

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 

H0

that the 

nested model with a smaller number of parameters is correct. The test statistic is

** *

* ** 2

[dim( ) dim( )]

2[ ( ) ( )] ~ ˆ ˆ

LRT   l l

−

= − −

γ γ

γ γ . Reject 

H0

at the 

 = .05

level.

2 - AIC selection criterion: compare the values of the AIC as obtained for the 

corresponding two models; 

3 - BIC selection criterion: compare the values of the BIC as obtained for the 

corresponding two models.
13

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 so far are summarized in table 1.
14

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.

M1

LRT AIC BIC

M2

LRT AIC BIC

M3

LRT 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 82 48.2 72.2 90 --- --- ---

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:

4 2 3

0 1 2 3 ( , ; ) logit( ) r ij ij ij ij ij p y x x y y γ = + + +     . 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 (2013) 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 

2



ˆ

and

3



ˆ

under the response models 

(1) ( , ) r ij ij p y x

and 

(2) ( , ) r ij ij p y x ,

(both assume NMAR nonresponse), when in fact the true response model is

(3) ( , ) r ij ij p y x

(MAR) or 

(1) ( , ) r ij ij p y x

(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 

(1) ( , ) r ij ij p y x , in 432 out of the 500 

simulations,

2



ˆ

was not significant at the 0.05 level. When testing the working 

response model 

(2) ( , ) r ij ij p y x , in 350 out of the 500 simulations, 

2



ˆ

was not 
15

significant at the 0.05 level, and in 398 simulations 

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 

(1) ( , ) r ij ij p y x

is correct and in 

all the 500 simulations, the estimator 

2



ˆ

was found significant. However, when 

testing the response model

(2) ( , ) r ij ij p y x , in 388 simulations the estimator

3



ˆ

was 

significant, even at the 0.01 level, and

2



ˆ

was significant in 498 simulations. 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 paper 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. 

REFERENCES

Cepillini, R., Siniscialco, M., and Smith, C.A.B. (1955). The estimation of gene 

frequencies in a random mating population. Annals of Human Genetics, 20, 97-

115.

Orchard, T., and Woodbury, M.A. (1972). A missing information principle: theory 

and application. Proceedings of the 6th Berkeley Symposium on Mathematical 

Statistics and Probability, 1, 697-715.
16

Pfeffermann, D. (2013). New Important Developments in Small Area Estimation. 

Statistical Science, 28, 40-68.

Pfeffermann, D., and V. Landsman (2011), “Are Private Schools Better than Public 

Schools? Appraisal for Ireland by Methods for Observational Studies,” Annals of 

Applied Statistics, 5, 1726–1751.

Pfeffermann, D. and Sikov N. (2011). Imputation and estimation under 

nonignorable nonresponse in household surveys with missing covariate 

information. Journal of Official Statistics, 27, 181–209.

Pfeffermann, D., and Sverchkov, M. (2007). Small-Area Estimation under 

Informative Probability Sampling of Areas and Within Selected Areas. Journal of 

the American Statistical Association, 102, 1427-1439.

Pfeffermann, D., and Sverchkov, M. (2009). Inference under Informative Sampling. 

In: Handbook of Statistics 29B; Sample Surveys: Inference and Analysis. Eds. D. 

Pfeffermann and C.R. Rao. North Holland. pp. 455-487.

Pfeffermann, D. and Sverchkov, M. (2019). Multivariate small area estimation 

under nonignorable nonresponse, Statistical Theory and Related Fields, 3, pp. 

213-223.

Rao, J.N.K., and Molina, I. (2015), Small Area Estimation, 2nd Edition, Wiley.

Riddles, K.M., Kim, J.K. and Im, J. (2016). A propensity-score adjustment method 

for nonignorable nonresponse. Journal of Survey Statistics and Methodology, 4, 

215-245.

Särndal, C.E. and Swensson B. (1987). A general view of estimation for two 

phases of selection with applications to two-phase sampling and nonresponse. 

International Statistical Review, 55, 279-294.

Sverchkov, M. (2008). A new approach to estimation of response probabilities 

when missing data are not missing at random. Joint Statistical Meetings,

Proceedings of the Section on Survey Research Methods, 867-874.
17

Sverchkov, M. (2013). Is it MAR or NMAR? 2013 JSM Meetings, Proceedings of 

the Section on Survey Methods Research, pp. 2307-2311

Sverchkov, M. (2022). An Algorithm for Small Area Estimation under Not Missing 

At Random Non-response. Joint Statistical Meetings, Proceedings of the Section 

on Survey Research Methods, pp. 1735-1745.

Sverchkov, M., and Pfeffermann, D. (2004). Prediction of Finite Population Totals 

Based on the Sample Distribution. Survey Methodology, 30, 79-92.

Sverchkov, M. and Pfeffermann, D. (2018). Small area estimation under 

informative sampling and not missing at random non-response. Journal of Royal 

Statistical Society, ser. A, 181, Part 4, pp. 981–1008.