Introduction to STDGLM

Introduction

The STDGLM package provides a framework for fitting spatio-temporal dynamic generalized linear models. These models are useful for analyzing data that varies over both space and time, allowing for the incorporation of spatial and temporal dependencies in the modeling process. The package provides functions for fitting these models, as well as tools for visualizing and interpreting the results.

Installation

You can install the package from GitHub using the following command:

if (!requireNamespace("devtools", quietly = TRUE)) {
  install.packages("devtools")
  require(devtools)
}
devtools::install_github("czaccard/STDGLM")

Run the following command to load the package:

library(STDGLM)

Quick Usage Example

data(ApuliaAQ)
p = length(unique(ApuliaAQ$AirQualityStation)) # 51
t = length(unique(ApuliaAQ$time))              # 365

# distance matrix
W = as.matrix(dist(cbind(ApuliaAQ$Longitude[1:p], ApuliaAQ$Latitude[1:p])))

# response variable: temperature
y = matrix(ApuliaAQ$CL_t2m, p, t)
# covariates (intercept + altitude)
X = array(1, dim = c(p, t, 2))
X[,,2] = matrix(ApuliaAQ$Altitude, p, t)

mod <- stdglm(y=y, X=X, W=W)

Detailed Explanation on supported STDGLMs

Let $p$ denote the number of spatial units (either georeferenced locations or areal units) where data is collected, and let $T$ denote the number of time points. As for the current version of the package (0.0.0.9000), only Gaussian outcomes are supported. Specifically only dynamic linear models (DLMs) with the following observation equation can be handled: $$y_{it} = \boldsymbol{x}_{it}' \boldsymbol{\beta}_{it} + \boldsymbol{z}_{it}' \boldsymbol{\gamma} + \epsilon_{it}, \quad \epsilon_{it} \sim N(0, \sigma_\epsilon^2)$$ where:

• $y_{it}$ is the response variable at spatial unit $i=1, \dots, p$ and time $t=1, \dots, T$,

• $\boldsymbol{x}_{it} = (x_{1,it}, \dots, x_{J,it})'$ is a $J$-dimensional ($J\ge 1$) vector of covariates at spatial unit $i$ at time $t$ (an intercept may or may not be included here),

• $\boldsymbol{\beta}_{it} = (\beta_{1,it}, \dots, \beta_{J,it})'$ is the state vector at time $t$ at spatial unit $i$,

• $\boldsymbol{z}_{it}$ is a $q$-dimensional vector of covariates whose effects are constant (an intercept may or may not be included here),

• $\boldsymbol{\gamma}$ is a vector of non-varying coefficients,

• $\epsilon_{it}$ is the observation error at time $t$ at spatial unit $i$.

The evolution of the state vector is described by the following state equation, which accounts for spatial correlations in the state vector: $$\boldsymbol{\beta}_{j, t} = \boldsymbol{F}_{j,t} \boldsymbol{\beta}_{j, t-1} + \boldsymbol{\eta}_{j,t}, \quad \boldsymbol{\eta}_{j,t} \sim N_p(0, \boldsymbol{\Sigma}_{\eta, j}), \quad j=1, \dots, J$$

where:

• $\boldsymbol{\beta}_{j,t} = (\beta_{j,1t}, \dots, \beta_{j,pt})'$ is the state vector related to the $j$-th covariate $x_{j,t} = (x_{j,1t}, \dots, x_{j,pt})'$, at time $t$ for all spatial units,

• $\boldsymbol{F}_{j,t} = \phi_j^{(\mathsf{T})} \boldsymbol{I}_p$ is a transition matrix,

• $\boldsymbol{\Sigma}_{\eta,j}$ is the covariance matrix of the state evolution error $\boldsymbol{\eta}_{j,t}$.

The transition matrix $\boldsymbol{F}_{j,t}, j=1, \dots, J,$ is assumed to be a scalar multiple of the identity matrix. The parameter $\phi_j^{(\mathsf{T})}$ controls the temporal autocorrelation of the state vector.

The state evolution covariance matrix $\boldsymbol{\Sigma}_{\eta,j}$ can be structured to reflect spatial relationships, e.g. by assuming an exponential covariance function if the data are point-referenced: $$(\boldsymbol{\Sigma}_{\eta,j})_{i\ell} = \rho_{j,2} \exp\left(-\frac{d_{i\ell}}{\rho_{j,2}}\right), \quad d_{i\ell} = \| \boldsymbol{s}_i - \boldsymbol{s}_\ell \|$$ where $\rho_{j,1}$ is the partial sill, $\rho_{j,2}$ is the range parameter, and $ d_{i}$ is the Euclidean distance between locations $\boldsymbol{s}_i$ and $\boldsymbol{s}_\ell$. At the moment, this is the only supported covariance structure for point-referenced data.

In this case, the evolution error $\boldsymbol{\eta}_{j,t}$ is assumed to be a zero-mean Gaussian process with exponential covariance matrix parameteterized by $\rho_{j,1}$ and $\rho_{j,2}$, which we will denote as $\boldsymbol{\eta}_{j,t} \sim GP(0, \rho_{j,1}, \rho_{j,2}; exp)$.

If the data are areal, a proper conditional autoregressive (PCAR) covariance structure is assumed: $$\boldsymbol{\Sigma}_{\eta,j} = \rho_{j,1} \left( \boldsymbol{D}_w - \rho_{j,2} \boldsymbol{W} \right)^{-1}$$ where $\boldsymbol{W}$ is a binary adjacency matrix, $\boldsymbol{D}_w$ is a diagonal matrix with row sums of $\boldsymbol{W}$ on the diagonal, and $\rho_{j,1}$ and $\rho_{j,2}$ are the conditional variance and autocorrelation parameters, respectively.

In this case, the evolution error $\boldsymbol{\eta}_{j,t}$ follows a zero-mean PCAR process, and we will denote this as $\boldsymbol{\eta}_{j,t} \sim PCAR(0, \rho_{j,1}, \rho_{j,2})$.

ANOVA Decomposition of the State Vector

The function stdglm allows for the decomposition of the state vector into components that can be interpreted as contributions from different sources of variability. Dropping the subscript $j$ for the sake of simplicity, the state vector $\beta_{it}$ is decomposed as follows: $$\beta_{it} = \overline{\beta} + \beta_{i}^{(\mathsf{S})} + \beta_{t}^{(\mathsf{T})} + \beta_{it}^{(\mathsf{ST})}$$ where:

• $\overline{\beta}$ is the overall mean effect,

• $\beta_{i}^{(\mathsf{S})}$ is the spatial effect at spatial unit $i$,

• $\beta_{t}^{(\mathsf{T})}$ is the temporal effect at time $t$,

• $\beta_{it}^{(\mathsf{ST})}$ is the interaction effect between space and time at spatial unit $i$ and time $t$.

Bayesian Hierarchical Structure

The Bayesian model is as follows, for $t=1, \dots, T$ and $j=1, \dots, J$: $$\begin{align*} y_{it} &\sim N(\boldsymbol{x}_{it}' \boldsymbol{\beta}_{it} + \boldsymbol{z}_{it}' \boldsymbol{\gamma}, \sigma_\epsilon^2) \\ \boldsymbol{\beta}_{j,t} &= \boldsymbol{1}_p \overline{\beta}_j + \boldsymbol{\beta}_j^{(\mathsf{S})} + \boldsymbol{1}_p \beta_{j,t}^{(\mathsf{T})} + \boldsymbol{\beta}_{j,t}^{(\mathsf{ST})} \\ \boldsymbol{\beta_j}^{(\mathsf{S})} &= (\beta_{j,1}^{(\mathsf{S})}, \dots, \beta_{j,p}^{(\mathsf{S})})' \sim GP(0, \rho_{j,1}^{(\mathsf{S})}, \rho_{j,2}^{(\mathsf{S})}; exp) \mbox{ or } PCAR(0, \rho_{j,1}^{(\mathsf{S})}, \rho_{j,2}^{(\mathsf{S})}) \\ \beta_{j,t}^{(\mathsf{T})} &\sim N(\phi_j^{(\mathsf{T})} \beta_{j,t-1}^{(\mathsf{T})}, V_{\beta,j}^{(\mathsf{T})}) \\ \boldsymbol{\beta}_{j,t}^{(\mathsf{ST})} &= (\beta_{j,1t}^{(\mathsf{ST})}, \dots, \beta_{j,pt}^{(\mathsf{ST})})' \sim GP(\phi_j^{(\mathsf{ST})} \boldsymbol{\beta}_{j,t-1}^{(\mathsf{ST})}, \rho_{j,1}^{(\mathsf{ST})}, \rho_{j,2}^{(\mathsf{ST})}; exp) \mbox{ or } PCAR(\phi_j^{(\mathsf{ST})} \boldsymbol{\beta}_{j,t-1}^{(\mathsf{ST})}, \rho_{j,1}^{(\mathsf{ST})}, \rho_{j,2}^{(\mathsf{ST})}) \end{align*}$$

The model is completed with the following priors (again, dropping the subscript $j$, since these are common across $j$): $$\begin{align*} \overline{\beta} &\sim N(0, V_\gamma) \\ \beta_{0}^{(\mathsf{T})} &\sim N(0, V_{\beta_0}) \\ \boldsymbol{\beta}_{0}^{(\mathsf{ST})} &\sim N_p(0, V_{\beta_0} \boldsymbol{I}_p) \\ \boldsymbol{\gamma} &\sim N_q(0, V_\gamma) \\ \sigma_\epsilon^2 &\sim IG(a_\epsilon, b_\epsilon) \\ \rho_1^{(\mathsf{S})} &\sim IG(0.01, 0.01) \\ \rho_2^{(\mathsf{S})} &\sim U(a_\rho, b_\rho) \\ \rho_1^{(\mathsf{ST})} &\sim IG(0.01, 0.01) \\ \rho_2^{(\mathsf{ST})} &\sim U(a_\rho, b_\rho) \\ \phi^{(\mathsf{T})} &\sim TN_{(-1, 1)}(0, 1) \\ \phi^{(\mathsf{ST})} &\sim TN_{(-1, 1)}(0, 1) \\ V_\beta^{(\mathsf{T})} &\sim IG(a^{(\mathsf{T})}, b^{(\mathsf{T})}) \end{align*}$$ where $TN_{(q, r)}(\mu, \sigma^2)$ denotes a normal distribution with mean $\mu$ and variance $\sigma^2$ truncated to the interval $(q, r)$. The hyperparameters $a_\rho$ and $b_\rho$ depend on the type of spatial data. If the data are point-referenced, they are set based on the minimum and maximum distances between points, respectively. If the data are areal, $a_\rho= 0.1$ and $b_\rho \rightarrow 1$.

Note that the spacetime-varying coefficients are assumed independent a priori across $j=1,\dots,J$.

Efficient Inference and Identifiability

To build an efficient sampler, the algorithm proposed by Chan and Jeliazkov (2009) is used in conjuction with sparse matrix techniques.

To make the model identifiable, some constraints are imposed on the varying parameters at each MCMC iteration:

• Set $\sum_{t=1}^{T} \beta_{t}^{(\mathsf{T})} = 0$.

• Set $\sum_{i=1}^{p} \beta_{i}^{(\mathsf{S})} = 0$.

• Set $\sum_{i=1}^{p} \beta_{i,t}^{(\mathsf{ST})} = 0$ for each $t=1,\dots,T$.

• Set $\sum_{t=1}^{T} \beta_{i,t}^{(\mathsf{ST})} = 0$ for each $i=1,\dots,p$.

References

Chan, Joshua CC, and Ivan Jeliazkov. 2009. “Efficient Simulation and Integrated Likelihood Estimation in State Space Models.” International Journal of Mathematical Modelling and Numerical Optimisation 1 (1-2): 101–20.