###################################################################################################################
###############################################################################################################
#### Bird Nest Survival Modesls 
##############################################################################################################
###################################################################################################################

library(tidyverse)
library(lme4)
library(car)
library(emmeans)
library(dplyr)
library(MuMIn)
library(jagsUI)
library(loo)

# KF Bayes Model test
# Set working directory & load data
setwd("D:/Studies/semo lead study pop model/oa analysis")
pbjags<- read.csv(file = 'Pb_AllSPecies_NS.csv', header=T, strip.white=T) 

#get rid of one cardinal nest with extreme lead concentration > 12000 and all EABL
pbjags = filter(pbjags, soilpb < 12000, species_code != 'EABL')

#clean up data
nest_id<-(pbjags$nest_id)
survive <- (pbjags$survive)
h <-(pbjags$h)
species<-(pbjags$species)
site<-(pbjags$site)
year<-(pbjags$year)

#log and scale soilpb
lspb<-log10(pbjags$soilpb+0.01)
slspb<-as.numeric(scale(lspb))

mean(lspb)
sd(lspb) 
summary(pbjags$soilpb)


#notes on category variables
#Bixby=1 Ftown =2 Jins =3 WashKing =4 AHaw =5 
#2016=1;2017=2;2018=3;
#FISP=1; EATO=2; INBU=3; NOCA=4; 


##########################################
# 4 open cup species
# Nest survival ~ intercept + FE Spp Pb Spp*PB + RE year site; 
#########################################

# Bundle data
jags.data <- list(survive = survive, 
                  year = year, 
                  site = site,
                  nest_id = as.numeric(nest_id),
                  nnests=length(levels(pbjags$nest_id)),
                  h = h, 
                  species = species,  
                  slspb = slspb, 
                  N = nrow(pbjags))



sink("NestSurv4.1.jags")
cat("
    # likelihood
    model{
    for (i in 1:N){
    survive[i]~dbern(p[i])
    p[i]<-z[i]^h[i]
    logit(z[i])<- a + s[site[i]] + yr[year[i]] + sp[species[i]] + pb*slspb[i] + interx[species[i]]*slspb[i]
    
    loglik.survive[i] <- logdensity.bern(survive[i],p[i])
      
      # ############ Goodness of fit tests ###################################
      e.survive[i]<- z[i]^h[i]                           # original model prediction
      E.survive[i]<- pow((survive[i]- e.survive[i]),2)/(e.survive[i]+0.5)
      survive.sim[i] ~ dbern(p[i])              # create new realization of model
      E.survive.sim[i]<- pow((survive.sim[i]-e.survive[i]),2)/(e.survive[i]+0.5)
      
    }
    fit.survive<- sum(E.survive[])
    fit.survive.sim<- sum(E.survive.sim[])
    
    # corner conditions for categorical site, year, species
    s[1]<-0
    yr[1]<-0
    sp[1]<-0
    interx[1]<-0
    
    # Priors
    a ~ dnorm(0, 0.0001)                                        # intercept
    pb ~ dnorm(0, 0.00001)                                      # FE lead
    for (j  in 2:4){ sp[j] ~ dnorm(0, 0.0001)}                  # FE species
    for (n in 2:4){ interx[n] ~ dnorm(0, 0.0001)}               # FE interaction
    for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
    for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
    tau~dgamma(0.01, 0.01)                                      # Survival precision
    mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
    tau.site~dgamma(0.01, 0.01)
    mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
    tau.year~dgamma(0.01, 0.01)

   #Derived these are maginal effects of PB on each species
    for (j  in 1:4){
    pb.sp[j]<-pb+interx[j]
    }

    }
    ",fill = TRUE)
sink()
dim


# Initial values
#inits <- function (){list(a=1, sp=c(NA,0,0,0), interx=c(NA,0,0,0), pb=0, yr=c(NA,0,0), s=c(NA, 0, 0, 0, 0), mu.site=0, mu.year=0, tau.site=1, tau.year=1)}
inits <- function (){list()}
# Parameters monitored
parameters <- c("sp", "pb", "interx", "mu.year", "mu.site", "tau", "a", "pb.sp","loglik.survive","fit.survive","fit.survive.sim")


# MCMC settings
ni <- 20000
nt <- 5
nb <- 5000
nc <- 3
# Call jags from R
NestSurv4.1 <- jags(jags.data, inits, parameters, "NestSurv4.1.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb,parallel = T)
print(NestSurv4.1)

save(NestSurv4.1, file='Pb_OCNests_NestSurv.RData')

# Check Bayes P-value
bpv.n <-mean(NestSurv4.1$sims.list$fit.survive.sim > NestSurv4.1$sims.list$fit.survive)
bpv.n
plot(NestSurv4.1$sims.list$fit.survive,NestSurv4.1$sims.list$fit.survive.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(NestSurv4.1$sims.list$fit.survive/NestSurv4.1$sims.list$fit.survive.sim )
c.hat


# Calculate WAIC & LOO Stats
loglik <- NestSurv4.1$sims.list$loglik.survive
loo.1<-loo(loglik, cores = 3)
print(loo.1)


##########################################
#Null model
# Nest survival ~ intercept + FE Spp RE year site; 
#########################################

# Bundle data
jags.data <- list(survive = survive, 
                  year = year, 
                  site = site,
                  nest_id = as.numeric(nest_id),
                  nnests=length(levels(pbjags$nest_id)),
                  h = h, 
                  species = species,  
                  slspb = slspb, 
                  N = nrow(pbjags))

sink("NestSurv4.0.jags")
cat("
    # likelihood
    model{
    for (i in 1:N){
    survive[i]~dbern(p[i])
    p[i]<-z[i]^h[i]
    logit(z[i])<- a + s[site[i]] + yr[year[i]] + sp[species[i]] 
    
      loglik.survive[i] <- logdensity.bern(survive[i],p[i])
    
      # ############ Goodness of fit tests ###################################
      e.survive[i]<- z[i]^h[i]                           # original model prediction
      E.survive[i]<- pow((survive[i]- e.survive[i]),2)/(e.survive[i]+0.5)
      survive.sim[i] ~ dbern(p[i])              # create new realization of model
      E.survive.sim[i]<- pow((survive.sim[i]-e.survive[i]),2)/(e.survive[i]+0.5)
      
    }
    
    fit.survive<- sum(E.survive[])
    fit.survive.sim<- sum(E.survive.sim[])
    
    # corner conditions for categorical site, year, species
    s[1]<-0
    yr[1]<-0
    sp[1]<-0
    
    # Priors
    a ~ dnorm(0, 0.0001)                                        # intercept
    
    for (j  in 2:4){ sp[j] ~ dnorm(0, 0.0001)}                  # FE species
    
    for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
    for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
    tau~dgamma(0.01, 0.01)                                      # Survival precision
    mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
    tau.site~dgamma(0.01, 0.01)
    mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
    tau.year~dgamma(0.01, 0.01)
    }
    ",fill = TRUE)
sink()
dim

# Initial values
inits <- function (){list(a=1, sp=c(NA,0,0,0), yr=c(NA,0,0), s=c(NA, 0, 0, 0, 0), mu.site=0, mu.year=0, tau.site=1, tau.year=1)}

# Parameters monitored
parameters <- c("sp", "mu.year", "mu.site", "tau", "a","loglik.survive","fit.survive","fit.survive.sim")

# MCMC settings
ni <- 20000
nt <- 5
nb <- 5000
nc <- 3
# Call jags from R
NestSurv4.0 <- jags(jags.data, inits, parameters, "NestSurv4.0.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb,parallel = T)
print(NestSurv4.0)

# Check Bayes P-value
bpv.n <-mean(NestSurv4.0$sims.list$fit.survive.sim > NestSurv4.0$sims.list$fit.survive)
bpv.n
plot(NestSurv4.0$sims.list$fit.survive,NestSurv4.0$sims.list$fit.survive.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(NestSurv4.0$sims.list$fit.survive/NestSurv4.0$sims.list$fit.survive.sim )
c.hat


# Calculate WAIC & LOO Stats
loglik <- NestSurv4.0$sims.list$loglik.survive
loo.0<-loo(loglik, cores = 3)
print(loo.0)

Model.compare<-loo_compare(loo.1,loo.0)
View(Model.compare)

###################################################################################################################
###### EABL Nest survival models 
################################################################################################################

library(tidyverse)
library(lme4)
library(car)
library(emmeans)
library(dplyr)
library(MuMIn)
library(jagsUI)
library(loo)

# KF Bayes Model test
# Set working directory & load data
setwd("D:/Studies/semo lead study pop model/oa analysis")
pbjags<- read.csv(file = 'Pb_AllSPecies_NS.csv', header=T, strip.white=T) 

#keep eabl observations
pbjags = filter(pbjags, species_code == 'EABL')

#clean up data
nest_id<-(pbjags$nest_id)
survive <- (pbjags$survive)
h <-(pbjags$h)
species<-(pbjags$species)
site<-(pbjags$site)
year<-(pbjags$year)

#log and scale soilpb
lspb<-log10(pbjags$soilpb+0.01)
slspb<-as.numeric(scale(lspb))

mean(lspb)
sd(lspb) 
summary(pbjags$soilpb)


#notes on category variables
#Bixby=1 Ftown =2 Jins =3 WashKing =4 AHaw =5 
#2016=1;2017=2;2018=3;
#FISP=1; EATO=2; INBU=3; NOCA=4; 


##########################################
# Nest survival ~ intercept + FE Spp Pb Spp*PB + RE year site; 
#########################################

# Bundle data
jags.data <- list(survive = survive, 
                  year = year, 
                  site = site,
                  nest_id = as.numeric(nest_id),
                  nnests=length(levels(pbjags$nest_id)),
                  h = h, 
                  slspb = slspb, 
                  N = nrow(pbjags))



sink("NestSurvEB.1.jags")
cat("
    # likelihood
    model{
    for (i in 1:N){
    survive[i]~dbern(p[i])
    p[i]<-z[i]^h[i]
    logit(z[i])<- a + s[site[i]] + yr[year[i]] + pb*slspb[i] 
    
    loglik.survive[i] <- logdensity.bern(survive[i],p[i])
      
      # ############ Goodness of fit tests ###################################
      e.survive[i]<- z[i]^h[i]                           # original model prediction
      E.survive[i]<- pow((survive[i]- e.survive[i]),2)/(e.survive[i]+0.5)
      survive.sim[i] ~ dbern(p[i])              # create new realization of model
      E.survive.sim[i]<- pow((survive.sim[i]-e.survive[i]),2)/(e.survive[i]+0.5)
      
    }
    fit.survive<- sum(E.survive[])
    fit.survive.sim<- sum(E.survive.sim[])
    
    # corner conditions for categorical site, year, species
    s[1]<-0
    yr[1]<-0
    
    # Priors
    a ~ dnorm(0, 0.0001)                                        # intercept
    pb ~ dnorm(0, 0.00001)                                      # FE lead
    for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
    for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
    tau~dgamma(0.01, 0.01)                                      # Survival precision
    mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
    tau.site~dgamma(0.01, 0.01)
    mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
    tau.year~dgamma(0.01, 0.01)

    }
    ",fill = TRUE)
sink()
dim


# Initial values
inits <- function (){list()}
# Parameters monitored
parameters <- c("pb", "mu.year", "mu.site", "tau", "a", "loglik.survive","fit.survive","fit.survive.sim")


# MCMC settings
ni <- 20000
nt <- 5
nb <- 5000
nc <- 3
# Call jags from R
NestSurvEB.1 <- jags(jags.data, inits, parameters, "NestSurvEB.1.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb,parallel = T)
print(NestSurvEB.1)

save(NestSurvEB.1, file='Pb_EABL_NestSurv.RData')

# Check Bayes P-value
bpv.n <-mean(NestSurvEB.1$sims.list$fit.survive.sim > NestSurvEB.1$sims.list$fit.survive)
bpv.n
plot(NestSurvEB.1$sims.list$fit.survive,NestSurvEB.1$sims.list$fit.survive.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(NestSurvEB.1$sims.list$fit.survive/NestSurvEB.1$sims.list$fit.survive.sim )
c.hat


# Calculate WAIC & LOO Stats
loglik <- NestSurvEB.1$sims.list$loglik.survive
loo.1<-loo(loglik, cores = 3)
print(loo.1)


##########################################
#Null model
# Nest survival ~ intercept + FE Spp RE year site; 
#########################################

# Bundle data
jags.data <- list(survive = survive, 
                  year = year, 
                  site = site,
                  nest_id = as.numeric(nest_id),
                  nnests=length(levels(pbjags$nest_id)),
                  h = h, 
                  slspb = slspb, 
                  N = nrow(pbjags))

sink("NestSurvEB.0.jags")
cat("
    # likelihood
    model{
    for (i in 1:N){
    survive[i]~dbern(p[i])
    p[i]<-z[i]^h[i]
    logit(z[i])<- a + s[site[i]] + yr[year[i]]  
    
      loglik.survive[i] <- logdensity.bern(survive[i],p[i])
    
      # ############ Goodness of fit tests ###################################
      e.survive[i]<- z[i]^h[i]                           # original model prediction
      E.survive[i]<- pow((survive[i]- e.survive[i]),2)/(e.survive[i]+0.5)
      survive.sim[i] ~ dbern(p[i])              # create new realization of model
      E.survive.sim[i]<- pow((survive.sim[i]-e.survive[i]),2)/(e.survive[i]+0.5)
      
    }
    
    fit.survive<- sum(E.survive[])
    fit.survive.sim<- sum(E.survive.sim[])
    
    # corner conditions for categorical site, year, species
    s[1]<-0
    yr[1]<-0
    
    # Priors
    a ~ dnorm(0, 0.0001)                                        # intercept
    
    for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
    for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
    tau~dgamma(0.01, 0.01)                                      # Survival precision
    mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
    tau.site~dgamma(0.01, 0.01)
    mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
    tau.year~dgamma(0.01, 0.01)
    }
    ",fill = TRUE)
sink()
dim

# Initial values
inits <- function (){list(a=1, yr=c(NA,0,0), s=c(NA, 0, 0, 0, 0), mu.site=0, mu.year=0, tau.site=1, tau.year=1)}

# Parameters monitored
parameters <- c("mu.year", "mu.site", "tau", "a","loglik.survive","fit.survive","fit.survive.sim")

# MCMC settings
ni <- 20000
nt <- 5
nb <- 5000
nc <- 3
# Call jags from R
NestSurvEB.0 <- jags(jags.data, inits, parameters, "NestSurvEB.0.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb,parallel = T)
print(NestSurvEB.0)

# Check Bayes P-value
bpv.n <-mean(NestSurvEB.0$sims.list$fit.survive.sim > NestSurvEB.0$sims.list$fit.survive)
bpv.n
plot(NestSurvEB.0$sims.list$fit.survive,NestSurvEB.0$sims.list$fit.survive.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(NestSurvEB.0$sims.list$fit.survive/NestSurvEB.0$sims.list$fit.survive.sim )
c.hat


# Calculate WAIC & LOO Stats
loglik <- NestSurvEB.0$sims.list$loglik.survive
loo.0<-loo(loglik, cores = 3)
print(loo.0)

Model.compare<-loo_compare(loo.1,loo.0)
View(Model.compare)


###################################################################################################################
###################################################################################################################
###### Bird Fledging Models
################################################################################################################
###################################################################################################################


library(tidyverse)
library(jagsUI)
library(loo)
library(ggplot2)
library(scales)


# Set working directory & load data
# setwd("D:/Studies/semo lead study pop model/oa analysis")

pbjags <- read.csv(file = 'DATA/Pb_OCNests_Repro.csv', header=T, strip.white=T)

#get rid of one cardinal nest with extreme lead concentration > 12000
pbjags = filter(pbjags, soilpb < 12000)

pbjags$lspb<-log10(pbjags$soilpb+0.01)
pbjags$slspb<-as.numeric(scale(pbjags$lspb))

### species 1=FISP 2=EATO 3=INBU4=NOCA
### site 1=Bixby 2=Ftown 3=Jins 4=WashKing 5=Hawn
 


##### Model 1 #############################################################
# Fledging Model
# glmm with gamma dist and log link
# fledged = intercept species lsoilpb species*lsoilpb;  random intercepts: site year
#################################################################

pbjagsf <- subset(pbjags, (fledged != "NA") & (fledged != "0"))
#results in sample of 136 nests
#mean(pbjags$lspb)
#sd(pbjags$lspb) 
#min(pbjags$slspb)
#max(pbjags$slspb) 

sink("NestFledge1.jags")
cat("
    # likelihood
    model{
      for (i in 1:N){
        fledged[i]  ~dgamma(shape[i],rate[i])
        shape[i] <- pow(f[i],2) / pow(sd,2)
        rate[i] <- f[i] / pow(sd,2)
        log(f[i])<- a + s[site[i]] + yr[year[i]] + sp[species[i]] + pb*lsoilpb[i] + interx[species[i]]*lsoilpb[i]
      
        # ############ Goodness of fit tests ###################################
        e.fledged[i]<- f[i]                           # original model prediction
        E.fledged[i]<- pow((fledged[i]- e.fledged[i]),2)/(e.fledged[i]+0.5)
        fledged.sim[i] ~ dgamma(shape[i],rate[i])              # create new realization of model
        E.fledged.sim[i]<- pow((fledged.sim[i]-e.fledged[i]),2)/(e.fledged[i]+0.5)
        
        loglik.fledged[i] <- logdensity.gamma(fledged[i],shape[i],rate[i])
  
      }

      fit.fledged <- sum(E.fledged[])
      fit.fledged.sim<- sum(E.fledged.sim[])
      
      sp[1]<-0
      interx[1]<-0
      s[1]<-0
      yr[1]<-0
     
      # Priors
      sd ~ dunif(0,10)
      a ~ dnorm(0, 0.0001)                                        # intercept
      pb ~ dnorm(0, 0.00001)                                      # FE lead
      for (j  in 2:4){ sp[j] ~ dnorm(0, 0.0001)}                  # FE species
      for (n in 2:4){ interx[n] ~ dnorm(0, 0.0001)}               # FE interaction
      for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
      for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
      tau~dgamma(0.01, 0.01)                                      # Survival precision
      mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
      tau.site~dgamma(0.01, 0.01)
      mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
      tau.year~dgamma(0.01, 0.01)
  
      #Derived these are maginal effects of PB on each species
      for (j  in 1:4){
        pb.sp[j]<-pb+interx[j]
      }
    }
    ",fill = TRUE)
sink()
dim

# Bundle data
fledge1.data <- list(fledged = (pbjagsf$fledged), year = (pbjagsf$YEAR), site = (pbjagsf$site), species = (pbjagsf$species),  lsoilpb = (pbjagsf$slspb), N = nrow(pbjagsf))

# Initial values
inits <- function (){list(a=1, sp=c(NA,0,0,0), interx=c(NA,0,0,0), pb=0, yr=c(NA,0,0), s=c(NA, 0, 0, 0, 0), mu.site=0, mu.year=0, tau.site=1, tau.year=1)}

# Parameters monitored
parameters <- c("sp", "pb", "interx", "mu.year", "mu.site", "tau", "a", "pb.sp","r","loglik.fledged","fit.fledged","fit.fledged.sim","fledged.sim")

# MCMC settings
ni <- 20000; nt <- 5; nb <- 5000; nc <- 3

# Call jags from R
NestFledge1 <- jags(fledge1.data, inits, parameters, "NestFledge1.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, parallel = T,codaOnly = "fledged.sim")

print(NestFledge1)
# plot(NestFledge)
save(NestFledge1, file='Pb_OCNests_Fledge.RData')

# Check Bayes P-value
bpv.n <-mean(NestFledge1$sims.list$fit.fledged.sim > NestFledge1$sims.list$fit.fledged)
bpv.n
plot(NestFledge1$sims.list$fit.fledged,NestFledge1$sims.list$fit.fledged.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(NestFledge1$sims.list$fit.fledged/NestFledge1$sims.list$fit.fledged.sim )
c.hat

#Plot data and simulated values
fit_df <- data.frame(y = c(c(fledge1.data$fledged), c(NestFledge1$sims.list$fledged.sim)),
                     data = c(rep("Observed", length(fledge1.data$fledged)),rep("Simulated",length(NestFledge1$sims.list$fledged.sim))))

# ggplot(fit_df, aes(x = y, fill = data)) + geom_histogram(position = "dodge")
ggplot(fit_df, aes(x=y,fill=data)) +  
  geom_histogram(aes(y=0.5*..density..),alpha=0.5,position='identity',binwidth=0.5)


# Calculate WAIC & LOO Stats
loglik <- NestFledge1$sims.list$loglik.fledged
loo.1<-loo(loglik, cores = 3)
print(loo.1)


##### Model 1 Null #############################################################
# Fledging Model
# glmm with poisson dist and log link
# fledged = intercept species;  random intercepts: site year
#################################################################


sink("NestFledge0.jags")
cat("
    # likelihood
    model{
      for (i in 1:N){
        fledged[i]  ~dgamma(shape[i],rate[i])
        shape[i] <- pow(f[i],2) / pow(sd,2)
        rate[i] <- f[i] / pow(sd,2)
        log(f[i]) <- a + s[site[i]] + yr[year[i]] + sp[species[i]]  
   
        # ############ Goodness of fit tests ###################################
        e.fledged[i]<- f[i]                           # original model prediction
        E.fledged[i]<- pow((fledged[i]- e.fledged[i]),2)/(e.fledged[i]+0.5)
        fledged.sim[i] ~ dgamma(shape[i],rate[i])              # create new realization of model
        E.fledged.sim[i]<- pow((fledged.sim[i]-e.fledged[i]),2)/(e.fledged[i]+0.5)
        
        loglik.fledged[i] <- logdensity.gamma(fledged[i],shape[i],rate[i])
  
      }
      
      fit.fledged <- sum(E.fledged[])
      fit.fledged.sim<- sum(E.fledged.sim[])
        
      sp[1]<-0
      s[1]<-0
      yr[1]<-0
      
      # Priors
      sd ~ dunif(0,10)
      a ~ dnorm(0, 0.0001)                                        # intercept
      for (j  in 2:4){ sp[j] ~ dnorm(0, 0.0001)}                  # FE species
      for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
      for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
      tau~dgamma(0.01, 0.01)                                      # Survival precision
      mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
      tau.site~dgamma(0.01, 0.01)
      mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
      tau.year~dgamma(0.01, 0.01)
    }
    ",fill = TRUE)
sink()
dim

# Bundle data
fledge0.data <- list(fledged = (pbjagsf$fledged), year = (pbjagsf$YEAR), site = (pbjagsf$site), species = (pbjagsf$species), N = nrow(pbjagsf))

# Initial values
inits <- function (){list(a=1, sp=c(NA,0,0,0), yr=c(NA,0,0), s=c(NA, 0, 0, 0, 0), mu.site=0, mu.year=0, tau.site=1, tau.year=1)}

# Parameters monitored
parameters <- c("sp", "mu.year", "mu.site", "tau", "a","r","loglik.fledged","fit.fledged","fit.fledged.sim","fledged.sim")

# MCMC settings
ni <- 20000; nt <- 5; nb <- 5000; nc <- 3

# Call jags from R
NestFledge0 <- jags(fledge0.data, inits, parameters, "NestFledge0.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, parallel = T,codaOnly = "fledged.sim")
print(NestFledge0)

# Check Bayes P-value
bpv.n <-mean(NestFledge0$sims.list$fit.fledged.sim > NestFledge0$sims.list$fit.fledged)
bpv.n
plot(NestFledge0$sims.list$fit.fledged,NestFledge0$sims.list$fit.fledged.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(NestFledge0$sims.list$fit.fledged/NestFledge0$sims.list$fit.fledged.sim )
c.hat

#Plot data and simulated values
fit_df <- data.frame(y = c(c(fledge1.data$fledged), c(NestFledge1$sims.list$fledged.sim)),
                     data = c(rep("Observed", length(fledge1.data$fledged)),rep("Simulated",length(NestFledge1$sims.list$fledged.sim))))

# ggplot(fit_df, aes(x = y, fill = data)) + geom_histogram(position = "dodge")
ggplot(fit_df, aes(x=y,fill=data)) +  
  geom_histogram(aes(y=0.5*..density..),alpha=0.5,position='identity',binwidth=0.5)

# Calculate WAIC & LOO Stats
loglik <- NestFledge0$sims.list$loglik.fledged
loo.0<-loo(loglik, cores = 3)
print(loo.0)



####################################################################################################################################
#### EABL Fledging model
###################################################################################################################################
library(tidyverse)
library(jagsUI)
library(loo)
library(ggplot2)
library(scales)

# Set working directory & load data
setwd("D:/Studies/semo lead study pop model/oa analysis")
eabl <- read.csv(file = "Pb_EABL_Repro.csv")
### note this uses updated soil PB from 2019

## Scale soil and blood Pb using updated soil pb 
eabl$slspb<-as.numeric(scale(eabl$lsoilpb))

#str(eabl)
#summary(eabl)

##### Model 2 #############################################################
# fledged = intercept scaled lsoil Pb;  random intercepts: site year
# glmm with gamma dist and log link
#################################################################

eabls <- subset(eabl, (slspb != "NA") & (fledged != "0"))

sink("eablfledge2.jags")
cat("
    # likelihood
    model{
      for (i in 1:N){
        fledged[i]  ~dgamma(shape[i],rate[i])
        shape[i] <- pow(f[i],2) / pow(sd,2)
        rate[i] <- f[i] / pow(sd,2)
    log(f[i])<- a + s[site[i]] + yr[year[i]] + pb*slspb[i]
    
        # ############ Goodness of fit tests ###################################
        e.fledged[i]<- f[i]                           # original model prediction
        E.fledged[i]<- pow((fledged[i]- e.fledged[i]),2)/(e.fledged[i]+0.5)
        fledged.sim[i] ~ dgamma(shape[i],rate[i])              # create new realization of model
        E.fledged.sim[i]<- pow((fledged.sim[i]-e.fledged[i]),2)/(e.fledged[i]+0.5)
        
        loglik.fledged[i] <- logdensity.gamma(fledged[i],shape[i],rate[i])
  
      }

      fit.fledged <- sum(E.fledged[])
      fit.fledged.sim<- sum(E.fledged.sim[])
    
    s[1]<-0
    yr[1]<-0
    sd ~ dunif(0,10)
    # Priors
    a ~ dnorm(0, 0.0001)                                        # intercept
    pb ~ dnorm(0, 0.00001)                                      # FE lead
    for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
    for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
    tau~dgamma(0.01, 0.01)                                      # Survival precision
    mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
    tau.site~dgamma(0.01, 0.01)
    mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
    tau.year~dgamma(0.01, 0.01)
    }
    ",fill = TRUE)
sink()
dim

# Bundle data
fledge.data <- list(fledged = (eabls$fledged), year = (eabls$year), site = (eabls$site), slspb = (eabls$slspb), N = nrow(eabls))

# Initial values
inits <- function (){list(a=1, pb=0, yr=c(NA,0,0), mu.site=0, mu.year=0, tau.site=1, tau.year=1)}

# Parameters monitored
parameters <- c("pb", "mu.year", "mu.site", "tau", "a", "r","loglik.fledged","fit.fledged","fit.fledged.sim","fledged.sim")

# MCMC settings
ni <- 20000
nt <- 5
nb <- 5000
nc <- 3

# Call jags from R
eablfledge2 <- jags(fledge.data, inits, parameters, "eablfledge2.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, parallel = T)

print(eablfledge2)
save(eablfledge2, file='Pb_EABLE_Fledge.RData')

# Check Bayes P-value
bpv.n <-mean(eablfledge2$sims.list$fit.fledged.sim > eablfledge2$sims.list$fit.fledged)
bpv.n
plot(eablfledge2$sims.list$fit.fledged,eablfledge2$sims.list$fit.fledged.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(eablfledge2$sims.list$fit.fledged/eablfledge2$sims.list$fit.fledged.sim )
c.hat

#Plot data and simulated values
fit_df <- data.frame(y = c(c(fledge.data$fledged), c(eablfledge2$sims.list$fledged.sim)),
                     data = c(rep("Observed", length(fledge.data$fledged)),rep("Simulated",length(eablfledge2$sims.list$fledged.sim))))

# ggplot(fit_df, aes(x = y, fill = data)) + geom_histogram(position = "dodge")
ggplot(fit_df, aes(x=y,fill=data)) +  
  geom_histogram(aes(y=0.5*..density..),alpha=0.5,position='identity',binwidth=0.5)


# Calculate WAIC & LOO Stats
loglik <- eablfledge2$sims.list$loglik.fledged
loo.1<-loo(loglik, cores = 3)
print(loo.1)



##### Model 2 NULL  #############################################################
# fledged = intercept ;  random intercepts: site year
# glmm with poisson dist and log link
#################################################################

sink("eablfledge2n.jags")
cat("
    # likelihood
    model{
      for (i in 1:N){
        fledged[i]  ~dgamma(shape[i],rate[i])
        shape[i] <- pow(f[i],2) / pow(sd,2)
        rate[i] <- f[i] / pow(sd,2)
    log(f[i])<- a + s[site[i]] + yr[year[i]] 
    
        # ############ Goodness of fit tests ###################################
        e.fledged[i]<- f[i]                           # original model prediction
        E.fledged[i]<- pow((fledged[i]- e.fledged[i]),2)/(e.fledged[i]+0.5)
        fledged.sim[i] ~ dgamma(shape[i],rate[i])              # create new realization of model
        E.fledged.sim[i]<- pow((fledged.sim[i]-e.fledged[i]),2)/(e.fledged[i]+0.5)
        
        loglik.fledged[i] <- logdensity.gamma(fledged[i],shape[i],rate[i])
  
      }

      fit.fledged <- sum(E.fledged[])
      fit.fledged.sim<- sum(E.fledged.sim[])
    
    s[1]<-0
    yr[1]<-0
    sd ~ dunif(0,10)    
    # Priors
    a ~ dnorm(0, 0.0001)                                        # intercept
    #pb ~ dnorm(0, 0.00001)                                      # FE lead
    for (k in 2:5){ s[k] ~ dnorm(mu.site, tau.site)}            # RE site
    for (m in 2:3){ yr[m] ~ dnorm(mu.year, tau.year)}           # RE year
    tau~dgamma(0.01, 0.01)                                      # Survival precision
    mu.site~dnorm(0, 0.00001)                                   # hyperpriors site
    tau.site~dgamma(0.01, 0.01)
    mu.year~dnorm(0, 0.00001)                                   # hyperpriors year
    tau.year~dgamma(0.01, 0.01)
    }
    ",fill = TRUE)
sink()
dim

# Bundle data
fledge.data <- list(fledged = (eabls$fledged), year = (eabls$year), site = (eabls$site), N = nrow(eabls))

# Initial values
inits <- function (){list(a=1, yr=c(NA,0,0), mu.site=0, mu.year=0, tau.site=1, tau.year=1)}

# Parameters monitored
parameters <- c("mu.year", "mu.site", "tau", "a", "r","loglik.fledged","fit.fledged","fit.fledged.sim","fledged.sim")

# MCMC settings
ni <- 100000
nt <- 5
nb <- 10000
nc <- 3

# Call jags from R
eablfledge2n <- jags(fledge.data, inits, parameters, "eablfledge2n.jags", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, parallel = T)

print(eablfledge2n)

# Check Bayes P-value
bpv.n <-mean(eablfledge2n$sims.list$fit.fledged.sim > eablfledge2n$sims.list$fit.fledged)
bpv.n
plot(eablfledge2n$sims.list$fit.fledged,eablfledge2n$sims.list$fit.fledged.sim,main=(paste0("Bayesian P-value: ",round(bpv.n,3))))
abline(0,1,col="red")

c.hat <- mean(eablfledge2n$sims.list$fit.fledged/eablfledge2n$sims.list$fit.fledged.sim )
c.hat

#Plot data and simulated values
fit_df <- data.frame(y = c(c(fledge.data$fledged), c(eablfledge2n$sims.list$fledged.sim)),
                     data = c(rep("Observed", length(fledge.data$fledged)),rep("Simulated",length(eablfledge2n$sims.list$fledged.sim))))

# ggplot(fit_df, aes(x = y, fill = data)) + geom_histogram(position = "dodge")
ggplot(fit_df, aes(x=y,fill=data)) +  
  geom_histogram(aes(y=0.5*..density..),alpha=0.5,position='identity',binwidth=0.5)


# Calculate WAIC & LOO Stats
loglik <- eablfledge2$sims.list$loglik.fledged
loo.1<-loo(loglik, cores = 3)
print(loo.1)








###################################################################################################################
#################################################################################################################################
###### Code to fit the IBM model across bird species 
#############################################################################################################################
###################################################################################################################


library(jagsUI)
library(ggplot2)
library(tidyr)
library(see)
library(ggridges)
#notes on category variables
#Bixby=1 Ftown =2 Jins =3 WashKing =4 AHaw =5 
#2016=1;2017=2;2018=3;
#FISP=1; EATO=2; INBU=3; NOCA=4; 

setwd("C:/Users/bonnott/OneDrive - University of Missouri/BONNOT/RESEARCH/USFWS/Pb")

# Set species list and species for looping and pulling posterior data
species.list <- c("FISP","EATO","INBU","NOCA","EABL")
PARMS <- read.csv("DATA/PB_SPECIES_NESTING_PARMS_Final.csv")

#Read data from Frank's initial models
pbjags<- read.csv(file = 'DATA/Pb_AllSpecies_NS.csv', header=T, strip.white=T) 

#get rid of one cardinal nest with extreme lead concentration > 12000
pbjags2 <- subset(pbjags, soilpb < 12000)
lspb<-log10(pbjags$soilpb+0.01)
slspb<-as.numeric(scale(lspb))
pb.mean <- mean(lspb)
pb.sd <- sd(lspb) 

#Load results from Frank's initial models
NS.post<-load("DATA/MANUSCRIPT/Pb_OCNests_NestSurv.RData")
FL.post<-load("DATA/MANUSCRIPT/Pb_OCNests_Fledge.RData")

iterations<-length(FL.post$sims.list$deviance)  #identify size of posterior sample to create simulations
slspb.sim<-c(25,50,100,200,400,800,1600,3200,4000)        #create the soil Pb levels to predict across

reps <- 100

for (species in 1:4){
  
  Prod <- array(0,dim=c(reps,length(slspb.sim),iterations))         #Initiate the annual productivity array that will hold annual productivity values
  Success<-array(0,dim=c(reps,length(slspb.sim),2))
  
  for (rep in 1:reps){
    
    #Start annual loop
    for (s in 1:length(slspb.sim)){
      
      pb.sim <- (log(slspb.sim[s]+0.1)-pb.mean)/pb.sd
      
      NS <- array(NA,dim=c(iterations,PARMS[1,species+1]))                   #Create an array to track Nest Status (NS) through breeding season
      NS[,1:3] <- 0                                                          #Set first three days of season to no activity (0) for nest building
      NS[,4]<-1                                                              #Start nesting activity on 4th day of season (mask filters out cells not within region)
      
      Prod.day<-matrix(0,iterations,PARMS[1,species+1])                        #Reset productivity and nest completions to zero each year
      Complete<-numeric(iterations)
      
      #Estimate logit of global-interactive nest survival model
      mu.NS <-  NS.post$sims.list$a + NS.post$sims.list$mu.site + NS.post$sims.list$mu.year + NS.post$sims.list$sp[,species] + NS.post$sims.list$pb*pb.sim + NS.post$sims.list$interx[,species]*pb.sim
      ds<-exp(mu.NS)/(1+exp(mu.NS))    #Transform NS_global from the logit scale
      mu.FL <- pmin(exp(FL.post$sims.list$a + FL.post$sims.list$mu.site + FL.post$sims.list$mu.year + FL.post$sims.list$sp[,species] + FL.post$sims.list$pb*pb.sim + FL.post$sims.list$interx[,species]*pb.sim),PARMS[13,species+1])
      
      # print(paste0("Mean =",mean(ds)))
      # hist(ds,main=paste0(slspb.sim[s]))
      # hist(mu.FL,main=paste0(slspb.sim[s]))
      
      
      #Start season loop; Iterate through  days from May 15th to August 30
      for(d in 5:PARMS[1,species+1]){
        
        #Draw nest fate and brood size replicate for each day
        FL <- pmin(rpois(iterations,mu.FL),PARMS[13,species+1])
        fate<-matrix(rbinom(iterations,1,ds))   #Derive the fate of cells based on daily survival with random draws from binomial dist. 
        
        #Calculate the status of each cell's nest for each day based on if/thens
        NS[,d] <- ifelse(NS[,d-1]==-9,                                                                                  #Is inactive?
                         -9,                                                                                                      #Yes -> Stay inactive
                         ifelse(NS[,d-1]>0,                                                                                       #No (active) -> Is it currently nesting?
                                ifelse(NS[,d-1]==PARMS[12,species+1],                                                                    #Yes -> Has nest survived for full length of days (considered successful)
                                       ifelse(Complete==2, -9, -2),                                                                              #Yes -> Is it the 2nd brood? Yes-> go inactive for season; No ->go into post-fledging period prior to 2nd attempt
                                       NS[,d-1]*fate+fate),                                                                                      #No -> Simulate fate of active nest for that day (e.g., a nest that survives day 5 will go to day 6, else it will go to 0)
                                ifelse(NS[,d-1]==-2,                                                                                     #No -> Is it in post-fledging status?
                                       ifelse(apply(NS[,max(1,d-PARMS[8,species+1]):d-1]==PARMS[12,species+1],c(1,2),sum)==1,                    #Yes ->Has is been that period for less than required days?
                                              -2,                                                                                                          #Yes -> Continue in post-fledging status
                                              ifelse(rbinom(iterations,1,PARMS[5,species+1])==1, 1, -9)),                                                  #No -> Draw according to DBROOD  f(x) ==1? Yes -> Start 2nd brood; No -> Go inactive
                                       ifelse(apply(NS[,max(1,d-PARMS[11,species+1]):d-1],c(1,2),max)>0,                                         #No (in post-failure status) -> Has it been less than minimum number of days since failure?
                                              0,                                                                                                           #Yes -> Continue in post-failure status 
                                              ifelse(d<(PARMS[1,species+1]-PARMS[12,species+1]), 
                                                     ifelse(rbinom(iterations,1,0.9)==1,1,-9),
                                                     -9)))))                                                #No -> Is there time to renest (>= Nesting cycle days till end of season)? Yes -> Start renest; No -> Go inactive
        
        
        Suc<-ifelse(NS[,d]==PARMS[12,species+1],1,0)                        #Identify all successful nest completions on that day
        Complete<-Complete+Suc                              #Update the number of broods completed
        
        #Estimate the brood sizes for completed nests on the day and add that to the cumulative productivity estimate
        Prod.day[,d]<-FL*Suc
        
      }     #Exit season loop
      
      #Write the cumulative productivity for the season into the annual productivity array
      Prod[rep,s,] <- apply(Prod.day,1,sum)
      Success[rep,s,1]<-sum(Complete>0)/iterations
      Success[rep,s,2]<-sum(Complete>1)/iterations
      
      
    }       #Exit the Pb level loop
    
  } #Exit replicate loop
  
  save(Prod,file=paste0("RESULTS/MANUSCRIPT/Pb_ProdSim_Prod_",species.list[species],"_2.RData"))
  save(Success,file=paste0("RESULTS/MANUSCRIPT/Pb_ProdSim_Success_",species.list[species],"_2.RData"))
  
}

#### For EABL ########################################################################
species <- 5


NS.post<-load("DATA/MANUSCRIPT/Pb_EABL_NestSurv.RData")
FL.post<-load("DATA/MANUSCRIPT/Pb_EABL_Fledge.RData")

iterations<-length(FL.post$sims.list$deviance)  #identify size of posterior sample to create simulations

Prod <- array(0,dim=c(reps,length(slspb.sim),iterations))         #Initiate the annual productivity array that will hold annual productivity values
Success<-array(0,dim=c(reps,length(slspb.sim),2))

for (rep in 1:reps){
  
  #Start annual loop
  for (s in 1:length(slspb.sim)){
    
    pb.sim <- (log(slspb.sim[s]+0.1)-pb.mean)/pb.sd
    
    NS <- array(NA,dim=c(iterations,PARMS[1,species+1]))                   #Create an array to track Nest Status (NS) through breeding season
    NS[,1:3] <- 0                                                          #Set first three days of season to no activity (0) for nest building
    NS[,4]<-1                                                              #Start nesting activity on 4th day of season (mask filters out cells not within region)
    
    Prod.day<-matrix(0,iterations,PARMS[1,species+1])                        #Reset productivity and nest completions to zero each year
    Complete<-numeric(iterations)
    
    #Estimate logit of global-interactive nest survival model
    mu.NS <-  NS.post$sims.list$a + NS.post$sims.list$mu.site + NS.post$sims.list$mu.year + NS.post$sims.list$pb*pb.sim
    ds<-exp(mu.NS)/(1+exp(mu.NS))    #Transform NS_global from the logit scale
    mu.FL <- pmin(exp(FL.post$sims.list$a + FL.post$sims.list$mu.site + FL.post$sims.list$mu.year + FL.post$sims.list$pb*pb.sim),PARMS[13,species+1])
    
    # print(paste0("Mean =",mean(ds)))
    # hist(ds,main=paste0(slspb.sim[s]))
    # hist(mu.FL,main=paste0(slspb.sim[s]))
    
    
    #Start season loop; Iterate through  days from May 15th to August 30
    for(d in 5:PARMS[1,species+1]){
      
      #Draw nest fate and brood size replicate for each day
      FL <- pmin(rpois(iterations,mu.FL),PARMS[13,species+1])
      fate<-matrix(rbinom(iterations,1,ds))   #Derive the fate of cells based on daily survival with random draws from binomial dist. 
      
      #Calculate the status of each cell's nest for each day based on if/thens
      NS[,d] <- ifelse(NS[,d-1]==-9,                                                                                  #Is inactive?
                       -9,                                                                                                      #Yes -> Stay inactive
                       ifelse(NS[,d-1]>0,                                                                                       #No (active) -> Is it currently nesting?
                              ifelse(NS[,d-1]==PARMS[12,species+1],                                                                    #Yes -> Has nest survived for full length of days (considered successful)
                                     ifelse(Complete==2, -9, -2),                                                                              #Yes -> Is it the 2nd brood? Yes-> go inactive for season; No ->go into post-fledging period prior to 2nd attempt
                                     NS[,d-1]*fate+fate),                                                                                      #No -> Simulate fate of active nest for that day (e.g., a nest that survives day 5 will go to day 6, else it will go to 0)
                              ifelse(NS[,d-1]==-2,                                                                                     #No -> Is it in post-fledging status?
                                     ifelse(apply(NS[,max(1,d-PARMS[8,species+1]):d-1]==PARMS[12,species+1],c(1,2),sum)==1,                    #Yes ->Has is been that period for less than required days?
                                            -2,                                                                                                          #Yes -> Continue in post-fledging status
                                            ifelse(rbinom(iterations,1,PARMS[5,species+1])==1, 1, -9)),                                                  #No -> Draw according to DBROOD  f(x) ==1? Yes -> Start 2nd brood; No -> Go inactive
                                     ifelse(apply(NS[,max(1,d-PARMS[11,species+1]):d-1],c(1,2),max)>0,                                         #No (in post-failure status) -> Has it been less than minimum number of days since failure?
                                            0,                                                                                                           #Yes -> Continue in post-failure status 
                                            ifelse(d<(PARMS[1,species+1]-PARMS[12,species+1]), 
                                                   ifelse(rbinom(iterations,1,0.9)==1,1,-9),
                                                   -9)))))                                                #No -> Is there time to renest (>= Nesting cycle days till end of season)? Yes -> Start renest; No -> Go inactive
      
      
      Suc<-ifelse(NS[,d]==PARMS[12,species+1],1,0)                        #Identify all successful nest completions on that day
      Complete<-Complete+Suc                              #Update the number of broods completed
      
      #Estimate the brood sizes for completed nests on the day and add that to the cumulative productivity estimate
      Prod.day[,d]<-FL*Suc
      
    }     #Exit season loop
    
    #Write the cumulative productivity for the season into the annual productivity array
    Prod[rep,s,] <- apply(Prod.day,1,sum)
    Success[rep,s,1]<-sum(Complete>0)/iterations
    Success[rep,s,2]<-sum(Complete>1)/iterations
    
    
  }       #Exit the Pb level loop
  
} #Exit replicate loop

save(Prod,file=paste0("RESULTS/MANUSCRIPT/Pb_ProdSim_Prod_",species.list[species],"_2.RData"))
save(Success,file=paste0("RESULTS/MANUSCRIPT/Pb_ProdSim_Success_",species.list[species],"_2.RData"))