R/sim2D_RandSet_HPPP.R

In jmleach-bst/sim2Dpredictr: Simulate Outcomes Using Spatially Dependent Design Matrices

#' Generate a Random Set Using a Poisson Process and Random Radii About Events
#'
#' A random set is generated by using a Poisson process in 2D space to choose
#' 'event' locations, about which a circle of random radius is 'drawn'. The 
#' union of the circles defines ultimately defines the set.
#'
#' @return A dataframe with columns for subject ID, x-coordinates, 
#' y-coordinates, and associated radii.
#'
#' @param xlim,ylim These are the 2D image limits. Defaults for both are 
#' \code{c(0, 1)}. It is not recommended to alter these arguments unless 
#' changing the limits has a specific practical utility.
#' @param N A scalar value determining the number of images to create.
#' @param radius.bounds A 2-element vector whose first and second entries 
#' determine the minimum and maximum radius sizes, respectively; these will 
#' be the bounds of the uniform distribution used to draw the radii. If 
#' \code{sub.area = TRUE}, then use \code{radius.bounds.min.sa} and 
#' \code{radius.bounds.max.sa}.
#' @param random.lambda \code{random.lambda = TRUE} allows the lambda 
#' (mean/intensity) parameter in the Poisson process to vary randomly by 
#' subject.
#' @param lambda A scalar value specifying the mean/intensity value of the 
#' Poisson process. If \code{random.lambda = FALSE} then this is the parameter
#' used to generate the binary image for each subject. If 
#' \code{random.lambda = TRUE}, then this is the mean parameter in the
#' distribution used to draw subject-specific lambda.
#' @param lambda.sd Only utilized when \code{random.lambda = TRUE}, and 
#' specifies the standard deviation in the distribution used to draw 
#' subject-specific lambda.
#' @param prior Only utilized when \code{random.lambda = TRUE}, and specifies
#' the distribution from which to draw the subject-specific lambda.
#' Options are \code{c("gaussian", "gamma")}.
#' @param lambda.bound Only utilized when \code{random.lambda = TRUE}, and 
#' allows the user to specify a lower and upper bound for the subject-specific
#' lambda; if the randomly selected value is outside of this range, then 
#' another draw is taken. This continues until a value is selected within the
#' specified bounds. If no bounds are desired then specify 
#' \code{lambda.bound = NULL}.
#' @param sub.area When \code{sub.area = TRUE}, a random sub-section of the 
#' image is chosen, within which the Poisson process is used to generate the
#' binary image.
#' @param min.sa,max.sa Only utilized when \code{sub.area = TRUE}, and
#' determines the width and height of the  minimum and maximum sub-areas; 
#' e.g., if \code{min.sa = c(0.1, 0.1)}, then the smallest possible random 
#' sub-area is a 0.1 x 0.1 square.
#' @param radius.bounds.min.sa,radius.bounds.max.sa Only utilized when 
#' \code{sub.area = TRUE}, and specifies \code{radius.bounds} for the minimum
#' and maximum sub-areas, respectively. This information is used to adaptively
#' alter the bounds in between the minimum and maximum sub-areas.
#' @param print.subj.sa,print.lambda,print.iter These arguments are either
#' \code{TRUE} or \code{FALSE}, and define print options for checking that the
#' function is working as the user intends. \code{print.subj.sa = TRUE} prints
#' the x-and y-limits for each subject's sub-area. \code{print.lambda = TRUE} 
#' prints each subject's mean and realized events; the means will be the same
#' unless \code{random.lambda = TRUE}, but the number of realized events will
#' always vary. \code{print.iter = TRUE} is only used when 
#' \code{random.lambda = TRUE} and \code{is.null(lambda.bound) = FALSE}, and
#' shows iterations for re-drawing when the randomly selected intensity is 
#' outside the specified bounds.
#' @importFrom stats rnorm rgamma rpois runif
#' @importFrom Rdpack reprompt
#' @references
#' \insertRef{Cressie+Wikle:2011}{sim2Dpredictr}
#' @export
sim2D_RandSet_HPPP <- function(N, xlim = c(0, 1), ylim = c(0, 1),
                               radius.bounds = c(0.02, 0.1),
                               lambda = 50, lambda.sd = 10, 
                               lambda.bound = NULL,
                               prior = "gamma", 
                               random.lambda = FALSE, 
                               sub.area = FALSE,
                               min.sa = c(0.1, 0.1), max.sa = c(0.3, 0.3),
                               radius.bounds.min.sa = c(0.02, 0.05),
                               radius.bounds.max.sa = c(0.08, 0.15),
                               print.subj.sa = FALSE, print.lambda = FALSE,
                               print.iter = FALSE){
  # store
  subj.event <- c()
  
  # store x- and y-coordinates for 'event' locations
  xcoord <- c()
  ycoord <- c()
  
  if (sub.area == TRUE) {
    
    if ( (max.sa[1] < min.sa[1]) | (max.sa[2] < min.sa[2]) ){
      stop("Must have max.sa[1] >= min.sa[1] and max.sa[2] >= min.sa[2]. \n ")
    }
    
    if (any(c(max.sa, min.sa) < 0)) {
      stop("Cannot have sub-area dimensions < 0. \n")
    }
    
    if ( max.sa[1] > abs(xlim[2] - xlim[1]) | max.sa[2] > abs(ylim[2] - ylim[1]) ) {
      stop("Sub-area bounds may not exceed image limits. \n")
    }
    
    # store sub-area limits
    xlim.sa <- c()
    ylim.sa <- c()
    
    ######################################################
    # The following code allows flexible 'radii.range'   #
    # based on random sub-area; i.e. the range of radius #
    # sizes will be proportional to the size of the      #
    # random sub-area.                                   #
    ######################################################
    
    # minimum possible sub-area
    A.min <- prod(min.sa)
    # maximum possible sub-area size
    A.max <- prod(max.sa)
    
    A <- matrix(c(1, 1, A.min, A.max), nrow = 2, ncol = 2)
    B.min <- solve(t(A) %*% A) %*% t(A) %*% c(radius.bounds.min.sa[1], 
                                              radius.bounds.min.sa[2])
    B.max <- solve(t(A) %*% A) %*% t(A) %*% c(radius.bounds.max.sa[1], 
                                              radius.bounds.max.sa[2])
    
    est.lm <- function(A,B){
      B[1] + A*B[2]
    }
  }
  
  # store the radii about each 'event'
  radii <- c()
  
  # store subject ID
  subj <- c()
  
  if (random.lambda == TRUE) {
    lambda.var <- lambda.sd^2
    
    # Convert mean, variance into alpha, beta parameters for gamma dist.
    gamma.alpha.beta <- function(mu, sigma2){
      alpha <-  (mu ^ 2) / sigma2
      beta <- sigma2 / mu
      param.gamma <- data.frame(alpha, beta)
      return(param.gamma)
    }
    gab <- gamma.alpha.beta(lambda,lambda.var)
    alpha <- gab$alpha
    beta <- gab$beta
  }
  
  # Initialization for subject ID; necessary for correctly tying coordinates
  # and radii to subjects. This will be updated after each iteration.
  old.subj <- 1
  
  # Begin generation of random binary images.
  for (i in 1:N){
    
    # Draw lambda from a Gaussian or Gamma Distribution.
    if (random.lambda == TRUE) {
      if (prior == "gaussian") {
        lambda.subj <- rnorm(1, mean = lambda,sd = lambda.sd)
        
        # If the random subject-specific lambda is outside the
        # user-specified bounds then re-draw.
        if (is.null(lambda.bound) == FALSE) {
          if (print.iter == TRUE & (lambda.subj > lambda.bound[2] |
                                    lambda.subj < lambda.bound[1])) {
            cat("Subject ", i, "lambda = ", lambda.subj, "REJECT!", "\n")
          }
          while (lambda.subj > lambda.bound[2] |
                 lambda.subj < lambda.bound[1]) {
            lambda.subj <- rnorm(1, mean = lambda, sd = lambda.sd)
            if (print.iter == TRUE & (lambda.subj > lambda.bound[2] |
                                      lambda.subj < lambda.bound[1])) {
              cat("Subject ", i, "lambda = ", lambda.subj, "REJECT!", "\n")
            }
          }
        }
      }
      if (prior=="gamma"){
        lambda.subj <- rgamma(1,alpha,scale=beta)
        
        # If the random subject-specific lambda is outside the
        # user-specified bounds then re-draw.
        if (is.null(lambda.bound) == FALSE) {
          if (print.iter == TRUE & (lambda.subj > lambda.bound[2] |
                                    lambda.subj < lambda.bound[1])) {
            cat("Subject", i, "lambda = ", lambda.subj, "REJECT!", "\n")
          }
          while (lambda.subj > lambda.bound[2] |
                 lambda.subj < lambda.bound[1]){
            lambda.subj <- rgamma(1,alpha,scale=beta)
            if (print.iter == TRUE & (lambda.subj > lambda.bound[2] |
                                      lambda.subj < lambda.bound[1])){
              cat("Subject", i, "lambda=", lambda.subj, "REJECT!", "\n")
            }
          }
        }
      }
    } else {
      # i.e. common lambda for all subjects.
      lambda.subj <- lambda
    }
    
    # Use the random lambda to draw the number of events from a 
    # Poisson distribution.
    subj.event[i] <- rpois(1, lambda.subj)
    if (print.lambda == TRUE){
      cat("Subject ",i, " draws with mean ",lambda.subj, "and obtains",
          subj.event[i], "events", "\n")
    }
    
    # For HPPP the number of events can be drawn
    # from a standard Poisson distribution (as above) and the locations of
    # the events will be uniformly distributed about the 2D area; i.e. we
    # can draw the coordinates from a 2D uniform distribution.
    
    # Use the entire area in assigning locations for events;
    # i.e. NO random sub-area.
    if (sub.area == FALSE) {
      subj[old.subj:(old.subj - 1 + subj.event[i])] <- i
      xcoord[old.subj:(old.subj - 1 + subj.event[i])] <- runif(subj.event[i], 
                                                               xlim[1], 
                                                               xlim[2])
      ycoord[old.subj:(old.subj - 1 + subj.event[i])] <- runif(subj.event[i], 
                                                               ylim[1],
                                                               ylim[2])
    }
    
    # Use only a random sub-area in assigning locations for events.
    if (sub.area == TRUE) {
      
      # Determine the random sub-area.
      bx <- xlim[2] - min.sa[1]
      by <- ylim[2] - min.sa[2]
      xlim.sa[1] <- runif(1, xlim[1], bx)
      ylim.sa[1] <- runif(1, ylim[1], by)
      xlim.sa[2] <- runif(1, xlim.sa[1] + min.sa[1],
                          min(xlim[2], xlim.sa[1] + max.sa[1]))
      ylim.sa[2] <- runif(1, ylim.sa[1] + min.sa[2],
                          min(ylim[2], ylim.sa[1] + max.sa[2]))
      
      # Assign locations to events
      subj[old.subj:(old.subj - 1 + subj.event[i])] <- i
      xcoord[old.subj:(old.subj - 1 + subj.event[i])] <- runif(subj.event[i],
                                                               xlim.sa[1],
                                                               xlim.sa[2])
      ycoord[old.subj:(old.subj - 1 + subj.event[i])] <- runif(subj.event[i],
                                                               ylim.sa[1],
                                                               ylim.sa[2])
      if (print.subj.sa == TRUE){
        cat("Subject", i, " has xlim.sa=", xlim.sa, 
            "and ylim.sa=", ylim.sa, "\n")
      }
    }
    
    #################################################################
    # Step 2: Observe the simplicity of generating a random radius. #
    #         It is the same no matter what variation of PPP used.  #
    #################################################################
    
    # The 'area' of the random sub-area.
    if (sub.area == TRUE) {
      obs.A <- (xlim.sa[2] - xlim.sa[1]) * (ylim.sa[2] - ylim.sa[1])
      radius.bounds <- c(est.lm(obs.A, B.min), est.lm(obs.A, B.max))
    }
    radii[old.subj:(old.subj - 1 + subj.event[i])] <- runif(subj.event[i],
                                                            radius.bounds[1],
                                                            radius.bounds[2])
    # ensure the next subject's data are stored in the correct location.
    old.subj <- old.subj + subj.event[i]
  }
  
  data <- data.frame(subj = subj, 
                     xcoord = xcoord, 
                     ycoord = ycoord, 
                     radii = radii)
  return(data)
}