Remove RandomFields from suggests

14 files changed, 755 insertions, 32 deletions

diff --git a/.local_checks/@making_AHMbook_package.R b/.local_checks/@making_AHMbook_package.R
index abd05ab..c60173c 100644
--- a/.local_checks/@making_AHMbook_package.R
+++ b/.local_checks/@making_AHMbook_package.R
@@ -28,7 +28,7 @@ unlink(list.files(pattern="Rplots.pdf", recursive=TRUE))
 system("R CMD build AHMbook")  # Produces the .tar.gz
 pkg <- "AHMbook_0.2.8.tar.gz"  # <-- fix version number here ################
 
-Sys.setenv("_R_CHECK_FORCE_SUGGESTS_" = FALSE) # Use if randomFields missing
+#Sys.setenv("_R_CHECK_FORCE_SUGGESTS_" = FALSE) # Use if randomFields missing
 
 # Pick one to check:
 ## on desktop
diff --git a/DESCRIPTION b/DESCRIPTION
index 8804d48..6cb761f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -5,7 +5,6 @@ Version: 0.2.9
 Date: 2023-08-21
 Depends: R (>= 2.10)
 Imports: grDevices, graphics, methods, stats, utils, plotrix, raster, sp, fields, coda, unmarked (>= 0.12.2), mvtnorm, spdep
-Suggests: RandomFields
 Authors@R: c(
   person("Marc", "Kéry", role="aut"),
   person("Andy", "Royle", role="aut"),
diff --git a/R/simDynoccSpatial_AHM2_9.R b/R/simDynoccSpatial_AHM2_9.R
index f49b377..8d71e2e 100644
--- a/R/simDynoccSpatial_AHM2_9.R
+++ b/R/simDynoccSpatial_AHM2_9.R
@@ -153,14 +153,14 @@ simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
   # 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 {
+  #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)
@@ -169,7 +169,7 @@ simDynoccSpatial <- function(side = 50, nyears = 10, nsurveys = 3,
     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
 
diff --git a/R/simExpCorrRF.R b/R/simExpCorrRF.R
index 225c263..6f20f8c 100644
--- a/R/simExpCorrRF.R
+++ b/R/simExpCorrRF.R
@@ -8,6 +8,7 @@
 # In DESCRIPTION file:
 # - add 'fields' to Imports
 # - move RandomFields from Imports to Suggests
+# - 8-21-2023: RandomFields completely removed as it does not seem to be coming back to CRAN
 
 # In NAMESPACE file:
 # - comment out or delete importFrom("RandomFields", "RFoptions", "RFsimulate", "RMexp")
@@ -36,12 +37,13 @@ y <- seq(1, size, step)
 # grid <- as.matrix(expand.grid(x,y))
 grid <- cbind(x = rep(x, each=size), y = y)
 
-if(requireNamespace("RandomFields", quietly=TRUE)) {
-  RandomFields::RFoptions(seed=seed)
-  field <- RandomFields::RFsimulate(RandomFields::RMexp(var = variance, scale = theta),
-      x=x, y=y, grid=TRUE)@data$variable1
-  RandomFields::RFoptions(seed=NA)
-} else {
+# RandomFields now completely removed from package
+#if(requireNamespace("RandomFields", quietly=TRUE)) {
+#  RandomFields::RFoptions(seed=seed)
+#  field <- RandomFields::RFsimulate(RandomFields::RMexp(var = variance, scale = theta),
+#      x=x, y=y, grid=TRUE)@data$variable1
+#  RandomFields::RFoptions(seed=NA)
+#} else {
  message("Using package 'fields' instead of 'RandomFields'; see help(simExpCorrRF).")
   if(!is.na(seed))
     set.seed(seed)  # Only for compatibility with RandomFields, better to set seed before calling simExpCommRF
@@ -50,7 +52,7 @@ if(requireNamespace("RandomFields", quietly=TRUE)) {
   if(inherits(tmp, "try-error"))
     stop("Simulation of random field failed.\nTry with smaller values for 'size' or 'theta'.")
   field <- as.vector(tmp * sqrt(variance))
-}
+#}
 
 # Plots
 # Correlation function
diff --git a/R/simPPe_AHM2_10.R b/R/simPPe_AHM2_10.R
index cd25af4..9e1d3d2 100644
--- a/R/simPPe_AHM2_10.R
+++ b/R/simPPe_AHM2_10.R
@@ -75,18 +75,18 @@ simPPe <- function(lscape.size = 150, buffer.width = 25, variance.X = 1, theta.X
   nsite <- length(mid.pt)^2
 
   # Simulate habitat covariate: a spatially correlated Gaussian random field (i.e., a Gaussian random field with negative exponential corr.)
-  if(requireNamespace("RandomFields", quietly=TRUE)) {
-    RandomFields::RFoptions(seed=NA)
-    field <- matrix(RandomFields::RFsimulate(RandomFields::RMexp(var = variance.X, scale = theta.X),
-      x=x, y=y, grid=TRUE)@data$variable1, ncol = lscape.size) # MVN r.v. with spatial correlation
-  } else {
+  #if(requireNamespace("RandomFields", quietly=TRUE)) {
+  #  RandomFields::RFoptions(seed=NA)
+  #  field <- matrix(RandomFields::RFsimulate(RandomFields::RMexp(var = variance.X, scale = theta.X),
+  #    x=x, y=y, grid=TRUE)@data$variable1, ncol = lscape.size) # MVN r.v. with spatial correlation
+  #} else {
     message("Using package 'fields' instead of 'RandomFields'; see help(simPPe).")
     obj <- circulantEmbeddingSetup(grid=list(x=x, y=y), Covariance="Exponential", aRange=theta.X)
     tmp <- try(circulantEmbedding(obj), silent=TRUE)
     if(inherits(tmp, "try-error"))
       stop("Simulation of random field failed.\nTry with smaller values for 'lscape.size' or 'theta.X'.")
     field <- matrix(tmp * sqrt(variance.X), ncol = lscape.size)
-  }
+  #}
 
   # --------------- Simulate points in the field --------------------
   #
diff --git a/inst/RandomFieldsSupport/simDynoccSpatial.R b/inst/RandomFieldsSupport/simDynoccSpatial.R
new file mode 100644
index 0000000..4bb53af
--- /dev/null
+++ b/inst/RandomFieldsSupport/simDynoccSpatial.R
@@ -0,0 +1,414 @@
+
+# 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)
+} # ------------------------------------------------------------------
diff --git a/inst/RandomFieldsSupport/simExpCorrRF.R b/inst/RandomFieldsSupport/simExpCorrRF.R
new file mode 100644
index 0000000..225c263
--- /dev/null
+++ b/inst/RandomFieldsSupport/simExpCorrRF.R
@@ -0,0 +1,81 @@
+# Define function for simulating spatially correlated random field
+# AHM2 - 9.2
+
+# Called by functions `simNmixSpatial` and `simOccSpatial`.
+
+# Modified to use package 'fields' if 'RandomFields' is not available.
+
+# In DESCRIPTION file:
+# - add 'fields' to Imports
+# - move RandomFields from Imports to Suggests
+
+# In NAMESPACE file:
+# - comment out or delete importFrom("RandomFields", "RFoptions", "RFsimulate", "RMexp")
+# - add: importFrom("fields", "circulantEmbeddingSetup", "circulantEmbedding")
+
+
+
+# ------------ Start of function definition ----------------
+simExpCorrRF <- function(variance = 1, theta = 1, size = 50, seed = NA, show.plots = TRUE){
+# Function creates Gaussian random field with negative
+# exponential correlation and visualizes correlation and random field
+#
+# Function arguments:
+# theta: parameter governing spatial correlation (=1/phi)
+# ("large theta means long range of correlation")
+# Note that RMexp is specified in terms of phi = 1/theta
+# variance: variance of field, set at 1
+# grid.size: Number of pixels in either direction
+# show.plot: if TRUE, plots of the data will be displayed;
+#  set to FALSE if you are running simulations or use inside of other fct's.
+
+# Generate correlated random variables in a square
+step <- 1
+x <- seq(1, size, step)
+y <- seq(1, size, step)
+# grid <- as.matrix(expand.grid(x,y))
+grid <- cbind(x = rep(x, each=size), y = y)
+
+if(requireNamespace("RandomFields", quietly=TRUE)) {
+  RandomFields::RFoptions(seed=seed)
+  field <- RandomFields::RFsimulate(RandomFields::RMexp(var = variance, scale = theta),
+      x=x, y=y, grid=TRUE)@data$variable1
+  RandomFields::RFoptions(seed=NA)
+} else {
+ message("Using package 'fields' instead of 'RandomFields'; see help(simExpCorrRF).")
+  if(!is.na(seed))
+    set.seed(seed)  # Only for compatibility with RandomFields, better to set seed before calling simExpCommRF
+  obj <- circulantEmbeddingSetup(grid=list(x=x, y=y), Covariance="Exponential", aRange=theta)
+  tmp <- try(circulantEmbedding(obj), silent=TRUE)
+  if(inherits(tmp, "try-error"))
+    stop("Simulation of random field failed.\nTry with smaller values for 'size' or 'theta'.")
+  field <- as.vector(tmp * sqrt(variance))
+}
+
+# Plots
+# Correlation function
+if(show.plots){
+  oldpar <- par(mfrow = c(1,2), mar = c(5,5,4,2), "cex.main") ; on.exit(par(oldpar))
+  tryPlot <- try( {
+    dis <- seq(0.01, 20, by = 0.01)
+    corr <- exp(-dis/theta)
+    plot(dis, corr, type = "l", xlab = "Distance", ylab = "Correlation", ylim = c(0,1), col = "blue", lwd = 2)
+    text(0.8*max(dis), 0.8, labels = paste("theta:", theta))
+
+    # Random field
+    # image(x, y, field,col=topo.colors(20), main = paste("Gaussian random field with \n negative exponential correlation (theta =", theta, ")"), cex.main = 1)
+    par(mar = c(3,2,5,1))
+    raster::plot(rasterFromXYZ(cbind(grid, field)), col=topo.colors(20),
+    main = paste("Gaussian random field with \n negative exponential correlation (theta =", theta, ")"), cex.main = 1, legend=FALSE, box=FALSE)
+    box()
+  }, silent = TRUE)
+  if(inherits(tryPlot, "try-error"))
+    tryPlotError(tryPlot)
+}
+
+# Output
+return(list(variance = variance, theta = theta, size = size, seed = seed,
+  field = field,
+  grid = grid))
+} # ------------ End of function definition ----------------
+
diff --git a/inst/RandomFieldsSupport/simPPe.R b/inst/RandomFieldsSupport/simPPe.R
new file mode 100644
index 0000000..cd25af4
--- /dev/null
+++ b/inst/RandomFieldsSupport/simPPe.R
@@ -0,0 +1,224 @@
+
+# AHM2 Chapter 10
+
+simPPe <- function(lscape.size = 150, buffer.width = 25, variance.X = 1, theta.X = 10, M = 250, beta = 1, quads.along.side = 6, show.plots = TRUE) {
+  #
+  # Name means 'SIMulate Point Pattern Educational version'
+
+  # Function to simulate spatial point pattern
+  # in a heterogeneous landscape simulated on a square grid.
+  # The study area ('core') is simulated inside of a larger landscape that
+  # includes a buffer. The size of the core is defined by the lscape.size
+  # minus twice the buffer.
+  # There is one habitat covariate X that affects the intensity
+  # of the points.
+  # X is spatially structured with negative exp. spatial autocorrelation;
+  # the parameters of the field can be chosen to create large
+  # 'islands' of similar values or no 'islands' at all,
+  # in which case the field is spatially unstructured.
+  #
+  # The intensity of STATIC points (e.g. home-range centers)
+  # may be inhomogeneous and affected by the coefficient beta,
+  # which is the log-linear effect of X.
+  #
+  # *** Function arguments ***
+  # lscape.size:  total side length of simulated landscape (length in units,
+  #               this includes a buffer)
+  # buffer.width: width of buffer around core study area
+  # variance.X:   variance of Gaussian random field (covariate X)
+  # theta.X:      scale parameter of correlation in field (must be >0)
+  # M:            Expected number of activity centers in core area
+  # beta:         coefficient of the habitat
+  # quads.along.side: number of quadrats along the side of the core area
+  #    (this is the parameter of the gridding process
+  #    and determines the size of the quadrats)
+
+  # -------------- Check and fix input -----------------------
+  buffer.width <- round(buffer.width[1])
+  stopifNegative(buffer.width)
+  quads.along.side <- round(quads.along.side[1])
+  stopifNegative(quads.along.side, allowZero=FALSE)
+  lscape.size <- round(lscape.size[1])
+  stopifnotGreaterthan(lscape.size, 2 * buffer.width + quads.along.side - 1)
+  variance.X <- variance.X[1]
+  stopifNegative(variance.X)
+  theta.X <- theta.X[1]
+  stopifNegative(theta.X)
+  M <- round(M[1])
+  stopifNegative(M)
+  # ------------------------------------------------------------
+
+  # --------------- Define basic geometry of simulation --------------
+  #
+  # Size of core study area (the 'core') and its proportion of total landscape area
+  size.core <- lscape.size - 2 * buffer.width
+  prop.core <- size.core^2/lscape.size^2 # ratio core area / total area
+
+  # Discrete approximation of total landscape
+  pixel.size <- 1      # length of side of square pixel of simulated landscape
+
+  # Coordinates (mid-points of basic pixel unit of simulation) and lscape
+  x <- seq(1, lscape.size, pixel.size)-0.5    # x coordinate of pixels
+  y <- seq(1, lscape.size, pixel.size)-0.5    # y coordinate of pixels
+  grid <- as.matrix(expand.grid(x,y))       # resulting grid
+
+  # Compute lambda of point pattern: limit of expected number of points per areal unit, when latter goes towards zero
+  lambda_pp <- M / size.core^2
+
+  # Define a core area in the middle of the square
+  # This core is then divided up to define a number of quadrats
+  # within which abundance and occurrence are measured
+  quad.size <- size.core / quads.along.side
+  breaks <- seq(buffer.width, size.core+buffer.width, by = quad.size) # boundaries of quadrats
+  mid.pt <- breaks[-length(breaks)] + 0.5 * quad.size # quadrat mid-points
+  core <- range(breaks)   # range of x and y coordinates in the core
+  nsite <- length(mid.pt)^2
+
+  # Simulate habitat covariate: a spatially correlated Gaussian random field (i.e., a Gaussian random field with negative exponential corr.)
+  if(requireNamespace("RandomFields", quietly=TRUE)) {
+    RandomFields::RFoptions(seed=NA)
+    field <- matrix(RandomFields::RFsimulate(RandomFields::RMexp(var = variance.X, scale = theta.X),
+      x=x, y=y, grid=TRUE)@data$variable1, ncol = lscape.size) # MVN r.v. with spatial correlation
+  } else {
+    message("Using package 'fields' instead of 'RandomFields'; see help(simPPe).")
+    obj <- circulantEmbeddingSetup(grid=list(x=x, y=y), Covariance="Exponential", aRange=theta.X)
+    tmp <- try(circulantEmbedding(obj), silent=TRUE)
+    if(inherits(tmp, "try-error"))
+      stop("Simulation of random field failed.\nTry with smaller values for 'lscape.size' or 'theta.X'.")
+    field <- matrix(tmp * sqrt(variance.X), ncol = lscape.size)
+  }
+
+  # --------------- Simulate points in the field --------------------
+  #
+  # Simulate binomial point process for activity centers
+  M2 <- round(M/prop.core) # Number of individuals in the total landscape
+                           #   (incl. buffer)
+
+  # Simulate point locations as function of habitat covariate x
+  probtemp <- exp(beta[1]*c(field)) # log-linear model for intensity on X
+  probs <- probtemp / sum(probtemp) # normalize to get probability of getting a point in a pixel
+  pixel.id <- sort(sample(1:lscape.size^2, M2 , replace=TRUE, prob=probs))
+
+  # Simulate locations randomly within the pixel (unlike sim.spatialDS)
+  u1 <- grid[pixel.id,1] + runif(M2, -pixel.size/2, pixel.size /2)
+  u2 <- grid[pixel.id,2] + runif(M2, -pixel.size/2, pixel.size /2)
+  u <- cbind(u1, u2)    # collect AC coordinates together
+
+
+  # ------ Summarization of point pattern within quadrats -------------
+  #
+  # This is INSIDE of the observation core
+  #
+  # Summarization for abundance (N) at every quadrat
+  Nac <- as.matrix(table(cut(u[,1], breaks=breaks),
+     cut(u[,2], breaks= breaks))) # quadrat-specific abundance for AC
+  E_N <- round(mean(Nac),2)       # average realized abundance per quadrat
+
+  # Summarization for presence/absence (z) at every quadrat
+  zac <- Nac   ;  zac[zac>1] <- 1 # quadrat-specific occurrence for AC
+  E_z <- round(mean(zac), 2)      # proportion occupied quadrats
+
+  # ------------------ Visualizations ---------------------------
+  if(show.plots) {
+    oldpar <- par(mfrow = c(1, 3), mar = c(4,2,5,2), cex.main = 1.8, cex.axis = 1.2) ; on.exit(par(oldpar))
+
+    tryPlot <- try( {
+      # *** Fig. 1: Original point pattern
+      # Random field of X with activity-centers overlaid
+      image(rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10),
+          main = "Point pattern with\ncore and buffer area",
+          xlab = "", ylab = "", axes = FALSE, asp = 1)
+      mtext(paste("Mean intensity (lambda) =", round(lambda_pp, 5)), side=1)
+      polygon(c(buffer.width, size.core+buffer.width, size.core+buffer.width, buffer.width),
+            c(buffer.width, buffer.width, size.core+buffer.width, size.core+buffer.width), 
+            lwd = 2, lty = 1)
+      points(u[,1], u[,2], pch=20, col='black', cex = 1.2)  # plot points
+      # points(u1, u2, pch=20, col='black', cex = 1.2)  # plot points
+
+      # *** Fig. 2: Show abundance and presence/absence in each quadrat on original landscape ***
+      # Covariate 1: the Gaussian random field with autocorrelation
+      # Reproduce random field with activity centers
+      image(rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10), main = "Abundance, N",
+          xlab = "", ylab = "", axes = FALSE, asp = 1)
+      mtext(paste0("Mean(N) = ", E_N, ", var(N) = ", round(var(c(Nac)), 2)), side=1)
+      polygon(c(buffer.width, size.core+buffer.width, size.core+buffer.width, buffer.width),
+          c(buffer.width, buffer.width, size.core+buffer.width, size.core+buffer.width), lwd = 2, lty = 1)
+      # Add activity centers
+      points(u[,1], u[,2], pch=20, col='black', cex = 1.2)  # plot points
+      # Overlay survey quadrats
+      for(i in 1:length(breaks)){
+         for(k in 1:length(breaks)){
+           segments(breaks[i], breaks[k], rev(breaks)[i], breaks[k])
+           segments(breaks[i], breaks[k], breaks[i], rev(breaks)[k])
+         }
+      }
+      # Print abundance into each quadrat
+      for(i in 1:length(mid.pt)){
+        for(k in 1:length(mid.pt)){
+          text(mid.pt[i], mid.pt[k], Nac[i,k], cex =4^(0.8-0.5*log10(quads.along.side)), col='red')
+        }
+      }
+
+      # Figure 3 for presence/absence of activity centers (= distribution)
+      # Reproduce random field with activity centers
+      image(rasterFromXYZ(cbind(grid, c(field))), col=topo.colors(10), main = "Occurrence, z",
+          xlab = "", ylab = "", axes = FALSE, asp = 1)
+      mtext(paste("Mean(z) =", E_z), side=1)
+      polygon(c(buffer.width, size.core+buffer.width, size.core+buffer.width, buffer.width),
+          c(buffer.width, buffer.width, size.core+buffer.width, size.core+buffer.width), lwd = 2, lty = 1)
+      # Add activity centers
+      points(u[,1], u[,2], pch=20, col='black', cex = 1.2)  # plot points
+      # Overlay quadrats
+      for(i in 1:length(breaks)){
+         for(k in 1:length(breaks)){
+           segments(breaks[i], breaks[k], rev(breaks)[i], breaks[k])
+           segments(breaks[i], breaks[k], breaks[i], rev(breaks)[k])
+         }
+      }
+      # Print presence/absence into each quadrat
+      for(i in 1:length(mid.pt)){
+        for(k in 1:length(mid.pt)){
+          text(mid.pt[i], mid.pt[k], zac[i,k], cex =4^(0.8-0.5*log10(quads.along.side)), col='red')
+        }
+      }
+
+      # Mike: Shade UNoccupied quadrats (which have abundance N = 0 or occurrence z = 0)
+      for(i in 1:(length(breaks)-1)){
+        for(k in 1:(length(breaks)-1)){
+          if(zac[i,k] == 1) # grey-out UNoccupied quads
+            next
+          polygon(c(breaks[i], breaks[i+1], breaks[i+1], breaks[i]),
+              c(breaks[k], breaks[k], breaks[k+1], breaks[k+1]),
+              col = adjustcolor("black", 0.6))
+        }
+      }
+    }, silent = TRUE)
+    if(inherits(tryPlot, "try-error"))
+      tryPlotError(tryPlot)
+  }
+  
+# Numerical output
+return(list(
+  # ----------------- arguments input -----------------------
+  grid.size = lscape.size, buffer.width = buffer.width, variance.X = variance.X, theta.X = theta.X, M = M, beta = beta, quads.along.side = quads.along.side,
+  # ---------------- generated values -------------------------
+  core = core,            # range of x and y coordinates in the 'core'
+  M2 = M2,                # Number of ACs in the total landscape (incl. buffer)
+  grid = grid,            # Coordinates of the centre of each pixel.
+  pixel.size = pixel.size,# 1; length of side of square pixel of landscape
+  size.core = size.core,  # the width = height of the core area
+  prop.core = prop.core,  # proportion of the landscape in the core
+  X = field,              # lscape.size x lscape.size matrix of covariate values for each pixel
+  probs = probs,          # corresponding matrix of probability of AC in pixel (sums to 1)
+  pixel.id = pixel.id,    # M2 vector, which pixel each AC is inside.
+  u = u,                  # M2 x 2 matrix, coordinates of each AC
+  nsite = nsite,          # total number of quadrats in the core
+  quad.size = quad.size,  # width = height of each quadrat
+  breaks = breaks,        # boundaries of each quadrat
+  mid.pt = mid.pt,        # mid=points of each quadrat
+  lambda_pp = lambda_pp,  # intensity of point pattern (ACs per unit area)
+  Nac = Nac,    # matrix, quads.along.side x quads.along.side, site-specific abundance of ACs
+  zac = zac,    # matrix, quads.along.side x quads.along.side, 0/1 occurrence
+  E_N = E_N,    # scalar, average realized abundance per quadrat.
+  E_z = E_z))   # scalar, average realized occupancy per quadrat.
+} # ------------ End of function definition --------------------
diff --git a/man/data_simDynoccSpatialData.Rd b/man/data_simDynoccSpatialData.Rd
index 8889feb..c7cb17d 100644
--- a/man/data_simDynoccSpatialData.Rd
+++ b/man/data_simDynoccSpatialData.Rd
@@ -14,7 +14,9 @@ The data analyzed in Kéry & Royle (2021) section 9.6.1.2 was generated using \c
 \details{
 The data set was generated on a machine with \pkg{RandomFields} installed, using the following code:
 
-\code{simDynoccSpatialData <- simDynoccSpatial(side = 30, nyears = 10, nsurveys = 3, mean.psi1 = 0.1,
+\code{
+  source(system.file("RandomFieldsSupport", "simDynoccSpatial.R", package="AHMbook"))
+  simDynoccSpatialData <- simDynoccSpatial(side = 30, nyears = 10, nsurveys = 3, mean.psi1 = 0.1,
   range.phi = c(0.5, 0.5), range.gamma = c(0.2, 0.2), range.p = c(0.4, 0.4),
   beta.Xautolog = c(1, 1), seed.XAC = 1, seed = 24, ask.plot = TRUE)}
 }
diff --git a/man/simDynoccSpatial.Rd b/man/simDynoccSpatial.Rd
index f47376d..ce2d8c6 100644
--- a/man/simDynoccSpatial.Rd
+++ b/man/simDynoccSpatial.Rd
@@ -21,7 +21,7 @@ Function to simulate detection/nondetection data under a general dynamic site-oc
 
 * 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.
 
-Previous versions used \pkg{RandomFields}, but that is not currently available on CRAN. If not available, \pkg{fields} will be used, but it cannot deal with large values of \code{side} and \code{theta.XAC}.
+Previous versions used \pkg{RandomFields}, but that is not currently available on CRAN. \pkg{fields} is now used instead, but it cannot deal with large values of \code{side} and \code{theta.XAC}. If you have \pkg{RandomFields} installed (perhaps by getting it from the CRAN archive), you can load a version of \code{simDynoccSpatial} that supports it with \code{source(system.file("RandomFieldsSupport", "simDynoccSpatial.R", package="AHMbook"))}.
 
 See also \code{\link{simDynoccSpatialData}}.
 }
diff --git a/man/simExpCorrRF.Rd b/man/simExpCorrRF.Rd
index 36d003f..e7fa0a7 100644
--- a/man/simExpCorrRF.Rd
+++ b/man/simExpCorrRF.Rd
@@ -8,7 +8,7 @@ Create a Gaussian random field with negative exponential correlation
 \description{
 Function creates Gaussian random field with negative exponential correlation and visualizes correlation and random field.
 
-Previous versions used \pkg{RandomFields}, but that is not currently available on CRAN. If not available, \pkg{fields} will be used, but it cannot deal with large values of \code{size} and \code{theta}.
+Previous versions used \pkg{RandomFields}, but that is not currently available on CRAN. \pkg{fields} is now used instead, but it cannot deal with large values of \code{size} and \code{theta}. If you have \pkg{RandomFields} installed (perhaps by getting it from the CRAN archive), you can load a version of \code{simExpCorrRF} that supports it with \code{source(system.file("RandomFieldsSupport", "simExpCorrRF.R", package="AHMbook"))}.
 }
 \usage{
 simExpCorrRF(variance = 1, theta = 1, size = 50, seed = NA, show.plots = TRUE)
diff --git a/man/simNmixSpatial.Rd b/man/simNmixSpatial.Rd
index 12a763a..b62fb15 100644
--- a/man/simNmixSpatial.Rd
+++ b/man/simNmixSpatial.Rd
@@ -137,7 +137,7 @@ str(dat<- simNmixSpatial(variance.RF = 0))
 str(dat<- simNmixSpatial(theta.RF = 0.0001))
 
 # Larger value of theta.RF gives larger 'islands'
-try(str(dat<- simNmixSpatial(theta.RF = 100))) # works with RandomFields
+#try(str(dat<- simNmixSpatial(theta.RF = 100))) # works with RandomFields
 
 # Truncate abundance in final plots to presence/absence
 str(dat<- simNmixSpatial(truncN = 0.5))
diff --git a/man/simOccSpatial.Rd b/man/simOccSpatial.Rd
index 8138981..5e89082 100644
--- a/man/simOccSpatial.Rd
+++ b/man/simOccSpatial.Rd
@@ -137,6 +137,6 @@ str(dat <- simOccSpatial(variance.RF = 0))
 str(dat <- simOccSpatial(theta.RF = 0.0001))
 
 # Larger value of theta.RF gives larger 'islands'
-try(str(dat <- simOccSpatial(theta.RF = 100)))  # Works with RandomFields
+#try(str(dat <- simOccSpatial(theta.RF = 100)))  # Works with RandomFields
 }
 }
diff --git a/man/simPPe.Rd b/man/simPPe.Rd
index d0b6cab..2e43932 100644
--- a/man/simPPe.Rd
+++ b/man/simPPe.Rd
@@ -14,7 +14,7 @@ The intensity of STATIC points (e.g. animal activity centers) may be inhomogeneo
 
 To recreate the data sets used in the book with R 3.6.0 or later, include \code{sample.kind="Rounding"} in the call to \code{set.seed}. This should only be used for reproduction of old results.
 
-Previous versions used \pkg{RandomFields}, but that is not currently available on CRAN. If not available, \pkg{fields} will be used, but it cannot deal with large values of \code{lscape.size} and \code{theta.X}.
+Previous versions used \pkg{RandomFields}, but that is not currently available on CRAN. \pkg{fields} is now used, but it cannot deal with large values of \code{lscape.size} and \code{theta.X}. If you have \pkg{RandomFields} installed (perhaps by getting it from the CRAN archive), you can load a version of \code{simPPe} that supports it with \code{source(system.file("RandomFieldsSupport", "simPPe.R", package="AHMbook"))}.
 
 }
 \usage{
@@ -82,8 +82,9 @@ Marc Kéry & Andy Royle.
 # Nice plot (produces the really nice Fig. 10-2 in the book)
 # RNGkind(sample.kind = "Rounding") # run this for R >= 3.6.0
 set.seed(117, kind="Mersenne-Twister")
-try(str(dat <- simPPe(lscape.size = 200, buffer.width = 25, variance.X = 1,
-  theta.X = 70, M = 200, beta = 1, quads.along.side = 6)))  # Fails if RandomFields is not available
+# Fails if RandomFields is not available
+#try(str(dat <- simPPe(lscape.size = 200, buffer.width = 25, variance.X = 1,
+#  theta.X = 70, M = 200, beta = 1, quads.along.side = 6)))
 
 str(dat <- simPPe(lscape.size = 200, buffer.width = 20, variance.X = 1,
   theta.X = 5, M = 250, beta = 1, quads.along.side = 6))