#' # Analysis of mountain goat data from Alaska and Washington
#' 
#' 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?  

# Clear out working memory 
  rm(list = ls())

#' ## Mountain goats in Alaska
#' 
# set global option for output when producing html file
  opts_chunk$set(comment=NA, fig.width=6, fig.height=6)

#+ warning=FALSE, message=FALSE
#' load libraries
  library(lme4)
  library(doBy)
  library(gplots)
  library(mosaic)
  library(AICcmodavg)
  library(MuMIn)
  
#' Function for getting SE
  SE<-function(x){sqrt(var(x)/length(x))} 
  
#' Read in data:
  mtg<-read.csv("../Data/MTG_Sight_Alaska.csv")
  
#' Look at number of observations for each animal
  table(mtg$ID)
  
#' Distribution of time since capture
  plot(density(mtg$Day.Since.Capture, from=0), xlab="Days since capture", ylab="Density", main="")  
  
#' Since so few animals were observed more than once, fit a glm (assuming independence), 
#' and then use bootstrap for standard errors if necessary (e.g., if
#' statistically significant when independence is assumed).
  fit.mtg<-glm(Seen.~Day.Since.Capture, family=binomial(), data=mtg)
  summary(fit.mtg) 
  
  
#' Could also consider controlling for area if pdetect depeneded on area.... 
  fit.mtg2<-glm(Seen.~Day.Since.Capture+Area2, family=binomial(), data=mtg)
  summary(fit.mtg2) 
   
  
#' Days since capture does not appear to be significant.  What about "years since capture" (since the number
#' of days since capture is bi-modal) 
  mtg$Itime<-I(mtg$Day.Since.Capture > 200) 
  fit.mtg3<-glmer(Seen.~Itime + (1|ID), family=binomial(), nAGQ=10, data=mtg)
  summary(fit.mtg3)
 
#'  Test and confidence interval for a difference in proportions 
  prop.test(Seen.~Itime, data=mtg)
  with(data=mtg,prop.table(table(Seen., Itime), margin=2))
  
  
#' ## Mountain Goats in WA
#'   
#' Read in data:
  mtgwa<-read.csv("../Data/WAgoats.csv")
 
#' Fit full model, with interaction between yr.s.cap & method, then model without the interaction
  fit.glmer1<-glmer(observed~grpsize+vegmidpt + terrainob+ method + yr.s.cap + method*yr.s.cap + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer1)
  
#' Drop interaction  
  fit.glmer2<-glmer(observed~grpsize+vegmidpt + method + terrainob+ yr.s.cap + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer2)
  
#' Drop other non-significant variables (control only for groupsize) to test for robustness.  First, include interaction
  fit.glmer3<-glmer(observed~grpsize+ method + yr.s.cap + method*yr.s.cap+(1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer3)
  
#' Now, drop interaction   
  fit.glmer3b<-glmer(observed~grpsize+ method + yr.s.cap +(1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer3b)
  
#' Now, method and yr.s.cap as only predictors (along with groupsize)  
  fit.glmer4<-glmer(observed~grpsize+ method + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer4)
  
  fit.glmer5<-glmer(observed~grpsize+  yr.s.cap + (1|id), family=binomial(), data=mtgwa, nAGQ=10 )  
  summary(fit.glmer5)


#' Yr.s.cap is not even close to significant in any of the models.  Construct plots to look at effect size.  
#' 
#' AK mountain goats: basic proportions and confidence intervals for proportion seen vs. time since cap
  mtg.ak<-summaryBy(Seen.~Itime, FUN=c(mean, SE),data=mtg) 
  
#' Model based predictions for Alaskan mtgs
  xmat<-matrix(rbind(c(1,0), c(1,1)), 2,2)
  mt.ak.prd<- matrix(NA,2,3)
  mt.ak.prd[,1]<-plogis(xmat%*%fixef(fit.mtg3))
  se.prd<-diag(xmat%*%vcov(fit.mtg3)%*% t(xmat)) 
  mt.ak.prd[,3]<- plogis(xmat%*%fixef(fit.mtg3)+1.96*se.prd)
  mt.ak.prd[,2]<- plogis(xmat%*%fixef(fit.mtg3)-1.96*se.prd)

#' Now for WA mtgs.  Use fit.glmer5, with grpsize and yr.s.cap (lowest AIC of models containing yr.s.cap)
#' Get predictions for median group size (which varied a bit by yr.s.cap, but not that much)
  tapply(mtgwa$grpsize, mtgwa$yr.s.cap, mean)
  tapply(mtgwa$grpsize, mtgwa$yr.s.cap, median)
  mgs<-median(mtgwa$grpsize)
  xmat<-matrix(rbind(c(1,mgs,0), c(1,mgs,1), c(1,mgs,2)), 3,3)
  mt.wa.prd<- matrix(NA,3,3)
  mt.wa.prd[,1]<-plogis(xmat%*%fixef(fit.glmer5))
  se.prd.wa<-diag(xmat%*%vcov(fit.glmer5)%*% t(xmat)) 
  mt.wa.prd[,3]<- plogis(xmat%*%fixef(fit.glmer5)+1.96*se.prd.wa)
  mt.wa.prd[,2]<- plogis(xmat%*%fixef(fit.glmer5)-1.96*se.prd.wa)

#' Now, get proportion observed as function of time since cap
  mtgwa$Iobs<-I(mtgwa$observed=="Y")
  mtg.wa.ob<-summaryBy(Iobs~yr.s.cap,FUN=c(mean, SE), data=mtgwa)
   
# Now, plot with helicopter/ground as well for washington data, using glmer3 (with interaction) between method & yr since cap
  xmat<-matrix(rbind(c(1,mgs,0,0,0),c(1,mgs,0,1,0),c(1,mgs,0,2,0), 
                    c(1,mgs,1,0,1),c(1,mgs,1,1,1),c(1,mgs,1,2,1)), 6,5)
  mt.wa.prd<- matrix(NA,6,3)
  mt.wa.prd[,1]<-plogis(xmat%*%fixef(fit.glmer3))
  se.prd.wa<-diag(xmat%*%vcov(fit.glmer3)%*% t(xmat)) 
  mt.wa.prd[,3]<- plogis(xmat%*%fixef(fit.glmer3)+1.96*se.prd.wa)
  mt.wa.prd[,2]<- plogis(xmat%*%fixef(fit.glmer3)-1.96*se.prd.wa)

#' # Plot results  
#+ fig.width=8, fig.height=6
par(mfrow=c(1,1), bty="L", cex.lab=1.8, cex=1, cex.axis=1.8, mar=c(6,6,4.1,2.1) , xpd=T, oma=c(0,1,0,0)) 
  plotCI(c(1,2), mt.ak.prd[,1], ui=mt.ak.prd[,3], li=mt.ak.prd[,2], xlab="Years since capture", ylab="Probability of detection", ylim=c(0,1), xlim=c(0,2), gap=0, pch=17, cex=2,xaxt="n")
  axis(side=1, c(0,1,2))
  plotCI(c(c(0,1,2)-.02, c(0,1,2)+.02), mt.wa.prd[,1], ui=mt.wa.prd[,3], li=mt.wa.prd[,2], xlab="", ylab="",add=T, ylim=c(0,1), xlim=c(0,2), gap=0, pch=c(rep(16,3), rep(1,3)),cex=2, xaxt="n")   
  legend("bottomright", pch=c(17,16,1), cex=2,lty=1, title="", legend=c("Alaska Helicopter/Helicopter", "Washington Ground/Helicopter", "Washington Helicopter/Helicopter"), bty="n")
  
  
#' For WA, where there were no missing values, also calculate model-averaged coefficients to get at effect sizes (use all models except 
#' those with interactions between method and years since capture; these models had larger AIC's anyway) 
  Cands <- list(fit.glmer2, fit.glmer3b, fit.glmer4, fit.glmer5)
   
#' model selection
  aictab(cand.set = Cands)

#' model average for difference between period 1 vs period 4
  summary(mod.avg(fit.glmer2, fit.glmer3b, fit.glmer4, fit.glmer5))
  
#' Effect sizes:
#'  
#'  - Yr.s.cap 
  c(exp(-0.161),  exp(-0.161-1.96*0.46395), exp(-0.161+1.96*0.46395))
  
#' - Capture mehtod  
  c(exp(0.1188361),exp(0.1188361-1.96*0.4996),exp(0.1188361+1.96*0.4996))
  
    
# Print out Session info
  sessionInfo()