R/sim_colext.r

In oliviergimenez/occuHMM: Fit occupancy data with hidden Markov models

#' Simulate data from a dynamic occupancy model
#'
#' This function allows you to simulate data from a dynamic occupancy model.
#' @param nb_sites number of sites (by default = 100)
#' @param nb_surveys number of surveys/visits (by default = 3)
#' @param nb_seasons number of seasons/years (by default = 10)
#' @param init_occ first-year occupancy probability (by default = 0.6)
#' @param det detection probability (by default = 0.7)
#' @param ext extinction probability (by default = 0.5)
#' @param col colonization probability (by default = 0.3)
#' @keywords multiseason
#' @export
#' @examples
#' vignette("colext_occuHMM")

sim_colext <- function(nb_sites=100,nb_surveys=3,nb_seasons=10, init_occ=0.6,det=0.7,ext=0.5,col=0.3){
# simulate data from dynamic occupancy model with constant parameters
# adapted from Kéry and Schaub book

# define various quantities
R = nb_sites # number of sites
J = nb_surveys # number of replicate surveys/visits
K = nb_seasons # number of years/seasons
psi1 = init_occ # occupancy prob in first year/season
p = det # detection prob
epsilon = ext # extinction prob
gamma = col # colonization prob

# gather parameters together
sim_param <- c(psi1,gamma,epsilon,p)

# pre-allocate memory
site <- 1:R # Sites
year <- 1:K # Years
psi <- rep(NA, K) # Occupancy probability
muZ <- z <- array(dim = c(R, K)) # Expected and realized occurrence
y <- array(NA, dim = c(R, J, K)) # Detection histories

# define state process
# first year/season
z[,1] <- rbinom(R, 1, psi1) # Initial occupancy state
# subsequent years/seasons
for(i in 1:R){ # Loop over sites
	for(k in 2:K){ # Loop over years
		muZ[k] <- z[i, k-1]*(1-epsilon) + (1-z[i, k-1])*gamma # Prob for occ.
		z[i,k] <- rbinom(1, 1, muZ[k])
		}
}

# define observation process
for(i in 1:R){
	for(k in 1:K){
		prob <- z[i,k] * p
		for(j in 1:J){
			y[i,j,k] <- rbinom(1, 1, prob)
			}
		}
}

# format data
yy <- matrix(y, R, J*K)
yy
}