Many data sources report related variables of interest that are also referenced over geographic regions
and time; however, there are relatively few general statistical methods that one can readily
use that incorporate these multivariate spatio-temporal dependencies. Additionally, many multivariate
spatio-temporal areal datasets are extremely high-dimensional, which leads to practical
issues when formulating statistical models. For example, we analyze Quarterly Workforce Indicators
(QWI) published by the US Census Bureau’s Longitudinal Employer-Household Dynamics
(LEHD) program. QWIs are available by different variables, regions, and time points, resulting
in millions of tabulations. Despite their already expansive coverage, by adopting a fully Bayesian
framework, the scope of the QWIs can be extended to provide estimates of missing values along
with associated measures of uncertainty. Motivated by the LEHD, and other applications in federal
statistics, we introduce the multivariate spatio-temporal mixed effects model (MSTM), which can
be used to efficiently model high-dimensional multivariate spatio-temporal areal datasets. The proposed
MSTM extends the notion of Moran’s I basis functions to the multivariate spatio-temporal
setting. This extension leads to several methodological contributions including extremely effective
dimension reduction, a dynamic linear model for multivariate spatio-temporal areal processes, and
the reduction of a high-dimensional parameter space using a novel parameter model.