Bayesian stratified sampling for 

establishment surveys with 

uncertain design parameters

Jonathan Mendelson1,2 and Michael Elliott3

FCSM, October 25, 2022

1

1 Bureau of Labor Statistics

2 University of Maryland

3 University of Michigan

Content does not represent BLS policy
Bottom line up front

• Imprecisely estimated survey design 

parameters could harm sample efficiency

• There is a Bayesian approach to sample 

design, which accounts for this

• We identify the Bayesian optimal design in a 

particular establishment survey context

– Outperforms Neyman-HT in simulation

– Performs similarly or better than the main 

model-assisted approach considered

2
I. Introduction

3
Introduction

• Sample designs often assume population 

characteristics are known.

• In practice, some are typically estimated.

• Example. Optimal STSRS for estimating finite 

population mean via separate ratio model

– Theory: nnh ∝ NNhSSddh/ cch (Cochran, 1977)

• SSdd is the stratum SD of a residual term

– Practice: nnh ∝ NNhSŜ

ddh/ cc

̂

h

• Typically, little attention is given to the effect of 

imperfect information on sample design.

4 In the above, SSdd

2 = NNh − 1 −1 ∑ii=1

NNh YYhii − YYh

XXh

XXhii

2

, and

NNh, nnh, cch are the population size, sample size, and per-unit cost in stratum h ∈ 1,... , HH .
Selected Bayesian design literature

• Bayesian optimal experimental design (Lindley, 1972) 

can be applied to STSRS sample allocation

– Flexible approach; accommodates uncertainty 

• Draper & Guttman (1968) consider continuous data

– Assumes use of pilot study data

– Special case leads approximately to Neyman allocation

– However, D&G assume fixed strata means and variances

• Rao & Ghangurde (1972) consider categorical data

– Assumes Dirichlet-multinomial model

– Applicability for continuous, skewed distributions?

5
Heteroscedasticity and design

• Consider {XXii, YYii; ii = 1, ... , NN}, where

– YYii = ββXXii + εεii

– EMM εεii = 0

– VarMM εεii = σσ2XXii

bb; known bb, XXii > 0

– Independent εεii’s

• “b” (coefficient of heteroscedasticity) can 

meaningfully affect optimal allocation

– PPS/GREG strategy: ππii ∝ XXii

bb/2 (e.g., SSW, 1992)

6
Heteroscedasticity, visualized

7

Data source: National Hospital Discharge Survey of 1968 (via PracTools)

• See Henry & Valliant (2009) for more real examples
Bayesian decision theory for optimal 

experimental design

• Lindley (1972) treats as a two-part decision:

– Choose the experiment, ee ∈ EE (e.g., ee = {nnh})

• This results in the sample (data), xx ∈ XX

– Translate the data into a terminal decision

• For example, compute estimate θθ� for parameter θθ ∈ Θ

(e.g., finite population mean)

• Define a loss function of the form LL θθ

�

, θθ, ee, xx

• Lindley suggests finding optimal θθ

�

, ee via

mm

ee �

XX

mm

θθ̂ �

Θ

LL θθ̂,θθ, ee, xx pp θθ xx, ee pp xx ee dd dd

8 In the above, pp(θθ|xx, ee) is the posterior pdf and pp(xx|ee) is a conditional pdf.
Our research

• We consider optimal STSRS design while 

accounting for heteroscedastic errors and 

uncertain design parameters

– We aim for weaker assumptions than some 

previous Bayesian work

– We accommodate uncertain design parameters 

via Bayesian decision theoretic formulation

9
II. Problem set-up and 

Bayesian analysis

10
Problem set-up

• Study design:

– Pilot is only used for designing main study

– Strata defined upfront

• Model: YYhii = ααhXXhii + εεhii, where εεhii ~

iiii NN 0, vvhXXhii

bb

– Known XXhii > 0; known bb

• Prior (diffuse): ππ ααh,

1

vvh

∝ ∏h=1

HH vvh

11

Pilot (D1) Main (D2)

Decide nnh
Overview: our Bayesian decision theoretic 

analysis for the finite population mean

Above, ee = nn 12 h , YY� = ∑h=1

HH ∑ii=1

NNh YYhii, and ss2h refers to the main study sample in stratum h.

1. Objective: LL YY�, ��

YY, ee,DDD = ��

YY − YY� 2

– Minimized when ��

YY = E YY� DDD, ee, bb

2. Posterior loss is Var YY� DDD, ee, bb

– Apply Ericson (1969) to obtain

3. Preposterior analysis: average over future data (DDD|DDD)

– Consider uncertainty with respect to:

• Posterior for parameters given pilot, ααh, vvh |DDD

• Sample indicators, {ss2h}

• Model uncertainty given above, DDD| ααh,vvh,DDD, ss2h

– Results provided in paper

4. Optimize via mathematical programming 
III. Simulation

13
Simulation design: compare strategies 

across a series of artificial populations

• “Strategy” denotes allocation + estimator

• We generated PP = 90 bivariate populations, and 

applied each strategy RR = 1000 times

• For population pp, simulation rr:

– Draw an equally allocated pilot sample of mm = 75 units

– For strategy aa:

• Allocate and draw a main study sample of nn = 500 units

• Obtain point estimate and 95% CI

• Compare strategies’ RMSE, bias, and CI coverage/width

– For instance: rr YY� pp,aa =

1

RR ∑rr=1

RR YY� pp,aa

rr − YY pp

2

14
We considered three size measures

15

MOS1 MOS

2 MOS

3

0 1 2 3

0

1

2

3

0

1

2

3

0

1

2

3

X

density

(Vertical lines denote strata boundaries)

Distributions of simulated stratified size measures, by MOS
We considered 30 structures for 

YYhii|XXhii~NN(ααhXXhii, vvhXXhii

bb )

16

• 5 levels of bb considered: bb ∈ 0, 0.5, 1, 1.5, 2

• 6 choices of ααh, vvh , where {vvh} were chosen 

as to approximately yield target correlations

Scenario αα1 αα2 αα3 αα4 αα5 ρρ1 ρρ2 ρρ3 ρρ4 ρρ5

1. Baseline 1 1 1 1 1 0.7 0.7 0.7 0.7 0.7

2. Lower correlations 1 1 1 1 1 0.5 0.5 0.5 0.5 0.5

3. Higher correlations 1 1 1 1 1 0.9 0.9 0.9 0.9 0.9

4. Increasing correlations, 

fixed slopes 1 1 1 1 1 0.5 0.6 0.7 0.8 0.9

5. Fixed correlations, 

decreasing slopes 1.4 1.2 1 0.8 0.6 0.7 0.7 0.7 0.7 0.7

6. Increasing correlations, 

decreasing slopes 1.4 1.2 1 0.8 0.6 0.5 0.6 0.7 0.8 0.9
We compared several strategies

• We focused on three main strategies:

– Neyman plug-in/HT estimator (N-HT)

– Cochran plug-in/separate ratio estimator (C-SR)

– Bayesian allocation/prediction estimator (B-P)

• We also considered three rule-of-thumb 

allocations suggested or implied by Cochran for 

different levels of bb (used SR estimator for each)

– nnh ∝ NNh

– nnh ∝ NNh XX�

h

– nnh ∝ NNhXX�

h

17
Simulation results: main strategies

• The three main strategies:

– were approximately unbiased; and

– had near-nominal coverage for 95% CIs.

• Therefore, we focused on analyzing RMSE

– Findings on RMSE were paralleled by analogous 

findings for CI relative width

18
B-P consistently outperformed N-HT

19

Relative increase in RMSE from N-HT versus B-P (among MOS1 pops)

Scenario for {ρρh, ααh}

1. Baseline

2. Lower 

corrs

3. Higher 

corrs 4. Inc. corrs

5. Dec. 

slopes

6. Inc. corrs, 

dec. slopes

b = 0 65% 32% 175% 72% 69% 77%

b = 0.5 44% 18% 138% 44% 51% 36%

b = 1 46% 24% 127% 35% 40% 28%

b = 1.5 52% 24% 139% 41% 51% 38%

b = 2 64% 44% 159% 75% 68% 62%

• Use of N-HT led to RMSE 11%–175% higher than B-P for 

individual populations (MOS1 pops displayed below)

• Results varied greatly by assumptions for ff(YYhii|XXhii)

– Compare 2nd and 3rd data columns below
B-P did about as well or better than C-SR, with 

marked differences across populations

20

Relative increase in RMSE from C-SR versus B-P (among MOS1 pops)

Scenario for {ρρh, ααh}

1. Baseline

2. Lower 

corrs

3. Higher 

corrs 4. Inc. corrs

5. Dec. 

slopes

6. Inc. corrs, 

dec. slopes

b = 0 19% 10% 18% 32% 29% 40%

b = 0.5 7% 3% 0% 2% 12% 7%

b = 1 6% 7% 3% 1% 4% 3%

b = 1.5 11% 9% 5% 5% 6% 5%

b = 2 23% 21% 23% 25% 24% 29%

• B-P showed the greatest advantage for a subset of MOS1 

scenarios (top and bottom rows below)

• In contrast, differences were fairly muted for most MOS2 

and MOS3 populations, which had less skewness
B-P sometimes produced more stable 

allocations than the main alternatives

• Differences in allocations’ stability were 

starkest for MOS1, bb = 2 pops, for instance:

21

Allocation summary statistics by allocation and stratum

Population 25 (MOS1, b=2, baseline ρρh, ααh scenario)

Neyman Cochran Bayesian

h E(nnh) ss nnh E(nnh) ss nnh E(nnh) ss nnh

1 149 48 135 50 112 19

2 90 20 92 22 98 17

3 76 16 80 19 85 15

4 85 17 85 19 91 15

5 101 22 107 24 114 19
Performance was mixed for rule of 

thumb strategies

• C-SR and B-P strategies, which incorporate 

pilot data for allocation, consistently did as 

well or better than the RT-SR strategies

– nnh ∝ NNhXX�h performed quite badly in some 

situations (e.g., RMSE 82%–204% higher than B-P 

for bb = 2, MOS1 populations)

– In contrast, nnh ∝ NNh XX�h had reasonable 

performance for a subset of the bb = 1 scenarios 

(depending on the ααh and ρρh)

22
IV. Application

23
We applied our methods to analyzing 

tax returns of public charities

• Source: IRS Form 990 data (National Center for 

Charitable Statistics [NCCS], Urban Institute)

– Analyzed 140,858 domestic operating public 

charities meeting inclusion criteria

– X = log revenue, 2008

– Y = log revenue, 2013

• Unstratified MCMC analysis yielded bb

�

= 0.55

and 95% CI of (0.25, 0.66)

24
NCCS application (continued)

• We formed 24 strata based on nonprofit 

sector (8 groups) by revenue class (3 groups)

• Methods paralleled earlier simulation

– RR = 10,000 equally allocated pilots of 360 units 

used to design main studies of 1,800 units

– Compared RMSE, relative bias, CI properties

25
C-SR and B-P again outperformed N-HT

• C-SR and B-P offered substantial reduction in 

RMSE than N-HT

• All three methods were approximately 

unbiased and had near-nominal CI coverage

26

Table. NCCS Simulation Results

Strategy Relative RMSE 1000*RelBias CI Coverage (%) 1000*CI RelWidth

N-HT 1.427 -0.00 94.7 6.48

C-SR 1.014 -0.01 94.8 4.57

B-P 1.000 -0.00 95.5 4.63

Note: RMSE is displayed relative to that of the B-P strategy.
V. Discussion

27
We provided a Bayesian approach to 

sample design for our problem

• We considered STSRS designs for establishments

– Allow for heteroscedastic errors  improved realism

– Problem formulated via Bayesian decision theory

– We derived the approximate expected posterior 

variance, which is then minimized

• We assessed performance via simulation to 

artificial and real data

– The proposed B-P strategy provided substantial gains 

versus design-based approach

– B-P strategy did as well or better than the main 

model-assisted approach considered

28
Potential future directions

• Consider other population structures, 

including those not following our model

• Compare to additional sampling strategies

• Extend to scenarios where “b” is unknown

• Consider other loss functions

• Identify other ways to express prior 

knowledge (e.g., in absence of pilot)

29
Jonathan Mendelson

PhD candidate, JPSM

Research Statistician, BLS (OSMR/BSRC)

mendelson.jonathan@bls.gov

30

Comments? Questions?