#Functions for simulating population growth and estimated realized growth
#Code by Subhash R. Lele and Mark L. Taper  12/12/2011 9:21:46 PM

library(plotrix) #used to plot confidence intervals

setwd("C:/YOUR_WORKING_DIRECTORY_HERE") #set the working directory

ProjectFemales=function(N0,T,sa,f){
#N0=initial number of females
#T=total number of transitions (years-1)
#sa=adult survival
#f=fertility=probability that an adult female who survives produces a female calf who survives to next census period
#reasonably, f < sa/2 to accommodate for sex ratio and greater mortality of calves
  vec.1=function(N){rep(1,N)}
  N.spr=list()
  for (t in 1:T){
    N.spr[[t]]=rbinom(n=N0,size=1,prob=sa)*(vec.1(N0)+rbinom(n=N0,size=1,prob=f))
    N0=sum(N.spr[[t]])
  }
  return(N.spr)
}

lambda.est=function(N.spr.vec,sa.smp.sz,R.smp.sz,sa.rb=0, cc.rb=0,cv2cw){
#N.spr=vector of 0s, 1s, and 2s representing cows who have died, cows who have
#survived and cows who have survied and produced female calves who have survived
  N.sa.smp=sample(N.spr.vec,size=sa.smp.sz,replace=F)  #sample from which to estimate sa
  sa.hat=mean(N.sa.smp>0)*(1+sa.rb)  #estimate of sa
  N.R.smp=sample(N.spr.vec,size=R.smp.sz,replace=F) #sample for R estimate
  cv=sum(N.R.smp>1) # number of calves in sample
  cw=sum(N.R.smp>0) # number of cows in sample
  cw.ob=cw*(1+cc.rb) #number of cows after apportionment of unknowns
  misc=rbinom(1,cv,cv2cw) #number of calves misclassified as cows
  R.hat=(cv-misc)/(cw.ob+cv) #estimate of R.hat with bias
  lambda.hat=sa.hat/(1-R.hat)
  return(lambda.hat)
}

lambda.tru.calc=function(N0,N.spr){
  T=length(N.spr)
  N=rep(NA,T+1)
  lambda.tru=rep(NA,T)
  N[1]=N0
  for (t in 1:T){
    N[t+1]=sum(N.spr[[t]])
    lambda.tru[t]=N[t+1]/N[t]
  }
  return(lambda.tru)
}

lambda.rel.bias=function(sab,f,sa.rb,cc.rb,cv2cw){
#calcualtes relative bias in lambda estimate given true values of sab,f
#and error quantities sa.rb, cc.rb, cv2cw
  cc=sab
  cf=sab*f
  lambda.true=sab/(1-(cf/(cc+cf)))
  lambda.est=sab*(1+sa.rb)/(1-((cf*(1-cv2cw))/(cc*(1+cc.rb)+cf)))
  lambda.rel.bias=(lambda.est/lambda.true)-1
  return(lambda.rel.bias)
}


Boot.RGrw=function(N.spr,bt.sz=2000,cvr=0.95,yr1=1993,
                   sa.smp.sz,R.smp.sz,sa.rb,cc.rb, cv2cw){
  T=length(N.spr)
  yr=(0:(T-1))+yr1
  center=rep(NA,T)
  UCL=rep(NA,T)
  LCL=rep(NA,T)
  BootRG=matrix(rep(NA,T*bt.sz),ncol=T)
  for (b in 1:bt.sz){
    BootRG[b,1]=lambda.est(N.spr[[1]],sa.smp.sz,R.smp.sz,sa.rb,cc.rb, cv2cw)
  }
  for (t in 2:T){
    for (b in 1:bt.sz){
      BootRG[b,t]=BootRG[b,t-1]*lambda.est(N.spr[[t]],sa.smp.sz,R.smp.sz,sa.rb,cc.rb, cv2cw)
    }
  }
  for (t in 1:T){
    RG.bvec=BootRG[,t]
    center[t]=mean(RG.bvec)
    UCL[t]=quantile(RG.bvec,1-((1-cvr)/2))
    LCL[t]=quantile(RG.bvec,(1-cvr)/2)
  }
  out=list(yr=yr,center=center,UCL=UCL,LCL=LCL)
  return(out)
}

N0=300     # initial true number of females in the population
T=16       # number of years overwhich to project
sa=0.866   # true adult survival
f=0.155    # true fertility
N.spr=ProjectFemales(N0,T,sa,f)  #projections of "true" population sizes
lam.tru=lambda.tru.calc(N0,N.spr) #calulation of vector of "true" lambdas
RG.tru=cumprod(lam.tru) #calculation of "true" realized population growth

B.RG.0.0.0=Boot.RGrw(N.spr,bt.sz=4000,cvr=0.95,yr1=1993, #simulation of Realized
                   sa.smp.sz=30,R.smp.sz=100,sa.rb=0,cc.rb=0, cv2cw=0)

B.RG.02.0.0=Boot.RGrw(N.spr,bt.sz=4000,cvr=0.95,yr1=1993,
                   sa.smp.sz=30,R.smp.sz=100,sa.rb=-0.02,cc.rb=.0, cv2cw=0.0)

B.RG.04.0.0=Boot.RGrw(N.spr,bt.sz=4000,cvr=0.95,yr1=1993,
                   sa.smp.sz=30,R.smp.sz=100,sa.rb=-0.04,cc.rb=.0, cv2cw=0.0)

B.RG.06.0.0=Boot.RGrw(N.spr,bt.sz=4000,cvr=0.95,yr1=1993,
                   sa.smp.sz=30,R.smp.sz=100,sa.rb=-0.06,cc.rb=.0, cv2cw=0.0)


pdf("FIGURE.1.0.remake.pdf")  #include to send figure to a pdf file


par(mfrow=c(2,2))
  plot(B.RG.0.0.0$yr,RG.tru,type="b",pch=19,ylim=c(0,1.75), col="blue", cex=1.25,
    main="Unbiased estimates of lambda",xlab="",ylab="Realized Growth")
  dispersion(B.RG.0.0.0$yr,B.RG.0.0.0$center,B.RG.0.0.0$UCL,B.RG.0.0.0$LCL,
     interval=F,arrow.cap=0)
  points(B.RG.0.0.0$yr,RG.tru,type="l",lty=1)
  abline(h=1,lty="dotted")
  
  plot(B.RG.02.0.0$yr,B.RG.02.0.0$center,type="b",pch=19,ylim=c(0,1.75),
    main="2% biased estimates of lambda",xlab="", col="blue", cex=1.25,ylab="")
  dispersion(B.RG.02.0.0$yr,B.RG.02.0.0$center,B.RG.02.0.0$UCL,B.RG.02.0.0$LCL,
       interval=F,arrow.cap=0)
  points(B.RG.0.0.0$yr,RG.tru,type="l",lty=1)
  abline(h=1,lty="dotted")
  
  plot(B.RG.04.0.0$yr,B.RG.04.0.0$center,type="b",pch=19,ylim=c(0,1.75),
    main="4% biased estimates of lambda",xlab="", col="blue", cex=1.25,ylab="Realized Growth")
  dispersion(B.RG.04.0.0$yr,B.RG.04.0.0$center,B.RG.04.0.0$UCL,B.RG.04.0.0$LCL,
       interval=F,arrow.cap=0)
  points(B.RG.0.0.0$yr,RG.tru,type="l",lty=1)
  abline(h=1,lty="dotted")
  
  plot(B.RG.06.0.0$yr,B.RG.06.0.0$center,type="b",pch=19,ylim=c(0,1.75),
    main="6% biased estimates of lambda",xlab="", col="blue", cex=1.25,ylab="")
  dispersion(B.RG.06.0.0$yr,B.RG.06.0.0$center,B.RG.06.0.0$UCL,B.RG.06.0.0$LCL,
      interval=F,arrow.cap=0)
  points(B.RG.0.0.0$yr,RG.tru,type="l",lty=1)
  abline(h=1,lty="dotted")
  

dev.off()  #returns output to console  need this on when pdf is on