## Chapter 5 Analysis
## Josephine Belia Smit

##Part 1: Activity analysis
##Part 2: Modeling event counts as a function of risk, diel period, group type, and camera placement relative to water and roads


##########################################
### Part 1: Activity Profile Analysis#####
##########################################

#Load packages
library(overlap) #version 0.3.3, published 20 May 2020
library(activity) #version 1.3, published 9 September 2019

#read in dataset
dat <- read.csv("All_Events.csv")
head(dat)

#event start times in 'sun time'
timeST <- dat$start_suntime

#extract data for RNP
rnpST <- timeST[dat$Grid == 1]

#extract data for MIO
mioST <- timeST[dat$Grid == 2]

#extract data for MBO
mboST <- timeST[dat$Grid == 3]

#extract data for RUI
ruiST <- timeST[dat$Grid == 4]

#calculate coefficient of overlapping for RNP & MIO 
min(length(rnpST), length(mioST))
rnpmioSTest <- overlapEst(rnpST, mioST, type="Dhat4")
rnpmioSTest
#0.7288065

#create figure
overlapPlot(rnpST, mioST, main="All Events RNP vs. MIO")
legend('topright', c("RNP", "MIO"), lty=c(1,2), col=c(1,4), bty='n')


#calculate coefficient of overlapping for RNP & MBO
min(length(rnpST), length(mboST))
rnpmboSTest <- overlapEst(rnpST, mboST, type="Dhat4")
rnpmboSTest

#create figure
overlapPlot(rnpST, mboST, main="All Events RNP vs. MBO")
legend('topright', c("RNP", "MBO"), lty=c(1,2), col=c(1,4), bty='n')

#overlap for RNP & RUI ST
min(length(rnpST), length(ruiST))
rnpruiSTest <- overlapEst(rnpST, ruiST, type="Dhat4")
rnpruiSTest

#create figure
overlapPlot(rnpST, ruiST, main="All Events RNP vs. RUI)")
legend('topright', c("RNP", "RUI"), lty=c(1,2), col=c(1,4), bty='n')


## Generate 95% confidence intervals for coefficient of overlapping

#RNP vs MIO
#Resample (smoothed bootstrap samples)
rnpSTboot <- resample(rnpST, 10000)
dim(rnpSTboot)
mioSTboot <- resample(mioST, 10000)
dim(mioSTboot)

#This produces matrices with a column for each bootstrap sample. The bootstrap sample
#size is the same as the original sample size. 
#Now we pass these two matrices to bootEst to generate estimates of overlap from each pair
#of samples.

#Bootstrapped mean 
rnpmioST <- bootEst(rnpSTboot, mioSTboot, type="Dhat4") 
(BSmean <- mean(rnpmioST)) 

bootCI(rnpmioSTest, rnpmioST)
#Use the basic0 output from bootCI as your confidence interval. 

#RNP vs MBO
#Resample
rnpSTboot <- resample(rnpST, 10000)
dim(rnpSTboot)
mboSTboot <- resample(mboST, 10000)
dim(mboSTboot)

#Bootstrapped mean 
rnpmboST <- bootEst(rnpSTboot, mboSTboot, type="Dhat4") 
(BSmean <- mean(rnpmboST)) 

bootCI(rnpmboSTest, rnpmboST)
#use basic0 


#RNP vs RUI
#Resample
rnpSTboot <- resample(rnpST, 10000)
dim(rnpSTboot)
ruiSTboot <- resample(ruiST, 10000)
dim(ruiSTboot)

#Bootstrapped mean 
rnpruiST <- bootEst(rnpSTboot, ruiSTboot, type="Dhat4") 
(BSmean <- mean(rnpruiST)) 

bootCI(rnpruiSTest, rnpruiST)
#use basic0 

#fit kernel density curves in the activity package and test whether activity curves are significantly different
rnpele <- timeST[dat$Grid == 1 & dat$Species == 'Elephant']
mioele <- timeST[dat$Grid == 2 & dat$Species == 'Elephant']
mboele <- timeST[dat$Grid == 3 & dat$Species == 'Elephant']
ruiele <- timeST[dat$Grid == 4 & dat$Species == 'Elephant']

rnpelefit <- fitact(rnpele)
mioelefit <- fitact(mioele)
mboelefit <- fitact(mboele)
ruielefit <- fitact(ruiele)

compareCkern(rnpelefit, mioelefit, reps = 1000)

compareCkern(rnpelefit, mboelefit, reps = 1000)

compareCkern(rnpelefit, ruielefit, reps = 1000)


##Compare CC and LB activity curves at low and high-risk grids

#coefficient of overlapping for CC & LB for RNP 
rnpccST <- timeST[dat$Grid == 1 & dat$Group == 'CC']
rnplbST <- timeST[dat$Grid == 1 & dat$Group == 'LB']

min(length(rnpccST), length(rnplbST))
rnpcclbSTest <- overlapEst(rnpccST, rnplbST, type="Dhat4")
rnpcclbSTest

overlapPlot(rnpccST, rnplbST, main="CC vs. LB in RNP")
legend('topright', c("CC", "LB"), lty=c(1,2), col=c(1,4), bty='n')

#coefficent of overlapping for CC & LB for MIO 
mioccST <- timeST[dat$Grid == 2 & dat$Group == 'CC']
miolbST <- timeST[dat$Grid == 2 & dat$Group == 'LB']

min(length(mioccST), length(miolbST))
miocclbSTest <- overlapEst(mioccST, miolbST, type="Dhat1")
miocclbSTest

overlapPlot(mioccST, miolbST, main="CC vs. LB in MIO")
legend('topright', c("CC", "LB"), lty=c(1,2), col=c(1,4), bty='n')

#coefficent of overlapping for CC & LB for MBO
mboccST <- timeST[dat$Grid == 3 & dat$Group == 'CC']
mbolbST <- timeST[dat$Grid == 3 & dat$Group == 'LB']

min(length(mboccST), length(mbolbST))
mbocclbSTest <- overlapEst(mboccST, mbolbST, type="Dhat1")
mbocclbSTest
 
overlapPlot(mboccST, mbolbST, main="CC vs. LB in MBO")
legend('topright', c("CC", "LB"), lty=c(1,2), col=c(1,4), bty='n')

#coefficent of overlapping for CC & LB for RUI ST
ruiccST <- timeST[dat$Grid == 4 & dat$Group == 'CC']
ruilbST <- timeST[dat$Grid == 4 & dat$Group == 'LB']

min(length(ruiccST), length(ruilbST))
ruicclbSTest <- overlapEst(ruiccST, ruilbST, type="Dhat4")
ruicclbSTest

overlapPlot(ruiccST, ruilbST, main="CC vs. LB in RUI")
legend('topright', c("CC", "LB"), lty=c(1,2), col=c(1,4), bty='n')


#calculate 95% confidence intervals for coefficient of overlapping for LB and CC comparisons

#CC vs LB RNP 
rnplbSTboot <- resample(rnplbST, 10000)
dim(rnplbSTboot)
rnpccSTboot <- resample(rnpccST, 10000)
dim(rnpccSTboot)
rnpcclbST <- bootEst(rnpccSTboot, rnplbSTboot, type="Dhat4") 
(BSmean <- mean(rnpcclbST)) 
bootCI(rnpcclbSTest, rnpcclbST)
#use basic0 

#CC vs LB MIO 
miolbSTboot <- resample(miolbST, 10000)
dim(miolbSTboot)
mioccSTboot <- resample(mioccST, 10000)
dim(mioccSTboot)
miocclbST <- bootEst(mioccSTboot, miolbSTboot, type="Dhat1") 
(BSmean <- mean(miocclbST)) 
bootCI(miocclbSTest, miocclbST)
#use basic0   

#CC vs LB MBO 
mbolbSTboot <- resample(mbolbST, 10000)
dim(mbolbSTboot)
mboccSTboot <- resample(mboccST, 10000)
dim(mboccSTboot)
mbocclbST <- bootEst(mboccSTboot, mbolbSTboot, type="Dhat1") 
(BSmean <- mean(mbocclbST)) 
bootCI(mbocclbSTest, mbocclbST)


#CC vs LB RUI
ruilbSTboot <- resample(ruilbST, 10000)
dim(ruilbSTboot)
ruiccSTboot <- resample(ruiccST, 10000)
dim(ruiccSTboot)
ruicclbST <- bootEst(ruiccSTboot, ruilbSTboot, type="Dhat4") 
(BSmean <- mean(ruicclbST)) 
bootCI(ruicclbSTest, ruicclbST)

# use activity package to test if CC and LB curves are significantly different as a function of risk

#extract data for CC and LB by grid
rnpccact <- timeST[dat$Grid == 1 & dat$Species == 'Elephant'& dat$Group == 'CC']
mioccact <- timeST[dat$Grid == 2 & dat$Species == 'Elephant'& dat$Group == 'CC']
mboccact <- timeST[dat$Grid == 3 & dat$Species == 'Elephant'& dat$Group == 'CC']
ruiccact <- timeST[dat$Grid == 4 & dat$Species == 'Elephant'& dat$Group == 'CC']

rnplbact <- timeST[dat$Grid == 1 & dat$Species == 'Elephant'& dat$Group == 'LB']
miolbact <- timeST[dat$Grid == 2 & dat$Species == 'Elephant'& dat$Group == 'LB']
mbolbact <- timeST[dat$Grid == 3 & dat$Species == 'Elephant'& dat$Group == 'LB']
ruilbact <- timeST[dat$Grid == 4 & dat$Species == 'Elephant'& dat$Group == 'LB']

rnplbfit <- fitact(rnplbact)
miolbfit <- fitact(miolbact)
mbolbfit <- fitact(mbolbact)
ruilbfit <- fitact(ruilbact)

rnpccfit <- fitact(rnpccact)
mioccfit <- fitact(mioccact)
mboccfit <- fitact(mboccact)
ruiccfit <- fitact(ruiccact)

#test if LB activity curves are significantly different between low and high risk grids
compareCkern(rnplbfit, miolbfit, reps = 1000)

compareCkern(rnplbfit, mbolbfit, reps = 1000)

compareCkern(rnplbfit, ruilbfit, reps = 1000)

#test if CC activity curves are significantly different between low and high risk grids
compareCkern(rnpccfit, mioccfit, reps = 1000)

compareCkern(rnpccfit, mboccfit, reps = 1000)

compareCkern(rnpccfit, ruiccfit, reps = 1000)

######################################
### Part 2: Modeling event counts ####
######################################

#load packages
library(lme4) #version 1.1.23, published 7 April 2020
library(MASS) #version 7.3.53, published 9 September 2020
library(DHARMa) #version 0.4.5, published 16 January 2022
library(performance) #version 0.9.0, published 30 March 2022
library(sjPlot) #version 2.8.11, published 7 August 2022
library(report)

#read in dataset
dat <- read.csv("Counts.csv")
head(dat)

##subset events for day and night diel periods
datdn <- subset(dat, Diel == "Day" | Diel == "Night")

#convert categorical variables to factors
datdn$Grid <- as.factor(datdn$Grid)  
datdn$Group <- as.factor(datdn$Group)
datdn$Diel <- as.factor(datdn$Diel)
datdn$Camera <- as.factor(datdn$Camera)


#relevel factors
datdn$Grid = relevel(datdn$Grid, ref="RNP LR")
datdn$Diel = relevel(datdn$Diel, ref="Day")
datdn$Group = relevel(datdn$Group, ref="Male")

#hours as offset due to variation in sampling effort by camera station
Hours <- datdn$Hours


#fit camera station as a random intercept due to repeated measures 

#null model: count ~ 1 + 1/camera
m3 <- glmer.nb(Count ~ 1 + (1|Camera) + offset(log(Hours)), data=datdn,
               control=glmerControl(optimizer="bobyqa",
                                    optCtrl=list(maxfun=1e5)))

summary(m3)
print(m3)
BIC(m3)
AIC(m3)
logLik(m3)
plot(m3)

#null model with diel and group as random slopes and camera as random intercept (road and water are property of station so cannot be fit as random slopes)
m3.1 <- glmer.nb(Count ~ 1 + (1 + Diel + Group|Camera) + offset(log(Hours)), data=datdn,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=1e5)))

summary(m3.1)
print(m3.1)
BIC(m3.1)
AIC(m3.1)
logLik(m3.1)
plot(m3.1)

#random-intercept only models

m4 <- update(m3, . ~ . + Grid*Diel + Grid*Road + Water)
summary(m4)
AIC(m4)

m5 <- update(m3, . ~ . + Grid*Road*Group)
summary(m5)
AIC(m5)

m6 <- update(m3, . ~ . + Water*Grid*Group)
summary(m6)

m7 <- update(m3, . ~ . + Diel*Grid*Group)
summary(m7)

m8 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Group)
summary(m8)

m9 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group)
summary(m9)
anova(m8, m9)

#print HTML table
tab_model(m9)

#plot coefficients
theme_set(theme_sjplot())
plot_model(m9, show.values = TRUE, value.offset = .5, value.size = 3, dot.size = 1.5, sort.est = TRUE, title="")

## fit top random intercept model with diel and group as random slopes for confirmatory purposes
m9.1 <- update(m3.1, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group)
summary(m9.1)
plot_model(m9.1, show.values = TRUE, value.offset = .5, value.size = 3, dot.size = 1.5, sort.est = TRUE, title="")
tab_model(m9.1)

#model fit
ED <- resid(m9, type="deviance")
plot(x=datdn$Water, y=ED, main="Deviance residuals vs Water") 
plot(x=datdn$Road, y=ED, main="Deviance residuals vs Road") 
plot(x=datdn$Grid, y=ED, main="Deviance residuals vs Grid") 
plot(x=datdn$Camera, y=ED, main="Deviance residuals vs Camera") 
plot(x=datdn$Group, y=ED, main="Deviance residuals vs Group") 
plot(x=datdn$Diel, y=ED, main="Deviance residuals vs Diel")

#model diagnostics for m9
#qqplot of random effect means should fall along straight line, checks the normality and heteroscedasticity assumptions
plot_model(m9, type="diag", show.values = TRUE)

#calculate simulated residuals (randomized quantile residuals)
simulationOutput <- simulateResiduals(fittedModel = m9, plot = T)

#plot simulated residuals against each predictor variable (for factors)
plotResiduals(simulationOutput, form = datdn$Grid) 
plotResiduals(simulationOutput, form = datdn$Group) 
plotResiduals(simulationOutput, form = datdn$Water) 
plotResiduals(simulationOutput, form = datdn$Road) 
plotResiduals(simulationOutput, form = datdn$Camera) 

#model diagnostics for m9.1
#qqplot of random effect means should fall along straight line, checks the normality and heteroscedasticity assumptions
plot_model(m9.1, type="diag", show.values = TRUE)

#calculate simulated residuals (randomized quantile residuals)
simulationOutput <- simulateResiduals(fittedModel = m9.1, plot = T)

#a qq-plot to detect overall deviations from the expected distribution
plotQQunif(simulationOutput) 

#plot simulated residuals against each predictor variable (for factors)
plotResiduals(simulationOutput, form = datdn$Grid) 
plotResiduals(simulationOutput, form = datdn$Group) 
plotResiduals(simulationOutput, form = datdn$Water) 
plotResiduals(simulationOutput, form = datdn$Road) 
plotResiduals(simulationOutput, form = datdn$Camera) 


#other random intercept only models
m10 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group + Group*Road)
summary(m10)

m11 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Diel*Group + Group*Water)
summary(m11)

m12 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group + Group*Diel*Grid)
summary(m12)

m13 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group + Group*Grid*Water)
summary(m13)

m14 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group + Group*Grid*Road)
summary(m14)

m15 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group + Group*Grid*Diel*Water)
summary(m15)

m16 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Water + Diel*Group + Group*Grid*Diel*Road)
summary(m16)

m17 <- update(m3, . ~ . + Diel*Grid*Group + Road*Grid*Group + Water*Group*Grid)
summary(m17)

m18 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Diel*Group)
summary(m18)

m19 <- update(m3, . ~ . + Diel*Grid + Road*Grid + Group)
summary(m19)

m20 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Grid*Group)
summary(m20)

m21 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Road*Water)
summary(m21)
anova(m9, m21) #m21 not better

m22 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Water*Diel)
summary(m22)
anova(m9, m24) #m22 not better

m23 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Road*Diel)
summary(m23)
anova(m9, m23) #m25 not better

m24 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Grid*Water)
summary(m24)
anova(m9, m24) #not better

m25 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Grid*Water*Diel)
summary(m25)
anova(m9, m25) #not better

m26 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Grid*Road*Diel)
summary(m26)
anova(m9, m26) #m28 not better

m27 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Grid*Water*Group)
summary(m27)

m28 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Grid*Diel*Group)
summary(m28)

m29 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Diel*Water*Group)
summary(m29)

m30 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Diel*Road*Group)
summary(m30)

m31 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Diel*Water*Road)
summary(m31)

m32 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Group*Water*Road)
summary(m32)

m33 <- update(m3, . ~ . + Diel*Grid + Water + Road*Grid + Diel*Group + Grid*Water*Road)
summary(m33)

m34 <- update(m3, . ~ . + Grid*Water*Road*Group)
summary(m34)