### Simulation study to compare the observation confirmation and site confirmation models
# rm(list = ls())

n.sims <- 1000			## Number of simulations (Be Careful: 1,000 simulations takes several hours)

### PARAMETERS
psi.list <- c(0.1, 0.9)		## Occupancy probability 
s11.list <- c(0.1, 0.9)		## Probability of >=1 observation-level true detection
s10.list <- c(0.1, 0.7)		## Probability of >=1 observation-level false positive

### SAMPLE CHARACTERSITICS
 K <- 30					## Number of survey occasions
 R <- 5000					## Number of sites
prop.conf.list <- c(0.05, 0.5)	## Proportion of sites where observations are confirmed


### Create a table to store simulation results
 TABLE <- matrix(NA, (2^4), 6)
colnames(TABLE) <- c("bias.psi.MM" , "bias.psi.OC" ,
				"error.psi.MM" , "error.psi.OC" ,
				"rmse.psi.MM" , "rmse.psi.OC" )

### Function for transformations
expit <- function(x) {1/(1 + exp(-x))}
logit <- function(p) {log(p/(1-p))}

### Likelihoods ###
##-------------------------------------------------##
###  Site  Confirmation Model -- Multiple Methods ###
##-------------------------------------------------##
SiteConf.MM <- function(beta,data,occ){
	psi = expit(beta[1])
	p11 = expit(beta[2])
	p10 = expit(beta[3])
 	r11 = expit(beta[4])			
	y = data[[1]]
	w = data[[2]]
 nLL = sum( -log( psi*dbinom(y,occ,p11) + (1-psi)*dbinom(y,occ,p10) ) ) +
		  sum( -log( psi*dbinom(w,occ,r11) + (1-psi)*dbinom(w,occ,0) ) )
return(nLL)
}

##----------------------------##
###  Obs. Confirmation Model ###
##----------------------------##
ObsConf <- function(beta,data,occ){
	psi = expit(beta[1])
	s11 = expit(beta[2])
	s10 = expit(beta[3])
  p11 <- s11 + s10 - (s11*s10)
  p10 <- s10
pr.X0.z1 <- (1-s11-s10+s11*s10)
pr.X1.z1 <- (s10-s11*s10)
pr.X2.z1 <- (s11-s11*s10)
pr.X3.z1 <- (s11*s10)
 pr.X0.z0 <- (1-s10)
 pr.X1.z0 <- s10
 pr.X2.z0 <- 0
 pr.X3.z0 <- 0
   pi.z1 <- c( pr.X0.z1 , pr.X1.z1 , pr.X2.z1 , pr.X3.z1 )  ;  sum(pi.z1)
   pi.z0 <- c( pr.X0.z0 , pr.X1.z0 , pr.X2.z0 , pr.X3.z0 )  ;  sum(pi.z0)
	y = data[[1]]
	w = data[[2]]
 nLL = sum( -log( psi*(p11^y)*((1-p11)^(occ-y)) + (1-psi)*(p10^y)*((1-p10)^(occ-y)) ) ) + 
		  sum( -log( psi*apply(t(pi.z1^t(w)), 1, prod) +  (1-psi)*apply(t(pi.z0^t(w)), 1, prod)  ) )
return(nLL)
}

###--------------------------------------------------------------------------------###

ct <- 0				## counter variable to keep track of iterations
time.start <- Sys.time()	## Start time

###---------------------------- START THE ITERATIONS ------------------------------###
###--------------------------------------------------------------------------------###

for (a in 1:length(psi.list)) 	{
for (b in 1:length(s11.list))		{
for (c in 1:length(s10.list)) 	{
for (f in 1:length(prop.conf.list)) {

  ct <- ct +1

### PARAMETERS
psi <- psi.list[a]
s11 <- s11.list[b]
s10 <- s10.list[c]

prop.conf <- prop.conf.list[f]	## Proportion of sites where data are confirmed (flase or true sample)

  I <- R*(1-prop.conf)		## nb of sites where cues where NOT confirmed
  J <- R*(prop.conf)		## nb of sites where cues where CONFIRMED 


###--------------------------------------------------------------------------------###

### POPULATION PARAMETERS

  p11 <- s11 + s10*(1-s11) ; p11
  p10 <- s10 ; p10

### Objects to save results MLE
 psi.MLE.MM <- NA
 p11.MLE.MM <- NA
 p10.MLE.MM <- NA

 psi.MLE.OC.Bin <- NA 
 s11.MLE.OC.Bin <- NA
 s10.MLE.OC.Bin <- NA

 psi.MLE.OC <- NA 
 s11.MLE.OC <- NA
 s10.MLE.OC <- NA

###--------------------------------------------------------------------------------###
###---------------------------- START THE SIMULATIONS -----------------------------###
###--------------------------------------------------------------------------------###

	for (s in 1:n.sims){

###----------------  SIMULATE THE DATA  ----------------###
### OCCUPANCY STATES
Z1 <- rbinom(I, 1, psi) ; Z1 ; length(Z1)
Z2 <- rbinom(J, 1, psi) ; Z2 ; length(Z2)

### DATA Y = non-confirmed data for I sites
 Y <- matrix(NA,I,K)
  for (i in 1:I){
 p <- Z1[i]*p11 + (1-Z1[i])*p10
Y[i,] <- rbinom(K, 1, p)	
   }							; head(Z1) ;  head(Y) ; dim(Y)
Y.sum <- apply(Y, 1, sum) ; head(Y.sum)

### Data W = confirmed data for J sites
 W.sum <- matrix(NA,J,4)
  for (j in 1:J){
pr.X0 <- Z2[j]*(1-s11-s10+s11*s10) + (1-Z2[j])*(1-s10)
pr.X1 <- Z2[j]*(s10-s11*s10) + (1-Z2[j])*s10
pr.X2 <- Z2[j]*(s11-s11*s10)
pr.X3 <- Z2[j]*(s11*s10)
  pi <- c( pr.X0 , pr.X1 , pr.X2 , pr.X3 )  ;  sum(pi)
W.sum[j,] <- as.vector(rmultinom(1,K,pi))
} 							; head(Z2)	;  head(W.sum) ; dim(W.sum)

###--------------------------------------------------------------------------------###
###--------------------------------------------------------------------------------###

## Site Confirmation Model -- Multiple METHODS ##
### First, 'simplify' the confirmed data for the SiteConf approach
 W.sum.MM <- W.sum
 W.sum.MM[,3] <- W.sum.MM[,3] + W.sum.MM[,4]  # Data CONFIRMED as TP <=> certain detection
 W.sum.MM <- W.sum.MM[,3]
data.MM <- list(Y.sum, W.sum.MM) ; length(data.MM)
beta.MM <- c(logit(psi), logit(p11), logit(p10), logit(s11)) + rnorm(4,0,sd=0.05)

  out.MM = optim(beta.MM, SiteConf.MM, data = data.MM, occ = K, hessian = T)  

est.MM <- round(expit(out.MM$par), 4)	

 psi.MLE.MM[s] <- est.MM[1]

###--------------------------------------------------------------------------------###
###--------------------------------------------------------------------------------###
## Obs. Confirmation model ##
data.OC <- list(Y.sum, W.sum)
beta.OC <- c(logit(psi), logit(s11), logit(s10)) + rnorm(3,0,sd=0.05) 

  out.OC = optim(beta.OC, ObsConf, data = data.OC, occ = K, hessian = T)  

est.OC <- round(expit(out.OC$par), 4)		

 psi.MLE.OC[s] <- est.OC[1]

##---------------------------------##
#  TIME and PROGRESS of SIMULATION  #
##---------------------------------##
 tot.sim <- (2^6)*n.sims
 avc <- (ct-1)*n.sims + s
if( round(avc/tot.sim, 3)-(avc/tot.sim) ==  0 ){ 
   print( paste((avc/tot.sim)*100, "% completed", sep = "") )
}else{}
flush.console()
Sys.sleep(0)

 } # end of loop over 's' -- simulations

###--------------------------------------------------------------------------------###
###----------------------------  END OF SIMULATIONS  ------------------------------###
###--------------------------------------------------------------------------------###

### Measures of estimator Performance -- Model Site Conf MM
 bias.psi.MM <- mean(psi.MLE.MM) - psi 
 error.psi.MM <- sqrt( sum( (psi.MLE.MM - mean(psi.MLE.MM))^2 ) / (length(psi.MLE.MM)-1) ) 
mse.psi.MM <- mean( (psi.MLE.MM - psi)^2 ) 
 rmse.psi.MM <- sqrt(mse.psi.MM) 

### Measures of estimator Performance -- Model Obs Conf (OC)
 bias.psi.OC <- mean(psi.MLE.OC) - psi 
 error.psi.OC <- sqrt( sum( (psi.MLE.OC - mean(psi.MLE.OC))^2 ) / (length(psi.MLE.OC)-1) )  
mse.psi.OC <- mean( (psi.MLE.OC - psi)^2 ) 
 rmse.psi.OC <- sqrt(mse.psi.OC) 

###--------------------------------------------------------------------------------###
TABLE[ct,] <- cbind(round(bias.psi.MM, 2), round(bias.psi.OC, 2), 
			round(error.psi.MM, 2), round(error.psi.OC, 2),
			round(rmse.psi.MM, 2), round(rmse.psi.OC, 2)
			)

}}} } ## END of ITERATIONS

### GET THE DURATION THAT IT TAKES TO SIMULATIONS TO COMPLETE
time.end <- Sys.time()
duration <- time.end - time.start 
duration

### LOOK AT RESULTS and the duration of simulations
TABLE
duration

### Save the rsults
# Create columns for variable parameters values
col <- matrix(NA, 2^4, 4)
colnames(col) <- c("psi","s11","s10","prop.conf")
col[,4] <- rep(prop.conf.list,2^3); length(col[,4])
col[,3] <- rep(c(rep(s10.list[1],2^1),rep(s10.list[2],2^1)) , 2^2) ; length(col[,3])
col[,2] <- rep(c(rep(s11.list[1],2^2),rep(s11.list[2],2^2)) , 2^1) ; length(col[,2])
col[,1] <- rep(c(rep(psi.list[1],2^3),rep(psi.list[2],2^3)) , 2^0) ; length(col[,1])

# Add the simulation resutl
TABLE.2 <- as.matrix(cbind(col, TABLE))
head(TABLE.2)

# Save the result as a .csv file
write.table(TABLE.2, "ResultSimul.csv", row.names = F, sep = ";")

###---------------------------------###   END of FILE   ###---------------------------------###