R/sim_colextfp.r

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

#' Simulate data from a dynamic occupancy model with false positives
#'
#' This function allows you to simulate data from a dynamic occupancy model with false positives.
#' Data coding is as follows: 0 is for non-detection, 1 is for detection with no ambiguity, 2 is for detection with false positives possibly occurring.
#' @param nb_sites number of sites (by default = 200)
#' @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 falsedet probability of false positive detection (by default = 0.1)
#' @param truedet probability of detection (by default = 0.4)
#' @param assign probability that a detection is classified as unambiguous (by default = 0.8)
#' @param ext extinction probability (by default = 0.5)
#' @param col colonization probability (by default = 0.3)
#' @keywords multiseason false positives
#' @export
#' @examples
#' vignette("colext_occuHMM")

sim_colextfp <- function(nb_sites=200,nb_surveys=3,nb_seasons=10,init_occ=0.6,falsedet=0.1,truedet=0.4,assign=0.8,ext=0.5,col=0.3){

# simulate data from dynamic occupancy model with constant parameters

# define various quantities
R = nb_sites
J = nb_surveys
K = nb_seasons
psi1 = init_occ # occupancy probability in first year/season
p10 = falsedet # p10
p11 = truedet # p11
delta = assign # delta
epsilon = ext # extinction
gamma = col # colonization

# gather parameters together
sim_param <- c(psi1,gamma,epsilon,p10,p11,delta)

# 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){ # loop sites
      for(k in 1:K){ # loop years
      	  for(j in 1:J){ # loop replicate surveys
				# if site unoccupied
      	  		if (z[i,k] == 0){
				   # step1: detection
					y_temp <- rbinom(1,1,p10)
      	  		   # step2: assignment
      	  		    if (y_temp == 0) y[i,j,k] <- 0
      	  		    if (y_temp == 1) y[i,j,k] <- 2
      	  		    } # end if
				# if site occupied
      	  		if (z[i,k] == 1){
				   # step1: detection
					y_temp <- rbinom(1,1,p11)
      	  		   # step2: assignment
      	  		    if (y_temp == 0) y[i,j,k] <- 0
      	  		    if (y_temp == 1) y[i,j,k] <- (1-rbinom(1,1,delta)) + 1
      	  		    } # end if
      	  		  } # end loop replicated surveys
      	  	} # end loop years
     } # end loop sites

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