setwd("C:/Users/g23b661/Desktop")

library(car)
library(ggplot2)

#enter the counts of unreplicated females wit cubs within the Yellowstone Demographic
#Monitoring Area, from the USFWS (2016) Conservation Plan (page 38)
#these are the basic population index data they use to track population growth
#and can be converted to units of total population size

year=seq(from=1984,to=2014,by=1)

# enter the counts of 'unreplicated females with cubs of the year' = FCOY in USFWS conservation plan
# within the GYE demographic monitoring area = DMA in USFWS conservation plan

N=c(23,18,28,17,21,16,25,38,40,21,23,38,37,39,37,36,54,48,59,47,58,31,45,54,57,40,53,48,57,60,59)

YGB <- data.frame(year,N)

# calculate annual growth rates r[a]
# recall that the intrinsic rate of increase (r)= ln(lambda), lambda = e^r
# and of course we can't get lambda for the last year
ln.lambda=log(N[-1]/N[-length(N)]) 

# the [-1] in the numerator drops the first observation and the [-length(N)] drops the last
# so that the value calculated is N next year divided by N this year
# then takes the log

# just renaming the growth rate variable to help keep things clear
r=ln.lambda

# examine the growth rates, r[a]
r

##############Summary of the data and inspection
summary(r)
hist(r)
################

# calculate mean & sd of log(lambda) values
# we need the geometric mean of the lambda values, which is
#equivalent to the arithmetic mean of r, as explained in class

geomean=mean(r)


#examine the geometric mean intrinsic rate of increase
geomean

# and the variance and standard deviation are calculated from the ln(lambda) values, 
# because the growth rates (lambdas) are lognormally distributed while the r values are
# normally distributed

var.r=var(r)
sd.r=sd(r)
var.r
sd.r


# examine the actual population dynamics that were observed


plot(N~year, xlab = "Year", ylab = "Population Index", pch = 19, ylim=c(0,100),
     type="b", main="Unreplicated Grizzly Females with Cubs within Yellowstone DMA ",cex.lab=1.5,
     cex=1,cex.axis=1.5)
############## 

# get some initial parameter estimates to use in nls()
# which is a nonlinear least squares regression function
# that can be used to fit a logistic (Pearl-Verhulst) growth model

coef(lm(logit(N/70)~year,data=YGB))

#and then fit the PV model using nls()

fcoy.growth <- nls(N~phi1/(1+exp(-(phi2+phi3*year))),
    start = list(phi1=70,phi2=-172,phi3=0.086), data= YGB, trace= TRUE)

summary(fcoy.growth)  

# ph1 = k is the carrying capacity, ph2 = r is the intrinsic rate of increase and ph3 = Nzero is the pop size in year zero
# which does not mean anything unless years would be started at year 0 by subtracting 1984 (or whatever
# the first year of the time series is)

#set parameters
phi1<-coef(fcoy.growth)[1]
phi2<-coef(fcoy.growth)[2]
phi3<-coef(fcoy.growth)[3]
x<-c(min(YGB$year):max(YGB$year)) #construct a range of x values bounded by the data
y<-phi1/(1+exp(-(phi2+phi3*x)))   #predicted pop size (in units of fcoy)
predict<-data.frame(x,y)          #create the prediction data frame

#plot
ygb.logistic <- ggplot(data=YGB,aes(x=year,y=N))+
  geom_point(color='blue',size=5)+theme_bw()+
  labs(x='Year',y='Population Index')+
  scale_x_continuous(breaks=c(1985,1990,1995,2000,2005,2010,2015))+
  scale_y_continuous(breaks=c(0,10,20,30,40,50,60,70,80))+
  theme(axis.text=element_text(size=18),axis.title=element_text(size=24))+
  geom_line(data=predict,aes(x=x,y=y), size=1)

ygb.logistic

ggsave(filename = "ygb_logistic_growth.jpg", plot= ygb.logistic, device= "jpeg", dpi=600)

#####################
# Fit the Pearl Verhulst equation as a state space model to partition out sampling error in the counts, 
# leaaving only process error in the variances of parameters in the growth model
#
# doing this with Bayesian model fitting, followung 
# Kery & Schaub 2012 Bayesian POpulation Analysis uing WInBUGS: a hierarchical perspective. Academic Press
#
# but using JAGS rather than BUGs, see link to code this was bsed on at
# https://www.vogelwarte.ch/de/projekte/publikationen/bpa/

library(lattice)
library(coda)
library(R2WinBUGS)
library(R2jags)

jags.dir <- "C:\\Program Files\\JAGS\\JAGS-4.3.0"
setwd("C:/Users/g23b661/Desktop")


n.years <-as.numeric(length(YGB$N))

# Specify model in JAGS language
sink("ssm.jags")
cat("
model { 
# Priors and constraints
N.est[1] ~ dunif(-500, 500)            # Prior for initial population size

mean.r ~ dunif(-10, 10)          # Prior for growth rate
mean.k ~ dunif(0,100)           # Prior for carrying capacity

#these vars not used this model
# sigma.k ~ dunif(0, 100)           # Prior for sd of k
# sigma2.k <- pow(sigma.k, 2)
# tau.k <- pow(sigma.k, -2)
# 
# sigma.r ~ dunif(0, 100)           # Prior for sd of r
# sigma2.r <- pow(sigma.r, 2)
# tau.r <- pow(sigma.r, -2)


sigma.obs ~ dunif(0, 100)           # Prior for sd of observation process
sigma2.obs <- pow(sigma.obs, 2)
tau.obs <- pow(sigma.obs, -2)

# Likelihood
# State process
for (t in 1:(T-1)){
   N.est[t+1] <- N.est[t]*exp(mean.r*(1-(N.est[t]/mean.k)))
   }
# Observation process
for (t in 1:T) {
   y[t] ~ dnorm(N.est[t], tau.obs)
   }
}
",fill = TRUE)
sink()

# Bundle data
jags.data <- list(y = YGB$N, T = n.years)

# Initial values
inits <- function(){list( mean.r = runif(1, -0.5, 0.5),
                          mean.k = runif(1,20,100),
                         sigma.obs = runif(1, 0, 10), N.est = c(runif(1, 20, 40), rep(NA, (n.years-1))))} 

# Parameters monitored
parameters <- c("mean.r", "mean.k",  "sigma2.obs", "N.est")

# MCMC settings
ni <- 25000
nt <- 3
nb <- 10000
nc <- 3

# Call JAGS from R (JRT <1 min)
ssm <- jags(jags.data, inits, parameters, "ssm.jags", n.chains = nc, n.thin = nt,
            n.iter = ni, n.burnin = nb, working.directory = getwd())
plot(ssm)
#traceplot(ssm)



# Define function to draw a graph to summarize results
graph.ssm <- function(ssm, N, y){
  fitted <- lower <- upper <- numeric()
  n.years <- length(y)
  for (i in 1:n.years){
    fitted[i] <- mean(ssm$BUGSoutput$sims.list$N.est[,i])
    lower[i] <- quantile(ssm$BUGSoutput$sims.list$N.est[,i], 0.025)
    upper[i] <- quantile(ssm$BUGSoutput$sims.list$N.est[,i], 0.975)}
  m1 <- min(c(y, fitted, N, lower))
  m2 <- max(c(y, fitted, N, upper))
  par(mar = c(4.5, 4, 1, 1), cex = 1.2)
  plot(0, 0, ylim = c(m1, m2), xlim = c(0.5, n.years), ylab = "Population Index (FCOY)", xlab = "Year", 
       las = 1, col = "black", type = "l", lwd = 2, frame = FALSE, axes = FALSE)
  axis(2, las = 1)
  axis(1, at = seq(0, n.years, 5), labels = seq(1983, 1983+n.years, 5))
  axis(1, at = 0:n.years, labels = rep("", n.years + 1), tcl = -0.25)
  polygon(x = c(1:n.years, n.years:1), y = c(lower, upper[n.years:1]), col = "gray90", border = "gray90")
  points(N, type = "b", col = "red", lwd = 2)
#  points(y, type = "l", col = "black", lwd = 2)
  points(fitted, type = "l", col = "blue", lwd = 2)
  legend(x = 1, y = m2, legend = c("Counts", "Bayesian Pearl-Verhulst \nState Space Model"), lty = c(1, 1, 1), lwd = c(2, 2, 2), col = c("red",  "blue"), bty = "n", cex = 1)
}

# Execute function: Produce figure 
graph.ssm(ssm, N, y)

#examine output
print(ssm,3)
median(ssm$BUGSoutput$sims.list$mean.k)
median(ssm$BUGSoutput$sims.list$mean.r)
# posterior CI for r and k, with histograms

hist(ssm$BUGSoutput$sims.list$mean.k, nclass = 30, col = "light gray", main = "", xlab = "Carrying Capacity (K)",
     ylab = "Frequency",
     xlim = c(20,100))
abline(v=median(ssm$BUGSoutput$sims.list$mean.k), lwd = 3)
text(x=80, y=1300, labels = "Median K = 62")

hist(ssm$BUGSoutput$sims.list$mean.r, nclass = 50, col = "light gray", main = "", xlab = "Intrinsic Rate of Increase (r)",
      ylab = "Frequency")
abline(v=median(ssm$BUGSoutput$sims.list$mean.r), lwd = 3)
text(x=0.25, y=700, labels = "Median r = 0.101")


plot(density(ssm$BUGSoutput$sims.list$mean.r), xlab = 'Intrinsic Rate of Increase (r)', lwd = 3, main =
       "Yellowstone Grizzly Bear 1984-2014", cex.lab = 1.3)
abline(v=median(ssm$BUGSoutput$sims.list$mean.r), lwd = 2, col = 'blue')


plot(density(ssm$BUGSoutput$sims.list$mean.k),  xlab = 'Carrying Capacity (K)', lwd = 3, main =
               "Yellowstone Grizzly Bear 1984-2014", cex.lab = 1.3)
abline(v=median(ssm$BUGSoutput$sims.list$mean.k), lwd = 2, col = 'blue')


#######################################################
#
#PVA using projection of density-dependent state-space model

nsims= 100
T = 50
N0 = 60
threshold = 20

N.agg <- matrix(NA, T, nsims)
N.agg[1,]=N0

#outer loop is stepping across repeated simulations
for (j in 1:nsims){      

  #inner loop is stepping across years
   for (t in 2:T){
     
     #draw values of r and K in the growth modelfrom their stochastic distributions  
     k.stoch <- sample(ssm$BUGSoutput$sims.list$mean.k, size = 1, replace = T)
     r.stoch <- sample(ssm$BUGSoutput$sims.list$mean.r, size = 1, replace = T)
    
      #project population size one time step with those values of r and k
     N.agg[t,j] <- N.agg[(t-1),j]*exp(r.stoch*(1-(N.agg[(t-1),j]/k.stoch)))
     
     #leav the loop if (pseudo) exntinction occurs
     if(N.agg[t,j]<=threshold) break
  
     } # close t

} # close j    

#end PVA

#plot results
plot(x=seq(from =1, to= T, by=1), y=N.agg[,1], 
     ylim = c(min(N.agg), (max(N.agg)+5)), xlab = 'Years', ylab = 'Population Size FCOY', cex.lab = 1.3)
for (i in 1:nsims){
lines(x=seq(from =1, to= T, by=1), y=N.agg[,i], col = i, lwd = 2)
    
}