R/Doptimality.R

In apear9/SSNdesign: An R Package for Optimal and Adaptive Sampling Designs on Stream Networks

#' Utility functions
#' 
#' @description
#'
#' Functions written in the pattern of the following \code{utility.function} can be used with \code{\link{optimiseSSNDesign}}. Off-the-shelf implementations are provided for the utility functions outlined in Falk et al. (2014) and modifications of these intended for creating adaptive designs.
#' 
#' @param ssn An object of class SpatialStreamNetwork
#' @param glmssn A fitted model object of class glmssn.
#' @param design.points A vector of pids corresponding to a set of observed sites in the obspoints slot of the SpatialStreamNetwork object.
#' @param prior.parameters A list of random functions that are parameterised in terms of n.draws.
#' @param n.draws A numeric scalar for the number of Monte Carlo draws to use when approximating the utility. 
#' @param extra.arguments A list of extra parameters that control the behaviour of the utility function. The distance matrices required to compute covariance matrices are also stored in this list. Note that these are generated inside \code{\link{optimiseSSNDesign}}.
#' @return A single number representing the expected utility for the design specified by \code{design.points}.
#' 
#' @details
#'
#' The following utility functions have been implemented.
#'
#' \itemize{
#' \item{DOptimality for optimal estimation of the fixed effects when the covariance parameters are known.}
#' \item{EDOptimality for optimal estimation of the fixed effects when the covariance parameters are not well known and must be estimated from the data.}
#' \item{CPOptimality for optimal estimation of the covariance parameters.}
#' \item{CPDOptimality for optimal estiamtion of both the covariance parameters and fixed effects parameters.}
#' \item{KOptimality for maximising the accuracy of kriging predictions when the covariance parameters are known.}
#' \item{EKOptimality for maximising the accuracy of kriging predictions when the covariance parameters are unknown and must be estimated fomr the data.}
#' \item{spaceFillingMaximin for constructing space filling designs that maximise the minimum distance between any pair of points in a design.}
#' \item{spaceFillingMorrisMitchell for constructing space filling designs using a modified version of the maximin utility function that also accounts for larger distances between pairs of points.}
#' \item{sequentialDOptimality for adaptively constructing designs that allow for optimal estimation of the fixed effects parameters when the covariance parameters are known.}
#' \item{sequentialEDOptimality for adaptively constructing designs that allow for optimal estimation of the fixed effects parameters when the covariance parameters are unknown and must be estimated from the data.}
#' \item{sequentialEDOptimality for adaptively constructing designs that allow for optimal estimation of the fixed effects parameters when the covariance parameters are unknown and must be estimated from the data.}
#' \item{sequentialCPOptimality for adaptively constructing designs that allow for optimal estimation of the covariance parameters.}
#' }
#'
#' When the returned value is \code{-1e9}, this indicates the utility function has failed to compute.
#' 
#' Do not worry about passing arguments to this function. All arguments are handled internally by \code{\link{optimiseSSNDesign}}. 
#' 
#' @examples 
#' 
#' \dontrun{
#' 
#' # Create stream network
#' s <- createSSN(100, systematicDesign(0.25), 
#' path = paste(tempdir(), "s.ssn", sep = "/"), importToR = TRUE)
#' createDistMat(s)
#' 
#' # Simulate data on network
#' s <- SimulateOnSSN(s, getSSNdata.frame(s), 
#' formula = ~1, coefficients = 1, CorParms = c(1,2,1,2,1,2,0.1),
#' addfunccol = "addfunccol")$ssn.object
#' 
#' # Fit a model to the simulated data
#' m <- glmssn(Sim_Values ~ 1, s, addfunccol = "addfunccol")
#' 
#' # Define the priors on the covariance parameters
#' p <- constructLogNormalCovPriors(m)
#' 
#' # Find the optimal design using D-optimality as the utility function
#' r <- optimiseSSNDesign(s, paste(tempdir(), "r.ssn", sep = "/"),
#'  m, 25, utility.function = DOptimality, prior.parameters = p)
#' 
#' # Plot result to check
#' plot(r$ssn.new, "Sim_Values")
#' 
#' }
#' 
#' @references 
#' 
#' Falk, M., McGree, J.M., and Pettit, A.N. (2014). Sampling designs on stream networks using the pseudo-Bayesian approach. \emph{Environmental and Ecological Statistics}, \emph{21}(4), 751-773.
#' 
#' @export
DOptimality <- function(ssn, glmssn, design.points, prior.parameters, n.draws, extra.arguments){
  
  ## Get data for design points. Note that design matrix is stored as obs.X in extra.arguments
  ind <- row.names(extra.arguments$obs.X) %in% design.points 
  X <- extra.arguments$obs.X[ind, ]
  Xt <- t(X)
  
  ## Get coordinates of points in case the covariance function includes a Euclidean dist-based component
  #  Note that the coordinates (with the same ordering as the design matrix) is stored as obs.C in extra.arguments
  ind.cds <- row.names(extra.arguments$obs.C) %in% design.points
  cds <- extra.arguments$obs.C[ind.cds, ] # names are already x and y
  
  ## Get distance matrices, etc.
  
  mat <- extra.arguments$Matrices.Obs
  ind.mat <- row.names(mat$d) %in% design.points
  mat$d <-  mat$d[ind.mat, ind.mat]
  mat$a <-  mat$a[ind.mat, ind.mat]
  mat$b <-  mat$b[ind.mat, ind.mat] 
  mat$w <-  mat$w[ind.mat, ind.mat]
  n.zero <- extra.arguments$net.zero.obs[ind.mat, ind.mat]
  
  ## Get other model parameters
  
  td <- glmssn$args$useTailDownWeight
  cm <- glmssn$args$CorModels
  un <- glmssn$args$use.nugget
  ua <- glmssn$args$use.anisotropy
  re <- glmssn$sampInfo$REs
  
  ## Perform MC simulations
  
  D <- vector("numeric", n.draws)
  
  for(i in 1:n.draws){
    
    theta.i <- prior.parameters[i, ]
    V <- SSN:::makeCovMat(
      theta.i, 
      mat$d, 
      mat$a, 
      mat$b, 
      mat$w, 
      n.zero, 
      cds[, "x"], 
      cds[, "y"], 
      cds[, "x"], 
      cds[, "y"], 
      td,
      cm, 
      un,
      ua,
      re
    )
    covbi <- Xt %*% solve(V) %*% X
    
    D[i] <- log(det(covbi)) 
    # Note that this is the same as -log(det(solve(covbi))) since the inverse of a determinant is the determinant of the inverse
    # Losing the inversion saves time, especially for large numbers of Monte Carlo draws.
    
  }
  
  if(any(is.infinite(D))){
    return(-1e9)
  } else {
    return(mean(D, na.rm = TRUE))
  }
  
}