#' # Analysis of moose data from MN
#' 
#' Programmer:  John Fieberg
#' 
#' R Code in support of Fieberg et al. Do Capture and Survey Methods Influence Whether Marked Animals are Representative of Unmarked Animals?  
#' 
#' Goal:  explore whether detection rates for individual animals changed as a function of time since they were captured. 
#' 
  rm(list = ls()) 

#+ warning=FALSE, message=FALSE
# Load libraries
  library(lme4)
  library(ggplot2)
  library(nlme)
  library(doBy)
  library(mosaic)
  
#' First read in detection data...
  sight.dat<-read.csv("../data/sightdat.csv")
   
#'  Evaluate whether there is a trend in the age-distrubion relative to year of capture
  boxplot(sight.dat$age~sight.dat$year, xlab="", ylab="Age distribution")  
  
#' Now, get years since capture
  sight.dat$yr.s.cap<-sight.dat$year-sight.dat$capyear

#' One observation was missing the date of capture, delete it:
  sight.dat<-sight.dat[is.na(sight.dat$capyear)!=T,]

#' Get number of individuals by capture year (use median or mean, below - they give the same answer as they should)
  capyrs<-summaryBy(capyear~id, data=sight.dat, FUN=c(median, mean))  
  table(capyrs$capyear.median)
    
#' Look at correlation between age and time since capture, and also the number of observations missing age 
  with(sight.dat[is.na(sight.dat$age)!=T,], cor.test(x=age, y=yr.s.cap))
  table(is.na(sight.dat$age))
  
#' Net guns were used in 2002, darting in all other years
  sight.dat$net<-0
  sight.dat$net[sight.dat$capyear==2002]<-1
  prop.table(table(sight.dat$net, sight.dat$observed), margin=1) # net gunners were more likely to be detected
  chisq.test(table(sight.dat$net, sight.dat$observed)) # not significant, but small sample size
  
#' Look at simple table and plot of detection probabilities as a function of time since capture
  temp.tab<-prop.table(table(sight.dat$observed, sight.dat$yr.s.cap), margin=2)
  temp.tab
 
#' Note:  a chi-squared test of Ho: detection probabilities were constant across years 
#' ends up being non-significant (but, the test is not geared towards alterantives where 
#' there is a monotonic increase/decrease).  Also, the test does not adjust for any
#' annual differences in detection that may be due to the distribution of animals in varying levels of cover/visual obstruction.  
  chisq.test(table(sight.dat$observed, sight.dat$yr.s.cap), simulate.p.value=T)
  
#' Detection models:  
#' 
#'   1. Full model with voc, time since capture, net gun (y/n), and age 
#'   2. ~ voc + time since capture
#'   3. ~ voc + net gun(Y/n)
#'   4. ~ yr.s. cap
#'   5. ~ net gun 
    sight.dat$id<-as.factor(sight.dat$id)
    fit.glmer1<-glmer(observed~voc+yr.s.cap + net + age + (1|id), family=binomial(), data=sight.dat, nAGQ=10 )
    summary(fit.glmer1)
    fit.glmer2<-glmer(observed~voc+yr.s.cap + (1|id), family=binomial(), data=sight.dat, nAGQ=10 )
    summary(fit.glmer2)
    fit.glmer3<-glmer(observed~voc+ net +(1|id), family=binomial(), data=sight.dat, nAGQ=10 )
    summary(fit.glmer3) 
    fit.glmer4<-glmer(observed~yr.s.cap+(1|id), family=binomial(), data=sight.dat, nAGQ=10 )
    summary(fit.glmer4) 
    fit.glmer5<-glmer(observed~ net +(1|id), family=binomial(), data=sight.dat, nAGQ=10 )
    summary(fit.glmer5) 
  
# Can't compare models with and without age using AIC because of missing age values.  Look at amount of missing data, below:
  table(is.na(sight.dat$age))
  
#' For plotting purposes fit a model that assumes independence (since re variance for animal is ~0)
  fit1<-glm(observed~voc+yr.s.cap, family=binomial(), data=sight.dat)
  summary(fit1)
  newdat<-data.frame(observe=1, voc=mean(sight.dat$voc), yr.s.cap=1:5)
  p1<-predict(fit1, newdat, se=T)
  
#' Odds of detection per unit increase in years since capture (to report effect size)
  exp(coef(fit1)[3])
  exp(coef(fit1)[3]+ 1.96*sqrt(vcov(fit1)[3,3]))
  exp(coef(fit1)[3]- 1.96*sqrt(vcov(fit1)[3,3]))
  
#Uncomment to produce single, 2 panel figure 
#par(mfrow=c(1,2), bty="L", cex.lab=1.8, cex.axis=1.8, mar=c(6,5,4.1,2.1) ) 

#+fig.height=6, fig.width=8, fig.keep="last", fig.show='asis'
  plot(as.numeric(colnames(temp.tab)), temp.tab[2,], xlab="Years since capture", ylab="Detection Probability", ylim=c(0,.8), main="", xlim=c(0.5, 5.5), pch=16)
  text(1:5, temp.tab[2,], paste("n = ", table(sight.dat$yr.s.cap), sep=""), pos=3, cex=1.4)
  lines(1:5,plogis(p1$fit))
  lines(1:5,plogis(p1$fit+1.96*p1$se.fit), lty=2)
  lines(1:5,plogis(p1$fit-1.96*p1$se.fit), lty=2)
  mtext(side=3, outer=F, cex=2, "A)")
  
#+fig.height=6, fig.width=8, fig.keep="last", fig.show='asis'
  boxplot(voc~yr.s.cap, data=sight.dat, ylab="Visual Obstruction", xlab="Years since capture")

# Alternative  plot

#+fig.height=6, fig.width=8, fig.keep="last", fig.show='asis'
  densityplot(~voc, groups=as.factor(yr.s.cap), data=sight.dat, xlab="Visual Obstruction", auto.key=TRUE) 
  mtext(side=3, outer=F, cex=2, "B)")
  
# Lets look at voc relative to time since capture
  fit.voc1<-lmer(voc~yr.s.cap+(1|id), data=sight.dat)
  summary(fit.voc1)
  fit.voc2<-lmer(voc~age+yr.s.cap+(1|id), data=sight.dat)
  summary(fit.voc2)
  
# Use lme to get approximate pvalue
   fit.voc.lme<-lme(voc~yr.s.cap, random=list(id=~1), data=sight.dat)
   summary(fit.voc.lme)
  
# Figure (for presentation)
  plot(as.numeric(colnames(temp.tab)), temp.tab[2,], xlab="Years since capture", ylab="Detection Probability", ylim=c(0,.8), main="", xlim=c(0.5, 5.5), pch=16)
  text(1:5, temp.tab[2,], paste("n = ", table(sight.dat$yr.s.cap), sep=""), pos=3, cex=1.4)
  lines(1:5,plogis(p1$fit))
  lines(1:5,plogis(p1$fit+1.96*p1$se.fit), lty=2)
  lines(1:5,plogis(p1$fit-1.96*p1$se.fit), lty=2)

  
# Print out Session info
  sessionInfo()