# SUPPLEMENT #1 - R code for estimating models of rare-event protected species bycatch using WinBUGS.
# Summer L. Martin (1,2), Stephen M. Stohs (2) & Jeffrey E. Moore (2)
# (1) Scripps Insititution of Oceanography and (2) NOAA Southwest Fisheries Science Center
# email: s2martin@ucsd.edu, summerLmartin@gmail.com
# Date: July 30, 2014

# WinBUGS code (run through R using R2WinBUGS package and Coda package) for estimating posterior distributions for model parameters
# for models M1 - M5z in Martin, Stohs & Moore manuscript on Bayesian estimation of rare-event bycatch 
# this code INCLUDES the raw data on observed sets and takes for 1990-2009 from the CA drift gillnet fishery from NOAA's observer program
# (takes shown here are for leatherback sea turtles and humpback whales only)
# therefore, other users can run this code for each model and reproduce our estimation results using this code
# provides parameter estimates and view posterior distribution summaries, DIC values, traceplots, R-hat, etc.
# requires WinBUGS to be installed (R calls WinBUGS to run models) and requires R packages below

library(R2WinBUGS)
library(coda)
library(stats)
#setwd("C:/Users/slm2.NMFS/Documents/Summer's documents/1 NOAA/3 Leatherback Bycatch Project/2 PAPER/PAPER doc and TeX files etc/WORKING DRAFT/SUBMISSION #4 to EA")


# ** Start MODEL 0 - NO EFFORT (number of sets) - did not use this model in paper*************************************************************************************************
# this "no-effort" model is WORSE (in terms of DIC) for leatherbacks (20 points higher than best model)
# and not too different (in terms of DIC) for humpbacks (1 point higher than best model, just like the others -- all pretty close)
sp <- "DC"     # pick a species: "DC" for leatherbacks
#sp <- "MN"    # pick a species: "MN" for humpbacks
# define data
if(sp=="DC") {ntakes <- c(1, 1, 5, 2, 1, 5, 2, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1)}    # leatherback takes 1990-2009
if(sp=="MN") {ntakes <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)}   # humpback takes 1990-2009 (run one species at a time)
# Define Poisson LIKELIHOOD function, pass to WinBUGS. Loops through years.
# purpose is to estimate the take rate theta. Bayesian estimation of parameters in log regression model.
likelihood <- function() {                  
  for (i in 1:20) {                       # loop thru 20 years, 1990-2009
    #lambda[i] <- theta                    # nsets is the number of observed sets in the current data year
    ntakes[i] ~ dpois(lambda)          # the nmbr of takes in the current data year has Poisson distribution
  }  
  lambda ~ dgamma(pr.th.a, pr.th.b)        # dgamma(alpha=takes,beta=sets); gamma keeps theta positive
}
# noninformative PRIOR for theta parameter - gives flat, diffuse, vague prior:
pr.th.a <- .001 #.001 #24       # alpha in gamma dist (takes)      
pr.th.b <- .001 #.001 #20     # beta in gamma dist (sets) ; using 100 here bc it works better in M1z to give comparable pD and DIC to other priors.                
curve(dgamma(x,pr.th.a,pr.th.b), xlim=c(0,2))  # plot prior
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
i.th <-  c(0,1,2)   # INITIAL starting values: for theta parameter - starting values for MCMC chains
data <- list("ntakes", "pr.th.a", "pr.th.b")   # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)                                        
initials <- list(list(lambda = i.th[1]),                 # 3 initial starting points (for the 3 chains) for theta prior, pull from gamma dist with set prior alpha and beta                                                                                                                                                    
                 list(lambda = i.th[2]),
                 list(lambda = i.th[3]))
parameters <- "lambda"
write.model(model=likelihood, con="model.bug")    # writes a text file that is called when "bugs" function calls WinBUGS
# Call WinBUGS using bugs() function 
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots          
#
# ** End MODEL M0 **********************************************************************************************************





# ** Start MODEL M1 *************************************************************************************************
#
sp <- "DC"     # pick a species: "DC" for leatherbacks
#sp <- "MN"    # pick a species: "MN" for humpbacks
# define data
nsets <- c(195, 477, 662, 762, 662, 587, 467, 748, 499, 528, 444, 323, 373, 295, 206, 228, 284, 158, 146, 108)  # observed sets 1990-2009
if(sp=="DC") {ntakes <- c(1, 1, 5, 2, 1, 5, 2, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1)}    # leatherback takes 1990-2009
if(sp=="MN") {ntakes <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)}   # humpback takes 1990-2009 (run one species at a time)
# Define Poisson LIKELIHOOD function, pass to WinBUGS. Loops through years.
# purpose is to estimate the take rate theta. Bayesian estimation of parameters in log regression model.
likelihood <- function() {                  
  for (i in 1:20) {                       # loop thru 20 years, 1990-2009
    lambda[i] <- theta * nsets[i]         # nsets is the number of observed sets in the current data year
    ntakes[i] ~ dpois(lambda[i])          # the nmbr of takes in the current data year has Poisson distribution
  }  
  theta ~ dgamma(pr.th.a, pr.th.b)        # dgamma(alpha=takes,beta=sets); gamma keeps theta positive
}
# noninformative PRIOR for theta parameter - gives flat, diffuse, vague prior:
pr.th.a <- 0       # alpha in gamma dist (takes)      
pr.th.b <- 100     # beta in gamma dist (sets) ; using 100 here bc it works better in M1z to give comparable pD and DIC to other priors.                
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
i.th <-  c(0.002,0.003,0.004)   # INITIAL starting values: for theta parameter - starting values for MCMC chains
data <- list("nsets", "ntakes", "pr.th.a", "pr.th.b")   # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)                                        
initials <- list(list(theta = i.th[1]),                 # 3 initial starting points (for the 3 chains) for theta prior, pull from gamma dist with set prior alpha and beta                                                                                                                                                    
                 list(theta = i.th[2]),
                 list(theta = i.th[3]))
parameters <- "theta"
write.model(model=likelihood, con="model.bug")    # writes a text file that is called when "bugs" function calls WinBUGS
# Call WinBUGS using bugs() function 
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots          
#
# ** End MODEL M1 **********************************************************************************************************







# ** Start MODEL M1z ******************************************************************************************************* 
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"    # pick a species: "MN" for humpbacks
# define data
nsets <- c(195, 477, 662, 762, 662, 587, 467, 748, 499, 528, 444, 323, 373, 295, 206, 228, 284, 158, 146, 108)  # observed sets 1990-2009
if(sp=="DC") {ntakes <- c(1, 1, 5, 2, 1, 5, 2, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1)}    # leatherback takes 1990-2009
if(sp=="MN") {ntakes <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)}   # humpback takes 1990-2009 (run one species at a time)
r <- c(rep(0,11), rep(1,9))     # r=0 = "no closure regulation in effect (1990-2000)"; r=1 = "regulation closure in effect (2001-2009)" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS. Loops through years.
# purpose is to estimate the take rate theta. Bayesian estimation of parameters in log regression model.
# gamma restricts theta to be positive since there can't be a negative take rate
likelihood <- function() {                      
  for (i in 1:20) {                                     # loop thru 20 seasons, 1990-2009
    p[i] <- p1*(1-r[i]) + p2*r[i]                       # r is regulation binary variable; 0 if before regul'n 1990-2000 (p1), 1 if after regul'n 2001-2009 (p2)    
    lambda[i] <- p[i] * theta * nsets[i] + (1-p[i])*0   # nsets is a vector with the number of observed sets per season
    ntakes[i] ~ dpois(lambda[i])                        # observed takes has Poisson distribution 
  }
  theta ~ dgamma(pr.th.a, pr.th.b)  # dgamma(alpha=takes,beta=sets); gamma keeps theta positive
  p1 ~ dbeta(pr.p.a, pr.p.b)        # beta dist constrains ZIP 'p' parameters to be between 0 and 1
  p2 ~ dbeta(pr.p.a, pr.p.b)        # prior for zero-infl poisson parameter (p gives prob that using regular poisson is sufficient over ZIP)  
}
# noninformative PRIOR for theta parameter - gives flat, diffuse, vague prior:
pr.th.a <- 0       # alpha in gamma dist (takes)      
pr.th.b <- 100     # beta in gamma dist (sets) ; using 100 here bc it works better in M1z to give comparable pD and DIC to other priors.                
# noninformative PRIORS for ZIP p parameters 
pr.p.a <- 1        # beta(1,1) is vague, diffuse, flat prior
pr.p.b <- 1        
# ************** INITIAL VALUES *************************
# INITIAL starting values: for theta parameter - starting values for MCMC chains
# turns out starting values are VERY important, particularly in this model, for getting the model to run
# making them realistic based on what the estimates are avoids the errors I kept getting 
# WinBUGS: "cannot sample from interval censored gamma full conditional theta"
# R Console: "Start, end and thin incompatible with data" 
if(sp=="DC") {i.th <- c(0.001, 0.0025, 0.004)}    
if(sp=="MN") {i.th <- c(0.0001, 0.0002, 0.0003)}  
# INITIAL starting values: for ZIP p parameters - starting values for MCMC chains
i.p1 <- c(.5, .55, .6)
i.p2 <- c(.5, .55, .6)
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("r", "nsets", "ntakes", "pr.th.a", "pr.th.b", "pr.p.a", "pr.p.b")         # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)                                        
initials <- list(list(theta = i.th[1], p1=i.p1[1], p2=i.p2[1]),        # 3 initial starting points (for the 3 chains) for theta prior, pull from gamma dist with set prior alpha and beta                                                                                                                                                    
                 list(theta = i.th[2], p1=i.p1[2], p2=i.p2[2]), 
                 list(theta = i.th[3], p1=i.p1[3], p2=i.p2[3])) 
parameters <- c("theta", "p1", "p2")
write.model(model=likelihood, con="model.bug")            # writes a text file that is called when "bugs" function calls WinBUGS 
# leatherbacks - Call WinBUGS using bugs() function 
if(sp=="DC") {bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                               n.chains=3, n.iter=200000, n.burnin=50000, n.sims=20000, debug=FALSE)}    # ... n.thin=5, n.sims=num.sims, took this out and put n.thin=5 instead
# humpbacks
if(sp=="MN") {bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                               n.chains=3, n.iter=200000, n.burnin=50000, n.sims=20000, debug=FALSE)} #n.sims=10000,    # ... n.thin=5, n.sims=num.sims, took this out and put n.thin=5 instead
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots          
#
# ** End MODEL M1z *****************************************************************************************






# ** Start MODEL M1r ****************************************************************************************
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"    # pick a species: "MN" for humpbacks
# define data
nsets <- c(195, 477, 662, 762, 662, 587, 467, 748, 499, 528, 444, 323, 373, 295, 206, 228, 284, 158, 146, 108)  # observed sets 1990-2009
if(sp=="DC") {ntakes <- c(1, 1, 5, 2, 1, 5, 2, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1)}    # leatherback takes 1990-2009
if(sp=="MN") {ntakes <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)}   # humpback takes 1990-2009 (run one species at a time)
r <- c(rep(0,11), rep(1,9))     # r=0 = "no closure regulation in effect (1990-2000)"; r=1 = "regulation closure in effect (2001-2009)" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                  
  for (i in 1:20) {                               # loop thru 20 seasons, 1990-2009
    log(theta[i]) <- beta0 + beta1*r[i]           # r is regulation binary variable; 0 if before regul'n 1990-2000 (p1), 1 if after regul'n 2001-2009 (p2)       
    lambda[i] <- theta[i] * nsets[i]              # number of observed sets in the current data season
    ntakes[i] ~ dpois(lambda[i])                  # the nmbr of takes in the current data season is poisson distributed  
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)         # beta priors: dnorm(mean,tau) where precision tau=1/(sd^2)=1/variance
  beta1 ~ dnorm(pr.b1.a, pr.b1.b)         # beta priors: dnorm(mean,tau) where precision tau=1/(sd^2)=1/variance
}
# noninformative PRIORS for beta0-beta1 parameters; NORMAL distribution (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b1.a <- 0
pr.b1.b <- .000001
# INITIAL starting values: for beta0, beta1, beta2, beta3 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0)  
i.b1 <-  c(-1,0,1)   
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("r", "nsets", "ntakes", "pr.b0.a", "pr.b0.b", "pr.b1.a", "pr.b1.b")         # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)                                        
initials <- list(list(beta0=i.b0[1], beta1=i.b1[1]),  # 3 lists of initials (for the 3 chains) initial starting point for beta priors
                 list(beta0=i.b0[2], beta1=i.b1[2]),
                 list(beta0=i.b0[3], beta1=i.b1[3]))                                                                                                                                                 
parameters <- c("beta0","beta1")                                                                
write.model(model=likelihood, con="model.bug")                            # writes a text file that is called when "bugs" function calls WinBUGS
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)  
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots          
#
# ** End MODEL M1r ******************************************************************************************






# ** START: DEFINE DATA FOR ALL REMAINING MODELS (M2 - M5z) ** #***********************************************
#     -- each model will use these data sets and combine the categories according to the model specifications
# create a matrix of observed sets in 4 data categories:
# Ctgy 1: OUTSIDE closure area / NOT DURING closure time
# Ctgy 2: INSIDE closure area / NOT DURING closure time
# Ctgy 3: OUTSIDE closure area / DURING closure time
# Ctgy 4: INSIDE closure area / DURING closure time
obsvd.sets.data <- matrix(data=c(c(79,214,221,294,266,189,163,249,184,178,181,122,165,119,124,112,130,63,53,42),       # OUTSIDE closure area / NOT DURING closure time
                                 c(15, 6, 41, 74, 4, 42, 81, 41, 19, 0, 2, 0, 4, 0, 0, 23, 14, 0, 4, 4),               # INSIDE closure area / NOT DURING closure time
                                 c(47,94,131,44,155,171,85,246,140,302,223,196,204,176,82,93,139,95,89,62),            # OUTSIDE closure area / DURING closure time
                                 c(54, 163, 269, 350, 237, 185, 138, 212, 156, 48, 38, 5, 0, 0, 0, 0, 1, 0, 0, 0)),    # INSIDE closure area / DURING closure time
                          nrow=20, ncol=4, byrow=FALSE)
# leatherback takes that occurred in the different time-area combinations
dc.take.data <- matrix(data=c(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),   # OUTSIDE closure area / NOT DURING closure time
                              c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),   # INSIDE closure area / NOT DURING closure time
                              c(0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1),   # OUTSIDE closure area / DURING closure time 
                              c(1, 1, 4, 2, 1, 4, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)),  # INSIDE closure area / DURING closure time
                       nrow=20, ncol=4, byrow=FALSE)
# humpback takes that occurred in the different time-area combinations
mn.take.data <- matrix(data=c(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),   # OUTSIDE closure area / NOT DURING closure time
                              c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),   # INSIDE closure area / NOT DURING closure time
                              c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0),   # OUTSIDE closure area / DURING closure time 
                              c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)),  # INSIDE closure area / DURING closure time
                       nrow=20, ncol=4, byrow=FALSE)
#
# ** END: DEFINE DATA FOR ALL REMAINING MODELS (M2-M5z)*******************************************************






# ** Start MODEL M2 (make sure to run data sets above M2 code first) *****************************************
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
# aggregate obsvd sets data into 2-column matrix for observed sets OUTSIDE and INSIDE the closure area 
nsets <- matrix(data=c(obsvd.sets.data[ ,1] + obsvd.sets.data[ ,3], 
                       obsvd.sets.data[ ,2] + obsvd.sets.data[ ,4]), nrow=20, ncol=2)      
# aggregate obsvd takes data into 2-column matrix for leatherback takes OUTSIDE and INSIDE the closure area
if(sp=="DC") {ntakes <- matrix(data=c(dc.take.data[ ,1] + dc.take.data[ ,3], 
                                      dc.take.data[ ,2] + dc.take.data[ ,4]), nrow=20, ncol=2)}
# 2-column matrix for humpback takes OUTSIDE and INSIDE the closure area (aggregate columns)
if(sp=="MN") {ntakes <- matrix(data=c(mn.take.data[ ,1] + mn.take.data[ ,3],
                                      mn.take.data[ ,2] + mn.take.data[ ,4]), nrow=20, ncol=2)}
a <- c(0,1)               # BINARY AREA variable: a=0 = "not inside closure area"; a=1 = "inside closure area" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 2 area categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                                      
  for (i in 1:20) {                              # loop thru 20 seasons, 1990-2009
    for (cg in 1:2){                             # cg is category 1:2; go through each category 1-2
      log(theta[i,cg]) <- beta0 + beta1*a[cg]    # a is area variable; 0 if not inside closure area, 1 if inside closure area      
      lambda[i,cg] <- theta[i,cg]*nsets[i,cg]
      ntakes[i,cg] ~ dpois(lambda[i,cg])
    } 
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta1 ~ dnorm(pr.b1.a, pr.b1.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
} 
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b1.a <- 0
pr.b1.b <- .000001
# INITIAL starting values: for beta0, beta1 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0)  
i.b1 <-  c(-1,0,1)   
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("a", "nsets", "ntakes",        # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)
             "pr.b0.a", "pr.b0.b", "pr.b1.a", "pr.b1.b")                                 
initials <- list(list(beta0=i.b0[1], beta1=i.b1[1]),  # 3 lists of initials (for the 3 chains) initial starting point for beta priors
                 list(beta0=i.b0[2], beta1=i.b1[2]),
                 list(beta0=i.b0[3], beta1=i.b1[3]))                                                                                                                                                 
parameters <- c("beta0","beta1")                                                                
write.model(model=likelihood, con="model.bug")      # writes a text file that is called when "bugs" function calls WinBUGS          
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)  
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots          
#
# ** End MODEL M2 ****************************************************************************************






# ** Start MODEL M2z (make sure to run data sets above M2 code first) ************************************
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
# aggregate obsvd sets data into 2-column matrix for observed sets OUTSIDE and INSIDE the closure area 
nsets <- matrix(data=c(obsvd.sets.data[ ,1] + obsvd.sets.data[ ,3], 
                       obsvd.sets.data[ ,2] + obsvd.sets.data[ ,4]), nrow=20, ncol=2)      
# aggregate obsvd takes data into 2-column matrix for leatherback takes OUTSIDE and INSIDE the closure area
if(sp=="DC") {ntakes <- matrix(data=c(dc.take.data[ ,1] + dc.take.data[ ,3], 
                                      dc.take.data[ ,2] + dc.take.data[ ,4]), nrow=20, ncol=2)}
# aggregate obsvd takes data into 2-column matrix for humpback takes OUTSIDE and INSIDE the closure area (aggregate columns)
if(sp=="MN") {ntakes <- matrix(data=c(mn.take.data[ ,1] + mn.take.data[ ,3],
                                      mn.take.data[ ,2] + mn.take.data[ ,4]), nrow=20, ncol=2)}
a <- c(0,1)               # BINARY AREA variable: a=0 = "not inside closure area"; a=1 = "inside closure area" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 2 area categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                                      
  for (i in 1:20) {                             # loop thru 20 seasons, 1990-2009     
    for (cg in 1:2){                            # cg is category 1:2; go through each category 1-2
      p[i,cg] <- p1*(1-a[cg]) + p2*a[cg]        # Zero-Inflated parameters (1-p)
      log(theta[i,cg]) <- beta0 + beta1*a[cg]   # a is area variable; 0 if not inside closure area, 1 if inside closure area      
      lambda[i,cg] <- p[i,cg]*theta[i,cg]*nsets[i,cg] + (1-p[i,cg])*0
      ntakes[i,cg] ~ dpois(lambda[i,cg])
    } 
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta1 ~ dnorm(pr.b1.a, pr.b1.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  p1 ~ dbeta(pr.p.a, pr.p.b)          # beta dist constrains ZIP 'p' parameters to be between 0 and 1
  p2 ~ dbeta(pr.p.a, pr.p.b)          # prior for zero-infl poisson parameter (p gives prob that using regular poisson is sufficient over ZIP)  
}   
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b1.a <- 0
pr.b1.b <- .000001
# noninformative PRIORS for ZIP p parameters 
pr.p.a <- 1        # beta(1,1) is vague, diffuse, flat prior
pr.p.b <- 1        
# INITIAL starting values: for beta0, beta1 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0)  
i.b1 <-  c(-1,0,1)   
# INITIAL starting values: for ZIP p parameters - starting values for MCMC chains
i.p1 <- c(0.45, 0.5, 0.55) 
i.p2 <- c(0.45, 0.5, 0.55) 
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("a", "nsets", "ntakes",        # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)
             "pr.b0.a", "pr.b0.b", "pr.b1.a", "pr.b1.b", "pr.p.a", "pr.p.b")                                 
initials <- list(list(beta0=i.b0[1], beta1=i.b1[1], p1=i.p1[1], p2=i.p2[1]),  # 3 lists of initials (for the 3 chains) initial starting point for beta priors
                 list(beta0=i.b0[2], beta1=i.b1[2], p1=i.p1[2], p2=i.p2[2]), 
                 list(beta0=i.b0[3], beta1=i.b1[3], p1=i.p1[3], p2=i.p2[3]))                                                                                                                                                
parameters <- c("beta0","beta1", "p1", "p2")                                                                
write.model(model=likelihood, con="model.bug")      # writes a text file that is called when "bugs" function calls WinBUGS          
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)     
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots          
#
# ** End MODEL M2 *************************************************************************************************






# ** Start MODEL M3 (make sure to run data sets above M2 code first) ***********************************************
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
# aggregate obsvd sets data into 2-column matrix for observed sets NOT DURING and DURING the closure time (Aug 15-Nov15) 
nsets <- matrix(data=c(obsvd.sets.data[ ,1] + obsvd.sets.data[ ,2], 
                       obsvd.sets.data[ ,3] + obsvd.sets.data[ ,4]), nrow=20, ncol=2)      
# aggregate obsvd takes data into 2-column matrix for leatherback takes NOT DURING and DURING the closure time (Aug 15-Nov15) 
if(sp=="DC") {ntakes <- matrix(data=c(dc.take.data[ ,1] + dc.take.data[ ,2], 
                                      dc.take.data[ ,3] + dc.take.data[ ,4]), nrow=20, ncol=2)}
# aggregate obvsd takes data into 2-column matrix for humpback takes NOT DURING and DURING the closure time (Aug 15-Nov15) 
if(sp=="MN") {ntakes <- matrix(data=c(mn.take.data[ ,1] + mn.take.data[ ,2],
                                      mn.take.data[ ,3] + mn.take.data[ ,4]), nrow=20, ncol=2)}
t <- c(0,1)               # BINARY TIME variable: t=0 = "not during closure time"; t=1 = "during closure time" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 2 time categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                                      # LOOP THROUGH CATEGORIES
  for (i in 1:20) {                                             # loop thru 20 seasons, 1990-2009
    for (cg in 1:2){                                      # cg is category 1:2; go through each category 1-2
      log(theta[i,cg]) <- beta0 + beta2*t[cg]             # 't' is time binary variable defined above      
      lambda[i,cg] <- theta[i,cg]*nsets[i,cg]
      ntakes[i,cg] ~ dpois(lambda[i,cg])
    }
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta2 ~ dnorm(pr.b2.a, pr.b2.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
}
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b2.a <- 0
pr.b2.b <- .000001
# INITIAL starting values: for beta0, beta2 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0)  
i.b2 <-  c(-1,0,1)   
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("t", "nsets", "ntakes",                             # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)
             "pr.b0.a", "pr.b0.b", "pr.b2.a", "pr.b2.b")                                 
initials <- list(list(beta0=i.b0[1], beta2=i.b2[1]),  # 3 lists of initials (for the 3 chains) initial starting point for beta priors
                 list(beta0=i.b0[2], beta2=i.b2[2]),
                 list(beta0=i.b0[3], beta2=i.b2[3]))                                                                                                                                                 
parameters <- c("beta0","beta2")                                                                
write.model(model=likelihood, con="model.bug")        # writes a text file that is called when 'bugs' function calls WinBUGS  
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)  
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots          
# ** End MODEL M3 ***************************************************************************************************






# ** Start MODEL M3z (make sure to run data sets above M2 code first) ***********************************************
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
# aggregate obsvd sets data into 2-column matrix for observed sets NOT DURING and DURING the closure time (Aug 15-Nov15) 
nsets <- matrix(data=c(obsvd.sets.data[ ,1] + obsvd.sets.data[ ,2], 
                       obsvd.sets.data[ ,3] + obsvd.sets.data[ ,4]), nrow=20, ncol=2)      
# aggregate obsvd takes data into 2-column matrix for leatherback takes NOT DURING and DURING the closure time (Aug 15-Nov15) 
if(sp=="DC") {ntakes <- matrix(data=c(dc.take.data[ ,1] + dc.take.data[ ,2], 
                                      dc.take.data[ ,3] + dc.take.data[ ,4]), nrow=20, ncol=2)}
# aggregate obvsd takes data into 2-column matrix for humpback takes NOT DURING and DURING the closure time (Aug 15-Nov15) 
if(sp=="MN") {ntakes <- matrix(data=c(mn.take.data[ ,1] + mn.take.data[ ,2],
                                      mn.take.data[ ,3] + mn.take.data[ ,4]), nrow=20, ncol=2)}
t <- c(0,1)               # BINARY TIME variable: t=0 = "not during closure time"; t=1 = "during closure time" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 2 time categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                                      # LOOP THROUGH CATEGORIES
  for (i in 1:20) {                                             # loop thru 20 seasons, 1990-2009
    for (cg in 1:2){                                    # cg is category 1:2; go through each category 1-2
      p[i,cg] <- p1*(1-t[cg]) + p2*t[cg]                # Zero-Inflated parameters (1-p)
      log(theta[i,cg]) <- beta0 + beta2*t[cg]           # 't' is time binary variable defined above         
      lambda[i,cg] <- p[i,cg]*theta[i,cg]*nsets[i,cg] + (1-p[i,cg])*0
      ntakes[i,cg] ~ dpois(lambda[i,cg])
    }
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta2 ~ dnorm(pr.b2.a, pr.b2.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  p1 ~ dbeta(pr.p.a, pr.p.b)        # beta dist constrains ZIP 'p' parameters to be between 0 and 1
  p2 ~ dbeta(pr.p.a, pr.p.b)        # prior for zero-infl poisson parameter (p gives prob that using regular poisson is sufficient over ZIP)  
}
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b2.a <- 0
pr.b2.b <- .000001
# noninformative PRIORS for ZIP p parameters 
pr.p.a <- 1        # beta(1,1) is vague, diffuse, flat prior
pr.p.b <- 1        
# INITIAL starting values: for beta0, beta1, beta2, beta3 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0)  
i.b2 <-  c(-1,0,1)   
# INITIAL starting values: for ZIP p parameters - starting values for MCMC chains
i.p1 <- c(0.45, 0.5, 0.55)
i.p2 <- c(0.45, 0.5, 0.55)
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("t", "nsets", "ntakes", 
             "pr.b0.a", "pr.b0.b", "pr.b2.a", "pr.b2.b", "pr.p.a", "pr.p.b")             # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)                    
initials <- list(list(beta0=i.b0[1], beta2=i.b2[1], p1=i.p1[1], p2=i.p2[1]),  # 3 lists of initials (for the 3 chains) initial starting point for beta priors
                 list(beta0=i.b0[2], beta2=i.b2[2], p1=i.p1[2], p2=i.p2[2]), 
                 list(beta0=i.b0[3], beta2=i.b2[3], p1=i.p1[3], p2=i.p2[3]))                                                                                                                                                 
parameters <- c("beta0","beta2", "p1", "p2")                                                                
write.model(model=likelihood, con="model.bug")        # writes a text file that is called when 'bugs' function calls WinBUGS  
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)  
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots   
#
# ** End MODEL M3z *************************************************************************************************







# ** Start MODEL M4 (make sure to run data sets above M2 code first) ***********************************************
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
# aggregate obsvd sets data into 2-column matrix for observed sets INSIDE/DURING closure and those NOT (INSIDE/DURING closure) 
nsets <- matrix(data=c(obsvd.sets.data[ ,1] + obsvd.sets.data[ ,2] + obsvd.sets.data[ ,3],
                        obsvd.sets.data[ ,4]), nrow=20, ncol=2)      
# aggregate obsvd takes data into 2-column matrix for leatherback takes 
if(sp=="DC") {ntakes <- matrix(data=c(dc.take.data[ ,1] + dc.take.data[ ,2] + dc.take.data[ ,3], 
                                      dc.take.data[ ,4]), nrow=20, ncol=2)}
# aggregate obvsd takes data into 2-column matrix for humpback takes 
if(sp=="MN") {ntakes <- matrix(data=c(mn.take.data[ ,1] + mn.take.data[ ,2] +  mn.take.data[ ,3],
                                      mn.take.data[ ,4]), nrow=20, ncol=2)}
at <- c(0,1)               # BINARY TIME-AREA variable: at=0 = "*NOT* INSIDE/DURING closure"; at=1 = "INSIDE/DURING closure" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 2 time-area categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                                      # LOOP THROUGH CATEGORIES
  for (i in 1:20) {                                             # loop thru 20 seasons, 1990-2009
    for (cg in 1:2){                                            # 'cg' is category 1:2 ; go through each category 1-2
      log(theta[i,cg]) <- beta0 + beta3*at[cg]                  # 'at' is time-area interaction variable; 0 if not inside&during, 1 if inside&during closure       
      lambda[i,cg] <- theta[i,cg]*nsets[i,cg]
      ntakes[i,cg] ~ dpois(lambda[i,cg])                                                                               
    }
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta3 ~ dnorm(pr.b3.a, pr.b3.b)             
}
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b3.a <- 0
pr.b3.b <- .000001
# INITIAL starting values: for beta0, beta1, beta2, beta3 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0)  
i.b3 <-  c(-1,0,1)   
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("at", "nsets", "ntakes",        # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)
             "pr.b0.a", "pr.b0.b", "pr.b3.a", "pr.b3.b")                                 
initials <- list(list(beta0=i.b0[1], beta3=i.b3[1]),  # 3 lists of initials (for the 3 chains) initial starting point for beta priors
                 list(beta0=i.b0[2], beta3=i.b3[2]),
                 list(beta0=i.b0[3], beta3=i.b3[3]))                                                                                                                                                 
parameters <- c("beta0","beta3")                                         # parameters for WinBUGS to monitor and save                        
write.model(model=likelihood, con="model.bug")                            # writes a text file that is called when "bugs" function calls WinBUGS
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)  
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots   
#
# ** End MODEL M4 *************************************************************************************************






# ** Start MODEL M4z (make sure to run data sets above M2 code first) ***********************************************
#
#sp <- "DC"     # pick a species: "DC" for leatherbacks
sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
# aggregate obsvd sets data into 2-column matrix for observed sets INSIDE/DURING closure and those NOT (INSIDE/DURING closure) 
nsets <- matrix(data=c(obsvd.sets.data[ ,1] + obsvd.sets.data[ ,2] + obsvd.sets.data[ ,3],
                       obsvd.sets.data[ ,4]), nrow=20, ncol=2)      
# aggregate obsvd takes data into 2-column matrix for leatherback takes 
if(sp=="DC") {ntakes <- matrix(data=c(dc.take.data[ ,1] + dc.take.data[ ,2] + dc.take.data[ ,3], 
                                      dc.take.data[ ,4]), nrow=20, ncol=2)}
# aggregate obvsd takes data into 2-column matrix for humpback takes 
if(sp=="MN") {ntakes <- matrix(data=c(mn.take.data[ ,1] + mn.take.data[ ,2] +  mn.take.data[ ,3],
                                      mn.take.data[ ,4]), nrow=20, ncol=2)}
at <- c(0,1)               # BINARY TIME-AREA variable: at=0 = "*NOT* INSIDE/DURING closure"; at=1 = "INSIDE/DURING closure" 
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 2 time-area categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                                      # LOOP THROUGH CATEGORIES
  for (i in 1:20) {                                             # loop thru 20 seasons, 1990-2009
    for (cg in 1:2){                                            # 'cg' is category 1:2 ; go through each category 1-2
      p[i,cg] <- p1*(1-at[cg]) + p2*at[cg]                      # let p vary by category, same as theta
      log(theta[i,cg]) <- beta0 + beta3*at[cg]                  # 'at' is time-area interaction variable; 0 if not inside&during, 1 if inside&during closure       
      lambda[i,cg] <- p[i,cg]*theta[i,cg]*nsets[i,cg] + (1-p[i,cg])*0
      ntakes[i,cg] ~ dpois(lambda[i,cg])                                                                               
    }
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta3 ~ dnorm(pr.b3.a, pr.b3.b)             
  p1 ~ dbeta(pr.p.a, pr.p.b)        # beta dist constrains ZIP 'p' parameters to be between 0 and 1
  p2 ~ dbeta(pr.p.a, pr.p.b)        # prior for zero-infl poisson parameter (p gives prob that using regular poisson is sufficient over ZIP)  
}
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b3.a <- 0
pr.b3.b <- .000001
# noninformative PRIORS for ZIP p parameters 
pr.p.a <- 1        # beta(1,1) is vague, diffuse, flat prior
pr.p.b <- 1        
# INITIAL starting values: for beta0, beta1, beta2, beta3 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0)  
i.b3 <-  c(-1,0,1)   
# INITIAL starting values: for ZIP p parameters - starting values for MCMC chains
i.p1 <- c(0.4, 0.5, 0.6)
i.p2 <- c(0.4, 0.5, 0.6)
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS)     
data <- list("at", "nsets", "ntakes",        # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS)
             "pr.b0.a", "pr.b0.b", "pr.b3.a", "pr.b3.b", "pr.p.a", "pr.p.b")                                  
initials <- list(list(beta0=i.b0[1], beta3=i.b3[1], p1=i.p1[1], p2=i.p2[1]),   # 3 lists of initials (for the 3 chains) initial starting point for beta priors
                 list(beta0=i.b0[2], beta3=i.b3[2], p1=i.p1[2], p2=i.p2[2]), 
                 list(beta0=i.b0[3], beta3=i.b3[3], p1=i.p1[3], p2=i.p2[3]))                                                                                                                                                  
parameters <- c("beta0","beta3", "p1", "p2")                              # parameters for WinBUGS to monitor and save                        
write.model(model=likelihood, con="model.bug")                            # writes a text file that is called when "bugs" function calls WinBUGS
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)  
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots   
#
# ** End MODEL M4z *************************************************************************************************






# ** Start MODEL M5 (make sure to run data sets above M2 code first) ***********************************************
#
#sp <- "DC"       # pick a species: "DC" for leatherbacks
sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
nsets <- obsvd.sets.data                # all four time-area categories are present in this model
if(sp=="DC") {ntakes <- dc.take.data}   # leatherback takes - all 4 categories
if(sp=="MN") {ntakes <- mn.take.data}   # humpback takes - all 4 categories
# Binary variables (a=AREA, t=TIME, at=AREA*TIME)
a <- c(0,1,0,1)               # Area Closure (inside?); categ 1-4: AreaNo, AreaYes, AreaNo, AreaYes
t <- c(0,0,1,1)               # Time Closure (during?); categ 1-4: TimeNo, TimeNo, TimeYes, TimeYes
at <- c(0,0,0,1)              # Intxn Area-Time Closure categ; 1-4: InteractionNo, InteractionNo, InteractionNo, InteractionYes
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 4 time-area categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                                          # LOOP THROUGHCATEGORIES
  for (i in 1:20) {                                                 # loop thru 20 seasons, 1990-2009
    for (cg in 1:4){                                                # cg is category 1:4 ;  go through each category 1-4
      log(theta[i,cg]) <- beta0 + beta1*a[cg] + beta2*t[cg] + beta3*at[cg]   # a, t, and at are vectors with just four 0s and 1s, so need to index w/ cg 1:4       
      lambda[i,cg] <- theta[i,cg]*nsets[i,cg]
      ntakes[i,cg] ~ dpois(lambda[i,cg])
    }
  }
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)   # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta1 ~ dnorm(pr.b1.a, pr.b1.b)
  beta2 ~ dnorm(pr.b2.a, pr.b2.b)
  beta3 ~ dnorm(pr.b3.a, pr.b3.b)
}
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b1.a <- 0
pr.b1.b <- .000001
pr.b2.a <- 0
pr.b2.b <- .000001
pr.b3.a <- 0
pr.b3.b <- .000001
# INITIAL starting values: for beta0, beta1, beta2, beta3 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0) 
i.b1 <-  c(-1,0,1)   
i.b2 <-  c(-1,0,1)   
i.b3 <-  c(-1,0,1)   
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS) 
data <- list("a", "t", "at", "nsets", "ntakes",        # specifies which data objects will be used when calling the "bugs" function (to call WinBUGS) 
             "pr.b0.a", "pr.b0.b", "pr.b1.a", "pr.b1.b","pr.b2.a", "pr.b2.b", "pr.b3.a", "pr.b3.b")                                 
initials <- list(list(beta0=i.b0[1], beta1=i.b1[1], beta2=i.b2[1], beta3=i.b3[1]),  # 3 lists of initials (for the 3 chains) initial starting point for beta priors, pull from normal dist with set mean and sd
                 list(beta0=i.b0[2], beta1=i.b1[2], beta2=i.b2[2], beta3=i.b3[2]),
                 list(beta0=i.b0[3], beta1=i.b1[3], beta2=i.b2[3], beta3=i.b3[3]))                                                                                                                                                 
parameters <- c("beta0","beta1", "beta2", "beta3")                                                                
write.model(model=likelihood, con="model.bug")                            # writes a text file that is called when "bugs" function calls WinBUGS    
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000, n.sims=20000, debug=FALSE)
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots   
#
# ** End MODEL M5 *************************************************************************************************





# ** Start MODEL M5z (make sure to run data sets above M2 code first) ***********************************************
#
sp <- "DC"       # pick a species: "DC" for leatherbacks
#sp <- "MN"      # pick a species: "MN" for humpbacks
# define data
nsets <- obsvd.sets.data                # all four time-area categories are present in this model
if(sp=="DC") {ntakes <- dc.take.data}   # leatherback takes - all 4 categories
if(sp=="MN") {ntakes <- mn.take.data}   # humpback takes - all 4 categories
# Binary variables (a=AREA, t=TIME, at=AREA*TIME)
a <-  c(0,1,0,1)               # Area Closure (inside?); categ 1-4: AreaNo, AreaYes, AreaNo, AreaYes
t <-  c(0,0,1,1)               # Time Closure (during?); categ 1-4: TimeNo, TimeNo, TimeYes, TimeYes
at <- c(0,0,0,1)              # Intxn Area-Time Closure categ; 1-4: InteractionNo, InteractionNo, InteractionNo, InteractionYes
NI <- c(1,0,0,0)              # outside both time and area closure
IA <- c(0,1,0,0)              # inside closure AREA only
IT <- c(0,0,1,0)              # inside closure TIME only
# Define Poisson LIKELIHOOD function, pass to WinBUGS, loop through years and 4 time-area categories of data
# priors for beta parameters have NORMAL distributions here with large standard deviations so they are noninformative, flat, diffuse, vague
likelihood <- function() {                # LOOP THROUGH CATEGORIES
  for (i in 1:20) {                       # loop thru 20 seasons, 1990-2009
    for (cg in 1:4){                                                         # cg is category 1:4 ; go through each category 1-4
      p[i,cg] <- p1*NI[cg] + p2*IA[cg] + p3*IT[cg] + p4*at[cg]               # one 'p' for each category
      log(theta[i,cg]) <- beta0 + beta1*a[cg] + beta2*t[cg] + beta3*at[cg]   # a, t, and 'at' are vectors with just four 0s and 1s, so need to index w/ cg 1:4     
      lambda[i,cg] <- p[i,cg]*theta[i,cg]*nsets[i,cg] + (1-p[i,cg])*0        
      ntakes[i,cg] ~ dpois(lambda[i,cg])                        
    } 
  }         
  beta0 ~ dnorm(pr.b0.a, pr.b0.b)     # beta priors: dnorm(mean,tau) where tau=1/(sd^2)
  beta1 ~ dnorm(pr.b1.a, pr.b1.b)
  beta2 ~ dnorm(pr.b2.a, pr.b2.b)
  beta3 ~ dnorm(pr.b3.a, pr.b3.b)
  p1 ~ dbeta(pr.p.a, pr.p.b)        # beta dist constrains ZIP 'p' parameters to be between 0 and 1
  p2 ~ dbeta(pr.p.a, pr.p.b)
  p3 ~ dbeta(pr.p.a, pr.p.b)                                 
  p4 ~ dbeta(pr.p.a, pr.p.b) 
}
# noninformative PRIORS for beta0-beta3 parameters; NORMAL distribution parameters (a=mean, b=tao precision)
# for dnorm(mean, tau): tau=.0001 is sd=100; tau=.000001 is sd=1000
pr.b0.a <- 0        
pr.b0.b <- .000001   # this precision is sd=1000    
pr.b1.a <- 0
pr.b1.b <- .000001
pr.b2.a <- 0
pr.b2.b <- .000001
pr.b3.a <- 0
pr.b3.b <- .000001
# noninformative PRIORS for ZIP p parameters 
pr.p.a <- 1        # beta(1,1) is vague, diffuse, flat prior
pr.p.b <- 1        
# INITIAL starting values: for beta0, beta1, beta2, beta3 parameters - starting values for MCMC chains
i.b0 <-  c(-15,-7,0) 
i.b1 <-  c(-1,0,1)   
i.b2 <-  c(-1,0,1)   
i.b3 <-  c(-1,0,1)   
# INITIAL starting values: for ZIP p parameters - starting values for MCMC chains
i.p1 <- c(0.4,0.5,0.6)
i.p2 <- c(0.4,0.5,0.6)
i.p3 <- c(0.4,0.5,0.6)
i.p4 <- c(0.4,0.5,0.6)
# Use MCMC to estimate model parameters using data and specified priors 
# R to WINBUGS: Settings in R to save parameters and set initial starting values (sends these to WinBUGS) 
data <- list("NI", "IA", "IT", "a", "t", "at", "nsets", "ntakes",
             "pr.b0.a", "pr.b0.b", "pr.b1.a", "pr.b1.b","pr.b2.a", "pr.b2.b", "pr.b3.a", "pr.b3.b", "pr.p.a", "pr.p.b")                                 
initials <- list(list(beta0=i.b0[1], beta1=i.b1[1], beta2=i.b2[1], beta3=i.b3[1], p1=i.p1[1], p2=i.p2[1], p3=i.p3[1], p4=i.p4[1]),   # 3 lists of initials (for the 3 chains) initial starting point for beta priors, pull from normal dist with set mean and sd
                 list(beta0=i.b0[2], beta1=i.b1[2], beta2=i.b2[2], beta3=i.b3[2], p1=i.p1[2], p2=i.p2[2], p3=i.p3[2], p4=i.p4[2]), 
                 list(beta0=i.b0[3], beta1=i.b1[3], beta2=i.b2[3], beta3=i.b3[3], p1=i.p1[3], p2=i.p2[3], p3=i.p3[3], p4=i.p4[3]))                                                                                                                                                
parameters <- c("beta0","beta1", "beta2", "beta3", "p1", "p2", "p3", "p4")                                                                
write.model(model=likelihood, con="model.bug")                            # writes a text file that is called when "bugs" function calls WinBUGS          
bugs.sim <- bugs(data=data, inits=initials, parameters.to.save=parameters, model.file="model.bug",
                 n.chains=3, n.iter=200000, n.burnin=50000,  n.sims=20000, debug=FALSE)  
print(bugs.sim, digits=5)               # prints summary of estimation (parameter estimates, pD, DIC)
mcmc.tracer <- as.mcmc.list(bugs.sim)   # list provides vector outputs from each chain (parameter posteriors + deviance)    
summary(mcmc.tracer)
traceplot(mcmc.tracer)        # plots parameters vs. iterations (convergence test) - if stable, hovering around mode with no trend, then indicates convergence
plot(mcmc.tracer)             # plots density of parameters (theta) and deviance, along with trace plots   
#
# ** End MODEL M5z *************************************************************************************************