dstm_ide: Integrodifference equation (IDE) model

Description

dstm_ide fits a type of dynamic spatio-temporal model called an integrodifference equation (IDE) model. It estimates a redistribution kernel---a probability distribution controlling diffusion across time and space. Currently, only Gaussian redistribution kernels are supported.

The process model is decomposed with an orthonormal basis function expansion (a Fourier series). It can then be estimated as a special case of a dynamic linear model (DLM), using the forward filtering backward sampling algorithm to estimate the state vector. The kernel parameters are estimated with a random walk Metropolis-Hastings update. The other parameters are estimated with conditionally conjugate updates.

Usage

dstm_ide(
  Y,
  locs = NULL,
  knot_locs = NULL,
  proc_error = "IW",
  J = 1L,
  n_samples = 1L,
  sample_sigma2 = TRUE,
  verbose = FALSE,
  params = NULL
)

Arguments

(numeric matrix) S by T data matrix containing response variable at S spatial locations and T time points. The t-th column (NOT row) corresponds to the t-th observation vector.

locs

(numeric matrix) S by 2 matrix containing the spatial locations of the observed data. The rows of `locs` correspond with the rows of `Y`.

knot_locs

(integer or numeric matrix) Knot locations for the spatially varying IDE model. The kernel parameters are estimated at these locations and then mapped to the spatial locations of the observed data via process convolution. If an integer is provided, then the knots are located on an equally spaced grid with dimension (`knot_locs`, `knot_locs`). If a matrix is provided, then each row of the matrix corresponds to a knot location. If NULL, then the standard (spatially constant) IDE is fit.

proc_error

(character string) Process error: "IW" (inverse-Wishart) or "Discount" (discount factor). "IW" is recommended because it is more computationally stable.

(integer) Extent of the Fourier approximation. The size of the state space is (2*`J` + 1)^2.

n_samples

(integer) Number of posterior samples to draw.

sample_sigma2

(logical) Whether to sample the variance of the iid observation error.

verbose

(logical) Whether to print additional information; e.g., iteration in sampling algorithm.

params

(list) List of hyperparameter values; see details.

Details

This section explains how to specify custom hyperparameters using the `params` argument. For each distribution referenced below, we use the scale parameterization found on the distribution's Wikipedia page. You may specify the following as named elements of the `params` list:

m_0: (numeric vector) The prior mean of the state vector at time zero (\(\theta_0\)).

C_0: (numeric matrix) The prior variance-covariance matrix of the state vector at time zero (\(\theta_0\)).

alpha_sigma2, beta_sigma2: (numeric scalars) The inverse-Gamma parameters (scale parameterization) of the prior distribution on the observation error (\(\sigma^2\)).

sigma2: (numeric scalar) The value to use for the observation error (\(\sigma^2\)) if `sample_sigma2` = FALSE.

alpha_lambda, beta_lambda: (numeric scalars) The inverse-Gamma parameters (scale parameterization) of the prior distribution on \(\lambda = (1 - \delta) / \delta\), where \(\delta\) is the discount factor.

scale_W: (numeric matrix) The scale matrix for the inverse-Wishart prior distribution on the variance-covariance matrix of the process error (`W`).

df_W: (numeric scalar) The degees of freedom for the inverse-Wishart prior distribution on the variance-covariance matrix of the process error (`W`).

L: (numeric scalar) The period of the Fourier series approximation. The spatial locations and knot locations are rescaled to range from -`L`/4 to `L`/4 because the Fourier decomposition assumes that the spatial surface is periodic. Regardless of the value of `L`, kernel parameter estimates are back-transformed to the original scale.

smoothing: (numeric scalar) Controls the degree of smoothing in the process convolution for models with spatially varying kernel parameters. The values in the process convolution matrix are proportional to exp(d/`smoothing`) where d is the distance between spatial locations before rescaling with `L`

mean_mu_kernel: (numeric vector) The mean of the normal prior distribution on `mu_kernel`, the mean of the redistribution kernel. In the spatially varying case, the prior distribution for `mu_kernel` is assumed to be the same at every knot location.

var_mu_kernel: (numeric matrix) The variance of the normal prior distribution on `mu_kernel`, the mean of the redistribution kernel.

scale_Sigma_kernel: (numeric matrix) The scale matrix for the inverse-Wishart prior distribution on `Sigma_kernel`, the variance-covariance matrix of the redistribution kernel.

df_Sigma_kernel: (numeric scalar) The degrees of freedom for the inverse-Wishart prior distribution on `Sigma_kernel`, the variance-covariance matrix of the redistribution kernel.

proposal_factor_mu: (numeric scalar) Controls the variance of the proposal distribution for `mu_kernel`. The proposals have a variance of `proposal_factor_mu`^2 * `var_mu_kernel`. `proposal_factor_mu` must generally be set lower for spatially varying models.

proposal_factor_Sigma: (numeric scalar) Controls the variance of the proposal distribution for `Sigma_kernel`. As is the case with `proposal_factor_mu`, a higher value corresponds to a higher variance. The degrees of freedom for the proposal distribution for `Sigma_kernel` is ncol(`locs`) + `df_Sigma_kernel` / `proposal_factor_Sigma`. `proposal_factor_Sigma` must generally be set lower for spatially varying models.

kernel_samples_per_iter: (numeric scalar) Number of times to update the kernel parameters per iteration of the sampling loop.

References

Cressie, N., and Wikle, C. K. (2011), Statistics for spatio-temporal data, John Wiley and Sons, New York, ISBN:978-0471692744.

Fruhwirth-Schnatter, S. (1994), “Data Augmentation and Dynamic Linear Models,” Journal of Time Series Analysis, 15, 183<U+2013>202, <doi:10.1111/j.1467-9892.1994.tb00184.x>.

Petris, G., Petrone, S., and Campagnoli, P. (2009), Dynamic Linear Models with R, useR!, Springer-Verlag, New York, ISBN:978-0387772370, <doi:10.1007/b135794>.

Wikle, C. K., and Cressie, N. (1999), “A dimension-reduced approach to space-time Kalman filtering,” Biometrika, 86, 815<U+2013>829, <https://www.jstor.org/stable/2673587>.

Wikle, C. K. (2002), “A kernel-based spectral model for non-Gaussian spatio-temporal processes,” Statistical Modelling, 2, 299<U+2013>314, <doi:10.1191/1471082x02st036oa>.

Examples

Run this code

# NOT RUN {
# Load example data
data("ide_standard", "ide_spatially_varying", "ide_locations")

# Basic IDE model with one kernel
mod <- dstm_ide(ide_standard, ide_locations)
predict(mod)
summary(mod)

# IDE model with spatially varying kernel
dstm_ide(ide_spatially_varying, ide_locations, knot_locs=4)

# Fix sigma2
dstm_ide(ide_standard, ide_locations, 
         sample_sigma2=FALSE, params=list(sigma2=1))

# Set proposal scaling factors, number of kernel updates per iteration,
# and prior distribution on kernel parameters
dstm_ide(ide_standard, ide_locations, 
         params=list(proposal_factor_mu=2, proposal_factor_Sigma=3,
                     kernel_updates_per_iter=2,
                     scale_Sigma_kernel=diag(2), df_Sigma_kernel=100,
                     mean_mu_kernel=c(0.2, 0.4), var_mu_kernel=diag(2)))
                     
# Set priors on state vector, process error, and observation error
J <- 1
P <- (2*J + 1)^2
dstm_ide(ide_standard, ide_locations,
         params=list(m_0=rep(1, P), C_0=diag(0.01, P),
                     alpha_sigma2=20, beta_sigma2=20,
                     scale_W=diag(P), df_W=100))
# }