Complex Survey Variance and Design Effects in R

using the Rstan and Survey packages

Matt Williams1 Terrance Savitsky2

1National Center for Science and Engineering Statistics

National Science Foundation

mrwillia@nsf.gov

2Office of Survey Methods Research

Bureau of Labor Statistics

Savitsky.Terrance@bls.gov

GASP

Sept 23, 2019

1 / 10
Population Inference from Complex Survey Samples

◮ Goal: perform inference about a finite population generated from an

unknown model, Pθ0

.

◮ Data: from under a complex sampling design distribution, Pν

◮ Probabilities of inclusion πi are often associated with the variable of

interest (purposefully)

◮ Sampling designs are “informative”: the balance of information in the

sample 6= balance in the population.

◮ Biased Estimation: estimate Pθ0 without accounting for Pν.

◮ Use inverse probability weights wi = 1/πi to mitigate bias.

◮ Incorrect Uncertainty Quantification:

◮ Failure to account for dependence induced by Pν leads to standard

errors and confidence intervals that are the wrong size.

2 / 10
Variance Estimation

◮ The de-facto approach:

◮ approximate sampling independence of the primary sampling units

(Heeringa et al. 2010).

◮ within-cluster dependence treated as nuisance

◮ Two common methods:

◮ Taylor linearization and replication based methods.

◮ A variety of implementations are available (Binder 1996, Rao et al.

1992).

3 / 10
Taylor Linearization

Let yij , Xij , and wij be the observed data for individual i in cluster j of

the sample. Assume the parameter θ is a vector of dimension d with

population model value θ0.

1. Approximate an estimate ˆθ, or a ‘residual’ (

ˆθ − θ0), as a weighted

sum:

ˆθ ≈

P

i,j wijzij (θ) where zij is a function evaluated at the

current values of yij , Xij , and ˆθ.

2. Compute the weighted components for each cluster (e.g., primary

sampling units (PSUs)): ˆθj =

P

i wijzij (θ).

3. Compute the variance \ between clusters:

V ar(

ˆθ) = 1

J−d

PJ

j=1(

ˆθ − ˆθj )(ˆθ − ˆθj )

T

4. For stratified designs, compute ˆθs and V ar\(

ˆ \ θs) within strata and sum

V ar(

ˆθ) = P

s

V ar\(

ˆθs).

4 / 10
Replication

Let yij , Xij , and wij be the observed data for individual i in cluster j of

the sample. Assume the parameter θ is a vector of dimension d with

population model value θ0.

1. Through randomization (bootstrap), leave-one-out (jackknife), or

orthogonal contrasts (balanced repeated replicates), create a set of K

replicate weights (wi)k for all i ∈ S and for every k = 1, . . . , K.

2. Each set of weights has a modified value (usually 0) for a subset of

clusters, and typically has a weight adjustment to the other clusters

to compensate: P

i∈S

(wi)k =

P

i∈S wi for every k.

3. Estimate ˆθk for each replicate k ∈ 1, . . . , K.

4. Compute the variance \ between replicates:

V ar(

ˆθ) = 1

K−d

PK

k=1(

ˆθ − ˆθk)(ˆθ − ˆθk)

T

.

5. For stratified designs, generate replicates such that each strata is

represented in every replicate.

5 / 10
Challenges

There are two notable trade-offs associated with these methods:

◮ Taylor linearization: value ˆθ computed once then used in a plug in for

zi(θ).

◮ Replication methods: estimate ˆθk computed K times.

◮ Sizable differences in computational effort

◮ Replication methods: no derivatives are needed.

◮ Taylor linearization: requires the calculation of a gradient to derive the

analytical form of the first order approximation zi(θ).

◮ This poses significant analytical challenges for all but the simplest

models.

6 / 10
Some Improvements

◮ Abstraction of Derivatives (less analytic work for Taylor Linearization)

◮ Recent advances in algorithmic differentiation (Margossian 2018),

allows us to specify the model as a log density but only treat the

gradient in the abstract without specifying it analytically.

◮ Implemented in Stan and Rstan (Carpenter 2015, Stan Development

Team 2016)

◮ Hybrid Approach or Taylor Linearization for replicate designs (less

computation for Replication approaches)

◮ Survey package (Lumley 2016) to calculate replication variance of

gradient

◮ Use plug in for θ, only estimate once

7 / 10
Example: Design Effect for Survey-Weighted Bayes

◮ Williams & Savitsky (2018): https://arxiv.org/abs/1807.11796

◮ Pseudo posterior ∝ Pseudo Likelihood × Prior

π

π

(λ|y, w ̃ ) ∝

"Yn

i=1

π (yi|λ)

w ̃i

#

π (λ)

◮ Variances Differ:

◮ Weighted MLE: H

−1

θ0

J

π

θ0H

−1

θ0

(Robust)

◮ Weighted Posterior: H

−1

θ0

(Model-Based)

◮ Adjust for Design Effect: R

−1

2 R1

◮ ˆθm ≡ sample pseudo posterior for m = 1, . . . , M draws with mean  ̄θ

◮ ˆθ

a

m =



ˆθm −  ̄θ



R

−1

2 R1 +  ̄θ

◮ where R′

1R1 = H

−1

θ0

J

π

θ0H

−1

θ0

◮ R′

2R2 = H

−1

θ0

8 / 10
R Code Schematic

Input R Code Output

Stan Model

sampling

(rstan)

svrepdesign

(survey)

Survey

Design

reps

 ̄θ

grad log prob

(rstan)

ˆθm

Hˆ

θ

withReplicates

(survey) Jˆπ

θ

aaply

(plyr)

ˆθ

a

m

9 / 10
References I

Binder, D. A. (1996), ‘Linearization methods for single phase and two-phase samples: a

cookbook approach’, Survey Methodology 22, 17–22.

Carpenter, B. (2015), ‘Stan: A probabilistic programming language’, Journal of Statistical

Software .

Heeringa, S. G., West, B. T. & Berglund, P. A. (2010), Applied Survey Data Analysis, Chapman

and Hall/CRC.

Lumley, T. (2016), ‘survey: analysis of complex survey samples’. R package version 3.32.

Margossian, C. C. (2018), ‘A review of automatic differentiation and its efficient

implementation’, CoRR abs/1811.05031.

URL: http://arxiv.org/abs/1811.05031

Rao, J. N. K., Wu, C. F. J. & Yue, K. (1992), ‘Some recent work on resampling methods for

complex surveys’, Survey Methodology 18, 209–217.

Stan Development Team (2016), ‘RStan: the R interface to Stan’. R package version 2.14.1.

URL: http://mc-stan.org/

Williams, M. R. & Savitsky, T. D. (2018), ‘Bayesian uncertainty estimation under complex

sampling’, arXiv preprint arXiv:1807.11796 .

10 / 10