### SETUP ### # clear workspace rm(list=ls()) # load packages library("caper") library("phytools") library("viridis") ### WEAVERBIRDS ### # read data weaverbird_data<-read.csv("weaverbird_data.csv", stringsAsFactors=T) str(weaverbird_data) # 59 species, 15 variables # create combined developmental period variable weaverbird_data$Dev_period<-weaverbird_data$Incubation_days+weaverbird_data$Nestling_days # descriptives/exploratory plots table(weaverbird_data$Tunnel, useNA="always") table(weaverbird_data$Attachment, useNA="always") table(weaverbird_data$Colonial, useNA="always") table(weaverbird_data$Thorns, useNA="always") table(weaverbird_data$Stinging, useNA="always") table(weaverbird_data$Raptors, useNA="always") table(weaverbird_data$Water, useNA="always") par(mfrow=c(2,4)) hist(weaverbird_data$Max_tunnel_length_cm, breaks=20, main="") hist(weaverbird_data$Height_m, breaks=20, main="") hist(weaverbird_data$Body_mass_g, breaks=20, main="") hist(weaverbird_data$Incubation_days, breaks=20, main="") hist(weaverbird_data$Nestling_days, breaks=20, main="") hist(weaverbird_data$Dev_period, breaks=20, main="") hist(abs(weaverbird_data$Latitude), breaks=20, main="") par(mfrow=c(2,4)) hist(log10(weaverbird_data$Max_tunnel_length_cm), breaks=20, main="") hist(log10(weaverbird_data$Height_m), breaks=20, main="") hist(log10(weaverbird_data$Body_mass_g), breaks=20, main="") hist(log10(weaverbird_data$Incubation_days), breaks=20, main="") hist(log10(weaverbird_data$Nestling_days), breaks=20, main="") hist(log10(weaverbird_data$Dev_period), breaks=20, main="") hist(log10(abs(weaverbird_data$Latitude)), breaks=20, main="") par(mfrow=c(1,1)) plot(weaverbird_data[,c("Max_tunnel_length_cm", "Height_m", "Body_mass_g", "Incubation_days", "Nestling_days", "Dev_period","Latitude")]) par(mfrow=c(2,2)) plot(Dev_period~Body_mass_g, weaverbird_data) plot(log10(Dev_period)~log10(Body_mass_g), weaverbird_data) plot(Dev_period~Max_tunnel_length_cm, weaverbird_data) plot(log10(Dev_period)~log10(Max_tunnel_length_cm), weaverbird_data) # read phylogeny weaverbird_tree<-read.nexus("deSilva2019.tre") str(weaverbird_tree) # 108 species par(mfrow=c(1,1)) plot(weaverbird_tree, cex=0.5) # combine data and tree into a comparative object for pgls analyses weaverbird_comp<-comparative.data(weaverbird_tree, weaverbird_data, names.col="Binomial", vcv=T, na.omit=F, warn.dropped=T) weaverbird_comp$dropped # all species in data match with phylogeny # pgls model: nest tunnels & developmental period weaverbird_tunnel_pgls_dev<-pgls(log10(Dev_period)~Tunnel+log10(Body_mass_g), data=weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_pgls_dev) length(fitted(weaverbird_tunnel_pgls_dev)) 10^(1.324648+(0.095569*mean(log10(weaverbird_tunnel_pgls_dev$data$data$Body_mass_g)))) # estimated dev period for non-tunnel builders 10^((1.324648+0.019675)+(0.095569*mean(log10(weaverbird_tunnel_pgls_dev$data$data$Body_mass_g)))) # estimated dev time for tunnel builders par(mfrow=c(2,2)) plot(weaverbird_tunnel_pgls_dev) # diagnostics # pgls model: nest tunnels & incubation period weaverbird_tunnel_pgls_inc<-pgls(log10(Incubation_days)~Tunnel+log10(Body_mass_g), data=weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_pgls_inc) length(fitted(weaverbird_tunnel_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_pgls_inc) # diagnostics # pgls model: nest tunnels & nestling period weaverbird_tunnel_pgls_nest<-pgls(log10(Nestling_days)~Tunnel+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_pgls_nest) length(fitted(weaverbird_tunnel_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_pgls_nest) # diagnostics # pgls model: tunnel length & developmental period weaverbird_tunnel_length_pgls_dev<-pgls(log10(Dev_period)~log10(Max_tunnel_length_cm)+log10(Body_mass_g), data=weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_length_pgls_dev) length(fitted(weaverbird_tunnel_length_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_length_pgls_dev) # diagnostics # pgls model: tunnel length & incubation period weaverbird_tunnel_length_pgls_inc<-pgls(log10(Incubation_days)~log10(Max_tunnel_length_cm)+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_length_pgls_inc) length(fitted(weaverbird_tunnel_length_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_length_pgls_inc) # diagnostics # pgls model: tunnel length & nestling period weaverbird_tunnel_length_pgls_nest<-pgls(log10(Nestling_days)~log10(Max_tunnel_length_cm)+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_length_pgls_nest) length(fitted(weaverbird_tunnel_length_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_length_pgls_nest) # diagnostics # pgls model: nest attachment & developmental period weaverbird_comp$data$Attachment<-relevel(weaverbird_comp$data$Attachment, "Supported") # change reference level weaverbird_attachment_pgls_dev<-pgls(log10(Dev_period)~Attachment+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_attachment_pgls_dev) length(fitted(weaverbird_attachment_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_attachment_pgls_dev) # pgls model: nest attachment & incubation period weaverbird_attachment_pgls_inc<-pgls(log10(Incubation_days)~Attachment+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_attachment_pgls_inc) length(fitted(weaverbird_attachment_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_attachment_pgls_inc) # pgls model: nest attachment & nestling period weaverbird_attachment_pgls_nest<-pgls(log10(Nestling_days)~Attachment+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_attachment_pgls_nest) length(fitted(weaverbird_attachment_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_attachment_pgls_nest) # pgls model: nest height & developmental period weaverbird_height_pgls_dev<-pgls(log10(Dev_period)~log10(Height_m)+log10(Body_mass_g), data=weaverbird_comp, lambda="ML") summary(weaverbird_height_pgls_dev) length(fitted(weaverbird_height_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_height_pgls_dev) # pgls model: nest height & incubation period weaverbird_height_pgls_inc<-pgls(log10(Incubation_days)~log10(Height_m)+log10(Body_mass_g), data=weaverbird_comp, lambda="ML") summary(weaverbird_height_pgls_inc) length(fitted(weaverbird_height_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_height_pgls_inc) # pgls model: nest height & nestling period weaverbird_height_pgls_nest<-pgls(log10(Nestling_days)~log10(Height_m)+log10(Body_mass_g), data=weaverbird_comp, lambda="ML") summary(weaverbird_height_pgls_nest) length(fitted(weaverbird_height_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_height_pgls_nest) # pgls model: coloniality & developmental period weaverbird_coloniality_pgls_dev<-pgls(log10(Dev_period)~Colonial+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_coloniality_pgls_dev) length(fitted(weaverbird_coloniality_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_coloniality_pgls_dev) # pgls model: coloniality & incubation period weaverbird_coloniality_pgls_inc<-pgls(log10(Incubation_days)~Colonial+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_coloniality_pgls_inc) length(fitted(weaverbird_coloniality_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_coloniality_pgls_inc) # pgls model: coloniality & nestling period weaverbird_coloniality_pgls_nest<-pgls(log10(Nestling_days)~Colonial+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_coloniality_pgls_nest) length(fitted(weaverbird_coloniality_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_coloniality_pgls_nest) # pgls model: nesting in thorns & developmental period weaverbird_thorns_pgls_dev<-pgls(log10(Dev_period)~Thorns+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_thorns_pgls_dev) length(fitted(weaverbird_thorns_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_thorns_pgls_dev) # pgls model: nesting in thorns & incubation period weaverbird_thorns_pgls_inc<-pgls(log10(Incubation_days)~Thorns+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_thorns_pgls_inc) length(fitted(weaverbird_thorns_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_thorns_pgls_inc) # pgls model: nesting in thorns & nestling period weaverbird_thorns_pgls_nest<-pgls(log10(Nestling_days)~Thorns+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_thorns_pgls_nest) length(fitted(weaverbird_thorns_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_thorns_pgls_nest) # pgls model: nesting over water & developmental period weaverbird_water_pgls_dev<-pgls(log10(Dev_period)~Water+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_water_pgls_dev) length(fitted(weaverbird_water_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_water_pgls_dev) # pgls model: nesting over water & incubation period weaverbird_water_pgls_inc<-pgls(log10(Incubation_days)~Water+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_water_pgls_inc) length(fitted(weaverbird_water_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_water_pgls_inc) # pgls model: nesting over water, nest tunnel & incubation period weaverbird_water_tunnel_pgls_inc<-pgls(log10(Incubation_days)~Water+Tunnel+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_water_tunnel_pgls_inc) length(fitted(weaverbird_water_tunnel_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_water_tunnel_pgls_inc) # pgls model: nesting over water & nestling period weaverbird_water_pgls_nest<-pgls(log10(Nestling_days)~Water+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_water_pgls_nest) length(fitted(weaverbird_water_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_water_pgls_nest) # pgls model: nesting near stinging insects & developmental period weaverbird_sting_pgls_dev<-pgls(log10(Dev_period)~Stinging+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_sting_pgls_dev) length(fitted(weaverbird_sting_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_sting_pgls_dev) # pgls model: nesting near stinging insects & incubation period weaverbird_sting_pgls_inc<-pgls(log10(Incubation_days)~Stinging+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_sting_pgls_inc) length(fitted(weaverbird_sting_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_sting_pgls_inc) # pgls model: nesting near stinging insects & nestling period weaverbird_sting_pgls_nest<-pgls(log10(Nestling_days)~Stinging+log10(Body_mass_g), data= weaverbird_comp, lambda="ML") summary(weaverbird_sting_pgls_nest) length(fitted(weaverbird_sting_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_sting_pgls_nest) # pgls model: nest tunnels & developmental period (controlling for latitude) weaverbird_tunnel_lat_pgls_dev<-pgls(log10(Dev_period)~Tunnel+log10(Body_mass_g)+abs(Latitude), data=weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_lat_pgls_dev) length(fitted(weaverbird_tunnel_lat_pgls_dev)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_lat_pgls_dev) # diagnostics # pgls model: nest tunnels & incubation period (controlling for latitude) weaverbird_tunnel_lat_pgls_inc<-pgls(log10(Incubation_days)~Tunnel+log10(Body_mass_g)+abs(Latitude), data= weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_lat_pgls_inc) length(fitted(weaverbird_tunnel_lat_pgls_inc)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_lat_pgls_inc) # diagnostics # pgls model: nest tunnels & nestling period (controlling for latitude) weaverbird_tunnel_lat_pgls_nest<-pgls(log10(Nestling_days)~Tunnel+log10(Body_mass_g)+abs(Latitude), data= weaverbird_comp, lambda="ML") summary(weaverbird_tunnel_lat_pgls_nest) length(fitted(weaverbird_tunnel_lat_pgls_nest)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_lat_pgls_nest) # diagnostics # pgls model: developmental period (intercept only, estimating lambda only) weaverbird_dev_lambda<-pgls(log10(Dev_period)~1, data=weaverbird_comp, lambda="ML") summary(weaverbird_dev_lambda) length(fitted(weaverbird_dev_lambda)) # pgls model: incubation period (intercept only, estimating lambda only) weaverbird_inc_lambda<-pgls(log10(Incubation_days)~1, data=weaverbird_comp, lambda="ML") summary(weaverbird_inc_lambda) length(fitted(weaverbird_inc_lambda)) # pgls model: nestling period (intercept only, estimating lambda only) weaverbird_nest_lambda<-pgls(log10(Nestling_days)~1, data=weaverbird_comp, lambda="ML") summary(weaverbird_nest_lambda) length(fitted(weaverbird_nest_lambda)) # pgls model: nest tunnels & developmental period (lambda fixed to 1) weaverbird_tunnel_pgls_dev_lfix<-pgls(log10(Dev_period)~Tunnel+log10(Body_mass_g), data=weaverbird_comp, lambda=1) summary(weaverbird_tunnel_pgls_dev_lfix) AIC(weaverbird_tunnel_pgls_dev, weaverbird_tunnel_pgls_dev_lfix) # pgls model: nest tunnels & incubation period (lambda fixed to 1) weaverbird_tunnel_pgls_inc_lfix<-pgls(log10(Incubation_days)~Tunnel+log10(Body_mass_g), data=weaverbird_comp, lambda=1) summary(weaverbird_tunnel_pgls_inc_lfix) AIC(weaverbird_tunnel_pgls_inc, weaverbird_tunnel_pgls_inc_lfix) # pgls model: nest tunnels & nestling period (lambda fixed to 1) weaverbird_tunnel_pgls_nest_lfix<-pgls(log10(Nestling_days)~Tunnel+log10(Body_mass_g), data= weaverbird_comp, lambda=1) summary(weaverbird_tunnel_pgls_nest_lfix) AIC(weaverbird_tunnel_pgls_nest, weaverbird_tunnel_pgls_nest_lfix) # pgls model: developmental period (intercept only, estimating delta + lambda) weaverbird_dev_delta<-pgls(log10(Dev_period)~1, data=weaverbird_comp, lambda="ML", delta="ML") summary(weaverbird_dev_delta) length(fitted(weaverbird_dev_delta)) # pgls model: incubation period (intercept only, estimating delta + lambda) weaverbird_inc_delta<-pgls(log10(Incubation_days)~1, data=weaverbird_comp, lambda="ML", delta="ML") summary(weaverbird_inc_delta) length(fitted(weaverbird_inc_delta)) # pgls model: nestling period (intercept only, estimating delta + lambda) weaverbird_nest_delta<-pgls(log10(Nestling_days)~1, data=weaverbird_comp, lambda="ML", delta="ML") summary(weaverbird_nest_delta) length(fitted(weaverbird_nest_delta)) # pgls model: nest tunnels & developmental period (estimating delta + lambda) weaverbird_tunnel_pgls_dev_delta<-pgls(log10(Dev_period)~Tunnel+log10(Body_mass_g), data=weaverbird_comp, lambda="ML", delta="ML") summary(weaverbird_tunnel_pgls_dev_delta) par(mfrow=c(2,2)) plot(weaverbird_tunnel_pgls_dev_delta) # diagnostics # pgls model: nest tunnels & incubation period (estimating delta + lambda) weaverbird_tunnel_pgls_inc_delta<-pgls(log10(Incubation_days)~Tunnel+log10(Body_mass_g), data=weaverbird_comp, lambda="ML", delta="ML") summary(weaverbird_tunnel_pgls_inc_delta) length(fitted(weaverbird_tunnel_pgls_inc_delta)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_pgls_inc_delta) # diagnostics # pgls model: nest tunnels & nestling period (estimating delta + lambda) weaverbird_tunnel_pgls_nest_delta<-pgls(log10(Nestling_days)~Tunnel+log10(Body_mass_g), data= weaverbird_comp, lambda="ML", delta="ML") summary(weaverbird_tunnel_pgls_nest_delta) length(fitted(weaverbird_tunnel_pgls_nest_delta)) par(mfrow=c(2,2)) plot(weaverbird_tunnel_pgls_nest_delta) # diagnostics # plot: nest tunnels & developmental period par(mfrow=c(1,3)) plot(log10(Dev_period)~log10(Body_mass_g), subset(weaverbird_tunnel_pgls_dev$data$data, Tunnel=="No"), pch=19, col=viridis(n=1, begin=1, end=1), ylab="Log10 developmental period", xlab="Log10 body mass", las=1, ylim=c(range(log10(weaverbird_tunnel_pgls_dev$data$data$Dev_period))), xlim=c(range(log10(weaverbird_tunnel_pgls_dev$data$data$Body_mass_g)))) # plot points for non-tunnel builders points(log10(Dev_period)~log10(Body_mass_g), subset(weaverbird_tunnel_pgls_dev$data$data, Tunnel=="Yes"), pch=19, col=viridis(n=1, begin=0.1, end=0.1)) # add points for tunnel builders summary(weaverbird_tunnel_pgls_dev) abline(1.324648, 0.095569, col=viridis(n=1, begin=1, end=1), lwd=2) # fit line for non-tunnel builders abline(1.324648+0.019675, 0.095569, col=viridis(n=1, begin=0.1, end=0.1), lwd=2) # fit line for tunnel builders legend("bottomright", pch=19, col=c(viridis(n=1, begin=0.1, end=0.1), viridis(n=1, begin=1, end=1)), legend=c("Tunnel", "No tunnel")) # plot: nest tunnels & incubation period plot(log10(Incubation_days)~log10(Body_mass_g), subset(weaverbird_tunnel_pgls_inc$data$data, Tunnel=="No"), pch=19, col=viridis(n=1, begin=1, end=1), ylab="Log10 incubation period", xlab="Log10 body mass", las=1, ylim=c(range(log10(weaverbird_tunnel_pgls_inc$data$data$Incubation_days))), xlim=c(range(log10(weaverbird_tunnel_pgls_inc$data$data$Body_mass_g)))) points(log10(Incubation_days)~log10(Body_mass_g), subset(weaverbird_tunnel_pgls_inc$data$data, Tunnel=="Yes"), pch=19, col=viridis(n=1, begin=0.1, end=0.1)) summary(weaverbird_tunnel_pgls_inc) abline(1.0604656, 0.0386005, col=viridis(n=1, begin=1, end=1), lwd=2) abline(1.0604656 + 0.0213149, 0.0386005, col=viridis(n=1, begin=0.1, end=0.1), lwd=2) # plot: nest tunnels & nestling period plot(log10(Nestling_days)~log10(Body_mass_g), subset(weaverbird_tunnel_pgls_nest$data$data, Tunnel=="No"), pch=19, col=viridis(n=1, begin=1, end=1), ylab="Log10 nestling period", xlab="Log10 body mass", las=1, ylim=c(range(log10(weaverbird_tunnel_pgls_nest$data$data$Nestling_days))), xlim=c(range(log10(weaverbird_tunnel_pgls_nest$data$data$Body_mass_g)))) points(log10(Nestling_days)~log10(Body_mass_g), subset(weaverbird_tunnel_pgls_nest$data$data, Tunnel=="Yes"), pch=19, col=viridis(n=1, begin=0.1, end=0.1)) summary(weaverbird_tunnel_pgls_nest) abline(1.006947, 0.135535, col=viridis(n=1, begin=1, end=1), lwd=2) abline(1.006947 + 0.017491, 0.135535, col=viridis(n=1, begin=0.1, end=0.1), lwd=2) # plot: tunnel length & developmental period par(mfrow=c(1,3)) plot(log10(Dev_period)~log10(Max_tunnel_length_cm), weaverbird_tunnel_length_pgls_dev$data$data, las=1, pch=19, col="grey", ylab="Log10 developmental period", xlab="Log10 max tunnel length") # plot: tunnel length & incubation period plot(log10(Incubation_days)~log10(Max_tunnel_length_cm), weaverbird_tunnel_length_pgls_inc$data$data, las=1, pch=19, col="grey", ylab="Log10 incubation period", xlab="Log10 max tunnel length") # plot: tunnel length & nestling period plot(log10(Nestling_days)~log10(Max_tunnel_length_cm), weaverbird_tunnel_length_pgls_nest$data$data, las=1, pch=19, col="grey", ylab="Log10 nestling period", xlab="Log10 max tunnel length") # plot: phylogeny + ancestral states reconstruction for nest tunnels weaverbird_ASR_data<-weaverbird_tunnel_pgls_dev$data$data$Tunnel # extract data for plotting names(weaverbird_ASR_data)<-rownames(weaverbird_tunnel_pgls_dev$data$data) # assign species names as rownames (to avoid risk of dissocation) weaverbird_ASR_tree<-weaverbird_tunnel_pgls_dev$data$phy weaverbird_ER_tunnel<-ace(weaverbird_ASR_data, weaverbird_ASR_tree, type="discrete", model="ER") # run ancestral states reconstruction (equal rates model) tunnel_cols<-setNames(c(viridis(n=1, begin=1, end=1), viridis(n=1, begin=0.1, end=0.1)), levels(weaverbird_ASR_data)) # assign colours for plot par(mfrow=c(1,1)) # see: http://www.phytools.org/Cordoba2017/ex/8/Anc-states-discrete.html par(mar=c(0,0,0,0)) plot(weaverbird_ASR_tree, cex=0.75, label.offset=0.01) # plot tree tiplabels(pie=to.matrix(weaverbird_ASR_data[weaverbird_ASR_tree$tip.label], levels(weaverbird_ASR_data)), piecol=tunnel_cols, cex=0.4) # add data as tip labels nodelabels(node=1:weaverbird_ASR_tree$Nnode+Ntip(weaverbird_ASR_tree), pie= weaverbird_ER_tunnel$lik.anc, piecol=tunnel_cols, cex=0.4) # add results of ancestral states model as node labels legend("topright", legend=c("Tunnel", "No tunnel"), pch=19, col=c(viridis(n=1, begin=0.1, end=0.1), viridis(n=1, begin=1, end=1))) # phylogenetic signal for tunnel presence (D statistic) weaverbird_D_data<-weaverbird_tunnel_pgls_dev$data$data weaverbird_D_data$Binomial<-rownames(weaverbird_D_data) weaverbird_D_data<-weaverbird_D_data[,c("Binomial", "Tunnel")] weaverbird_D_phy<-weaverbird_tunnel_pgls_dev$data$phy tunnel_D<-phylo.d(weaverbird_D_data, weaverbird_D_phy, names.col=Binomial, binvar=Tunnel) ### ICTERIDS ### # load data icterid_data<-read.csv("icterid_data.csv", stringsAsFactors=T) str(icterid_data) # 54 species, 14 variables # create combined developmental period variable icterid_data$Dev_period<-icterid_data$Incubation_days+icterid_data$Nestling_days # descriptives/exploratory plots table(icterid_data$Nest_type, useNA="always") table(icterid_data$Colonial, useNA="always") table(icterid_data$Thorns, useNA="always") table(icterid_data$Stinging, useNA="always") table(icterid_data$Raptors, useNA="always") table(icterid_data$Water, useNA="always") par(mfrow=c(2,4)) hist(icterid_data$Max_length_cm, breaks=20, main="") hist(icterid_data$Height_m, breaks=20, main="") hist(icterid_data$Body_mass_g, breaks=20, main="") hist(icterid_data$Incubation_days, breaks=20, main="") hist(icterid_data$Nestling_days, breaks=20, main="") hist(icterid_data$Dev_period, breaks=20, main="") hist(abs(icterid_data$Latitude), breaks=20, main="") par(mfrow=c(2,4)) hist(log10(icterid_data$Max_length_cm), breaks=20, main="") hist(log10(icterid_data$Height_m), breaks=20, main="") hist(log10(icterid_data$Body_mass_g), breaks=20, main="") hist(log10(icterid_data$Incubation_days), breaks=20, main="") hist(log10(icterid_data$Nestling_days), breaks=20, main="") hist(log10(icterid_data$Dev_period), breaks=20, main="") hist(log10(abs(icterid_data$Latitude)), breaks=20, main="") par(mfrow=c(1,1)) plot(icterid_data[,c("Max_length_cm", "Height_m", "Body_mass_g", "Incubation_days", "Nestling_days", "Dev_period", "Latitude")]) par(mfrow=c(2,2)) plot(Dev_period~Body_mass_g, icterid_data) plot(log10(Dev_period)~log10(Body_mass_g), icterid_data) plot(Dev_period~Max_length_cm, icterid_data) plot(log10(Dev_period)~log10(Max_length_cm), icterid_data) # read tree icterid_tree<-read.nexus("Powell2014.nex") str(icterid_tree) # 118 species icterid_tree<-multi2di(icterid_tree) # randomly resolve polytomies par(mfrow=c(1,1)) plot(icterid_tree, cex=0.5) # combine data and tree into a comparative object for pgls analyses icterid_comp<-comparative.data(icterid_tree, icterid_data, names.col="Binomial", vcv=T, na.omit=F, warn.dropped=T) icterid_comp$dropped # all species in data match with phylogeny # pgls model: nest type & developmental period icterid_comp$data$Nest_type<-relevel(icterid_comp$data$Nest_type, "Supported") # change reference level icterid_nest_type_pgls_dev<-pgls(log10(Dev_period)~Nest_type+log10(Body_mass_g), data=icterid_comp, lambda="ML") summary(icterid_nest_type_pgls_dev) length(fitted(icterid_nest_type_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_nest_type_pgls_dev) 10^(1.121135+(0.163288*mean(log10(icterid_nest_type_pgls_dev$data$data$Body_mass_g)))) # estimated dev period for supported nest-builders 10^((1.121135+0.073340)+(0.163288*mean(log10(icterid_nest_type_pgls_dev$data$data$Body_mass_g)))) # estimated dev period for suspended nest-builders 10^((1.121135+0.135994)+(0.163288*mean(log10(icterid_nest_type_pgls_dev$data$data$Body_mass_g)))) # estimated dev period for pendulous nest-builders # pgls model: nest type & incubation period icterid_nest_type_pgls_inc<-pgls(log10(Incubation_days)~Nest_type+log10(Body_mass_g), data=icterid_comp, lambda="ML") summary(icterid_nest_type_pgls_inc) length(fitted(icterid_nest_type_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_nest_type_pgls_inc) # pgls model: nest type & nestling period icterid_nest_type_pgls_nest<-pgls(log10(Nestling_days)~Nest_type+log10(Body_mass_g), data=icterid_comp, lambda="ML") summary(icterid_nest_type_pgls_nest) length(fitted(icterid_nest_type_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_nest_type_pgls_nest) # pgls model: nest length & developmental period icterid_nest_length_pgls_dev<-pgls(log10(Dev_period)~log10(Max_length_cm)+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_length_pgls_dev) length(fitted(icterid_nest_length_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_nest_length_pgls_dev) # pgls model: nest length & incubation period icterid_nest_length_pgls_inc<-pgls(log10(Incubation_days)~log10(Max_length_cm)+log10(Body_mass_g), data=icterid_comp, lambda="ML") summary(icterid_nest_length_pgls_inc) length(fitted(icterid_nest_length_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_nest_length_pgls_inc) # pgls model: nest length & nestling period icterid_nest_length_pgls_nest<-pgls(log10(Nestling_days)~log10(Max_length_cm)+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_length_pgls_nest) length(fitted(icterid_nest_length_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_nest_length_pgls_nest) # pgls model: nest height & developmental period icterid_nest_height_pgls_dev<-pgls(log10(Dev_period)~log10(Height_m+1)+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_height_pgls_dev) length(fitted(icterid_nest_height_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_nest_height_pgls_dev) # pgls model: nest height & incubation period icterid_nest_height_pgls_inc<-pgls(log10(Incubation_days)~log10(Height_m+1)+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_height_pgls_inc) length(fitted(icterid_nest_height_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_nest_height_pgls_inc) # pgls model: nest height & nestling period icterid_nest_height_pgls_nest<-pgls(log10(Nestling_days)~log10(Height_m+1)+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_height_pgls_nest) length(fitted(icterid_nest_height_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_nest_height_pgls_nest) # pgls model: nest height, type & developmental period icterid_nest_type_height_pgls_dev<-pgls(log10(Dev_period)~log10(Height_m+1)+Nest_type+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_type_height_pgls_dev) length(fitted(icterid_nest_type_height_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_nest_type_height_pgls_dev) # pgls model: nest height, type & incubation period icterid_nest_type_height_pgls_inc<-pgls(log10(Incubation_days)~log10(Height_m+1)+Nest_type+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_type_height_pgls_inc) length(fitted(icterid_nest_type_height_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_nest_type_height_pgls_inc) # pgls model: nest height, type & nestling period icterid_nest_type_height_pgls_nest<-pgls(log10(Nestling_days)~log10(Height_m+1)+Nest_type+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_nest_type_height_pgls_nest) length(fitted(icterid_nest_type_height_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_nest_type_height_pgls_nest) # pgls model: nest type & developmental period (controlling for latitude) icterid_comp$data$Nest_type<-relevel(icterid_comp$data$Nest_type, "Supported") # change reference level icterid_nest_type_lat_pgls_dev<-pgls(log10(Dev_period)~Nest_type+log10(Body_mass_g)+abs(Latitude), data=icterid_comp, lambda="ML") summary(icterid_nest_type_lat_pgls_dev) length(fitted(icterid_nest_type_lat_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_nest_type_lat_pgls_dev) # pgls model: nest type & incubation period (controlling for latitude) icterid_nest_type_lat_pgls_inc<-pgls(log10(Incubation_days)~Nest_type+log10(Body_mass_g)+abs(Latitude), data= icterid_comp, lambda="ML") summary(icterid_nest_type_lat_pgls_inc) length(fitted(icterid_nest_type_lat_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_nest_type_lat_pgls_inc) # pgls model: nest type & nestling period (controlling for latitude) icterid_nest_type_lat_pgls_nest<-pgls(log10(Nestling_days)~Nest_type+log10(Body_mass_g)+abs(Latitude), data= icterid_comp, lambda="ML") summary(icterid_nest_type_lat_pgls_nest) length(fitted(icterid_nest_type_lat_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_nest_type_lat_pgls_nest) # plot: nest type & developmental period par(mfrow=c(1,3)) plot(log10(Dev_period)~log10(Body_mass_g), subset(icterid_nest_type_pgls_dev$data$data, Nest_type=="Supported"), pch=19, col=viridis(n=1, begin=1, end=1), ylab="Log10 developmental period", xlab="Log10 body mass", las=1, ylim=c(range(log10(icterid_nest_type_pgls_dev $data$data$Dev_period))), xlim=c(range(log10(icterid_nest_type_pgls_dev $data$data$Body_mass_g)))) # plot points for supported nest builders points(log10(Dev_period)~log10(Body_mass_g), subset(icterid_nest_type_pgls_dev$data$data, Nest_type=="Suspended"), pch=19, col=viridis(n=1, begin=0.5, end=0.5)) # add points for suspended nest builders points(log10(Dev_period)~log10(Body_mass_g), subset(icterid_nest_type_pgls_dev$data$data, Nest_type=="Pendulous"), pch=19, col=viridis(n=1, begin=0.1, end=0.1)) # add points for pendulous nest builders summary(icterid_nest_type_pgls_dev) abline(1.121135, 0.163288,col=viridis(n=1, begin=1, end=1), lwd=2) # fit line for supported nest builders abline(1.121135+0.073340, 0.163288,col=viridis(n=1, begin=0.5, end=0.5), lwd=2) # fit line for suspended nest builders abline(1.121135+0.135994, 0.163288,col=viridis(n=1, begin=0.1, end=0.1), lwd=2) # fit line for pendulous nest builders legend("topleft", pch=19, col=c(viridis(n=1, begin=0.1, end=0.1), viridis(n=1, begin=0.5, end=0.5), col=viridis(n=1, begin=1, end=1)), legend=c("Pendulous", "Suspended", "Supported")) # plot: nest type & incubation period plot(log10(Incubation_days)~log10(Body_mass_g), subset(icterid_nest_type_pgls_inc$data$data, Nest_type=="Supported"), pch=19, col=viridis(n=1, begin=1, end=1), ylab="Log10 incubation period", xlab="Log10 body mass", las=1, ylim=c(range(log10(icterid_nest_type_pgls_inc $data$data$Incubation_days))), xlim=c(range(log10(icterid_nest_type_pgls_inc$data$data$Body_mass_g)))) points(log10(Incubation_days)~log10(Body_mass_g), subset(icterid_nest_type_pgls_inc$data$data, Nest_type=="Suspended"), pch=19, col=viridis(n=1, begin=0.5, end=0.5)) points(log10(Incubation_days)~log10(Body_mass_g), subset(icterid_nest_type_pgls_inc$data$data, Nest_type=="Pendulous"), pch=19, col=viridis(n=1, begin=0.1, end=0.1)) summary(icterid_nest_type_pgls_inc) abline(0.949080, 0.099200,col=viridis(n=1, begin=1, end=1), lwd=2) abline(0.949080 + 0.026996, 0.099200,col=viridis(n=1, begin=0.5, end=0.5), lwd=2) abline(0.949080 + 0.062251, 0.099200,col=viridis(n=1, begin=0.1, end=0.1), lwd=2) # plot: nest type & nestling period plot(log10(Nestling_days)~log10(Body_mass_g), subset(icterid_nest_type_pgls_nest$data$data, Nest_type=="Supported"), pch=19, col=viridis(n=1, begin=1, end=1), ylab="Log10 nestling period", xlab="Log10 body mass", las=1, ylim=c(range(log10(icterid_nest_type_pgls_nest$data$data$Nestling_days))), xlim=c(range(log10(icterid_nest_type_pgls_nest $data$data$Body_mass_g)))) points(log10(Nestling_days)~log10(Body_mass_g), subset(icterid_nest_type_pgls_nest$data$data, Nest_type=="Suspended"), pch=19, col=viridis(n=1, begin=0.5, end=0.5)) points(log10(Nestling_days)~log10(Body_mass_g), subset(icterid_nest_type_pgls_nest$data$data, Nest_type=="Pendulous"), pch=19, col=viridis(n=1, begin=0.1, end=0.1)) summary(icterid_nest_type_pgls_nest) abline(0.742685, 0.193606,col=viridis(n=1, begin=1, end=1), lwd=2) abline(0.742685 + 0.114934, 0.193606,col=viridis(n=1, begin=0.5, end=0.5), lwd=2) abline(0.742685 + 0.223820, 0.193606,col=viridis(n=1, begin=0.1, end=0.1), lwd=2) # plot: nest length & developmental period par(mfrow=c(1,3)) plot(log10(Dev_period)~log10(Max_length_cm), icterid_nest_length_pgls_dev$data$data, las=1, pch=19, col="grey", ylab="Log10 developmental period", xlab="Log10 max nest length") summary(icterid_nest_length_pgls_dev) icterid_nest_length_fitted<-1.127868+(log10(icterid_nest_length_pgls_dev$data$data$Max_length_cm)*0.061136)+(mean(log10(icterid_nest_length_pgls_dev$data$data$Body_mass_g))*0.163559) points(icterid_nest_length_fitted~log10(icterid_nest_length_pgls_dev$data$data$Max_length_cm), type='l', lty=1, lwd=2) # plot: nest length & incubation period plot(log10(Incubation_days)~log10(Max_length_cm), icterid_nest_length_pgls_inc$data$data, las=1, pch=19, col="grey", ylab="Log10 incubation period", xlab="Log10 max nest length") # plot: nest length & nestling period plot(log10(Nestling_days)~log10(Max_length_cm), icterid_nest_length_pgls_nest$data$data, las=1, pch=19, col="grey", ylab="Log10 nestling period", xlab="Log10 max nest length") # pgls model: colonality & developmental period icterid_coloniality_pgls_dev<-pgls(log10(Dev_period)~Colonial+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_coloniality_pgls_dev) length(fitted(icterid_coloniality_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_coloniality_pgls_dev) # pgls model: colonality & incubation period icterid_coloniality_pgls_inc<-pgls(log10(Incubation_days)~Colonial+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_coloniality_pgls_inc) length(fitted(icterid_coloniality_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_coloniality_pgls_inc) # pgls model: colonality & nestling period icterid_coloniality_pgls_nest<-pgls(log10(Nestling_days)~Colonial+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_coloniality_pgls_nest) length(fitted(icterid_coloniality_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_coloniality_pgls_nest) # pgls model: nesting in thorns & developmental period icterid_thorns_pgls_dev<-pgls(log10(Dev_period)~Thorns+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_thorns_pgls_dev) length(fitted(icterid_thorns_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_thorns_pgls_dev) # pgls model: nesting in thorns & incubation period icterid_thorns_pgls_inc<-pgls(log10(Incubation_days)~Thorns+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_thorns_pgls_inc) length(fitted(icterid_thorns_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_thorns_pgls_inc) # pgls model: nesting in thorns & nestling period icterid_thorns_pgls_nest<-pgls(log10(Nestling_days)~Thorns+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_thorns_pgls_nest) length(fitted(icterid_thorns_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_thorns_pgls_nest) # pgls model: nesting over water & developmental period icterid_water_pgls_dev<-pgls(log10(Dev_period)~Water+log10(Body_mass_g), data=icterid_comp, lambda="ML") summary(icterid_water_pgls_dev) length(fitted(icterid_water_pgls_dev)) par(mfrow=c(2,2)) plot(icterid_water_pgls_dev) # pgls model: nesting over water & incubation period icterid_water_pgls_inc<-pgls(log10(Incubation_days)~Water+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_water_pgls_inc) length(fitted(icterid_water_pgls_inc)) par(mfrow=c(2,2)) plot(icterid_water_pgls_inc) # pgls model: nesting over water & nestling period icterid_water_pgls_nest<-pgls(log10(Nestling_days)~Water+log10(Body_mass_g), data= icterid_comp, lambda="ML") summary(icterid_water_pgls_nest) length(fitted(icterid_water_pgls_nest)) par(mfrow=c(2,2)) plot(icterid_water_pgls_nest) # plot traits on the phylogeny + ancestral states icterid_ASR_data<-icterid_nest_type_pgls_dev$data$data$Nest_type names(icterid_ASR_data)<-rownames(icterid_nest_type_pgls_dev$data$data) icterid_ASR_tree<-icterid_nest_type_pgls_dev$data$phy icterid_ER_nest<-ace(icterid_ASR_data, icterid_ASR_tree, type="discrete", model="ER") nest_cols<-setNames(c(viridis(n=1, begin=1, end=1), viridis(n=1, begin=0.1, end=0.1), viridis(n=1, begin=0.5, end=0.5)), levels(icterid_ASR_data)) par(mfrow=c(1,1)) par(mar=c(0,0,0,0)) plot(icterid_ASR_tree, cex=0.75, label.offset=0.015) tiplabels(pie=to.matrix(icterid_ASR_data[icterid_ASR_tree$tip.label], levels(icterid_ASR_data)), piecol= nest_cols, cex=0.4) nodelabels(node=1: icterid_ASR_tree$Nnode+Ntip(icterid_ASR_tree), pie=icterid_ER_nest$lik.anc, piecol=nest_cols, cex=0.4) legend("bottomright", legend=c("Pendulous", "Suspended", "Supported"), pch=19, col=c(viridis(n=1, begin=0.1, end=0.1), viridis(n=1, begin=0.5, end=0.5), viridis(n=1, begin=1, end=1))) # phylogenetic signal for nest type (D statistic) icterid_D_data<-icterid_nest_type_pgls_dev$data$data icterid_D_data$Binomial<-rownames(icterid_D_data) icterid_D_data <-icterid_D_data[,c("Binomial", "Nest_type")] icterid_D_data$Nest_type_binary<-ifelse(icterid_D_data$Nest_type=="Supported", "No", "Yes") icterid_D_phy<-icterid_nest_type_pgls_dev$data$phy nest_type_D<-phylo.d(icterid_D_data, icterid_D_phy, names.col=Binomial, binvar=Nest_type_binary)