### Investigate likelihood convergence ###

		###---------------------------------------------------------------###
		###		STEP 1: Simulate Data (Obs. Confirmation Design)	###
		###---------------------------------------------------------------###

# rm(list = ls())
set.seed(1584)

### PARAMETERS
psi <- 0.8	# occupancy probability 
s1 <- 0.7	# probability of >=1 observation-level true detection (unconditional on site status)
s0 <- 0.1	# 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 <- 250	# number of sites where observations confirmed
R <- I + J	# number total of sites

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

### 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

# Organize the confirmed data for the Site Confirmation Model
V.MM <- V[,3] + V[,4]  # Group all observations CONFIRMED as True Positives <=> Unambiguous detections


		###---------------------------------------------------------------###
		###		STEP 2: Implement the Site Confirmation Model		###
		###---------------------------------------------------------------###

### Site Confirmation Model  -- Multiple Detection Methods
## Funtcions
expit <- function(x) {1/(1 + exp(-x))}
logit <- function(p) {log(p/(1-p))}

## Likelihood
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]]
	v = data[[2]]
 nLL = sum( -log( psi*dbinom(y,occ,p11) + (1-psi)*dbinom(y,occ,p10) ) ) +
		  sum( -log( psi*dbinom(v,occ,r11) + (1-psi)*dbinom(v,occ,0) ) )
return(nLL)
}

## Data
dataSC <- list(Y.sum, V.MM)
occSC <- K

####### COMPARISON #######
### 1. Initial values : SET 1 (close to real values)
betaSC <- c(1, 1, -2, 1)

### RUN the Model
  outSC.MM = optim(betaSC, SiteConf.MM, data = dataSC, occ = occSC, hessian = T)  
### Comparison real values and estimated values 
 round(expit(betaSC), 4)			### Initial values for (in order) psi, p11, p10, r11
 round(expit(outSC.MM$par), 4)		### Estimated parameters are (in order) psi, p11, p10, r11
 c(psi, p11, p10, s1)				### Real values -- Note that r11 = s1 (= Probability of observing a true positive, given that the site is occupied)

### Deviance: -2LogLik
outSC.MM$value*2
###-------------------------------------------------------------------------------------------------------###
### NOTE: Estimated and real values are pretty close -- Seems to have converged to the right answer
###-------------------------------------------------------------------------------------------------------###

### 2. Initial values : SET 2 (far from real values)
betaSC.2 <- c(-1, -1, 2, -1)

### RUN the Model
  outSC.MM.2 = optim(betaSC.2, SiteConf.MM, data = dataSC, occ = occSC, hessian = T)  
### Comparison real values and estimated values  (initial values also shown)
 round(expit(betaSC.2), 4)			### Initial values for (in order) psi, p11, p10, r11
 round(expit(outSC.MM.2$par), 4)		### Estimated parameters are (in order) psi, p11, p10, r11
 c(psi, p11, p10, s1)				### Real values -- Note that r11 = s1 (= Probability of observing a true positive, given that the site is occupied)

### Deviance: -2LogLik
outSC.MM.2$value*2

###-------------------------------------------------------------------------------------------------------###
### NOTE: psi is biased low, and detection estimates ar off. There was a convergence issue.
### To further convince yourself, try with different real values and initial values.
###-------------------------------------------------------------------------------------------------------###


		###---------------------------------------------------------------###
		###		STEP 3: Implement the Obs Confirmation Model		###
		###---------------------------------------------------------------###

## Likelihood
ObsConf <- function(beta,data,occ){
	psi = expit(beta[1])
	s1 = expit(beta[2])
	s0 = expit(beta[3])
	  p11 <- s1 + s0 - (s1*s0)
	  p10 <- s0
  pi.z1 <- c( (1-s1-s0+s1*s0) , (s0-s1*s0) , (s1-s1*s0), (s1*s0) )
  pi.z0 <- c( (1-s0) , s0 , 0 , 0 )
	y = data[[1]]
	v = 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(v)), 1, prod) +  (1-psi)*apply(t(pi.z0^t(v)), 1, prod)  ) )
return(nLL)
}

## Data
dataOC <- list(Y.sum, V)
occOC <- K

####### COMPARISON #######
### 1. Initial values : SET 1 (close to real values)
betaOC <- c(1, 1, -2)

### RUN the Model
  outOC = optim(betaOC, ObsConf, data = dataOC, occ = occOC, hessian = T)  

### Comparison real values and estimated values 
 round(expit(betaOC), 4)		### Initial values for (in order) psi, s1, s0
 round(expit(outOC$par), 4)		### Estimated parameters are (in order) psi, s1, s0
 c(psi, s1, s0)				### Real values -- Note that r11 = s1 (= Probability of observing a true positive, given that the site is occupied)

### Deviance: -2LogLik
outOC$value*2
###-------------------------------------------------------------------------------------------------------###
### NOTE: Estimated and real values are pretty close -- Seems to have converged to the right answer
### To further convince yourself, try again with different real values and initial values.
###-------------------------------------------------------------------------------------------------------###

### 2. Initial values : SET 2 (far from real values)
betaOC.2 <- c(-1, -1, 2)

### RUN the Model
  outOC.2 = optim(betaOC.2, ObsConf, data = dataOC, occ = occOC, hessian = T)  

### Comparison real values and estimated values 
 round(expit(betaOC.2), 4)		### Initial values for (in order) psi, s1, s0
 round(expit(outOC.2$par), 4)		### Estimated parameters are (in order) psi, s1, s0
 c(psi, s1, s0)				### Real values -- Note that r11 = s1 (= Probability of observing a true positive, given that the site is occupied)

### Deviance: -2LogLik
outOC.2$value*2
###-------------------------------------------------------------------------------------------------------###
### NOTE: Despite the different starting values, the response stays the same (same estimates as above)
### This model appears much more robust to intial values than the Site Confirmation model
###-------------------------------------------------------------------------------------------------------###

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