Compare evolutionary rates
fevicallmass<-read.csv("/home/m/Dropbox/Projects-Manuscripts/Manuscripts/Parrot Calls-Phylo paper/Results/call mass and fevi parrot phylo 9.9.14.csv",row.names=1,header=T)
# str(fevicallmass)
callparam2<-fevicallmass[,c(1,4:9,13:15)]
Arini51sp<-read.tree("/home/m/Dropbox/Projects-Manuscripts/Manuscripts/Parrot Calls-Phylo paper/Trees/Arini51spsep14.tre")
library(phytools)
#######################
#modified by DCA
library(ape)
###DCA: Read csv this way, and is already a matrix
#callparam<-read.csv("Call parameters parrots.csv")
# callparammat<-read.csv("Call parameters parrots.csv",header=TRUE,row.names=1)
#make it a matrix with species names in rows
#callparammat<-as.matrix(callparam[,2:9])
#rownames(callparammat)<-callparam$X
#scale due to different units
#scaledcallparam<-scale(callparammat)
#Note: In my paper, I recommend log transformation for scale issues (see also Felsenstein)
scaledcallparam<-log(callparam2)
Arini51sp<-read.tree("/home/m/Dropbox/Projects-Manuscripts/Manuscripts/Parrot Calls-Phylo paper/Trees/Arini51spsep14.tre")
#read tree
# Arini51sp<-read.tree("Arini51spsep14.tre")
# plot(Arini51sp)
#check match among data and tree
library(geiger)
name.check(Arini51sp, scaledcallparam)
## [1] "OK"
#compare evo rates
source('/home/m/Dropbox/Documentos R/CompareRates.multTrait function.R')
cr<-CompareRates.multTrait(Arini51sp,scaledcallparam,TraitCov=F)
colnames(cr$Robs)<-rownames(cr$Robs)<-colnames(scaledcallparam)
findCI<-function(phy,x){
library(numDeriv)
x<-as.matrix(x)
N<-nrow(x)
p<-ncol(x)
C<-vcv.phylo(phy)
C<-C[rownames(x),rownames(x)]
a.obs<-colSums(solve(C))%*%x/sum(solve(C))
D<-matrix(0,N*p,p)
for(i in 1:(N*p)) for(j in 1:p) if((j-1)*N<i&&i<=j*N) D[i,j]=1.0
y<-as.matrix(as.vector(x))
one<-matrix(1,N,1)
R.obs<-t(x-one%*%a.obs)%*%solve(C)%*%(x-one%*%a.obs)/N
lik<-function(R.obs,y,D,a.obs,C,N,p){
logL1<--t(y-D%*%t(a.obs))%*%solve(kronecker(R.obs,C))%*%(y-D%*%t(a.obs))/2-N*p*log(2*pi)/2-determinant(kronecker(R.obs,C))$modulus[1]/2
}
H<-hessian(lik,R.obs,y=y,D=D,a.obs=a.obs,C=C,N=N,p=p) #Hessian of likelihood fxn with parameters
std.err.all<-sqrt(diag(solve(-H))) #standard errors for parameters (NOTE: from full R matrix)
std.err.all[is.na(std.err.all)]<-0.0
std.err<-diag(matrix(std.err.all,nrow=p))/sqrt(N) #divide sqrt(N) for std. err
R.val<-diag(R.obs)
error <- qnorm(0.975)*std.err #Compute 95% CI and compare across traits
CI.l <- R.val-error; CI.u <- R.val+error #lower and upper CI
return(list(R=R.obs,CI.Lo=CI.l,CI.Hi=CI.u))
}
a<-findCI(Arini51sp,scaledcallparam[,1:5])
b<-findCI(Arini51sp,scaledcallparam[,6:10])
# c<-findCI(Arini51sp,scaledcallparam[,1:10])
sigma2<-c(diag(a$R),diag(b$R))
upsig<-c(a$CI.Lo,b$CI.Lo)
losig<-c(a$CI.Hi,b$CI.Hi)
sigma<-data.frame(names(sigma2),sigma2,upsig,losig)
combs<-combn(c(1:ncol(scaledcallparam)),2)
ahi<-list()
vars<-numeric()
Robs1<-numeric()
Robs2<-numeric()
P<-numeric()
AICobs<-numeric()
AICcon<-numeric()
messa<-numeric()
mat<-data.frame(matrix(nrow=ncol(scaledcallparam),ncol=ncol(scaledcallparam)))
for(i in 1:ncol(combs))
{ ahi<-CompareRates.multTrait2(Arini51sp,scaledcallparam[combs[,i]],TraitCov=T)
vars[i]<-paste(names(diag(ahi$Robs)),collapse ="-")
Robs1[i]<-diag(ahi$Robs)[1]
Robs2[i]<-diag(ahi$Robs)[2]
P[i]<-ahi$Prob
AICobs[i]<-ahi$AICc.obs
AICcon[i]<-ahi$AICc.constrained
messa[i]<-ahi$optimmessage
mat[combs[1,i],combs[2,i]]<-ahi$Prob
mat[combs[2,i],combs[1,i]]<-ahi$Prob
}
mat<-as.matrix(mat)
diag(mat)<-rep(0,ncol(mat))
mat<-as.data.frame(mat)
mat[mat>.05]<-NA
colnames(mat)<-rownames(mat)<-colnames(scaledcallparam)
# mat
# results<-data.frame(vars,Robs1,Robs2,P,AICobs,AICcon,messa)
# results[results$P<.05,]
rates<-data.frame(colnames(scaledcallparam)[which(colnames(scaledcallparam) %in% rownames(cr$Robs))],diag(cr$Robs))
colnames(rates)[1]<-"vars"
#
# #calcuate Confidence intervals for sigma2 (evolutionary rate)
# CIsigma2<-as.data.frame(matrix(nrow=ncol(scaledcallparam),ncol=3))
#
# for(i in 1:ncol(scaledcallparam))
# {aa<-ace(scaledcallparam[,i],Arini51sp)
# CIsigma2[i,2:3]<-aa$sigma2
# CIsigma2[i,1]<-colnames(scaledcallparam)[i]}
rates<-merge(rates,sigma,by.x="vars",by.y="names.sigma2.")
rates<-rates[,c(1,3:5)]
rates$vars=as.character(rates$vars)
rates$vars[grep("evi",rates$vars)]<-"Evi"
rates<-merge(rates,physigres[,c(1,5)],by.x="vars",by.y="Parameter",all=F)
colnames(rates)[3:5]<-c("Uplim","Lowlim","phylosig")
rates$phylosig[rates$phylosig<0.05]<-"yes"
rates$phylosig[rates$phylosig!="yes"]<-"no"
# rates[order(-rates[,2]),]
# rates<-rates[c(1,5,4,2,3,6,7),]
rates<-rates[c(9,6,3,7,5,1,10,8,2,4),]
rownames(rates)<-1:10
library(RColorBrewer)
# dev.off()
plot(rownames(rates),rates[,2],xaxt ="n",pch=20,cex=2,col="white",xlab="Acoustic parameters or predictors",ylab="Evolutionary rate",ylim=c(0,0.052))
for(i in 1:nrow(rates)){
lines(rep(i,2),rates[i,3:4])
lines(c(i+.1,i-.1),rep(rates[i,3],2))
lines(c(i+.1,i-.1),rep(rates[i,4],2))
points(1:nrow(rates),rates[,2],pch=20,cex=1)
if(rates$phylosig[i]=="yes") points(i,rates[i,2],pch=20,cex=0.5)}
axis(1,at=c(1:nrow(rates)),adj=1, tick=T,labels=rates[,1])
# rates[order(rates[,2]),]
# mat[order(rates[,2]),order(rates[,2])]
ley<-c("a","b","bc","c","c","cd","d","d","de","e")
ley<-data.frame(rates$vars[order(rates[,2])],ley)
rates<-merge(rates,ley,by.x="vars",by.y="rates.vars.order.rates...2...",sort = F)
# rates2<-data.frame(rates[order(rates[,2]),],ley)
# rates2<-rates2[order(as.numeric(rownames(rates2))),]
# rates2<-rates2[c(1,5,4,2,3,6:10),]
# str(rates)
text(1:nrow(rates),rates[,3]+0.006,labels=rates[,6])
points(1.02,0.059,pch=20)
text(1,0.059,"=phylogenetic signal",cex=0.8,pos=4)
abline(v=c(3.5,7.5))
text(c(2,5.5,9),rep(0.04,3),labels=c("Correlated evolution to morphology","Not correlated","Eco/morpho pedictors"))
# leye<-data.frame(rates[order(rates[,2]),1]
# ,character(nrow(rates)))
# CompareRates.multTrait(Arini51sp,scaledcallparam[,c(1,3)],TraitCov=T)
# CompareRates.multTrait(Arini51sp,scaledcallparam[,3:4],TraitCov=T)
# #Here it blows up. Try different optimization
# CompareRates.multTrait(Arini51sp,scaledcallparam[,1:4],TraitCov=T)
#
# source('CompareRatesAmongTraits2.r')
# CompareRates.multTrait2(Arini51sp,scaledcallparam[,1:4],TraitCov=T)
# #TAKES A LONG TIME to converge. BUT, results were gibberish. (Rconstrained matrix wasn't constrained...)
#
#
# CompareRates.multTrait2(Arini51sp,scaledcallparam,TraitCov=T)
#
# write.csv(callparam2,"/home/m/Dropbox/Projects-Manuscripts/Manuscripts/Parrot Calls-Phylo paper/Results/Call parameters parrots.csv")
#
#
# #read data
# callparam<-read.csv("Call parameters parrots.csv")
#
# #make it a matrix with species names in rows
# callparammat<-as.matrix(callparam[,2:9])
# rownames(callparammat)<-callparam$X
#
# #scale due to different units
# scaledcallparam<-scale(callparammat)
#
# #read tree
# Arini51sp<-read.tree("Arini51spsep14.tre")
#
# #check match among data and tree
# library(geiger)
# name.check(Arini51sp, scaledcallparam)
#
# #compare evo rates
# CompareRates.multTrait2(Arini51sp,scaledcallparam)
#
#
# CompareRates.multTrait(Arini51sp,as.matrix(sortedvars[,c(1:8)]*100))
#
# CompareRates.multTrait(Arini51sp,as.matrix(sortedvars[,c(9:11)]*10))
#
# CompareRates.multTrait(Arini51sp,as.matrix(sortedvars[,c(9:10)]*10))
# CompareRates.multTrait(Arini51sp,as.matrix(sortedvars[,c(10:11)]*10))
# CompareRates.multTrait(Arini51sp,as.matrix(sortedvars[,c(9,11)]*10))
Divergence times??
Sympatric species (overlap of smallest range species >20%)
## [1] 2
Allopatric species (overlap of smallest range species <20%)
## [1] 15
