library(phytools)
library(paleotree)

###Extant birds

#Data preparation
data.ind=read.csv("raw data.csv")
bm= tapply(data.ind$bodymass, data.ind$species, mean)
kl= tapply(data.ind$keel, data.ind$species, mean)
il= tapply(data.ind$ilium, data.ind$species, mean)

lg.bm=log10(bm)
lg.kl=log10(kl)
lg.il=log10(il)

#Phylogeny

trees=read.nexus("24721.tre")

#phylogenetic size-correction

phyres=vector("list",length=1000)

for (i in 1:1000){
	phyres[[i]]<-phyl.resid(trees[[i]],x=lg.bm,Y=cbind(lg.kl,lg.il),method="lambda")
}

#Estimating evolutionary rate matrices

RR=vector("list",length=1000)

for (i in 1:1000){
  RR[[i]]<-evol.vcv(trees[[i]],phyres[[i]]$resid)
}

#

params=c()
for(i in 1:1000){
params=rbind(params,
data.frame(row.names=NULL,tree=i,k1=RR[[i]]$k1,logL1=RR[[i]]$logL1,r=cov2cor(RR[[i]]$R.single)[1,2]))
}
        
params$AICc=-2*params$logL1+2*params$k1+2*params$k1*(params$k1+1)/(137-params$k1-1)
params$delta=params$AICc-min(params$AICc)
params$w=exp(-0.5*params$delta)/sum(exp(-0.5*params$delta))

#model averaged correlation cofficient, 95% confidence interval of correlation. Fischer Transformation.

CIcorr<-function(corr,w,n){
  z<-qnorm(0.975)
  Fr<-atanh(corr)
  Fr.ave<-weighted.mean(Fr,w)
  err<-sqrt(1/(n-3)+(Fr-Fr.ave)^2)
  se.ave<-weighted.mean(err,w)
  corr.ave<-tanh(Fr.ave)
  LL<-tanh(Fr.ave-z*se.ave)
  UL<-tanh(Fr.ave+z*se.ave)
  x<-c(corr.ave,LL,UL)
  return(x)
}

r.ave=CIcorr(params$r,params$w,137)

#Plot Figure 2

order=read.csv("order.csv")
m=which(params$AICc==min(params$AICc))
dataset=data.frame(species=row.names(phyres[[m]]$resid),phyres[[m]]$resid)
dataset=merge(dataset,order,by="species")

library(ggplot2)
library(RColorBrewer)

p<-ggplot(data=dataset,aes(x=lg.kl,y=lg.il,colour=Order,shape=Order,label=species))+geom_point(size=3)+
  scale_colour_manual(values=brewer.pal(12,"Paired")[c(2,4,12,8,6,10,2,4,12,8,6,10,2,4,12,8,6,10)])+scale_shape_manual(values=c(19,19,19,19,19,19,17,17,17,17,17,17,18,18,18,18,18,18))+
  xlab("lg(keel length) (size-corrected)")+ylab("lg(ilium length) (size-corrected)")

#show names
q<-ggplot(data=dataset,aes(x=lg.kl,y=lg.il,colour=Order,shape=Order,label=species))+geom_point(size=3)+
  scale_colour_manual(values=brewer.pal(12,"Paired")[c(2,4,12,8,6,10,2,4,12,8,6,10,2,4,12,8,6,10)])+scale_shape_manual(values=c(19,19,19,19,19,19,17,17,17,17,17,17,18,18,18,18,18,18))+
  geom_text(aes(label=species))+xlab("lg(keel length) (size-corrected)")+ylab("lg(ilium length) (size-corrected)")

###Mesozoic birds

tree.eb=read.tree("early birds.tre")
eb=read.csv("early birds.csv")
fl.eb=eb$fl=log10(eb$fl)
kl.eb=eb$kl=log10(eb$kl)
il.eb=eb$il=log10(eb$il)
names(fl.eb)=eb$species
names(kl.eb)=eb$species
names(il.eb)=eb$species

ages=eb[,6:7]
row.names(ages)=eb$species

eb.trees=timePaleoPhy(tree.eb,ages,type="equal",vartime=2,ntrees=1000,dateTreatment="minMax")

eb.data.sc=vector("list",length=1000)
eb.ev=vector("list",length=1000)

for(i in 1:1000){
  eb.data.sc[[i]]=phyl.resid(eb.trees[[i]],x=fl.eb,Y=cbind(kl.eb,il.eb),method="BM")
  eb.ev[[i]]=evol.vcv(eb.trees[[i]],eb.data.sc[[i]]$resid)
  
}

eb.params=c()
for(i in 1:1000){
  eb.params=rbind(eb.params,data.frame(tree=i,logL1=eb.ev[[i]]$logL1,k1=eb.ev[[i]]$k1,r=cov2cor(eb.ev[[i]]$R.single)[1,2]))
}
eb.params$AICc=-2*eb.params$logL1+2*eb.params$k1+2*eb.params$k1*(eb.params$k1+1)/(10-eb.params$k1-1)
eb.params$delta=eb.params$AICc-min(eb.params$AICc)
eb.params$w=exp(-0.5*eb.params$delta)/sum(exp(-0.5*eb.params$delta))

r.ave.eb=CIcorr(eb.params$r,eb.params$w,10)


#Plot Figure 3

n=which(eb.params$AICc==min(eb.params$AICc))
phylomorphospace(eb.trees[[n]],eb.data.sc[[n]]$resid,label="horizontal",node.by.map=TRUE,
                 xlab="lg(keel length) (size corrected)",ylab="lg(ilium length) (size-corrected)")

###DONE