## Code for statistical analysis in 
## 'Lee & Cowlishaw. 2017. Switching spatial scale reveals dominance-dependent social foraging tactics in a wild primate.' 

## load packages
library(lme4)
library(colorspace)

## inverse logit function (for figs)
logit <- function(x) 1/(1+(exp(-x)))

## load required data files (available in Dryad)
AllData <- read.csv("ForagingDecisions_PJ.csv", header=T, sep=";")
DiffData <- read.csv("DominanceDifferences_PJ.csv", header = T, sep=";")
CompEx <- read.csv("CompetitiveExclusion_PJ.csv", header=T, sep=";")

###########################################################################################
## Analysis 1 -- Dominance and social foraging decisions at different spatial scales
## uses AllData
###########################################################################################

## Variable explanations for those not clear:
## ExploitNumeric -- binary indicator of produce (0) or join (1)
## SpatialScale -- patch (A) or sub-patch (B)

## make categorical variables non-numeric
AllData$Year <- as.factor(AllData$Year)
AllData$AgeClassBinary <- factor(c("J","A")[AllData$AgeClassBinary])

## build maximal model
All_Model_Max <- glmer(ExploitNumeric ~ SpatialScale * DominanceRank * AgeClassBinary + Troop + Year + (1|FocalID) + (1|FollowID), data=AllData, family = binomial, control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e8)))
summary(All_Model_Max)

All_Model_2 <- glmer(ExploitNumeric ~ SpatialScale * DominanceRank + DominanceRank:AgeClassBinary + SpatialScale:AgeClassBinary + AgeClassBinary + Troop + Year + (1|FocalID) + (1|FollowID), data=AllData, family = binomial, control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
summary(All_Model_2) ## final model

## LR tests
drop1(All_Model_Max, ~ SpatialScale:DominanceRank:AgeClassBinary, test="Chisq")
drop1(All_Model_2, test="Chisq")

## so keep All_Model_2

## Figure 1
newdat_All <- expand.grid(
  SpatialScale = c(0,1)
  ,DominanceRank = seq(0,1,0.01)
  ,AgeClassBinary = c(0,1)
  ,Troop = c(0,1)
  , Year = c(0,1)
  ,ExploitNumeric=0
)

mm_All <- model.matrix(terms(All_Model_2),newdat_All)
newdat_All$ExploitNumeric <- mm_All %*% fixef(All_Model_2)
pvar1_All <- diag(mm_All %*% tcrossprod(vcov(All_Model_2),mm_All))
newdat_All <- data.frame(
  newdat_All
  , plo = newdat_All$ExploitNumeric-2*sqrt(pvar1_All)
  , phi = newdat_All$ExploitNumeric+2*sqrt(pvar1_All)
)

AllMean <- tapply(newdat_All$ExploitNumeric, as.factor(newdat_All$SpatialScale):as.factor(newdat_All$AgeClassBinary):as.factor(newdat_All$DominanceRank), mean)
AllHi <- tapply(newdat_All$plo, as.factor(newdat_All$SpatialScale):as.factor(newdat_All$AgeClassBinary):as.factor(newdat_All$DominanceRank), mean)
AllLo <- tapply(newdat_All$phi, as.factor(newdat_All$SpatialScale):as.factor(newdat_All$AgeClassBinary):as.factor(newdat_All$DominanceRank), mean)

## for 95% CI polygons
All_P_Adult_x <- c(seq(0,1,0.01),rev(seq(0,1,0.01)))
All_P_Adult_y <- c(logit(AllLo[1:101]),rev(logit(AllHi[1:101])))

All_P_Juv_x <- c(seq(0,1,0.01),rev(seq(0,1,0.01)))
All_P_Juv_y <- c(logit(AllLo[102:202]),rev(logit(AllHi[102:202])))

All_SP_Adult_x <- c(seq(0,1,0.01),rev(seq(0,1,0.01)))
All_SP_Adult_y <- c(logit(AllLo[203:303]),rev(logit(AllHi[203:303])))

All_SP_Juv_x <- c(seq(0,1,0.01),rev(seq(0,1,0.01)))
All_SP_Juv_y <- c(logit(AllLo[304:404]),rev(logit(AllHi[304:404])))

## plot
par(mfrow=c(1,2))
plot(AllData$ExploitNumeric ~ AllData$DominanceRank, type='n', bty='n', xlab="Dominance rank", ylab="Probability of joining", ylim=c(0.3,1), cex.lab=1.4, cex.axis=1.4)
polygon(All_P_Adult_x,All_P_Adult_y,col=rainbow_hcl(10, alpha=0.8)[5], border=NA)
polygon(All_P_Juv_x,All_P_Juv_y,col=rainbow_hcl(10, alpha=0.8)[9], border=NA)
lines(seq(0,1,0.01),logit(AllMean[1:101]))#patch & adult
lines(seq(0,1,0.01),logit(AllMean[102:202]), lty=2)# patch & juve
text(0.05,1,"a", cex=1.4)

plot(AllData$ExploitNumeric ~ AllData$DominanceRank, type='n', bty='n', xlab="Dominance rank", ylab="", ylim=c(0,0.2), cex.lab=1.4, cex.axis=1.4)
polygon(All_SP_Adult_x,All_SP_Adult_y,col=rainbow_hcl(10, alpha=0.4)[5], border=NA)
polygon(All_SP_Juv_x,All_SP_Juv_y,col=rainbow_hcl(10, alpha=0.4)[9], border=NA)
lines(seq(0,1,0.01),logit(AllMean[203:303]), lty=1)# subpatch & adult
lines(seq(0,1,0.01),logit(AllMean[304:404]), lty=2)# subpatch & juve
text(0.05,0.2,"b", cex=1.4)

###########################################################################################
## Analysis 2 -- Social constraints on joining behaviour at different spatial scales
## uses DiffData
###########################################################################################

## Variable explanations for those not clear:
## SpatialScale -- patch (A) or sub-patch (B)

## make categorical variables non-numeric
DiffData$AgeClassBinary <- factor(c("J","A")[DiffData$AgeClassBinary])

## observed values are DiffMeanRank, calculated as the rank of the joiner minus the mean rank of the joined individual(s)

## to generate random distributions, need to randomly resample mean rank of joined individual from all possible values
## so make a sample to randomly draw from, representing all possible ranks (including own, as could scrounge from one +1 rank and one -1 rank) for each troop year combo

## import ranks for each year-troop combo for random sampling
DomRanks <- read.csv("Ranks_for_rand_sampling_PJ.csv", header=T, sep=";")

## generate random distributions
DiffData$TroopYear <- as.numeric(as.factor(paste(DiffData$Troop, DiffData$Year)))
nsims <- 10000

## function for sampling that ignores NAs
sample_dom <- function(x){
  sample(x[!is.na(x)], size = 1)
}

## 10000 reps for each observation, drawn randomly from dominance hierarchy for each troop-year combo (as hierarchies different)
Rand_dom <- replicate(10000,sapply(DiffData$TroopYear, function(x) sample_dom(DomRanks[,x])))

## Rank difference between observed joiner and each randomly selected joined individual
Rand_dom_diff <- DiffData$DominanceRank-Rand_dom

## now generate random distributions for proportion of obs <=0 for rank diff for each spatial scale / age class combination
## separately for different scale-age combos
Rand_dom_diff.prop.negatives.AP <- apply(Rand_dom_diff[DiffData$AgeClassBinary=="A" & DiffData$SpatialScale=="A",], 2, function(x) sum(x<=0)/length(x))
Rand_dom_diff.prop.negatives.JP <- apply(Rand_dom_diff[DiffData$AgeClassBinary=="J" & DiffData$SpatialScale=="A",], 2, function(x) sum(x<=0)/length(x))
Rand_dom_diff.prop.negatives.ASP <- apply(Rand_dom_diff[DiffData$AgeClassBinary=="A" & DiffData$SpatialScale=="B",], 2, function(x) sum(x<=0)/length(x))
Rand_dom_diff.prop.negatives.JSP <- apply(Rand_dom_diff[DiffData$AgeClassBinary=="J" & DiffData$SpatialScale=="B",], 2, function(x) sum(x<=0)/length(x))

## Figure 2
## plot distributions, with 95% TIs and observed values. Sample sizes at top of each distribution
par(mfrow=c(2,2),mar=c(4.5,4.5,2,1))
## adult patch
plot(density(Rand_dom_diff.prop.negatives.AP), xlab='', main="N = 313", xlim=c(0.3,0.5), bty="n", cex.lab=1.8, cex.axis=1.6, type='n')
polygon(density(Rand_dom_diff.prop.negatives.AP), col=rainbow_hcl(10, alpha=0.8)[5], lwd=1.5)
abline(v=sort(Rand_dom_diff.prop.negatives.AP)[length(Rand_dom_diff.prop.negatives.AP)*0.975], col="black", lty=2, lwd=1.5) # upper CI
abline(v=sort(Rand_dom_diff.prop.negatives.AP)[length(Rand_dom_diff.prop.negatives.AP)*0.025], col="black", lty=2, lwd=1.5) # lower CI
abline(v=sum(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="A" & DiffData$SpatialScale=="A"]<=0)/length(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="A" & DiffData$SpatialScale=="A"]), col="black", lwd=1.5) # observed value

## adult sub-patch
plot(density(Rand_dom_diff.prop.negatives.ASP), xlab='', main="N = 154", xlim=c(0,0.5), bty="n", cex.lab=1.8, cex.axis=1.6, type='n', ylab='')
polygon(density(Rand_dom_diff.prop.negatives.ASP), col=rainbow_hcl(10, alpha=0.4)[5], lwd=1.5)
abline(v=sort(Rand_dom_diff.prop.negatives.ASP)[length(Rand_dom_diff.prop.negatives.ASP)*0.975], col="black", lty=2, lwd=1.5) # upper CI
abline(v=sort(Rand_dom_diff.prop.negatives.ASP)[length(Rand_dom_diff.prop.negatives.ASP)*0.025], col="black", lty=2, lwd=1.5) # lower CI
abline(v=sum(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="A" & DiffData$SpatialScale=="B"]<=0)/length(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="A" & DiffData$SpatialScale=="B"]), col="black", lwd=1.5) # observed value

## juvenile patch
plot(density(Rand_dom_diff.prop.negatives.JP), xlab='Proportion of events with subordinate joiner', main="N = 349", xlim=c(0.5,0.8), bty="n", cex.lab=1.8, cex.axis=1.6, type='n')
polygon(density(Rand_dom_diff.prop.negatives.JP), col=rainbow_hcl(10, alpha=0.8)[9], lwd=1.5)
abline(v=sort(Rand_dom_diff.prop.negatives.JP)[length(Rand_dom_diff.prop.negatives.JP)*0.975], col="black", lty=2, lwd=1.5) # upper CI
abline(v=sort(Rand_dom_diff.prop.negatives.JP)[length(Rand_dom_diff.prop.negatives.JP)*0.025], col="black", lty=2, lwd=1.5) # lower CI
abline(v=sum(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="J" & DiffData$SpatialScale=="A"]<=0)/length(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="J" & DiffData$SpatialScale=="A"]), col="black", lwd=1.5) # observed value

## juvenile sub-patch
plot(density(Rand_dom_diff.prop.negatives.JSP), xlab='Proportion of events with subordinate joiner', main="N = 147", xlim=c(0.4,0.9), bty="n", cex.lab=1.8, cex.axis=1.6, type='n', ylab="")
polygon(density(Rand_dom_diff.prop.negatives.JSP), col=rainbow_hcl(10, alpha=0.4)[9], lwd=1.5)
abline(v=sort(Rand_dom_diff.prop.negatives.JSP)[length(Rand_dom_diff.prop.negatives.JSP)*0.975], col="black", lty=2, lwd=1.5) # upper CI
abline(v=sort(Rand_dom_diff.prop.negatives.JSP)[length(Rand_dom_diff.prop.negatives.JSP)*0.025], col="black", lty=2, lwd=1.5) # lower CI
abline(v=sum(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="J" & DiffData$SpatialScale=="B"]<=0)/length(DiffData$DiffMeanRank[DiffData$AgeClassBinary=="J" & DiffData$SpatialScale=="B"]), col="black", lwd=1.5) # observed value

###########################################################################################
## Analysis 3 -- Competitive exclusion at different spatial scales
## uses CompEx
###########################################################################################

## Variable explanations for those not clear:
## Supp -- binary indicator of whether the joined individual was competitively excluded (1) or not (0)
## SpatialScale -- patch (A) or sub-patch (B)

## make categorical variables non-numeric
CompEx$AgeClassBinary <- factor(c("J","A")[CompEx$AgeClassBinary])

levels(CompEx$Supp) <- c(0,1)

## build maximal model 
m.CompEx <- glmer(Supp ~ SpatialScale * AgeClassBinary + Troop + (1|FocalID) + (1|FollowID), data=CompEx, family = binomial)
summary(m.CompEx)

m.CompEx2 <- glmer(Supp ~ SpatialScale + AgeClassBinary + Troop + (1|FocalID) + (1|FollowID), data=CompEx, family = binomial)
summary(m.CompEx2) #final

## LR tests
drop1(m.CompEx, test="Chisq")
drop1(m.CompEx2, test="Chisq")