Small Area Modeling for Survey Data with Smoothed Error Covariance 

Structure via Generalized Design Effects 

A.C. Singh1, R.E. Folsom, Jr2., and A.K. Vaish2 

Statistics Canada1 and RTI International2 

A.C. Singh, 16-Q, R.H. Coats, 120 Parkdale Ave, Ottawa, ON K1A 0T6 avi.singh@statcan.ca

Abstract 

We consider the problem of specifying the design-based error covariance structure in small area modeling with survey data 

because of too many unknown design parameters. As a compromise, it is customary to treat the estimated covariance as 

known although undesirable because of its instability, i.e., relative variance of this estimate could be high, and as a result, it 

could cause serious underestimation of variance of estimates. To alleviate this problem, one can either model the estimated 

error covariance structure in addition to the main task of modeling the estimated small area estimates (SAEs), or smooth the 

estimated covariance by only specifying its mean function. We advocate smoothing over modeling the error covariance 

because of strong simplifying assumptions needed in modeling the error covariance structure of estimated covariance of 

SAEs, and because practitioners, in general, prefer to take the least assumptions in modeling path philospphy. In practice, 

smoothing may work quite well as supported by the recent work of Singh, You, and Mantel (2005) who used the well known 

property of constant deff (design effect) over suitable subgroups of estimates to smooth the error covariance in modeling 

direct estimates from the cross-sectional and longitudinal data obtained from the Canadian Labour Force Survey. In this 

application, however, cross-sectionally the small areas were treated as strata, i.e., cross-sectionally, the error covariance was 

diagonal. This is rather restrictive because small areas or domains need not be strata. Moreover, using optimal estimating 

functions (EFs) of Godambe and Thompson (1986) for estimating model parameters, it can be shown that even when small 

areas are strata, additional summary statistics (in particular, the estimated domain count if it is not constant) for each domain 

or area should be used besides the direct SAEs, and thus the error covariance becomes necessarily nondiagonal. Use of EFs 

was also advocated by Singh, Folsom, and Vaish (2002, 2003) to define approximate EF-based Gaussian likelihood (EFGL) 

for obtaining efficient SAEs under a HB (hierarchical Bayes) framework. In this paper, we propose the idea of using g-deff 

(generalized design effect), defined earlier by Rao and Scott (1981) in the context of categorical data analysis, to deal with 

the case of nondiagonal error covariance. Simulation results for SAEs based on a linear mixed model show that the EFGL 

method with the proposed smoothing works quite well in general, and provides improved coverage of confidence intervals in 

particular.

Key Words: Estimating Functions; Generalized Deffs; Ignorable and nonignorable designs; unstable estimated error 

covarianc 
1. Introduction 

In any modeling problem, one needs to specify the mean function and the error covariance structure under a semi-parametric 

approach based on first two moments. In the case of modeling direct domain total estimates { } of small areas 

(see e.g., Fay and Herriot, 1979) from surveys, special problems arise in specifying the error covariance structure because 

it depends on too many unknown quantities (e.g., the first and second order inclusion probabilities for the nonsampled units) 

that are not estimable. However, since an approximately unbiased estimate of may be available, it is customary (see 

Rao, 2003, pp. 76) for such models to treat as known. For large samples, such an assumption is, of course, commonly 

made. However, for small samples , this is clearly not desirable because like the direct point estimates , the direct 

variance estimates are also subject to instability, and so treating as known may cause serious underestimation of 

variance of SAEs. As an alternative, one may want to model as well in addition to which renders the SAE problem 

even more compounded as variance of now involves unknown third and fourth order sample inclusion probabilities. 

Under certain simplifying assumptions about the design, and the superpopulation model, the problem can, however, be 

simplified, see e.g., the use of Wishart distribution for by Otto and Bell (1995). Nevertheless, in practice, it would be 

desirable to make rather weak assumptions about that seem plausible, because practitioners, in general, prefer to take the 

path of least assumptions in modeling. 

,

: 1,..., y d td D =

Vt

ˆVt Vt

ˆVt

d

n y d,

t

( )

ˆVt d ( )ˆVt d

ˆVt y d,

t

ˆVt

ˆVt

ˆVt

In trying to model under weak assumptions, it may be useful to first observe that the problem of modeling the mean of 

is lot simpler than that of the error covariance structure, i.e., modeling the variance of . If the goal is to find a stable or 

smoothed estimate of without incurring too much bias, then one can simply specify the mean function for and leave 

the covariance function unspecified, and then estimate parameters of the mean function by analogy with ordinary least 

squares. The smoothed estimate is treated as known or may be allowed to depend on mean parameters of the model for 

. It is advocated in this paper that smoothing might provide a reasonable practical compromise between the two extremes 

of modeling or no modeling of . 

( )

ˆVt d

ˆVt

ˆVt

Vt

 Vt

ˆVt

Vt



y d,

t

ˆVt

In this paper, we propose the use of generalized design effect (g-deff) considered earlier by Rao and Scott (1981) in the 

context of categorical analysis of survey data. G-deffs are defined as the eigen-values ( λi

’s) of the matrix product 

where is the error covariance under a suitable working assumption such as that of simple random sampling design 

or that the design is ignorable for the superpopulation model under consideration. In the case of diagonal (i.e., when the 

small areas can be treated as strata), Singh, You and Mantel (2005) used the well-known property of constant deff 

(

* 1 Vt− Vt* Vt

ˆVt

* V V td td () () ) over suitable subgroups of direct estimates to smooth by assuming that is stable. More specifically, in the ˆVt* ˆVt
context of the Canadian Labour Force Survey with partial overlapping samples over six months for each panel, they 

considered the combined longitudinal and cross-sectional data set of direct small area estimates for which the matrix now 

becomes block diagonal where each block is a 6x6 matrix corresponding to each area. Now under the weak and rather 

plausible assumption of constant deff over six months for each area, and common fixed-lag correlations over areas and time, 

they got improved results over the alternative smoothing methods described in You, Rao, and Gambino (1999). However, the 

above method has some limitations: first, the nonnegativity of the smoothing matrix is not ensured when fixed lag 

correlations are replaced by averages over time and areas; second, for the more general case of areas not necessarily being 

strata as in the case of demographic domains, the matrix becomes non-block-diagonal and it is not clear how the simple 

assumption of constant correlations can be applied to smooth the off-diagonal terms. In all these cases, however, the proposed 

method of g-deff can be used wherein the eigen-values are averaged over suitable subgroups to reduce their instability. 

ˆVt

ˆVt

For the problem of SAE modeling with direct estimates { ,

: 1,..., y d td D = }, one can think of an implicit underlying 

superpopulation model at an aggregate level for the outcome variable for each unit k in domain ; here the aggregate 

level model signifies that all the model covariates are at the area or domain ( d ) level. Conditional on the covariates, the 

model under the joint randomization of design (

y d

π ) and superpopulation model (ξ , this is conditional on any random effects 

in the model) can be expressed for d D =1,..., as 

 

, , ,, ,

1

,, ,

,

yd yd yd yd d yd

y d x d d d dk x d d y d k

t T e T NA

AA N A β η ,

ε βη μ

−

=+ =

= ++ ≈ += ′ ′ ∑

 (1.1) 

 

2

, ( )

2 1

~ (0, ); ~ (0, )

~ (0, ), 0

yd t d d

dk d dk k

e NV N

N

η

ε

η σ

εσ ε

− ∑ ≈

 (1.2) 

where is the finite population total of for area ; is the known total count for area d which we will also denote as 

, the total for the count variable taking value 1 or 0 depending on whether the unit is in domain or not; is of course 

the whole population count;

Ty d,

y d Nd

Tc d,

c d N

Ax,d

is the vector of domain- d averages for the p -vector of auxiliary variables x ; β is a 

p -vector of fixed parameters, and ηd

is the random effect corresponding to domain assumed to be normal with mean 

and variance 

d

0

2 σ η

; is the sampling or observation error in the estimated domain total and assumed to be normal with 

mean 0 and variance , and independent of 

y d,

e

Vt d( ) ηd

 but may or may not be independent over d , and dk ε is the superpopulation 

model error whose contribution to the small area total or average is negligible compared to other terms because is 

assumed to be large. The model error 

Nd

dk ε has mean 0 and variance 2 σ ε

but is not necessarily assumed to be normal although 

is assumed to be approximately normal due to central limit theorem provided is sufficiently large. Note that y d,

t

d

n μ y d,

is the 

superpopulation mean parameter of y given Ax,d

and is the limit of the corresponding finite population parameter as 

. 

Ay d,

Nd → ∞
Now, the problem of smoothing gets more complicated if one considers the above superpopulation model explicitly and 

the issue of optimal sample summary statistics to be used for estimating model parameters. In fact, the usual choice of 

summary statistics { } is somewhat ad hoc. Using Godambe and Thompson (1986), it can be shown (see next 

section) using census estimating functions that for aggregate level modeling, the optimal set of summary statistics is 

{

ˆVt

,

: 1,..., y d td D =

, ,

( , ) : 1,..., d y d c d ψ = = tt d D }. Here it is assumed that is not constant (i.e., not equal to ) as is typically the case in 

practice, else the usual set of summary statistics {

c d,

t Nd

,

: 1,..., y d td D = } becomes optimal. Thus the error covariance for each 

becomes a 2x2 matrix and the resulting overall matrix 

( )

ˆVψ d

d Vˆ

ψ necessarily nondiagonal requiring a more general smoothing 

method such as the proposed one. 

Further, in the case of unit-level superpopulation model, the problem of smoothing Vˆ

ψ becomess considerably more involved 

because of additional summary statistics that enter in the picture from the theory of EFs. Note that Singh, Folsom, and Vaish 

(2002, 2003) proposed estimating function-based Gaussian likelihood (EFGL) methodology in a hierarchical Bayes 

framework to generalize the usual SAE modeling of Fay and Herriot (1979) at the aggregate level to unit level modeling in 

the interest of efficiency gains obtained from using more detailed information. However, they didn’t specify the summary 

statistics explicitly. In this paper, the vector ψ of optimal summary statistics is specified explicitly which makes it 

computationally easier to use g-deff for smoothing the error covariance structure Vˆ

ψ

 under unit-level modeling. 

In Section 2, we consider first the general problem of modeling with survey data, and then the choice of appropriate finite 

population parameters in the context of SAE modeling and the corresponding sample summary statistics. We consider both 

aggregate and unit level modeling and identify the problem of instability of Vˆ

ψ

in each case for the corresponding vector ψ . 

In the next section 3, we review the application of usual deff in the generalized variance function modeling, and its use for 

smoothing when it is diagonal. We next describe the proposed method of g-deff for smoothing nondiagonal or ˆVt

ˆVt Vˆ

ψ

. In 

Section 4, we present a simulation study for the linear mixed superpopulation model and show how EFGL-smothed compares 

with EFGL-unsmoothed. In the simulation study, the case of aggregate level modeling was, however, not considered 

although it follows readily from the general case by replacing the unit-level covariate value dk x with the domain level average 

value Ax,d

for each unit k d in domain . Finally, Section 5 contains summary and some concluding remarks. 

2. Modeling with Survey Data 

The difficulty in modeling survey data is well known in view of the two basic results of Godambe (1955, 1966): first, the 

nonexistence of uniformly minimum variance unbiased estimate in a suitable linear class which causes difficulty in using a 

semiparametric approach, and second, likelihood being flat for the unseen (i.e., nonselected population units) given the seen 

and thus making it difficult to use the likelihood approach either in a frequentist or a Bayesian framework. The main reason 
underlying these problems that distinguish survey modeling from mainstream statistics is that there are too many finite 

population parameters to cope with if one identifies each unit’s characteristic as a parameter of interest. The reality is that in 

practice we are not interested in characteristics at the unit level, but instead we need to define suitable finite population 

quantities corresponding to a group of units. For the difficult but realistic problem of nonignorable designs for a given 

superpopulation model, this can be done using census EFs (i.e., assuming the sample is the census, see, e.g., Binder, 1983) 

which depend, of course, on the model parameters to be estimated. Alternatively, Pfeffermann and Sverchkov (1999, 2003) 

proposed an elegant use of weighted likelihood approach for survey data. However, they require modeling of sampling 

weights which may not be desirable for practitioners who prefer the route of least assumptions in modeling as mentioned 

earlier. Using the idea of EFs, Godambe and Thompson (1986) showed how optimal EFs can be obtained under joint πξ -

randomization. Thus for the aggregate level superpopulation model conditional on covariates Ax,d

and random effects ηd

 dk x d d dk ,

y A = ++ ′ β η ε (2.1)

 the optimal EFs for (, ) η β are 

 

() , , , ,

,, , , , ,

( )( )

() (

d dk x d d dk y d x d d c d k

,

)

x d dk x d d dk x d y d x d x d d c d d k d d

yA w t A t

A y A w At A A t

η

β

φ β η β η

φ β η

= −− = ′ ′ − +

= − ′ ′ − = − +

∑

∑∑ ∑ ∑ β η

 (2.2) 

where are the sampling weights, and for each domain , varies from 1 to . It follows that the summary statistics 

{

wdk d k dn

, ,

( , ) : 1,..., d y d c d ψ = = tt d D ′ } where , corresponding to the finite population quantities 

{ } form an optimal set for estimating the parameters( ,

, ,

ˆ

,

d d n n

y d dk dk c d dk d k k

t y = = ∑ ∑ w t w = N

, ,

( , ) : 1,..., TT d D yd cd = η β ) . This is an interesting finding which came 

about by working with EFs under the superpopulation model. In contrast, with the traditional area-level model, one starts 

directly with as an estimate of the parameter of interest y d,

t

,, ,

( T NA N A yd d yd d xd d = ≈ + ′ β η ) , in which case it is easy to 

overlook the utility of the other summary statistic because the domain-level covariates and c d,

t Nd Ax,d

are assumed to be 

known. 

The aggregate-level model (1.1) can be generalized for the enlarged summary statistics ψ d

by introducing unknown 

incidental parameters defined as the limit of d

p N N d

as N → ∞ . The model with all the parameters is given by 

 (2.3)

, , ,

, , ,

yd d yd yd

cd d cd cd d cd

t Np e

t Np eNp e

μ

μ

≈ +

≈ + = + ,

where 1 ,, ,

, μ yd xd d cd =+ = A′ βημ . The sampling errors for in equations (1.1) and (2.3) do not change because y d,

t Np N d d ≈ . 

Observe that the introduction of incidental parameters (although not of direct interest) is bit of an artifact but has the 

desired effect of engaging in estimating the parameters 

d

p

c d,

t (, ) β ηd

as implied by EFs. However, unlike EFs, in the model 

(2.3), the data is separable from the parameters in that the mean function on the right hand side consists of only nonrandom 

covariates and unknown parameters. As a result, it becomes easier to apply the proposed method of g-deff for smoothing 

Vˆ

ψ unlike Vˆφ

corresponding to (2.2). Note that Vˆ

φ

involves the mean parameters and so smoothing would be required at each 
cycle of MCMC in using EFGL. With model (2.3), the smoothing of Vˆ

ψ

is done only once before starting the cycles of 

MCMC. However, unlike EFGL with φ , here the model’s mean function is nonlinear in the enlarged parameter set 

{ , , , 1,..., } d d β η pd D = which would be computationally more demanding than the linear case. However, in the MCMC 

cycles for finding conditional posteriors, model mean turns out to be linear in (, ) β ηd

given and so no new calculations 

are required except that there would be more steps in each MCMC cycle because of incidental parameters. 

d

p

Next for unit-level modeling, (2.3) can be modified as follows: 

 

( ) , , ,

,

( )

( )

d dk dk d dk y d x d c d d k

dk dk dk d dk xy xx x d d d k d

yx w t t t

xy x w t t t

η

β

φ β η β η

φ β η ′β η

= − ′ ′ − = − −

= − ′ ′ − = − +

∑

∑∑ ∑ (2.4) 

The optimal summary statistics ψ for estimating (, ) β ηd

are obtained as 

,,,

{ , , ,,} y d x d c d xy xx t t t tt

′ with somewhat self-explanatory 

new notations for certain estimated population totals. In this case, the unit-level model for this new set of summary statistics 

can be expressed with the new set of incidental parameters 

,

{, , } d x d xx p μ μ ′

as : 

 

, , ,

, ,

, , ,

, ,

( )

yd d yd yd

cd d cd

xd d xd xd

xy xx d x d d x d

xx xx xx

t Np e

t Np e

t Np e

tN p e

tN e

μ

μ

μβ μη

μ

′

′ ′′

≈ +

≈ +

≈ +

≈+ +

= +

∑ y d

 (2.5)

where, as defined earlier, μ y d,

is the limit of as and Ay d, Nd → ∞ μ x,dis similarly defined. As in the case of aggregate-level, 

although the model (2.5) is now nonlinear because of incidental parameters, model mean turns out to be linear for computing 

conditional posteriors, and so no new computational complexity is involved except that there are more steps in MCMC cycles 

because of extra parameters. It may be noted that for unit-level modeling, the covariates dk x ’s are typically taken as 

categorical (such as demographic group indicators) because the domain totals Tx,d

or averages Ax,drequired for each variable 

may only be available in practice for domain counts or proportions for each covariate category. 

Often, in practice, one may be interested in nonlinear modeling. Although it seems difficult to define optimal summary 

statistics for general nonlinear models at the unit level, it is indeed possible with categorical x 's . Suppose, l denotes the 

cell obtained by completely cross-classifying all the th l x -categories, and denotes the known population count for -

subdomain. It can be shown from EFs that the optimal summary statistics 

Ndl dl

ψ now are but there are too many 

or because the total number of cross-classified l -categories can be very large. While it may be tempting to use all 

these summary statistics by introducing a large number of corresponding incidental parameters , it is best to avoid 

modeling with very unstable summary statistics. A compromise may be to collapse over to reduce the set of summary 

statistics. In doing so, we will only need to introduce unknown incidental parameters as before, but we will also need to 

assume that the values of relative to are known. This is a rather innocuous assumption because 

, ,

{ , ,} y d c dl xy ttt

c dl ,

t

ˆNdl

dl p

c dl ,

t l

d

p

dl p d

p l d p (

dl d = p p ) can 
be well approximated by the corresponding known finite population quantity N N dl d . Now that the modeling problem is 

formulated in terms of the summary statistics ψ , we can turn to the question of smoothing Vˆ

ψ before model fitting using the 

EFGL method mentioned earlier. 

3. Proposed Method of Smoothing Vˆ

ψ using g-deff

First we review the use of Deff for generalized variance function (GVF) modeling. 

3.1 Deff for GVF modeling 

Suppose Vψ

is a D×1 diagonal matrix, and let dγ be the design effect for the domain-d SAE, i.e., d d V V ψ ψ () () γ

∗

= d

. In view 

of the fact that the design effects d

γ are often approximately constant over a suitable set of statistics, a simple type of GVF 

modeling assumes that the mean of is proportional to where the covariate is taken to be approximately 

known. As mentioned earlier, is estimated under the assumption of simple random sampling or that the design is 

ignorable for the model and that the estimate is assumed to be stable. So for the GVF model with , 

( )

ˆVψ d Vψ ( ) d

∗ Vψ ( ) d∗

Vψ ( ) d

∗

( )

ˆVψ d

∗

d D =1,...,

 (3.1)

() () ,

ˆV Ve ψ ψ d d γ

∗

= + v d

where the error covariance structure is not specified, the mean parameter γ can be estimated by the average of ˆ

d

γ over 

domains where () ()

ˆ ˆ ˆ

d d V V ψ ψ

γ

∗

= d d. The smoothed estimate of is then obtained as . If it does not seem 

reasonable to assume that the deff is constant over all areas, then one can compute separate estimates of 

Vψ ( ) d () ()

ˆ ˆ V V ψ ψ d

γ

∗  =

γ for subgroups of 

areas for which the constant deff assumption seems plausible. Note that may depend on the mean parameters Vψ ( ) d

∗

{, } β ηdof the small area model as in the case of binary data, and then the smoothed variance estimate can also be allowed to 

depend on unknown mean parameters. When dealing with discrete data, this is clearly a desirable feature. 

3.2 GVF-type Modeling for Non-diagonal Vˆ

ψ

The above GVF modeling to smooth Vˆ

ψ

is not applicable when Vˆ

ψ

is nondiagonal because the concept of deff is not defined 

for off-diagonal terms of covariances. In this case, using the concept of g-deff (see e.g., Rao and Scott, 1981) defined as 

eigen-values ( , 1,..., k λ k K = ) of * 1 Vψ

−

Vψ

, K being the dimension of vector ψ , we can write a GVF-type model as 

 (3.2)

ˆ

( ) ( ) () j k k j k

vec V vec q q vec e ψ = ∑ ∑ λ ′ + v

where ‘vec’ notation is used to signify that columns of the matrix are stacked one above the other, the second sum is over the 

th j subgroup of summary statitics, and the first sum is over all subgroups. The above model is motivated by the matrix result 

for a pair of real symmetric matrices with at least one of them being positive definite (cf: C.R. Rao, 1973, pp.41) which states 

that there exists a nonsingular matrix Q such that (here the pair of matrices are Vψ

and Vψ

∗

), 
 

K

kkk k

K

k k k

V Q Q qq

V QQ q q

ψ

ψ

λ

∗

=Λ =′ ′

= =′ ′

∑

∑

 (3.3)

where ( ) k Λ = diag λ , and is the column of . It follows from (3.3) that if k

q

th k Q Vˆ

ψ

∗

as well as for each selected 

subgroup of areas are stable, then the instability of 

Kj

k k k ∑ q q′

Vˆ

ψ

can be overcome by smoothing the estimated eigen-values ˆλk over 

subgroups. Note that the eigen-values are nonnegative for a nonnegative definite real symmetric matrix. Thus as with deff, if 

the estimated g-deffs are averaged over suitable subgroups ( j =1,..., J ) as defined by the GVF-type model (3.2), the 

smoothed estimate of Vψ

is obtained as 

 

 

 

ˆ

(

j k j k

Vψ λ )

k = ∑ ∑ q q′



 (3.4)

In the actual application of g-deff for smoothing, it may be better to avoid the effect of possibly different scaling of direct

estimates on eigen-values. To do this, smooth parts of Vˆ

ψ

expressed as a sandwich given by 

() ()

ˆˆ ˆ ( ) ( ˆ) ( V diag V R diag V ψ ψ = k ρ ) ψ k

 (3.5)

where R( ) ρˆ is the correlation matrix in the middle, and ρˆ denotes the estimated vector of correlations. Now we 

can use the deff-based GVF (equation 3.1) to smooth the left and right parts of the sandwich, while for the middle part, the gdeff based GVF (equation 3.2) can be used to extract the smoothed correlation vector 

K K( 1− ) / 2

ρ from the corresponding smoothed 

matrix R(this need not be a correlation matrix), and then obtainVψ

 as 

() () ( ) ( ) ( V diag V R diag V  ψ ψ = k ρ ) ψ k

 (3.6)

Some additional benefits of using the above sandwich form for smoothing are that it allows for more freedom in choosing 

parameters for smoothing because we have K deffs and K eigen-values based on the matrix R( ) ρˆ as opposed to only 

K eigen-values or g-deffs based on the matrix Vˆ

ψ

 , and that the two outside parts of the sandwich can be allowed to depend 

on model mean parameters as mentioned earlier under section 3.1. Note that in the g-deff based smoothing (3.4), it is not 

clear how one can make the eigen-value or the eigen-vector depend on unknown mean parameters. 

3.3 Aggregate-level vs. Unit-level Modeling

Finding eigen-values and eigen-vectors for g-deff based smoothing could be computationally difficult if the dimension K is 

large. For the aggregate-level modeling, Vˆ

ψ

is typically block-diagonal, and so the proposed method of g-deff is not too 

complicated computationally. However, for unit-level modeling, it remains nondiagonal because of summary statistics 

(obtained from EFs for fixed parameters) that aggregate over all areas. In these situations, one can take advantage of certain

patterns that typically arise in Vˆ

ψ

. Observe that the small areas or domains can often be grouped in practice into strata or 

superstrata, and then Vˆ

ψ

can be partitioned as 
 

 

ˆ

A B

V

B C ψ

⎛

= ⎜

′

⎝ ⎠

⎞

⎟

 (3.7)

where A is a high dimensional block diagonal matrix corresponding to domain or strata-level summary statistics, and C is a 

low dimensional matrix corresponding to summary statistics aggregated over all domains. Now decompose Vˆ

ψ

as 

1 1

1

,

1

A B AO I

H H H

B C O C BA B BA I

− −

−

⎛⎞⎛ ⎞ ⎛

= ′

⎜⎟⎜ ⎟ ⎜ ⎝⎠⎝ ⎠ ⎝ ′ ′ − −

O

−

⎞

= ⎟

′

⎠

 (3.8) 

where is expected to be stable, and the effect of instability of the matrix 1 C BA B−

− ′ H on Vˆ

ψ

is expected to be subsumed in 

that of the matrix A because covariance of the transformed vector Hψ turns out to be block-diag . So it is 

probably sufficient to just smooth

1

{, } AC BA B−

− ′

A to obtain A

 (using the g-deff idea) which being block-diagonal is not computationally

demanding, and then obtain Vψ



from the decomposition in (3.8) wherein the matrices H and are not modified or 

smoothed. 

1 C BA B−

− ′

4. Simulation Study

We design our study along the lines of Singh, Folsom, and Vaish (2003) which is based on Pfeffermann et al. (1998). 

Consider a universe of d = 1, ,  D strata (small areas) where D = 100 and let denote the number of population members 

in stratum- . In this simulation experiment, we set where is a constant and is obtained by 

truncating at 

Nd

d

*

0

(1 exp( )) NN u d = + d N0

*

d

u

~ (0, 0.2) du N ± 0.2 . For simplicity, we consider a single covariate super-population linear mixed 

model dk dk d dk 0 1 y x = β + βη ε + + where 0 β = 0.5 , 1 β = 1, ~ (0, 0.2) ηd N , ~ (0, 1) ε dk N , and 1, , . The covariate d

k =  N

dk d dk x = υ + δ where ~ (0, 0.1) υd N and ~ (0,1) δ dk N . We generate M = 150 population level data sets with common dk x

and where ’s are generated using =3000. Note that the substratum sizes vary over the 150 populations. We selected 

a sample from each of these populations so that the design was nonignorable. To select a sample with nonignorable design, 

we stratify the stratum- d population into two substrata 

Nd Nd N0

Ωd + with 0 dk ε > and Ωd − with 0 dk ε ≤ . Let , Nd + Nd − denote the 

sizes of these substrata and , denote the sizes of the simple random samples selected without replacement from these 

strata, respectively. Note that the substratum sizes vary across populations. We have and 

where . For 150 populations, we generate the corresponding 150 samples. In our simulation 

experiment, , so that we have common sample sizes from each area. 

d

n + dn −

100

1

d d

N N== ∑

100

1

d d

n

=

= ∑ n dd d +

nn n = − +

( , ) (2,10), (5,15) d d n n + − =

The results in this paper for the EFGL-unsmoothed method are not directly comparable to those in the Singh, Folsom, and 

Vaish (2003) paper because here we are using the version of EFGL based on the summary statistics ψ which consists of 

; is excluded because it is constant and equal to for the stratified SRS design. Note that the simpler 

case of aggregate-level modeling was not considered in this limited simulation study. Also, we didn’t consider the sandwich 

2 , ,

{ , ,,} y d x d xy xt t tt c d,

t Nd
-type g-deff method (3.5) for smoothing Vˆ

ψ but used the initial g-deff method (3.4). To avoid computational complexity, we 

used the smoothing based on the decomposition (3.8) where A is a block diagonal matrix of 100 blocks, each of dimension 

2x2 corresponding to each small area, and is only a 2x2 matrix. In order to use EFGL under a HB framework, 

customary priors for 

1 C BA B−

− ′

β and 2 σ η were chosen as (estimation of the parameter 2 σ ε was not considered as it was not needed) 

 

 

( ) ( ) 0

2

~ ,~ 2 0 U R IG , β ση ν σ

2

2 η (4.1)

where , a small positive number, and 

0

2

.01 ση = 0 ν = 3 , the smallest possible integral value greater than 2 which gives a 

practically meaningful interpretation of the posterior mean of 2

[ η σ η] obtained as a weighted average 

( ( ) ( ) ) ( )

0

2 22

0

| 2 d

d

E DD D ση η ν σ ν η

⎡ ⎤ = + − +

⎣ ⎦ ∑ 0

2 η − (4.2) 

The results from the simulation study are presented in Tables 1 and 2. In terms of point estimation and Standard deviation, 

both methods (EFGL-u for unsmoothed, and EFGL-s for smoothed) perform very similarly as seen from Table 1. However, 

in terms of coverage probabilities and length of prediction intervals, the unsmoothed method EFGL-u starts showing 

deterioration as seen from Table 2. It is expected that variance of SAEs would suffer from downward bias due to 

underestimation of their variances by regarding Vˆ

ψ or its smoothed version Vψ

 as known. However, it is interesting to note 

that for smaller sample sizes, there is marked difference between the coverage probability and interval width between the 

smoothed and unsmoothed methods and smoothing is seen to correct in the right direction. Nevertheless, when sample sizes 

for small areas are rather low as seen in the lower half of Table 2, even after smoothing, there is significant undercoverage of 

prediction intervals. This suggests that Vψ



is still not sufficiently smooth and perhaps collapsing of small areas with similar 

random effects as suggested in Singh, Folsom and Vaish (2003) might alleviate this problem. In future, this and other issues 

related to the performance of smoothing methods for other choices of area sample sizes will be investigated further. 

5. Summary and Concluding Remarks 

In this paper we considered the problem of smoothing the error covariance in small area modeling of survey data. This 

problem has been around for quite some time. The proposed method of smoothing the error covariance based on deff and gdeff provides a simpler alternative to other methods including modeling the error covariance, and seemed to perform well in a 

limited simulation study. An important consideration in using the proposed method is that similar to GVF modeling, the 

underlying assumptions are quite mild; this is likely to be attractive to practitioners. An interesting finding of the simulation 

study was that although smoothing seemed to help correct underestimation of variance of SAEs, it may not be sufficient for 

areas with very small sample sizes. In future, it would be interesting to see if suitable collapsing of areas would help to 

provide a sufficiently stable smoothed estimate of error covariance as well as help in justifying the normal approximation for 

summary statistics. We also note that the basic idea of covariance smoothing proposed in this paper is applicable to other 

SAE problems as well involving spatial and temporal modeling. 
Finally, it is interesting to note that the use of GVF was also identified by Eltinge et. al (2002) in trying to stabilize the 

covariance matrix in the context of quadratic-score type tests for model checks in analysing survey data, and considered their 

relationships with Rao-Scott corrections to chi-square tests. However, they didn’t consider the use of g-deff in GVF-type 

modeling that seems to provide a simple but important approach for the problem of SAE considered in this paper. 

Acknowledgments: The first author’s research was supported in part by a grant from Natural Sciences and Engineering 

Research Council of Canada held at Carleton University, Ottawa under an adjunct research professorship. 

References 

Binder, D.A. (1983), “On the Variances of Asymptotically Normal Estimators from Complex Surveys,” International 

Statistical Review, 51, pp. 279-292. 

Eltinge, J.L., Cho, M.J., and Hinrichs, P. (2002). Use of Generalized variance functions in Multivariate Analysis, 

Proceedings of American Statistical association, Section on Survey Research Methods. 

Fay, R.E. and Herriot, R.A. (1979), “Estimates of Income for Small Places: An Application of James-Stein Procedures to 

Census Data,” Journal of the American Statistical Association, 74, pp. 269-277. 

Godambe, V.P. (1955). A unified theory of sampling from finite populations. JRSS (B), 17, 269-278. 

Godambe, V.P. (1960). A new approach to sampling from finite populations I, II (with discussion), JRSS (B), 28, 310-328. 

Godambe and Thompson, M.E, (1986), “Parameters of Super population and Survey Population, Their Relationship and 

Estimation,” International Statistical Review, 54, pp. 127-38. 

Otto, M.C., and Bell, W.R. (1995), “Sampling Error Modeling of Poverty and Income Statistics for States,” Proceedings of 

the Government Statistics Section, American Statistical Association, pp. 160-165. 

Pfeffermann, D. and Sverchkov, M. (1999). Parametric and semi-parametric estimation of regression models fitted to survey 

data. Sankhya: The Indian Journal of Statistics, Ser. B, 61, 166-186. 

Pfeffermann, D. and Sverchkov, M. (2003). Fitting generalized linear models under informative probability sampling. In: 

Analysis of Survey Data, eds. C.J. Skinner, and R.L. Chambers, New York: Wiley, 175-195.

Pfeffermann, D., Skinner, C.J., Holmes, D.J., Goldstein, H., and Rasbash, J. (1998), “Weighting for Unequal Selection 

Probabilities in Multilevel Models,” Journal of the Royal Statistical Society, B, 60, pp. 23-40 

Rao, C.R. (1973). Linear Statistical Inference and Its Applications, Second Ed., New York: John Wiley. 

Rao, J.N.K. (2003). Small Area Estimation, New York: John Wiley 

Rao, J.N.K., and Scott, A.J. (1981) The analysis of categorical data from complex sample surveys: Chi-squared tests for 

goodness-of-fit and independence in two-way tables, Journal of the American Statistical Association, 76 221-230. 

Singh, A.C., Folsom, R.E., Jr., and Vaish, A.K. (2002). Estimating function-based approach to hierarchical Bayes Small Area 

Estimation from survey data, International Conference on Recent Advances in Survey Sampling, , in honor of J.N.K. Rao’s 

65th birthday, July 10-13 
Singh, A.C., Folsom, R.E., Jr., and Vaish, A.K. (2003). Hierarchical Bayes small area estimation for survey data by EFGL: 

the method of estimating function-based Gaussian Likelihood (with discussion), FCSM Statistiscal Policy Working paper 

#36, pp. 47-74 (www.fcsm.gov/working-papers/) 

Singh, A.C., You, Y., and Mantel, H.J. (2005). Use of generalized design effect for variance function modeling in small area 

estimation from survey data, Statistical Society of Canada, 33rd Annual Meeting, Saskatoon, Saskatchewan, June 12-15. 

You, Y, Rao, J.N.K., and Gambino, J.G. (2003). Model-based unemployment rate estimation for the Canadian Labour Force 

Survey: A hierarchical Bayes approach, Survey Methodology, 29, 25-32. 
Table 1: Average Posterior Mean and Standard Deviation for Model Parameters: Nonignorable Sample

Design 

( , ) (5,15 d d n n + − = )

Average Posterior Mean Average Posterior 

Parameter Standard Deviation 

(True 

Value) EFGL-u EFGL-s EFGL-u EFGL-s 

β 0 (0.5) 0.4927 0.4933 0.0478 0.0477 

β1 (1.0) 0.9984 0.9983 0.1225 0.1226 

2 σ η (0.2) 0.1970 0.1963 0.0320 0.0320 

( , ) (2,10 d d n n + − = )

Average Posterior Mean Average Posterior 

Parameter Standard Deviation 

(True 

Value) EFGL-u EFGL-s EFGL-u EFGL-s 

β 0 (0.5) 0.4816 0.4790 0.0505 0.0504 

β1 (1.0) 0.9852 0.9839 0.1011 0.1030 

2 σ η (0.2) 0.2048 0.1980 0.0347 0.0345 

Table 2: 95% Coverage Probability and Ratio of Prediction Interval (PI) Widths: Nonignorable Sample 

Design 

( , ) (5,15 d d n n + − = )

Percentiles and Coverage Probability Ratio of Average PI Widths 

Means 

over Small Areas EFGL-u EFGL-s EFGL-s/EFGL-u 

75% 0.9333 0.9533 1.02 

Mean 0.9183 0.9386 1.01 

25% 0.9000 0.9267 1.00 

( , ) (2,10 d d n n + − = )

Percentiles and Coverage Probability Ratio of Average PI Widths 

Means 

over Small Areas EFGL-u EFGL-s EFGL-s/EFGL-u 

75% 0.8433 0.8833 1.06 

Mean 0.8252 0.8659 1.06 

25% 0.8067 0.8467 1.04