R/simconf.inla.R

In davidbolin/excursions: Excursion Sets and Contour Credibility Regions for Random Fields

## simconf.inla.R
##
##   Copyright (C) 2014 David Bolin, Finn Lindgren
##
##   This program is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   This program is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
##   You should have received a copy of the GNU General Public License
##   along with this program.  If not, see <http://www.gnu.org/licenses/>.

#' Simultaneous confidence regions for latent Gaussian models
#'
#' `simconf.inla` is used for calculating simultaneous confidence regions
#' for latent Gaussian models estimated using `INLA`.
#'
#' @param result.inla Result object from INLA call.
#' @param stack The stack object used in the INLA call.
#' @param name The name of the component for which to do the calculation. This
#' argument should only be used if a stack object is not provided, use the tag
#' argument otherwise.
#' @param tag The tag of the component in the stack for which to do the calculation.
#' This argument should only be used if a stack object is provided, use the name
#' argument otherwise.
#' @param ind If only a part of a component should be used in the calculations, this
#' argument specifies the indices for that part.
#' @param alpha Error probability for the region.
#' @param method Method for handling the latent Gaussian structure:
#' \describe{
#' \item{'EB' }{Empirical Bayes (Gaussian approximation of posterior).}
#' \item{'NI' }{Numerical integration (Calculation based on the Gaussian mixture
#' approximation of the posterior, as calculated by INLA).}}
#' @param n.iter Number or iterations in the MC sampler that is used for approximating
#' probabilities. The default value is 10000.
#' @param verbose Set to TRUE for verbose mode (optional).
#' @param link Transform output to the scale of the data using the link function as defined in
#' the model estimated with INLA (default FALSE).
#' @param max.threads Decides the number of threads the program can use. Set to 0 for
#' using the maximum number of threads allowed by the system (default).
#' @param compressed If INLA is run in compressed mode and a part of the linear
#' predictor is to be used, then only add the relevant part. Otherwise the
#' entire linear predictor is added internally (default TRUE).
#' @param seed Random seed (optional).
#' @param inla.sample Set to TRUE if inla.posterior.sample should be used for the MC
#' integration.
#'
#' @return An object of class "excurobj" with elements
#' \item{a }{The lower bound.}
#' \item{b }{The upper bound.}
#' \item{a.marginal }{The lower bound for pointwise confidence bands.}
#' \item{b.marginal }{The upper bound for pointwise confidence bands.}
#' @export
#' @details See [simconf()] for details.
#'
#'
#' @note This function requires the `INLA` package, which is not a CRAN package.
#' See <https://www.r-inla.org/download-install> for easy installation instructions.
#' @author David Bolin \email{davidbolin@@gmail.com}
#' @references Bolin et al. (2015) *Statistical prediction of global sea level
#' from global temperature*, Statistica Sinica, vol 25, pp 351-367.
#'
#' Bolin, D. and Lindgren, F. (2018), *Calculating Probabilistic Excursion Sets and Related Quantities Using excursions*, Journal of Statistical Software, vol 86, no 1, pp 1-20.
#' @seealso [simconf()], [simconf.mc()], [simconf.mixture()]
#' @examples
#' \dontrun{
#' if (require.nowarnings("INLA")) {
#'   n <- 10
#'   x <- seq(0, 6, length.out = n)
#'   y <- sin(x) + rnorm(n)
#'   mu <- 1:n
#'   result <- inla(y ~ 1 + f(mu, model = "rw2"),
#'     data = list(y = y, mu = mu), verbose = FALSE,
#'     control.compute = list(
#'       config = TRUE,
#'       return.marginals.predictor = TRUE
#'     ),
#'     num.threads = "1:1"
#'   )
#'
#'   res <- simconf.inla(
#'     result, name = "mu", alpha = 0.05,
#'     max.threads = 1, num.threads = "1:1")
#'
#'   plot(result$summary.random$mu$mean, ylim = c(-2, 2))
#'   lines(res$a)
#'   lines(res$b)
#'   lines(res$a.marginal, col = "2")
#'   lines(res$b.marginal, col = "2")
#' }
#' }
#'
simconf.inla <- function(result.inla,
                         stack,
                         name = NULL,
                         tag = NULL,
                         ind = NULL,
                         alpha,
                         method = "NI",
                         n.iter = 10000,
                         verbose = FALSE,
                         link = FALSE,
                         max.threads = 0,
                         compressed = TRUE,
                         seed = NULL,
                         inla.sample = TRUE) {
  if (!requireNamespace("INLA", quietly = TRUE)) {
    stop("This function requires the INLA package (see www.r-inla.org/download-install)")
  }
  if (missing(result.inla)) {
    stop("Must supply INLA result object")
  }

  if (missing(alpha)) {
    stop("Must supply error probability alpha")
  }

  if (result.inla$.args$control.compute$config == FALSE) {
    stop("INLA result must be calculated using control.compute$config=TRUE")
  }

  n <- length(result.inla$misc$configs$config[[1]]$mean)

  if (!missing(ind)) {
    ind <- private.as.vector(ind)
  }


  # Get indices for the component of interest in the configs
  tmp <- inla.output.indices(result.inla,
                             name = name, stack = stack,
                             tag = tag, compressed = compressed
  )
  ind.stack <- tmp$index
  result.inla.orig <- result.inla
  if (tmp$result.updated) {
      result.inla <- tmp$result
      ind.stack.original <- tmp$index.original
  } else {
      ind.stack.original <- ind.stack
  }
  n <- length(result.inla$misc$configs$config[[1]]$mean)
  n.out <- length(ind.stack)
  ind.int <- seq_len(n.out)
  # ind is assumed to contain indices within the component of interest
  if (!missing(ind) && !is.null(ind)) {
      ind.int <- ind.int[ind]
      ind.stack <- ind.stack[ind]
      ind.stack.original <- ind.stack.original[ind]
  }
  ind <- ind.stack
  ind.original <- ind.stack.original
  
  links <- NULL
  if (link) {
    links <- result.inla$misc$linkfunctions$names[
      result.inla$misc$linkfunctions$link
    ]
    links <- links[ind]
  }

  n.theta <- result.inla$misc$configs$nconfig

  if (method == "EB") {
    for (i in 1:n.theta) {
      config <- private.get.config(result.inla, i)
      if (config$lp == 0) {
        break
      }
    }
    res <- simconf(
      alpha = alpha, mu = config$mu, Q = config$Q,
      vars = config$vars, n.iter = n.iter, ind = ind,
      verbose = verbose, max.threads = max.threads, seed = seed
    )
    res$meta$call <- match.call()
    return(private.simconf.link(res, links, link))
  } else if (method == "NI") {
    mu <- Q <- vars <- list()
    w <- rep(0, n.theta)
    for (i in 1:n.theta) {
      config <- private.get.config(result.inla, i)
      mu[[i]] <- config$mu
      Q[[i]] <- config$Q
      vars[[i]] <- config$vars
      w[i] <- config$lp
    }
    w <- exp(w) / sum(exp(w))
    if (inla.sample) {
      num.threads <- INLA::inla.getOption("num.threads")
      if (max.threads > 0) {
        if (is.character(num.threads)) {
          num.threads <- as.numeric(strsplit(num.threads, ":")[[1]])
        }
        num.threads <- min(max.threads, num.threads[1])
      }
      s <- suppressWarnings(INLA::inla.posterior.sample(n.iter, result.inla.orig,
                                                        use.improved.mean = FALSE,
                                                        skew.corr = FALSE,
                                                        num.threads = num.threads))
      samp <- matrix(0, n.iter, length(ind))

      for (i in seq_len(n.iter)) {
        samp[i, ] <- s[[i]]$latent[ind]
      }

      mu.m <- matrix(0, n.theta, length(ind))
      sd.m <- matrix(0, n.theta, length(ind))
      for (k in seq_len(n.theta)) {
        mu.m[k, ] <- mu[[k]][ind]
        sd.m[k, ] <- sqrt(vars[[k]][ind])
      }
      limits <- c(-1000, 1000)
      a.marg <- sapply(seq_len(length(ind)), function(i) {
        Fmix_inv(
          p = alpha / 2,
          mu = mu.m[, i], sd = sd.m[, i],
          w = w, br = limits
        )
      })

      b.marg <- sapply(seq_len(length(ind)), function(i) {
        Fmix_inv(
          p = 1 - alpha / 2,
          mu = mu.m[, i], sd = sd.m[, i],
          w = w, br = limits
        )
      })
      while (min(a.marg) == limits[1] || max(b.marg) == limits[2]) {
        limits <- 2 * limits
        a.marg <- sapply(seq_len(length(ind)), function(i) {
          Fmix_inv(
            p = alpha / 2,
            mu = mu.m[, i], sd = sd.m[, i],
            w = w, br = limits
          )
        })

        b.marg <- sapply(seq_len(length(ind)), function(i) {
          Fmix_inv(
            p = -alpha / 2,
            mu = mu.m[, i], sd = sd.m[, i],
            w = w, br = limits
          )
        })
      }

      r.o <- optimize(fmix.samp.opt,
        interval = c(0, alpha), mu = mu.m, alpha = alpha,
        sd = sd.m, w = w, limits = limits, samples = samp,
        verbose = verbose
      )

      a <- sapply(seq_len(length(ind)), function(i) {
        Fmix_inv(
          p = r.o$minimum / 2,
          mu = mu.m[, i], sd = sd.m[, i],
          w = w, br = limits
        )
      })

      b <- sapply(seq_len(length(ind)), function(i) {
        Fmix_inv(
          p = 1 - r.o$minimum / 2,
          mu = mu.m[, i], sd = sd.m[, i],
          w = w, br = limits
        )
      })

      res <- list(
        a = a, b = b, a.marginal = a.marg, b.marginal = b.marg,
        mean = mu, vars = vars
      )

      res$meta <- list(
        calculation = "simconf",
        alpha = alpha,
        n.iter = n.iter,
        ind = ind,
        call = match.call()
      )
      class(res) <- "excurobj"

      return(private.simconf.link(res, links, link))
    } else {
      res <- simconf.mixture(
        alpha = alpha, mu = mu, Q = Q, vars = vars,
        w = w, n.iter = n.iter, ind = ind, verbose = verbose,
        max.threads = max.threads, seed = seed
      )
      res$meta$call <- match.call()
      return(private.simconf.link(res, links, link))
    }
  } else {
    stop("Method must be EB or NI")
  }
}