#------------------------------------------------------------------------------------------------------------------
# Harmon et al. 2020. Seasonal Patterns in Nest Survival of a subtropical wading bird, the Hawaiian Stilt
#------------------------------------------------------------------------------------------------------------------

#Load data
log.exp<-read.csv("C:/Users/kmcor/OneDrive/Graduate School/Papers/Nesting Time Paper/Nesting Time Paper R Datasheets/Logistic Exposure_HI Stilt Nesting Data_Clean.csv", stringsAsFactors = F)

log.exp$Year<-as.factor(log.exp$Year)

#--------------------------------------------------------------------------------------------------------------------

##Get summary of habitat variables at nests
library("psych")

#Distance to water
describe(log.exp$DistWatm, type=2)
#Vegetation height
describe(log.exp$VegHtcm, type=2)

#------------------------------------------------------------------------------------------------------------------------

##Changes in nest-site variables over time

#Linear regression for veg height over time
vegheight.lm<-lm(VegHtcm~FirstFound, data=log.exp)
summary(vegheight.lm)


#Linear regression for distance to water over time
distwater.lm<-lm(DistWatm~FirstFound, data=log.exp)
summary(distwater.lm)

#----------------------------------------------------------------------------------------------------------------------------------------

##Test for effects of other potential predictors
install.packages("car")
library("car")

#Obtain DSR for nests with and without cameras
surv.cam<-glm(Survive~Camera, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)
summary(surv.cam) #Not signifcantly different 
pred.data <- data.frame(Camera=unique(log.exp$Camera))
logit.dsr.pred.cam <- predict(surv.cam,
                              newdata=pred.data, se.fit=TRUE)
dsr.cam <- cbind(pred.data,
                 logit.dsr=logit.dsr.pred.cam$fit,
                 logit.dsr.se=logit.dsr.pred.cam$se.fit)
dsr.cam$dsr <- expit(dsr.cam$logit.dsr)
dsr.cam$dsr.se <- dsr.cam$logit.dsr.se* dsr.cam$dsr *
  (1-dsr.cam$dsr)
dsr.cam

#Null model
surv.null<-glm(Survive~1, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)
summary(surv.null)
dsr.null<-expit(coef(surv.null))
dsr.se.null<-arm::se.coef(surv.null)*dsr.null*(1-dsr.null)

#Compare camera model with null model
Anova(surv.cam) #Camera variable does not have signifcant impact

#Analysis of Deviance Table (Type II tests)
#Response: Survive
#LR Chisq Df Pr(>Chisq)
#Camera    1.218  1     0.2697

#-----------------------------------------------------------------
#Test for effects of island
log.exp.island<-read.csv("C:/Users/kmcor/OneDrive/Graduate School/Papers/Nesting Time Paper/Nesting Time Paper R Datasheets/Logistic Exposure_HI Stilt Nesting Data_Islands.csv", stringsAsFactors = F)
log.expand.island <- expand.nest.data(log.exp.island)
table(log.exp.island$Island)

surv.island<-glm(Survive~Island, family=binomial(logexp(na.omit(log.expand.island$Exposure))), data=log.expand.island)
summary(surv.island) #Not significantly different 
pred.data.island <- data.frame(Island=unique(log.exp.island$Island))
logit.dsr.pred.island <- predict(surv.island,
                                 newdata=pred.data.island, se.fit=TRUE)
dsr.island <- cbind(pred.data.island,
                    logit.dsr=logit.dsr.pred.island$fit,
                    logit.dsr.se=logit.dsr.pred.island$se.fit)
dsr.island$dsr <- expit(dsr.island$logit.dsr)
dsr.island$dsr.se <- dsr.island$logit.dsr.se* dsr.island$dsr *
  (1-dsr.island$dsr)
dsr.island

#Compare island model with null model
Anova(surv.island) #Island model performs slightly better than null model

#Analysis of Deviance Table (Type II tests)
#Response: Survive
#LR Chisq Df Pr(>Chisq)
#Island  0.32456  1     0.5689

#-----------------------------------------------------------------
#Test for effects of mammal-exclusion fence
log.exp.fence<-read.csv("C:/Users/kmcor/OneDrive/Graduate School/Papers/Nesting Time Paper/Nesting Time Paper R Datasheets/Logistic Exposure_HI Stilt Nesting Data_Fence.csv", stringsAsFactors = F)
log.exp.fence$Year <- as.factor(log.exp.fence$Year)
log.expand.fence <- expand.nest.data(log.exp.fence)

surv.fence<-glm(Survive~Year, family=binomial(logexp(na.omit(log.expand.fence$Exposure))), data=log.expand.fence)
summary(surv.fence) #Not signifcantly different 
pred.data.fence <- data.frame(Year=unique(log.exp.fence$Year))
logit.dsr.pred.fence<- predict(surv.fence,
                               newdata=pred.data.fence, se.fit=TRUE)
dsr.fence<- cbind(pred.data.fence,
                  logit.dsr=logit.dsr.pred.fence$fit,
                  logit.dsr.se=logit.dsr.pred.fence$se.fit)
dsr.fence$dsr <- expit(dsr.fence$logit.dsr)
dsr.fence$dsr.se <- dsr.fence$logit.dsr.se* dsr.fence$dsr *
  (1-dsr.fence$dsr)
dsr.fence

#Compare island model with null model
Anova(surv.fence) #Model including year performs better than null model
#Analysis of Deviance Table (Type II tests)
#Response: Survive
#LR Chisq Df Pr(>Chisq)    
#Year   15.231  1  9.515e-05 ***


#-----------------------------------------------------------------------------------------------------------------------------------------

#Test for effects of wetland and year
surv.wetland<-glm(Survive~Wetland, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)
surv.year<-glm(Survive~Year, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)
Anova(surv.wetland)
Anova(surv.year)

#----------------------------------------------------------------------------

##Logistic Exposure Model for Nest Survival

Look for correlations among variables
install.packages("Hmisc")
library(Hmisc)

stilt_mat <- as.matrix(log.exp[, c(9,10)])
rcorr(stilt_mat, type="pearson")


#Expand data to create daily survival records
expand.nest.data <- function( nestdata, 
                              FirstFound="FirstFound", 
                              LastChecked="LastChecked",
                              LastPresent="LastPresent",
                              Fate="Fate",
                              AgeDay1="AgeDay1"){
  require(plyr)
  nestdata2 <- plyr::adply(nestdata, 1, function(x){
    # expand for days when nest is known to be alive
    # create data frame because tibbles can cause problems
    x <- as.data.frame(x)
    #browser()
    orig.x <- x
    x <- x[rep(1,x[,LastPresent]-x[,FirstFound]),]  
    #browser()
    if(nrow(x)>0){
      x$Day <- x[,FirstFound][1]:(x[,LastPresent][1]-1)
      x$Exposure<- 1
      x$Survive <- 1
      #browser()
      if(AgeDay1 %in% names(orig.x)) x$NestAge <- orig.x[,AgeDay1] + x$Day -1 
    }
    
    # now for the final record where the nest fails in interval
    # We define the date of the last interval as the midpoint of the interval
    # We define the nestage in the last interval as the age in the midpoint of the interval as well.
    #browser()
    if(orig.x[,LastChecked] > orig.x[,LastPresent]){
      x2 <- orig.x
      x2$Exposure <- orig.x[,LastChecked] - orig.x[,LastPresent]
      x2$Survive  <- 1-orig.x[,Fate]
      x2$Day      <- x2[,LastPresent] + x2$Exposure/2
      if(AgeDay1 %in% names(orig.x)) x2$NestAge  <- orig.x[,AgeDay1] + x2$Day -1
      x <- rbind(x, x2)
      #browser()
    }
    # remove any records with 0 exposure
    x <- x[ !x$Exposure==0,]
    # return the expanded data set
    x
  })
  nestdata2}



#Expand the data    
log.expand <- expand.nest.data(log.exp)
head(log.expand)


#modified link function
logexp<-function(days=1)
{
  linkfun <- function(mu) qlogis(mu^(1/days))
  linkinv <-function(eta) plogis(eta)^days
  mu.eta<-function(eta) days*plogis(eta)^(days-1) *
    .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
  valideta<-function(eta) TRUE
  link<-paste("logexp(",deparse(substitute(days)),")",sep="")
  structure(list(linkfun=linkfun,linkinv=linkinv,
                 mu.eta=mu.eta,valideta=valideta,name=link),
            class="link-glm")
}


#Construct models
surv1<-glm(Survive~Day + I(Day^2) + VegHtcm + DistWatm + Day*VegHtcm + Day*DistWatm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv2<-glm(Survive~Day +  Day*VegHtcm + Day*DistWatm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv3<-glm(Survive~Day +  DistWatm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv4<-glm(Survive~Day +  VegHtcm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv5<-glm(Survive~Day + DistWatm + Day*DistWatm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv6<-glm(Survive~Day + VegHtcm + Day*VegHtcm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv7<-glm(Survive~Day + I(Day^2), family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv8<-glm(Survive~Day, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv9<-glm(Survive~VegHtcm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv10<-glm(Survive~DistWatm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv11<-glm(Survive~DistWatm + VegHtcm, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)

surv12<-glm(Survive~1, family=binomial(logexp(na.omit(log.expand$Exposure))), data=log.expand)


#Compare models
install.packages("MuMIn")
library("MuMIn")
model.sel(surv1, surv2, surv3, surv4, surv5, surv6, surv7, surv8, surv9, surv10, surv11, surv12)
#------------------------------------------------------------------------------

##Plot top model results

#Plot DSR Day + Veg Height

#Make list of date values based on exposure data for 2018 and 2019
day.surv6 <- data.frame(Day=seq(min(log.expand$Day), max(log.expand$Day), length.out = 153), VegHtcm=mean(log.expand$VegHtcm))
veg.surv6<-data.frame(VegHtcm=seq(0, 70, by = 5), Day=mean(log.expand$Day))

##Model predictions for graphing
pred.date.surv6 <- predict(surv6, day.surv6, re.form=NA, se.fit=TRUE, type="link")
pred.veg.surv6 <- predict(surv6, veg.surv6, re.form=NA, se.fit=TRUE, type="link")


#Predicted values and 95% CI (to be used in new prediction dataframe)
date.values.surv6 <- pred.date.surv6$fit #predicted values
date.low.surv6 <- pred.date.surv6$fit - pred.date.surv6$se.fit*1.44 #85% CI lower bound
date.high.surv6 <- pred.date.surv6$fit + pred.date.surv6$se.fit*1.44 #85% CI upper bound

veg.values.surv6 <- pred.veg.surv6$fit #predicted values
veg.low.surv6 <- pred.veg.surv6$fit - pred.veg.surv6$se.fit*1.44 #85% CI lower bound
veg.high.surv6 <- pred.veg.surv6$fit + pred.veg.surv6$se.fit*1.44 #85% CI upper bound


#Convert predicted values to scale of DSR (to be used in new prediction dataframe)
date.values.surv6.scale <- exp(1)^date.values.surv6/(1+exp(1)^date.values.surv6)
date.low.surv6.scale <- exp(1)^date.low.surv6/(1+exp(1)^date.low.surv6)
date.high.surv6.scale <- exp(1)^date.high.surv6/(1+exp(1)^date.high.surv6)

veg.values.surv6.scale <- exp(1)^veg.values.surv6/(1+exp(1)^veg.values.surv6)
veg.low.surv6.scale <- exp(1)^veg.low.surv6/(1+exp(1)^veg.low.surv6)
veg.high.surv6.scale <- exp(1)^veg.high.surv6/(1+exp(1)^veg.high.surv6)


#New prediction dataframe for date
date.surv6.df <- data.frame(day.surv6, date.values.surv6.scale, date.low.surv6.scale, date.high.surv6.scale)
veg.surv6.df <- data.frame(veg.surv6, veg.values.surv6.scale, veg.low.surv6.scale, veg.high.surv6.scale)


#Create column for Constance Survival
date.surv6.df$constant <- rep("Constant Survival",length(date.surv6.df$Day))
colnames(date.surv6.df) <- c("Day", "VegHtcm", "DSR","DSRLow","DSRHigh","Constant Survival")
date.surv6.df$`Constant Survival`  #not sure what this does? From Lindsey's code

veg.surv6.df$constant <- rep("Constant Survival",length(veg.surv6.df$VegHtcm))
colnames(veg.surv6.df) <- c("VegHtcm", "Day", "DSR","DSRLow","DSRHigh","Constant Survival")
veg.surv6.df$`Constant Survival`  #not sure what this does? From Lindsey's code


#Plot the predicted DSR for date and veg height
library("ggplot2")

ggplot(data=date.surv6.df, aes(x=Day, y=DSR)) +
  geom_line(size=0.75) +
  scale_x_continuous(limits = c(1, 155)) + 
  scale_y_continuous(limits = c(0.90, 1)) +
  geom_ribbon(aes(ymin=DSRLow, ymax=DSRHigh),alpha=0.15) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))


ggplot(data=veg.surv6.df, aes(x=VegHtcm, y=DSR)) +
  geom_line(size=0.75) +
  scale_x_continuous(limits = c(0, 70)) + 
  scale_y_continuous(limits = c(0.97, 1)) +
  geom_ribbon(aes(ymin=DSRLow, ymax=DSRHigh),alpha=0.15) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))


#Plot interaction between veg height and day of nesting season
install.packages("interplot")
library("interplot")

interplot(m = surv6, var1 = "Day", var2 = "VegHtcm", ci = .85, point = T) +
  theme(axis.text.x  = element_text(angle=90))


day.sim <- seq(from = 1, to = 153, by = 1) #list of dates of nesting season
eff.veg <- coef(surv6)[3] + coef(surv6)[4] * day.sim #calculate effect of veg height for different days of nesting season
eff.data<-data.frame(day.sim, eff.veg) #create data frame of days and effects of veg based on day
eff.data$se.eff.veg <- sqrt(vcov(surv6)[3, 3] + 
                              eff.data$day.sim^2 * vcov(surv6)[4, 4] + 
                              2 * eff.data$day.sim * vcov(surv6)[3, 4]) #get standard errors for each variable of model
eff.data$CIlow <- eff.data$eff.veg - eff.data$se.eff.veg*1.44 #85% CI lower bound
eff.data$CIhigh <- eff.data$eff.veg + eff.data$se.eff.veg*1.44 #85% CI upper bound


ggplot(data=eff.data, aes(x=day.sim, y=eff.veg)) +
  geom_line(size=0.75) +
  geom_ribbon(aes(ymin=CIlow, ymax=CIhigh),alpha=0.15) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))