correlation: Correlation factor structures in generic model components

Description

Element 'factor' of a model component created using function gen is a formula composed of several possible terms described below. It is used to derive a (typically sparse) precision matrix for a set of coefficients, and possibly a matrix representing a set of linear constraints to be imposed on the coefficient vector.

iid(f): Independent effects corresponding to the levels of factor f.
RW1(f, circular=FALSE, w=NULL): First-order random walk over the levels of factor f. The random walk can be made circular and different (fixed) weights can be attached to the innovations. If specified, w must be a positive numeric vector of length one less than the number of factor levels. For example, if the levels correspond to different times, it would often be reasonable to choose w proportional to the reciprocal time differences. For equidistant times there is generally no need to specify w.
RW2(f): Second-order random walk.
AR1(f, phi, w=NULL, control=NULL): First-order autoregressive correlation structure among the levels of f. Argument phi can be a single numerical value of the autoregressive parameter, or an appropriate prior specification if phi should be inferred. If not supplied, a uniform prior on (-1, 1] is assumed. For irregularly spaced AR(1) processes weights can be specified, in the same way as for RW1.
season(f, period): Dummy seasonal with period period.
spatial(f, graph, snap, queen): CAR spatial correlation. Argument graph can either be an object of (S4) class SpatialPolygonsDataFrame or an object of (S3) class sf. The latter can be obtained, e.g., by reading in a shape file using function st_read. Arguments snap and queen are passed to poly2nb, which computes a neighbours list. Alternatively, a neighbours list object of class nb can be passed directly as argument graph.
spline(f, knots, degree): P-splines, i.e. penalized B-splines structure over the domain of a quantitative variable f. Arguments knots and degree are passed to splineDesign. If knots is a single value it is interpreted as the number of knots, otherwise as a vector of knot positions. By default 40 equally spaced knots are used, and a degree of 3.
custom(f, D=NULL, Q=NULL, R=NULL, derive.constraints=NULL): Either a custom precision or incidence matrix associated with factor f can be passed to argument Q or D. Optionally a constraint matrix can be supplied as R, or constraints can be derived from the null space of the precision matrix by setting derive.constraints=TRUE.

Usage

iid(name)
RW1(name, circular = FALSE, w = NULL)
RW2(name)
AR1(name, phi, w = NULL, control = NULL)
season(name, period)
spatial(
  name,
  graph = NULL,
  snap = sqrt(.Machine$double.eps),
  queen = TRUE,
  poly.df = NULL,
  derive.constraints = FALSE
)
spline(name, knots, degree)
custom(name, D = NULL, Q = NULL, R = NULL, derive.constraints = NULL)

Arguments

name: name of a variable, unquoted.
circular: whether the random walk is circular.
w: a vector of weights.
phi: prior distribution, or fixed value, for an autoregressive parameter. The default is a uniform prior over the interval [-1, 1]. A single numeric value is interpreted as a fixed value, corresponding to a degenerate prior, which can also be specified as pr_fixed(value). Alternatively, link{pr_truncnormal} can be used to specify a truncated normal prior.
control: options for Metropolis-Hastings sampling from the conditional posterior for an autoregressive parameter. These options can be set using function set_MH. Supported proposal types are "TN" and "RWTN". By default an independence truncated normal proposal (type="TN"), or a random walk truncated normal proposal (type="RWTN") with adaptive scale initialised at 0.025 is used, depending on whether the specified random effects' distribution is Gaussian or non-Gaussian.
period: a positive integer specifying the seasonal period.
graph: either a spatial object of class SpatialPolygons, sf, sfc, or a neighbours list of class nb.
snap: passed to poly2nb. Ignored if graph is a neighbours list.
queen: passed to poly2nb. Ignored if graph is a neighbours list.
poly.df: a spatial data frame. DEPRECATED, use argument graph instead.
derive.constraints: whether to derive the constraint matrix for an IGMRF model component numerically from the precision matrix. The use of derive.constraints in function spatial is DEPRECATED, as it is no longer needed.
knots: passed to splineDesign.
degree: passed to splineDesign.
D: custom incidence matrix.
Q: custom precision matrix.
R: custom restriction matrix.

References

B. Allevius (2018). On the precision matrix of an irregularly sampled AR(1) process. arXiv:1801.03791.

H. Rue and L. Held (2005). Gaussian Markov Random Fields. Chapman & Hall/CRC.

Examples

Run this code

# \donttest{
# example of CAR spatial random effects
if (requireNamespace("sf")) {
  # 1. load a shape file of counties in North Carolina
  nc <- sf::st_read(system.file("shape/nc.shp", package="sf"))
  # 2. generate some data according to a model with a few regression
  # effects, as well as spatial random effects
  gd <- generate_data(
    ~ reg(~ AREA + BIR74, prior=pr_normal(precision=1), name="beta") +
      gen(factor = ~ spatial(NAME, graph=nc), name="vs"),
    sigma.mod = pr_invchisq(df=10, scale=1),
    data = nc
  )
  # add the generated target variable and the spatial random effects to the
  # spatial dataframe object
  nc$y <- gd$y
  nc$vs_true <- gd$pars$vs
  # 3. fit a model to the generated data, and see to what extent the
  #    parameters used to generate the data, gd$pars, are reproduced
  sampler <- create_sampler(
    y ~ reg(~ AREA + BIR74, prior=pr_normal(precision=1), name="beta") +
    gen(factor = ~ spatial(NAME, graph=nc), name="vs"),
    data=nc
  )
  sim <- MCMCsim(sampler, store.all=TRUE, n.iter=600, n.chain=2, verbose=FALSE)
  (summ <- summary(sim))
  nc$vs <- summ$vs[, "Mean"]
  plot(nc[c("vs_true", "vs")])
  plot(gd$pars$vs, summ$vs[, "Mean"]); abline(0, 1, col="red")
}
# }

Run the code above in your browser using DataLab