library("lme4")
library(MuMIn)
library(multcomp)
library("ggplot2")

#Logistic function from Ben Bolker and Mark Herzog code available here https://rpubs.com/bbolker/logregexp
logexp <- function(exposure = 1) {
  linkfun <- function(mu) qlogis(mu^(1/exposure))
  linkinv <- function(eta) plogis(eta)^exposure
  logit_mu_eta <- function(eta) {
    ifelse(abs(eta)>30,.Machine$double.eps,
           exp(eta)/(1+exp(eta))^2)
  }
  mu.eta <- function(eta) {
    exposure * plogis(eta)^(exposure-1) *
      logit_mu_eta(eta)
  }
  valideta <- function(eta) TRUE
  link <- paste("logexp(", deparse(substitute(exposure)), ")",
                sep="")
  structure(list(linkfun = linkfun, linkinv = linkinv,
                 mu.eta = mu.eta, valideta = valideta,
                 name = link),
            class = "link-glm")
}

#Place the file "Data.csv" in your working drive 
DATA <- read.csv("Data.csv", sep=";")

mod_glmer<-glmer(Fate~Species+Day+Height+BandingWeek+Zone+BandingWeek*Zone+BandingWeek*Species+Zone*Species+Height*Species+Height*Zone+(1|Year/TreeID), family=binomial(link=logexp(DATA$t)), data=DATA)
options(na.action="na.fail")
modelselection<-dredge(mod_glmer,rank="AICc") #CPU time 52min
modelselection

################################
# Best model : m1<-glmer(Fate~Species+Day+Height+BandingWeek+Zone+BandingWeek*Zone+BandingWeek*Species+Zone*Species+Height*Species+(1|Year/TreeID), family=binomial(link=logexp(DATA$t)), data=DATA)
##################################
m1<-get.models(modelselection,1)[[1]]
summary(m1)

DATA$predm1 <-predict(m1, type="response") 
DATA$Height<- factor(DATA$Height, levels = c("low", "med", "high"))

#Figure Height effect per Species from Model 1
a<-data.frame(Species=gl(3,3,labels=c("BUBIBI","EGRGAR","PLEFAL")),Height=rep(unique(DATA$Height),3),SampleSize=c(5,71,48,23,69,7,15,64,34))

FigHeight<-ggplot(DATA, aes(x=Species, y=predm1, fill=Height))+
  geom_boxplot(position=position_dodge(0.6))+
  xlab("Species")+
  ylab("Daily survival rate")+
  scale_fill_brewer() + 
  theme_classic()+
  scale_x_discrete(labels=c("Cattle Egret","Little Egret","Glossy Ibis"))+
  coord_cartesian(ylim=c(0.7, 1))+
  geom_text(data=a,aes(y=1.01,label=SampleSize),size=3,position=position_dodge(width=0.6))
FigHeight

#Figure Zone effect per species from Model 1
aggregate(predm1~Species+Zone, data=DATA,FUN=mean) #for Discussion on the 4% reduced survival of Glossy Ibis's DSR
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") #color blind friendly palette

FigZone<-ggplot(DATA, aes(x=Species, y=predm1, fill=Zone))+
  geom_boxplot(position=position_dodge(0.6))+
  xlab("Species")+
  ylab("Daily survival rate")+
  scale_fill_manual(values=cbPalette)+
  scale_x_discrete(labels=c("Cattle Egret","Little Egret","Glossy Ibis"))+
  coord_cartesian(ylim=c(0.75, 1))
FigZone

####################
# Height effect and Species effect from Model 3 with no interaction (most parsimonious)
####################

m3<-get.models(modelselection,3)[[1]]
summary(m3)
summary(glht(m3, mcp(Height="Tukey")))
boxplot(predict(m3,type="response")~Height, data=DATA, ylab="Chick survival")
summary(glht(m3, mcp(Species="Tukey")))
boxplot(predict(m3,type="response")~Species, data=DATA, ylab="Chick survival")

###################################
#Banding week effect from Model 4
##################################
m4<-get.models(modelselection,4)[[1]]
summary(m4)
boxplot(predict(m4,type="response")~BandingWeek, data=DATA, ylab="Chick survival")
DATA$predm4 <-predict(m4, type="response") 
aggregate(predm4~BandingWeek, data=DATA,FUN=mean)
aggregate(predm4~BandingWeek, data=DATA,FUN=sd)

############################
#Complementary analysis
#Editor asked for differences among species in the two separate dataset (control vs. disturb zones). 
########################################################
#Control zone dataset
DATACont<-DATA[DATA$Zone=="Control",]
ModCont<-glmer(Fate~Species+Day+Height+Species*Height+(1|Year/TreeID), family=binomial(link=logexp(DATACont$t)), data=DATACont)
modelselectionCont<-dredge(ModCont,rank="AICc") 
modelselectionCont
m2Cont<-get.models(modelselectionCont,2)[[1]] #first model did not include the variable "Species"
summary(m2Cont)#model nearly unidentifiable
summary(glht(m2Cont, mcp(Species="Tukey"))) 
plot(predict(m2Cont,type="response")~Species, data=DATACont, ylab="Chick survival")
# same survival for three species in control zone

#Disturbed zone dataset
DATADist<-DATA[DATA$Zone=="Disturbed",]
ModDist<-glmer(Fate~Species+Day+Height+Species*Height+(1|Year/TreeID), family=binomial(link=logexp(DATADist$t)), data=DATADist)
modelselectionDist<-dredge(ModDist,rank="AICc") 
modelselectionDist
m1Dist<-get.models(modelselectionDist,1)[[1]]
summary(m1Dist)
summary(glht(m1Dist, mcp(Species="Tukey")))
plot(predict(m1Dist,type="response")~Species, data=DATADist, ylab="Chick survival")
# Lower survival for Glossy ibis