rm(list=ls())
library(lme4)
library(car)

master_root <- '~/Apes/Julia/newdata/'

# Load data
group <- "north"
dpn <- read.csv(paste("~/Apes/Julia/newdata/masterfile_meat_eating_",group,"_August_2015_presence_prepped.csv",sep=""),header=T)
group <- "south"
dps <- read.csv(paste("~/Apes/Julia/newdata/masterfile_meat_eating_",group,"_September_2015_presence_prepped.csv",sep=""),header=T)

dp <- rbind(dpn,dps)
length(unique(dp$EventID))

ab <- sapply(dpn$EventID,FUN=function(x){y=strsplit(as.character(x),
                                                    split='-')[[1]][1]})
ab <- as.numeric(ab)
plot(ab)
setdiff(1:max(ab),ab)


# MODELS PART I PROBABILITY OF BEING AT A MEAT EATING EVENT
# Presence data

pres.in.party.temp <- NA
pres.in.party.temp[dp$pres.in.party==0] <- "no"
pres.in.party.temp[dp$pres.in.party==1] <- "yes"
dp$pres.in.party <- as.factor(pres.in.party.temp)

estrus.temp <- NA
estrus.temp[dp$estrus==0] <- "no"
estrus.temp[dp$estrus==1] <- "yes"
dp$estrus <- as.factor(estrus.temp)
dp$estrus.yes <- 0
dp$estrus.yes[which(dp$estrus=="yes")] <- 1

is.lac.temp <- NA
is.lac.temp[dp$is.lac==0] <- "no"
is.lac.temp[dp$is.lac==1] <- "yes"
dp$is.lac <- as.factor(is.lac.temp)
dp$lac.yes <- 0
dp$lac.yes[which(dp$is.lac=="yes")] <- 1

dp$age <- round(dp$age/365)
dp$age.z <- as.numeric(scale(dp$age))
dp$nfemales.z <- as.numeric(scale(dp$nfemales))
dp$nmales.z <- as.numeric(scale(dp$nmales))
dp$communitysize.z <- as.numeric(scale(dp$communitysize))
dp$FAI.z <- as.numeric(scale(dp$FAI))

dp$rank <- relevel(dp$rank,ref="low")
dp$rank.middle <- 0
dp$rank.high <- 0
dp$rank.middle[which(as.character(dp$rank)=="middle")] <- 1
dp$rank.high[which(as.character(dp$rank)=="high")] <- 1
dp$rank.middle.FAI.z <- 0
dp$rank.middle.FAI.z[which(dp$rank.middle==1)] <- dp$FAI.z[which(dp$rank.middle==1)]
dp$rank.high.FAI.z <- 0
dp$rank.high.FAI.z[which(dp$rank.high==1)] <- dp$FAI.z[which(dp$rank.high==1)]
dp$rank.middle.nfemales.z <- 0
dp$rank.middle.nfemales.z[which(dp$rank.middle==1)] <- dp$nfemales.z[which(dp$rank.middle==1)]
dp$rank.high.nfemales.z <- 0
dp$rank.high.nfemales.z[which(dp$rank.high==1)] <- dp$nfemales.z[which(dp$rank.high==1)]
dp$rank.middle.group <- 0
dp$rank.middle.group[which(dp$group=="south")] <- 1
dp$rank.high.group <- 0
dp$rank.high.group[which(dp$group=="south")] <- 1

dp$RDate <- as.Date(paste0(dp$year,'/',dp$month,'/',dp$day),"%Y/%m/%d")


# nfemales and community size covary
xx <- glm(pres.in.party~communitysize.z+estrus+rank+age.z+nfemales.z+nmales.z+is.lac+group+FAI.z,family=binomial, data=dp)
vif(xx)
plot(dp$communitysize,dp$nfemales)

yy <- glm(pres.in.party~communitysize.z+estrus+rank.middle+rank.high+age.z+nfemales.z+nmales.z+is.lac+group+FAI.z,family=binomial, data=dp)

xx <- glm(pres.in.party~communitysize.z+estrus+rank+age.z+nfemales.z+nmales.z+is.lac+group+FAI.z+rank*FAI.z+rank*group,family=binomial, data=dp)
summary(xx)
#with(dp,boxplot(age~EventID))

# Summarize the data
length(unique(dp$EventID)) # Number of unique events
tot.rank.by.event <- tapply(dp$rank,dp$EventID,table)
sum(sapply(tot.rank.by.event,FUN=function(x){sum(names(tot.rank.by.event[[1]])=='low')>0}))
sum(sapply(tot.rank.by.event,FUN=function(x){sum(names(tot.rank.by.event[[1]])=='middle')>0}))
sum(sapply(tot.rank.by.event,FUN=function(x){sum(names(tot.rank.by.event[[1]])=='high')>0}))

# Part 1- Presence models

# Part 1A- Presence models with nfmales instead of community size
# Null model with only random effects to see which ones have estimated sd's away from zero
# For a given event, there will be no variation in FAI, so don't include this or its interactions in random slopes
#   for a given event
presence.re.structure.form <-as.formula(pres.in.party~
                                          1+(1|FID)+
                                          (0+age.z|FID)+
                                          (0+nfemales.z|FID)+
                                          (0+nmales.z|FID)+
                                          (0+FAI.z|FID)+
                                          (1|EventID)+
                                          (0+rank.middle|EventID)+
                                          (0+rank.high|EventID)+
                                          (0+age.z|EventID)+
                                          (0+lac.yes|EventID)+
                                          (0+estrus.yes|EventID)) 

pm.null <- glmer(presence.re.structure.form,
                 control=glmerControl(optimizer = "bobyqa"),
                 family=binomial,data=dp)
summary(pm.null)
# Singular fit, removing random slopes categorical pv's seems to fix this

presence.re.structure.form <-as.formula(pres.in.party~
                                          1+(1|FID)+
                                          (0+age.z|FID)+
                                          (0+nfemales.z|FID)+
                                          (0+nmales.z|FID)+
                                          (0+FAI.z|FID)+
                                          (1|EventID)+
                                          #(0+rank.middle|EventID)+
                                          #(0+rank.high|EventID)+
                                          (0+age.z|EventID))#+
                                          #(0+lac.yes|EventID)+
                                          #(0+estrus.yes|EventID)) 

pm.null <- glmer(presence.re.structure.form,
                 control=glmerControl(optimizer = "bobyqa"),
                 family=binomial,data=dp)
summary(pm.null)

# Estimates all away from zero.
pm.full.formula <- as.formula(pres.in.party~
                                rank+
                                is.lac+
                                estrus.yes+
                                age.z+
                                nmales.z+
                                nfemales.z+
                                FAI.z+
                                group+
                                rank*nfemales.z+
                                rank*FAI.z+
                                rank*group+
                                rank*estrus.yes+
                                (1|FID)+
                                (0+age.z|FID)+
                                (0+nfemales.z|FID)+
                                (0+nmales.z|FID)+
                                (0+FAI.z|FID)+
                                (1|EventID)+
                                (0+age.z|EventID))

pm.full <- glmer(pm.full.formula,family=binomial,
                 control=glmerControl(optimizer = "bobyqa"),data=dp)
summary(pm.full)
(pm.full.vs.null.sig.test<-anova(pm.full,pm.null,test="Chisq")) # Full model better than no fixed effects at all
pm.full.vs.null.sig.test$Chisq

# Test if interactions are needed one by one. LRTs suggest they are not
pm.no.rank.nfemales.inter.formula <- update(pm.full.formula,.~.-rank:nfemales.z)
pm.no.rank.nfemales.inter <- glmer(pm.no.rank.nfemales.inter.formula,
                                   control=glmerControl(optimizer = "bobyqa"),
                                   family=binomial,data=dp)
(pm.full.vs.no.rank.nfemales.inter.sig.test <- anova(pm.full,pm.no.rank.nfemales.inter,test="Chisq"))

pm.no.rank.FAI.inter.formula <- update(pm.full.formula,.~.-rank:FAI.z)
pm.no.rank.FAI.inter <- glmer(pm.no.rank.FAI.inter.formula,
                              control=glmerControl(optimizer = "bobyqa"),
                              family=binomial,data=dp)
(pm.full.vs.no.rank.FAI.inter.sig.test <- anova(pm.full,pm.no.rank.FAI.inter,test="Chisq"))

pm.no.rank.group.inter.formula <- update(pm.full.formula,.~.-rank:group)
pm.no.rank.group.inter <- glmer(pm.no.rank.group.inter.formula,
                                control=glmerControl(optimizer = "bobyqa"),
                                family=binomial,data=dp)
(pm.full.vs.no.rank.group.inter.sig.test<-anova(pm.full,pm.no.rank.group.inter,test="Chisq"))

pm.no.rank.estrus.inter.formula <- update(pm.full.formula,.~.-rank:estrus.yes)
pm.no.rank.estrus.inter <- glmer(pm.no.rank.estrus.inter.formula,
                                 control=glmerControl(optimizer = "bobyqa"),
                                 family=binomial,data=dp)
(pm.full.vs.no.rank.estrus.inter.sig.test<-anova(pm.full,pm.no.rank.estrus.inter,test="Chisq"))

# Fit model without interactions so terms involved in interactions are directly interpretable
pm.no.inter.formula <- update(pm.full.formula,.~.-rank:nfemales.z-rank:FAI.z-rank:group-rank:estrus.yes)
pm.no.inter <- glmer(pm.no.inter.formula,
                     control=glmerControl(optimizer = "bobyqa"),
                     family=binomial,data=dp)
# Model without interactions supported over null model
(pm.no.inter.vs.null.sig.test<-anova(pm.no.inter,pm.null,test="Chisq"))
# Model without interactions not worse than one with interactions
anova(pm.full,pm.no.inter,test="Chisq")

# LRTs for each individual term
pm.terms.LRT.test <- drop1(pm.no.inter,test="Chisq")

# Test for residual autocorrelation to address whether time since last recorded 
# event of an individual being present has an influence.
r <- resid(pm.no.inter)
rn <- r[which(dp$group=="north")]
rs <- r[which(dp$group=="south")]

dp$resid <- r

uFID <- unique(dp$FID)
pres.cor.store <- data.frame(FID=uFID,cor=NA,pval=NA,acf=NA)
for(ui in 1:length(uFID)){
  dp.sub <- subset(dp,FID==uFID[ui])
  diff.time <- as.numeric(difftime(dp.sub$RDate[-1],dp.sub$RDate[-nrow(dp.sub)],
                                   units='days'))
  if(length(diff.time)>3){
    cor.temp <- lm(dp.sub$resid[-1]~diff.time)
    pres.cor.store$cor[ui] <- coef(cor.temp)[2]
    pres.cor.store$pval[ui] <- summary(cor.temp)$coefficients[2,4]
    cor.temp.acf <- acf(dp.sub$resid[-1],ci.type='ma')
    pres.cor.store.alt$acf[ui] <- max(cor.temp.acf$acf[-1])
  }  
}
pres.cor.store

# Make table summarizing results
export.table <- summary(pm.no.inter)$coefficients[,1:3]
export.table <- rbind(export.table,summary(pm.full)$coefficients[-(1:nrow(summary(pm.no.inter)$coefficients)),1:3])
export.table <- as.data.frame(export.table)
export.table$LRT <- NA
export.table$Df <- NA
export.table$p.val <- NA

rownames(pm.terms.LRT.test)[which(rownames(pm.terms.LRT.test)=="rank")] <- "rankhigh"
rownames(pm.terms.LRT.test)[which(rownames(pm.terms.LRT.test)=="is.lac")] <- "is.lacyes"
rownames(pm.terms.LRT.test)[which(rownames(pm.terms.LRT.test)=="group")] <- "groupsouth"
index.temp <- match(rownames(pm.terms.LRT.test),rownames(export.table))
index.temp <- index.temp[is.finite(index.temp)]

export.table$LRT[index.temp] <- pm.terms.LRT.test$LRT[-1]
export.table$Df[index.temp] <- pm.terms.LRT.test$Df[-1]
export.table$p.val[index.temp] <- pm.terms.LRT.test$"Pr(Chi)"[-1]

f.extract<-function(x){
	LRT=-2*(x$logLik[1]-x$logLik[2])
	Df=x$Df[2]-x$Df[1]
	P.value=pchisq(LRT,df=Df,lower.tail=F)
	return(c(LRT,Df,P.value))
}

f.extract(pm.full.vs.null.sig.test)

#LRT.values.temp <- f.extract(pm.no.inter.vs.null.sig.test)
#export.table["(Intercept)",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(pm.full.vs.no.rank.group.inter.sig.test)
export.table["rankhigh:groupsouth",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(pm.full.vs.no.rank.nfemales.inter.sig.test)
export.table["rankhigh:nfemales.z",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(pm.full.vs.no.rank.FAI.inter.sig.test)
export.table["rankhigh:FAI.z",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(pm.full.vs.no.rank.estrus.inter.sig.test)
export.table["rankhigh:estrus.yes",c("LRT","Df","p.val")] <- LRT.values.temp

export.table <- round(export.table,3)

export.table <- data.frame(Term=rownames(export.table),export.table)
rownames(export.table) <- NULL
pm.export.table <- export.table
write.csv(format(pm.export.table,digits=3),
          file=file.path(master_root,"Table-meat-presence-model-results-with-rank-estrus-interaction-Jan-2019.csv"),
          row.names=F)

# Part 1B- Presence model using community size instead of number of females
pm.null.alt <- glmer(pres.in.party~
                       (1|FID)+
                       (0+age.z|FID)+
                       (0+communitysize.z|FID)+
                       (0+nmales.z|FID)+
                       (0+FAI.z|FID)+
                       (1|EventID)+
                       (0+age.z|EventID),
                     family=binomial,control=glmerControl(optimizer = "bobyqa"),data=dp)
summary(pm.null.alt)
# Random effects SD estimates away from zero

pm.full.formula.alt <- as.formula(pres.in.party~
                                    rank+
                                    is.lac+
                                    estrus.yes+
                                    age.z+
                                    nmales.z+
                                    communitysize.z+
                                    FAI.z+
                                    group+
                                    rank*communitysize.z+
                                    rank*FAI.z+
                                    rank*group+
                                    rank*estrus.yes+
                                    (1|FID)+
                                    (0+age.z|FID)+
                                    (0+communitysize.z|FID)+
                                    (0+nmales.z|FID)+
                                    (0+FAI.z|FID)+
                                    (1|EventID)+
                                    (0+age.z|EventID))

pm.full.alt <- glmer(pm.full.formula.alt,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=dp)
summary(pm.full.alt)
(pm.full.vs.null.sig.test.alt<-anova(pm.full.alt,pm.null.alt,test="Chisq")) 
# Full model better than no fixed effects at all

# Test if interactions are needed one by one. LRTs suggest they are not
pm.no.rank.communitysize.inter.formula.alt <- update(pm.full.formula.alt,.~.-rank:communitysize.z)
pm.no.rank.communitysize.inter.alt <- glmer(pm.no.rank.communitysize.inter.formula.alt,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=dp)
(pm.full.vs.no.rank.communitysize.inter.sig.test.alt <- anova(pm.full.alt,pm.no.rank.communitysize.inter.alt,
                                                              test="Chisq"))

pm.no.rank.FAI.inter.formula.alt <- update(pm.full.formula.alt,.~.-rank:FAI.z)
pm.no.rank.FAI.inter.alt <- glmer(pm.no.rank.FAI.inter.formula.alt,
                                  family=binomial,
                                  control=glmerControl(optimizer = "bobyqa"),data=dp)
(pm.full.vs.no.rank.FAI.inter.sig.test.alt <- anova(pm.full.alt,
                                                    pm.no.rank.FAI.inter.alt,
                                                    test="Chisq"))

pm.no.rank.group.inter.formula.alt <- update(pm.full.formula.alt,.~.-rank:group)
pm.no.rank.group.inter.alt <- glmer(pm.no.rank.group.inter.formula.alt,
                                    family=binomial,
                                    control=glmerControl(optimizer = "bobyqa"),data=dp)
(pm.full.vs.no.rank.group.inter.sig.test.alt<-anova(pm.full.alt,
                                                    pm.no.rank.group.inter.alt,
                                                    test="Chisq"))

pm.no.rank.estrus.inter.formula.alt <- update(pm.full.formula.alt,.~.-rank:estrus.yes)
pm.no.rank.estrus.inter.alt <- glmer(pm.no.rank.estrus.inter.formula.alt,
                                     family=binomial,
                                     control=glmerControl(optimizer = "bobyqa"),data=dp)
(pm.full.vs.no.rank.estrus.inter.sig.test.alt<-anova(pm.full.alt,
                                                     pm.no.rank.estrus.inter.alt,
                                                     test="Chisq"))

# Fit model without interactions so terms involved in interactions are directly interpretable
pm.no.inter.formula.alt <- update(pm.full.formula.alt,.~.-rank:communitysize.z-rank:FAI.z-rank:group-rank:estrus.yes)
pm.no.inter.alt <- glmer(pm.no.inter.formula.alt,
                         family=binomial,
                         control=glmerControl(optimizer = "bobyqa"),
                         data=dp)
(pm.no.inter.vs.null.sig.test.alt<-anova(pm.no.inter.alt,pm.null.alt,test="Chisq"))  # Model without interactions supported over null model

anova(pm.full.alt,pm.no.inter.alt,test="Chisq")

# LRTs for each individual term
pm.terms.LRT.test.alt <- drop1(pm.no.inter.alt,test="Chisq")

# Test for residual autocorrelation to address whether time since last recorded event of an individual being present has an influence.
r.alt <- resid(pm.no.inter.alt)
rn.alt <- r.alt[which(dp$group=="north")]
rs.alt <- r.alt[which(dp$group=="south")]

dp$resid.alt <- r.alt

pres.cor.store.alt <- data.frame(FID=uFID,cor=NA,pval=NA,acf=NA)
for(ui in 1:length(uFID)){
  dp.sub <- subset(dp,FID==uFID[ui])
  diff.time <- as.numeric(difftime(dp.sub$RDate[-1],dp.sub$RDate[-nrow(dp.sub)],
                                   units='days'))
  if(ui %in% c(1,26,51)){
    quartz()
    par(mfrow=c(5,5),mar=c(3,3,2,1))
  }
  if(length(diff.time)>3){
    cor.temp.alt <- lm(dp.sub$resid.alt[-1]~diff.time)
    pres.cor.store.alt$cor[ui] <- coef(cor.temp)[2]
    pres.cor.store.alt$pval[ui] <- summary(cor.temp)$coefficients[2,4]
    cor.temp.acf <- acf(dp.sub$resid.alt[-1],ci.type='ma')
    pres.cor.store.alt$acf[ui] <- max(cor.temp.acf$acf[-1])
  } 
}




# Make table summarizing results
export.table.alt <- summary(pm.no.inter.alt)$coefficients[,1:3]
export.table.alt <- rbind(export.table.alt,summary(pm.full.alt)$coefficients[-(1:nrow(summary(pm.no.inter.alt)$coefficients)),1:3])
export.table.alt <- as.data.frame(export.table.alt)
export.table.alt$LRT <- NA
export.table.alt$Df <- NA
export.table.alt$p.val <- NA

rownames(pm.terms.LRT.test.alt)[which(rownames(pm.terms.LRT.test.alt)=="rank")] <- "rankhigh"
rownames(pm.terms.LRT.test.alt)[which(rownames(pm.terms.LRT.test.alt)=="is.lac")] <- "is.lacyes"
rownames(pm.terms.LRT.test.alt)[which(rownames(pm.terms.LRT.test.alt)=="group")] <- "groupsouth"
index.temp.alt <- match(rownames(pm.terms.LRT.test.alt),rownames(export.table.alt))

export.table.alt$LRT[index.temp.alt[is.finite(index.temp.alt)]] <- pm.terms.LRT.test.alt$LRT[which(is.finite(index.temp.alt))]
export.table.alt$Df[index.temp.alt[is.finite(index.temp.alt)]] <- pm.terms.LRT.test.alt$Df[which(is.finite(index.temp.alt))]
export.table.alt$p.val[index.temp.alt[is.finite(index.temp.alt)]] <- pm.terms.LRT.test.alt[,"Pr(Chi)"][which(is.finite(index.temp.alt))]

f.extract<-function(x){
  LRT=-2*(x$logLik[1]-x$logLik[2])
  Df=x$Df[2]-x$Df[1]
  P.value=pchisq(LRT,df=Df,lower.tail=F)
  return(c(LRT,Df,P.value))
}

f.extract(pm.full.vs.null.sig.test.alt)

#LRT.values.temp.alt <- f.extract(pm.no.inter.vs.null.sig.test.alt)
#export.table.alt["(Intercept)",c("LRT","Df","p.val")] <- LRT.values.temp.alt

LRT.values.temp.alt <- f.extract(pm.full.vs.no.rank.group.inter.sig.test.alt)
export.table.alt["rankhigh:groupsouth",c("LRT","Df","p.val")] <- LRT.values.temp.alt

LRT.values.temp.alt <- f.extract(pm.full.vs.no.rank.communitysize.inter.sig.test.alt)
export.table.alt["rankhigh:communitysize.z",c("LRT","Df","p.val")] <- LRT.values.temp.alt

LRT.values.temp.alt <- f.extract(pm.full.vs.no.rank.FAI.inter.sig.test.alt)
export.table.alt["rankhigh:FAI.z",c("LRT","Df","p.val")] <- LRT.values.temp.alt

LRT.values.temp.alt <- f.extract(pm.full.vs.no.rank.estrus.inter.sig.test.alt)
export.table.alt["rankhigh:estrus.yes",c("LRT","Df","p.val")] <- LRT.values.temp.alt

export.table.alt <- round(export.table.alt,3)

export.table.alt <- data.frame(Term=rownames(export.table.alt),export.table.alt)
rownames(export.table.alt) <- NULL
pm.export.table.alt <- export.table.alt
write.csv(format(pm.export.table.alt,digits=3),
          file=file.path(master_root,'Table-meat-presence-model-results-with-rank-estrus-interaction-Jan-2019-alt.csv'),
          row.names=F)

# PART 2 PROBABILITY OF GETTING MEAT GIVEN PRESENCE AT MEAT EATING EVENT
# Sharing data
group <- "north"
den <- read.csv(paste("~/Apes/Julia/newdata/masterfile_meat_eating_",group,"_August_2015_meat_prepped.csv",sep=""),header=T)
group <- "south"
des <- read.csv(paste("~/Apes/Julia/newdata/masterfile_meat_eating_",group,"_September_2015_meat_prepped.csv",sep=""),header=T)
de <- rbind(den,des)

# Number of observed meat eating events: 827
length(unique(de$EventID))

# 54 meat sharing events with no adult females present. Remove
# or 42???? 42 with no females at all, 54 with no adult females
#sum(sapply(tapply(de$sex,de$EventID,table),FUN=function(x){x['F']})==0)

ab <- tapply(de$rank,de$EventID,FUN=function(x){sum(as.numeric(!is.na(x)))})
(ab=ab[which(ab==0)])  # Meat sharing events with no adult females present. Remove?
length(ab)  # 54 events with no adult females are revmoved

head(which(de$sex=='F' & is.na(de$rank)))


# Use 773 remaining events for analyses
nrow(de)-length(which(!is.finite(match(de$EventID,names(ab)))))
de <- de[which(!is.finite(match(de$EventID,names(ab)))),]
de$EventID <- as.character(de$EventID)
ab <- tapply(de$rank,de$EventID,FUN=function(x){sum(as.numeric(!is.na(x)))})
(ab <- ab[which(ab==0)])  

# Only adult females used in analysis
de <- subset(de,sex=="F" & age>(13*365)) 
de$sex <- as.character(de$sex)
unique(de$sex)

meat.temp <- NA
meat.temp[de$meat==0] <- "no"
meat.temp[de$meat==1] <- "yes"
de$meat <- as.factor(meat.temp)

estrus.temp <- NA
estrus.temp[de$estrus==0] <- "no"
estrus.temp[de$estrus==1] <- "yes"
de$estrus <- as.factor(estrus.temp)
de$estrus.yes <- 0
de$estrus.yes[which(de$estrus=="yes")] <- 1

is.lac.temp <- NA
is.lac.temp[de$is.lac==0] <- "no"
is.lac.temp[de$is.lac==1] <- "yes"
de$is.lac <- as.factor(is.lac.temp)
de$lac.yes <- 0
de$lac.yes[which(de$is.lac=="yes")] <- 1

de$age <- round(de$age/365)
de$age.z <- as.numeric(scale(de$age))
de$nfemales.z <- as.numeric(scale(de$nfemales))
de$nmales.z <- as.numeric(scale(de$nmales))
de$communitysize.z <- as.numeric(scale(de$communitysize))
de$FAI.z <- as.numeric(scale(de$FAI))

de$rank <- relevel(de$rank,ref="low")
de$rank.middle <- 0
de$rank.high <- 0
de$rank.middle[which(as.character(de$rank)=="middle")] <- 1
de$rank.high[which(as.character(de$rank)=="high")] <- 1

de$rank.middle.FAI.z <- 0
de$rank.middle.FAI.z[which(de$rank.middle==1)] <- de$FAI.z[which(de$rank.middle==1)]
de$rank.high.FAI.z <- 0
de$rank.high.FAI.z[which(de$rank.high==1)] <- de$FAI.z[which(de$rank.high==1)]

de$rank.middle.nfemales.z <- 0
de$rank.middle.nfemales.z[which(de$rank.middle==1)] <- de$nfemales.z[which(de$rank.middle==1)]
de$rank.high.nfemales.z <- 0
de$rank.high.nfemales.z[which(de$rank.high==1)] <- de$nfemales.z[which(de$rank.high==1)]

de$rank.middle.communitysize.z <- 0
de$rank.middle.communitysize.z[which(de$rank.middle==1)] <- de$communitysize.z[which(de$rank.middle==1)]
de$rank.high.communitysize.z <- 0
de$rank.high.communitysize.z[which(de$rank.high==1)] <- de$communitysize.z[which(de$rank.high==1)]

de$RDate <- as.Date(paste0(de$year,'/',de$month,'/',de$day),"%Y/%m/%d")

# Stats
table(de$meat)
with(de,table(meat,rank))
ab <- with(de,table(meat,rank))
ab <- ab[,c("low","middle","high")]
barplot(ab)
barplot(ab/rbind(apply(ab,2,sum),apply(ab,2,sum)))

range(dp$RDate)

length(unique(dp$FID))
length(unique(de$who))


ab <- tapply(dp$pres.in.party,dp$FID,table)
which(sapply(ab,FUN=function(x){x['yes']})==0) #Individuals never eating meat

# Number of events with at least one individual of a given rank present
ab <- tapply(de$rank,de$EventID,table)
sum(sapply(ab,FUN=function(x){x['low']})>0)
sum(sapply(ab,FUN=function(x){x['middle']})>0)
sum(sapply(ab,FUN=function(x){x['high']})>0)

# Proportion of times each rank got meat given the rank class was present
ab <- subset(de,rank=='low')
length(unique(ab$EventID)) # Number of events at least one low ranking individual present
sum(tapply(ab$meat,ab$EventID,FUN=function(x){sum(x=='yes')})>0)
sum(tapply(ab$meat,ab$EventID,FUN=function(x){sum(x=='yes')})>0)/length(unique(ab$EventID))

ab <- subset(de,rank=='middle')
length(unique(ab$EventID)) # Number of events at least one low ranking individual present
sum(tapply(ab$meat,ab$EventID,FUN=function(x){sum(x=='yes')})>0)
sum(tapply(ab$meat,ab$EventID,FUN=function(x){sum(x=='yes')})>0)/length(unique(ab$EventID))

ab <- subset(de,rank=='high')
length(unique(ab$EventID)) # Number of events at least one low ranking individual present
sum(tapply(ab$meat,ab$EventID,FUN=function(x){sum(x=='yes')})>0)
sum(tapply(ab$meat,ab$EventID,FUN=function(x){sum(x=='yes')})>0)/length(unique(ab$EventID))



# Meat events on different days
length(unique(de$EventID))
de$RDate <- as.Date(paste(de$year,de$month,de$day,sep="-"),format="%Y-%m-%d")
length(unique(de$RDate))
1-(length(unique(de$EventID))-length(unique(de$RDate)))/length(unique(de$EventID))  # Percent of days
unique.RDate <- tapply(de$RDate,de$group,FUN=function(x){unique(x)})

inter.event.times.north <- as.numeric(difftime(unique.RDate$north[-1],unique.RDate$north[-length(unique.RDate$north)]))
inter.event.times.south <- as.numeric(difftime(unique.RDate$south[-1],unique.RDate$south[-length(unique.RDate$south)]))
mean(inter.event.times.north)
mean(inter.event.times.south)
range(inter.event.times.north)
range(inter.event.times.south)

xx <- glm(meat~rank+is.lac+estrus+age.z+nfemales.z+nmales.z+group+communitysize.z+FAI.z+rank*communitysize.z+rank*FAI.z,family=binomial,data=de)
yy <- glm(meat~rank.middle+rank.high+is.lac+estrus+age.z+nfemales.z+nmales.z+group+communitysize.z+FAI.z+rank.middle.communitysize.z+rank.high.communitysize.z+rank.middle.FAI.z+rank.high.FAI.z,family=binomial, data=de)

with(de,boxplot(age~who))
with(de,boxplot(nfemales~who))
with(de,boxplot(nmales~who))
with(de,boxplot(FAI~who))

range(tapply(de$rank,de$EventID,FUN=function(x){length(unique(x))}))
tapply(de$rank,de$EventID,FUN=function(x){length(unique(x))})
ab=subset(de,EventID=="95-south")

# Sharing Models 
xx <- glm(meat~rank+is.lac+estrus+age.z+nfemales.z+nmales.z+group+
            communitysize.z+FAI.z,family=binomial,data=de)
vif(xx)
plot(de$communitysize,de$nfemales)

em.null.all.slopes <- glmer(meat~
                              (1|who)+
                              (0+age.z|who)+
                              (0+nfemales.z|who)+
                              (0+nmales.z|who)+
                              (0+FAI.z|who)+
                              (0+communitysize.z|who)+
                              (1|EventID)+
                              (0+rank.middle|EventID)+
                              (0+rank.high|EventID)+
                              (0+age.z|EventID)+
                              (0+lac.yes|EventID)+
                              (0+estrus.yes|EventID),
                            family=binomial,
                            control=glmerControl(optimizer = "bobyqa"),data=de)
summary(em.null.all.slopes)
# Several r.e. sd estimates near zero and singular fit
# Remove: estrus.yes, lac.yes, rank (high and middle together), age within event
# Any random slopes lead to singular fit when main effects are included so just 
# have random intercepts

em.null <- glmer(meat~
                   (1|who)+
                   #(0+age.z|who)+
                   #(0+nfemales.z|who)+
                   #(0+nmales.z|who)+
                   #(0+FAI.z|who)+
                   #(0+communitysize.z|who)+
                   (1|EventID),
                 control=glmerControl(optimizer = "bobyqa"),
                 family=binomial,data=de)
summary(em.null)
# RE sd estimatates all away from zero

em.full.formula <- as.formula(meat~
	rank+
	is.lac+
	estrus.yes+
	age.z+
	nmales.z+
	nfemales.z+
	FAI.z+
	communitysize.z+
	group+
	rank*nfemales.z+
	rank*FAI.z+
	rank*group+
	rank*estrus.yes+
	  (1|who)+
	  #(0+age.z|who)+
	  #(0+nfemales.z|who)+
	  #(0+nmales.z|who)+
	  #(0+FAI.z|who)+
	  #(0+communitysize.z|who)+
	  (1|EventID))

em.full <- glmer(em.full.formula,
                 control=glmerControl(optimizer = "bobyqa"),
                 family=binomial,data=de)

summary(em.full)
(em.full.vs.null.sig.test<-anova(em.full,em.null,test="Chisq")) # Full model better than no fixed effects at all
em.full.vs.null.sig.test$Chisq

# Test if interactions are needed one by one. LRTs suggest they are not
em.no.rank.nfemales.inter.formula <- update(em.full.formula,.~.-rank:nfemales.z)
em.no.rank.nfemales.inter <- glmer(em.no.rank.nfemales.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de)
(em.full.vs.no.rank.nfemales.inter.sig.test <- anova(em.full,em.no.rank.nfemales.inter,test="Chisq"))

em.no.rank.FAI.inter.formula <- update(em.full.formula,.~.-rank:FAI.z)
em.no.rank.FAI.inter <- glmer(em.no.rank.FAI.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de)
(em.full.vs.no.rank.FAI.inter.sig.test <- anova(em.full,em.no.rank.FAI.inter,test="Chisq"))

em.no.rank.group.inter.formula <- update(em.full.formula,.~.-rank:group)
em.no.rank.group.inter <- glmer(em.no.rank.group.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de)
(em.full.vs.no.rank.group.inter.sig.test<-anova(em.full,em.no.rank.group.inter,test="Chisq"))

em.no.rank.estrus.inter.formula <- update(em.full.formula,.~.-rank:estrus.yes)
em.no.rank.estrus.inter <- glmer(em.no.rank.estrus.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de)
(em.full.vs.no.rank.estrus.inter.sig.test<-anova(em.full,em.no.rank.estrus.inter,test="Chisq"))

# Fit model without interactions so terms involved in interactions are directly interpretable
em.no.inter.formula <- update(em.full.formula,.~.-rank:nfemales.z-rank:FAI.z-rank:group-rank:estrus.yes)
em.no.inter <- glmer(em.no.inter.formula,
                     control=glmerControl(optimizer = "bobyqa"),
                     family=binomial,data=de)
# No interactions better than null
(em.no.inter.vs.null.sig.test<-anova(em.no.inter,em.null,test="Chisq"))  # Model without interactions supported over null model
em.no.inter.vs.null.sig.test$Chisq
# No interactions not worse than one with interactions
anova(em.full,em.no.inter,test="Chisq")

# LRTs for each individual term
em.terms.LRT.test <- drop1(em.no.inter,test="Chisq")

# Residual analysis
de$resid <- resid(em.no.inter)

em.uwho <- unique(de$who)
em.cor.store <- data.frame(who=em.uwho,cor=NA,pval=NA,acf=NA)
for(ui in 1:length(em.uwho)){
  de.sub <- subset(de,who==em.uwho[ui])
  diff.time <- as.numeric(difftime(de.sub$RDate[-1],de.sub$RDate[-nrow(de.sub)],
                                   units='days'))
  if(ui %in% c(1,26,51)){
    quartz()
    par(mfrow=c(5,5),mar=c(3,3,2,1))
  }
  if(length(diff.time)>3){
    cor.temp <- lm(de.sub$resid[-1]~diff.time)
    em.cor.store$cor[ui] <- coef(cor.temp)[2]
    em.cor.store$pval[ui] <- summary(cor.temp)$coefficients[2,4]
    cor.temp.acf <- acf(de.sub$resid[-1],ci.type='ma')
    em.cor.store$acf[ui] <- max(cor.temp.acf$acf[-1])
  }  
}

# Make table summarizing results
export.table <- summary(em.no.inter)$coefficients[,1:3]
export.table <- rbind(export.table,summary(em.full)$coefficients[-(1:nrow(summary(em.no.inter)$coefficients)),1:3])
export.table <- as.data.frame(export.table)
export.table$LRT <- NA
export.table$Df <- NA
export.table$p.val <- NA

rownames(em.terms.LRT.test)[which(rownames(em.terms.LRT.test)=="rank")] <- "rankhigh"
rownames(em.terms.LRT.test)[which(rownames(em.terms.LRT.test)=="is.lac")] <- "is.lacyes"
rownames(em.terms.LRT.test)[which(rownames(em.terms.LRT.test)=="group")] <- "groupsouth"
index.temp <- match(rownames(em.terms.LRT.test),rownames(export.table))
index.temp <- index.temp[is.finite(index.temp)]
                         
export.table$LRT[index.temp] <- em.terms.LRT.test$LRT[-1]
export.table$Df[index.temp] <- em.terms.LRT.test$Df[-1]
export.table$p.val[index.temp] <- em.terms.LRT.test$"Pr(Chi)"[-1]

f.extract<-function(x){
	LRT=-2*(x$logLik[1]-x$logLik[2])
	Df=x$Df[2]-x$Df[1]
	P.value=pchisq(LRT,df=Df,lower.tail=F)
	return(c(LRT,Df,P.value))
}
f.extract(em.full.vs.null.sig.test)

#LRT.values.temp <- f.extract(em.no.inter.vs.null.sig.test)
#export.table["(Intercept)",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(em.full.vs.no.rank.group.inter.sig.test)
export.table["rankhigh:groupsouth",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(em.full.vs.no.rank.nfemales.inter.sig.test)
export.table["rankhigh:nfemales.z",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(em.full.vs.no.rank.FAI.inter.sig.test)
export.table["rankhigh:FAI.z",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(em.full.vs.no.rank.estrus.inter.sig.test)
export.table["rankhigh:estrus.yes",c("LRT","Df","p.val")] <- LRT.values.temp

export.table <- round(export.table,3)

export.table <- data.frame(Term=rownames(export.table),export.table)
rownames(export.table) <- NULL
em.export.table <- export.table
write.csv(format(em.export.table,digits=3),
          file=file.path(master_root,"Table-meat-access-model-results-Jan-2019.csv"),
          row.names=F)


# PART 3
# Proportion of events with at least one of each rank present
ab <- table(de$Event,de$rank)
apply(ab,2,median)
apply(ab,2,range)

ef <- apply(ab,1,FUN=function(x){sum(x==0)})  # Count the number of missing ranks
# ef <- apply(ab,1,FUN=function(x){y=x; y[1]==y[2] & y[1]==y[3]}) Number of events for which the number of individuals in each rank class are equal
sum(ef>=1)   # Number of events with at least one rank not being present: 309 have at least one rank not present
length(unique(de$EventID))-sum(ef>=1)  # 464 events with at least one individual from each of the three rank classes present

# subset of meat eating events where at least one of each rank is present
all.ranks.pres.events <- names(ef[which(ef==0)])
de.rel <- de[is.finite(match(de$EventID,all.ranks.pres.events)),]
length(unique(de.rel$EventID))

de.rel$age.z <- as.numeric(scale(de.rel$age))
de.rel$nfemales.z <- as.numeric(scale(de.rel$nfemales))
de.rel$nmales.z <- as.numeric(scale(de.rel$nmales))
de.rel$communitysize.z <- as.numeric(scale(de.rel$communitysize))
de.rel$FAI.z <- as.numeric(scale(de.rel$FAI))

de.rel$rank.middle.FAI.z <- 0
de.rel$rank.middle.FAI.z[which(de.rel$rank.middle==1)] <- de.rel$FAI.z[which(de.rel$rank.middle==1)]
de.rel$rank.high.FAI.z <- 0
de.rel$rank.high.FAI.z[which(de.rel$rank.high==1)] <- de.rel$FAI.z[which(de.rel$rank.high==1)]

de.rel$rank.middle.nfemales.z <- 0
de.rel$rank.middle.nfemales.z[which(de.rel$rank.middle==1)] <- de.rel$nfemales.z[which(de.rel$rank.middle==1)]
de.rel$rank.high.nfemales.z <- 0
de.rel$rank.high.nfemales.z[which(de.rel$rank.high==1)] <- de.rel$nfemales.z[which(de.rel$rank.high==1)]

de.rel$rank.middle.communitysize.z <- 0
de.rel$rank.middle.communitysize.z[which(de.rel$rank.middle==1)] <- de.rel$communitysize.z[which(de.rel$rank.middle==1)]
de.rel$rank.high.communitysize.z <- 0
de.rel$rank.high.communitysize.z[which(de.rel$rank.high==1)] <- de.rel$communitysize.z[which(de.rel$rank.high==1)]

# Fit models
em.rel.null.all.slopes <- glmer(meat~
	(1|who)+(0+age.z|who)+(0+nfemales.z|who)+(0+nmales.z|who)+(0+FAI.z|who)+(0+communitysize.z|who)+
	(1|EventID)+(0+rank.middle|EventID)+(0+rank.high|EventID)+(0+age.z|EventID)+(0+lac.yes|EventID)+(0+estrus.yes|EventID),
	family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de.rel)
summary(em.rel.null.all.slopes)
# Several r.e. sd estimates near zero. Remove: estrus.yes, lac.yes, rank (high and middle together), age within event

em.rel.null <- glmer(meat~
	(1|who)+(0+age.z|who)+(0+nfemales.z|who)+(0+nmales.z|who)+(0+FAI.z|who)+(0+communitysize.z|who)+(1|EventID),
	family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de.rel)
summary(em.rel.null)
# RE sd estimatates all away from zero

em.rel.full <- glmer(em.full.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de.rel)
summary(em.rel.full)
(em.rel.full.vs.null.sig.test<-anova(em.rel.full,em.rel.null,test="Chisq")) # Full model better than no fixed effects at all

# Test if interactions are needed one by one. LRTs suggest they are not
em.rel.no.rank.nfemales.inter <- glmer(em.no.rank.nfemales.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de.rel)
(em.rel.full.vs.no.rank.nfemales.inter.sig.test <- anova(em.rel.full,em.rel.no.rank.nfemales.inter,test="Chisq"))

em.rel.no.rank.FAI.inter <- glmer(em.no.rank.FAI.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de.rel)
(em.rel.full.vs.no.rank.FAI.inter.sig.test <- anova(em.rel.full,em.rel.no.rank.FAI.inter,test="Chisq"))

em.rel.no.rank.group.inter <- glmer(em.no.rank.group.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de.rel)
(em.rel.full.vs.no.rank.group.inter.sig.test<-anova(em.rel.full,em.rel.no.rank.group.inter,test="Chisq"))

# Fit model without interactions so terms involved in interactions are directly interpretable
em.rel.no.inter <- glmer(em.no.inter.formula,family=binomial,control=glmerControl(optimizer = "bobyqa"),data=de.rel)
summary(em.rel.no.inter)
(em.rel.no.inter.vs.null.sig.test<-anova(em.rel.no.inter,em.rel.null,test="Chisq"))  # Model without interactions supported over null model

# LRTs for each individual term
em.rel.terms.LRT.test <- drop1(em.rel.no.inter,test="Chisq")

# Make table summarizing results
export.table <- summary(em.rel.no.inter)$coefficients[,1:3]
export.table <- rbind(export.table,summary(em.rel.full)$coefficients[-(1:nrow(summary(em.rel.no.inter)$coefficients)),1:3])
export.table <- as.data.frame(export.table)
export.table$LRT <- NA
export.table$Df <- NA
export.table$p.val <- NA

rownames(em.rel.terms.LRT.test)[which(rownames(em.rel.terms.LRT.test)=="rank")] <- "rankhigh"
rownames(em.rel.terms.LRT.test)[which(rownames(em.rel.terms.LRT.test)=="is.lac")] <- "is.lacyes"
index.temp <- match(rownames(em.rel.terms.LRT.test),rownames(export.table))

export.table$LRT[index.temp[is.finite(index.temp)]] <- em.rel.terms.LRT.test$LRT[which(is.finite(index.temp))]
export.table$Df[index.temp[is.finite(index.temp)]] <- em.rel.terms.LRT.test$Df[which(is.finite(index.temp))]
export.table$p.val[index.temp[is.finite(index.temp)]] <- em.rel.terms.LRT.test[,"Pr(Chi)"][which(is.finite(index.temp))]

f.extract<-function(x){
	LRT=-2*(x$logLik[1]-x$logLik[2])
	Df=x$Df[2]-x$Df[1]
	P.value=pchisq(LRT,df=Df,lower.tail=F)
	return(c(LRT,Df,P.value))
}

#LRT.values.temp <- f.extract(em.rel.no.inter.vs.null.sig.test)
#export.table["(Intercept)",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(em.rel.full.vs.no.rank.group.inter.sig.test)
export.table["rankhigh:groupsouth",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(em.rel.full.vs.no.rank.nfemales.inter.sig.test)
export.table["rankhigh:nfemales.z",c("LRT","Df","p.val")] <- LRT.values.temp

LRT.values.temp <- f.extract(em.rel.full.vs.no.rank.FAI.inter.sig.test)
export.table["rankhigh:FAI.z",c("LRT","Df","p.val")] <- LRT.values.temp

export.table <- round(export.table,3)

export.table <- data.frame(Term=rownames(export.table),export.table)
rownames(export.table) <- NULL
em.rel.export.table <- export.table
write.csv(em.rel.export.table,
          file=file.path(master_root,
                         "Table-meat-access-all-ranks-present-model-results-Jan-2019.csv"),
          row.names=F)


# Figure 1- Percent of events where individuals are there by rank, and accessed by rank

# Proportion of events with each rank there
de.unique.event <- unique(de$EventID)
de.ab <- as.list(tapply(de$rank,de$EventID,unique))
de.low.per <- sum(unlist(lapply(de.ab,FUN=function(x){sum(x=="low")})))/length(de.unique.event)
de.mid.per <- sum(unlist(lapply(de.ab,FUN=function(x){sum(x=="middle")})))/length(de.unique.event)
de.high.per <- sum(unlist(lapply(de.ab,FUN=function(x){sum(x=="high")})))/length(de.unique.event)

# Proportion of events where each rank accessed meat
de.unique.event <- unique(de$EventID)
low.and.accessed <- rep(NA,length(de.unique.event))
mid.and.accessed <- rep(NA,length(de.unique.event))
high.and.accessed <- rep(NA,length(de.unique.event))
for(i in 1:length(de.unique.event)){
	dt <- subset(de,EventID==de.unique.event[i])
	index.temp <- which(dt$rank=="high" & dt$meat=="yes")
	high.and.accessed[i] <- 0*(length(index.temp)==0)+1*(length(index.temp)!=0)
	index.temp <- which(dt$rank=="middle" & dt$meat=="yes")
	mid.and.accessed[i] <- 0*(length(index.temp)==0)+1*(length(index.temp)!=0)
	index.temp <- which(dt$rank=="low" & dt$meat=="yes")
	low.and.accessed[i] <- 0*(length(index.temp)==0)+1*(length(index.temp)!=0)
}

percent.at.and.accessed <- data.frame(Event=c("Present","","","Accessed","",""),rank=rep(c("Low","Middle","High"),2),
	Percent=c(de.low.per,de.mid.per,de.high.per,sum(low.and.accessed)/length(de.unique.event),sum(mid.and.accessed)/length(de.unique.event),sum(high.and.accessed)/length(de.unique.event)))

write.csv(percent.at.and.accessed,
          file="~/Apes/Julia/newdata/Table-meat-presence-and-access-by-rank-data-summary-Jan-2019.csv",row.names=F)

save(list=ls(),file="~/Apes/Julia/newdata/meat-eating-models-Jan-2019.RData")

#load("~/Apes/Julia/newdata/meat-eating-models-Jan-2019.RData)
# Diagnostics
ranef.diagn.plot<-function(model.res, QQ=F){
  old.par = par(no.readonly = TRUE)
	n.plots=sum(unlist(lapply(ranef(model.res), length)))
	x=ifelse(n.plots%%2==1,n.plots+1,n.plots)
	xmat=outer(1:x, 1:(1+x/2), "*")-n.plots
	colnames(xmat)=1:(1+x/2)
	rownames(xmat)=1:x
	xmat=as.data.frame(as.table(xmat))
	xmat=subset(xmat, as.numeric(as.character(xmat$Var1))>=as.numeric(as.character(xmat$Var2)) & xmat$Freq>=0)
	sum.diff=as.numeric(as.character(xmat$Var1))-as.numeric(as.character(xmat$Var2))+xmat$Freq
	xmat=xmat[sum.diff==min(sum.diff),]
	xmat=xmat[which.min(xmat$Freq),]
	par(mfrow=c(xmat$Var2, xmat$Var1))
	par(mar=c(rep(2, 3), 1))
	par(mgp=c(1, 0.5, 0))
	for(i in 1:length(ranef(model.res))){
		to.plot=ranef(model.res)[[i]]
		for(k in 1:ncol(to.plot)){
			if(QQ){
				qqnorm(to.plot[,k], main="", tcl=-0.25, xlab="", ylab="", pch=19, col=grey(0.5, alpha=0.75))
				qqline(to.plot[,k])
			}else{
				hist(to.plot[,k], main="", tcl=-0.25, xlab="", ylab="")
			}
			mtext(text=paste(c(names(ranef(model.res)[i]), colnames(to.plot)[k]), collapse=", "), side=3)
		}
	}
	par(old.par)
}
diagnostics.plot<-function(mod.res, col=grey(level=0.25, alpha=0.5)){
  old.par = par(no.readonly = TRUE)
  par(mfrow=c(2, 2))
  par(mar=c(3, 3, 1, 0.5))
  hist(residuals(mod.res), probability=T, xlab="", ylab="", main="")
  mtext(text="histogram of residuals", side=3, line=0)
  x=seq(min(residuals(mod.res)), max(residuals(mod.res)), length.out=100)
  lines(x, dnorm(x, mean=0, sd=sd(residuals(mod.res))))
  qqnorm(residuals(mod.res), main="", pch=19)
  qqline(residuals(mod.res))
  mtext(text="qq-plot of residuals", side=3, line=0)
  plot(fitted(mod.res), residuals(mod.res), pch=19, col=col)
  abline(h=0, lty=2)
  mtext(text="residuals against fitted values", side=3, line=0)
  par(old.par)
}
ranef.diagn.plot(pm.no.inter)
diagnostics.plot(pm.no.inter)  # Probably not useful for logistic regression