Resolve issue with "All nests scenario" coding in PVA.

R/01-data-manipulation.R

@@ -730,7 +730,6 @@ inits.func1 <- function (){
   r = mean(outp$r),
   N = outp$N[1:7,1:13,1:2, 
               sample(seq(1, 5200, by=400), 1, replace = F)] # sample from inits of chains that worked
-
   )}
 
 

R/06-ipm-simp.R → R/02-ipm-simp.R

@@ -1,13 +1,12 @@
-## ---- ipm1 --------
+## ---- ipm-simp --------
 library('nimble')
 library('parallel')
 library ('coda')
-#load("data/data.rdata")
 load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
 
-#################################
-# The model
-####################################################
+#**********************
+#* Parameter descriptions
+#**********************
 # Mark-resight-recovery data
 #   Observations (po) = y  
 #     1 seen first-year (age=0, just before 1st b-day)
@@ -33,13 +32,27 @@ load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
 #   pFY: resight probability first-years
 #   pA: resight probability nonbreeders
 #   pB: resight probability breeders
+#   lmu.prod: mean productivity (males and females) per territory on the log scale
+#   sds: standard deviations for multivariate normal random effects for time and site
+#   R: correlation coefficients for multivariate normal random effects for time and site
+#   lambda: population growth rate (derived)
+#   extinct: binary indicator of extirpation at a site
+#   gamma: coefficient of effect from nest treatments
+#   betas: coefficient of effect from translocations
+#   deltas: coefficient of effect from survey effort
+#   mus: overall means for survival, recruitment, and detections
+#   r and rr: "r" parameter for negative binomial distribution
+#             also described as omega in manuscript
 
+#**********************
+#* Model code
+#**********************
 mycode <- nimbleCode(
   {
     ###################################################
     # Priors and constraints
     ###################################################
-    # survival, recruitment, and detection can be correlated
+    # Survival, recruitment, and detection can be correlated
     for (k in 1:8){ 
       betas[k] ~ dunif(-20, 20)  # prior for coefficients
     } # k
@@ -49,51 +62,40 @@ mycode <- nimbleCode(
         mus[j,s] ~ dbeta(1,1) # prior for means
       }}     # m population #s sex #h hacked 
 
-    # Temporal random effects and correlations among all sites
-    for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation
-    # estimated using the multivariate normal distribution
+    # Temporal random effects and correlations among all sites, synchrony
+    for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
     R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
     Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
     U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
-    # multivariate normal for temporal variance
-    for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
+    # Multivariate normal for temporal variance
+    for (t in 1:nyr){ 
       eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
                           cholesky = U[1:p, 1:p], prec_param = 0)
     }
 
     # Temporal random effects and correlations between sites
-    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation
-    # estimated using the multivariate normal distribution
+    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
     R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
     Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
     U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
-    # multivariate normal for temporal variance
-    for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
+    # multivariate normal for temporal and site variance
+    for (t in 1:nyr){ 
       for (s in 1:nsite){
         eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
                                cholesky = U2[1:p2, 1:p2], prec_param = 0)
       } } # s t 
-    #######################
-    # Derived params
-    #######################
-    # for (s in 1:nsite){
-    #   for (t in 1:(nyr-1)){
-    #     lambda[t, s] <-  Ntot[t+1, s]/(Ntot[t, s])
-    #     loglambda[t, s] <- log(lambda[t, s])
-    #   }} #t
-
+
     ###############################
-    # Likelihood for fecundity
+    # Likelihood for productivity
     ###############################
-    # Priors for number of fledglings
-    # and nest success
+    # Priors for productivity
     for (s in 1:nsite){
       lmu.f[s] ~ dnorm(0, sd=5)
     } # s
     gamma ~ dunif(-20, 20)
     rr ~ dexp(0.05)
 
-    # Fecundity       
+    # Productivity likelihood      
     for (k in 1:nnest){
       f[k] ~ dnegbin(ppp[k], rr)
       ppp[k] <- rr/(rr+mu.f[k])
@@ -102,7 +104,10 @@ mycode <- nimbleCode(
         eps[9, year.nest[k] ] + 
         eta[9, site.nest[k], year.nest[k] ] 
     } # k
-    # derive yearly brood size for population model
+    # Derive yearly brood size for population model
+    # Need to reorder because nimble doesn't 
+    # handle nonconsecutive indices
+    # yrind.nest is a matrix of indices for each site
     for (t in 1:nyr){
       for (s in 1:nsite){
         for (xxx in 1:nest.end[t,s]){
@@ -127,132 +132,16 @@ mycode <- nimbleCode(
     ################################
     # Likelihood for counts
     ################################
-    # # Abundance for year=1
-    # for (v in 1:7){ 
-    #   for (s in 1:nsite){
-    #     # subtract one because to allow dcat to include zero
-    #     N[v, 1, s] <- N2[v, 1, s] - 1 
-    #     N2[v, 1, s] ~ dcat(pPrior[1:s.end[s], s]) 
-    #   }} # s t
-    # # Abundance for years > 1
-    # for (t in 1:(nyr-1)){
-    #   for (s in 1:nsite){
-    #     # Number of wild born juvs
-    #     N[1, t+1, s] ~ dpois( (NFY[t, s]*mn.phiFY[t, s]*mn.psiFYB[t, s] + # first year breeders
-    #                              NF[t, s]*mn.phiA[t, s]*mn.psiAB[t, s] + # nonbreeders to breeders
-    #                              NB[t, s]*mn.phiB[t, s]*(1-mn.psiBA[t, s])) # breeders remaining
-    #                           *mn.f[t+1, s] ) # end Poisson
-    #     # Abundance of nonbreeders
-    #     ## Second year nonbreeders
-    #     N[2, t+1, s] ~ dbin(mn.phiFY[t, s]*(1-mn.psiFYB[t, s]), NFY[t, s]) # Nestlings to second year nonbreeders
-    #     ## Adult nonbreeders
-    #     N[3, t+1, s] ~ dbin(mn.phiA[t, s]*(1-mn.psiAB[t, s]), NF[t, s]) # Nonbreeders to nonbreeders
-    #     N[4, t+1, s] ~ dbin(mn.phiB[t, s]*mn.psiBA[t, s], NB[t, s]) # Breeders to nonbreeders
-    #     # Abundance of breeders
-    #     ## Second year breeders
-    #     N[5, t+1, s] ~ dbin(mn.phiFY[t, s]*mn.psiFYB[t, s], NFY[t, s]) # Nestlings to second year breeders
-    #     ## Adult breeders
-    #     N[6, t+1, s] ~ dbin(mn.phiA[t, s]*mn.psiAB[t, s], NF[t, s]) # Nonbreeder to breeder
-    #     N[7, t+1, s] ~ dbin(mn.phiB[t, s]*(1-mn.psiBA[t, s]), NB[t, s]) # Breeder to breeder
-    #   }} # s t
-    # 
-    # # Observation process    
-    # for (t in 1:nyr){
-    #   for (s in 1:nsite){
-    #     counts[2, t, s] ~ dpois(NF[t, s]) # nonbreeding adult females 
-    #     counts[3, t, s] ~ dpois(NB[t, s]) # breeding females 
-    #     # constrain N1+hacked.counts to be >=0
-    #     constraint_data[t, s] ~ dconstraint( (N[1, t, s] + hacked.counts[t, s]) >= 0 ) # Transfers translocated first-year females
-    #     NFY[t, s] <- N[1, t, s] + hacked.counts[t, s] # Transfers translocated first-year females
-    #     NF[t, s] <- sum(N[2:4, t, s]) # number of adult nonbreeder females
-    #     NB[t, s] <- sum(N[5:7, t, s]) # number of adult breeder females
-    #     Ntot[t, s] <- sum(N[1:7, t, s]) # total number of females
-    #   }# s
-    #   # First-years at different sites have different distributions
-    #   # for better model fit
-    #   counts[1, t, 1] ~ dpois(N[1, t, 1]) # doesn't have any zeroes so poisson
-    #   counts[1, t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
-    #   pp[t] <- r/(r+N[1, t, 2])
-    # } # t
-    # r ~ dexp(0.05)
-    # ###################
-    # # Assess GOF of the state-space models for counts
-    # # Step 1: Compute statistic for observed data
-    # # Step 2: Use discrepancy measure: mean absolute error
-    # # Step 3: Use test statistic: number of turns
-    # ###################
-    # for (t in 1:nyr){
-    #   c.repFY[t, 1] ~ dpois( N[1, t, 1] )
-    #   c.repFY[t, 2] ~ dnegbin( pp[t], r )
-    #   for (s in 1:nsite){
-    #     c.expB[t, s] <- NB[t, s]  # expected counts adult breeder
-    #     c.expA[t, s] <- NF[t, s] # nonbreeder
-    #     c.expFY[t, s] <- N[1, t, s]
-    #     c.obsB[t, s] <- counts[3, t, s]
-    #     c.obsA[t, s] <- counts[2, t, s]
-    #     c.obsFY[t, s] <- counts[1, t, s]  # first year
-    #     c.repB[t, s] ~ dpois(NB[t,s] ) # simulated counts
-    #     c.repA[t, s] ~ dpois(NF[t,s] )
-    #     dssm.obsB[t, s] <- abs( ( (c.obsB[t, s]) - (c.expB[t, s]) ) / (c.obsB[t, s]+0.001)  )
-    #     dssm.obsA[t, s] <- abs( ( (c.obsA[t, s]) - (c.expA[t, s]) ) / (c.obsA[t, s]+0.001)  )
-    #     dssm.obsFY[t, s] <- abs( ( (c.obsFY[t, s]) - (c.expFY[t, s]) ) / (c.obsFY[t, s]+0.001)  )
-    #     dssm.repB[t, s] <- abs( ( (c.repB[t, s]) - (c.expB[t, s]) ) / (c.repB[t, s]+0.001) )
-    #     dssm.repA[t, s] <- abs( ( (c.repA[t, s]) - (c.expA[t, s]) ) / (c.repA[t, s]+0.001) )
-    #     dssm.repFY[t, s] <- abs( ( (c.repFY[t, s]) - (c.expFY[t, s]) ) / (c.repFY[t, s]+0.001) )
-    #   }} # t
-    # dmape.obs[1] <- sum(dssm.obsB[1:nyr, 1:nsite])
-    # dmape.obs[2] <- sum(dssm.obsA[1:nyr, 1:nsite])
-    # dmape.obs[3] <- sum(dssm.obsFY[1:nyr, 1:nsite])
-    # dmape.rep[1] <- sum(dssm.repB[1:nyr, 1:nsite])
-    # dmape.rep[2] <- sum(dssm.repA[1:nyr, 1:nsite])
-    # dmape.rep[3] <- sum(dssm.repFY[1:nyr, 1:nsite])
-    # tvm.obs[1] <- sd(c.obsB[1:nyr, 1:nsite])^2/mean(c.obsB[1:nyr, 1:nsite])
-    # tvm.obs[2] <- sd(c.obsA[1:nyr, 1:nsite])^2/mean(c.obsA[1:nyr, 1:nsite])
-    # tvm.obs[3] <- sd(c.obsFY[1:nyr, 1:nsite])^2/mean(c.obsFY[1:nyr, 1:nsite])
-    # tvm.rep[1] <- sd(c.repB[1:nyr, 1:nsite])^2/mean(c.repB[1:nyr, 1:nsite])
-    # tvm.rep[2] <- sd(c.repA[1:nyr, 1:nsite])^2/mean(c.repA[1:nyr, 1:nsite])
-    # tvm.rep[3] <- sd(c.repFY[1:nyr, 1:nsite])^2/mean(c.repFY[1:nyr, 1:nsite])
-    # # # Test statistic for number of turns
-    # for (s in 1:nsite){
-    #   for (t in 1:(nyr-2)){
-    #     tt1.obsB[t, s] <- step(c.obsB[t+2, s] - c.obsB[t+1, s])
-    #     tt2.obsB[t, s] <- step(c.obsB[t+1, s] - c.obsB[t, s])
-    #     tt3.obsB[t, s] <- equals(tt1.obsB[t, s] + tt2.obsB[t, s], 1)
-    #     tt1.obsA[t, s] <- step(c.obsA[t+2, s] - c.obsA[t+1, s])
-    #     tt2.obsA[t, s] <- step(c.obsA[t+1, s] - c.obsA[t, s])
-    #     tt3.obsA[t, s] <- equals(tt1.obsA[t, s] + tt2.obsA[t, s], 1)
-    #     tt1.obsFY[t, s] <- step(c.obsFY[t+2, s] - c.obsFY[t+1, s])
-    #     tt2.obsFY[t, s] <- step(c.obsFY[t+1, s] - c.obsFY[t, s])
-    #     tt3.obsFY[t, s] <- equals(tt1.obsFY[t, s] + tt2.obsFY[t, s], 1)
-    #   }} # t
-    # tturn.obs[1] <- sum(tt3.obsB[1:(nyr-2), 1:nsite])
-    # tturn.obs[2] <- sum(tt3.obsA[1:(nyr-2), 1:nsite])
-    # tturn.obs[3] <- sum(tt3.obsFY[1:(nyr-2), 1:nsite])
-    # for (s in 1:nsite){
-    #   for (t in 1:(nyr-2)){
-    #     tt1.repB[t, s] <- step(c.repB[t+2, s] - c.repB[t+1, s])
-    #     tt2.repB[t, s] <- step(c.repB[t+1, s] - c.repB[t, s])
-    #     tt3.repB[t, s] <- equals(tt1.repB[t, s] + tt2.repB[t, s], 1)
-    #     tt1.repA[t, s] <- step(c.repA[t+2, s] - c.repA[t+1, s])
-    #     tt3.repA[t, s] <- equals(tt1.repA[t, s] + tt2.repA[t, s], 1)
-    #     tt2.repA[t, s] <- step(c.repA[t+1, s] - c.repA[t, s])
-    #     tt1.repFY[t, s] <- step(c.repFY[t+2, s] - c.repFY[t+1, s])
-    #     tt2.repFY[t, s] <- step(c.repFY[t+1, s] - c.repFY[t, s])
-    #     tt3.repFY[t, s] <- equals(tt1.repFY[t, s] + tt2.repFY[t, s], 1)
-    #   }} # t
-    # tturn.rep[1] <- sum(tt3.repB[1:(nyr-2), 1:nsite])
-    # tturn.rep[2] <- sum(tt3.repA[1:(nyr-2), 1:nsite])
-    # tturn.rep[3] <- sum(tt3.repFY[1:(nyr-2), 1:nsite])
-    # 
+# Omitted. See IPM or PVA.
+
     ################################
     # Likelihood for survival
     ################################ 
-    # Calculate an averages for sites
-    # each year for integration
+    # Calculate yearly averages for sites for integration
     for (t in 1:nyr){
       for (s in 1:nsite){
         for (xxxx in 1:surv.end[t,s]){
-          # need to reorder because nimble doesn't 
+          # Need to reorder because nimble doesn't 
           # handle nonconsecutive indices
           # yrind.surv is a matrix of indices for each site
           phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
@@ -358,6 +247,9 @@ mycode <- nimbleCode(
     } #i
   } )
 
+#**********************
+#* Function to run model in NIMBLE
+#**********************
 run_ipm <- function(info, datl, constl, code){
   library('nimble')
   library('coda')
@@ -431,7 +323,9 @@ run_ipm <- function(info, datl, constl, code){
   return(post)
 } # run_ipm function end
 
-
+#*****************
+#* Run chains in parallel
+#*****************
 this_cluster <- makeCluster(10)
 post <- parLapply(cl = this_cluster, 
                   X = par_info,