Hierarchical models in the production of

official statistics: a discussion of some

practical aspects

Andreea L. Erciulescu

Westat

GASP, Washington, DC

September 23, 2019

The views presented in this paper are those of the author(s) and do not represent

the views of any Government Agency/Department or Westat

Improving Lives Through Research 1
Why hierarchical models in the production of official

statistics?

◮ need to account for various sources of error

◮ sampling error, measurement error, linking error

◮ need for a transparent, reproducible and validated process

◮ analytic techniques

◮ need for measures of uncertainty

◮ entire distribution

Improving Lives Through Research 2
Hierarchical model - general form

Data model: [y|θ, ψ]

Process model: [θ|ψ]

Parameter model: [ψ]

Joint distribution

[y, θ, ψ] = [y|θ, ψ][θ|ψ][ψ]

Predictive distribution

[θ|y, ψ]

Three cases

◮ known ψ

◮ unknown, fixed ψ

◮ unknown, random ψ

Improving Lives Through Research 3
Why model-based small area estimation (SAE) in the

production of official statistics?

◮ need to account for various sources of error

◮ need for a transparent, reproducible and validated process

◮ need for measures of uncertainty

◮ integration of data from multiple sources

◮ observed survey data within area + auxiliary data within area +

information across all areas

◮ growing demand for granular statistics

◮ quantities of interest: i.e. totals, means, ratios

◮ ‘area’ / ‘domain’: i.e. geography, socio-economic status, occupation

◮ increasing costs of data collection

◮ ‘small’: amount of survey data available for estimation within a given

area (realized sample size as small as zero)

Improving Lives Through Research 4
Small domain hierarchical model - general form

Sampling model: [yi|θi, xi, ψθ,i, ψy,i]

Linking model: [θi|xi, ψθ]

Parameter model: [ψθ,i, ψy,i]

Joint distribution

[yi, θi, xi, ψθ, ψy,i] = [yi|θi, xi, ψθ, ψy,i][θi|xi, ψθ][ψθ, ψy,i]

Predictive distribution

[θi|yi, xi, ψθ, ψy,i]

◮ typically known xi

, ψy,i (area-level models)

◮ still three cases

◮ known ψθ

◮ unknown, fixed ψθ

◮ unknown, random ψθ

Improving Lives Through Research 5
Examples of model-based SAE in government programs

◮ Chilean Ministerio de Desarrollo, World Bank: poverty mapping

◮ Casas Cordero Valencia et al. (2016); Elbers et al. (2003)

◮ Census Bureau: income and poverty measures

◮ Bell et al. (2016)

◮ National Agricultural Statistics Service:

◮ cash rental rates, Erciulescu et al. (2018)

◮ crops production, Erciulescu et al. (2019a)

◮ agricultural labor wages, Erciulescu (2018)

◮ Organisation for Economic Co-operation and Development: adult

competency

◮ Krenzke et al. (2019)

◮ Bureau of Labor Statistics: employee compensation components

◮ Erciulescu and Opsomer (2019)

Improving Lives Through Research 6
A small domain hierarchical Bayes multivariate model

for employee compensation components*

Sampling Model: yi|θi ∼ N(θi, Σei)

Linking Model: θi|(β, Σv ) ∼ N(xiβ, Σv )

Independent priors: π(β, Σv ) = π(β)π(Σv )

◮ domain i, cross-tabulation of census divisions, 6-digit SOC system

codes, work levels, binary characteristics

◮ θi

, quantities of interest, wage and benefits

◮ yi

, direct survey estimates, domain-level wage and benefits direct

survey estimates

◮ xi

, known covariates, selected using sample domains definitions

◮ β, regression coefficients

◮ Σv , linking model variance-covariance matrix

◮ Σei, known survey variance-covariances matrices

* Erciulescu and Opsomer, 2019

Improving Lives Through Research 7
A common approach for fit and prediction

Fit

◮ Markov chain Monte Carlo (MCMC): multiple chains, iterations,

burn-in, thining; keep R samples for inference

Prediction

◮ in-sample domain i

◮ R samples θiζ , ζ = 1, · · · , R, from [θi|(yi, ψθ, ψy,i)]

◮ not-in-sample domain i

′

:

◮ generate new R samples θi

′

ζ

, ζ = 1, · · · , R, from [θi|xi, ψθ,ζ ]

◮ small domain posterior means: R

−1 PR

ζ=1

θiζ

◮ small domain posterior variances: R

−1 PR

ζ=1



θiζ − R

−1 PR

ζ=1

θiζ

2

◮ small domain posterior p quantiles: θi(p)

, θi = (θi1, · · · , θiR)

Practical challenges: number of MCMC samples, large number of domains (in-sample

and not-in-sample), variable selection

Improving Lives Through Research 8
Prediction for other functions

For example,

◮ sum of θi components (total employee compensation, proportions of

adult literacy/numeracy scores in prespecified ranges)

◮ inference using R samples θ

S

iζ

:= θiζ,1 + θiζ,2, ζ = 1, · · · , R

◮ ratio of θi components (labor wage as ratio of total income to total

hours worked, crop yield as the ratio of total production to total

harvested acres)

◮ inference using R samples θ

R

iζ

:= θiζ,1/θiζ,2, ζ = 1, · · · , R

◮ aggregates of θi for pre-defined domains (2-digits SOC codes,

county-agricultural district-state-census division-nation)

◮ inference using R samples

θ

D

iζ

:= P

i∈pre-defined domain θiζ , ζ = 1, · · · , R

Improving Lives Through Research 9
Examples of model validation

Internal

◮ mixing and convergence diagnostics: Rˆ, MC effective sample size,

autocorrelation

◮ residuals diagnostics: unconditional and conditional

◮ posterior predictive checks: indicator, correlation, deviance, residuals

◮ alternative model specifications: prior distributions

External

◮ predictions versus direct estimates

◮ predictions for not-in-sample domains versus predictions for

in-sample domains

◮ cross-validation

Practical challenges: autocorrelation, cross-validation, visualization, storage

Improving Lives Through Research 10
Simulation study: Data generation model

yi|θi ∼ N(θi, Σei)

θi|(β, Σv ) ∼ N(xiβ, Σv )

π(β, Σv ) = π(β)π(Σv )

◮ i = 1, · · · , m

◮ yi = (yi,1, yi,2)

◮ xi

, two identical rows xi,row = (xi,1, xi,2, xi,3, ..., , xi,(p+1))

◮ xi,1 = 1, xi,2 ∼ Beta(2, 4), xi,k ∼ N(μxk , σ2

xk ), k = 3, ...,(p + 1)

◮ μxk ∼ Unif (1, 50), σxk ∼ Unif (1, 10), k = 3, ...,(p + 1)

◮ β = (β1, β2), βj = (βj,1, βj,2, βj,3, ..., , βj,(p+1)), j = 1, 2

◮ β1,1 = 1, β1,2 = 10, β1,k ∼ Unif (1, 5) + N(0, 1), k = 3, ...,(p + 1)

◮ β2,1 = 1, β2,2 = 10, β2,k ∼ Unif (1, 2) + N(0, 1), k = 3, ...,(p + 1)

◮ Σv ∼ inverse − Wishart(3, I2)

◮ diag(Σei) = (σ

2

ei,1

, σ2

ei,2

), σ

2

ei,j = exp(log(xi,jβj) + N(0, 1)), j = 1, 2

◮ cor(yi,1, yi,2) ∼ Unif (0.3, 0.9)

Improving Lives Through Research 11
Simulation study: Fitted models

yi|θi ∼ N(θi, Σei)

θi|(β, Σv ) ∼ N(xiβ, Σv ) ↔ θi|(β, vi) = xiβ + vi, vi|Σv ∼ N(0, Σv )

π(β, Σv ) = π(β)π(Σv )

Software

◮ R STAN, R JAGS

Bayesian specification

◮ hierarchical Bayes: π(β) = N(0, 104

) and

Σv ∼ Inverse − Wishart(3, I)

◮ empirical Bayes: π(β1) = N(0, 104

); β−1 set equal to the least

squares estimates, based on a multiple linear regression model and

Σv ∼ Inverse − Wishart(3, I)

Practical challenges: prior distributions for the linking model variance-covariance

components (Inverse-Gamma/Uniform/Cauchy/F; Inverse-Wishart/LKJ)

Improving Lives Through Research 12
Simulation study: Fitted models specifications

Table 1: Model Specifications

Parameters Model Software # MC samples / chain

(m, p) (start, burn-in, thin)

(100, 10) S1 STAN (3000, 1000, 10)

J1 JAGS (3000, 1000, 10)

J1l JAGS (30000, 10000, 10)

(1000, 100) S2 STAN (3000, 1000, 10)

J2 JAGS (3000, 1000, 10)

J2l JAGS (30000, 10000, 10)

(10000, 100) S3 STAN (3000, 1000, 10)

J3 JAGS (3000, 1000, 10)

J3l JAGS (30000, 10000, 10)

Improving Lives Through Research 13
Simulation study: JAGS HB model specification

Improving Lives Through Research 14
Simulation study: JAGS EB model specification

Improving Lives Through Research 15
Simulation study: STAN HB model specification

Improving Lives Through Research 16
Simulation study: STAN EB model specification

Improving Lives Through Research 17
Simulation study: Computational time results

Table 2: Computational Time Summaries

Model Time (seconds)

HB EB

S1 525 284

J1 5 2

J1l 35 11

S2 32744 14702

J2 1438 24

J2l 14227 109

S3 354416 76992

J3 18210 273

J3l Too long ...

Improving Lives Through Research 18
Simulation study: Convergence results

Table 3: Convergence Results Summaries HB

Model Approach Rˆ MC Effective Sample Size

(min, median, max) (min, median, max)

S1 HB (0.995, 1.000, 1.072) (44, 571, 825)

EB (0.995, 1.005, 1.086) (35, 342, 810)

J1 HB (0.999, 1.004, 1.070) (44, 600, 600)

EB (0.999, 1.023, 1.216) (14, 440, 600)

J1l HB (1.001, 1.001, 1.030) (130, 5800, 6000)

EB (1.001, 1.003, 1.014) (190, 3400, 7500)

S2 HB (0.995, 0.999, 1.289) (14, 580, 908)

EB (0.995, 0.999, 1.262) (9, 583, 802)

J2 HB (0.999, 1.010, 1.818) (6, 600, 600)

EB (1.037, 1.116, 2.029) (5, 79, 600)

J2l HB (1.001, 1.002, 1.140) (33, 6000, 6000)

EB (1.001, 1.002, 1.138) (59, 3500, 7500)

S3 HB (0.995, 1.000, 4.097) (2, 582, 989)

EB (0.995, 1.003, 2.427) (3, 503, 1109)

J3 HB (0.999, 1.008, 2.118) (5, 600, 600)

EB (0.999, 1.043, 3.003) (1, 600, 600)

Improving Lives Through Research 19
Example output from the small domain hierarchical Bayes

multivariate model for employee compensation components*

◮ J1/2/3-type model (m = 16, 107; p = 18)

◮ 16,107 survey estimates and in-sample model predictions, and 556,221

not-in-sample model predictions

* Erciulescu and Opsomer, 2019

Improving Lives Through Research 20
Final thoughts...

◮ a flavor of innovation in official statistics programs

◮ existing tools and potential for development of novel ones

◮ software: more than just R JAGS and R STAN

◮ alternative sampling methods

◮ practical challenges: time + number of MCMC samples, large

number of domains (in-sample and not-in-sample), variable

selection, prior distributions for the linking model variance-covariance

components, autocorrelation, cross-validation, visualization, storage

◮ beyond small area estimation - for example, bridging models for the

Association of Fish and Wildlife Agencies (Erciulescu et al. 2019b)

Improving Lives Through Research 21
Selected references

Bell, W. R., Basel W. W. , and Maples J. J. (2016). “An Overview of the U.S. Census Bureau’s Small Area

Income and Poverty Estimates Program.” Analysis of Poverty Data by Small Area Estimation, M. Pratesi

(Ed.). West Sussex: Wiley Sons, Inc., 349-77.

Casas Cordero Valencia, C., Encina J., and Lahiri P. (2016). “Poverty Mapping in Chilean Comunas.” Analysis of

Poverty Data by Small Area Estimation, M. Pratesi (Ed.). West Sussex: Wiley Sons, Inc., 379-403.

Elbers, C., Lanjouw, J. and Lanjouw, P. (2003). “Micro-Level Estimation of Poverty and Inequality.”

Econometrica, 71, 1, 355-364. http://www.jstor.org/stable/3082050.

Erciulescu, A. L. (2018). “Transparent and Reproducible Research in Agricultural Official Statistics.” Government

Advances in Statistical Programming Workshop. http://washstat.org/presentations/20181024/Erciulescu.pdf.

Erciulescu, A. L., Berg E. J., Cecere W., and Ghosh, M. (2018). “A Bivariate Hierarchical Bayesian Model for

Estimating Cropland Cash Rental Rates at the County Level.” Survey Methodology, Statistics Canada,

Catalogue No. 12-001-X, 45, 2.

https://www150.statcan.gc.ca/n1/pub/12-001-x/2019002/article/00001-eng.htm.

Erciulescu, A. L., Cruze N. B., and Nandram, B. (2019a). “Model-Based County Level Crop Estimates

Incorporating Auxiliary Sources of Information.” Journal of the Royal Statistical Society: Series A (Statistics in

Society), 182, 1, 283-303. doi:10.1111/rssa.12390.

Erciulescu, A.L., and Opsomer, J.D. (2019). “A model-based approach to predict employee compensation

components.” Joint Statistical Meetings, Denver, CO, July 2019.

https://ww2.amstat.org/meetings/jsm/2019/onlineprogram/AbstractDetails.cfm?abstractid=300325.

Erciulescu, A.L., Opsomer, J.D., Breidt, J., and Morganstein, D. (2019b). “The data science challenge of bridging

two complex surveys.” New England Statistics Symposium.

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 (366a). Taylor Francis:

269-77. doi:10.1080/01621459.1979.10482505.

Krenzke T., Mohadjer L., Li J., Erciulescu A. L., Fay R., Ren W., Van de Kreckhove W., Li L., Rao J. N. K.

(forthcoming) “Program for the International Assessment of Adult Competencies. State and County Indirect

Estimation Methodology.” United States Department of Education. NCES2019012.

Lewandowski, D., Kurowicka, D., Joe, H. (2009). “Generating random correlation matrices based on vines and

extended onion method.” Journal of Multivariate Analysis, 100, 9, 1989-2001. doi: 10.1016/j.jmva.2009.04.008.

Improving Lives Through Research 22
Thank you!

AndreeaErciulescu@westat.com

Improving Lives Through Research 23