correlation: Correlation factor structures in generic model components

In mcmcsae: Markov Chain Monte Carlo Small Area Estimation

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.

splines(f, knots, degree)

P-splines, i.e. penalised 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
)

splines(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.03791v2.

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

Examples

## Not run: 
# 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"),
    family=f_gaussian(var.prior = 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
  )
  # increase burnin and n.iter below to improve MCMC convergence
  sim <- MCMCsim(sampler, store.all=TRUE, burnin=100, n.iter=200, 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")
}

## End(Not run)

mcmcsae documentation built on Sept. 12, 2025, 10:19 a.m.