Open Population Hierarchical Abundance Model with Removal Sampling

DM_Multinomial

# Model
cat("
    model{
    # Abundance Priors
    a.N ~ dnorm(0, 0.01)
    b.elev ~ dnorm(0, 0.01)
    #b.elev2 ~ dnorm(0, 0.01)
    #b.slope ~ dnorm(0, 0.01)
    b.drainage ~ dnorm(0, 0.01)
    b.hardwood ~ dnorm(0, 0.01)
    b.softwood ~ dnorm(0, 0.01)
    b.herb ~ dnorm(0, 0.01)
    b.cwd ~ dnorm(0, 0.01)
    b.litter ~ dnorm(0, 0.01)
    b.trap ~ dnorm(0, 0.01)
    b.stems ~ dnorm(0, 0.01)
    b.year2 ~ dnorm(0, 0.01)
    b.year3 ~ dnorm(0, 0.01)
    gam1 ~ dnorm(0, 0.01)
    
    # Random priors
    for(i in 1:nsites){
    b.site[i] ~ dnorm(0, tau.site)
    gam0[i] ~ dnorm(0, tau.gam0)
    }
    
    # Hyperpriors for random effects
    sigma.site ~ dunif(0, 10)
    tau.site <- 1/(sigma.site*sigma.site)
    
    sigma.gam0 ~ dunif(0, 5)
    tau.gam0 <- 1/(sigma.gam0*sigma.gam0)
    
    # Detection Priors
    p0 ~ dunif(-4, 4)
    #logitp0 <- log(p0/(1-p0))
    p.precip ~ dunif(-3, 3)
    p.trap ~ dunif(-3, 3)
    # p.day ~ dunif(-3, 3)

# First Closed Period
for(k in 1:1){
  for(i in 1:nplots){
    fac[i,1,k] <- 1
    p[i,1,k] <- p0 + p.precip*precip[i,1] + p.trap*traptype[i]
    mu[i,1,k] <- p[i,1,k]
    for(j in 2:8){
      # Detection for closed period 1
      p[i,j,k] <- p0 + p.precip*precip[i,j] + p.trap*traptype[i]
      fac[i,j,k] <- fac[i,j-1,k] * (1 - p[i,j-1,k])
      mu[i,j,k] <- p[i,j,k] * fac[i,j,k]
    }

      mu0[i,k] <- sum(mu[i,1:8,k])
      pi0[i,k] <- 1 - mu0[i,k]
      pcap[i,k] <- 1 - pi0[i,k] # redundant with mu0

    for(j in 1:8){
      muc[i,j,k] <- mu[i,j,k]/pcap[i,k]
    }
    
    
    # Multinomial Logit???
    mu0[i,k] <- sum(exp(mu[i,1:8,k]))
    for(j in 1:8){
      muc[i,j,k] <- exp(mu[i,j,k])/mu0[i,j,k]
    }
    
    # Abundance for closed period 1
    log(lambda[i])<- a.N + b.site[site[i]] + b.drainage*drainage[i] + b.softwood*softwood[i] + b.herb*herb[i] + b.cwd*cwd[i] + b.litter*litter[i] + b.stems*stems[i] + b.trap*traptype[i] + b.year2*year2[i] + b.year3*year3[i] + b.elev*elev[i] + b.hardwood*hardwood[i]#  + b.elev2*elev[i]^2  + b.slope*slope[i] 
    N[i,k] ~ dpois(lambda[i])
    ncap[i,k] ~ dbin(pcap[i,k],N[i,k])
    y[i,1:8] ~ dmulti(muc[i,1:8,k], ncap[i,k])
    
    # Open Period Change in Abundance
    gamma[i] <- min(exp(gam0[site[i]] + gam1*N[i,1]), 500)
  }
}

# Second Closed Period
for(k in 2:2){
  for(i in 1:nplots){
    fac[i,1,k] <- 1
    p[i,1,k] <- p0 + p.precip*precip[i,9] + p.trap*traptype[i]
    mu[i,1,k] <- p[i,1,k]
    for(j in 2:8){
      # Detection for closed period 1
      p[i,j,k] <- p0 + p.precip*precip[i,j+8] + p.trap*traptype[i]
      fac[i,j,k] <- fac[i,j-1,k] * (1 - p[i,j-1,k])
      mu[i,j,k] <- p[i,j,k] * fac[i,j,k]
    }

      mu0[i,k] <- sum(mu[i,1:8,k])
      pi0[i,k] <- 1 - mu0[i,k]
      pcap[i,k] <- 1 - pi0[i,k]

    for(j in 1:8){
      muc[i,j,k] <- mu[i,j,k]/pcap[i,k]
    }
    
    # Multinomial Logit???
    mu0[i,k] <- sum(exp(mu[i,1:8,k]))
    for(j in 1:8){
      muc[i,j,k] <- exp(mu[i,j,k])/mu0[i,j,k]
    }
    
    N[i,k] ~ dpois(gamma[i])
    ncap[i,k] ~ dbin(pcap[i,k],N[i,k])
    y[i,9:16] ~ dmulti(muc[i,1:8,k], ncap[i,k])
  }
}

    ### compute a bunch of fit-related stuff for Bayesian
    ### p-values.....
    
    for(i in 1:nplots){
      y.fit[i,1:8] ~ dmulti(muc[i,1:8,1],ncap[i,1])
      y.fit[i,9:16] ~ dmulti(muc[i,1:8,2],ncap[i,2])
      for(k in 1:2){
        ncap.fit[i,k] ~ dbin(pcap[i,k],N[i,k])
      }
        for(t in 1:8){
          e1[i,t]<- muc[i,t,1]*ncap[i,1]
          resid1[i,t]<- pow(pow(y[i,t],0.5)-pow(e1[i,t],0.5),2)
          resid1.fit[i,t]<- pow(pow(y.fit[i,t],0.5) - pow(e1[i,t],0.5),2)

          e1[i,t+8]<- muc[i,t,2]*ncap[i,2]
          resid1[i,t+8]<- pow(pow(y[i,t+8],0.5)-pow(e1[i,t+8],0.5),2)
          resid1.fit[i,t+8]<- pow(pow(y.fit[i,t+8],0.5) - pow(e1[i,t+8],0.5),2)
        }
      
      e2[i,1] <- N[i,1]*lambda[i]
      e2[i,2] <- N[i,2]*gamma[i]
      for(k in 1:2){
        resid2[i,k]<- pow(pow(ncap[i,k],0.5) - pow(e2[i,k],0.5),2)
        resid2.fit[i,k]<- pow(pow(ncap.fit[i,k],0.5) - pow(e2[i,k],0.5),2)
      }
    }
    ft1.data<- sum(resid1[,])
    ft1.post<- sum(resid1.fit[,])
    
    ft2.data<- sum(resid2[,])
    ft2.post<- sum(resid2.fit[,])
}
", fill=TRUE, file='DM3.txt')