Spatio-temporal change of support modeling with R

Abstract

Spatio-temporal change of support methods are designed for statistical analysis on spatial and temporal domains which can differ from those of the observed data. Previous work introduced a parsimonious class of Bayesian hierarchical spatio-temporal models, which we refer to as STCOS, for the case of Gaussian outcomes. Application of STCOS methodology from this literature requires a level of proficiency with spatio-temporal methods and statistical computing which may be a hurdle for potential users. The present work seeks to bridge this gap by guiding readers through STCOS computations. We focus on the R computing environment because of its popularity, free availability, and high quality contributed packages. The stcos package is introduced to facilitate computations for the STCOS model. A motivating application is the American Community Survey (ACS), an ongoing survey administered by the U.S. Census Bureau that measures key socioeconomic and demographic variables for various populations in the United States. The STCOS methodology offers a principled approach to compute model-based estimates and associated measures of uncertainty for ACS variables on customized geographies and/or time periods. We present a detailed case study with ACS data as a guide for change of support analysis in R, and as a foundation which can be customized to other applications.

Computational details and proofs

We will make use of the following well-known property in several places.

Property 1

If ${\varvec{A}} \in \mathbb {R}^{m \times k}$, ${\varvec{B}} \in \mathbb {R}^{k \times l}$, ${\varvec{C}} \in \mathbb {R}^{l \times n}$, then $\text {vec}({\varvec{A}} {\varvec{B}} {\varvec{C}}) = ({\varvec{C}}^\top \otimes {\varvec{A}}) \text {vec}({\varvec{B}})$.

The following proposition gives the explicit solution to the minimization problem stated in (3.10). Bradley et al. (2015a) considers a similar problem featuring a more general objective function but assuming that the columns of ${\varvec{S}}$ are orthonormal. Higham (1988) gives a general discussion of problems involving Frobenius and 2-norm distance minimization.

Proposition 1

(Frobenius Norm Minimization) Suppose ${\varvec{S}} \in \mathbb {R}^{n \times r}$ has rank r and ${\varvec{\varSigma }} \in \mathbb {R}^{n \times n}$ is positive definite. The minimizer ${\varvec{X}} \in \mathbb {R}^{r \times r}$ of $\Vert {\varvec{\varSigma }} - {\varvec{S}} {\varvec{X}} {\varvec{S}}^\top \Vert _\text {F}$ is ${\varvec{X}} = ({\varvec{S}}^\top {\varvec{S}})^{-1} {\varvec{S}}^\top {\varvec{\varSigma }} {\varvec{S}} ({\varvec{S}}^\top {\varvec{S}})^{-1}$.

Proof

Using Property 1, we have

$$\begin{aligned} \Vert {\varvec{\varSigma }} - {\varvec{S}} {\varvec{X}} {\varvec{S}}^\top \Vert _{\text {F}}^2&= \text {vec}\left[ {\varvec{\varSigma }} - {\varvec{S}} {\varvec{X}} {\varvec{S}}^\top \right] ^\top \text {vec}\left[ {\varvec{\varSigma }} - {\varvec{S}} {\varvec{X}} {\varvec{S}}^\top \right] \nonumber \\&= \left[ \text {vec}({\varvec{\varSigma }}) - \text {vec}({\varvec{S}} {\varvec{X}} {\varvec{S}}^\top ) \right] ^\top \left[ \text {vec}({\varvec{\varSigma }}) - \text {vec}({\varvec{S}} {\varvec{X}} {\varvec{S}}^\top ) \right] \nonumber \\&= \left[ \text {vec}({\varvec{\varSigma }}) - ({\varvec{S}} \otimes {\varvec{S}}) \text {vec}({\varvec{X}}) \right] ^\top \left[ \text {vec}({\varvec{\varSigma }}) - ({\varvec{S}} \otimes {\varvec{S}}) \text {vec}({\varvec{X}}) \right] \nonumber \\&= \Vert \text {vec}({\varvec{\varSigma }}) - ({\varvec{S}} \otimes {\varvec{S}}) \text {vec}({\varvec{X}}) \Vert _2^2, \end{aligned}$$

(A.1)

where the norm on the last line is the usual 2-norm on $\mathbb {R}^{n^2}$. We recognize the expression in (A.1) as a standard least squares minimization whose solution is

$$\begin{aligned} \text {vec}({\varvec{X}})&= [({\varvec{S}} \otimes {\varvec{S}})^\top ({\varvec{S}} \otimes {\varvec{S}})]^{-1} ({\varvec{S}} \otimes {\varvec{S}})^\top \text {vec}({\varvec{\varSigma }}) \\&= [({\varvec{S}}^\top \otimes {\varvec{S}}^\top ) ({\varvec{S}} \otimes {\varvec{S}})]^{-1} ({\varvec{S}}^\top \otimes {\varvec{S}}^\top ) \text {vec}({\varvec{\varSigma }}) \\&= [ {\varvec{S}}^\top {\varvec{S}} \otimes {\varvec{S}}^\top {\varvec{S}} ]^{-1} \text {vec}({\varvec{S}}^\top {\varvec{\varSigma }} {\varvec{S}}) \\&= [({\varvec{S}}^\top {\varvec{S}})^{-1} \otimes ({\varvec{S}}^\top {\varvec{S}})^{-1}] \text {vec}({\varvec{S}}^\top {\varvec{\varSigma }} {\varvec{S}}) \\&= \text {vec}\left[ ({\varvec{S}}^\top {\varvec{S}})^{-1} {\varvec{S}}^\top {\varvec{\varSigma }} {\varvec{S}} ({\varvec{S}}^\top {\varvec{S}})^{-1} \right] . \end{aligned}$$

Therefore, the minimizer is ${\varvec{X}} = ({\varvec{S}}^\top {\varvec{S}})^{-1} {\varvec{S}}^\top {\varvec{\varSigma }} {\varvec{S}} ({\varvec{S}}^\top {\varvec{S}})^{-1}$, as desired. $\square$

Remark 1

(MLE Computation) To compute the MLE for the STCOS model, we first note that the likelihood, excluding the parameter model, is

$$\begin{aligned} f({\varvec{Z}} \mid {\varvec{\mu }}_B, \sigma _K^2, \sigma _\xi ^2)&= \int \phi ({\varvec{Z}} \mid {\varvec{H}} {\varvec{\mu }}_B + {\varvec{S}} {\varvec{\eta }}, \sigma _\xi ^2 {\varvec{I}} + {\varvec{V}}) \cdot \phi ({\varvec{\eta }} \mid {\varvec{0}}, \sigma _K^2 {\varvec{K}}) d {\varvec{\eta }} \\&= \phi ({\varvec{Z}} \mid {\varvec{H}} {\varvec{\mu }}_B, {\varvec{\varDelta }}) \\&= (2 \pi )^{-N/2} |{\varvec{\varDelta }}|^{-1/2} \exp \left\{ -\frac{1}{2} ({\varvec{Z}} - {\varvec{H}} {\varvec{\mu }}_B)^\top {\varvec{\varDelta }}^{-1} ({\varvec{Z}} - {\varvec{H}} {\varvec{\mu }}_B) \right\} , \end{aligned}$$

where ${\varvec{\varDelta }} = \sigma _\xi ^2 {\varvec{I}} + {\varvec{V}} + \sigma _K^2 {\varvec{S}} {\varvec{K}} {\varvec{S}}^\top$. Given $\sigma _K^2$ and $\sigma _\xi ^2$, the likelihood is maximized by the weighted least squares estimator $\hat{{\varvec{\mu }}}_B = ({\varvec{H}}^\top {\varvec{\varDelta }}^{-1} {\varvec{H}})^{-1} {\varvec{H}}^\top {\varvec{\varDelta }}^{-1} {\varvec{Z}}$. To estimate the unknown $\sigma _K^2$ and $\sigma _\xi ^2$, we carry out numerical maximization on the partially maximized log-likelihood

$$\begin{aligned} \ell (\sigma _K^2, \sigma _\xi ^2) = -\frac{N}{2} \log (2 \pi ) -\frac{1}{2} \log |{\varvec{\varDelta }}| -\frac{1}{2} ({\varvec{Z}} - {\varvec{H}} \hat{{\varvec{\mu }}}_B)^\top {\varvec{\varDelta }}^{-1} ({\varvec{Z}} - {\varvec{H}} \hat{{\varvec{\mu }}}_B). \end{aligned}$$

To enforce the constraints that $\sigma _K^2 > 0$ and $\sigma _\xi ^2 > 0$, we optimize over $(\vartheta _1, \vartheta _2) \in \mathbb {R}^2$ and take $\sigma _K^2 = \exp (\vartheta _1)$, $\sigma _\xi ^2 = \exp (\vartheta _2)$.