2024 update to NARW population model, including calf integration

NARW_JS_2model.R

@@ -1,8 +1,9 @@
 # code for fitting Jolly Seber model to North Atlantic right whale sightings (1990-)
 
 # Load the data/inputs----
+Year.Max <- 2023
 #source("NARW_JS_1data.R")
-load("NARW_JS_inputs.Rdata")
+load(paste0("./data/NARW_JS_inputs_",Year.Max,".Rdata"))
 library(parallel)
 library(MCMCvis)
 library(mcmcplots)
@@ -12,7 +13,7 @@ library(nimble)
 obs <- whale.data$y[,-1]
 obs[obs==2] <- 0 #turn to 0/1s
 obs <- obs[which(apply(obs,1,sum) > 0),]  #limit to known individuals
-write.csv(obs,paste0("sightings_NARW_1990-",Year.Max,".csv"),row.names = F)
+write.csv(obs,paste0("./data/sightings_NARW_1990-",Year.Max,".csv"),row.names = F)
 
 # Parameters monitored
 params <- c(
@@ -22,7 +23,7 @@ params <- c(
   "pi.female", "psi",             
   "sigma.phi", "sigma.p.i", "sigma.p.t",
   "epsilon.phi.t","epsilon.p.t",
-  "gamma", "entry", "N", "NF","NM", "B", "Nd","aliveT"
+  "gamma", "entry", "N", "NF","NM", "NadF", "B", "Nd","aliveT"
   )
 
 # multiple chain inits
@@ -40,7 +41,8 @@ run_MCMC_allcode <- function(info, whale.data, whale.constants, params,
                              n.iter, n.burnin, n.thin){
   library(nimble)
 
-  source("Pace2017_NARW_model.txt")
+  source("NARW_model_w_calfs.txt")
+  #source("Pace2017_NARW_model.txt")
 
   # testing only, not when running cluster
   # n.iter=2000; n.burnin=1000; n.thin=1
@@ -96,7 +98,8 @@ save(outL, run_MCMC_allcode, whale.data, whale.constants, #whale.tbl,
 # Examine output----
 load(file = paste0("./out/outL_NARW_1990-",Year.Max,".Rdata"))
 
-param.rm <- paste(c("N","B","entry","epsilon","gamma","aliveT"),collapse="|")
+param.rm <- paste(c("N", "NF","NM", "NadF","Nd", "B",
+                    "entry","epsilon","gamma","aliveT"),collapse="|")
 
 MCMCsummary(outL, round = 3, 
             excl = colnames(outL[[1]])[grep(pat = param.rm,colnames(outL[[1]]))])

NARW_model_w_calfs.txt

@@ -0,0 +1,115 @@
+
+whale.code.pop <- nimbleCode( {
+
+#model {
+
+  pi.female ~ dbeta(5, 5)      ### probability of individual being female
+
+  a0 ~ dlogis(0,1)             ### logit-scale intercept for sighting prob (p)
+  a.female ~ dunif(-5, 5)
+
+  b0 ~ dlogis(0,1)             ### logit-scale intercept for survival (phi)
+  b.female ~ dunif(-5, 5)
+  b.age ~ dunif(-5, 5)
+  b.regime ~ dunif(-5, 5)
+
+  sigma.p.t ~ dunif(0, 10)       ###  prior sd of random year effect on p
+  sigma.p.i ~ dunif(0, 10)       ###  prior sd of random individual effect on p
+  sigma.phi ~ dunif(0, 10)       ###  prior sd of random year effect on phi
+
+  phi_calf ~ dbeta(1, 1)
+
+  for (i in 1:(M-1)) {
+    epsilon.p.i[i] ~ dnorm(0, sd=sigma.p.i) ### random individual effect on p
+  }
+  for (t in 1:(n.occasions - 2)) {
+    epsilon.phi.t[t] ~ dnorm(0, sd=sigma.phi)  ### random year effect on phi
+    epsilon.p.t[t] ~ dnorm(0, sd=sigma.p.t)    ### random year effect on p
+  }
+  #### sum-to-zero epsilons
+  epsilon.p.i[M] <- -sum(epsilon.p.i[1:(M-1)])  
+  epsilon.phi.t[n.occasions-1] <- -sum(epsilon.phi.t[1:(n.occasions-2)])
+  epsilon.p.t[n.occasions-1] <- -sum(epsilon.p.t[1:(n.occasions-2)])
+
+  # Prior for entry probabilities
+  gamma[1] ~ dunif(0, 1)
+  for (t in 2:(n.occasions - 1)) {
+    gamma[t] <- (calves[t]*phi_calf) / (M*prod(1-gamma[1:(t-1)]))                  
+  } #t
+
+  #########  Probability models
+  for (i in 1:M) {
+    Sex[i] ~ dbern(pi.female)
+    for (t in 1:(n.occasions - 1)) {
+      logit(phi[i,t]) <- b0 + b.age * Age[i,t] + b.female * Sex[i] * Adult[i,t] + b.regime *
+        Regime[t] + epsilon.phi.t[t]
+      logit(p[i,t]) <- a0 + a.female * (Sex[i]) + epsilon.p.t[t] + epsilon.p.i[i]
+    } #t
+  } #i
+  ### Define state-transition and observation matrices
+  for (i in 1:M) {
+    # Given state S(t), define probabilities of S(t+1) 
+    for (t in 1:(n.occasions - 1)) {
+      ps[1,i,t,1] <- 1 - gamma[t]
+      ps[1,i,t,2] <- gamma[t]
+      ps[1,i,t,3] <- 0
+      ps[2,i,t,1] <- 0
+      ps[2,i,t,2] <- phi[i,t]
+      ps[2,i,t,3] <- 1 - phi[i,t]
+      ps[3,i,t,1] <- 0
+      ps[3,i,t,2] <- 0
+      ps[3,i,t,3] <- 1
+      # Given state S(t), define probabilities of O(t)
+      po[1,i,t,1] <- 0
+      po[1,i,t,2] <- 1
+      po[2,i,t,1] <- p[i,t]*(1-equals(z[i,n.occasions - 1], 1))
+      po[2,i,t,2] <- 1 - (p[i,t]*(1-equals(z[i,n.occasions - 1], 1)))
+      po[3,i,t,1] <- 0
+      po[3,i,t,2] <- 1
+    } #t
+  } #i
+
+  # Likelihood
+  for (i in 1:M) {
+    # Define latent state at first occasion.
+    z[i, 1] ~ dcat(px[1:3])   # all M individuals are known z=1 at t=1
+    for (t in 2:n.occasions) {
+      # State process: draw S(t) given S(t-1)
+      z[i,t] ~ dcat(ps[z[i,t-1],i,t-1,1:3])
+      # Observation process: draw O(t) given S(t)
+      y[i,t] ~ dcat(po[z[i,t],i,t-1,1:2])
+    } #t
+  } #i
+  # Calculate derived population parameters
+  # for (t in 1:(n.occasions - 1)) {
+  #   qgamma[t] <- 1 - gamma[t]
+  # }
+  # cprob[1] <- gamma[1]
+  # for (t in 2:(n.occasions - 1)) {
+  #   cprob[t] <- gamma[t] * prod(qgamma[1:(t-1)])
+  # } #t
+  # psi <- sum(cprob[1:(n.occasions-1)])            # Inclusion probability
+  # for (t in 1:(n.occasions - 1)) {
+  #   entry[t] <- cprob[t] / psi      # Entry probability
+  # } #t
+  for (i in 1:M) {
+    for (t in 2:n.occasions) {
+      al[i,t-1] <- equals(z[i,t], 2)
+      dead[i,t-1] <- equals(z[i,t], 3) * equals(z[i,t-1], 2)
+    } #t
+    for (t in 1:(n.occasions - 1)) {
+      d[i,t] <- equals(z[i,t] - al[i,t], 0)
+    } #t
+    aliveT[i] <- al[i,(n.occasions-1)]
+  } #i
+  for (t in 1:(n.occasions - 1)) {
+    N[t] <- sum(al[1:M,t])
+    NF[t] <- sum(al[1:M,t] * Sex[1:M])
+    NadF[t] <- sum(al[1:M,t] * Sex[1:M] * Proven[1:M,t]) 
+    NM[t] <- sum(al[1:M,t] * (1-Sex[1:M]))
+    Nd[t] <- sum(dead[1:M,t]) 
+    B[t] <- sum(d[1:M,t])         # Number of entries
+  } #t
+}
+)
+

README.md

@@ -1,5 +1,5 @@
 # Population size estimation for North Atlantic right whales
 
-This repository contains the data and scripts for fitting a Jolly-Seber model to North Atlantic right whale (NARW) sightings, using the general framework described by [Pace et al. (2017)](https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.3406).
+This repository contains the data and scripts for fitting a Jolly-Seber model to North Atlantic right whale (NARW) sightings, using the general framework described by [Pace et al. (2017)](https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.3406) and further refined by [Linden (2024)](https://www.biorxiv.org/content/10.1101/2024.10.11.617830v1).
 
 *This repository is a scientific product and is not official communication of the National Oceanic and Atmospheric Administration, or the United States Department of Commerce. All NOAA GitHub project code is provided on an 'as is' basis and the user assumes responsibility for its use. Any claims against the Department of Commerce or Department of Commerce bureaus stemming from the use of this GitHub project will be governed by all applicable Federal law. Any reference to specific commercial products, processes, or services by service mark, trademark, manufacturer, or otherwise, does not constitute or imply their endorsement, recommendation or favoring by the Department of Commerce. The Department of Commerce seal and logo, or the seal and logo of a DOC bureau, shall not be used in any manner to imply endorsement of any commercial product or activity by DOC or the United States Government.*