Minor changes and some additional annotations

JamieBrusa · Dec 29, 2023 · ef80ce0 · ef80ce0
1 parent 3f92b40
commit ef80ce0
Show file tree

Hide file tree

Showing 3 changed files with 33 additions and 19 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/NegBinomBreedMolt.txt b/NegBinomBreedMolt.txt
@@ -15,6 +15,7 @@ model{
 	}
 
 
+#set up priors
 tau.eps.year ~ dgamma(1,1) 
 tau.eps.pair ~ dgamma(1,1)
 tau.eps.psu ~ dgamma(1,1)
@@ -44,14 +45,15 @@ beta.chl~dnorm(0, 0.1)
 beta.BM~dnorm(0,0.1)
 beta.offshore~dnorm(0,0.1)
 alpha~dnorm(0,.1)
-sigma0~dunif(0, 7)
+sigma0~dunif(0,7)
 sigma.eps.year <- 1/sqrt(tau.eps.year)
 sigma.eps.pair <- 1/sqrt(tau.eps.pair)
 sigma.eps.psu <- 1/sqrt(tau.eps.psu)
 
 r.N~dunif(0,100)
 
 
+#the observation model
 
 		for (j in 1:nsites){
 		log(sigma[j]) <- sigma0 + beta.bss.1*BSS.1[j] + beta.bss.2*BSS.2[j] + beta.bss.3*BSS.3[j] + 
@@ -67,27 +69,32 @@ r.N~dunif(0,100)
 
  	pcap[j]<-sum(f[j,1:nG]) 
 
+#the observed counts   	    
+y[j]~ dbin(pcap[j],N[j])
+
+#the abundance model
 
 	    lambda[j] <- exp(alpha + beta.ho*haul.count[j] + beta.NPGO*NPGO[j] + beta.Depth*MeanDepth[j] + 
 beta.upwell*Upwell[j] +
  eps.year[Year[j]] + beta.ST*ST[j] + beta.FT*FT[j] + beta.shore7*shore7[j] + beta.shore5*shore5[j] + 
                beta.shore9A*shore9A[j] + beta.shore1A*shore1A[j] + beta.shore2A*shore2A[j] + eps.PSU[SP[j]] +
 beta.shore6A*shore6A[j] + beta.shore4*shore4[j] + beta.shore6D*shore6D[j] + beta.shore8A*shore8A[j] + beta.sal*sal[j] + 
 beta.sst*sst[j] + beta.chl*chl[j] + beta.river*River[j] + beta.BM*breedmolt[j] + beta.offshore*Offshore[j] + log(area[j]))
-     	    y[j]~ dbin(pcap[j],N[j])
-#	    N[j]~dpois(lambda[j])
+
+
 
 
-###stuff for negbin
+
+##Poisson-gamma formulation of the negative binomial
      	    N[j]~dpois(lambda.star[j])
 	    lambda.star[j]<-lambda[j]*rho.N[j]
 	    rho.N[j]~dgamma(r.N,r.N)
 
-	###create replicate abundances
+###create replicate abundances for Bayesian p-value analysis
 	Nnew[j]~dpois(lambda.star[j])
 #	Nnew[j]~dpois(lambda[j])
 
-	### create fit statistic 1
+### create fit statistic 1
 	FT1[j]<-pow(sqrt(N[j])-sqrt(lambda[j]),2)
 	FT1new[j]<-pow(sqrt(Nnew[j])-sqrt(lambda[j]),2)
 }
@@ -97,15 +104,15 @@ beta.sst*sst[j] + beta.chl*chl[j] + beta.river*River[j] + beta.BM*breedmolt[j] +
 	T1newp<-sum(FT1new[1:nsites])
 
 
-Bp.N<-sum(T1newp)>sum(T1p)
+        Bp.N<-sum(T1newp)>sum(T1p)
 
+##distance classes for observation model and Bayesian p-value fit statistics
       for(i in 1:nind){
       dclass[i] ~ dcat(fct[site[i], 1:nG]) 
       dclassnew[i] ~ dcat(fct[site[i], 1:nG]) 
       Tobsp[i]<- pow(1- sqrt(fct[site[i], dclass[i]]),2)
       Tobspnew[i]<- pow(1- sqrt(fct[site[i], dclassnew[i]]),2)
      }
-
 Bp.Obs<-sum(Tobspnew[1:nind])>sum(Tobsp[1:nind])
 
 

diff --git a/SS HDS model final.R b/SS HDS model final.R
@@ -80,12 +80,10 @@ multiseal_df$dclass<-with(multiseal_df,ifelse(Perp_Dist<=dist.breaks[1],1,ifelse
                                                                                  12))))))))))))
 
 
-##negative binomial dispersion parameter in the Poisson-Gamma formulation
-r <- 1
 
 
 ###########Full model##############
-#Make an array with 1 species x 129 transects x 12 distance classes and populate with N
+#Make an array with 1 species x 129 transects x 12 distance classes and populate with N (number of observations)
 transect <- unique(multiseal_df$StratPSUSegYJ)
 nSites <- length(unique(multiseal_df$StratPSUSegYJ))
 distclass <- unique(multiseal_df$dclass)
@@ -323,7 +321,7 @@ dataCovs<-list(nG=nG, xg=dist.breaks[-1]-13.5, nsites=nSites, Year = Year, pair
 ### initial values for N
 N.in<-t(seals.sum)+3
 
-initsCovs<-function(){list(N=as.vector(N.in), alpha = runif(1,1, 3), sigma0 = runif(1, 0, 7),
+initsCovs<-function(){list(N=as.vector(N.in), alpha = runif(1,1, 3), sigma0 = runif(1, 5, 6),
                            beta.bss.1=rnorm(1, 0, 0.1),
                            beta.bss.2=rnorm(1, 0, 0.1),
                            beta.bss.3=rnorm(1, 0, 0.1),
@@ -359,15 +357,28 @@ setwd(here())
 modelFileCovs='NegBinomBreedMolt.txt'
 
 
+
+
+
+# Set up to run chains in parallel
+library(parallel)
+cl <- makeCluster(3)
+start.time=Sys.time()
+
 covs.mod<-run.jags(model = modelFileCovs,
                    monitor = params.Covs,
                    data = dataCovs,
                    n.chains = 3,
-                   adapt = 2000,
                    burnin = 10000,
                    sample = 50000,
+                   adapt = 2000,
                    inits = initsCovs,
-                   thin = 1)
+                   thin = 1, method="rjparallel", cl=cl)
+
+stopCluster(cl)
+
+end.time=Sys.time()
+end.time-start.time
 
 
 #Diagnostics
@@ -386,9 +397,5 @@ ggs_density(covs.mod_ggs, c("beta.b"))
 ggs_density(covs.mod_ggs, c("alpha"))
 ggs_density(covs.mod_ggs, c("N"))
 ggs_density(covs.mod_ggs, c("sigma"))
-
-
-
-
-
+ggs_density(covs.mod_ggs, c("Bp.N"))