R/spde_smooth.R

In dill/gamUtils: Some extra utilities for generalized additive modelling

#' Matérn via SPDEs in GAMs
#'
#' Uses constructors from package R-INLA to build the stochastic partial differential equation Matérn model of Lindgren et al (2011). Details also in Miller et al (2019).
#'
#' For one dimensional smooths a mesh will be autogenerated using [`INLA::inla.mesh.1d`] with `k` evenly spaces basis functions (degree 2 B-splines). For other cases you need to supply `xt$mesh` which should be an object generated by [`INLA::inla.mesh`].
#'
#' If you want to compare "smoothing parameter" estimates with those produced by `R-INLA` then you need to add the `control = gam.control(scalePenalty = FALSE)` option to your [`mgcv::gam`] call.
#'
#' @note The `INLA` package is required to use this function see [the R-INLA website](http://r-inla.org) for details on installation.
#'
#' @inheritParams mgcv::smooth.construct.bs.smooth.spec
#'
#' @references
#' Miller, D. L., Glennie, R., & Seaton, A. E. (2019). Understanding the Stochastic Partial Differential Equation Approach to Smoothing. Journal of Agricultural, Biological and Environmental Statistics. https://doi.org/10.1007/s13253-019-00377-z
#' 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(4), 423–498. https://doi.org/10.1111/j.1467-9868.2011.00777.x
#' @importFrom mgcv smooth.construct
# @importFrom INLA inla.mesh.1d inla.spde.make.A inla.mesh.fem inla.mesh
#' @export
#' @author David L Miller
#' @examples
#' \dontrun{
#' library(mgcv)
#' set.seed(2) ## simulate some data... 
#' dat <- gamSim(1,n=400,dist="normal",scale=2)
#'
#' # fit model, ensuring smoothing parameters are comparable
#' # to those generated by INLA
#' b <- gam(y~s(x2, bs="spde", k=50),data=dat,
#'          control = gam.control(scalePenalty = FALSE))
#' # look at results
#' summary(b)
#' plot(b)
#'
#' # get hyperparameter estimates (as defined in Lindgren et al. 2011)
#' tau <- b$sp[1]
#' kappa <- b$sp[2]
#' # compute correlation range (rho) and marginal variance (sigma)
#' rho <- sqrt(8*1.5) / kappa
#' }
smooth.construct.spde.smooth.spec <- function(object, data, knots){

  # make sure that INLA is installed
  if (!requireNamespace("INLA", quietly = TRUE)) {
    stop("Package \"INLA\" is needed for this function to work. Please install it.",
      call. = FALSE)
  }

  # observation locations
  dim <- length(object$term)
  if (dim > 2 | dim < 1) stop("SPDE Matern can only be fit in 1D or 2D.")
  if (dim == 1) {
    x <- data[[object$term]]
  } else {
    x <- matrix(0, nrow = length(data[[1]]), ncol = 2)
    x[,1] <- data[[object$term[1]]]
    x[,2] <- data[[object$term[2]]]
  }
  # setup mesh or accept user mesh
  if (is.null(object$xt)) {
    if (dim == 1) {
      t <- seq(min(x), max(x), len=object$bs.dim)
      mesh <- INLA::inla.mesh.1d(loc=t, degree=2, boundary="free")
    } else {
      stop("For 2D, mesh must be supplied as argument xt$mesh in s(...,xt = )")
    }
  } else {
    if (class(object$xt$mesh) != "inla.mesh") stop("xt must be NULL or an inla.mesh object")
    mesh <- object$xt$mesh
  }
  # model matrix: projects parameters to observation locations on mesh 
  object$X <- as.matrix(INLA::inla.spde.make.A(mesh, x))
  # compute finite element matrices used as smoothing penalty matrices 
  inlamats <- INLA::inla.mesh.fem(mesh)
  object$S <- list()
  object$S[[1]] <- as.matrix(inlamats$c1)
  object$S[[2]] <- 2 * as.matrix(inlamats$g1)
  object$S[[3]] <- as.matrix(inlamats$g2)
  # L is a matrix with a column for each smoothing parameter (tau, kappa)
  # and a row for each smoothing matrix (c1, g1, g2).
  # The (i,j)^th entry of L contains the power that smoothing parameter i
  # is computed to before being multiplied by smoothing matrix j.
  # E.g. If (1, 2) has value 4, then smoothing parameter 2 (kappa) is taken
  # to the power 4 before being multiplied by smoothing matrix 1 (c1): i.e. kappa^4*c1
  # All of these computations for each element of L are then summed to create a single
  # smoothing matrix.
  object$L <- matrix(c(2,2,2,4,2,0), ncol = 2)
  # Rank is the basis dimension, it is repeated for each smoothing matrix
  object$rank <- rep(mesh$n, 3)
  object$knots <- mesh$n
  # As kappa > 0, the null space of the Matern SPDE is empty
  object$null.space.dim <- 0
  # Save the mesh
  object$mesh <- mesh
  object$df <- ncol(object$X)     # maximum DoF
  # set constraint matrix to have no cols for consistency with INLA
  object$C <- matrix(0, nrow=0, ncol=object$df)
  # Give object a class
  class(object) <- "spde.smooth"
  return(object)
}

# Prediction function for the `spde' smooth class
#' @rdname smooth.construct.spde.smooth.spec
# @importFrom INLA inla.spde.make.A
#' @export
#' @importFrom mgcv Predict.matrix
Predict.matrix.spde.smooth <- function(object, data){

  # make sure that INLA is installed
  if (!requireNamespace("INLA", quietly = TRUE)) {
    stop("Package \"INLA\" is needed for this function to work. Please install it.",
      call. = FALSE)
  }

  dim <- length(object$term)
  if (dim > 2 | dim < 1) stop("SPDE Matern can only be fit in 1D or 2D.")
  if (dim == 1) {
    x <- data[[object$term]]
  } else {
    x <- matrix(0, nrow = length(data[[1]]), ncol = 2)
    x[,1] <- data[[object$term[1]]]
    x[,2] <- data[[object$term[2]]]
  }
  Xp <- INLA::inla.spde.make.A(object$mesh, x)
  return(as.matrix(Xp))
}