### SIMULATE a DATASET under the OBSERVATION CONFIRMATION DESIGN
rm(list = ls())	# Clean the workspace
set.seed(367)

### PARAMETERS
psi <- 0.6	# occupancy probability 
s1 <- 0.7	# probability of >=1 observation-level true detection (unconditional on site status)
s0 <- 0.2	# probability of >=1 observation-level false positive (unconditional on site status)

  p11 <- s1 + s0 - (s1*s0)	# Define site-level detection probability
  p10 <- s0				# Define observation-level false positive probability

### SAMPLE CHARACTERSITICS
K <- 5	# number of survey occasions
I <- 500	# number of sites where observations are not confirmed
J <- 300	# number of sites where observations confirmed
R <- I + J	# number total of sites

### SIMULATE OCCUPANCY STATES
Z1 <- rbinom(I, 1, psi)	  # for sites I sites where observations are not confirmed
Z2 <- rbinom(J, 1, psi)   # for sites J sites where observations are confirmed

### SIMULATE AMBIGUOUS-DETECTION/NON-DETECTION DATA Y (site-level data, for sites I sites where observations are not confirmed)
 Y <- matrix(NA,I,K)			### Create a matrix to store detection/non-detection data 
  for (i in 1:I){				### Loop over all I sites 'i'
 p <- Z1[i]*p11 + (1-Z1[i])*p10	### 'Binomial' parameter, conditional on site occupancy status 
Y[i,] <- rbinom(K, 1, p)		### Site-level data Y [0/1]
   }						### End Loop over sites
Y.sum <- apply(Y, 1, sum) 		### Summarize the data as 'number of detection for each site'

### SIMULATE Confirmed Observation (for sites J sites)
 V <- matrix(NA,J,4)					### Create a matrix to store the confirmed observation data 
  for (j in 1:J){							### Loop over all J sites 'j'
pr.X0 <- Z2[j]*(1-s1-s0+s1*s0) + (1-Z2[j])*(1-s0)	### Define the probability of observation-state v = 0 (non-detection)
pr.X1 <- Z2[j]*(s0-s1*s0) + (1-Z2[j])*s0			### Define the probability of observation-state v = 1 (observation-level false positives only)
pr.X2 <- Z2[j]*(s1-s1*s0)					### Define the probability of observation-state v = 2 (observation-level true detections only)
pr.X3 <- Z2[j]*(s1*s0)						### Define the probability of observation-state v = 3 (both true and false observation-level detections were found)
 pi <- c(pr.X0,pr.X1,pr.X2,pr.X3)				### Define the vector of probabilities of the Categorical distribution
V[j,] <- as.vector(rmultinom(1,K,pi))			### Site-level data V [mutli-state]
}									### End Loop over sites

### Save the data used in the model
data.ObsConf <- list(Y.sum, V) 		### Observation DATA
occ.ObsConf <- K					### number of survey occasions
 TrueValues.ObsConf <- cbind(psi, s1, s0)	### Save the true parameter values for comparison wih estimates

# Save these obsjects as a 'RData' file
save(file = "DataObsConf.RData", list = c("data.ObsConf", "occ.ObsConf", "TrueValues.ObsConf"))