## Population size estimates for paper "Status, genetic diversity,
## habitat suitability and conservation implications for a fragmented 
## Asian elephant (Elephas maximus) population in Cambodia"  

install.packages("capwire")
library(capwire)

## Create model dataframe. Requires a two column dataframe 
## 1. Capture class is the number of times and one individual has been captured
## 2. Number of Individuals 

## Prey Lang Population Estimate
## Unique IDs split into Prey Lang and Preah Roka/Chhaeb individuals using the genetic analysis report
## UG17 1 sample removed due to anomaly
PLcounts <- c(8,3,11,1,3,6,2,1,2,2,1,2,1,2,1,1,3,1,2,1,1,3) ## Counts retrieved 
## Sample count from Table 2 Genetics analysis report
PreyLang <- buildClassTable(PLcounts) ## Counts aggregated into capture classes to create final dataframe for the model
PreyLang

## Fitting the Equal Capture Model
## Assumes all individuals have an equal chance of being captured
res.ecmPL <- fitEcm(data=PreyLang, max.pop=100)
## Maximum population is much larger than reasonably expect population to be to generate upper bounds for model
res.ecmPL
##Confidence intervals are estimated using parametric bootstrapping
boot.ecmPL <- bootstrapCapwire(x=res.ecmPL, bootstraps = 1000, CI= c(0.025, 0.975))
boot.ecmPL

## Fitting the Two-Innate Rates Model (TIRM) to estimate population size
res.tirmPL <- fitTirm(data=PreyLang, max.pop=100)
res.tirmPL


## Model selection
## To select between the ECM and the TIRM models use a likelihood ratio test.
## To evaluate the significance of likelihood ration, generate a null 
# distribution for the LR by simulating data under the ECM using the
# maximum likelihood pop est. from fitting the ECM.

lik.ratio <- lrtCapwire(ecm=res.ecmPL, tirm = res.tirmPL, bootstraps = 100)
lik.ratio

## Reject null hypothesis that there is no population heterogeneity, use TIRM model
## Perform parametric bootstraps to obtain 95% confidence intervals for est. using TIRM

boot.tirmPL <- bootstrapCapwire(x=res.tirmPL, bootstraps = 1000, CI= c(0.025, 0.975))
boot.tirmPL

## PREY LANG POP ESTIMATE = 28 CONFIDENCE INTERVALS 2.5% = 23 97.5% = 38


###############Preah Roka/Chhaeb Population Estimate############################

PRCcounts <- c(8,5,7,1,7,4,8,1,1,1,1,1,7) #UG23 1 sampled removed due to anomaly

PreahRokaCh <- buildClassTable(PRCcounts)
PreahRokaCh
res.ecmPRC <- fitEcm(data= PreahRokaCh, max.pop = 100)
res.ecmPRC
boot.ecmPRC <- bootstrapCapwire(x=res.ecmPRC, bootstraps = 1000, CI= c(0.025, 0.975))
boot.ecmPRC

## Fitting the Two-Innate Rates Model (TIRM) to estimate population size

res.tirmPRC <- fitTirm(data=PreahRokaCh, max.pop=100)
res.tirmPRC


## Model selection
## To select between the ECM and the TIRM models use a likelihood ratio test.
## To evaluate the significance of likelihood ration, generate a null 
# distribution for the LR by simulating data under the ECM using the
# maximum likelihood pop est. from fitting the ECM.

lik.ratio <- lrtCapwire(ecm=res.ecmPRC, tirm = res.tirmPRC, bootstraps = 100)
lik.ratio

boot.tirmPRC <- bootstrapCapwire(x=res.tirmPRC, bootstraps = 1000, CI= c(0.025, 0.975))
boot.tirmPRC

## PREAH ROKA/CHHAEB POP ESTIMATE = 18 CONFIDENCE INTERVALS 2.5% = 12, 97.5% = 21

## Graph final results
install.packages(ggplot2)
library(ggplot2)

popest <- c(31,20)
PAs <- c("Prey Lang", "Preah Roka/Chhaeb")
df <- data.frame(PAs, popest)
df

ymin <- c(24,13)
ymax <- c(41,22)

g1 <- ggplot(df, aes(x= PAs, y= popest)) + geom_point() + theme_bw() + 
  theme(panel.grid.major.x = element_blank(), 
        panel.grid.major.y = element_line(colour = "black",
                                          linetype = "dotted")) +
  labs(x= "Population", y = "Population Estimate") +
  geom_errorbar(aes(ymin = ymin, ymax = ymax), width=0.2)
#geom_errorbar(data= CI,mapping = aes(x= PAs, ymin= 13, ymax=22))
g1