library(MASS) 
library(lme4)
library(reshape2)
library(ggplot2)
library(bbmle)
library(MuMIn)
library(ggpubr)
library(dplyr)
library(dfpk)
library(RColorBrewer)
library(arm)

###  Data = Logistic_Exposure_1.xlsx #####
dat2 <- subset(peerj_62023_Logistic_Exposure,
               select=c(NestID,DaysPast,Year,Fate,FolDen,TrHt,NestHt,DistEdg,TrHtavg,JunCt))

dat2S <- subset(dat2,DaysPast>0)

dat2<-dat2S[complete.cases(dat2S), ]

### Preliminary Data Plotting
mdat <- melt(dat2,id.var=1:4)

ggplot(mdat,aes(x=value,y=Fate))+
  geom_point()+
  facet_wrap(~variable,scales="free")+
  geom_smooth(method="loess")+
  scale_y_continuous(expand=c(0,0), limits=c(0.5,1.0), oob=scales::rescale_none)
###Juniper Count is unimodal, will use I(JunCt^2)

###Check correlation coefficients of independent variables
cor(dat2[,5:10])                #nothing over 0.7
###Convert numerical to factors
dat2$Year<-as.factor(dat2$Year)
dat2$Fate<-as.factor(dat2$Fate)

#####################################################
# For Logistic Exposure Models we need to source a custom link function. Courtesy of Ben Bolker. 
# Note that this version of the logexp link function will allow the use of
# lme4 (glmer) to fit mixed effect models as well.  

logexp <- function(exposure = 1) { 
  linkfun <- function(mu) qlogis(mu^(1/exposure)) 
  ## FIXME: is there some trick we can play here to allow 
  ## evaluation in the context of the 'data' argument? 
  linkinv <- function(eta) plogis(eta)^exposure 
  logit_mu_eta <- function(eta) { 
    ifelse(abs(eta)>30,.Machine$double.eps, 
           exp(eta)/(1+exp(eta))^2) 
    ## OR .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats") 
  } 
  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")
} 

###Standardize continuous variables so that they're all on similar scales
normFunc <- function(x){(x-mean(x, na.rm = T))/sd(x, na.rm = T)}
dat2[5:10]<-apply(dat2[5:10],2, normFunc)

### Make Global Model to Dredge
globalmod<-glmer(Fate~TrHtavg*NestHt+FolDen*DistEdg+TrHt+JunCt+I(JunCt^2)+NestHt*FolDen+(1|Year), family=binomial(link=logexp(dat2$DaysPast)), data=dat2,na.action = na.pass)
summary(globalmod)
###No highly correlated independent variables (> 0.7)

### Dredge Global Model and Export as CSV
msubset <- expression(  dc(JunCt, `I(JunCt^2)`)) ###makes it so I(JunCt) only included in models with JunCt
dm<-dredge(globalmod, beta = F, evaluate = T, rank = AICc,subset=msubset)

dm

#Pull Out Competitive Models
m1 <- glmer(Fate ~ FolDen*DistEdg+(1|Year), family=binomial(link=logexp(dat2$DaysPast)), data=dat2) 
summary(m1)
confint.merMod(m1,"FolDen",method="profile")
confint.merMod(m1,"DistEdg",method="profile")
confint.merMod(m1,"FolDen:DistEdg",method="profile")

m2 <- glmer(Fate ~ FolDen*DistEdg+NestHt*TrHtavg+(1|Year), family=binomial(link=logexp(dat2$DaysPast)), data=dat2) 
summary(m2)
confint.merMod(m2,"FolDen",method="profile")
confint.merMod(m2,"DistEdg",method="profile")
confint.merMod(m2,"FolDen:DistEdg",method="profile")
confint.merMod(m2,"NestHt",method="profile")
confint.merMod(m2,"TrHtavg",method="profile")
confint.merMod(m2,"NestHt:TrHtavg",method="profile")

m3 <- glmer(Fate ~ FolDen*DistEdg+JunCt+I(JunCt^2)+(1|Year), family=binomial(link=logexp(dat2$DaysPast)), data=dat2) 
summary(m3)
confint.merMod(m3,"FolDen",method="profile")
confint.merMod(m3,"DistEdg",method="profile")
confint.merMod(m3,"FolDen:DistEdg",method="profile")
confint.merMod(m3,"JunCt",method="profile")
confint.merMod(m3,"I(JunCt^2)",method="profile") ### 95% CI overlap with Zero

NullM <- glmer(Fate ~ 1+(1|Year), family=binomial(link=logexp(dat2$DaysPast)), data=dat2) 
summary(NullM)
#Odds Ratios
cc <- confint(m1,parm="beta_")  ## slow (~ 11 seconds)
ctab <- cbind(est=fixef(m1),cc)
rtab <- exp(ctab)
print(rtab,digits=3)

cc <- confint(m2,parm="beta_")  ## slow (~ 11 seconds)
ctab <- cbind(est=fixef(m2),cc)
rtab <- exp(ctab)
print(rtab,digits=3)



######### Revert Back to Unstandardized Logistic_Exposure for DNS Estimate/Plotting ############
dat2 <- subset(peerj_62023_Logistic_Exposure,
               select=c(NestID,DaysPast,Year,Fate,TrHt,FolDen,NestHt,DistEdg,TrHtavg,FolDenavg,JunCt))

dat2S <- subset(dat2,DaysPast>0)

dat2<-dat2S[complete.cases(dat2S), ]
dat2$Year<-as.factor(dat2$Year)
dat2$Fate<-as.factor(dat2$Fate)

m1<-glm(Fate ~ FolDen*DistEdg, family=binomial(link=logexp(dat2$DaysPast)),data=dat2)
m2<-glm(Fate ~ FolDen*DistEdg+NestHt*TrHtavg, family=binomial(link=logexp(dat2$DaysPast)),data=dat2)


###Model Averaging Top Models

top.models<-list(m1,m2)
x<-model.avg(top.models)
summary(x)

###DNS prediction Using Means of Independent Variables
ndata <- expand.grid(DistEdg = mean(dat2$DistEdg),FolDen=mean(dat2$FolDen),NestHt=mean(dat2$NestHt),TrHtavg=mean(dat2$TrHtavg))

preds <- predict(x,newdata=ndata, se.fit=TRUE)
preds

plottable <- with(preds, 
                  
                  data.frame(ndata, 
                             
                             fit = invlogit(fit),
                             
                             lwr = invlogit(fit + 2*se.fit),
                             
                             upr = invlogit(fit - 2*se.fit)))
###DNS prediction
plottable$fit

###Cumulative Nest survival prediction
(plottable$fit)^32.4


#DNS Figure 2 and 3
ndata <- expand.grid(JunCt = seq(5,30,5),DistEdg = seq(1, 250, 25),FolDen=seq(1,4,1),NestHt=seq(1,4,0.5),TrHtavg=seq(2,5,1),JunCt=seq(1,34,5))

preds <- predict(m1,newdata=ndata, se.fit=TRUE)

plottable <- with(preds, 
                  
                  data.frame(ndata, 
                             
                             fit = invlogit(fit),
                             
                             lwr = invlogit(fit + 2*se.fit),
                             
                             upr = invlogit(fit - 2*se.fit)) )
x <- plottable$FolDen
plottable$FD_3group <-
  case_when(x == 1 ~ "0-25%",
            x == 2 ~ "25-50%",
            x ==  3 ~ "50-75%",
            x == 4 ~ "75-100%")

ip1<-ggplot(plottable, aes(x = DistEdg, y = fit, col=FD_3group)) + 
  geom_smooth(aes(group=FD_3group),method=loess)+
  labs(x = "Distance From Edge (cm)", y = "Daily Nest Survival Probability",col = "Foliage Density")+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"), 
        axis.title.y = element_text( colour="black", size = 12),
        axis.title.x = element_text( colour="black", size = 12),
        axis.text.x=element_text(colour="black", size= 10,margin=unit(c(0.3,0.1,0.1,0.1), "cm")), 
        axis.text.y=element_text( colour="black", size = 10,margin=unit(c(0.1,0.3,0.1,0.1), "cm")),
        plot.margin = unit(c(0.,0.0,0,0),"in"),
        axis.ticks.length=unit(-0.15, "cm"),
        legend.text = element_text(colour="black", size = 10))

ip1

ggsave(file="Figure2.tiff", ip1,units="in", width=7, height=3.5, dpi=600, compression = 'lzw')

preds2 <- predict(m2,newdata=ndata, se.fit=TRUE)

plottable2 <- with(preds2, 
                  
                  data.frame(ndata, 
                             
                             fit = invlogit(fit),
                             
                             lwr = invlogit(fit + 2*se.fit),
                             
                             upr = invlogit(fit - 2*se.fit)) )

ip2<-ggplot(plottable2, aes(x = NestHt, y = fit, col=TrHtavg)) + 
  geom_smooth(aes(group=TrHtavg),method=loess)+
  labs(x = "Nest Height (m)", y = "Daily Nest Survival Probability",col = "Average Tree Height (m)")+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"), 
        axis.title.y = element_text( colour="black", size = 12),
        axis.title.x = element_text( colour="black", size = 12),
        axis.text.x=element_text(colour="black", size= 10,margin=unit(c(0.3,0.1,0.1,0.1), "cm")), 
        axis.text.y=element_text( colour="black", size = 10,margin=unit(c(0.1,0.3,0.1,0.1), "cm")),
        plot.margin = unit(c(0.,0.0,0,0),"in"),
        axis.ticks.length=unit(-0.15, "cm"),
        legend.text = element_text(colour="black", size = 10))

ip2

ggsave(file="Figure3.tiff", ip2,units="in", width=7, height=3.5, dpi=600, compression = 'lzw')

### Bonus Plot of Optimal Juniper Density
m3<-glm(Fate ~ FolDen*DistEdg+JunCt+I(JunCt^2), family=binomial(link=logexp(dat2$DaysPast)),data=dat2)

preds <- predict(m3,newdata=ndata, se.fit=TRUE)

plottable <- with(preds, 
                  
                  data.frame(ndata, 
                             
                             fit = invlogit(fit),
                             
                             lwr = invlogit(fit + 2*se.fit),
                             
                             upr = invlogit(fit - 2*se.fit)) )

ip3<-ggplot(plottable, aes(x = JunCt, y = fit)) + 
  stat_smooth(method = "lm", formula = y ~ x + I(x^2))+
  labs(x = "Juniper Count (25-m radius)", y = "Daily Nest Survival Probability",col = "Foliage Density")+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        panel.background = element_blank(),axis.line = element_line(colour = "black"), 
        axis.title.y = element_text( colour="black", size = 10),
        axis.title.x = element_text( colour="black", size = 10),
        axis.text.x=element_text(colour="black", size= 8,margin=unit(c(0.3,0.1,0.1,0.1), "cm")), 
        axis.text.y=element_text( colour="black", size = 8,margin=unit(c(0.1,0.3,0.1,0.1), "cm")),
        plot.margin = unit(c(0.,0.0,0,0),"in"),
        axis.ticks.length=unit(-0.15, "cm"),
        legend.text = element_text(colour="black", size = 8))

ip3

###Optimum Juniper Density at the 25-m radius scale is ~15-20% of the available cover.