R/uhcdenscalc.R

In aaarchmiller/uhcplots: Used-Habitat Calibration Plots (UHC Plots)

#' uhcdenscalc
#'
#' \code{uhcdenscalc} calculates kernal density estimates to create UHC plots.
#'
#' \code{uhcdenscalc} calculates density estimates for the environmental
#' characteristics associated with the observed and available locations in the
#' test data and also those associated with the randomly chosen locations
#' generated by the \code{uhcsim} or \code{uhcsimstrat} functions.
#'
#' @param rand_sims The characteristics associated with the simulated locations
#' in the test data set
#' @param dat The characteristics associated with observed locations in the
#' test data set
#' @param avail The characteristics associated with available locations in the
#' test data set
#' @param gridsize The size of grid for estimating the kernel density estimate
#' (kde). Default value is 500.
#'
#' @return A list of density estimates associated with the simulated locations,
#' (densrand), the observed locations (densdat), and the available locations
#' (densavail)
#'
#' @seealso Full archive of the data and code necessary to replicate the
#' manuscript at \url{http://doi.org/10.13020/D6T590}.
#'
#' @examples
#' # Simulate training data for the non-linear example
#' nonlinear.train <- uhcdatasimulator(nused = 100,
#'     navail = 10000,
#'     betas = c(2,-1),
#'     ntemp = 1000000,
#'     example = "non-linear")
#'
#' # Simulate test data for the non-linear example
#' nonlinear.test <- uhcdatasimulator(nused = 100,
#'     navail = 10000,
#'     betas = c(2,-1),
#'     ntemp = 1000000,
#'     example = "non-linear")
#'
#' # Fit GLM with quadratic relationship
#' train.correct <- glm(y~temp + I(temp^2),
#'    family = binomial,
#'    data = nonlinear.train)
#'
#' # Fit GLM with linear (misspecified) relationship
#' train.misspec <- glm(y~temp,
#'    family = binomial,
#'    data = nonlinear.train)
#'
#' # Simulate data for quadratic model
#' xhat.correct <- uhcsim(nsims = 1000,
#'    nused_test = 100,
#'    xmat = model.matrix(y~temp + I(temp^2), data = nonlinear.test)[,-1],
#'    fit_rsf = train.correct,
#'    z = as.matrix(nonlinear.test[,"temp"]))
#'
#' # Simulate data for linear (misspecified) model
#' xhat.misspec <- uhcsim(nsims = 1000,
#'    nused_test = 100,
#'    xmat = as.matrix(model.matrix(y~temp, data = nonlinear.test)[,2]),
#'    fit_rsf = train.misspec,
#'    z = as.matrix(nonlinear.test[,"temp"]))
#'
#' # Get density estimates for quadratic model
#' denshats.correct <- uhcdenscalc(rand_sims = xhat.correct[,,1],
#'    dat = subset(nonlinear.test, y==1, select="temp"),
#'    avail = subset(nonlinear.test, y==0, select="temp"))
#'
#' # Get density estimates for linear (misspecified) model
#' denshats.misspec <- uhcdenscalc(rand_sims = xhat.misspec[,,1],
#'    dat = subset(nonlinear.test, y==1, select="temp"),
#'    avail = subset(nonlinear.test, y==0, select="temp"))
#' @export
uhcdenscalc <- function(rand_sims, dat, avail, gridsize=500){

  # Density plot of x for each simulated data set & also for the test data set
  rand_sims <- as.matrix(rand_sims)
  dat <- as.matrix(dat)
  avail <- as.matrix(avail)
  alldat <- c(rand_sims,dat, avail) # combine data to get reasonable axis limits
  range.x <- c(min(alldat), max(alldat)) # to get reasonable axis limits
  nsims <- nrow(rand_sims)

  # to hold density estimates fU^(z) for all nsims datasets
  densrand <- matrix(NA, nsims,gridsize)
  # observed data density, fu(z)
  densdat <- KernSmooth::bkde(as.vector(dat), range.x=range.x, gridsize=gridsize)
  for(i in 1:nsims){ # densities for simulated data sets
    temp <- KernSmooth::bkde(rand_sims[i,], range.x=range.x, gridsize=gridsize)
    densrand[i,] <- temp$y #estimated density
  }
  # density of available points
  densavail <- KernSmooth::bkde(as.vector(avail),range.x=range.x, gridsize=gridsize)
  out <- list(densdat=densdat, densrand=densrand, densavail=densavail)
  return(out)
}