model{
  # Pr(Y=0|x,z)

   for(i in 1:Nz){
    zeros[i] <- 0
    zeros[i] ~ dpois(mu[i])
    mu[i] <- -log(1-p0[i] + p0[i]*exp(-lambda1[i]))

    p0[i]<- p[i]    
    lambda1[i]<- lambda[i]*offset[i]
    logit(p[i]) <- b0.site[Site[i]] + b0.frt*Frt[i] + b0.pco*PCO[i] + b0.pca1*PCA1[i] +                 	                                                  
                       b0.pca2*PCA2[i] +b0.pca3*PCA3[i]+b0.lat*STX[i]+b0.long*STY[i]+eps0[Month[i], Sp[i]]*xi0
    log(lambda[i])<-b.site[Site[i]] + b.frt*Frt[i] + b.pco*PCO[i] + b.pca1*PCA1[i] + b.pca2*PCA2[i] +                                                   
                             b.pca3*PCA3[i]+b.lat*STX[i]+b.long*STY[i]+eps[Month[i], Sp[i]]*xi  
}
 
  # Pr(Y>0|x,z)

   for(i in (Nz+1):N){
    zeros[i] <- 0
    zeros[i] ~ dpois(mu[i])
    mu[i] <- -(log(p[i])-lambda1[i] + y[i]*log(lambda1[i])-logfact(y[i]))

    p1[i]<- p[i] 
    lambda1[i]<- lambda[i]*offset[i]
    logit(p[i]) <- b0.site[Site[i]]  + b0.frt*Frt[i] + b0.pco*PCO[i] + b0.pca1*PCA1[i] +                                                                      
                       b0.pca2*PCA2[i] +b0.pca3*PCA3[i]+b0.lat*STX[i]+b0.long*STY[i]+eps0[Month[i], Sp[i]]*xi0  
    log(lambda[i])<-b.site[Site[i]] + b.frt*Frt[i] + b.pco*PCO[i] + b.pca1*PCA1[i] + b.pca2*PCA2[i]+                                                 
                             b.pca3*PCA3[i]+b.lat*STX[i]+b.long*STY[i]+eps[Month[i], Sp[i]]*xi 
}
  

  # Vague Priors for model coefficients
   
   for(i in 1:n.cat){
    b0.site[i] ~ dnorm(0.0,0.001)
    b.site[i] ~ dnorm(0.0,0.001)
   }

    b.frt ~ dnorm(0.0,0.0001)
    b.pco ~ dnorm(0.0,0.0001)
    b.pca1 ~ dnorm(0.0,0.0001)
    b.pca2 ~ dnorm(0.0,0.0001)
    b.pca3 ~ dnorm(0.0,0.0001)
    b0.frt ~ dnorm(0.0,0.0001)
    b0.pco ~ dnorm(0.0,0.0001)
    b0.pca1 ~ dnorm(0.0,0.0001)
    b0.pca2 ~ dnorm(0.0,0.0001)
    b0.pca3 ~ dnorm(0.0,0.0001)
    b.lat ~ dnorm(0.0,0.0001)
    b.long ~ dnorm(0.0,0.0001)
    b0.lat ~ dnorm(0.0,0.0001)
    b0.long ~ dnorm(0.0,0.0001)

  # Month and species random effect

    for(m in 1:n.mo){
     for(s in 1:n.sp){
      eps0[m, s] ~ dnorm(0.0,tau0.mo)
      eps[m, s] ~ dnorm(0.0,tau.mo)
     }
  }

 # Setting up scale for half-cauchy prior on random effects 

   xi0~dnorm(0, tau0.xi)
   tau0.xi<-pow(prior.scale, -2)       #prior scale is sd of y (observations)
   xi~dnorm(0, tau.xi)
   tau.xi<-pow(prior.scale, -2) 

  # Vague Priors for precision

    tau0.mo ~ dgamma(0.5,0.5)       # chi^2 with 1 df (Gelman and Hill, pg. 430)
    tau.mo ~ dgamma(0.5,0.5)
    sigma0.mo<-abs(xi0)/sqrt(tau0.mo)      #sd of random effect
    sigma.mo<-abs(xi)/sqrt(tau.mo)           #sd of random effect
  
  # Prediction of parameters

    for(j in 1:n.cat){

     # probability
      post.logitp[j] <- b0.site[j] + b0.frt*avgfruit + b0.pco*avgpco + b0.pca1*avgpca1+                                                                 
                              b0.pca2*avgpca2 + b0.pca3*avgpca3 +b0.lat*avgstx+b0.long*avgsty
        post.p[j] <- 1/(1+exp(-post.logitp[j]))

     # mean
      post.loglam[j] <- b.site[j] + b.frt*avgfruit + b.pco*avgpco + b.pca1*avgpca1+                                                                     
                                 b.pca2*avgpca2 + b.pca3*avgpca3  +b.lat*avgstx+b.long*avgsty
       post.lambda[j] <- exp(post.loglam[j])

     # mean abundance
      post.abund[j]<- post.p[j]*post.lambda[j]    
    }
}
pardpardeftab720ri0qlqnatural

f1fs24 cf0 



}