## Chapter 6 Analysis Script
## Josephine B. Smit

#Part 1: Activity pattern figures and analyses
#Part 2: Modeling the effect of risk, sex, season, and diel period on elephant event counts
#Part 3: Group size tests

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

##https://www.rdocumentation.org/packages/overlap/versions/0.3.2
#overlap version 0.3.2, published 3 May 2018
#activity version 1.3.2, published 25 April 2022

#load the overlap and activity packages
library(overlap)
library(activity)

#read in wet season data, times already converted to 'sun time' in radians
wet <- read.csv("wetsuntime.csv")

timeST <- wet$suntime

#extract the data for each risk level
lowwet <- timeST[wet$Risk == "Low"]
highparkwet <- timeST[wet$Risk == "Highpark"]
highvillagewet <- timeST[wet$Risk == "Highvillage"]

#Comparison: all events low and highpark
overlapPlot(lowwet, highparkwet, main="")
legend('topleft', c("Low-risk park", "High-risk park"), lty=c(1,2), col=c(1,4), bty='n')

#Comparison: all events low and highvillage
overlapPlot(lowwet, highvillagewet, main="")
legend('topleft', c("Low-risk park", "High-risk village"), lty=c(1,2), col=c(1,4), bty='n')


##Read in dry season data

dry <- read.csv("drysuntime.csv")

timeST <- dry$suntime

#extract the data for all elephant events
eledry <- timeST[dry$X01_Species == "01_Elephant"]

#extract the data for cow-calf groups
cclowdry <- timeST[dry$X02_Group_type == "02_CC" & dry$Risk == "Low"]
cchighpdry <- timeST[dry$X02_Group_type == "02_CC" & dry$Risk == "Highpark"]
cchighvdry <- timeST[dry$X02_Group_type == "02_CC" & dry$Risk == "Highvillage"]

#extract the data for lone bull groups
lblowdry <- timeST[dry$X02_Group_type == "02_LB" & dry$Risk == "Low"]
lbhighpdry <- timeST[dry$X02_Group_type == "02_LB" & dry$Risk == "Highpark"]
lbhighvdry <- timeST[dry$X02_Group_type == "02_LB" & dry$Risk == "Highvillage"]

#extract the data for bull groups
bglowdry <- timeST[dry$X02_Group_type == "02_BG" & dry$Risk == "Low"]
bghighpdry <- timeST[dry$X02_Group_type == "02_BG" & dry$Risk == "Highpark"]
bghighvdry <- timeST[dry$X02_Group_type == "02_BG" & dry$Risk == "Highvillage"]

#extract the data for mixed groups
mlowdry <- timeST[dry$X02_Group_type == "02_M" & dry$Risk == "Low"]
mhighpdry <- timeST[dry$X02_Group_type == "02_M" & dry$Risk == "Highpark"]
mhighvdry <- timeST[dry$X02_Group_type == "02_M" & dry$Risk == "Highvillage"]


#extract the data for each risk level
lowdry <- timeST[dry$Risk == "Low"]
highparkdry <- timeST[dry$Risk == "Highpark"]
highvillagedry <- timeST[dry$Risk == "Highvillage"]


#Comparison: low and highpark (all elephants) dry season
overlapPlot(lowdry, highparkdry, main="")
legend('topleft', c("Low-risk park", "High-risk park"), lty=c(1,2), col=c(1,4), bty='n')

#Comparison: low and highvillage (all elephants) dry season
overlapPlot(lowdry, highvillagedry, main="")
legend('topleft', c("Low-risk park", "High-risk village"), lty=c(1,2), col=c(1,4), bty='n')

#Comparison: wet and dry season (all elephants, low risk only)
overlapPlot(lowdry, lowwet, main="Low Risk Wet vs Dry season")
legend('topleft', c("Dry", "Wet"), lty=c(1,2), col=c(1,4), bty='n')

#Comparison: wet and dry season (all elephants, high risk park only)
overlapPlot(highparkdry, highparkwet, main="High Park Wet vs Dry season")
legend('topleft', c("Dry", "Wet"), lty=c(1,2), col=c(1,4), bty='n')

#Comparison: wet and dry season (all elephants, high risk village only)
overlapPlot(highvillagedry, highvillagewet, main="High village Wet vs Dry season")
legend('topleft', c("Dry", "Wet"), lty=c(1,2), col=c(1,4), bty='n')

#Comparison: lone bull low vs high park & high village
overlapPlot(lblowdry, lbhighpdry, main="")
legend('topleft', c("Low-risk park", "High-risk park"), lty=c(1,2), col=c(1,4), bty='n')

overlapPlot(lblowdry, lbhighvdry, main="")
legend('topleft', c("Low-risk park", "High-risk village"), lty=c(1,2), col=c(1,4), bty='n')

#Comparison: bull groups low vs high park & high village
overlapPlot(bglowdry, bghighpdry, main="")
legend('topleft', c("Low-risk park", "High-risk park"), lty=c(1,2), col=c(1,4), bty='n')

overlapPlot(bglowdry, bghighvdry, main="")
legend('topleft', c("Low-risk park", "High-risk village"), lty=c(1,2), col=c(1,4), bty='n')


#Comparison: cow-calf groups low vs high park & high village

overlapPlot(cclowdry, cchighpdry, main="")
legend('topleft', c("Low-risk park", "High-risk park"), lty=c(1,2), col=c(1,4), bty='n')

overlapPlot(cclowdry, cchighvdry, main="")
legend('topleft', c("Low-risk park", "High-risk village"), lty=c(1,2), col=c(1,4), bty='n')


#Comparison: mixed groups low vs high park & high village

overlapPlot(mlowdry, mhighpdry, main="")
legend('topleft', c("Low-risk park", "High-risk park"), lty=c(1,2), col=c(1,4), bty='n')

overlapPlot(mlowdry, mhighvdry, main="")
legend('topleft', c("Low-risk park", "High-risk village"), lty=c(1,2), col=c(1,4), bty='n')


##calculate coefficient of overlapping
#if one or both of the sample sizes is below 75 events, use Dhat1
#if both sample sizes are 75 events or more, use Dhat4

#calculate coefficient of overlapping between low and highpark activity curves, dry, all elephants
min(length(lowdry), length(highparkdry))
lowhighparkdryest <- overlapEst(lowdry, highparkdry, type="Dhat4")
lowhighparkdryest
#0.6159376 

#calculate coefficient of overlapping between low and highvillage activity curves, dry, all elephants
min(length(lowdry), length(highvillagedry))
lowhighvillagedryest <- overlapEst(lowdry, highvillagedry, type="Dhat4")
lowhighvillagedryest
#0.5743572 

#calculate coefficient of overlapping between low and highpark activity curves, wet, all elephants
min(length(lowwet), length(highparkwet))
lowhighparkwetest <- overlapEst(lowwet, highparkwet, type="Dhat1")
lowhighparkwetest
#0.8801125 

#calculate coefficient of overlapping between low and highvillage activity curves, wet, all elephants
min(length(lowwet), length(highvillagewet))
lowhighvillagewetest <- overlapEst(lowwet, highvillagewet, type="Dhat1")
lowhighvillagewetest
#0.6889812 


#coefficient of overlapping between low and highpark activity curves, dry, cow-calf
min(length(cclowdry), length(cchighpdry))
cclowhighparkdryest <- overlapEst(cclowdry, cchighpdry, type="Dhat1")
cclowhighparkdryest
#0.6321221

#coefficient of overlapping between low and highvillage activity curves, dry, cow-calf
min(length(cclowdry), length(cchighvdry))
cclowhighvillagedryest <- overlapEst(cclowdry, cchighvdry, type="Dhat1")
cclowhighvillagedryest
#0.6031421

#coefficient of overlapping between low and highpark activity curves, dry, lone bull
min(length(lblowdry), length(lbhighpdry))
lblowhighparkdryest <- overlapEst(lblowdry, lbhighpdry, type="Dhat1")
lblowhighparkdryest
#0.6086988 

#coefficient of overlapping between low and highvillage activity curves, dry, lone bull
min(length(lblowdry), length(lbhighvdry))
lblowhighvillagedryest <- overlapEst(lblowdry, lbhighvdry, type="Dhat1")
lblowhighvillagedryest
#0.4912552 

#coefficient of overlapping between low and highpark activity curves, dry, bull group
min(length(bglowdry), length(bghighpdry))
bglowhighparkdryest <- overlapEst(bglowdry, bghighpdry, type="Dhat1")
bglowhighparkdryest
#0.7272095

#coefficient of overlapping between low and highvillage activity curves, dry, bull group
min(length(bglowdry), length(bghighvdry))
bglowhighvillagedryest <- overlapEst(bglowdry, bghighvdry, type="Dhat1")
bglowhighvillagedryest
#0.5387675

#coefficient of overlapping between low and highpark activity curves, dry, mixed group
min(length(mlowdry), length(mhighpdry))
mlowhighparkdryest <- overlapEst(mlowdry, mhighpdry, type="Dhat1")
mlowhighparkdryest
#0.5137704 

#coefficient of overlapping between low and highvillage activity curves, dry, mixed group
min(length(mlowdry), length(mhighvdry))
mlowhighvillagedryest <- overlapEst(mlowdry, mhighvdry, type="Dhat1")
mlowhighvillagedryest
#0.5313888

#coefficient of overlap between wet and dry season, low risk, all elephants
min(length(lowwet), length(lowdry)) 
lowdrywetest <- overlapEst(lowdry, lowwet, type="Dhat4")
lowdrywetest
#0.8941837 

#coefficient of overlapping between wet and dry season, highvillage, all elephants
min(length(highvillagedry), length(highvillagewet))
highvillagedrywetest <- overlapEst(highvillagedry, highvillagewet, type="Dhat1")
highvillagedrywetest
#0.7131617

#coefficient of overlapping between wet and dry season, highpark, all elephants
min(length(highparkdry), length(highparkwet))
highparkdrywetest <- overlapEst(highparkdry, highparkwet, type="Dhat1")
highparkdrywetest
#0.61992 

##calculate confidence intervals for coefficients of overlapping

#wet vs dry, low risk, all elephants
#Resample (smoothed bootstrap samples)
lowwetboot <- resample(lowwet, 10000)
dim(lowwetboot)

lowdryboot <- resample(lowdry, 10000)
dim(lowdryboot)

#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.

#Calculate bootstrapped mean (takes a bit of time)
lowwetdry <- bootEst(lowwetboot, lowdryboot, type="Dhat4") 
lowwetdry

#generate CI
bootCI(lowdrywetest, lowwetdry)
#basic0 0.8356586 0.9464675


#low vs highpark, dry season, all elephants
#Resample (smoothed bootstrap samples)
highparkdryboot <- resample(highparkdry, 10000)
dim(highparkdryboot)

lowdryboot <- resample(lowdry, 10000)
dim(lowdryboot)

#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.

#Calculate bootstrapped mean (takes a bit of time)
lowhighparkdry <- bootEst(lowdryboot, highparkdryboot, type="Dhat4") 
lowhighparkdry

#generate CI
bootCI(lowhighparkdryest, lowhighparkdry)
#basic0  0.5417343 0.6920289


#low vs high village, dry season, all elephants
#Resample (smoothed bootstrap samples)
highvillagedryboot <- resample(highvillagedry, 10000)
dim(highvillagedryboot)

lowdryboot <- resample(lowdry, 10000)
dim(lowdryboot)

#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.

#Calculate bootstrapped mean (takes a bit of time)
lowhighvillagedry <- bootEst(lowdryboot, highvillagedryboot, type="Dhat4") 
lowhighvillagedry

#generate CI
bootCI(lowhighvillagedryest, lowhighvillagedry)
#basic0 0.5232067 0.6254414


#low vs high park, wet season, all elephants
#Resample (smoothed bootstrap samples)
lowwetboot <- resample(lowwet, 10000)
dim(lowwetboot)

highparkwetboot <- resample(highparkwet, 10000)
dim(highparkwetboot)

#Calculate bootstrapped mean (takes a bit of time)
lowhighparkwet <- bootEst(lowwetboot, highparkwetboot, type="Dhat4") 
lowhighparkwet

#generate CI
bootCI(lowhighparkwetest, lowhighparkwet)
#basic0 0.7903913 0.9552299


#low bs high village, wet season, all elephants
#Resample (smoothed bootstrap samples)
lowwetboot <- resample(lowwet, 10000)
dim(lowwetboot)

highvillagewetboot <- resample(highvillagewet, 10000)
dim(highvillagewetboot)

#Calculate bootstrapped mean (takes a bit of time)
lowhighvillagewet <- bootEst(lowwetboot, highvillagewetboot, type="Dhat4") 
lowhighvillagewet

#generate CI
bootCI(lowhighvillagewetest, lowhighvillagewet)
#basic0 0.5758664 0.7921402


#CC low vs highpark, dry season
#Resample (smoothed bootstrap samples)
cchighparkdryboot <- resample(cchighpdry, 10000)
dim(cchighparkdryboot)

cclowdryboot <- resample(cclowdry, 10000)
dim(cclowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
cclowhighparkdry <- bootEst(cclowdryboot, cchighparkdryboot, type="Dhat1") 

#generate CI
bootCI(cclowhighparkdryest, cclowhighparkdry)
#basic0 0.5240949 0.7387547


#CC low vs highvillage, dry season
#Resample (smoothed bootstrap samples)
cchighvillagedryboot <- resample(cchighvdry, 10000)
dim(cchighvillagedryboot)

cclowdryboot <- resample(cclowdry, 10000)
dim(cclowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
cclowhighvillagedry <- bootEst(cclowdryboot, cchighvillagedryboot, type="Dhat4") 

#generate CI
bootCI(cclowhighvillagedryest, cclowhighvillagedry)
#basic0 0.5033231 0.6986600


#LB low vs highpark, dry season
#Resample (smoothed bootstrap samples)
lbhighparkdryboot <- resample(lbhighpdry, 10000)
dim(lbhighparkdryboot)

lblowdryboot <- resample(lblowdry, 10000)
dim(lblowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
lblowhighparkdry <- bootEst(lblowdryboot, lbhighparkdryboot, type="Dhat1") 

#generate CI
bootCI(lblowhighparkdryest, lblowhighparkdry)
#basic0 0.4471988 0.7681741


#LB low vs highvillage, dry season
#Resample (smoothed bootstrap samples)
lbhighvillagedryboot <- resample(lbhighvdry, 10000)
dim(lbhighvillagedryboot)

lblowdryboot <- resample(lblowdry, 10000)
dim(lblowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
lblowhighvillagedry <- bootEst(lblowdryboot, lbhighvillagedryboot, type="Dhat1") 

#generate CI
bootCI(lblowhighvillagedryest, lblowhighvillagedry)
#basic0 0.3750924 0.6086714


#BG low vs highpark, dry season
#Resample (smoothed bootstrap samples)
bghighparkdryboot <- resample(bghighpdry, 10000)
dim(bghighparkdryboot)

bglowdryboot <- resample(bglowdry, 10000)
dim(bglowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
bglowhighparkdry <- bootEst(bglowdryboot, bghighparkdryboot, type="Dhat1") 

#generate CI
bootCI(bglowhighparkdryest, bglowhighparkdry)
#basic0 0.5157146 0.9108317


#BG low vs highvillage, dry season
#Resample (smoothed bootstrap samples)
bghighvillagedryboot <- resample(bghighvdry, 10000)
dim(bghighvillagedryboot)

bglowdryboot <- resample(bglowdry, 10000)
dim(bglowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
bglowhighvillagedry <- bootEst(bglowdryboot, bghighvillagedryboot, type="Dhat1") 

#generate CI
bootCI(bglowhighvillagedryest, bglowhighvillagedry)
#basic0 0.3597479 0.7029986


#M low vs highpark, dry season
#Resample (smoothed bootstrap samples)
mhighparkdryboot <- resample(mhighpdry, 10000)
dim(mhighparkdryboot)

mlowdryboot <- resample(mlowdry, 10000)
dim(mlowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
mlowhighparkdry <- bootEst(mlowdryboot, mhighparkdryboot, type="Dhat1") 

#generate CI
bootCI(mlowhighparkdryest, mlowhighparkdry)
#basic0 0.3452578 0.6742305


#M low vs highvillage, dry season
#Resample (smoothed bootstrap samples)
mhighvillagedryboot <- resample(mhighvdry, 10000)
dim(mhighvillagedryboot)

mlowdryboot <- resample(mlowdry, 10000)
dim(mlowdryboot)

#Calculate bootstrapped mean (takes a bit of time)
mlowhighvillagedry <- bootEst(mlowdryboot, mhighvillagedryboot, type="Dhat1") 

#generate CI
bootCI(mlowhighvillagedryest, mlowhighvillagedry)
#basic0 0.3733463 0.6843581


#Highpark, wet vs dry season
#Resample (smoothed bootstrap samples)
highparkdryboot <- resample(highparkdry, 10000)
dim(highparkdryboot)

highparkwetboot <- resample(highparkwet, 10000)
dim(highparkwetboot)

#Calculate bootstrapped mean (takes a bit of time)
highparkwetdry <- bootEst(highparkdryboot, highparkwetboot, type="Dhat1") 

#generate CI
bootCI(highparkdrywetest, highparkwetdry)
#basic0 0.5073883 0.7275011

#Highvillage, wet vs dry season
#Resample (smoothed bootstrap samples)
highvillagedryboot <- resample(highvillagedry, 10000)
dim(highvillagedryboot)

highvillagewetboot <- resample(highvillagewet, 10000)
dim(highvillagewetboot)

#Calculate bootstrapped mean (takes a bit of time)
highvillagewetdry <- bootEst(highvillagedryboot, highvillagewetboot, type="Dhat1") 

#generate CI
bootCI(highvillagedrywetest, highvillagewetdry)
#basic0 0.6056527 0.8095193


### Use Activity package to determine whether activity curves are significantly different

#dry season
timeST <- dry$suntime
alowdry <- timeST[dry$Risk == "Low"]
ahighpdry <- timeST[dry$Risk == "Highpark"]
ahighvdry <- timeST[dry$Risk == "Highvillage"]

#wet season
timeST <- wet$suntime
alowwet <- timeST[wet$Risk == "Low"]
ahighpwet <- timeST[wet$Risk == "Highpark"]
ahighvwet <- timeST[wet$Risk == "Highvillage"]


#extract the data for each group type
acclowdry <- timeST[dry$Risk == "Low" & dry$X02_Group_type == "02_CC"]
acchighpdry <- timeST[dry$Risk == "Highpark" & dry$X02_Group_type == "02_CC"]
acchighvdry <- timeST[dry$Risk == "Highvillage" & dry$X02_Group_type == "02_CC"]

alblowdry <- timeST[dry$Risk == "Low" & dry$X02_Group_type == "02_LB"]
albhighpdry <- timeST[dry$Risk == "Highpark" & dry$X02_Group_type == "02_LB"]
albhighvdry <- timeST[dry$Risk == "Highvillage" & dry$X02_Group_type == "02_LB"]

abglowdry <- timeST[dry$Risk == "Low" & dry$X02_Group_type == "02_BG"]
abghighpdry <- timeST[dry$Risk == "Highpark" & dry$X02_Group_type == "02_BG"]
abghighvdry <- timeST[dry$Risk == "Highvillage" & dry$X02_Group_type == "02_BG"]

amlowdry <- timeST[dry$Risk == "Low" & dry$X02_Group_type == "02_M"]
amhighpdry <- timeST[dry$Risk == "Highpark" & dry$X02_Group_type == "02_M"]
amhighvdry <- timeST[dry$Risk == "Highvillage" & dry$X02_Group_type == "02_M"]


#Fit and plot kernel density curve to radian time-of-day data  

alowdryfit <- fitact(alowdry)
plot(alowdryfit, main="All elephants Low Dry")

ahighpdryfit <- fitact(ahighpdry)
plot(ahighpdryfit, main="All elephants Highpark Dry")

ahighvdryfit <- fitact(ahighvdry)
plot(ahighvdryfit, main="All elephants Highvillage Dry")

acclowdryfit <- fitact(acclowdry)
plot(acclowdryfit, main="CC elephants Low Dry")

acchighpdryfit <- fitact(acchighpdry)
plot(acchighpdryfit, main="CC elephants Highpark Dry")

acchighvdryfit <- fitact(acchighvdry)
plot(acchighvdryfit, main="CC elephants Highvillage Dry")

alblowdryfit <- fitact(alblowdry)
plot(alblowdryfit, main="LB elephants Low Dry")

albhighpdryfit <- fitact(albhighpdry)
plot(albhighpdryfit, main="LB elephants Highpark Dry")

albhighvdryfit <- fitact(albhighvdry)
plot(albhighvdryfit, main="LB elephants Highvillage Dry")

abglowdryfit <- fitact(abglowdry)
plot(abglowdryfit, main="BG elephants Low Dry")

abghighpdryfit <- fitact(abghighpdry)
plot(abghighpdryfit, main="BG elephants Highpark Dry")

abghighvdryfit <- fitact(abghighvdry)
plot(abghighvdryfit, main="BG elephants Highvillage Dry")

amlowdryfit <- fitact(amlowdry)
plot(amlowdryfit, main="M elephants Low Dry")

amhighpdryfit <- fitact(amhighpdry)
plot(amhighpdryfit, main="M elephants Highpark Dry")

amhighvdryfit <- fitact(amhighvdry)
plot(amhighvdryfit, main="M elephants Highvillage Dry")

alowwetfit <- fitact(alowwet)
plot(alowwetfit, main="All elephants Low Wet")

ahighpwetfit <- fitact(ahighpwet)
plot(ahighpwetfit, main="All elephants Highpark Wet")

ahighvwetfit <- fitact(ahighvwet)
plot(ahighvwetfit, main="All elephants Highvillage Wet")

#Conduct randomisation test for the probability that two sets of circular observations come from the same distribution.
#Output: A named 4-element vector: obs = observed overlap index; null = mean null overlap index; seNull =standard error of the null distribution; pNull = probability observed index arose by chance

#calculate whether there is a statistically significant difference between the low and high park dry season curves
compareCkern(alowdryfit, ahighpdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.61583187 0.90956509 0.02337311 0.00000000 
#obs = observed coefficient of overlapping
#pNull = probability that the observed coefficient of overlapping occurred by chance

#low dry and high village
compareCkern(alowdryfit, ahighvdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.57417551 0.93446485 0.01764156 0.00000000 
#obs = observed coefficient of overlapping
#pNull = probability that the observed coefficient of overlapping occurred by chance


#CC low dry and high park
compareCkern(acclowdryfit, acchighpdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.6286585 0.8777427 0.0328042 0.0000000

#CC low dry and high village
compareCkern(acclowdryfit, acchighvdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.60444360 0.87105774 0.03004453 0.00000000 

#LB low dry and high park
compareCkern(alblowdryfit, albhighpdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.58835770 0.80915583 0.05689589 0.00100000 

#LB low dry and high village
compareCkern(alblowdryfit, albhighvdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.48535215 0.86038290 0.04366492 0.00000000 


#BG low dry and high park
compareCkern(abglowdryfit, abghighpdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.74997509 0.71001530 0.09208705 0.65000000

#BG low dry and high village
compareCkern(abglowdryfit, abghighvdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.53174343 0.82208090 0.06079477 0.00100000


#M low dry and high park
compareCkern(amlowdryfit, amhighpdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.51045075 0.79088262 0.06257136 0.00000000 

#M low dry and high village
compareCkern(amlowdryfit, amhighvdryfit, reps = 1000)
#obs       null     seNull      pNull 
#0.52520571 0.80812978 0.06149135 0.00000000 

#low, wet vs dry season
compareCkern(alowdryfit, alowwetfit, reps = 1000)
#obs       null     seNull      pNull 
#0.89413723 0.93132378 0.01802315 0.04000000
#obs = observed coefficient of overlapping
#pNull = probability that the observed coefficient of overlapping occurred by chance

#highpark, wet vs dry season
compareCkern(ahighpdryfit, ahighpwetfit, reps = 1000)
#obs       null     seNull      pNull 
#0.61631467 0.88003633 0.03574822 0.00000000 

#highvillage, wet vs dry season
compareCkern(ahighvdryfit, ahighvwetfit, reps = 1000)
#obs       null     seNull      pNull 
#0.72404741 0.82827418 0.04499953 0.01400000 


#wet season: low vs highpark
compareCkern(alowwetfit, ahighpwetfit, reps = 1000)
#obs       null     seNull      pNull 
#0.88298008 0.89168490 0.03228321 0.35700000

#wet season: low vs highvillage
compareCkern(alowwetfit, ahighvwetfit, reps = 1000)
#obs       null     seNull      pNull 
#0.68770158 0.85848796 0.04155448 0.00100000 


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

##GLMM to model elephant event counts at water sources as a function of: risk (low-risk Park, high-risk Park, high-risk village), season(dry/wet), sex (male/female) and diel period (night/day)

#Response variable is Count = number of camera trap events detected, summed by camera trap station, season, risk and group type
#Station is specified as a random intercept term to account for the fact that I have repeated measures from camera trap stations
#Sex and season are specified as random slopes (risk is not, as risk is a property of station)
#Used Negative Binomial distribution because of overdispersion of Count data

#load packages
library(lme4) #version 1.1.29, published 7 April 2022
library(DHARMa) #version 0.4.5, published 16 January 2022
library(MASS) #version 7.3.56, published 23 March 2022
library(performance) #version 0.9.0, published 30 March 2022
library(sjPlot) #version 2.8.10, published 26 November 2021
library(report) #version 0.5.1, published 22 February 2022

#read in dataset of event counts by station, diel period, sex, and season
dat <- read.csv("Counts.csv")
head(dat)


#convert categorical variables to factors
dat$Risk <- as.factor(dat$Risk)  
dat$Season <- as.factor(dat$Season)
dat$Sex <- as.factor(dat$Sex)
dat$Diel <- as.factor(dat$Diel)
dat$Station <- as.factor(dat$Station)

#relevel factors 
dat$Risk = relevel(dat$Risk, ref="Low")
dat$Diel = relevel(dat$Diel, ref="Day")
dat$Sex = relevel(dat$Sex, ref="Male")
dat$Season = relevel(dat$Season, ref="Dry")

#offset variable: sampling effort per diel period per camera station in hours to account for differences in sampling effort between camera stations and diel period
Hours <- dat$Hours
summary(Hours)

#null model
m0.1 <- glmer.nb(Count ~ 1 + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m0.1)
print(m0.1)
logLik(m0.1)
AICc(m0.1)
#AICc 790.0


#Global model 
m1.1 <- glmer.nb(Count ~ Risk*Diel*Sex + Risk*Diel*Season + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.1)
AICc(m1.1) #789.5
logLik(m1.1)


##drop 3-way interaction with sex
m1.2 <- glmer.nb(Count ~ Risk*Diel*Season + Risk*Sex + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.2)
AICc(m1.2) #781.2
logLik(m1.2)


##drop 2-way interaction with sex
m1.3 <- glmer.nb(Count ~ Risk*Diel*Season + Sex + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.3)
AICc(m1.3) #780.8
logLik(m1.3)


##drop sex 
m1.4 <- glmer.nb(Count ~ Risk*Diel*Season + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.4)
AICc(m1.4) #778.4
logLik(m1.4)


##2-way interactions 
m1.5 <- glmer.nb(Count ~ Risk*Diel + Risk*Season + Risk*Sex + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.5)
AICc(m1.5) #778.2


##2-way interactions minus 2-way interaction with sex
m1.6 <- glmer.nb(Count ~ Risk*Diel + Risk*Season + Sex + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.6)
AICc(m1.6) #777.9


##2-way interactions minus sex
m1.7 <- glmer.nb(Count ~ Risk*Diel + Risk*Season + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.7) #775.7
AICc(m1.7)


##risk*diel interaction + season 
m1.8 <- glmer.nb(Count ~ Risk*Diel + Season + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                 control=glmerControl(optimizer="bobyqa",
                                      optCtrl=list(maxfun=3e5)))
summary(m1.8) 
AICc(m1.8) #780.7


##risk*season interaction + diel
m1.9 <- glmer.nb(Count ~ Risk*Season + Diel + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                  control=glmerControl(optimizer="bobyqa",
                                       optCtrl=list(maxfun=3e5)))
summary(m1.9) 
AICc(m1.9) #778.2

##risk*season interaction + diel
m1.10 <- glmer.nb(Count ~ Risk*Season + Diel + Sex + (1 + Diel + Sex + Season |Station) + offset(log(Hours)), data=dat,
                  control=glmerControl(optimizer="bobyqa",
                                       optCtrl=list(maxfun=3e5)))
summary(m1.10) 
AICc(m1.10) #780.4


#model comparisons
anova(m1.1, m1.7) #m1.1 (global model) is not better than m1.7 (lowest AICc)
anova(m1.7, m1.6) #m1.7 is better than m1.6
anova(m1.7, m1.8) #m1.7 is better than simpler m1.8
anova(m1.7, m1.9) #m1.7 is better than simpler m1.9
anova(m1.7, m1.10) #m1.7 is better than simpler m1.10

#print plots and table for top-ranked model
tab_model(m1.7)

theme_set(theme_sjplot())
plot_model(m1.7, show.values = TRUE, value.offset = .5, value.size = 3, dot.size = 1.5, sort.est = TRUE, title="")

#marginal effects plot
plot_model(m1.7, type = "pred", terms = c("Risk", "Diel", "Season"))

#model diagnostics
plot_model(m1.7, type="diag", show.values = TRUE)
check_model(m1.7)
check_overdispersion(m1.7) #No overdispersion detected

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

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

plotResiduals(simulationOutput) 

#plot simulated residuals against each predictor variable (for factors)
plotResiduals(simulationOutput, form = dat$Risk) #NS
plotResiduals(simulationOutput, form = dat$Diel) #NS
plotResiduals(simulationOutput, form = dat$Sex) #NS
plotResiduals(simulationOutput, form = dat$Season) #NS
plotResiduals(simulationOutput, form = dat$Station) #NS

report(m1.7)

####################################
## Part 3: Group size tests #######
####################################

###Group size figures and Kolmogorov-Smirnov tests

library("dplyr"), published 28 April 2022 

#read in group size data
all <- read.csv(file = 'All_Size.csv')


### Wet season tests ####

###CC low vs high risk

ccwetlow <- all[all$Season == "Wet" & all$Group == "CC" & all$Risk == "Low",]
ccwethp <- all[all$Season == "Wet" & all$Group == "CC" & all$Risk == "Highpark",]

ccwetlow <- ccwetlow[,4]
ccwethp <- ccwethp[,4]
ccwetlow<- as.numeric(ccwetlow)
ccwethp<- as.numeric(ccwethp)


ks.test(ccwetlow, ccwethp)
#Two-sample Kolmogorov-Smirnov test

#data:  ccwetlow and ccwethp
#D = 0.26778, p-value = 0.1084



#bull groups low vs high risk (excluding lone bulls)

bwetlow <- all[all$Season == "Wet" & all$Group == "BG" & all$Risk == "Low",]
bwethp <- all[all$Season == "Wet" & all$Group == "BG" & all$Risk == "Highpark",]
bwethv <- all[all$Season == "Wet" & all$Group == "BG" & all$Risk == "Highvillage",]

bwetlow <- bwetlow[,4]
bwethp <- bwethp[,4]
bwethv <- bwethv[,4]
bwetlow<- as.numeric(bwetlow)
bwethp<- as.numeric(bwethp)
bwethv<- as.numeric(bwethv)

ks.test(bwetlow, bwethp)
#Two-sample Kolmogorov-Smirnov test

#data:  bwetlow and bwethp
#D = 0.18919, p-value = 0.7112
#alternative hypothesis: two-sided

ks.test(bwetlow, bwethv)
#data:  bwetlow and bwethv
#D = 0.56892, p-value = 3.109e-05


#mixed groups low vs high risk

mwetlow <- all[all$Season == "Wet" & all$Group == "M" & all$Risk == "Low",]
mwethp <- all[all$Season == "Wet" & all$Group == "M" & all$Risk == "Highpark",]

mwetlow <- mwetlow[,4]
mwethp <- mwethp[,4]
mwetlow<- as.numeric(mwetlow)
mwethp<- as.numeric(mwethp)


ks.test(mwetlow, mwethp)
#Two-sample Kolmogorov-Smirnov test

#data:  mwetlow and mwethp
#D = 0.15588, p-value = 0.8915
#alternative hypothesis: two-sided



### Dry season tests ####

###CC low vs high risk

ccdrylow <- all[all$Season == "Dry" & all$Group == "CC" & all$Risk == "Low",]
ccdryhp <- all[all$Season == "Dry" & all$Group == "CC" & all$Risk == "Highpark",]
ccdryhv <- all[all$Season == "Dry" & all$Group == "CC" & all$Risk == "Highvillage",]

ccdrylow <- ccdrylow[,4]
ccdryhp <- ccdryhp[,4]
ccdryhv <- ccdryhv[,4]
ccdrylow<- as.numeric(ccdrylow)
ccdryhp<- as.numeric(ccdryhp)
ccdryhv<- as.numeric(ccdryhv)

ks.test(ccdrylow, ccdryhp)
#Two-sample Kolmogorov-Smirnov test

#D = 0.11798, p-value = 0.5157


ks.test(ccdrylow, ccdryhv)
#Two-sample Kolmogorov-Smirnov test

#D = 0.13889, p-value = 0.187


summary(ccdrylow)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#1.000   3.000   5.000   5.902   7.250  35.000 

summary(ccdryhp)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#1.000   3.000   5.000   5.046   6.000  12.000 

summary(ccdryhv)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#1.000   4.000   6.000   6.667   9.000  18.000 

###Bull groups low vs high risk (excluding lone bulls)

bdrylow <- all[all$Season == "Dry" & all$Group == "BG" & all$Risk == "Low",]
bdryhp <- all[all$Season == "Dry" & all$Group == "BG" & all$Risk == "Highpark",]
bdryhv <- all[all$Season == "Dry" & all$Group == "BG" & all$Risk == "Highvillage",]

bdrylow <- bdrylow[,4]
bdryhp <- bdryhp[,4]
bdryhv <- bdryhv[,4]
bdrylow<- as.numeric(bdrylow)
bdryhp<- as.numeric(bdryhp)
bdryhv<- as.numeric(bdryhv)

ks.test(bdrylow, bdryhp)
#Two-sample Kolmogorov-Smirnov test

#data:  bdrylow and bdryhp
#D = 0.15789, p-value = 0.5835

ks.test(bdrylow, bdryhv)
#Two-sample Kolmogorov-Smirnov test

#data:  bdrylow and bdryhv
#D = 0.57727, p-value = 2.295e-06

summary(bdrylow)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#2.000   2.000   2.000   2.842   3.000   6.000

summary(bdryhp)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#2.000   2.000   2.000   2.462   3.000   3.000 

summary(bdryhv)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#2.00    3.00    6.00    8.32   11.00   36.00 

#mixed groups low vs high risk

mdrylow <- all[all$Season == "Dry" & all$Group == "M" & all$Risk == "Low",]
mdryhp <- all[all$Season == "Dry" & all$Group == "M" & all$Risk == "Highpark",]
mdryhv <- all[all$Season == "Dry" & all$Group == "M" & all$Risk == "Highvillage",]

mdrylow <- mdrylow[,4]
mdryhp <- mdryhp[,4]
mdryhv <- mdryhv[,4]
mdrylow<- as.numeric(mdrylow)
mdryhp<- as.numeric(mdryhp)
mdryhv<- as.numeric(mdryhv)

ks.test(mdrylow, mdryhp)
#Two-sample Kolmogorov-Smirnov test

#D = 0.19908, p-value = 0.3209
#alternative hypothesis: two-sided

ks.test(mdrylow, mdryhv)
#Two-sample Kolmogorov-Smirnov test

#D = 0.20276, p-value = 0.3317

summary(mdrylow)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#3.00    5.00   10.00   11.71   16.00   26.00 

summary(mdryhp)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#2.00    5.00    8.00   11.39   12.50   62.00 

summary(mdryhv)
#Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#2.00    7.00    9.00   10.58   12.00   48.00 


#seasonal comparisons for low risk site only

#CC
ks.test(ccwetlow, ccdrylow)
#Two-sample Kolmogorov-Smirnov test

#data:  ccwetlow and ccdrylow
#D = 0.17, p-value = 0.04728

#BG (excluding LB)
ks.test(bwetlow, bdrylow)
#Two-sample Kolmogorov-Smirnov test

#data:  bwetlow and bdrylow
#D = D = 0.076814, p-value = 0.8755

#M
ks.test(mwetlow, mdrylow)
#Two-sample Kolmogorov-Smirnov test

#data:  mwetlow and mdrylow
#D = 0.16807, p-value = 0.3408