MakeIntegrationStack.R

#' Function to create stack for integration points
#'
#' @param mesh INLA mesh (e.g., \code{$mesh} of object generated by \code{\link{MakeSpatialRegion}}).
#' @param data \code{data.frame} with at least columns "X" and "Y" (or \code{coordnames}) giving locations, plus optional other covariate. Alternatively, a \code{SpatialPointsDataFrame} object.
#' @param area Area around each integration point (e.g., the \code{$w} element of the object generated by \code{\link{MakeSpatialRegion}}).
#' @param tag name for the stack.
#' @param coordnames Names for coordinates in data (if \code{data} is not a \code{SpatialPointsDataFrame}).
#' @param InclCoords should coordinates be included in data? (defaults to \code{FALSE}).
#' @return An INLA stack for integration
#'
#' @export
#' @import INLA
MakeIntegrationStack <- function(mesh, data, area, tag='mesh',
                                 coordnames=c("X","Y"), InclCoords=FALSE) {
  if(class(data)!="SpatialPointsDataFrame" & !all(coordnames%in%names(data))) stop("Coordinates not in the data")

  if(class(data)=="SpatialPointsDataFrame") {
    coordnames <- colnames(data@coords)
    Names <- colnames(data@data)
  } else {
    Names  <- names(data)[!names(data)%in%coordnames]
  }
  if(InclCoords) Names  <- c(coordnames, Names)

  Points <- cbind(c(mesh$loc[,1]), c(mesh$loc[,2]))
  colnames(Points) <- coordnames
  NearestCovs <- GetNearestCovariate(points=Points, covs=data)
  NearestCovs$Intercept <- rep(1,nrow(NearestCovs))
  if(InclCoords) {
    NearestCovs@data[,colnames(NearestCovs@coords)] <- NearestCovs@coords
  }

  # Projector matrix for integration points.
  projmat.ip <- Matrix::Diagonal(mesh$n, rep(1, mesh$n))  # from mesh to integration points

  stk.ip <- inla.stack(data=list(resp=cbind(rep(0,mesh$n), NA), e=area),
                       A=list(1,projmat.ip), tag=tag,
                       effects=list(NearestCovs@data, list(i=1:mesh$n)))
  return(stk.ip)
}