## xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
# R code for analyses in:  
# Aggressive signaling among competing species of birds
## xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

signal.dat<-read.csv(file.choose(), stringsAsFactors = FALSE,strip.white = TRUE, na.strings = c("NA","") )#use dialog box to import data file "aggressive_signals.csv"

names(signal.dat) # lists all variable names



############################################################
####Postures used in aggressive interspecific encounters####
############################################################

### summarize data for all
summary(signal.dat$substrate) #not considered a component of posture, but used in description of data
summary(signal.dat$body.orientation)
summary(signal.dat$head.position)
summary(signal.dat$wing.position)
summary(signal.dat$shoulder.position)
summary(signal.dat$bill.position)
summary(signal.dat$tail.position)
summary(signal.dat$feet.position)
summary(signal.dat$closest.point)

#random effects
str(signal.dat$animal)
str(signal.dat$focal.interaction)
signal.dat$focal.interaction<-as.factor(signal.dat$focal.interaction) # make factor

library(dplyr)
library(ggplot2)


#####posture: body.orientation####

##make long dataset for MCMCglmm analysis
names(signal.dat)
body.orientation.dat<-signal.dat[,c(1:5,12,55,56)]
levels(as.factor(body.orientation.dat$body.orientation))
#create new binary columns for each body orientation category
body.orientation.dat$above<-NA
body.orientation.dat$feet.forward<-NA
body.orientation.dat$aa.forward.lowered<-NA
body.orientation.dat$forward.normal<-NA
body.orientation.dat$forward.upright<-NA
body.orientation.dat$side.oriented<-NA
body.orientation.dat$upside.down<-NA
for(i in 1:nrow(body.orientation.dat)){
  if(body.orientation.dat$body.orientation[i]=="above"){body.orientation.dat$above[i]<-1} else{body.orientation.dat$above[i]<-0}
  if(body.orientation.dat$body.orientation[i]=="feet.forward"){body.orientation.dat$feet.forward[i]<-1} else{body.orientation.dat$feet.forward[i]<-0}
  if(body.orientation.dat$body.orientation[i]=="forward.lowered"){body.orientation.dat$aa.forward.lowered[i]<-1} else{body.orientation.dat$aa.forward.lowered[i]<-0}
  if(body.orientation.dat$body.orientation[i]=="forward.normal"){body.orientation.dat$forward.normal[i]<-1} else{body.orientation.dat$forward.normal[i]<-0}
  if(body.orientation.dat$body.orientation[i]=="forward.upright"){body.orientation.dat$forward.upright[i]<-1} else{body.orientation.dat$forward.upright[i]<-0}
  if(body.orientation.dat$body.orientation[i]=="side.oriented"){body.orientation.dat$side.oriented[i]<-1} else{body.orientation.dat$side.oriented[i]<-0}
  if(body.orientation.dat$body.orientation[i]=="upside.down"){body.orientation.dat$upside.down[i]<-1} else{body.orientation.dat$upside.down[i]<-0}
}

#makde dataset long for analysis
library(tidyr)
names(body.orientation.dat)
body.orientation.long<-gather(body.orientation.dat,body.orientation,posture.assumed,c(9:15))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
body.orientation.long$posture.assumed.bv[body.orientation.long$posture.assumed > 0] <- "yes"
body.orientation.long$posture.assumed.bv[body.orientation.long$posture.assumed == 0] <- "no"
body.orientation.long$posture.assumed.bv <- as.factor(body.orientation.long$posture.assumed.bv)

k <- length(levels(body.orientation.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
body.orientation.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ body.orientation, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = body.orientation.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
body.orientation.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ body.orientation, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = body.orientation.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
body.orientation.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ body.orientation, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = body.orientation.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
body.orientation.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = body.orientation.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(body.orientation.model.MCMCglmm$Sol,body.orientation.model.MCMCglmm.2$Sol, body.orientation.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(body.orientation.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(body.orientation.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(body.orientation.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- body.orientation.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(body.orientation.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
body.orientation.model.MCMCglmm$DIC
body.orientation.model.MCMCglmm.null$DIC

summary(body.orientation.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#forward.lowered
forward.lowered.est<-inv.logit(mean(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])))
forward.lowered.est
forward.lowered.ci<-inv.logit(HPDinterval(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])))
forward.lowered.ci
#forward.upright
forward.upright.est<-inv.logit(mean(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"]) + body.orientation.model.MCMCglmm$Sol[,"body.orientationforward.upright"]))
forward.upright.est
forward.upright.ci<-inv.logit(HPDinterval(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])+ body.orientation.model.MCMCglmm$Sol[,"body.orientationforward.upright"])) 
forward.upright.ci
#forward.normal
forward.normal.est<-forward.normal.est<-inv.logit(mean(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"]) + body.orientation.model.MCMCglmm$Sol[,"body.orientationforward.normal"])) 
forward.normal.est
forward.normal.ci<-inv.logit(HPDinterval(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])+ body.orientation.model.MCMCglmm$Sol[,"body.orientationforward.normal"]))
forward.normal.ci
#feet.forward
feet.forward.est<-inv.logit(mean(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"]) + body.orientation.model.MCMCglmm$Sol[,"body.orientationfeet.forward"])) 
feet.forward.est
feet.forward.ci<-inv.logit(HPDinterval(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])+ body.orientation.model.MCMCglmm$Sol[,"body.orientationfeet.forward"])) 
feet.forward.ci
#side.oriented
side.oriented.est<-inv.logit(mean(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"]) + body.orientation.model.MCMCglmm$Sol[,"body.orientationside.oriented"])) 
side.oriented.est
side.oriented.ci<-inv.logit(HPDinterval(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])+ body.orientation.model.MCMCglmm$Sol[,"body.orientationside.oriented"])) 
side.oriented.ci
#above
above.est<-inv.logit(mean(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"]) + body.orientation.model.MCMCglmm$Sol[,"body.orientationabove"]))
above.est
above.ci<-inv.logit(HPDinterval(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])+ body.orientation.model.MCMCglmm$Sol[,"body.orientationabove"])) 
above.ci
#upside.down
upside.down.est<-inv.logit(mean(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"]) + body.orientation.model.MCMCglmm$Sol[,"body.orientationupside.down"])) 
upside.down.est
upside.down.ci<-inv.logit(HPDinterval(mcmc(body.orientation.model.MCMCglmm$Sol[,"(Intercept)"])+ body.orientation.model.MCMCglmm$Sol[,"body.orientationupside.down"])) 
upside.down.ci

##insert figure
l.body.orientation <- list(
  a = list(body.orientation2= "aa.forward.lowered", posture.assumed = forward.lowered.est, lower = forward.lowered.ci[1] ,upper = forward.lowered.ci[2]) 
  , b = list(body.orientation2= "b.forward.upright", posture.assumed =forward.upright.est, lower = forward.upright.ci[1],upper =  forward.upright.ci[2])
  , c = list(body.orientation2= "c.forward.normal", posture.assumed = forward.normal.est, lower=forward.normal.ci[1], upper= forward.normal.ci[2])
  , d = list(body.orientation2= "d.feet.forward", posture.assumed = feet.forward.est, lower=feet.forward.ci[1], upper= feet.forward.ci[2])
  , e = list(body.orientation2= "e.side.oriented", posture.assumed = side.oriented.est, lower=side.oriented.ci[1], upper= side.oriented.ci[2])
  , f = list(body.orientation2= "f.above", posture.assumed = above.est, lower=above.ci[1], upper= above.ci[2])
  , g = list(body.orientation2= "g.upside.down", posture.assumed = upside.down.est, lower=upside.down.ci[1], upper= upside.down.ci[2])
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
body.orientation.df <- ldply (l.body.orientation, data.frame) #create dataframe with estimates
detach(package:plyr)

body.orientation.long["body.orientation2"]<-NA #new column to adjust order of x-axis
for (i in 1:nrow(body.orientation.long)){
  if (body.orientation.long$body.orientation[i]=="aa.forward.lowered") {body.orientation.long$body.orientation2[i]<-"aa.forward.lowered"}
  if (body.orientation.long$body.orientation[i]=="forward.upright") {body.orientation.long$body.orientation2[i]<-"b.forward.upright"}
  if (body.orientation.long$body.orientation[i]=="forward.normal") {body.orientation.long$body.orientation2[i]<-"c.forward.normal"}
  if (body.orientation.long$body.orientation[i]=="feet.forward") {body.orientation.long$body.orientation2[i]<-"d.feet.forward"}
  if (body.orientation.long$body.orientation[i]=="side.oriented") {body.orientation.long$body.orientation2[i]<-"e.side.oriented"}
  if (body.orientation.long$body.orientation[i]=="above") {body.orientation.long$body.orientation2[i]<-"f.above"}
  if (body.orientation.long$body.orientation[i]=="upside.down") {body.orientation.long$body.orientation2[i]<-"g.upside.down"}
}

library(ggplot2)
##figure
p<-ggplot(data=body.orientation.long,aes(body.orientation2,posture.assumed))
p2<-p+geom_jitter(data=body.orientation.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=body.orientation2,y=posture.assumed),data=body.orientation.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=body.orientation2,ymin=lower,ymax=upper),width=0.1,data=body.orientation.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)



#####posture: head.position####

##make long dataset for MCMCglmm analysis
names(signal.dat)
head.position.dat<-signal.dat[,c(1:5,13,55,56)]
levels(as.factor(head.position.dat$head.position))
head.position.dat<-head.position.dat[-c(head.position.dat$head.position=="not.visible"),] #remove interactions in which head position was not visible
nrow(head.position.dat) #336 interactions remain

#create new binary columns for each body orientation category
head.position.dat$aa.forward.lowered<-NA
head.position.dat$b.forward.upright<-NA
head.position.dat$c.forward.normal<-NA
head.position.dat$d.side.oriented<-NA
head.position.dat$e.held.back.upward<-NA

for(i in 1:nrow(head.position.dat)){
  if(head.position.dat$head.position[i]=="forward.lowered"){head.position.dat$aa.forward.lowered[i]<-1} else{head.position.dat$aa.forward.lowered[i]<-0}
  if(head.position.dat$head.position[i]=="forward.upright"){head.position.dat$b.forward.upright[i]<-1} else{head.position.dat$b.forward.upright[i]<-0}
  if(head.position.dat$head.position[i]=="forward.normal"){head.position.dat$c.forward.normal[i]<-1} else{head.position.dat$c.forward.normal[i]<-0}
  if(head.position.dat$head.position[i]=="side.oriented"){head.position.dat$d.side.oriented[i]<-1} else{head.position.dat$d.side.oriented[i]<-0}
  if(head.position.dat$head.position[i]=="held.back.upward"){head.position.dat$e.held.back.upward[i]<-1} else{head.position.dat$e.held.back.upward[i]<-0}
}

#makde dataset long for analysis
library(tidyr)
names(head.position.dat)
head.position.long<-gather(head.position.dat,head.position,posture.assumed,c(9:13))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose())#load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
head.position.long$posture.assumed.bv[head.position.long$posture.assumed > 0] <- "yes"
head.position.long$posture.assumed.bv[head.position.long$posture.assumed == 0] <- "no"
head.position.long$posture.assumed.bv <- as.factor(head.position.long$posture.assumed.bv)

k <- length(levels(head.position.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
head.position.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ head.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = head.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
head.position.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ head.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = head.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
head.position.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ head.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = head.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
head.position.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = head.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(head.position.model.MCMCglmm$Sol,head.position.model.MCMCglmm.2$Sol, head.position.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(head.position.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(head.position.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(head.position.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- head.position.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(head.position.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
head.position.model.MCMCglmm$DIC
head.position.model.MCMCglmm.null$DIC

summary(head.position.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#forward.lowered
forward.lowered.est<-inv.logit(mean(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
forward.lowered.est
forward.lowered.ci<-inv.logit(HPDinterval(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
forward.lowered.ci
#forward.upright
forward.upright.est<-inv.logit(mean(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"]) + head.position.model.MCMCglmm$Sol[,"head.positionb.forward.upright"])) 
forward.upright.est
forward.upright.ci<-inv.logit(HPDinterval(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"])+ head.position.model.MCMCglmm$Sol[,"head.positionb.forward.upright"]))
forward.upright.ci
#forward.normal
forward.normal.est<-forward.normal.est<-inv.logit(mean(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"]) + head.position.model.MCMCglmm$Sol[,"head.positionc.forward.normal"])) 
forward.normal.est
forward.normal.ci<-inv.logit(HPDinterval(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"])+ head.position.model.MCMCglmm$Sol[,"head.positionc.forward.normal"]))
forward.normal.ci
#side.oriented
side.oriented.est<-inv.logit(mean(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"]) + head.position.model.MCMCglmm$Sol[,"head.positiond.side.oriented"]))
side.oriented.est
side.oriented.ci<-inv.logit(HPDinterval(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"])+ head.position.model.MCMCglmm$Sol[,"head.positiond.side.oriented"]))
side.oriented.ci
#held.back.upward
held.back.upward.est<-inv.logit(mean(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"]) + head.position.model.MCMCglmm$Sol[,"head.positione.held.back.upward"])) 
held.back.upward.est
held.back.upward.ci<-inv.logit(HPDinterval(mcmc(head.position.model.MCMCglmm$Sol[,"(Intercept)"])+ head.position.model.MCMCglmm$Sol[,"head.positione.held.back.upward"]))
held.back.upward.ci

##insert figure
l.head.position <- list(
  a = list(head.position= "aa.forward.lowered", posture.assumed = forward.lowered.est, lower = forward.lowered.ci[1] ,upper = forward.lowered.ci[2]) 
  , b = list(head.position= "b.forward.upright", posture.assumed =forward.upright.est, lower = forward.upright.ci[1],upper =  forward.upright.ci[2])
  , c = list(head.position= "c.forward.normal", posture.assumed = forward.normal.est, lower=forward.normal.ci[1], upper= forward.normal.ci[2])
  , d = list(head.position= "d.side.oriented", posture.assumed = side.oriented.est, lower=side.oriented.ci[1], upper= side.oriented.ci[2])
  , e = list(head.position= "e.held.back.upward", posture.assumed = held.back.upward.est, lower=held.back.upward.ci[1], upper= held.back.upward.ci[2])
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
head.position.df <- ldply (l.head.position, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=head.position.long,aes(head.position,posture.assumed))
p2<-p+geom_jitter(data=head.position.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=head.position,y=posture.assumed),data=head.position.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=head.position,ymin=lower,ymax=upper),width=0.1,data=head.position.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)



####posture: wing.position####

##make long dataset for MCMCglmm analysis
names(signal.dat)
wing.position.dat<-signal.dat[,c(1:5,14,55,56)]
levels(as.factor(wing.position.dat$wing.position))
nrow(wing.position.dat) #337 interactions

#create new binary columns for each body orientation category
wing.position.dat$aa.closed.held.slightly.out<-NA
wing.position.dat$b.closed.flat<-NA
wing.position.dat$c.closed.raised.off.back<-NA
wing.position.dat$d.flapping<-NA
wing.position.dat$e.spread.outward<-NA
wing.position.dat$f.partially.spread<-NA
wing.position.dat$g.raised.upward<-NA
wing.position.dat$h.soaring.gliding<-NA

for(i in 1:nrow(wing.position.dat)){
  if(wing.position.dat$wing.position[i]=="closed.held.slightly.out"){wing.position.dat$aa.closed.held.slightly.out[i]<-1} else{wing.position.dat$aa.closed.held.slightly.out[i]<-0}
  if(wing.position.dat$wing.position[i]=="closed.flat"){wing.position.dat$b.closed.flat[i]<-1} else{wing.position.dat$b.closed.flat[i]<-0}
  if(wing.position.dat$wing.position[i]=="closed.raised.off.back"){wing.position.dat$c.closed.raised.off.back[i]<-1} else{wing.position.dat$c.closed.raised.off.back[i]<-0}
  if(wing.position.dat$wing.position[i]=="flapping"){wing.position.dat$d.flapping[i]<-1} else{wing.position.dat$d.flapping[i]<-0}
  if(wing.position.dat$wing.position[i]=="spread.outward"){wing.position.dat$e.spread.outward[i]<-1} else{wing.position.dat$e.spread.outward[i]<-0}
  if(wing.position.dat$wing.position[i]=="partially.spread"){wing.position.dat$f.partially.spread[i]<-1} else{wing.position.dat$f.partially.spread[i]<-0}
  if(wing.position.dat$wing.position[i]=="raised.upward"){wing.position.dat$g.raised.upward[i]<-1} else{wing.position.dat$g.raised.upward[i]<-0}
  if(wing.position.dat$wing.position[i]=="soaring.gliding"){wing.position.dat$h.soaring.gliding[i]<-1} else{wing.position.dat$h.soaring.gliding[i]<-0}
  
}

#makde dataset long for analysis
library(tidyr)
names(wing.position.dat)
wing.position.long<-gather(wing.position.dat,wing.position,posture.assumed,c(9:ncol(wing.position.dat)))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
wing.position.long$posture.assumed.bv[wing.position.long$posture.assumed > 0] <- "yes"
wing.position.long$posture.assumed.bv[wing.position.long$posture.assumed == 0] <- "no"
wing.position.long$posture.assumed.bv <- as.factor(wing.position.long$posture.assumed.bv)

k <- length(levels(wing.position.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
wing.position.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ wing.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
wing.position.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ wing.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
wing.position.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ wing.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
wing.position.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(wing.position.model.MCMCglmm$Sol,wing.position.model.MCMCglmm.2$Sol, wing.position.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(wing.position.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(wing.position.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(wing.position.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- wing.position.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(wing.position.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
wing.position.model.MCMCglmm$DIC
wing.position.model.MCMCglmm.null$DIC

summary(wing.position.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)

#closed.held.slightly.out
closed.held.slightly.out.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
closed.held.slightly.out.est
closed.held.slightly.out.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
closed.held.slightly.out.ci
#closed.flat
closed.flat.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]) + wing.position.model.MCMCglmm$Sol[,"wing.positionb.closed.flat"])) 
closed.flat.est
closed.flat.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"])+ wing.position.model.MCMCglmm$Sol[,"wing.positionb.closed.flat"]))
closed.flat.ci
#closed.raised.off.back
closed.raised.off.back.est<-closed.raised.off.back.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]) + wing.position.model.MCMCglmm$Sol[,"wing.positionc.closed.raised.off.back"]))
closed.raised.off.back.est
closed.raised.off.back.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"])+ wing.position.model.MCMCglmm$Sol[,"wing.positionc.closed.raised.off.back"]))
closed.raised.off.back.ci
#flapping
flapping.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]) + wing.position.model.MCMCglmm$Sol[,"wing.positiond.flapping"]))
flapping.est
flapping.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"])+ wing.position.model.MCMCglmm$Sol[,"wing.positiond.flapping"]))
flapping.ci
#spread.outward
spread.outward.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]) + wing.position.model.MCMCglmm$Sol[,"wing.positione.spread.outward"])) 
spread.outward.est
spread.outward.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"])+ wing.position.model.MCMCglmm$Sol[,"wing.positione.spread.outward"]))
spread.outward.ci
#partially.spread
partially.spread.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]) + wing.position.model.MCMCglmm$Sol[,"wing.positionf.partially.spread"]))
partially.spread.est
partially.spread.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"])+ wing.position.model.MCMCglmm$Sol[,"wing.positionf.partially.spread"]))
partially.spread.ci
#raised.upward
raised.upward.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]) + wing.position.model.MCMCglmm$Sol[,"wing.positiong.raised.upward"]))
raised.upward.est
raised.upward.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"])+ wing.position.model.MCMCglmm$Sol[,"wing.positiong.raised.upward"])) 
raised.upward.ci
#soaring.gliding
soaring.gliding.est<-inv.logit(mean(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"]) + wing.position.model.MCMCglmm$Sol[,"wing.positionh.soaring.gliding"])) 
soaring.gliding.est
soaring.gliding.ci<-inv.logit(HPDinterval(mcmc(wing.position.model.MCMCglmm$Sol[,"(Intercept)"])+ wing.position.model.MCMCglmm$Sol[,"wing.positionh.soaring.gliding"])) 
soaring.gliding.ci

##insert figure
l.wing.position <- list(
  a = list(wing.position= "aa.closed.held.slightly.out", posture.assumed = closed.held.slightly.out.est, lower = closed.held.slightly.out.ci[1] ,upper = closed.held.slightly.out.ci[2]) 
  , b = list(wing.position= "b.closed.flat", posture.assumed =closed.flat.est, lower = closed.flat.ci[1],upper =  closed.flat.ci[2])
  , c = list(wing.position= "c.closed.raised.off.back", posture.assumed = closed.raised.off.back.est, lower=closed.raised.off.back.ci[1], upper= closed.raised.off.back.ci[2])
  , d = list(wing.position= "d.flapping", posture.assumed = flapping.est, lower=flapping.ci[1], upper= flapping.ci[2])
  , e = list(wing.position= "e.spread.outward", posture.assumed = spread.outward.est, lower=spread.outward.ci[1], upper= spread.outward.ci[2])
  , f = list(wing.position= "f.partially.spread", posture.assumed = partially.spread.est, lower=partially.spread.ci[1], upper= partially.spread.ci[2])
  , g = list(wing.position= "g.raised.upward", posture.assumed = raised.upward.est, lower=raised.upward.ci[1], upper= raised.upward.ci[2])
  , h = list(wing.position= "h.soaring.gliding", posture.assumed = soaring.gliding.est, lower=soaring.gliding.ci[1], upper= soaring.gliding.ci[2])
  
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
wing.position.df <- ldply (l.wing.position, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=wing.position.long,aes(wing.position,posture.assumed))
p2<-p+geom_jitter(data=wing.position.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=wing.position,y=posture.assumed),data=wing.position.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=wing.position,ymin=lower,ymax=upper),width=0.1,data=wing.position.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)


####wing.position.2: compare closed vs open wings####
wing.position.dat$wings.open<-NA
wing.position.dat$wings.closed<-NA
for(i in 1:nrow(wing.position.dat)){
  if(wing.position.dat$wing.position[i]=="closed.held.slightly.out"){(wing.position.dat$wings.open[i]<-0)&(wing.position.dat$wings.closed[i]<-1)}
  if(wing.position.dat$wing.position[i]=="closed.flat"){(wing.position.dat$wings.open[i]<-0)&(wing.position.dat$wings.closed[i]<-1)}
  if(wing.position.dat$wing.position[i]=="closed.raised.off.back"){(wing.position.dat$wings.open[i]<-0)&(wing.position.dat$wings.closed[i]<-1)}
  if(wing.position.dat$wing.position[i]=="flapping"){(wing.position.dat$wings.open[i]<-1)&(wing.position.dat$wings.closed[i]<-0)}
  if(wing.position.dat$wing.position[i]=="spread.outward"){(wing.position.dat$wings.open[i]<-1)&(wing.position.dat$wings.closed[i]<-0)}
  if(wing.position.dat$wing.position[i]=="partially.spread"){(wing.position.dat$wings.open[i]<-1)&(wing.position.dat$wings.closed[i]<-0)}
  if(wing.position.dat$wing.position[i]=="raised.upward"){(wing.position.dat$wings.open[i]<-1)&(wing.position.dat$wings.closed[i]<-0)}
  if(wing.position.dat$wing.position[i]=="soaring.gliding"){(wing.position.dat$wings.open[i]<-1)&(wing.position.dat$wings.closed[i]<-0)}
}

library(tidyr)
names(wing.position.dat)
wing.position.dat.sub<-wing.position.dat[,c(1:8,17,18)]
wing.position.long2<-gather(wing.position.dat.sub,wings.open,posture.assumed,c(9,10))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
wing.position.long2$posture.assumed.bv[wing.position.long2$posture.assumed > 0] <- "yes"
wing.position.long2$posture.assumed.bv[wing.position.long2$posture.assumed == 0] <- "no"
wing.position.long2$posture.assumed.bv <- as.factor(wing.position.long2$posture.assumed.bv)

k <- length(levels(wing.position.long2$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
wings.open.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ wings.open, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long2, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
wings.open.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ wings.open, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long2, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
wings.open.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ wings.open, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long2, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
wings.open.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = wing.position.long2, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(wings.open.model.MCMCglmm$Sol,wings.open.model.MCMCglmm.2$Sol, wings.open.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(wings.open.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(wings.open.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(wings.open.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- wings.open.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(wings.open.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
wings.open.model.MCMCglmm$DIC
wings.open.model.MCMCglmm.null$DIC

summary(wings.open.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)

#wings.closed
wings.closed.est<-inv.logit(mean(mcmc(wings.open.model.MCMCglmm$Sol[,"(Intercept)"]))) 
wings.closed.est
wings.closed.ci<-inv.logit(HPDinterval(mcmc(wings.open.model.MCMCglmm$Sol[,"(Intercept)"]))) 
wings.closed.ci
#wings.open
wings.open.est<-inv.logit(mean(mcmc(wings.open.model.MCMCglmm$Sol[,"(Intercept)"]) + wings.open.model.MCMCglmm$Sol[,"wings.openwings.open"])) 
wings.open.est
wings.open.ci<-inv.logit(HPDinterval(mcmc(wings.open.model.MCMCglmm$Sol[,"(Intercept)"])+ wings.open.model.MCMCglmm$Sol[,"wings.openwings.open"]))
wings.open.ci



####posture: shoulder.position####

##make long dataset for MCMCglmm analysis
names(signal.dat)
shoulder.position.dat<-signal.dat[,c(1:5,15,55,56)]
levels(as.factor(shoulder.position.dat$shoulder.position))
shoulder.position.dat<-shoulder.position.dat[-c(shoulder.position.dat$shoulder.position=="not.visible"),] #remove interactions in which position was not visible
nrow(shoulder.position.dat) #336 interactions remain

#create new binary columns for each body orientation category
shoulder.position.dat$aa.underwing.forward<-NA
shoulder.position.dat$b.wing.closed.shoulder.visible<-NA
shoulder.position.dat$c.wing.horizontal<-NA
shoulder.position.dat$d.wing.closed.shoulder.concealed<-NA
shoulder.position.dat$e.upperwing.forward<-NA


for(i in 1:nrow(shoulder.position.dat)){
  if(shoulder.position.dat$shoulder.position[i]=="underwing.forward"){shoulder.position.dat$aa.underwing.forward[i]<-1} else{shoulder.position.dat$aa.underwing.forward[i]<-0}
  if(shoulder.position.dat$shoulder.position[i]=="wing.closed.shoulder.visible"){shoulder.position.dat$b.wing.closed.shoulder.visible[i]<-1} else{shoulder.position.dat$b.wing.closed.shoulder.visible[i]<-0}
  if(shoulder.position.dat$shoulder.position[i]=="wing.horizontal"){shoulder.position.dat$c.wing.horizontal[i]<-1} else{shoulder.position.dat$c.wing.horizontal[i]<-0}
  if(shoulder.position.dat$shoulder.position[i]=="wing.closed.shoulder.concealed"){shoulder.position.dat$d.wing.closed.shoulder.concealed[i]<-1} else{shoulder.position.dat$d.wing.closed.shoulder.concealed[i]<-0}
  if(shoulder.position.dat$shoulder.position[i]=="upperwing.forward"){shoulder.position.dat$e.upperwing.forward[i]<-1} else{shoulder.position.dat$e.upperwing.forward[i]<-0}
  
}

#makde dataset long for analysis
library(tidyr)
names(shoulder.position.dat)
shoulder.position.long<-gather(shoulder.position.dat,shoulder.position,posture.assumed,c(9:13))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
shoulder.position.long$posture.assumed.bv[shoulder.position.long$posture.assumed > 0] <- "yes"
shoulder.position.long$posture.assumed.bv[shoulder.position.long$posture.assumed == 0] <- "no"
shoulder.position.long$posture.assumed.bv <- as.factor(shoulder.position.long$posture.assumed.bv)

k <- length(levels(shoulder.position.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
shoulder.position.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ shoulder.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = shoulder.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
shoulder.position.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ shoulder.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = shoulder.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
shoulder.position.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ shoulder.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = shoulder.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
shoulder.position.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = shoulder.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(shoulder.position.model.MCMCglmm$Sol,shoulder.position.model.MCMCglmm.2$Sol, shoulder.position.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(shoulder.position.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(shoulder.position.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(shoulder.position.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- shoulder.position.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(shoulder.position.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
shoulder.position.model.MCMCglmm$DIC
shoulder.position.model.MCMCglmm.null$DIC

summary(shoulder.position.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#underwing.forward
underwing.forward.est<-inv.logit(mean(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
underwing.forward.est
underwing.forward.ci<-inv.logit(HPDinterval(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
underwing.forward.ci
#wing.closed.shoulder.visible
wing.closed.shoulder.visible.est<-inv.logit(mean(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"]) + shoulder.position.model.MCMCglmm$Sol[,"shoulder.positionb.wing.closed.shoulder.visible"])) 
wing.closed.shoulder.visible.est
wing.closed.shoulder.visible.ci<-inv.logit(HPDinterval(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"])+ shoulder.position.model.MCMCglmm$Sol[,"shoulder.positionb.wing.closed.shoulder.visible"]))
wing.closed.shoulder.visible.ci
#wing.horizontal
wing.horizontal.est<-wing.horizontal.est<-inv.logit(mean(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"]) + shoulder.position.model.MCMCglmm$Sol[,"shoulder.positionc.wing.horizontal"])) 
wing.horizontal.est
wing.horizontal.ci<-inv.logit(HPDinterval(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"])+ shoulder.position.model.MCMCglmm$Sol[,"shoulder.positionc.wing.horizontal"]))
wing.horizontal.ci
#wing.closed.shoulder.concealed
wing.closed.shoulder.concealed.est<-inv.logit(mean(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"]) + shoulder.position.model.MCMCglmm$Sol[,"shoulder.positiond.wing.closed.shoulder.concealed"])) 
wing.closed.shoulder.concealed.est
wing.closed.shoulder.concealed.ci<-inv.logit(HPDinterval(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"])+ shoulder.position.model.MCMCglmm$Sol[,"shoulder.positiond.wing.closed.shoulder.concealed"]))
wing.closed.shoulder.concealed.ci
#upperwing.forward
upperwing.forward.est<-inv.logit(mean(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"]) + shoulder.position.model.MCMCglmm$Sol[,"shoulder.positione.upperwing.forward"])) 
upperwing.forward.est
upperwing.forward.ci<-inv.logit(HPDinterval(mcmc(shoulder.position.model.MCMCglmm$Sol[,"(Intercept)"])+ shoulder.position.model.MCMCglmm$Sol[,"shoulder.positione.upperwing.forward"]))
upperwing.forward.ci


##insert figure
l.shoulder.position <- list(
  a = list(shoulder.position= "aa.underwing.forward", posture.assumed = underwing.forward.est, lower = underwing.forward.ci[1] ,upper = underwing.forward.ci[2]) 
  , b = list(shoulder.position= "b.wing.closed.shoulder.visible", posture.assumed =wing.closed.shoulder.visible.est, lower = wing.closed.shoulder.visible.ci[1],upper =  wing.closed.shoulder.visible.ci[2])
  , c = list(shoulder.position= "c.wing.horizontal", posture.assumed = wing.horizontal.est, lower=wing.horizontal.ci[1], upper= wing.horizontal.ci[2])
  , d = list(shoulder.position= "d.wing.closed.shoulder.concealed", posture.assumed = wing.closed.shoulder.concealed.est, lower=wing.closed.shoulder.concealed.ci[1], upper= wing.closed.shoulder.concealed.ci[2])
  , e = list(shoulder.position= "e.upperwing.forward", posture.assumed = upperwing.forward.est, lower=upperwing.forward.ci[1], upper= upperwing.forward.ci[2])
 
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
shoulder.position.df <- ldply (l.shoulder.position, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=shoulder.position.long,aes(shoulder.position,posture.assumed))
p2<-p+geom_jitter(data=shoulder.position.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=shoulder.position,y=posture.assumed),data=shoulder.position.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=shoulder.position,ymin=lower,ymax=upper),width=0.1,data=shoulder.position.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)



####posture: bill.position####

##make long dataset for MCMCglmm analysis
names(signal.dat)
bill.position.dat<-signal.dat[,c(1:5,16,55,56)]
levels(as.factor(bill.position.dat$bill.position))
bill.position.dat<-bill.position.dat[-c(bill.position.dat$bill.position=="not.visible"),] #remove interactions in which position was not visible
nrow(bill.position.dat) #336 interactions remain

#create new binary columns for each body orientation category
bill.position.dat$aa.open.forward<-NA
bill.position.dat$b.closed.forward<-NA
bill.position.dat$c.open.side<-NA
bill.position.dat$d.closed.side<-NA
bill.position.dat$e.closed.downward<-NA
bill.position.dat$f.open.downward<-NA
bill.position.dat$g.closed.upward<-NA

for(i in 1:nrow(bill.position.dat)){
  if(bill.position.dat$bill.position[i]=="open.forward"){bill.position.dat$aa.open.forward[i]<-1} else{bill.position.dat$aa.open.forward[i]<-0}
  if(bill.position.dat$bill.position[i]=="closed.forward"){bill.position.dat$b.closed.forward[i]<-1} else{bill.position.dat$b.closed.forward[i]<-0}
  if(bill.position.dat$bill.position[i]=="open.side"){bill.position.dat$c.open.side[i]<-1} else{bill.position.dat$c.open.side[i]<-0}
  if(bill.position.dat$bill.position[i]=="closed.side"){bill.position.dat$d.closed.side[i]<-1} else{bill.position.dat$d.closed.side[i]<-0}
  if(bill.position.dat$bill.position[i]=="closed.downward"){bill.position.dat$e.closed.downward[i]<-1} else{bill.position.dat$e.closed.downward[i]<-0}
  if(bill.position.dat$bill.position[i]=="open.downward"){bill.position.dat$f.open.downward[i]<-1} else{bill.position.dat$f.open.downward[i]<-0}
  if(bill.position.dat$bill.position[i]=="closed.upward"){bill.position.dat$g.closed.upward[i]<-1} else{bill.position.dat$g.closed.upward[i]<-0}
}

#makde dataset long for analysis
library(tidyr)
names(bill.position.dat)
bill.position.long<-gather(bill.position.dat,bill.position,posture.assumed,c(9:15))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
bill.position.long$posture.assumed.bv[bill.position.long$posture.assumed > 0] <- "yes"
bill.position.long$posture.assumed.bv[bill.position.long$posture.assumed == 0] <- "no"
bill.position.long$posture.assumed.bv <- as.factor(bill.position.long$posture.assumed.bv)

k <- length(levels(bill.position.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
bill.position.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ bill.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = bill.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
bill.position.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ bill.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = bill.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
bill.position.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ bill.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = bill.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
bill.position.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = bill.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(bill.position.model.MCMCglmm$Sol,bill.position.model.MCMCglmm.2$Sol, bill.position.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(bill.position.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(bill.position.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(bill.position.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- bill.position.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(bill.position.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
bill.position.model.MCMCglmm$DIC
bill.position.model.MCMCglmm.null$DIC

summary(bill.position.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)

#open.forward
open.forward.est<-inv.logit(mean(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
open.forward.est
open.forward.ci<-inv.logit(HPDinterval(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
open.forward.ci
#closed.forward
closed.forward.est<-inv.logit(mean(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]) + bill.position.model.MCMCglmm$Sol[,"bill.positionb.closed.forward"])) 
closed.forward.est
closed.forward.ci<-inv.logit(HPDinterval(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"])+ bill.position.model.MCMCglmm$Sol[,"bill.positionb.closed.forward"]))
closed.forward.ci
#open.side
open.side.est<-open.side.est<-inv.logit(mean(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]) + bill.position.model.MCMCglmm$Sol[,"bill.positionc.open.side"])) 
open.side.est
open.side.ci<-inv.logit(HPDinterval(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"])+ bill.position.model.MCMCglmm$Sol[,"bill.positionc.open.side"])) 
open.side.ci
#closed.side
closed.side.est<-inv.logit(mean(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]) + bill.position.model.MCMCglmm$Sol[,"bill.positiond.closed.side"]))
closed.side.est
closed.side.ci<-inv.logit(HPDinterval(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"])+ bill.position.model.MCMCglmm$Sol[,"bill.positiond.closed.side"])) 
closed.side.ci
#closed.downward
closed.downward.est<-inv.logit(mean(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]) + bill.position.model.MCMCglmm$Sol[,"bill.positione.closed.downward"])) 
closed.downward.est
closed.downward.ci<-inv.logit(HPDinterval(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"])+ bill.position.model.MCMCglmm$Sol[,"bill.positione.closed.downward"])) 
closed.downward.ci
#open.downward
open.downward.est<-inv.logit(mean(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]) + bill.position.model.MCMCglmm$Sol[,"bill.positionf.open.downward"])) 
open.downward.est
open.downward.ci<-inv.logit(HPDinterval(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"])+ bill.position.model.MCMCglmm$Sol[,"bill.positionf.open.downward"])) 
open.downward.ci
#closed.upward
closed.upward.est<-inv.logit(mean(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"]) + bill.position.model.MCMCglmm$Sol[,"bill.positiong.closed.upward"])) 
closed.upward.est
closed.upward.ci<-inv.logit(HPDinterval(mcmc(bill.position.model.MCMCglmm$Sol[,"(Intercept)"])+ bill.position.model.MCMCglmm$Sol[,"bill.positiong.closed.upward"])) 
closed.upward.ci


##insert figure
l.bill.position <- list(
  a = list(bill.position= "aa.open.forward", posture.assumed = open.forward.est, lower = open.forward.ci[1] ,upper = open.forward.ci[2]) 
  , b = list(bill.position= "b.closed.forward", posture.assumed =closed.forward.est, lower = closed.forward.ci[1],upper =  closed.forward.ci[2])
  , c = list(bill.position= "c.open.side", posture.assumed = open.side.est, lower=open.side.ci[1], upper= open.side.ci[2])
  , d = list(bill.position= "d.closed.side", posture.assumed = closed.side.est, lower=closed.side.ci[1], upper= closed.side.ci[2])
  , e = list(bill.position= "e.closed.downward", posture.assumed = closed.downward.est, lower=closed.downward.ci[1], upper= closed.downward.ci[2])
  , f = list(bill.position= "f.open.downward", posture.assumed = open.downward.est, lower=open.downward.ci[1], upper= open.downward.ci[2])
  , g = list(bill.position= "g.closed.upward", posture.assumed = closed.upward.est, lower=closed.upward.ci[1], upper= closed.upward.ci[2])
  
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
bill.position.df <- ldply (l.bill.position, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=bill.position.long,aes(bill.position,posture.assumed))
p2<-p+geom_jitter(data=bill.position.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=bill.position,y=posture.assumed),data=bill.position.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=bill.position,ymin=lower,ymax=upper),width=0.1,data=bill.position.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)



####posture: tail.position####

##make long dataset for MCMCglmm analysis
names(signal.dat)
tail.position.dat<-signal.dat[,c(1:5,17,55,56)]
levels(as.factor(tail.position.dat$tail.position))
tail.position.dat<-tail.position.dat[-c(tail.position.dat$tail.position=="not.visible"),] #remove interactions in which position was not visible
nrow(tail.position.dat) #336 interactions remain

#create new binary columns for each body orientation category
tail.position.dat$aa.trailing.notfanned<-NA
tail.position.dat$b.trailing.fanned<-NA
tail.position.dat$c.down.fanned<-NA
tail.position.dat$d.down.notfanned<-NA
tail.position.dat$e.raised.fanned<-NA
tail.position.dat$f.partlyraised.notfanned<-NA
tail.position.dat$g.partlyraised.fanned<-NA
tail.position.dat$h.raised.notfanned<-NA
tail.position.dat$i.side.fanned<-NA
tail.position.dat$j.side.notfanned<-NA

for(i in 1:nrow(tail.position.dat)){
  if(tail.position.dat$tail.position[i]=="trailing.notfanned"){tail.position.dat$aa.trailing.notfanned[i]<-1} else{tail.position.dat$aa.trailing.notfanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="trailing.fanned"){tail.position.dat$b.trailing.fanned[i]<-1} else{tail.position.dat$b.trailing.fanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="down.fanned"){tail.position.dat$c.down.fanned[i]<-1} else{tail.position.dat$c.down.fanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="down.notfanned"){tail.position.dat$d.down.notfanned[i]<-1} else{tail.position.dat$d.down.notfanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="raised.fanned"){tail.position.dat$e.raised.fanned[i]<-1} else{tail.position.dat$e.raised.fanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="partlyraised.notfanned"){tail.position.dat$f.partlyraised.notfanned[i]<-1} else{tail.position.dat$f.partlyraised.notfanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="partlyraised.fanned"){tail.position.dat$g.partlyraised.fanned[i]<-1} else{tail.position.dat$g.partlyraised.fanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="raised.notfanned"){tail.position.dat$h.raised.notfanned[i]<-1} else{tail.position.dat$h.raised.notfanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="side.fanned"){tail.position.dat$i.side.fanned[i]<-1} else{tail.position.dat$i.side.fanned[i]<-0}
  if(tail.position.dat$tail.position[i]=="side.notfanned"){tail.position.dat$j.side.notfanned[i]<-1} else{tail.position.dat$j.side.notfanned[i]<-0}
  
}

#makde dataset long for analysis
library(tidyr)
names(tail.position.dat)
tail.position.long<-gather(tail.position.dat,tail.position,posture.assumed,c(9:ncol(tail.position.dat)))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
tail.position.long$posture.assumed.bv[tail.position.long$posture.assumed > 0] <- "yes"
tail.position.long$posture.assumed.bv[tail.position.long$posture.assumed == 0] <- "no"
tail.position.long$posture.assumed.bv <- as.factor(tail.position.long$posture.assumed.bv)

k <- length(levels(tail.position.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
tail.position.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ tail.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = tail.position.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
tail.position.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ tail.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = tail.position.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)
tail.position.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ tail.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = tail.position.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
tail.position.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = tail.position.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(tail.position.model.MCMCglmm$Sol,tail.position.model.MCMCglmm.2$Sol, tail.position.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(tail.position.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(tail.position.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(tail.position.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- tail.position.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(tail.position.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
tail.position.model.MCMCglmm$DIC
tail.position.model.MCMCglmm.null$DIC

summary(tail.position.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#trailing.notfanned
trailing.notfanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
trailing.notfanned.est
trailing.notfanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
trailing.notfanned.ci
#trailing.fanned
trailing.fanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positionb.trailing.fanned"])) 
trailing.fanned.est
trailing.fanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positionb.trailing.fanned"]))
trailing.fanned.ci
#down.fanned
down.fanned.est<-down.fanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positionc.down.fanned"]))
down.fanned.est
down.fanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positionc.down.fanned"]))
down.fanned.ci
#down.notfanned
down.notfanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positiond.down.notfanned"]))
down.notfanned.est
down.notfanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positiond.down.notfanned"]))
down.notfanned.ci
#raised.fanned
raised.fanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positione.raised.fanned"]))
raised.fanned.est
raised.fanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positione.raised.fanned"]))
raised.fanned.ci
#partlyraised.notfanned
partlyraised.notfanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positionf.partlyraised.notfanned"]))
partlyraised.notfanned.est
partlyraised.notfanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positionf.partlyraised.notfanned"]))
partlyraised.notfanned.ci
#partlyraised.fanned
partlyraised.fanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positiong.partlyraised.fanned"])) 
partlyraised.fanned.est
partlyraised.fanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positiong.partlyraised.fanned"]))
partlyraised.fanned.ci
#raised.notfanned
raised.notfanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positionh.raised.notfanned"])) 
raised.notfanned.est
raised.notfanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positionh.raised.notfanned"]))
raised.notfanned.ci
#side.fanned
side.fanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positioni.side.fanned"]))
side.fanned.est
side.fanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positioni.side.fanned"]))
side.fanned.ci
#side.notfanned
side.notfanned.est<-inv.logit(mean(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"]) + tail.position.model.MCMCglmm$Sol[,"tail.positionj.side.notfanned"]))
side.notfanned.est
side.notfanned.ci<-inv.logit(HPDinterval(mcmc(tail.position.model.MCMCglmm$Sol[,"(Intercept)"])+ tail.position.model.MCMCglmm$Sol[,"tail.positionj.side.notfanned"]))
side.notfanned.ci


##insert figure
l.tail.position <- list(
  a = list(tail.position= "aa.trailing.notfanned", posture.assumed = trailing.notfanned.est, lower = trailing.notfanned.ci[1] ,upper = trailing.notfanned.ci[2]) 
  , b = list(tail.position= "b.trailing.fanned", posture.assumed =trailing.fanned.est, lower = trailing.fanned.ci[1],upper =  trailing.fanned.ci[2])
  , c = list(tail.position= "c.down.fanned", posture.assumed = down.fanned.est, lower=down.fanned.ci[1], upper= down.fanned.ci[2])
  , d = list(tail.position= "d.down.notfanned", posture.assumed = down.notfanned.est, lower=down.notfanned.ci[1], upper= down.notfanned.ci[2])
  , e = list(tail.position= "e.raised.fanned", posture.assumed = raised.fanned.est, lower=raised.fanned.ci[1], upper= raised.fanned.ci[2])
  , f = list(tail.position= "f.partlyraised.notfanned", posture.assumed = partlyraised.notfanned.est, lower=partlyraised.notfanned.ci[1], upper= partlyraised.notfanned.ci[2])
  , g = list(tail.position= "g.partlyraised.fanned", posture.assumed = partlyraised.fanned.est, lower=partlyraised.fanned.ci[1], upper= partlyraised.fanned.ci[2])
  , h = list(tail.position= "h.raised.notfanned", posture.assumed = raised.notfanned.est, lower=raised.notfanned.ci[1], upper= raised.notfanned.ci[2])
  , i = list(tail.position= "i.side.fanned", posture.assumed = side.fanned.est, lower=side.fanned.ci[1], upper= side.fanned.ci[2])
  , j = list(tail.position= "j.side.notfanned", posture.assumed = side.notfanned.est, lower=side.notfanned.ci[1], upper= side.notfanned.ci[2])
  
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
tail.position.df <- ldply (l.tail.position, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=tail.position.long,aes(tail.position,posture.assumed))
p2<-p+geom_jitter(data=tail.position.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=tail.position,y=posture.assumed),data=tail.position.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=tail.position,ymin=lower,ymax=upper),width=0.1,data=tail.position.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)


####posture: feet position####
##make long dataset for MCMCglmm analysis
names(signal.dat)
feet.position.dat<-signal.dat[,c(1:5,18,55,56)]
levels(as.factor(feet.position.dat$feet.position))
feet.position.dat<-feet.position.dat[-c(feet.position.dat$feet.position=="not.visible"),] #remove interactions in which position was not visible
nrow(feet.position.dat) #336 interactions remain

#create new binary columns for each body orientation category
feet.position.dat$aa.on.substrate<-NA
feet.position.dat$b.extended<-NA
feet.position.dat$c.tucked.up<-NA
feet.position.dat$d.partly.extended<-NA
feet.position.dat$e.hanging<-NA


for(i in 1:nrow(feet.position.dat)){
  if(feet.position.dat$feet.position[i]=="on.substrate"){feet.position.dat$aa.on.substrate[i]<-1} else{feet.position.dat$aa.on.substrate[i]<-0}
  if(feet.position.dat$feet.position[i]=="extended"){feet.position.dat$b.extended[i]<-1} else{feet.position.dat$b.extended[i]<-0}
  if(feet.position.dat$feet.position[i]=="tucked.up"){feet.position.dat$c.tucked.up[i]<-1} else{feet.position.dat$c.tucked.up[i]<-0}
  if(feet.position.dat$feet.position[i]=="partly.extended"){feet.position.dat$d.partly.extended[i]<-1} else{feet.position.dat$d.partly.extended[i]<-0}
  if(feet.position.dat$feet.position[i]=="hanging"){feet.position.dat$e.hanging[i]<-1} else{feet.position.dat$e.hanging[i]<-0}
}

#makde dataset long for analysis
library(tidyr)
names(feet.position.dat)
feet.position.long<-gather(feet.position.dat,feet.position,posture.assumed,c(9:ncol(feet.position.dat)))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
feet.position.long$posture.assumed.bv[feet.position.long$posture.assumed > 0] <- "yes"
feet.position.long$posture.assumed.bv[feet.position.long$posture.assumed == 0] <- "no"
feet.position.long$posture.assumed.bv <- as.factor(feet.position.long$posture.assumed.bv)

k <- length(levels(feet.position.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
feet.position.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ feet.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = feet.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
feet.position.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ feet.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = feet.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
feet.position.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ feet.position, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = feet.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
feet.position.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = feet.position.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(feet.position.model.MCMCglmm$Sol,feet.position.model.MCMCglmm.2$Sol, feet.position.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(feet.position.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(feet.position.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(feet.position.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- feet.position.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(feet.position.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
feet.position.model.MCMCglmm$DIC
feet.position.model.MCMCglmm.null$DIC

summary(feet.position.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#on.substrate
on.substrate.est<-inv.logit(mean(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
on.substrate.est
on.substrate.ci<-inv.logit(HPDinterval(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"]))) 
on.substrate.ci
#extended
extended.est<-inv.logit(mean(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"]) + feet.position.model.MCMCglmm$Sol[,"feet.positionb.extended"])) 
extended.est
extended.ci<-inv.logit(HPDinterval(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"])+ feet.position.model.MCMCglmm$Sol[,"feet.positionb.extended"]))
extended.ci
#tucked.up
tucked.up.est<-tucked.up.est<-inv.logit(mean(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"]) + feet.position.model.MCMCglmm$Sol[,"feet.positionc.tucked.up"]))
tucked.up.est
tucked.up.ci<-inv.logit(HPDinterval(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"])+ feet.position.model.MCMCglmm$Sol[,"feet.positionc.tucked.up"]))
tucked.up.ci
#partly.extended
partly.extended.est<-inv.logit(mean(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"]) + feet.position.model.MCMCglmm$Sol[,"feet.positiond.partly.extended"]))
partly.extended.est
partly.extended.ci<-inv.logit(HPDinterval(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"])+ feet.position.model.MCMCglmm$Sol[,"feet.positiond.partly.extended"]))
partly.extended.ci
#hanging
hanging.est<-inv.logit(mean(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"]) + feet.position.model.MCMCglmm$Sol[,"feet.positione.hanging"]))
hanging.est
hanging.ci<-inv.logit(HPDinterval(mcmc(feet.position.model.MCMCglmm$Sol[,"(Intercept)"])+ feet.position.model.MCMCglmm$Sol[,"feet.positione.hanging"]))
hanging.ci

##create figure
l.feet.position <- list(
  a = list(feet.position= "aa.on.substrate", posture.assumed = on.substrate.est, lower = on.substrate.ci[1] ,upper = on.substrate.ci[2]) 
  , b = list(feet.position= "b.extended", posture.assumed =extended.est, lower = extended.ci[1],upper =  extended.ci[2])
  , c = list(feet.position= "c.tucked.up", posture.assumed = tucked.up.est, lower=tucked.up.ci[1], upper= tucked.up.ci[2])
  , d = list(feet.position= "d.partly.extended", posture.assumed = partly.extended.est, lower=partly.extended.ci[1], upper= partly.extended.ci[2])
  , e = list(feet.position= "e.hanging", posture.assumed = hanging.est, lower=hanging.ci[1], upper= hanging.ci[2])
 
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
feet.position.df <- ldply (l.feet.position, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=feet.position.long,aes(feet.position,posture.assumed))
p2<-p+geom_jitter(data=feet.position.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=feet.position,y=posture.assumed),data=feet.position.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=feet.position,ymin=lower,ymax=upper),width=0.1,data=feet.position.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)



####posture: closest.point####
##make long dataset for MCMCglmm analysis
names(signal.dat)
closest.point.dat<-signal.dat[,c(1:5,19,55,56)]
levels(as.factor(closest.point.dat$closest.point))
nrow(closest.point.dat) #337 interactions

#create new binary columns for each body orientation category
closest.point.dat$aa.bill<-NA
closest.point.dat$b.bill.wing<-NA
closest.point.dat$c.bill.feet<-NA
closest.point.dat$d.bill.neck<-NA
closest.point.dat$e.bill.forehead<-NA
closest.point.dat$f.bill.face<-NA
closest.point.dat$g.bill.head<-NA
closest.point.dat$h.bill.throat.neck<-NA
closest.point.dat$i.bill.tail<-NA
closest.point.dat$j.feet<-NA
closest.point.dat$k.wing<-NA
closest.point.dat$l.breast<-NA
closest.point.dat$m.feet.tail<-NA
closest.point.dat$n.side<-NA

for(i in 1:nrow(closest.point.dat)){
  if(closest.point.dat$closest.point[i]=="bill"){closest.point.dat$aa.bill[i]<-1} else{closest.point.dat$aa.bill[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.wing"){closest.point.dat$b.bill.wing[i]<-1} else{closest.point.dat$b.bill.wing[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.feet"){closest.point.dat$c.bill.feet[i]<-1} else{closest.point.dat$c.bill.feet[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.neck"){closest.point.dat$d.bill.neck[i]<-1} else{closest.point.dat$d.bill.neck[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.forehead"){closest.point.dat$e.bill.forehead[i]<-1} else{closest.point.dat$e.bill.forehead[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.face"){closest.point.dat$f.bill.face[i]<-1} else{closest.point.dat$f.bill.face[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.head"){closest.point.dat$g.bill.head[i]<-1} else{closest.point.dat$g.bill.head[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.throat.neck"){closest.point.dat$h.bill.throat.neck[i]<-1} else{closest.point.dat$h.bill.throat.neck[i]<-0}
  if(closest.point.dat$closest.point[i]=="bill.tail"){closest.point.dat$i.bill.tail[i]<-1} else{closest.point.dat$i.bill.tail[i]<-0}
  if(closest.point.dat$closest.point[i]=="feet"){closest.point.dat$j.feet[i]<-1} else{closest.point.dat$j.feet[i]<-0}
  if(closest.point.dat$closest.point[i]=="wing"){closest.point.dat$k.wing[i]<-1} else{closest.point.dat$k.wing[i]<-0}
  if(closest.point.dat$closest.point[i]=="breast"){closest.point.dat$l.breast[i]<-1} else{closest.point.dat$l.breast[i]<-0}
  if(closest.point.dat$closest.point[i]=="feet.tail"){closest.point.dat$m.feet.tail[i]<-1} else{closest.point.dat$m.feet.tail[i]<-0}
  if(closest.point.dat$closest.point[i]=="side"){closest.point.dat$n.side[i]<-1} else{closest.point.dat$n.side[i]<-0}
}

#makde dataset long for analysis
library(tidyr)
names(closest.point.dat)
closest.point.long<-gather(closest.point.dat,closest.point,posture.assumed,c(9:ncol(closest.point.dat)))

##analysis controlling for phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary response variable
closest.point.long$posture.assumed.bv[closest.point.long$posture.assumed > 0] <- "yes"
closest.point.long$posture.assumed.bv[closest.point.long$posture.assumed == 0] <- "no"
closest.point.long$posture.assumed.bv <- as.factor(closest.point.long$posture.assumed.bv)

k <- length(levels(closest.point.long$posture.assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
closest.point.model.MCMCglmm <- MCMCglmm(posture.assumed.bv ~ closest.point, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = closest.point.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
closest.point.model.MCMCglmm.2 <- MCMCglmm(posture.assumed.bv ~ closest.point, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = closest.point.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)
closest.point.model.MCMCglmm.3 <- MCMCglmm(posture.assumed.bv ~ closest.point, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = closest.point.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
closest.point.model.MCMCglmm.null <- MCMCglmm(posture.assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = closest.point.long, nitt=5000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(closest.point.model.MCMCglmm$Sol,closest.point.model.MCMCglmm.2$Sol, closest.point.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(closest.point.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(closest.point.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(closest.point.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- closest.point.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(closest.point.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
closest.point.model.MCMCglmm$DIC
closest.point.model.MCMCglmm.null$DIC

summary(closest.point.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#bill
bill.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]))) 
bill.est
bill.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]))) 
bill.ci
#bill.wing
bill.wing.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointb.bill.wing"])) 
bill.wing.est
bill.wing.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointb.bill.wing"]))
bill.wing.ci
#bill.feet
bill.feet.est<-bill.feet.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointc.bill.feet"]))
bill.feet.est
bill.feet.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointc.bill.feet"])) 
bill.feet.ci
#bill.neck
bill.neck.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointd.bill.neck"]))
bill.neck.est
bill.neck.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointd.bill.neck"])) 
bill.neck.ci
#bill.forehead
bill.forehead.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointe.bill.forehead"])) 
bill.forehead.est
bill.forehead.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointe.bill.forehead"]))
bill.forehead.ci
#bill.face
bill.face.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointf.bill.face"])) 
bill.face.est
bill.face.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointf.bill.face"])) 
bill.face.ci
#bill.head
bill.head.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointg.bill.head"])) 
bill.head.est
bill.head.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointg.bill.head"])) 
bill.head.ci
#bill.throat.neck
bill.throat.neck.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointh.bill.throat.neck"])) 
bill.throat.neck.est
bill.throat.neck.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointh.bill.throat.neck"])) 
bill.throat.neck.ci
#bill.tail
bill.tail.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointi.bill.tail"])) 
bill.tail.est
bill.tail.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointi.bill.tail"])) 
bill.tail.ci
#feet
feet.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointj.feet"])) 
feet.est
feet.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointj.feet"])) 
feet.ci
#wing
wing.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointk.wing"]))
wing.est
wing.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointk.wing"]))
wing.ci
#breast
breast.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointl.breast"]))
breast.est
breast.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointl.breast"]))
breast.ci
#feet.tail
feet.tail.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointm.feet.tail"])) 
feet.tail.est
feet.tail.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointm.feet.tail"])) 
feet.tail.ci
#side
side.est<-inv.logit(mean(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"]) + closest.point.model.MCMCglmm$Sol[,"closest.pointn.side"])) 
side.est
side.ci<-inv.logit(HPDinterval(mcmc(closest.point.model.MCMCglmm$Sol[,"(Intercept)"])+ closest.point.model.MCMCglmm$Sol[,"closest.pointn.side"])) 
side.ci


##insert figure
l.closest.point <- list(
  a = list(closest.point= "aa.bill", posture.assumed = bill.est, lower = bill.ci[1] ,upper = bill.ci[2]) 
  , b = list(closest.point= "b.bill.wing", posture.assumed =bill.wing.est, lower = bill.wing.ci[1],upper =  bill.wing.ci[2])
  , c = list(closest.point= "c.bill.feet", posture.assumed = bill.feet.est, lower=bill.feet.ci[1], upper= bill.feet.ci[2])
  , d = list(closest.point= "d.bill.neck", posture.assumed = bill.neck.est, lower=bill.neck.ci[1], upper= bill.neck.ci[2])
  , e = list(closest.point= "e.bill.forehead", posture.assumed = bill.forehead.est, lower=bill.forehead.ci[1], upper= bill.forehead.ci[2])
  , f = list(closest.point= "f.bill.face", posture.assumed = bill.face.est, lower = bill.face.ci[1] ,upper = bill.face.ci[2]) 
  , g = list(closest.point= "g.bill.head", posture.assumed =bill.head.est, lower = bill.head.ci[1],upper =  bill.head.ci[2])
  , h = list(closest.point= "h.bill.throat.neck", posture.assumed = bill.throat.neck.est, lower=bill.throat.neck.ci[1], upper= bill.throat.neck.ci[2])
  , i = list(closest.point= "i.bill.tail", posture.assumed = bill.tail.est, lower=bill.tail.ci[1], upper= bill.tail.ci[2])
  , j = list(closest.point= "j.feet", posture.assumed = feet.est, lower=feet.ci[1], upper= feet.ci[2])
  , k = list(closest.point= "k.wing", posture.assumed = wing.est, lower=wing.ci[1], upper= wing.ci[2])
  , l = list(closest.point= "l.breast", posture.assumed = breast.est, lower = breast.ci[1] ,upper = breast.ci[2]) 
  , m = list(closest.point= "m.feet.tail", posture.assumed =feet.tail.est, lower = feet.tail.ci[1],upper =  feet.tail.ci[2])
  , n = list(closest.point= "n.side", posture.assumed = side.est, lower=side.ci[1], upper= side.ci[2])
  
)  #list containing lists of estimates and confidence intervals of each predictor

library (plyr)
closest.point.df <- ldply (l.closest.point, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=closest.point.long,aes(closest.point,posture.assumed))
p2<-p+geom_jitter(data=closest.point.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=closest.point,y=posture.assumed),data=closest.point.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=closest.point,ymin=lower,ymax=upper),width=0.1,data=closest.point.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)




#######################################################################
####Body regions highlighted in aggressive interspecific encounters####
#######################################################################

names(signal.dat) # lists all variable names

consensus.regions<-signal.dat[,c(1:5,29:45,55,56)]
names(consensus.regions)
#find which regions are agreed upon as most highlighted by at least 2 (>1) observers
for(i in 1:nrow(consensus.regions)){
  if(consensus.regions$mouth[i]>1){consensus.regions$mouth[i]<-1}else{consensus.regions$mouth[i]<-0}
  if(consensus.regions$bill[i]>1){consensus.regions$bill[i]<-1}else{consensus.regions$bill[i]<-0}
  if(consensus.regions$face.throat[i]>1){consensus.regions$face.throat[i]<-1}else{consensus.regions$face.throat[i]<-0}
  if(consensus.regions$breast[i]>1){consensus.regions$breast[i]<-1}else{consensus.regions$breast[i]<-0}
  if(consensus.regions$belly[i]>1){consensus.regions$belly[i]<-1}else{consensus.regions$belly[i]<-0}
  if(consensus.regions$sides[i]>1){consensus.regions$sides[i]<-1}else{consensus.regions$sides[i]<-0}
  if(consensus.regions$crown[i]>1){consensus.regions$crown[i]<-1}else{consensus.regions$crown[i]<-0}
  if(consensus.regions$nape.back[i]>1){consensus.regions$nape.back[i]<-1}else{consensus.regions$nape.back[i]<-0}
  if(consensus.regions$shoulders[i]>1){consensus.regions$shoulders[i]<-1}else{consensus.regions$shoulders[i]<-0}
  if(consensus.regions$upperwings[i]>1){consensus.regions$upperwings[i]<-1}else{consensus.regions$upperwings[i]<-0}
  if(consensus.regions$underwings[i]>1){consensus.regions$underwings[i]<-1}else{consensus.regions$underwings[i]<-0}
  if(consensus.regions$uppertail[i]>1){consensus.regions$uppertail[i]<-1}else{consensus.regions$uppertail[i]<-0}
  if(consensus.regions$uppertail.coverts[i]>1){consensus.regions$uppertail.coverts[i]<-1}else{consensus.regions$uppertail.coverts[i]<-0}
  if(consensus.regions$undertail[i]>1){consensus.regions$undertail[i]<-1}else{consensus.regions$undertail[i]<-0}
  if(consensus.regions$undertail.coverts[i]>1){consensus.regions$undertail.coverts[i]<-1}else{consensus.regions$undertail.coverts[i]<-0}
  if(consensus.regions$tarsal.feathers[i]>1){consensus.regions$tarsal.feathers[i]<-1}else{consensus.regions$tarsal.feathers[i]<-0}
  if(consensus.regions$legs.feet[i]>1){consensus.regions$legs.feet[i]<-1}else{consensus.regions$legs.feet[i]<-0}
}


#make long dataset for analysis (omitting regions that were never highlighted: uppertail.coverts,undertain.coverts,tarsal.feathers)
library(tidyr)
names(consensus.regions)
regions.long.all<-gather(consensus.regions,region,highlighted,c(6:22))

library(dplyr)

for (i in 1:nrow(regions.long.all)){
  if (regions.long.all$region[i]=="face.throat") {regions.long.all$region[i]<-"aa.face.throat"}
  if (regions.long.all$region[i]=="mouth") {regions.long.all$region[i]<-"b.mouth"}
  if (regions.long.all$region[i]=="underwings") {regions.long.all$region[i]<-"c.underwings"}
  if (regions.long.all$region[i]=="bill") {regions.long.all$region[i]<-"d.bill"}
  if (regions.long.all$region[i]=="breast") {regions.long.all$region[i]<-"e.breast"}
  if (regions.long.all$region[i]=="nape.back") {regions.long.all$region[i]<-"f.nape.back"}
  if (regions.long.all$region[i]=="legs.feet") {regions.long.all$region[i]<-"g.legs.feet"}
  if (regions.long.all$region[i]=="crown") {regions.long.all$region[i]<-"h.crown"}
  if (regions.long.all$region[i]=="undertail") {regions.long.all$region[i]<-"i.undertail"}
  if (regions.long.all$region[i]=="upperwings") {regions.long.all$region[i]<-"j.upperwings"}
  if (regions.long.all$region[i]=="belly") {regions.long.all$region[i]<-"k.belly"}
  if (regions.long.all$region[i]=="uppertail") {regions.long.all$region[i]<-"l.uppertail"}
  if (regions.long.all$region[i]=="shoulders") {regions.long.all$region[i]<-"m.shoulders"}
  if (regions.long.all$region[i]=="sides") {regions.long.all$region[i]<-"n.sides"}
  if (regions.long.all$region[i]=="uppertail.coverts") {regions.long.all$region[i]<-"o.uppertail.coverts"}
  if (regions.long.all$region[i]=="undertail.coverts") {regions.long.all$region[i]<-"p.undertail.coverts"}
  if (regions.long.all$region[i]=="tarsal.feathers") {regions.long.all$region[i]<-"q.tarsal.feathers"}
}

regions.long<-subset(regions.long.all,region=="aa.face.throat"|region=="b.mouth"|region=="c.underwings"|region=="d.bill"|region=="e.breast"|region=="f.nape.back"|region=="g.legs.feet"|region=="h.crown"|region=="i.undertail"|region=="j.upperwings"|region=="k.belly"|region=="l.uppertail"|region=="m.shoulders"|region=="n.sides")

##phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary highlighted variable
regions.long$highlighted.bv[regions.long$highlighted > 0] <- "highlighted"
regions.long$highlighted.bv[regions.long$highlighted == 0] <- "aa.not.highlighted"
regions.long$highlighted.bv <- as.factor(regions.long$highlighted.bv)

k <- length(levels(regions.long$highlighted.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
regions.model.MCMCglmm <- MCMCglmm(highlighted.bv ~ region, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = regions.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
regions.model.MCMCglmm.2 <- MCMCglmm(highlighted.bv ~ region, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = regions.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
regions.model.MCMCglmm.3 <- MCMCglmm(highlighted.bv ~ region, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = regions.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
regions.model.MCMCglmm.null <- MCMCglmm(highlighted.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = regions.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(regions.model.MCMCglmm$Sol,regions.model.MCMCglmm.2$Sol, regions.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(regions.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(regions.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(regions.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- regions.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(regions.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
regions.model.MCMCglmm$DIC
regions.model.MCMCglmm.null$DIC

summary(regions.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#face.throat.bill
face.throat.bill.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]))) 
face.throat.bill.est
face.throat.bill.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])))
face.throat.bill.ci
#belly
belly.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionk.belly"]))
belly.est
belly.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionk.belly"]))
belly.ci
#bill
bill.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regiond.bill"]))
bill.est
bill.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regiond.bill"])) 
bill.ci
#breast
breast.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regione.breast"])) 
breast.est
breast.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regione.breast"])) 
breast.ci
#crown
crown.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionh.crown"]))
crown.est
crown.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionh.crown"])) 
crown.ci
#legs.feet
legs.feet.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regiong.legs.feet"])) 
legs.feet.est
legs.feet.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regiong.legs.feet"]))
legs.feet.ci
#mouth
mouth.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionb.mouth"])) 
mouth.est
mouth.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionb.mouth"])) 
mouth.ci
#nape.back
nape.back.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionf.nape.back"]))
nape.back.est
nape.back.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionf.nape.back"]))
nape.back.ci
#shoulders
shoulders.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionm.shoulders"])) 
shoulders.est
shoulders.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionm.shoulders"])) 
shoulders.ci
#sides
sides.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionn.sides"])) 
sides.est
sides.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionn.sides"]))
sides.ci
#undertail
undertail.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regioni.undertail"])) 
undertail.est
undertail.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regioni.undertail"]))
undertail.ci
#underwings
underwings.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionc.underwings"])) 
underwings.est
underwings.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionc.underwings"])) 
underwings.ci
#uppertail
uppertail.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionl.uppertail"])) 
uppertail.est
uppertail.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionl.uppertail"]))
uppertail.ci
#upperwings
upperwings.est<-inv.logit(mean(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"]) + regions.model.MCMCglmm$Sol[,"regionj.upperwings"])) 
upperwings.est
upperwings.ci<-inv.logit(HPDinterval(mcmc(regions.model.MCMCglmm$Sol[,"(Intercept)"])+ regions.model.MCMCglmm$Sol[,"regionj.upperwings"])) 
upperwings.ci


##insert figure
l.regions <- list(
  a = list(region= "aa.face.throat", highlighted = face.throat.bill.est, lower = face.throat.bill.ci[1] ,upper = face.throat.bill.ci[2]) 
  , b = list(region= "b.mouth", highlighted =mouth.est, lower = mouth.ci[1],upper =  mouth.ci[2])
  , c = list(region= "c.underwings", highlighted = underwings.est, lower=underwings.ci[1], upper= underwings.ci[2])
  , d = list(region= "d.bill", highlighted = bill.est, lower=bill.ci[1], upper= bill.ci[2])
  , e = list(region= "e.breast", highlighted = breast.est, lower=breast.ci[1], upper= breast.ci[2])
  , f = list(region= "f.nape.back", highlighted = nape.back.est, lower=nape.back.ci[1], upper= nape.back.ci[2])
  , g = list(region= "g.legs.feet", highlighted = legs.feet.est, lower=legs.feet.ci[1], upper= legs.feet.ci[2])
  , h = list(region= "h.crown", highlighted = crown.est, lower=crown.ci[1], upper= crown.ci[2])
  , i = list(region= "i.undertail", highlighted = undertail.est, lower=undertail.ci[1], upper= undertail.ci[2])
  , j = list(region= "j.upperwings", highlighted = upperwings.est, lower=upperwings.ci[1], upper= upperwings.ci[2])
  , k = list(region= "k.belly", highlighted = belly.est, lower=belly.ci[1], upper= belly.ci[2])
  , l = list(region= "l.uppertail", highlighted = uppertail.est, lower=uppertail.ci[1], upper= uppertail.ci[2])
  , m = list(region= "m.shoulders", highlighted = shoulders.est, lower=shoulders.ci[1], upper= shoulders.ci[2])
  , n = list(region= "n.sides", highlighted = sides.est, lower=sides.ci[1], upper= sides.ci[2])
  , o = list(region= "o.uppertail.coverts", highlighted = 0, lower = 0,upper = 0) 
  , p = list(region= "p.undertail.coverts", highlighted =0, lower = 0,upper = 0)
  , q = list(region= "q.tarsal.feathers", highlighted = 0, lower = 0,upper = 0)
)  #list containing lists of estimates and confidence intervals of each treatment

library (plyr)
regions.df <- ldply (l.regions, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=regions.long.all,aes(region,highlighted))
p2<-p+geom_jitter(data=regions.long.all,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=region,y=highlighted),data=regions.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=region,ymin=lower,ymax=upper),width=0.1,data=regions.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)


##print phylogeny figure with data
#for species represented in more than one interaction, this shows the most common result for each body region
library(ggtree)
library(tidyr)
library(dplyr)

fig.dat.reg<-consensus.regions
str(fig.dat.reg)

((unique(fig.dat.reg$animal))) #164
fig.dat.reg.bars<-data.frame(matrix(ncol=1,nrow=164))
fig.dat.reg.bars["animal"]<-(unique(fig.dat.reg$animal))
fig.dat.reg.bars["face.throat"]<-NA
fig.dat.reg.bars["mouth"]<-NA
fig.dat.reg.bars["underwings"]<-NA
fig.dat.reg.bars["bill"]<-NA
fig.dat.reg.bars["breast"]<-NA
fig.dat.reg.bars["nape.back"]<-NA
fig.dat.reg.bars["legs.feet"]<-NA
fig.dat.reg.bars["crown"]<-NA
fig.dat.reg.bars["undertail"]<-NA
fig.dat.reg.bars["upperwings"]<-NA
fig.dat.reg.bars["belly"]<-NA
fig.dat.reg.bars["uppertail"]<-NA
fig.dat.reg.bars["shoulders"]<-NA
fig.dat.reg.bars["sides"]<-NA
fig.dat.reg.bars["uppertail.coverts"]<-NA
fig.dat.reg.bars["undertail.coverts"]<-NA
fig.dat.reg.bars["tarsal.feathers"]<-NA

face.throat.ag<-aggregate(face.throat~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
mouth.ag<-aggregate(mouth~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
underwings.ag<-aggregate(underwings~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
bill.ag<-aggregate(bill~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
breast.ag<-aggregate(breast~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
nape.back.ag<-aggregate(nape.back~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
legs.feet.ag<-aggregate(legs.feet~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
crown.ag<-aggregate(crown~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
undertail.ag<-aggregate(undertail~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
upperwings.ag<-aggregate(upperwings~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
belly.ag<-aggregate(belly~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
uppertail.ag<-aggregate(uppertail~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
shoulders.ag<-aggregate(shoulders~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
sides.ag<-aggregate(sides~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
uppertail.coverts.ag<-aggregate(uppertail.coverts~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
undertail.coverts.ag<-aggregate(undertail.coverts~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)
tarsal.feathers.ag<-aggregate(tarsal.feathers~animal,data=fig.dat.reg,FUN=mean,na.action=na.omit)

for(i in 1:nrow(fig.dat.reg.bars)){
  
  tryCatch({
    if(face.throat.ag$face.throat[which(face.throat.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$face.throat[i]<-"highlighted"} else{fig.dat.reg.bars$face.throat[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(mouth.ag$mouth[which(mouth.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$mouth[i]<-"highlighted"} else{fig.dat.reg.bars$mouth[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(underwings.ag$underwings[which(underwings.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$underwings[i]<-"highlighted"} else{fig.dat.reg.bars$underwings[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(bill.ag$bill[which(bill.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$bill[i]<-"highlighted"} else{fig.dat.reg.bars$bill[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(breast.ag$breast[which(breast.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$breast[i]<-"highlighted"} else{fig.dat.reg.bars$breast[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(nape.back.ag$nape.back[which(nape.back.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$nape.back[i]<-"highlighted"} else{fig.dat.reg.bars$nape.back[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(legs.feet.ag$legs.feet[which(legs.feet.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$legs.feet[i]<-"highlighted"} else{fig.dat.reg.bars$legs.feet[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(crown.ag$crown[which(crown.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$crown[i]<-"highlighted"} else{fig.dat.reg.bars$crown[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(undertail.ag$undertail[which(undertail.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$undertail[i]<-"highlighted"} else{fig.dat.reg.bars$undertail[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(upperwings.ag$upperwings[which(upperwings.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$upperwings[i]<-"highlighted"} else{fig.dat.reg.bars$upperwings[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(belly.ag$belly[which(belly.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$belly[i]<-"highlighted"} else{fig.dat.reg.bars$belly[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(uppertail.ag$uppertail[which(uppertail.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$uppertail[i]<-"highlighted"} else{fig.dat.reg.bars$uppertail[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(shoulders.ag$shoulders[which(shoulders.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$shoulders[i]<-"highlighted"} else{fig.dat.reg.bars$shoulders[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(sides.ag$sides[which(sides.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$sides[i]<-"highlighted"} else{fig.dat.reg.bars$sides[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(uppertail.coverts.ag$uppertail.coverts[which(uppertail.coverts.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$uppertail.coverts[i]<-"highlighted"} else{fig.dat.reg.bars$uppertail.coverts[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(undertail.coverts.ag$undertail.coverts[which(undertail.coverts.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$undertail.coverts[i]<-"highlighted"} else{fig.dat.reg.bars$undertail.coverts[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(tarsal.feathers.ag$tarsal.feathers[which(tarsal.feathers.ag$animal==fig.dat.reg.bars$animal[i])]>=0.5){fig.dat.reg.bars$tarsal.feathers[i]<-"highlighted"} else{fig.dat.reg.bars$tarsal.feathers[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
}

str(fig.dat.reg.bars)
fig.dat.reg.bars<-fig.dat.reg.bars[,c(2:ncol(fig.dat.reg.bars))]
library(tibble)
fig.dat.reg.bars<-fig.dat.reg.bars%>%column_to_rownames('animal')
str(fig.dat.reg.bars)

p <- ggtree(tree) + 
  geom_tiplab(fontface='italic',size=2, align=TRUE, linesize=.5) + 
  theme_tree2(legend.position="none")
gheatmap(p, fig.dat.reg.bars, offset=50, width=1, 
         colnames=FALSE, legend_title="highlighted") +
  scale_x_ggtree() + 
  scale_fill_manual(breaks=c("highlighted", "not.highlighted"), 
                    values=c("black", "lightgray"), name="highlighted")





#################################################################
####Colors highlighted in aggressive interspecific encounters####
#################################################################

names(signal.dat) # lists all variable names

consensus.colours<-signal.dat[,c(1:5,46:53,55,56)]
names(consensus.colours)

#make long dataset for analysis
library(tidyr)
names(consensus.colours)
colours.long.all<-gather(consensus.colours,colours,highlighted.all,c(red.pink.orange.yellow:gray))

#omit NAs from dataset - examine only colours and colour groups that are present in each species
colours.long<-colours.long.all[which(is.na(colours.long.all$highlighted.all)==FALSE),]
nrow(colours.long)

#find which regions are agreed upon as most highlighted by at least 2 observers (>1)
for(i in 1:nrow(colours.long)){
  if(colours.long$highlighted.all[i]>1){colours.long$highlighted[i]<-1}else{colours.long$highlighted[i]<-0}
}  

library(dplyr)

for (i in 1:nrow(colours.long)){
  if(colours.long$colours[i]=="red.pink.orange.yellow"){colours.long$colours[i]<-"aa.red.pink.orange.yellow"}
  if (colours.long$colours[i]=="blue.green.violet") {colours.long$colours[i]<-"b.blue.green.violet"}
  if (colours.long$colours[i]=="black") {colours.long$colours[i]<-"c.black"}
  if (colours.long$colours[i]=="dark.white.contrast") {colours.long$colours[i]<-"d.dark.white.contrast"}
  if (colours.long$colours[i]=="brown.beige") {colours.long$colours[i]<-"e.brown.beige"}
  if (colours.long$colours[i]=="gray") {colours.long$colours[i]<-"f.gray"}
  if (colours.long$colours[i]=="white") {colours.long$colours[i]<-"g.white"}
  if (colours.long$colours[i]=="rufous.chestnut") {colours.long$colours[i]<-"h.rufous.chestnut"}
}


##phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary highlighted variable
colours.long$highlighted.bv[colours.long$highlighted > 0] <- "highlighted"
colours.long$highlighted.bv[colours.long$highlighted == 0] <- "aa.not.highlighted"
colours.long$highlighted.bv <- as.factor(colours.long$highlighted.bv)

k <- length(levels(colours.long$highlighted.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
colours.model.MCMCglmm <- MCMCglmm(highlighted.bv ~ colours, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = colours.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
colours.model.MCMCglmm.2 <- MCMCglmm(highlighted.bv ~ colours, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = colours.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
colours.model.MCMCglmm.3 <- MCMCglmm(highlighted.bv ~ colours, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = colours.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
colours.model.MCMCglmm.null <- MCMCglmm(highlighted.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = colours.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)


gelman.diag(list(colours.model.MCMCglmm$Sol,colours.model.MCMCglmm.2$Sol, colours.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(colours.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(colours.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(colours.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- colours.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(colours.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
colours.model.MCMCglmm$DIC
colours.model.MCMCglmm.null$DIC

summary(colours.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#red.pink.orange.yellow
red.pink.orange.yellow.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]))) 
red.pink.orange.yellow.est
red.pink.orange.yellow.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])))
red.pink.orange.yellow.ci
#blue.green.violet
blue.green.violet.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]) + colours.model.MCMCglmm$Sol[,"coloursb.blue.green.violet"]))
blue.green.violet.est
blue.green.violet.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])+ colours.model.MCMCglmm$Sol[,"coloursb.blue.green.violet"]))
blue.green.violet.ci
#black
black.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]) + colours.model.MCMCglmm$Sol[,"coloursc.black"]))
black.est
black.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])+ colours.model.MCMCglmm$Sol[,"coloursc.black"])) 
black.ci
#dark.white.contrast
dark.white.contrast.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]) + colours.model.MCMCglmm$Sol[,"coloursd.dark.white.contrast"])) 
dark.white.contrast.est
dark.white.contrast.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])+ colours.model.MCMCglmm$Sol[,"coloursd.dark.white.contrast"])) 
dark.white.contrast.ci
#brown.beige 
brown.beige.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]) + colours.model.MCMCglmm$Sol[,"colourse.brown.beige"]))
brown.beige.est
brown.beige.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])+ colours.model.MCMCglmm$Sol[,"colourse.brown.beige"])) 
brown.beige.ci
#gray
gray.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]) + colours.model.MCMCglmm$Sol[,"coloursf.gray"])) 
gray.est
gray.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])+ colours.model.MCMCglmm$Sol[,"coloursf.gray"]))
gray.ci
#white
white.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]) + colours.model.MCMCglmm$Sol[,"coloursg.white"])) 
white.est
white.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])+ colours.model.MCMCglmm$Sol[,"coloursg.white"])) 
white.ci
#rufous.chestnut
rufous.chestnut.est<-inv.logit(mean(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"]) + colours.model.MCMCglmm$Sol[,"coloursh.rufous.chestnut"]))
rufous.chestnut.est
rufous.chestnut.ci<-inv.logit(HPDinterval(mcmc(colours.model.MCMCglmm$Sol[,"(Intercept)"])+ colours.model.MCMCglmm$Sol[,"coloursh.rufous.chestnut"]))
rufous.chestnut.ci

##insert figure
l.colours <- list(
  a = list(colours= "aa.red.pink.orange.yellow", highlighted = red.pink.orange.yellow.est, lower = red.pink.orange.yellow.ci[1] ,upper = red.pink.orange.yellow.ci[2]) 
  , b = list(colours= "b.blue.green.violet", highlighted =blue.green.violet.est, lower = blue.green.violet.ci[1],upper =  blue.green.violet.ci[2])
  , c = list(colours= "c.black", highlighted = black.est, lower=black.ci[1], upper= black.ci[2])
  , d = list(colours= "d.dark.white.contrast", highlighted = dark.white.contrast.est, lower=dark.white.contrast.ci[1], upper= dark.white.contrast.ci[2])
  , e = list(colours= "e.brown.beige", highlighted = brown.beige.est, lower=brown.beige.ci[1], upper= brown.beige.ci[2])
  , f = list(colours= "f.gray", highlighted = gray.est, lower=gray.ci[1], upper= gray.ci[2])
  , g = list(colours= "g.white", highlighted = white.est, lower=white.ci[1], upper= white.ci[2])
  , h = list(colours= "h.rufous.chestnut", highlighted = rufous.chestnut.est, lower=rufous.chestnut.ci[1], upper= rufous.chestnut.ci[2])
)  #list containing lists of estimates and confidence intervals of each treatment

library (plyr)
colours.df <- ldply (l.colours, data.frame) #create dataframe with estimates
detach(package:plyr)

library(ggplot2)
##figure
p<-ggplot(data=colours.long,aes(colours,highlighted))
p2<-p+geom_jitter(data=colours.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=colours,y=highlighted),data=colours.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=colours,ymin=lower,ymax=upper),width=0.1,data=colours.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)


##print phylogeny figure with data
#for species represented in more than one interaction, this shows the most common result for each colour
library(ggtree)
library(ape)
library(tidyr)
library(dplyr)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre

names(consensus.colours)
fig.dat.col<-consensus.colours
str(fig.dat.col)

((unique(fig.dat.col$animal))) #164
fig.dat.col.bars<-data.frame(matrix(ncol=1,nrow=164))
fig.dat.col.bars["animal"]<-(unique(fig.dat.col$animal))
fig.dat.col.bars["red.pink.orange.yellow"]<-NA
fig.dat.col.bars["blue.green.violet"]<-NA
fig.dat.col.bars["dark.white.contrast"]<-NA
fig.dat.col.bars["black"]<-NA
fig.dat.col.bars["brown.beige"]<-NA
fig.dat.col.bars["gray"]<-NA
fig.dat.col.bars["white"]<-NA
fig.dat.col.bars["rufous.chestnut"]<-NA

red.pink.orange.yellow.ag<-aggregate(red.pink.orange.yellow~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)
blue.green.violet.ag<-aggregate(blue.green.violet~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)
dark.white.contrast.ag<-aggregate(dark.white.contrast~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)
black.ag<-aggregate(black~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)
brown.beige.ag<-aggregate(brown.beige~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)
gray.ag<-aggregate(gray~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)
white.ag<-aggregate(white~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)
rufous.chestnut.ag<-aggregate(rufous.chestnut~animal,data=fig.dat.col,FUN=mean,na.action=na.omit)

for(i in 1:nrow(fig.dat.col.bars)){
  
  tryCatch({
    if(red.pink.orange.yellow.ag$red.pink.orange.yellow[which(red.pink.orange.yellow.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$red.pink.orange.yellow[i]<-"highlighted"} else{fig.dat.col.bars$red.pink.orange.yellow[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(blue.green.violet.ag$blue.green.violet[which(blue.green.violet.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$blue.green.violet[i]<-"highlighted"} else{fig.dat.col.bars$blue.green.violet[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(dark.white.contrast.ag$dark.white.contrast[which(dark.white.contrast.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$dark.white.contrast[i]<-"highlighted"} else{fig.dat.col.bars$dark.white.contrast[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(black.ag$black[which(black.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$black[i]<-"highlighted"} else{fig.dat.col.bars$black[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(brown.beige.ag$brown.beige[which(brown.beige.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$brown.beige[i]<-"highlighted"} else{fig.dat.col.bars$brown.beige[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(gray.ag$gray[which(gray.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$gray[i]<-"highlighted"} else{fig.dat.col.bars$gray[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(white.ag$white[which(white.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$white[i]<-"highlighted"} else{fig.dat.col.bars$white[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
  tryCatch({
    if(rufous.chestnut.ag$rufous.chestnut[which(rufous.chestnut.ag$animal==fig.dat.col.bars$animal[i])]>=0.5){fig.dat.col.bars$rufous.chestnut[i]<-"highlighted"} else{fig.dat.col.bars$rufous.chestnut[i]<-"not.highlighted"}
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
}

str(fig.dat.col.bars)
fig.dat.col.bars<-fig.dat.col.bars[,c(2:ncol(fig.dat.col.bars))]
library(tibble)
fig.dat.col.bars<-fig.dat.col.bars%>%column_to_rownames('animal')
str(fig.dat.col.bars)

p <- ggtree(tree) + 
  geom_tiplab(fontface='italic',size=2, align=TRUE, linesize=.5) + 
  theme_tree2(legend.position="none")
gheatmap(p, fig.dat.col.bars, offset=30, width=0.5, 
         colnames=FALSE, legend_title="color highlighted") +
  scale_x_ggtree() + 
  scale_fill_manual(breaks=c("highlighted", "not.highlighted"), 
                    values=c("red", "lightgray"), name="highlighted")




##############################################################################
####Similarlity to signals used in aggressive encounters with conspecifics####
##############################################################################

names(signal.dat)

#we have scored similarity for 307 interactions
#253 are "same"

similarity.dat<-signal.dat[,c(1:5,20,55,56)]
names(similarity.dat)

#remove unknown species from dataframe
similarity.dat.2<-subset(similarity.dat,relation.to.intra.signal=="same"|relation.to.intra.signal=="similar"|relation.to.intra.signal=="similar to congener"|relation.to.intra.signal=="different")
nrow(similarity.dat.2) #307 interactions

for(i in 1:nrow(similarity.dat.2)){
  if(similarity.dat.2$relation.to.intra.signal[i]=="same"){similarity.dat.2$a.same[i]<-1}else{similarity.dat.2$a.same[i]<-0}
  if(similarity.dat.2$relation.to.intra.signal[i]=="similar"){similarity.dat.2$b.similar[i]<-1}else{similarity.dat.2$b.similar[i]<-0}
  if(similarity.dat.2$relation.to.intra.signal[i]=="similar to congener"){similarity.dat.2$c.similar.to.congener[i]<-1}else{similarity.dat.2$c.similar.to.congener[i]<-0}
  if(similarity.dat.2$relation.to.intra.signal[i]=="different"){similarity.dat.2$d.different[i]<-1}else{similarity.dat.2$d.different[i]<-0}
}


#make long dataset for analysis
library(tidyr)
names(similarity.dat.2)
similarity.long<-gather(similarity.dat.2,similarity,assumed,c(9:12))

##phylogeny
library(ape)
tree<-read.tree(file.choose()) #load phylogeny: tree_1.tre
is.ultrametric(tree)

library(MCMCglmm)

#create categorical variable for binary highlighted variable
similarity.long$assumed.bv[similarity.long$assumed > 0] <- "yes"
similarity.long$assumed.bv[similarity.long$assumed == 0] <- "no"
similarity.long$assumed.bv <- as.factor(similarity.long$assumed.bv)

k <- length(levels(similarity.long$assumed.bv))
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))
prior.model1 <- list(R=list(fix=1, V=(1/k)*(I+J), nu=k-1), G=list(G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1),G1=list(V=diag(k-1), nu=k-1)))
similarity.model.MCMCglmm <- MCMCglmm(assumed.bv ~ similarity, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = similarity.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

# model checks (see details and references in Supplemental Materials and Methods)
similarity.model.MCMCglmm.2 <- MCMCglmm(assumed.bv ~ similarity, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = similarity.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)
similarity.model.MCMCglmm.3 <- MCMCglmm(assumed.bv ~ similarity, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = similarity.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

#null model for comparison
similarity.model.MCMCglmm.null <- MCMCglmm(assumed.bv ~ 1, random =~ animal + focal.interaction + signaling.IOC, pedigree = tree, family="categorical", data = similarity.long, nitt=2000000, burnin=10000, thin=100, prior=prior.model1)

gelman.diag(list(similarity.model.MCMCglmm$Sol,similarity.model.MCMCglmm.2$Sol, similarity.model.MCMCglmm.3$Sol)) # check for model convergence across the 3 runs; report results from first model; upper limits of potential scale reduction factors (psrf) should approach 1
library(lattice)
print(xyplot(similarity.model.MCMCglmm$Sol)) # trace plots should show no trends and look like white noise, with rapid and featureless variation
round(sort(effectiveSize(similarity.model.MCMCglmm$Sol))) # effective sample sizes should be >200
geweke.diag(similarity.model.MCMCglmm$Sol) # values should fall within the 95% confidence intervals of the standard normal (i.e., |x| < 1.96)
vv <- similarity.model.MCMCglmm$VCV
print(xyplot(vv)) # trace plots should show no trends
effectiveSize(vv) # effective sample sizes should be >200
print(densityplot(similarity.model.MCMCglmm$Sol)) # density plots of the posterior distributions - should be generally symmetrical and normal
print(densityplot(vv))

#model comparison
similarity.model.MCMCglmm$DIC
similarity.model.MCMCglmm.null$DIC

summary(similarity.model.MCMCglmm) 

##back-transform model estimates and CIs
library(boot)
#same
same.est<-inv.logit(mean(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"]))) 
same.est
same.ci<-inv.logit(HPDinterval(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"])))
same.ci
#similar
similar.est<-inv.logit(mean(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"]) + similarity.model.MCMCglmm$Sol[,"similarityb.similar"]))
similar.est
similar.ci<-inv.logit(HPDinterval(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"])+ similarity.model.MCMCglmm$Sol[,"similarityb.similar"]))
similar.ci
#similar.to.congener
similar.to.congener.est<-inv.logit(mean(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"]) + similarity.model.MCMCglmm$Sol[,"similarityc.similar.to.congener"]))
similar.to.congener.est
similar.to.congener.ci<-inv.logit(HPDinterval(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"])+ similarity.model.MCMCglmm$Sol[,"similarityc.similar.to.congener"])) 
similar.to.congener.ci
#different
different.est<-inv.logit(mean(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"]) + similarity.model.MCMCglmm$Sol[,"similarityd.different"])) 
different.est
different.ci<-inv.logit(HPDinterval(mcmc(similarity.model.MCMCglmm$Sol[,"(Intercept)"])+ similarity.model.MCMCglmm$Sol[,"similarityd.different"])) 
different.ci


##create figure
l.similarity <- list(
  a = list(similarity= "a.same", assumed = same.est, lower = same.ci[1] ,upper = same.ci[2]) 
  , b = list(similarity= "b.similar", assumed =similar.est, lower = similar.ci[1],upper =  similar.ci[2])
  , c = list(similarity= "c.similar.to.congener", assumed = similar.to.congener.est, lower=similar.to.congener.ci[1], upper= similar.to.congener.ci[2])
  , d = list(similarity= "d.different", assumed = different.est, lower=different.ci[1], upper= different.ci[2])
)  #list containing lists of estimates and confidence intervals of each treatment

library (plyr)
similar.df <- ldply (l.similarity, data.frame) #create dataframe with estimates
detach(package:plyr)


library(ggplot2)
##figure
p<-ggplot(data=similarity.long,aes(similarity,assumed))
p2<-p+geom_jitter(data=similarity.long,width=0.25,height=0.025,color="gray",pch=1)
p3<-p2+geom_point(aes(x=similarity,y=assumed),data=similar.df,size=3,inherit.aes=FALSE)
p4<-p3+geom_errorbar(aes(x=similarity,ymin=lower,ymax=upper),width=0.1,data=similar.df)
p6<-p4+scale_x_discrete(guide = guide_axis(angle=45)) +theme_light() +theme(axis.text.x=element_text(size=rel(1.5)),axis.text.y=element_text(size=rel(1.5)))
print(p6)