glmmfields: Fit a spatiotemporal random fields GLMM

Description Usage Arguments Details Examples

Description

Fit a spatiotemporal random fields model that optionally uses the MVT distribution instead of a MVN distribution to allow for spatial extremes through time. It is also possible to fit a spatial random fields model without a time component.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
glmmfields(formula, data, lon, lat, time = NULL, nknots = 15L,
  prior_gp_theta = half_t(3, 0, 5), prior_gp_sigma = half_t(3, 0, 5),
  prior_sigma = half_t(3, 0, 5), prior_rw_sigma = half_t(3, 0, 5),
  prior_intercept = student_t(3, 0, 10), prior_beta = student_t(3, 0,
  3), prior_phi = student_t(1000, 0, 0.5), fixed_df_value = 1000,
  fixed_phi_value = 0, estimate_df = FALSE, estimate_ar = FALSE,
  family = gaussian(link = "identity"),
  covariance = c("squared-exponential", "exponential", "matern"),
  matern_kappa = 0.5, algorithm = c("sampling", "meanfield"),
  year_re = FALSE, nb_lower_truncation = 0,
  control = list(adapt_delta = 0.9), save_log_lik = FALSE,
  df_lower_bound = 2, cluster = c("pam", "kmeans"), ...)

Arguments

formula

The model formula.

data

A data frame.

lon

A character object giving the name of the longitude column.

lat

A character object giving the name of the latitude column.

time

A character object giving the name of the time column. Leave as NULL to fit a spatial GLMM without a time element.

nknots

The number of knots to use in the predictive process model. Smaller values will be faster but may not adequately represent the shape of the spatial pattern.

prior_gp_theta

The prior on the Gaussian Process scale parameter. Must be declared with half_t(). Here, and throughout, priors that are normal or half-normal can be implemented by setting the first parameter in the half-t or student-t distribution to a large value. E.g. something greater than 100.

prior_gp_sigma

The prior on the Gaussian Process eta parameter. Must be declared with half_t().

prior_sigma

The prior on the observation process scale parameter. Must be declared with half_t(). This acts as a substitute for the scale parameter in whatever observation distribution is being used. E.g. the CV for the Gamma or the dispersion parameter for the negative binomial.

prior_rw_sigma

The prior on the standard deviation parameter of the random walk process (if specified). Must be declared with half_t().

prior_intercept

The prior on the intercept parameter. Must be declared with student_t().

prior_beta

The prior on the slope parameters (if any). Must be declared with student_t().

prior_phi

The prior on the AR parameter. Must be declared with student_t().

fixed_df_value

The fixed value for the student-t degrees of freedom parameter if the degrees of freedom parameter is fixed in the MVT. If the degrees of freedom parameter is estimated then this argument is ignored. Must be 1 or greater. Very large values (e.g. the default value) approximate the normal distribution. If the value is >=1000 then a true MVN distribution will be fit.

fixed_phi_value

The fixed value for temporal autoregressive parameter, between random fields at time(t) and time(t-1). If the phi parameter is estimated then this argument is ignored.

estimate_df

Logical: should the degrees of freedom parameter be estimated?

estimate_ar

Logical: should the AR (autoregressive) parameter be estimated? Here, this refers to a autoregressive process in the evolution of the spatial field through time.

family

Family object describing the observation model. Note that only one link is implemented for each distribution. Gamma, negative binomial (specified via nbinom2() as nbinom2(link = "log"), and Poisson must have a log link. Binomial must have a logit link. Also implemented is the lognormal (specified via lognormal() as lognormal(link = "log"). Besides the negative binomial and lognormal, other families are specified as shown in family.

covariance

The covariance function of the Gaussian Process. One of "squared-exponential", "exponential", or "matern".

matern_kappa

Optional parameter for the Matern covariance function. Optional values are 1.5 or 2.5. Values of 0.5 are equivalent to exponential.

algorithm

Character object describing whether the model should be fit with full NUTS MCMC or via the variational inference mean-field approach. See rstan::vb(). Note that the variational inference approach should not be trusted for final inference and is much more likely to give incorrect inference than MCMC.

year_re

Logical: estimate a random walk for the time variable? If TRUE, then no fixed effects (B coefficients) will be estimated. In this case, prior_intercept will be used as the prior for the initial value in time.

nb_lower_truncation

For NB2 only: lower truncation value. E.g. 0 for no truncation, 1 for 1 and all values above. Note that estimation is likely to be considerably slower with lower truncation because the sampling is not vectorized. Also note that the log likelihood values returned for estimating quantities like LOOIC will not be correct if lower truncation is implemented.

control

List to pass to rstan::sampling(). For example, increase adapt_delta if there are warnings about divergent transitions: control = list(adapt_delta = 0.99). By default, glmmfields sets adapt_delta = 0.9.

save_log_lik

Logical: should the log likelihood for each data point be saved so that information criteria such as LOOIC or WAIC can be calculated? Defaults to FALSE so that the size of model objects is smaller.

df_lower_bound

The lower bound on the degrees of freedom parameter. Values that are too low, e.g. below 2 or 3, it might affect chain convergence. Defaults to 2.

cluster

The type of clustering algorithm used to determine the knot locations. "pam" = cluster::pam(). The "kmeans" algorithm will be faster on larger datasets.

...

Any other arguments to pass to rstan::sampling().

Details

Note that there is no guarantee that the default priors are reasonable for your data. Also, there is no guarantee the default priors will remain the same in future versions. Therefore it is important that you specify any priors that are used in your model, even if they replicate the defaults in the package. It is particularly important that you consider that prior on gp_theta since it depends on the distance between your location points. You may need to scale your coordinate units so they are on a ballpark range of 1-10 by, say, dividing the coordinates (say in UTMs) by several order of magnitude.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# Spatiotemporal example:
set.seed(1)
s <- sim_glmmfields(n_draws = 12, n_knots = 12, gp_theta = 1.5,
gp_sigma = 0.2, sd_obs = 0.2)
print(s$plot)
# options(mc.cores = parallel::detectCores()) # for parallel processing
# should use 4 or more chains for real model fits
m <- glmmfields(y ~ 0, time = "time",
 lat = "lat", lon = "lon", data = s$dat,
 nknots = 12, iter = 1000, chains = 2, seed = 1)

# Spatial example (with covariates) from the vignette and customizing
# some priors:
set.seed(1)
N <- 100 # number of data points
temperature <- rnorm(N, 0, 1) # simulated temperature data
X <- cbind(1, temperature) # design matrix
s <- sim_glmmfields(n_draws = 1, gp_theta = 1.2, n_data_points = N,
  gp_sigma = 0.3, sd_obs = 0.1, n_knots = 12, obs_error = "gamma",
  covariance = "squared-exponential", X = X,
  B = c(0.5, 0.2)) # B represents our intercept and slope
d <- s$dat
d$temperature <- temperature
library(ggplot2)
ggplot(s$dat, aes(lon, lat, colour = y)) +
  viridis::scale_colour_viridis() +
  geom_point(size = 3)
m_spatial <- glmmfields(y ~ temperature, data = d, family = Gamma(link = "log"),
  lat = "lat", lon = "lon", nknots = 12, iter = 2000, chains = 2,
  prior_beta = student_t(100, 0, 1), prior_intercept = student_t(100, 0, 5),
  control = list(adapt_delta = 0.95))

glmmfields documentation built on May 18, 2019, 9 a.m.