Hummer color evolution dec-14

http://marceloarayasalas.weebly.com/

Organize data

setwd("/home/m/Dropbox/Projects-Manuscripts/Projects/Hummer song evol/Archivos Paulina/datos")

maxir<-read.csv("pmaxbcgnrsininfrared2-2.csv")

colnames(maxir)[4]<-"wl"
colnames(maxir)[2]<-"wl.smooth"
# str(maxir)
range(maxir$wl)

## [1] 196.183 616.168

range(maxir$wl.smooth)

## [1] 195.8077 612.4068

gen<-as.factor(sapply(strsplit(colnames(maxir), ".",fixed=T), "[",1))
sp<-as.factor(sapply(strsplit(colnames(maxir), ".",fixed=T), "[",2))

species<-paste(gen,sp)[5:ncol(maxir)]  

library(pavo)

## Loading required package: rgl

x<-0
for(i in unique(maxir$Body)) {
  x<-x+1
  sub<-maxir[maxir$Body==i,c(4,5:ncol(maxir))]
sub<-as.rspec(sub,whichwl=1,lim = c(200,600),interp = F)
  if(x==1) aggdata<-aggspec(sub,paste(species,i,sep=".")) else {
    aggsub<-aggspec(sub,paste(species,i,sep="."))
    aggdata<-merge(aggdata,aggsub, by = "wl")
    }
}

aggdata<-as.rspec(aggdata,lim = c(300,700))

## wavelengths found in column 1

## Warning in as.rspec(aggdata, lim = c(300, 700)): Specified wavelength
## limits outside of actual data. Check 'lim' argument.

aggdata<-procspec(aggdata,opt = "smooth")

## processing options applied:
##  smoothing spectra with a span of 0.25 
##

range(aggdata$wl)

## [1] 300 700

convert to avian visual model

vis.aggdata <- vismodel(aggdata, visual='avg.uv',relative = F,bkg = maxir[maxir$Body=="c","white"])

tcsdata<-tcs(vis.aggdata)

## Warning in tcs(vis.aggdata): Quantum catch are not relative, and have been
## transformed

tcsdata<-data.frame("refl"=log(summary(aggdata)$B2),tcsdata)

belly<-tcsdata[grep(".b$",rownames(tcsdata)),]
colnames(belly)<-paste(colnames(belly),"b",sep=".")

crown<-tcsdata[grep(".c$",rownames(tcsdata)),]
colnames(crown)<-paste(colnames(crown),"c",sep=".")

gorgette<-tcsdata[grep(".g$",rownames(tcsdata)),]
colnames(gorgette)<-paste(colnames(gorgette),"g",sep=".")

neck<-tcsdata[grep(".n$",rownames(tcsdata)),]
colnames(neck)<-paste(colnames(neck),"n",sep=".")

rump<-tcsdata[grep(".r$",rownames(tcsdata)),]
colnames(rump)<-paste(colnames(rump),"r",sep=".")


# tcsplot(belly,col="red")
# tcspoints(neck,col="blue")
# tcspoints(rump,col="green")
# 
#put all data toghether
tcssortdata<-cbind(belly,crown,gorgette,neck,rump)


# voloverlap(neck,crown)
# dev.off()
# plot(aggdata,type="overlay")
# summary(tcsdata)
# str(aggdata)  
# colnames(aggdata)  

# sumdata<-summary(aggdata)
# # str(sumdata)
# sumdataF<-sumdata[sapply(sumdata,function(x) any(!is.na(x)))]
# str(sumdataF)

# xx<-cor(tcsdata,use="complete.obs",method="spearman")
# xx<-round(xx,2)
# xx[xx<0.7]<-"NA"
# xx<-as.matrix(xx)
# apply(xx,2,function(x) abs(length(x[x=="NA"])-(length(x)-1)))
#

Subset tree and data and change spp names

library(phytools)

## Loading required package: ape
## Loading required package: maps

library(geiger)
library(ape)

hum294<-read.tree("/home/m/Dropbox/Projects-Manuscripts/Projects/Hummer song evol/hum294.tre")
# hum294$tip.label
# hum63<-read.tree("/home/m/Dropbox/Projects-Manuscripts/Projects/Hummer song evol/Trochilidae63spp.fixed.names.tre")


library(gdata)

## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
## 
## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
## 
## Attaching package: 'gdata'
## 
## The following object is masked from 'package:stats':
## 
##     nobs
## 
## The following object is masked from 'package:utils':
## 
##     object.size

sp.names = read.xls("/home/m/Dropbox/Projects-Manuscripts/Projects/Hummer song evol/treelabels check sep-14.xlsx", sheet = "treelabels", header = TRUE)
# str(sp.names)
# sp.names$treelabs
# hum294$tip.label
# sp.names$Species.name

dat<-merge(as.data.frame(hum294$tip.label),sp.names,by.x=1,by.y=1,all.x=T,sort = F)
# View(dat)
colnames(dat)[1]<-"labs"
dtips<-as.character(dat$labs[grep("rop",dat$Drop.tip)])

sp291<-droplevels(dat[grep("rop",dat$Drop.tip,invert = T),])
# str(sp291)
hum291<-drop.tip(hum294,dtips)
hum286<-drop.tip(hum291,which(duplicated(sp291$Species.name)))
sp286<-sp291[which(!duplicated(sp291$Species.name)),]
# sp.names291<-sp.names[match(hum291$tip.label,sp.names291$treelabs),]
# # data.frame(sp.names291$treelabs[order(sp.names291$treelabs)],hum291$tip.label[order(hum291$tip.label)])
# data.frame(sp.names291$treelabs,hum291$tip.label)
# sp.names$treelabs[match(hum291$tip.label,as.character(sp.names291$treelabs))]
# hum291$tip.label[match(hum291$tip.label,as.character(sp.names291$treelabs))]

# which(sp.names291$treelabs=="Florisuga.mellivora.9646")
# View(sp.names291)

hum286$tip.label<-gsub(" ","_",sp286$Species.name)

tcssortdata<-cbind(belly,crown,gorgette,neck,rump)

gen<-as.factor(sapply(strsplit(rownames(tcssortdata), ".",fixed=T), "[",1))
sp<-as.factor(sapply(strsplit(rownames(tcssortdata), ".",fixed=T), "[",2))

rownames(tcssortdata)<-paste(gen,sp,sep="_")


# setdiff(hum286$tip.label,rownames(tcssortdata))
# setdiff(rownames(tcssortdata),hum286$tip.label)

# write.csv(setdiff(rownames(tcssortdata),hum286$tip.label),"/home/m/Dropbox/Projects-Manuscripts/Projects/Hummer song evol/delete.csv",row.names=F)
nn<-read.csv("/home/m/Dropbox/Projects-Manuscripts/Projects/Hummer song evol/delete.csv")

nn$nn.<-gsub(" ","_",nn$nn.)

# rownames(tcssortdata)

for(i in 1:nrow(nn)) rownames(tcssortdata)[rownames(tcssortdata)==nn$x[i]]<-as.character(nn$nn.[i])

rownames(tcssortdata)<-gsub(" ","_",rownames(tcssortdata))

# setdiff(rownames(tcssortdata),hum286$tip.label)


# name.check(hum286, rownames(tcssortdata))

hum59<-drop.tip(hum286,setdiff(hum286$tip.label,rownames(tcssortdata)))
# nrow(tcssortdata)
# str(sp286)
sp286$Species.name<-gsub(" ","_",sp286$Species.name,fixed=T)

EVI59X<-sp286[sp286$Species.name %in% rownames(tcssortdata),c("Species.name","map")]
EVI59<-as.numeric(as.character(droplevels(EVI59X[,2])))
names(EVI59)<-EVI59X$Species.name

name.check(hum59, tcssortdata)

## [1] "OK"

PCA

b<-tcssortdata[,grep(".b$",colnames(tcssortdata))]
# xx<-cor(b,use="complete.obs",method="spearman")
# xx<-round(xx,2)
# xx[xx<0.7]<-"NA"
# xx<-as.matrix(xx)
# a<-apply(xx,2,function(x) abs(length(x[x=="NA"])-(length(x)-1)))
# b1<-b[,a<2]
# bpca<-phyl.pca(hum59,b1, method="lambda", mode="corr")$S
bpcaX<-princomp(log(b[,-1]+1.1+abs(min(b))),cor=TRUE)
summary(bpcaX)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     2.7178456 2.0683035 1.5711543 1.24420595 0.53685027
## Proportion of Variance 0.4616678 0.2673675 0.1542829 0.09675303 0.01801301
## Cumulative Proportion  0.4616678 0.7290353 0.8833181 0.98007116 0.99808418
##                             Comp.6       Comp.7       Comp.8       Comp.9
## Standard deviation     0.157472153 0.0742433496 1.813132e-02 3.734995e-03
## Proportion of Variance 0.001549842 0.0003445047 2.054654e-05 8.718868e-07
## Cumulative Proportion  0.999634018 0.9999785228 9.999991e-01 9.999999e-01
##                             Comp.10      Comp.11      Comp.12      Comp.13
## Standard deviation     8.806252e-04 4.027925e-04 4.670313e-05 8.208851e-06
## Proportion of Variance 4.846879e-08 1.014011e-08 1.363239e-10 4.211577e-12
## Cumulative Proportion  1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
##                             Comp.14      Comp.15 Comp.16
## Standard deviation     4.166254e-06 4.272254e-07       0
## Proportion of Variance 1.084854e-12 1.140760e-14       0
## Cumulative Proportion  1.000000e+00 1.000000e+00       1

#number of PCs
npc<-ncol(bpcaX$scores)
bpca<-cbind("refl"=b[,1],bpcaX$scores[,1:npc])


c<-tcssortdata[,grep(".c$",colnames(tcssortdata))]
# xx<-cor(c,use="complete.obs",method="spearman")
# xx<-round(xx,2)
# xx[xx<0.7]<-"NA"
# xx<-as.matrix(xx)
# a<-apply(xx,2,function(x) abs(length(x[x=="NA"])-(length(x)-1)))
# c1<-c[,a<2]
# cpca<-phyl.pca(hum59,c1, method="lambda", mode="corr")$S
cpcaX<-princomp(log(c[,-1]+1.1+abs(min(c))),cor=TRUE)
summary(cpcaX)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     2.7178456 2.0683035 1.5711543 1.24420595 0.53685027
## Proportion of Variance 0.4616678 0.2673675 0.1542829 0.09675303 0.01801301
## Cumulative Proportion  0.4616678 0.7290353 0.8833181 0.98007116 0.99808418
##                             Comp.6       Comp.7       Comp.8       Comp.9
## Standard deviation     0.157472153 0.0742433496 1.813132e-02 3.734995e-03
## Proportion of Variance 0.001549842 0.0003445047 2.054654e-05 8.718868e-07
## Cumulative Proportion  0.999634018 0.9999785228 9.999991e-01 9.999999e-01
##                             Comp.10      Comp.11      Comp.12      Comp.13
## Standard deviation     8.806252e-04 4.027925e-04 4.670313e-05 8.208851e-06
## Proportion of Variance 4.846879e-08 1.014011e-08 1.363239e-10 4.211577e-12
## Cumulative Proportion  1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
##                             Comp.14      Comp.15 Comp.16
## Standard deviation     4.166254e-06 4.272254e-07       0
## Proportion of Variance 1.084854e-12 1.140760e-14       0
## Cumulative Proportion  1.000000e+00 1.000000e+00       1

cpca<-cbind("refl"=c[,1],cpcaX$scores[,1:npc])

g<-tcssortdata[,grep(".g$",colnames(tcssortdata))]
# xx<-cor(g,use="complete.obs",method="spearman")
# xx<-round(xx,2)
# xx[xx<0.7]<-"NA"
# xx<-as.matrix(xx)
# a<-apply(xx,2,function(x) abs(length(x[x=="NA"])-(length(x)-1)))
# g1<-g[,a<1]
# gpca<-phyl.pca(hum59,g1, method="lambda", mode="corr")$S
gpcaX<-princomp(log(g[,-1]+1.1+abs(min(g))),cor=TRUE)
summary(gpcaX)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3    Comp.4     Comp.5
## Standard deviation     2.8024962 2.0435268 1.4250381 1.2618368 0.50803531
## Proportion of Variance 0.4908741 0.2610001 0.1269208 0.0995145 0.01613124
## Cumulative Proportion  0.4908741 0.7518742 0.8787950 0.9783095 0.99444077
##                             Comp.6      Comp.7       Comp.8       Comp.9
## Standard deviation     0.250256447 0.142718237 0.0764970337 9.026928e-03
## Proportion of Variance 0.003914268 0.001273031 0.0003657373 5.092839e-06
## Cumulative Proportion  0.998355037 0.999628067 0.9999938047 9.999989e-01
##                             Comp.10      Comp.11      Comp.12      Comp.13
## Standard deviation     4.078333e-03 9.367555e-04 3.501412e-04 7.511454e-05
## Proportion of Variance 1.039550e-06 5.484443e-08 7.662430e-09 3.526371e-10
## Cumulative Proportion  9.999999e-01 1.000000e+00 1.000000e+00 1.000000e+00
##                             Comp.14      Comp.15      Comp.16
## Standard deviation     1.109324e-05 6.200831e-06 4.332712e-09
## Proportion of Variance 7.691252e-12 2.403144e-12 1.173275e-18
## Cumulative Proportion  1.000000e+00 1.000000e+00 1.000000e+00

gpca<-cbind("refl"=g[,1],gpcaX$scores[,1:npc])

n<-tcssortdata[,grep(".n$",colnames(tcssortdata))]
# xx<-cor(n,use="complete.obs",method="spearman")
# xx<-round(xx,2)
# xx[xx<0.7]<-"NA"
# xx<-as.matrix(xx)
# a<-apply(xx,2,function(x) abs(length(x[x=="NA"])-(length(x)-1)))
# n1<-n[,a<1]
# npca<-phyl.pca(hum59,n1, method="lambda", mode="corr")$S
npcaX<-princomp(log(n[,-1]+1.1+abs(min(n))),cor=TRUE)
summary(npcaX)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     2.7003861 2.2394951 1.4163609 1.19143033 0.49179983
## Proportion of Variance 0.4557553 0.3134586 0.1253799 0.08871914 0.01511669
## Cumulative Proportion  0.4557553 0.7692139 0.8945938 0.98331297 0.99842966
##                             Comp.6       Comp.7       Comp.8       Comp.9
## Standard deviation     0.148123430 0.0543423301 1.422046e-02 5.355980e-03
## Proportion of Variance 0.001371284 0.0001845681 1.263884e-05 1.792908e-06
## Cumulative Proportion  0.999800947 0.9999855148 9.999982e-01 9.999999e-01
##                             Comp.10      Comp.11      Comp.12      Comp.13
## Standard deviation     7.730077e-04 4.487827e-04 2.354174e-04 5.887804e-06
## Proportion of Variance 3.734630e-08 1.258787e-08 3.463835e-09 2.166640e-12
## Cumulative Proportion  1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
##                             Comp.14      Comp.15 Comp.16
## Standard deviation     2.823682e-06 2.076119e-06       0
## Proportion of Variance 4.983239e-13 2.693918e-13       0
## Cumulative Proportion  1.000000e+00 1.000000e+00       1

npca<-cbind("refl"=n[,1],npcaX$scores[,1:npc])

r<-tcssortdata[,grep(".r$",colnames(tcssortdata))]
# xx<-cor(r,use="complete.obs",method="spearman")
# xx<-round(xx,2)
# xx[xx<0.7]<-"NA"
# xx<-as.matrix(xx)
# a<-apply(xx,2,function(x) abs(length(x[x=="NA"])-(length(x)-1)))
# r1<-r[,a<2]
# rpca<-phyl.pca(hum59,r1[,c(1,3,4)], method="lambda", mode="corr")$S
rpcaX<-princomp(log(r[,-1]+1.1+abs(min(r))),cor=TRUE)
summary(cpcaX)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     2.7178456 2.0683035 1.5711543 1.24420595 0.53685027
## Proportion of Variance 0.4616678 0.2673675 0.1542829 0.09675303 0.01801301
## Cumulative Proportion  0.4616678 0.7290353 0.8833181 0.98007116 0.99808418
##                             Comp.6       Comp.7       Comp.8       Comp.9
## Standard deviation     0.157472153 0.0742433496 1.813132e-02 3.734995e-03
## Proportion of Variance 0.001549842 0.0003445047 2.054654e-05 8.718868e-07
## Cumulative Proportion  0.999634018 0.9999785228 9.999991e-01 9.999999e-01
##                             Comp.10      Comp.11      Comp.12      Comp.13
## Standard deviation     8.806252e-04 4.027925e-04 4.670313e-05 8.208851e-06
## Proportion of Variance 4.846879e-08 1.014011e-08 1.363239e-10 4.211577e-12
## Cumulative Proportion  1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
##                             Comp.14      Comp.15 Comp.16
## Standard deviation     4.166254e-06 4.272254e-07       0
## Proportion of Variance 1.084854e-12 1.140760e-14       0
## Cumulative Proportion  1.000000e+00 1.000000e+00       1

rpca<-cbind("refl"=r[,1],rpcaX$scores[,1:npc])

name.check(hum59, bpca)

## [1] "OK"

Multivariate phylogenetic signal

source("/home/m/Dropbox/Documentos R/Functions/Kmult test for phylo signal in mutivariate data.R")

PCs<-list(bpca,cpca,gpca,npca,rpca)
library(parallel)

physig<-mclapply(PCs,Test.Kmult,phy=hum59)

psr<-data.frame(c("belly","crown","gorget","neck","rump"),t(matrix(unlist(physig),ncol = 5,nrow = 2)))

colnames(psr)<-c("Body part","physig","p-value")

library(knitr)
kable(psr)

Body part	physig	p-value
belly	0.5187670	0.002
crown	0.5187670	0.001
gorget	0.3281156	0.027
neck	0.3590322	0.019
rump	0.3633574	0.013

Multivariate PGLS on PCs and log(reflectance)

evi<-data.frame(names(EVI59),EVI59)
EVI59Y<-merge(as.data.frame(rownames(bpca)),evi,1,1,sort = F)
l<-lapply(list("bpca","cpca","gpca","npca","rpca"),function(x) return(as.formula(paste(x,paste("~","EVI59")))))


source("/home/m/Dropbox/Documentos R/Functions/Dpgls multivariate pgls.R")
EVI59<-log(EVI59Y[,2]+1.1+abs(min(EVI59Y[,2])))
names(EVI59)<-EVI59Y[,1]
# D.pgls((rpca[,-1])~EVI59,phy=hum59)
results<-mclapply(l,D.pgls,phy=hum59)

setNames(results,c("belly","crown","gorget","neck","rump"))

## $belly
##          df     SS.obs       MS        F P.val        Rsq
## EVI59     1   1.551032 1.551032 0.883381 0.952 0.01526139
## Residual 57 100.080089 1.755791       NA    NA         NA
## 
## $crown
##          df     SS.obs       MS        F P.val        Rsq
## EVI59     1   1.551032 1.551032 0.883381 0.956 0.01526139
## Residual 57 100.080089 1.755791       NA    NA         NA
## 
## $gorget
##          df    SS.obs        MS        F P.val        Rsq
## EVI59     1  14.34771 14.347714 5.911857 0.644 0.09397048
## Residual 57 138.33550  2.426939       NA    NA         NA
## 
## $neck
##          df     SS.obs       MS        F P.val        Rsq
## EVI59     1   3.849519 3.849519 1.649793 0.906 0.02812956
## Residual 57 133.000072 2.333335       NA    NA         NA
## 
## $rump
##          df     SS.obs       MS        F P.val        Rsq
## EVI59     1   3.452901 3.452901 1.487881 0.947 0.02543914
## Residual 57 132.278924 2.320683       NA    NA         NA

Analisis using chroma,hue and brightness

# str(aggdata)
CHBdata<-summary(aggdata)[,c("B2","S8","H2")]
# rownames(CHBdata)

belly<-CHBdata[grep(".b$",rownames(CHBdata)),]
colnames(belly)<-paste(colnames(belly),"b",sep=".")

crown<-CHBdata[grep(".c$",rownames(CHBdata)),]
colnames(crown)<-paste(colnames(crown),"c",sep=".")

gorgette<-CHBdata[grep(".g$",rownames(CHBdata)),]
colnames(gorgette)<-paste(colnames(gorgette),"g",sep=".")

neck<-CHBdata[grep(".n$",rownames(CHBdata)),]
colnames(neck)<-paste(colnames(neck),"n",sep=".")

rump<-CHBdata[grep(".r$",rownames(CHBdata)),]
colnames(rump)<-paste(colnames(rump),"r",sep=".")

# 
#put all data toghether
chbsortdata<-cbind(belly,crown,gorgette,neck,rump)

gen<-as.factor(sapply(strsplit(rownames(chbsortdata), ".",fixed=T), "[",1))
sp<-as.factor(sapply(strsplit(rownames(chbsortdata), ".",fixed=T), "[",2))

rownames(chbsortdata)<-paste(gen,sp,sep="_")

nn<-read.csv("/home/m/Dropbox/Projects-Manuscripts/Projects/Hummer song evol/delete.csv")

nn$nn.<-gsub(" ","_",nn$nn.)


for(i in 1:nrow(nn)) rownames(chbsortdata)[rownames(chbsortdata)==nn$x[i]]<-as.character(nn$nn.[i])

rownames(chbsortdata)<-gsub(" ","_",rownames(chbsortdata))

name.check(hum59, chbsortdata)

## [1] "OK"

library(nlme)
pa<-corPagel(1,hum59)
x<-chbsortdata[,1]
names(x)<-rownames(chbsortdata)

# summary(gls(x~EVI59,correlation=pa,method="REML"))
hum59$edge.length<-hum59$edge.length*100

results<-try(sapply(chbsortdata[,],function(x) {r1<-summary(gls(log(x)~EVI59,correlation=pa,method="REML"))
return(round(as.data.frame(r1$tTable)[2,4],3))}),silent=T)

names(results)<-gsub("B2","Brightness",names(results))
names(results)<-gsub("S8","Chroma",names(results))
names(results)<-gsub("H2","Hue",names(results))
names(results)<-gsub(".b$",".Belly",names(results),fixed = F)
names(results)<-gsub(".n$",".Neck",names(results),fixed = F)
names(results)<-gsub(".c$",".Crown",names(results),fixed = F)
names(results)<-gsub(".r$",".Rump",names(results),fixed = F)
names(results)<-gsub(".g$",".Gorget",names(results),fixed = F)

print("P-values")

## [1] "P-values"

results

##  Brightness.Belly      Chroma.Belly         Hue.Belly  Brightness.Crown 
##             0.076             0.492             1.000             0.076 
##      Chroma.Crown         Hue.Crown Brightness.Gorget     Chroma.Gorget 
##             0.492             1.000             0.157             0.117 
##        Hue.Gorget   Brightness.Neck       Chroma.Neck          Hue.Neck 
##             0.000             0.006             0.037             0.108 
##   Brightness.Rump       Chroma.Rump          Hue.Rump 
##             0.000             0.046             0.984

Hummer color evolution dec-14

Marcelo Araya-Salas

12/27/2014

Organize data

convert to avian visual model

Subset tree and data and change spp names

PCA

Multivariate phylogenetic signal

Multivariate PGLS on PCs and log(reflectance)

Analisis using chroma,hue and brightness