#JAGS code to run the Bayesain hierarchical species distribution model (BHSDM) #descirbed in Sigourney et al. #Data are fit to line transect data on fin whales from four distinct #surveys conducted from the summer of 2010 to the fall of 2013. #Shipboard and aerial surveys were conducted by the Northeast Fisheries Science (NEFSC) and # the Southeast Fisheries Science Center (SEFSC). # For surveys 1=NEFSC and 2=SEFSC and for platforms 1=Shipboard and 2=aerial (e.g. 1_2 means NEFSC_Aerial) #Libraries library(mgcv) library(rjags) ################################### #Import data files ################################### Dat_Sightings<-read.csv("Dat_Sightings_Fin_Whales.csv") ncols_Sights<-dim(Dat_Sightings)[2] Dat_Sightings<-Dat_Sightings[,2:ncols_Sights] Final_Cov_Mat<-read.csv("Final_Cov_Mat_Fin_Whales.csv") ncols_Covs<-dim(Final_Cov_Mat)[2] Final_Cov_Mat<- Final_Cov_Mat[,2:ncols_Covs] #################################### #Prepare Data inputs into model #################################### #NE Ship Dat_Sightings_1_1<-subset(Dat_Sightings, Survey==1 & Platform==1) N_Sights_1_1<-dim(Dat_Sightings_1_1)[1] ones.dist_1_1<-array(NA,N_Sights_1_1) for (i in 1:N_Sights_1_1){ ones.dist_1_1[i]<-1 } #SE Ship Dat_Sightings_2_1<-subset(Dat_Sightings, Survey==2 & Platform==1) N_Sights_2_1<-dim(Dat_Sightings_2_1)[1] ones.dist_2_1<-array(NA,N_Sights_2_1) for (i in 1:N_Sights_2_1){ ones.dist_2_1[i]<-1 } #NE Aerial Dat_Sightings_1_2<-subset(Dat_Sightings, Survey==1 & Platform==2) N_Sights_1_2<-dim(Dat_Sightings_1_2)[1] ones.dist_1_2<-array(NA,N_Sights_1_2) for (i in 1:N_Sights_1_2){ ones.dist_1_2[i]<-1 } #SE Aerial Dat_Sightings_2_2<-subset(Dat_Sightings, Platform==2 & Survey==2) N_Sights_2_2<-dim(Dat_Sightings_2_2)[1] ones.dist_2_2<-array(NA,N_Sights_2_2) for (i in 1:N_Sights_2_2){ ones.dist_2_2[i]<-1 } #Truncation distances (meters) W_1_1=5000 #NE Shipboard W_2_1=8840 #SE Shipboard W_1_2=1000 #NE Aerial W_2_2=600 #SE Aerial #Subset out fin whale sightings with observations of group size Dat_Sightings_Fin<-subset(Dat_Sightings, Species==1) N.size<-dim(Dat_Sightings_Fin)[1] gs<-Dat_Sightings_Fin$Group.Size #Observed group sizes of fin whales #Total observed grid cells N_t<-dim(Final_Cov_Mat)[1] #Results from g(0) to develop informative priors #NE Surveys g0_est_NE=0.503 g0_se_NE=0.08 g0_CV_NE=g0_se_NE/g0_est_NE # g0_CV_NE=0.168 c_g0_NE=(g0_est_NE*(1-g0_est_NE)/(g0_est_NE*g0_CV_NE)^2)-1 a_g0_NE=c_g0_NE*g0_est_NE b_g0_NE=c_g0_NE*(1-g0_est_NE) #SE Surveys g0_est_SE=0.898 g0_se_SE=0.10 g0_CV_SE=g0_se_SE/g0_est_SE # g0_CV_SE=0.091 c_g0_SE=(g0_est_SE*(1-g0_est_SE)/(g0_est_SE*g0_CV_SE)^2)-1 a_g0_SE=c_g0_SE*g0_est_SE b_g0_SE=c_g0_SE*(1-g0_est_SE) a_g0<-c(a_g0_NE, a_g0_SE) b_g0<-c(b_g0_NE, b_g0_SE) #Information on surface availability to develop informative priors Av_est=0.374 Av_CV=0.336 c_Av=((Av_est*(1-Av_est))/((Av_est*Av_CV)^2))-1 a_Av=c_Av*Av_est b_Av=c_Av*(1-Av_est) zeros<-matrix(0,nrow=N_t, ncol=1) Platform<-Final_Cov_Mat$Platform Survey<-Final_Cov_Mat$Survey Indicator_1_1<-as.numeric(Survey==1)*as.numeric(Platform==1) #NE Shipboard Indicator_2_1<-as.numeric(Survey==2)*as.numeric(Platform==1) #SE Shipboard Indicator_1_2<-as.numeric(Survey==1)*as.numeric(Platform==2) #NE Aerial Indicator_2_2<-as.numeric(Survey==2)*as.numeric(Platform==2) #SE Aerial #To integrate hazard rate model using Heun's method steps<-50 #For NE Shipboard y.step_all_1_1<-seq(0.01,W_1_1/1000, length=steps) #Convert to KM y.step_1_1<-y.step_all_1_1[1:(steps-1)] y.step_plus1_1_1<-y.step_all_1_1[2:steps] step.size_1_1<-(W_1_1/1000)/(steps-1) #Convert to KM #For NE Aerial y.step_all_1_2<-seq(0.01,W_1_2/1000, length=steps) #Convert to KM y.step_1_2<-y.step_all_1_2[1:(steps-1)] y.step_plus1_1_2<-y.step_all_1_2[2:steps] step.size_1_2<-(W_1_2/1000)/(steps-1) #Convert to KM #Inputs for GAM model using jagam function from the mgcv package jd <- jagam(NSights~s(CHL, bs="ts", k=5)+s(BTEMP,bs="ts", k=5)+s(DIST125,bs="ts", k=5)+s(LAT,bs="ts", k=5), data= Final_Cov_Mat, family = poisson(link=log), file = "Fin_jagam") #Inputs for smooth functions from jagam X<-jd$jags.data$X S1<-jd$jags.data$S1 S2<-jd$jags.data$S2 S3<-jd$jags.data$S3 S4<-jd$jags.data$S4 # S5<-jd$jags.data$S5 zero<-jd$jags.data$zero #GAM Inits b<-jd$jags.ini$b lambda<-jd$jags.ini$lambda #Format Data for CPG distribution Sightings=Final_Cov_Mat[,25] NonZeroObs<-which(Sightings>0) Nnonzero<-length(NonZeroObs) ZeroObs<-which(Sightings==0) Nzero<-length(ZeroObs) sink("Fin_jagam_CPG.bugs") cat(" model{ #Parameters for integrating half-normal detection function (adapted from Moore and Barlow 2011) pi <- 3.141593 # approx constant used in the error function (erf) for estimating esw, # described above a <- 0.147 # priors for fixed effect parms for half-normal detection parm sigma ##### PRIORS ###### # MRDS COMPONENT OF MODEL FOR SURFACE DETECTABILYTY ##### #Priors on MR Parameters (just shipboard surveys) #NE Ship b1.0_1 ~ dnorm(0,0.0001) b.dist_1~ dnorm(0,0.0001) #SE Ship b1.0_2 ~ dnorm(0,0.0001) b2.0_2 ~ dnorm(0,0.0001) b.dist1_2~ dnorm(0,0.0001) b.dist2_2~ dnorm(0,0.0001) b.beauf_2~ dnorm(0,0.0001) #Priors on DS Parameters #Scale Parameters #Hazard Rate b.df.0_1_1 ~ dunif(0,10) #NE Ship b.df.0_1_2 ~ dunif(0,10)#NE Aerial #Half-Normal b.df.0_2_1 ~ dnorm(0,0.0001) #SE Ship b.df.0_2_2 ~ dnorm(0,0.0001) #SE Aerial #Shape parameters for hazard rate function b_1_1 ~ dunif(0,10) #NE Shipboard b_1_2 ~ dunif(0,10) #NE Aerial #Beaufort effect b.df.beaufort_2_1 ~ dnorm(0,0.0001) #SE Ship b.df.beaufort_2_2 ~ dnorm(0,0.0001) #SE Aerial #Prior for power parameter of CPG likelihood nu <- (2 - zeta)/(zeta - 1) zeta ~ dunif(1.01,2) #Prior for average group size b0.size~dnorm(0, 0.001) log(lambda.size)<-b0.size #Informed Priors for g(0) for aerial surveys for (c in 1:2){ g.knot[c]~dbeta(a_g0[c], b_g0[c]) } #Informed Prior for surafce avaialbility Av~dbeta(a_Av, b_Av) #Group size model for (g in 1:N.size){ gs[g]~dpois(lambda.size) } #NE Shipboard for(i in 1:N_Sights_1_1){ ############################################################################################## # HAZARD RATE KEY ############################################################################################ mu.df_1_1[i] <- b.df.0_1_1 esw_pred_1_1[i,1]<-0 #Numerically integrate hazard rate function for (c in 1:(steps-1)){ dx1_1_1[i,c]<- 1-exp(-((y.step_1_1[c]/mu.df_1_1[i])^(-b_1_1))) dx2_1_1[i,c]<- 1-exp(-((y.step_plus1_1_1[c]/mu.df_1_1[i])^(-b_1_1))) esw_pred_1_1[i,c+1]<-esw_pred_1_1[i,c]+0.5*step.size_1_1*(dx1_1_1[i,c]+dx2_1_1[i,c]) } esw_1_1[i]<-esw_pred_1_1[i,steps] #Ones Trick L.f0_1_1[i] <- (1-exp(-((dist_1_1[i]/ mu.df_1_1[i])^(-b_1_1)))) * 1/esw_1_1[i] #hazard rate likelihood phi_1_1[i]<- L.f0_1_1[i]/C ones.dist_1_1[i] ~ dbern(phi_1_1[i]) #NE Ship MR component logit(p1_1[i])<-b1.0_1+b.dist_1*dist_1_1[i] logit(p2_1[i])<-b1.0_1+b.dist_1*dist_1_1[i] logit(p1.0_1[i])<-b1.0_1 logit(p2.0_1[i])<-b1.0_1 p.dot_1[i]<-p1_1[i]+p2_1[i]-p1_1[i]*p2_1[i] p0_1[i]<-p1.0_1[i]+p2.0_1[i]-p1.0_1[i]*p2.0_1[i] pi_1.1[i]<-p1_1[i]*(1-p2_1[i])/p.dot_1[i] #y10 pi_1.2[i]<-(1-p1_1[i])*p2_1[i]/p.dot_1[i] #y01 pi_1.3[i]<-p1_1[i]*p2_1[i]/p.dot_1[i] #y11 Obs_1_1[i,1]~dbern(pi_1.1[i]) Obs_1_1[i,2]~dbern(pi_1.2[i]) Obs_1_1[i,3]~dbern(pi_1.3[i]) } mean.esw_1_1 <- mean(esw_1_1[]) mean.p0_1_1 <- mean(p0_1[]) Pcds_1_1<-mean.esw_1_1/W_1_1 P.det_1_1<- mean.p0_1_1*Pcds_1_1 #SE Shipboard for(i in 1:N_Sights_2_1){ ############################################################################ #USE HALF-NORMAL KEY ########################################################################### mu.df_2_1[i] <- b.df.0_2_1+b.df.beaufort_2_1*Beauf_2_1[i] ########## # estimate of sd and var, given coeffs above sig.df_2_1[i] <- exp(mu.df_2_1[i]) sig2.df_2_1[i] <- sig.df_2_1[i]*sig.df_2_1[i] # estimate effective strip width... # ...(esw, which is also the integral from 0 to W of g(y)) # argument of 'error func' in the integrand for detect func g(y) x_2_1[i] <-(W_2_1/sqrt(2*sig2.df_2_1[i])) # approximation of the 'error function' erf_2_1[i] <- sqrt(1-exp(-x_2_1[i]*x_2_1[i]*(4/pi + a*x_2_1[i]*x_2_1[i])/(1 + a*x_2_1[i]*x_2_1[i]))) # effective strip width for Shipboard esw_2_1[i] <- sqrt(pi * sig2.df_2_1[i] / 2) * erf_2_1[i] #Ones Trick L.f0_2_1[i] <- exp(-dist_2_1[i]*dist_2_1[i] / (2*sig2.df_2_1[i])) * 1/esw_2_1[i] #half-normal likelihood phi_2_1[i]<- L.f0_2_1[i]/C ones.dist_2_1[i] ~ dbern(phi_2_1[i]) #SE Ship MR component logit(p1_2[i])<-b1.0_2+b.dist1_2*dist_1_2[i]+b.beauf_2*Beauf_2_1[i] logit(p2_2[i])<-b2.0_2+b.dist2_2*dist_1_2[i]+b.beauf_2*Beauf_2_1[i] logit(p1.0_2[i])<-b1.0_2+b.beauf_2*Beauf_2_1[i] logit(p2.0_2[i])<-b2.0_2+b.beauf_2*Beauf_2_1[i] p.dot_2[i]<-p1_2[i]+p2_2[i]-p1_2[i]*p2_2[i] p0_2[i]<-p1.0_2[i]+p2.0_2[i]-p1.0_2[i]*p2.0_2[i] pi_2.1[i]<-p1_2[i]*(1-p2_2[i])/p.dot_2[i] #y10 pi_2.2[i]<-(1-p1_2[i])*p2_2[i]/p.dot_2[i] #y01 pi_2.3[i]<-p1_2[i]*p2_2[i]/p.dot_2[i] #y11 Obs_2_1[i,1]~dbern(pi_2.1[i]) Obs_2_1[i,2]~dbern(pi_2.2[i]) Obs_2_1[i,3]~dbern(pi_2.3[i]) } mean.esw_2_1 <- mean(esw_2_1[]) mean.p0_2_1 <- mean(p0_2[]) Pcds_2_1<-mean.esw_2_1/W_2_1 P.det_2_1<- mean.p0_2_1*Pcds_2_1 #NE Aerial for(i in 1:N_Sights_1_2){ ############################################################################################## # USE HAZARD RATE KEY ############################################################################################ mu.df_1_2[i] <- b.df.0_1_2 esw_pred_1_2[i,1]<-0 #Numerically integrate hazard rate function for (c in 1:(steps-1)){ dx1_1_2[i,c]<- 1-exp(-((y.step_1_2[c]/mu.df_1_2[i])^(-b_1_2))) dx2_1_2[i,c]<- 1-exp(-((y.step_plus1_1_2[c]/mu.df_1_2[i])^(-b_1_2))) esw_pred_1_2[i,c+1]<-esw_pred_1_2[i,c]+0.5*step.size_1_2*(dx1_1_2[i,c]+dx2_1_2[i,c]) } esw_1_2[i]<-esw_pred_1_2[i,steps] #Ones Trick L.f0_1_2[i] <- (1-exp(-((dist_1_2[i]/ mu.df_1_2[i])^(-b_1_2)))) * 1/esw_1_2[i] #hazard rate likelihood phi_1_2[i]<- L.f0_1_2[i]/C ones.dist_1_2[i] ~ dbern(phi_1_2[i]) } mean.esw_1_2 <- mean(esw_1_2[]) Pcds_1_2<-mean.esw_2_1/W_1_2 P.det_1_2<- g.knot[1]*Pcds_1_2 #SE Aerial for(i in 1:N_Sights_2_2){ ############################################################################ #USE HALF-NORMAL KEY ########################################################################### mu.df_2_2[i] <- b.df.0_2_2+b.df.beaufort_2_2*Beauf_2_2[i] ########## # estimate of sd and var, given coeffs above sig.df_2_2[i] <- exp(mu.df_2_2[i]) sig2.df_2_2[i] <- sig.df_2_2[i]*sig.df_2_2[i] # estimate effective strip width... # ...(esw, which is also the integral from 0 to W of g(y)) # argument of 'error func' in the integrand for detect func g(y) x_2_2[i] <-( W_2_2/sqrt(2*sig2.df_2_2[i])) # approximation of the 'error function' erf_2_2[i] <- sqrt(1-exp(-x_2_2[i]*x_2_2[i]*(4/pi + a*x_2_2[i]*x_2_2[i])/(1 + a*x_2_2[i]*x_2_2[i]))) # effective strip width for Shipboard esw_2_2[i] <- sqrt(pi * sig2.df_2_2[i] / 2) * erf_2_2[i] #Ones Trick L.f0_2_2[i] <- exp(-dist_2_2[i]*dist_2_2[i] / (2*sig2.df_2_2[i])) * 1/esw_2_2[i] #half-normal likelihood phi_2_2[i]<- L.f0_2_2[i]/C ones.dist_2_2[i] ~ dbern(phi_2_2[i]) } mean.esw_2_2 <- mean(esw_2_2[]) Pcds_2_2<-mean.esw_2_2/W_2_2 P.det_2_2<- g.knot[2]*Pcds_2_2 ################################################################################################################### #Priors for GAM smooth terms #Priors for CPG parameters gamma1 ~ dunif(-100,100) gammaCHL ~ dunif(-100,100) gammaBTEMP ~ dunif(-100,100) gammaDIST125 ~ dunif(-100,100) gammaLAT ~ dunif(-100,100) b.plusgamma[1]<-b[1]+gamma1 b.plusgamma[2:5]<-b[2:5]+gammaCHL b.plusgamma[6:9]<-b[6:9]+gammaBTEMP b.plusgamma[10:13]<-b[10:13]+gammaDIST125 b.plusgamma[14:17]<-b[14:17]+gammaLAT b.minusgamma[1]<-b[1]-gamma1 b.minusgamma[2:5]<-b[2:5]-gammaCHL b.minusgamma[6:9]<-b[6:9]-gammaBTEMP b.minusgamma[10:13]<-b[10:13]-gammaDIST125 b.minusgamma[14:17]<-b[14:17]-gammaLAT eta <- X %*% b ## linear predictor eta.plusgamma<-X %*% b.plusgamma eta.minusgamma<-X %*% b.minusgamma ## GAM priors from jagam function in mgcv for (i in 1:1) { b[i] ~ dnorm(0,0.0021) } ## prior for s(CHL)... K1 <- S1[1:4,1:4] * lambda[1] b[2:5] ~ dmnorm(zero[2:5],K1) ## prior for s(BTEMP)... K2 <- S2[1:4,1:4] * lambda[2] b[6:9] ~ dmnorm(zero[6:9],K2) ## prior for s(DIST125)... K3 <- S3[1:4,1:4] * lambda[3] b[10:13] ~ dmnorm(zero[10:13],K3) ## prior for s(LAT)... K4 <- S4[1:4,1:4] * lambda[4] b[14:17] ~ dmnorm(zero[14:17],K4) ## smoothing parameter priors CHECK... for (i in 1:4) { lambda[i] ~ dgamma(.05,.005) rho[i] <- log(lambda[i]) } #Habitat Loop for (t in 1:N_t) { ############################################################################## #NE Shipboard ############ # USE HAZARD RATE KEY ###### new_mu.df_1_1[t] <- b.df.0_1_1 new_esw_pred_1_1[t,1]<-0 #Numerically integrate hazard rate function for (b in 1:(steps-1)){ new_dx1_1_1[t,b]<- 1-exp(-((y.step_1_1[b]/new_mu.df_1_1[t])^(-b_1_1))) new_dx2_1_1[t,b]<- 1-exp(-((y.step_plus1_1_1[b]/new_mu.df_1_1[t])^(-b_1_1))) new_esw_pred_1_1[t,b+1]<-new_esw_pred_1_1[t,b]+0.5*step.size_1_1*(new_dx1_1_1[t,b]+new_dx2_1_1[t,b]) } new_esw_1_1[t]<-new_esw_pred_1_1[t,steps] #Calculate g.knot NE new_f0_1_1[t] <- 1/new_esw_1_1[t] P.0.1_1_1[t]<-(exp(b1.0_1))/(1+exp(b1.0_1)) #NE Ship P.0.2_1_1[t]<-(exp(b1.0_1))/(1+exp(b1.0_1)) #NE Ship P_0_1_1[t]<- P.0.1_1_1[t]+ P.0.2_1_1[t]- P.0.1_1_1[t]*P.0.2_1_1[t] #Availability Av_1_1[t]<-1 #NE Shipboard P_t_1_1[t]<- P_0_1_1[t]*new_esw_1_1[t]/W_1_1* Av_1_1[t] #Overall detection for transect t ################################################################################################################ #SE Shipboard ################ # USE HALF-NORMAL RATE KEY ####################### new_mu.df_2_1[t] <- (b.df.0_2_1+ b.df.beaufort_2_1*Beaufort[t]) new_sig.df_2_1[t] <- exp(new_mu.df_2_1[t]) new_sig2.df_2_1[t] <- new_sig.df_2_1[t]*new_sig.df_2_1[t] new_x_2_1[t] <- (W_2_1/sqrt(2*new_sig2.df_2_1[t])) new_erf_2_1[t] <- sqrt(1-exp(-new_x_2_1[t]*new_x_2_1[t]*(4/pi + a*new_x_2_1[t]*new_x_2_1[t])/(1 + a*new_x_2_1[t]*new_x_2_1[t]))) new_esw_2_1[t] <- sqrt(pi * new_sig2.df_2_1[t] / 2) * new_erf_2_1[t] #Calculate g.knot SE Shipboard P.0.1_2_1[t]<-(exp(b1.0_2+b.beauf_2*Beaufort[t]))/(1+exp(b1.0_2+b.beauf_2*Beaufort[t])) #SE Ship P.0.2_2_1[t]<-(exp(b2.0_2+b.beauf_2*Beaufort[t]))/(1+exp(b2.0_2+b.beauf_2*Beaufort[t])) #SE Ship P_0_2_1[t]<- P.0.1_2_1[t]+ P.0.2_2_1[t]- P.0.1_2_1[t]*P.0.2_2_1[t] #Surface Availability Shipboard Av_2_1[t]<-1 #SE Shipboard #Overall detection for transect t SE Shipboard P_t_2_1[t]<- P_0_2_1[t]*new_esw_2_1[t]/W_2_1* Av_2_1[t] #################################################################################################################### # NE Aerial ###################### # USE HAZARD RATE KEY ##################### new_mu.df_1_2[t] <- b.df.0_1_2 new_esw_pred_1_2[t,1]<-0 #Numerically integrate hazard rate function for (b in 1:(steps-1)){ new_dx1_1_2[t,b]<- 1-exp(-((y.step_1_2[b]/new_mu.df_1_2[t])^(-b_1_2))) new_dx2_1_2[t,b]<- 1-exp(-((y.step_plus1_1_2[b]/new_mu.df_1_2[t])^(-b_1_2))) new_esw_pred_1_2[t,b+1]<-new_esw_pred_1_2[t,b]+0.5*step.size_1_2*(new_dx1_1_2[t,b]+new_dx2_1_2[t,b]) } new_esw_1_2[t]<-new_esw_pred_1_2[t,steps] #Surface Availability Aerial Av_1_2[t]<-Av #Overall detection for transect t for NE Aerial P_t_1_2[t]<- g.knot[1]*new_esw_1_2[t]/W_1_2* Av_1_2[t] ################################################################ ############################################################################################## # USE HALF-NORMAL RATE KEY ############################################################################################ new_mu.df_2_2[t] <- b.df.0_2_2 + b.df.beaufort_2_2*Beaufort[t] new_sig.df_2_2[t] <- exp(new_mu.df_2_2[t]) new_sig2.df_2_2[t] <- new_sig.df_2_2[t]*new_sig.df_2_2[t] new_x_2_2[t] <- (W_2_2/sqrt(2*new_sig2.df_2_2[t])) new_erf_2_2[t] <- sqrt(1-exp(-new_x_2_2[t]*new_x_2_2[t]*(4/pi + a*new_x_2_2[t]*new_x_2_2[t])/(1 + a*new_x_2_2[t]*new_x_2_2[t]))) new_esw_2_2[t] <- sqrt(pi * new_sig2.df_2_2[t] / 2) * new_erf_2_2[t] new_f0_2_2[t] <- 1/new_esw_2_2[t] #Surface Availability Aerial Av_2_2[t]<-Av P_t_2_2[t]<- g.knot[2]*new_esw_2_2[t]/W_2_2* Av_2_2[t] #Overall detection for transect t #Search Effort shipboard surveys effort_offset_1_1[t]<-(Effort[t]*2*W_1_1)/Area[t] effort_offset_2_1[t]<-(Effort[t]*2*W_2_1)/Area[t] #Search Effort aerial surveys effort_offset_1_2[t]<-(Effort[t]*2*W_1_2)/Area[t] effort_offset_2_2[t]<-(Effort[t]*2*W_2_2)/Area[t] P_t[t]<-P_t_1_1[t]*Indicator_1_1[t]+P_t_1_2[t]*Indicator_1_2[t]+P_t_2_1[t]*Indicator_2_1[t]+P_t_2_2[t]*Indicator_2_2[t] effort_offset[t]<-effort_offset_1_1[t]*Indicator_1_1[t]+effort_offset_1_2[t]*Indicator_1_2[t]+effort_offset_2_1[t]*Indicator_2_1[t]+ effort_offset_2_2[t]*Indicator_2_2[t] D[t]<-Sightings[t]/(P_t[t]*effort_offset[t])*(lambda.size+1) mu[t] <- exp(eta[t]) ## expected response mu_g[t] <- exp(eta.plusgamma[t]/2) mu_p[t] <- exp(eta.minusgamma[t]/2) eff.mu_p[t]<-mu_p[t]*P_t[t]*effort_offset[t] eff.mu[t]<-mu[t]*P_t[t]*effort_offset[t] ################################################################ } #Compound Poisson-Gamma (CPG) Likelihood for (i in 1:Nnonzero){ Sightings[NonZeroObs[i]] ~ dgamma(psi[NonZeroObs[i]] * nu , nu / mu_g[NonZeroObs[i]]) psi[NonZeroObs[i]] ~ dpois(eff.mu_p[NonZeroObs[i]]) } for (i in 1:Nzero){ Sightings[ZeroObs[i]] ~ dpois(eff.mu_p[ZeroObs[i]]) psi[ZeroObs[i]] ~ dpois(eff.mu_p[ZeroObs[i]]) } } # end model ",fill=TRUE) sink() # Bundle data Model.data <- list( N_Sights_1_1=dim(Dat_Sightings_1_1)[1], #NE Shipboard N_Sights_2_1=dim(Dat_Sightings_2_1)[1], #SE Shipboard N_Sights_1_2=dim(Dat_Sightings_1_2)[1], #NE Aerial N_Sights_2_2=dim(Dat_Sightings_2_2)[1], #SE Aerial Sightings=Sightings, Obs_1_1=Dat_Sightings_1_1[,10:12], #Double platform observer data NE Ship Obs_2_1=Dat_Sightings_2_1[,10:12], #Double platform observer data SE Ship C=10000, #For implememnting "Ones Trick" in JAGS ones.dist_1_1=ones.dist_1_1, ones.dist_1_2=ones.dist_1_2, ones.dist_2_1=ones.dist_2_1, ones.dist_2_2=ones.dist_2_2, W_1_1=W_1_1/1000, #Truncation distance NE Ship. Convert to KM. W_2_1=W_2_1/1000, #Truncation distance SE Ship. Convert to KM. W_1_2=W_1_2/1000, #Truncation distance NE Aerial. Convert to KM. W_2_2=W_2_2/1000, #Truncation distance SE Aerial. Convert to KM. N_t=N_t, #Number of unique grid cells gs=gs-1, #Observed group sizes (subtract 1 for truncated Poisson) N.size=N.size, dist_1_1=c(Dat_Sightings_1_1[,6])/1000+0.0001, #column with distance data NE Ship. Add a little to avoid zeros. dist_2_1=c(Dat_Sightings_2_1[,6])/1000+0.0001, #column with distance data SE Ship. Add a little to avoid zeros. dist_1_2=c(Dat_Sightings_1_2[,6])/1000+0.0001, #column with distance data NE Aerial. Add a little to avoid zeros. dist_2_2=c(Dat_Sightings_2_2[,6])/1000+0.0001, #column with distance data SE Aerial. Add a little to avoid zeros. Beauf_2_1=c(Dat_Sightings_2_1[,7]), Beauf_2_2=c(Dat_Sightings_2_2[,7]), Area=Final_Cov_Mat[,1], Effort=Final_Cov_Mat[,8], Beaufort=Final_Cov_Mat[,9], Indicator_1_1=Indicator_1_1, Indicator_1_2=Indicator_1_2, Indicator_2_1=Indicator_2_1, Indicator_2_2=Indicator_2_2, a_g0=a_g0, b_g0=b_g0, a_Av=a_Av, b_Av=b_Av, #Inputs for integrating hazard rate function steps=steps, y.step_1_1=y.step_1_1, y.step_plus1_1_1=y.step_plus1_1_1, y.step_1_2=y.step_1_2, y.step_plus1_1_2=y.step_plus1_1_2, step.size_1_1=step.size_1_1, step.size_1_2=step.size_1_2, #GAM inputs X=X, S1=S1, S2=S2, S3=S3, S4=S4, zero=zero, #CPG inputs NonZeroObs=NonZeroObs, Nnonzero=Nnonzero, ZeroObs=ZeroObs, Nzero=Nzero ) # Initial values inits <- function(){list(b.df.0_1_1 = runif(1,1,3), b.df.0_1_2 = runif(1,1,3), b.df.0_2_1 = runif(1,1,3), b.df.0_2_2 = runif(1,1,3), b.df.beaufort_2_1 = runif(1,-1,1), b.df.beaufort_2_2 = runif(1,-1,1), b_1_1 = runif(1,0,1), b_2_1 = runif(1,0,1), b1.0_1 = runif(1,0,1), b1.0_2 = runif(1,0,1), b2.0_2 = runif(1,0,1), b.dist_1= runif(1,0,1), b.dist1_2= runif(1,0,1), b.dist2_2= runif(1,0,1), b.beauf_2=runif(1,0,1), psi=as.numeric(Sightings != 0), gammaBTEMP=0, gammaCHL=0, gammaLAT=0, gammaDIST125=0, #GAM inits b=b, lambda=lambda ) } # Define parameters to be monitored # varsToMonitor<-c("b.df.0_1_1", "b_1_1", "b.dist_1", "b1.0_1", "b.df.0_1_2", "b.df.beaufort_2_1", "b1.0_2","b2.0_2","b.dist1_2","b.dist2_2","b.beauf_2", "b.df.0_2_1", "b_2_1", "b.df.0_2_2", "b.df.beaufort_2_2", "mean.esw_1_1", "mean.esw_1_2","mean.esw_2_1","mean.esw_2_2", "P_0_1_1","P_0_1_2","P_0_2_1","P_0_2_2","P_t_1_1","P_t_1_2","P_t_2_1","P_t_2_2", "mu", "eta","y.rep", "fit.pred", "fit.true", "lambda.size", "b0.size", "lambda", "b" , "rho", "g.knot", "Av", "gamma", "fit.pred", "fit.true", "D", "r", "eff.mu", "lam", "mean.p0_1_1", "Pcds_1_1", "P.det_1_1", "mean.p0_1_2","Pcds_1_2", "P.det_1_2","Pcds_2_1", "P.det_2_1","Pcds_2_2", "P.det_2_2") require(rjags) nitt=2000 burnin=1000 chains=1 thin=2 jm <- jags.model("Fin_jagam_CPG.bugs", data = Model.data, inits = inits, n.adapt = burnin, n.chains = chains) sam <- jags.samples(jm, varsToMonitor, n.iter = nitt, thin = thin) save(sam, file="Fin_jagam_CPG_output.RData")