Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!
Get Started for Free

IMRF

Interpolated Markov Random Field models


Description

IMRF is a syntax to specify random-effect terms of the forms considered by Lindgren et al. (2011) or Nychka et al. (2015, 2019). For example, using IMRF with its model argument provides good approximations of random effects with Matern correlation structure with fixed smoothness<2.

The random effects considered here all involve a multivariate Gaussian random effect over a lattice, from which the random-effect value in any spatial position is determined by interpolation of values on the lattice. IMRF stands for Interpolated Markov Random Field because the specific process considered on the lattice is currently known as a Gaussian Markov Random Field (see the Details for further information). Lindgren et al. considered irregular lattices that can be specified by the model argument, while Nychka et al. considered regular grids that can be specified by the other arguments.

The multIMRF syntax implements the multiresolution model of Nychka et al. Any multIMRF term in a formula is immediately converted to IMRF terms. This has distinct implications for controlling the parameters of these or other random effects in the model by init or fixed values: see Details if you need such control.

Usage

# IMRF( 1 | <coordinates>, model, nd, m, no, ce, ...) 
# multIMRF( 1 | <coordinates>, levels, margin, coarse=10L, 
#            norm=TRUE, centered=TRUE )

Arguments

model

An inla.spde2 object as produced by INLA::inla.spde2.matern or INLA::inla.spde2.pcmatern (see Examples below, and https://www.r-inla.org for further information).

levels

integer; Number of levels in the hierarchy, i.e. number of component IMRFs.

margin, m

integer; width of the margin, as a number of additional grid points on each side (applies to all levels of the hierarchy).

coarse

integer; number of grid points (excluding the margins) per dimension for the coarsest IMRF. The number of grids steps nearly doubles with each level of the hierarchy (see Details).

nd

integer; number of grid steps (excluding the margins) per dimension for the given IMRF.

norm, no

Boolean; whether to apply normalization (see Details), or not.

centered, ce

Boolean; whether to center the grid in all dimensions, or not.

...

Not documented, for programming purposes

Details

Gaussian Markov Random Field (MRF) and conditional autoregressive models are essentially the same thing, apart from details of specification. adjacency and AR1 random effects can be seen as specific MRFs. The common idea is the Markov-like property that the distribution of each element b_i of the random-effect b, given values of a few specific elements (the “neighbours” of i), is independent of other elements (i.e., of non-neighbours). The non-zero non-diagonal elements of a precision matrix characterize the neighbours.

Given the inferred vector b of values of the MRF on the lattice, the interpolation of the MRF in any focal point is of the form Ab where each row of A weights the values of b according to the position of the focal point relative to the vertices of the lattice. Following the original publications, for regular grids (NULL model), the weights are computed as <Wendland function>(<scaled Euclidean distances between focal point and vertices>); and for grids given by model=<inla.spde2 object>, the non-zero weights are the barycentric coordinates of the focal point in the enclosing triangle from the mesh triangulation (points from outside the mesh would have zero weights, so the predicted effect Ab=0).

The IMRF model defines both a lattice in space, the precision matrix for a Gaussian MRF over this lattice, and the A weights. The full specification of the MRF on irregular lattices is complex. The κ parameter considered by spaMM is the κ considered by Lindgren et al. The α argument of the INLA::inla.spde2.matern controls the smoothness of the approximated Matern model, as α=ν + d/2) where d is the dimension of the space. Correlation models created by INLA::inla.spde2.pcmatern are handled so as to give the same results as when INLA::inla.spde2.matern is used with the same mesh and alpha argument (thus, the extra functionalities of “pcmatern are ignored).

Not all options of the INLA functions may be compatible or meaningful when used with spaMM (only the effects of alpha and cutoff have been checked). The correlation models thus defined are fitted by the same methods as other models in spaMM. spaMM does not call INLA code except optionally INLA::inla.spde.make.A (if available to the R session) to construct the A matrix.

For the MRFs on default regular grids (missing model argument), the precision matrix is defined (up to a variance parameter) as M'M where the diagonal elements m_{ii} of M are 4+κ^2 and the m_{ij} for the four nearest neighbours are -1 (note that M'M involves more than these four neighbours).

The precision matrix defined in this way is the inverse of an heteroscedastic covariance matrix C, but by default a normalization is applied so that the random effect is homoscedastic. As for other random effects, the variance is further controlled by a multiplicative factor λ. The normalization is as follows: the design matrix of the random effect term is viewed as WAL where W is a diagonal normalization matrix, A is the above-described weight matrix, and L is a “square root” of C. If no normalization is applied, the covariance matrix of the random effect is of the form λALL'A', which is heteroscedastic; λ may then be quite different from the marginal variance of the random effect, and is difficult to describe in a simple way. Hence, by default, W is defined such that WALL'A'W' has unit diagonal; then, λ is the marginal variance of the random effect.

By default (meaning in particular that model is not used to specify a lattice defined by the INLA procedures), the IMRF lattice is rectangular (currently the only option) and is made of a core lattice, to which margins of margin steps are added on each side. The core lattice is defined as follows: in each of the two spatial dimensions, the range of axial coordinates is determined. The largest range is divided in nd-1 steps, determining nd values and step length L. The other range is divided in steps of the same length L. If it extends over (say) 2.5 L, a grid of 2 steps and 3 values is defined, and by default centered on the range (the extreme points therefore typically extend slightly beyond the grid, within the first of the additional steps defined by the margin; if not centered, the grid start from the lower coordinate of the range).

multIMRF implements multilevel IMRFs. It defines a sequence of IMRFs, with progressively finer lattices, a common κ value hy_kap for these IMRFs, and a single variance parameter hy_lam that determines λ values decreasing by a factor of 4 for successive IMRF terms. By default, each component IMRF is normalized independently as described above (as in Nychka et al. 2019), and hy_lam is the sum of the variances of these terms (e.g., if there are three levels and hy_lam=1, the successive variances are (1,1/4,1/16)/(21/16) ). The nd of the first IMRF is set to the coarse value, and its lattice is defined accordingly. If coarse=4 and margin=5, a grid of 14 coordinates is therefore defined over the largest range. In the second IMRF, the grid spacing is halved, so that new steps are defined halfway between the previous ones (yielding a grid of 27 step in the widest range). The third IMRF proceeds from the second in the same way, and so on.

To control initial or fixed values of multIMRF κ and variance parameters, which are hyper-parameter controlling several IMRF terms, the hyper syntax shown in the Examples should be used. hyper is a nested list whose possible elements are named "1", "2", ... referring to successive multIMRF terms in the input formula, not to successive random effect in the expanded formula with distinct IMRF terms (see Examples). But the different IMRF terms should be counted as distinct random effects when controlling other parameters (e.g., for fixing the variances of other random effects).

References

D. Nychka, S. Bandyopadhyay, D. Hammerling, F. Lindgren, S. Sain (2015) A multiresolution gaussian process model for the analysis of large spatial datasets. Journal of Computational and Graphical Statistics 24 (2), 579-599. doi: 10.1080/10618600.2014.914946

D. Nychka, D. Hammerling, Mitchel. Krock, A. Wiens (2018) Modeling and emulation of nonstationary Gaussian fields. Spat. Stat. 28: 21-38. doi: 10.1016/j.spasta.2018.08.006

Lindgren F., Rue H., Lindström J. (2011) An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73: 423-498. doi: 10.1111/j.1467-9868.2011.00777.x

Examples

if (spaMM.getOption("example_maxtime")>6) {

data("blackcap") ## toy examples; but IMRF may be useful only for much larger datasets
## and when using the 'cutoff' parameter of INLA::inla.mesh.2d()

########################## Irregular lattice specified by 'model':
#
data("small_spde") ## load object of class 'inla.spde2', created and saved by :
  # spd <- sp::SpatialPointsDataFrame(coords = blackcap[, c("longitude", "latitude")],
  #                            data = blackcap)
  # small_mesh <- INLA::inla.mesh.2d(loc = INLA::inla.mesh.map(sp::coordinates(spd)), 
  #                           max.n=100, # only for demonstration purposes
  #                           max.edge = c(3, 20)) 
  # small_spde <- INLA::inla.spde2.matern(small_mesh)
  # save(small_spde, file="small_spde.RData", version=2)
#  
fit_SPDE <- fitme(migStatus ~ means + IMRF(1|longitude+latitude, model=small_spde), 
                  data=blackcap)
                  
########################## Regular lattices:   
# 
#Using 'hyper' to control fixed hyper-parameters
#
(mrf <- fitme(migStatus ~ 1 + (1|pos) + 
                          multIMRF(1|longitude+latitude,margin=5,levels=2), 
              data=blackcap, fixed=list(phi=1,lambda=c("1"=0.5),
              hyper=list("1"=list(hy_kap=0.1,hy_lam=1)))) )
              
# Using 'hyper' to control initial hyper-parameters
#
(mrf <- fitme(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2),
                data=blackcap, method="ML", fixed =list(phi=1),
                init=list(hyper=list("1"=list(hy_kap=0.1,hy_lam=1)))) )
                
# *Independent* IMRF terms with default rectangular lattice (often giving dubious results)
#
(mrf <- fitme(migStatus ~ 1 + IMRF(1|longitude+latitude,margin=5, nd=4L)
                              + IMRF(1|longitude+latitude,margin=5, nd=7L),
          data=blackcap,  
          fixed=list(phi=1,lambda=c(1/4,1/16),
                       corrPars=list("1"=list(kappa=0.1),"2"=list(kappa=0.1)))))
                    
}

spaMM

Mixed-Effect Models, with or without Spatial Random Effects

v3.10.0
CeCILL-2
Authors
François Rousset [aut, cre, cph] (<https://orcid.org/0000-0003-4670-0371>), Jean-Baptiste Ferdy [aut, cph], Alexandre Courtiol [aut] (<https://orcid.org/0000-0003-0637-2959>), GSL authors [ctb] (src/gsl_bessel.*)
Initial release
2022-02-06

We don't support your browser anymore

Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.