inst/RandomFieldsSupport/simDynoccSpatial.R

In AHMbook: Functions and Data for the Book 'Applied Hierarchical Modeling in Ecology' Vols 1 and 2

# 1. Define an R function to generate dynamic presence/absence systems with 'space'

# Code to define a function for simulating data.
# AHM2 - 9.6.1.1

#### Uses 'fields' instead of 'RandomFields' ####

simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
      mean.psi1 = 0.4, beta.Xpsi1 = 0,
      range.phi = c(0.8, 0.8), beta.Xphi = 0,
      range.gamma = c(0.1, 0.1), beta.Xgamma = 0,
      range.p = c(0.4, 0.4), beta.Xp = 0,
      theta.XAC = 5000, beta.XAC = c(0, 0, 0, 0), beta.Xautolog = c(0, 0),
      trend.sd.site = c(0, 0), trend.sd.survey = c(0, 0),
      seed.XAC = NA, seed = NULL, show.plots= TRUE, ask.plot = TRUE, verbose=TRUE) {
  #
  #      Written by Marc Kéry, 2014-2018
  #
  # Function to simulate detection/nondetection data in a square area
  #   under a very general dynamic site-occ model, including the
  #   following effects:
  # (1) annual variation in the probabilities of patch persistence,
  #   colonization and detection can be specified by the bounds of a
  #   uniform distribution.
  # (2) one site-, site/year-, and site/year/rep-specific covariate
  #   is allowed to affect the probabilities of occupancy
  #   (beta.Xpsi1 for site-covariate), colonisation/persistence
  #   (beta.Xgamma, beta.Xphi, for yearly site-covariate), and
  #   detection (beta.Xp for observational covariate), respectively.
  # (3) a single, spatially structured covariate for habitat suitability
  #   may affect all parameters via coefficient beta.XAC (for a
  #   biologically reasonable way, choose coefficients with the same sign
  #   for all 4 (mediated by underlying density).
  #   That spatial covariate is simulated as a Gaussian random field
  #   with negative exponential correlation function with
  #   'range parameter' theta.XAC
  # (4) autologistic effects (beta.Xautolog) in persistence and colonization
  #   probability can be chosen, which fits a logistic regression of
  #   these parameters on the proportion of occupied neighbouring cells
  #   (in a queen's or 2nd order neighbourhood) during the previous time step
  # (5) Additional detection heterogeneity can be introduced
  #   at the site- or the individual survey level, with the possibility of a
  #   temporal trend in this heterogeneity. For instance, an annual trend in
  #   detection heterogeneity at the site or the survey level is specified by
  #   the value in the first and the last year.
  #   Hence, trend.sd.site = c(0, 1) will result in a linear trend in
  #   the magnitude of site heterogeneity in detection from 0 in the
  #   first year to 1 in the last year.
  #
  #
  # Function arguments:
  # -------------------
  #
  # *** Design of study and basic 'magnitude' of parameters ***
  # side – side length of square simulation area. Therefore,
  #     the number of sites, or cells, M = side^2
  # nsurveys – Number of replicate surveys within a 'season', year or primary period
  # nyears – Number of years (or 'seasons')
  # mean.psi1 – intercept of occupancy probability in year 1
  # range.phi and range.gamma – bounds of uniform distribution from which
  #   annual intercepts for persistence (phi) and colonisation (gamma)
  #   are drawn
  # range.p – same for detection probability p
  #
  #
  # *** Covariates ***
  # beta.Xpsi1: coefficient of a site covariate in psi1
  # beta.Xphi: coefficient of a site/year covariate in phi
  # beta.Xgamma: coefficient of a site/year covariate in gamma
  # beta.Xp: coefficient of a site/year/rep covariate in p
  #
  #
  # *** Parameters governing the spatial correlations ***
  # theta.XAC: 'range parameter' of a covariate with exponential
  #   spatial correlation (i.e., a Gaussian random field is used as an
  #   environmental covariate). NOTE: if you want to set to zero the effects
  #   of this spatially autocorrelated variable, you CANNOT
  #   set theta.XAC=0 because this breaks the function,
  #   nor can you simply choose a very small value.
  #   Instead you MUST set the elements of coefficients vector beta.XAC
  #   to zero.
  # beta.XAC: vector of coefficients of that field for the 4 model params:
  #   psi1, phi, gamma, and p (in that order)
  # beta.Xautolog – vector of coefficients of autologistic covariate
  #   in the following order: persistence (phi), colonization (gamma).
  #   Autocovariate is computed at every season as the proportion of
  #   occupied cells in a queen's neighbourhood around each cell.
  #
  #
  # *** Detection heterogeneity ***
  # trend.sd.site: range of year-specific values of SD of Gaussian
  #   random site effects in p: c(1,1) specifies constant value of 1
  #   for all years, while c(0,1) specifies linear increase over the years
  #   from 0 to 1.
  # trend.sd.survey: range of year-specific values of standard deviation
  #   of Gaussian random survey effects in p: specification as
  #   for trend.sd.site
  #
  # *** Graphics control and other ***
  # seed – allows to 'fix' the simulation such that it becomes reproducible
  # ask.plot – if TRUE permits to browse through plots (otherwise if FALSE)

  if(FALSE) {x <- NULL; rm(x)}  # Stops R CMD check choking on 'curve'.

  # Checks and fixes for input data -----------------------------
  side <- round(side[1])
  nyears <- round(nyears[1])
  stopifnotGreaterthan(nyears, 1)
  nsurveys <- round(nsurveys[1])
  stopifnotProbability(mean.psi1)
  stopifnotProbability(range.phi) # bounds
  stopifnotProbability(range.gamma) # bounds
  stopifnotProbability(range.p) # bounds
  stopifNegative(theta.XAC, allowZero=FALSE)
  stopifnotLength(beta.XAC, 4)
  stopifnotLength(beta.Xautolog, 2)
  stopifnotLength(trend.sd.site, 2) # trend
  stopifNegative(trend.sd.site)
  stopifnotLength(trend.sd.survey, 2) # trend
  stopifNegative(trend.sd.survey)
  # ----------------------------------------------------------------
  # Restore graphical settings on exit -----------------------------
  if(show.plots) {
    oldpar <- par("mfrow", "mar", "cex.main", "cex.lab", "cex.axis")
    oldAsk <- devAskNewPage(ask = ask.plot && dev.interactive(orNone=TRUE))
    on.exit({par(oldpar); devAskNewPage(oldAsk)})
  }
  # ----------------------------------------------------------------

  # Create grid
  xcoord <- 1:side
  ycoord <- 1:side
  grid <- as.matrix(expand.grid(x=xcoord, y=ycoord))
  M <- side^2             # Total number of cells or sites

  # Compute adjacency matrix for grid
  neigh <- spdep::dnearneigh(as.matrix(grid), d1 = 0, d2 = sqrt(2) * 1 + 0.01)
  winnb <- spdep::nb2WB(neigh)        # Function to get CAR ingredients for BUGS
  nneigh <- winnb$num          # number of neighbours
  amatrix <- spdep::nb2mat(neigh)
  amatrix[amatrix > 0] <- 1    # Neighbours get a 1, non-neighbours a 0

  # Set up arrays needed
  site <- 1:M                                      # Sites
  year <- 1:nyears                                 # Years
  prob <- array(dim = c(side, side))               # p matrix
  psi <- muZ <- z <- array(dim = c(side, side, nyears))      # Occupancy, occurrence
  phi <- gamma <- array(NA, dim = c(side, side, (nyears-1))) # Survival, colonisation
  Xauto <- array(NA, dim = c(side, side, nyears))            # Autocovariate
  y <- p <- array(NA, dim = c(side, side, nsurveys, nyears)) # Det. histories and p

  # Create values of 1 spatially autocorrelated covariate XAC
  # Generate correlated random variables in a square
  
  if(requireNamespace("RandomFields", quietly=TRUE)) {
    RandomFields::RFoptions(seed=seed.XAC)     # Default NA; 88 gives cool pattern
    XAC <- matrix(RandomFields::RFsimulate(RandomFields::RMexp(var = 1, scale = theta.XAC),
        x=xcoord, y=ycoord, grid=TRUE)@data$variable1,
        ncol = side, byrow = TRUE)  # variance 1
    if(!is.na(seed.XAC))
      RandomFields::RFoptions(seed=NA)
  } else {
    message("Using package 'fields' instead of 'RandomFields'; see help(simDynoccSpatial).")
    if(!is.na(seed.XAC))
      set.seed(seed.XAC)
    obj <- circulantEmbeddingSetup(grid=list(x=xcoord, y=ycoord), Covariance="Exponential", aRange=theta.XAC)
    tmp <- try(circulantEmbedding(obj), silent=TRUE)
    if(inherits(tmp, "try-error"))
      stop("Simulation of random field failed.\nTry with smaller values for 'side' or 'theta.XAC'.")
    XAC <- matrix(tmp, ncol = side, byrow = TRUE)
  }
 
  set.seed(seed=seed)  # Default NULL; do this AFTER RFsimulate

  # Create four spatially unstructured covariates
  # Site covariate for psi1
  Xpsi1 <- matrix(runif(M, -2, 2), ncol = side)

  # Yearly-site covariates for phi and gamma
  Xphi <- Xgamma <- array(runif(M*nyears, -2, 2), dim = c(side, side, nyears))

  # Observational covariate for p
  Xp <- array(runif(M*nsurveys*nyears,-2,2), dim=c(side, side,nsurveys,nyears))

  # Draw values of baseline levels of the main parameters
  # (i.e., draw year effects if any)
  mean.phi <- runif(n = nyears-1, min = min(range.phi), max = max(range.phi))
  mean.gamma <- runif(n = nyears-1, min = min(range.gamma), max = max(range.gamma))
  mean.p <- runif(n = nyears, min = min(range.p), max = max(range.p))



  # (a) Simulate state process parameters: initial state (first year)
  psi[,,1] <- plogis(qlogis(mean.psi1) + beta.Xpsi1 * Xpsi1 +
     beta.XAC[1] * XAC) # psi1

  # (b) Simulate state in first year
  z[,,1] <- rbinom(M, 1, psi[,,1])  # Initial occurrence state

  # Compute value of autocovariate after first year = proportion of neighbours occupied
  # first vectorize and then put into matrix again
  Xautovec <- amatrix %*% c(z[,,1])
  Xauto[,,1] <- matrix(Xautovec/nneigh, ncol = side) # Put back in matrix by column again
  
  # Do the pre-loop plots
  # ---------------------
  if(show.plots) {
    tryPlot <- try( {
      # Plot effects of autocovariate on (year-specific) phi and gamma
      par(mfrow = c(1, 2))
      curve(plogis(qlogis(mean.phi[1]) + beta.Xautolog[1] * x), 0, 1, 
          main = "Persistence: \nphi ~ Year + Autocovariate", xlab = "Autocov. (prop. occupied neighb.)",
          ylab = "phi", ylim = c(0,1), frame = FALSE)
      for(k in 2:(nyears-1)){
         curve(plogis(qlogis(mean.phi[k])+beta.Xautolog[1]*x),0,1,add=TRUE)
      }

      curve(plogis(qlogis(mean.gamma[1]) + beta.Xautolog[2] * x), 0, 1, 
          main = "Colonization: \ngamma ~ Year + Autocovariate",
          xlab = "Autocovariate (prop. occupied neighb.)", ylab = "gamma", ylim = c(0,1), frame = FALSE)
      for(k in 2:(nyears-1)){
         curve(plogis(qlogis(mean.gamma[k])+beta.Xautolog[2]*x),0,1,add=TRUE)
      }

      # Simulate true system dynamics
      par(mfrow = c(2,2), mar = c(5,4,5,2), cex.main = 1.3, cex.lab = 1.5, cex.axis = 1.2)

      # Plot random field covariate XAC
      # rows are in x, columns in y direction
      image(1:side, 1:side, XAC, col=topo.colors(100), 
          main = paste("Gaussian random field XAC with \n neg. exponential correlation (range =", theta.XAC, ")"), 
          xlab = 'x', ylab = 'y')

      image(1:side, 1:side, psi[,,1], col=topo.colors(100),
          main = paste("Initial occupancy probability"), xlab = 'x', ylab = 'y')

      image(1:side, 1:side, z[,,1], col=c("white", "black"), 
          main = paste("Initial presence/absence (true system state z):\n black = occupied, white = unoccupied"),
          xlab = 'x', ylab = 'y')
      abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")

      image(1:side, 1:side, Xauto[,,1], col=topo.colors(100), 
          main = "Autocovariate between year 1 and year 2", xlab = 'x', ylab = 'y')
    }, silent = TRUE)
    if(inherits(tryPlot, "try-error")) {
      show.plots <- FALSE
      tryPlotError(tryPlot)
    }
  }

  # (c) Simulate state process parameters: time steps 2:nyears
  for(k in 2:nyears){
    par(mfrow = c(2,2), mar = c(5,4,5,2), cex.main = 1.3, cex.lab = 1.5, cex.axis = 1.2)
    if(verbose)
      cat(paste("** Year", k, "**\n"))
    # Compute colonisation and extinction parameters and plot
    phi[,,k-1] <- plogis(qlogis(mean.phi[k-1]) + beta.Xphi * Xphi[,,k-1] + beta.XAC[2] * XAC +
      beta.Xautolog[1] * Xauto[,,k-1])
    gamma[,,k-1] <- plogis(qlogis(mean.gamma[k-1]) + beta.Xgamma * Xgamma[,,k-1] + beta.XAC[3] * XAC +
      beta.Xautolog[2] * Xauto[,,k-1])

    # Compute latent states and plot
    muZ[,,k] <- z[,,k-1]*phi[,,k-1] + (1-z[,,k-1])*gamma[,,k-1]
    z[,,k] <- rbinom(M, 1, muZ[,,k])

    # Compute autocovariate and plot
    Xautovec <- amatrix %*% c(z[,,k])
    Xauto[,,k] <- matrix(Xautovec/nneigh, ncol = side) # re-assemble by column
    
    # Do the in-loop plots
    # --------------------
    if(show.plots) {
      tryPlot <- try( {
        image(1:side, 1:side, phi[,,k-1], col=topo.colors(100),
              main = paste("Persistence between year", k-1, "and year", k), xlab = 'x', ylab = 'y')
        image(1:side, 1:side, gamma[,,k-1], col=topo.colors(100),
              main = paste("Colonization between year", k-1, "and year", k), xlab = 'x', ylab = 'y')
        image(1:side, 1:side, z[,,k], col=c("white", "black"),
            main = paste('Presence/absence (z) in year', k, ':\n black = occupied, white = unoccupied'),
            xlab = 'x', ylab = 'y')
        abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")
        image(1:side, 1:side, Xauto[,,k], col=topo.colors(100),
            main = paste("Autocovariate between year", k, "and year", k+1), xlab = 'x', ylab = 'y')
      }, silent = TRUE)
      if(inherits(tryPlot, "try-error")) {
        show.plots <- FALSE
        tryPlotError(tryPlot)
      }
    }

  }

  # (d) Observation process parameters
  # First choose values of annual SD of p random effects
  sd.site <- seq(from = trend.sd.site[1], to = trend.sd.site[2],
     length.out = nyears)
  sd.survey <- seq(from = trend.sd.survey[1], to = trend.sd.survey[2],
     length.out = nyears)
  for(k in 1:nyears){
    # Site random effects
    eps1 <- matrix(rnorm(n = M, mean = 0, sd = sd.site[k]), ncol = side)
    # Survey random eff.
    eps2 <- rnorm(n = nsurveys, mean = 0, sd = sd.survey[k])
    for(j in 1:nsurveys){
      p[,,j,k] <- plogis(qlogis(mean.p[k]) + beta.Xp * Xp[,,j,k] + beta.XAC[4] * XAC + eps1[,] + eps2[j])
    }
  }

  # Simulate actual observation process (also updating entire grid in one go)
  for(k in 1:nyears){                 # Loop over years
    for(j in 1:nsurveys){               # Loop over replicates
      prob <- z[,,k] * p[,,j,k]  # zero out p for unoccupied sites
      y[,,j,k] <- rbinom(M, 1, prob)
    # image(1:side, 1:side, y[,,j,k], main = paste("Year", k, "and rep", j))
        # Look at clumped pattern in y
    }
  }

  # Derived quantities
  # Compute annual population occupancy
  for (k in 2:nyears){
     psi[,,k] <- psi[,,k-1]*phi[,,k-1] + (1-psi[,,k-1])*gamma[,,k-1]
  }
  mean.psi <- apply(psi, 3, mean)              # Average psi over sites

  # Compute true and observed number of occupied sites
  zobs <- apply(y, c(1,2,4), max)
  nocc <- apply(z, 3, sum)
  nocc.obs <- apply(zobs, 3, sum)

  # Do the post-loop plots
  # ----------------------
  psi.app <- apply(apply(y, c(1,2,4), max), 3, mean)
  if(show.plots) {
    tryPlot <- try( {
      # (4) More plots comparing true states and observations
      # Plot realised and apparent occupancy
      par(mfrow = c(1,1))
      plot(year, apply(z, 3, mean), type = "l", xlab = "Year", ylab = "Occupancy or Detection prob.", col = "red", xlim = c(0,nyears+1), ylim = c(0,1), lwd = 2, lty = 1, frame.plot = FALSE, las = 1)
      lines(year, mean.p, type = "l", col = "red", lwd = 2, lty = 2)
      lines(year, psi.app, type = "l", col = "black", lwd = 2)
      text(0.8*nyears, 0.1, labels = "red solid - true occupancy prob.\n red dashed - detection prob.\n black - observed proportion occupied", cex = 1)

      # Plots comparing true and observed latent states
      par(mfrow = c(2,2), mar = c(5,4,5,2), cex.main = 1.3, cex.lab = 1.5, cex.axis = 1.2)
      for(k in 1:nyears){
        image(1:side, 1:side, z[,,k], col=c("white", "black"), main = paste('Presence/absence (z) in year', k), xlab = 'x', ylab = 'y')
        abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")
        image(1:side, 1:side, zobs[,,k], col=c("white", "black"), main = paste('Ever_detected (zobs) in year', k), xlab = 'x', ylab = 'y')
        abline(h = 0:side+0.5, v = 0:side+0.5, col = "lightgrey")
      }
    }, silent = TRUE)
    if(inherits(tryPlot, "try-error"))
      tryPlotError(tryPlot)
  }
  
  # Compute values of naive autocovariate (observed prop.
  #     of occupied neighbouring cells)
  Xautoobs <- array(NA, dim = dim(zobs))
  for(k in 1:nyears){            # Loop over years
    for(i1 in 1:side){      # Loop over one side (hopefully X)
      for(i2 in 1:side){    # Loop over other side (hopefully Y)
        i1.start <- max(1,(i1-1))
        i1.end <- min(side,(i1+1))
        i2.start <- max(1,(i2-1))
        i2.end <- min(side,(i2+1))
        Xautoobs[i1,i2,k] <- (sum(zobs[i1.start:i1.end,i2.start:i2.end,k])- zobs[i1,i2,k])  / (length(zobs[i1.start:i1.end,i2.start:i2.end,k]) - 1)
      }
    }
  }

  out <- list(
    # ----------------- values input -----------------------
    side=side, nyears=nyears, nsurveys=nsurveys, mean.psi1=mean.psi1, beta.Xpsi1=beta.Xpsi1,
    range.phi=range.phi, beta.Xphi=beta.Xphi,
    range.gamma=range.gamma, beta.Xgamma=beta.Xgamma,
    range.p=range.p, beta.Xp=beta.Xp,
    theta.XAC=theta.XAC, beta.XAC= beta.XAC, beta.Xautolog=beta.Xautolog,
    trend.sd.site=trend.sd.site, trend.sd.survey=trend.sd.survey,
    seed=seed, seed.XAC = seed.XAC,
    # ----------------- values generated --------------------
    M=M,                # total number of pixels in the study area
    grid=grid,          # 2-column matrix, the x and y coordinates of the pixels
    amatrix = amatrix,  # MxM matrix, [i,j] = 1 if cell i and j are neighbours, 0 otherwise
    Xpsi1 = Xpsi1,      # side x side matrix, value of covariate affecting initial occupancy (psi1)
    Xphi = Xphi,        # side x side x nyears array, value of covariate affecting persistence (phi)
    Xgamma = Xgamma,    # side x side x nyears array, value of covariate affecting colonisation (gamma)
    Xp = Xp,            # side x side x nsurveys x nyears array, value of covariate affecting detection (p)
    XAC=XAC,            # side x side matrix, the spatially correlated covariate
    Xauto = Xauto,      # side x side x nyears array, autocovariate, proportion of neighbouring cells occupied
    Xautoobs = Xautoobs,# side x side x nyears array, observed autocovariate, proportion of neighbouring cells where species detected
    sd.site=sd.site,    # vector nyears, year-specific values of SD of Gaussian random site effects in p
    sd.survey=sd.survey,# vector nyears, year-specific values of SD of Gaussian random survey effects in p
    mean.phi=mean.phi,  # vector nyears-1, year-specific intercept of persistence on probability scale
    mean.gamma=mean.gamma,# vector nyears-1, year-specific intercept of colonisation on probability scale
    mean.p=mean.p,      # vector nyears, year-specific intercept of detection probability on probability scale
    psi=psi,            # side x side x nyears array, probability of occupancy of cell
    mean.psi=mean.psi,  # vector nyears, mean occupancy over all cells
    psi.app=psi.app,    # vector nyears, apparent occupancy, proportion of cells where species detected
    z=z,                # side x side x nyears array, true occupancy status of each cell in each year (1 if occupied)
    zobs=zobs,          # side x side x nyears array, observed occupancy status of each cell in each year (1 if detected)
    nocc = nocc,        # vector nyears, the true number of cells occupied each year
    nocc.obs = nocc.obs,# vector nyears, the number of cells where species detected each year
    phi=phi,            # side x side x nyears-1 array, probability of persistence in each interval between years
    gamma=gamma,        # side x side x nyears-1 array, probability of colonisation in each interval between years
    p=p,                # side x side x nsurveys x nyears array, probability of detection
    y = y)             # side x side x nsurveys x nyears array, detection history, 1 if species detected.

  # Add an unmarked data frame object
  out$umf <- conv2UM(out)

  return(out)
} # ------------------------------------------------------------------