# F&WL 501 - lab 13 - count-based PVA
# Data source:
# Interagency Grizzly Bear Study Team - Annual Report 2008
# http://www.nrmsc.usgs.gov/files/norock/products/IGBST/2008report_1-15-2010.pdf

year=seq(1983,2008)

N=c(19, 22, 18, 28, 17, 21, 18, 25, 38, 41, 21, 23, 43, 38, 
39, 37, 36, 51, 48, 58, 46, 58, 31, 45, 53, 56)

ln.lambda=log(N[-1]/N[-length(N)])
ln.lambda[length(N)]=NA

# put various values of interest together in data frame
dat=data.frame(cbind(year,N,ln.lambda))
dat

# calculate mean & sd of log(lambda) values
gmean=mean(ln.lambda, na.rm=T)
var.ln.lambda=var(ln.lambda, na.rm=T)
sd.ln.lambda=sqrt(var.ln.lambda)

c(gmean, sd.ln.lambda)

# plot population trajectory
plot(year,N,'o',ylim=c(0,60))
lines(year,N[1]*exp(gmean)^(0:(length(N)-1)),'l',lty=3,col='blue')

# simulate population trajectories

N0=N[length(N)]  # start at most recent value
Ne=25             #  quasi-extinction threshold
time.horizon=25  # number of years to project forward
num.sims=5000   # number of simulations to do

# storage container for population sizes
pop=matrix(NA,time.horizon,num.sims)
pop[1,]=N0

# loop to create population sizes
for (sim in 1:num.sims){
  for (t in 2:time.horizon) {
    lambda=exp(rnorm(1,gmean,sd.ln.lambda))
    pop[t,sim]=pop[(t-1),sim]*lambda
    if(pop[t,sim]<=Ne) break  # stop trajectory if go <Ne
}}

# summarize population size information through time
# calculate SD and CI's on log-scale, i.e., lognormal distribution
pop.mean=apply(pop,1,mean, na.rm=T)
log.pop.sd  =apply(log(pop),1,sd, na.rm=T)
pop.lcl =exp(log(pop.mean)-1.96*log.pop.sd)
pop.ucl =exp(log(pop.mean)+1.96*log.pop.sd)

# make plots
windows()
par(mfrow=c(2,1))
plot(1:time.horizon,log(pop.mean),'l',ylim=c(0,max(log(pop.ucl))),
  xlab='Years into the future', ylab='log(N)')
lines(1:time.horizon,log(pop.lcl),'l',lty=3)
lines(1:time.horizon,log(pop.ucl),'l',lty=3)

plot(1:time.horizon,pop.mean,'l',ylim=c(0,max(pop.ucl)),
  xlab='Years into the future', ylab='N')
lines(1:time.horizon,pop.lcl,'l',lty=3)
lines(1:time.horizon,pop.ucl,'l',lty=3)

# work with quasi-extinction
# evaluate populations for quasi-extinction through time
# copy the pop sizes over to new variable
pop.quasi.ext=pop 

# replace NA values with the value for Ne - needed for later calculation
pop.quasi.ext=ifelse(is.na(pop.quasi.ext),Ne,pop.quasi.ext) 

# substitute a value of 1 if pop is quasi.ext, 0 if not
pop.quasi.ext=ifelse(pop.quasi.ext<=Ne,1,0)

# calculate the ppn that is quasi-ext at each time step
#  by taking the row sums (# of 1's) and dividing by num. sims
ppn.quasi.ext=apply(pop.quasi.ext,1,sum)/num.sims 

cbind(t=1:time.horizon,ppn.quasi.ext)

windows()
par(mfrow=c(1,1))
plot(1:time.horizon,ppn.quasi.ext,'l',ylim=c(0,1))

## diffusion-approximation = analytical solution based on continuous time model
#require(popbio)
#ex=extCDF(gmean, var.ln.lambda, N[length(N)], Ne=Ne, tmax = time.horizon)
#ex
#points(1:time.horizon,ex,'o',col='blue')
## look at results from both
#out=data.frame(cbind(t=1:time.horizon,ppn.qe.sim=ppn.quasi.ext,ppn.qe.diff.approx=ex))
#round(out,2)