# Syntax employed in analyses of Rezende et al. (2017) Shrinking dinosaurs and the evolution of endothermy in birds # Cross-validation of the tip data and phylogenetic structure reported by Lee et al. (2014) and Benson et al. (2014) # Ancestral estimates could not be compared because they were not estimated in Benson et al. (2014) (Roger Benson, pers. comm.) # ----------------------------------------------------------------------------------------------------------------------------------- # Comparison between body mass estimates obtained by Lee and Benson # Whereas Lee employs femur length to estimate body mass, Benson employ equations based onlong bone circumference data <- read.table("Lee vs Benson.txt",header=TRUE) data$mass_lee <- 10^(-6.288 + 3.222*log10(data$femur_lee)) # There is a typo in Lee's dataset for Piatnitzkysaurus floresi (femur length = 152 instead of 552) # This species is not included in the Dataset by Rezende et al, but was fixed here to compare mass estimates obtained with different methods piat <- 10^(-6.288 + 3.222*log10(552)) # Plot data plot(data$mass_benson,data$mass_lee,log="xy",ylab="Mass Lee (kg)",xlab="Mass Benson (kg)",lwd=1.5); abline(0,1) points(data$mass_benson[286],piat,col="red",lwd=1.5) arrows(data$mass_benson[286], data$mass_lee[286],data$mass_benson[286],piat) text(0.5,1000,"R-square (correcting typo) = 0.9778") data$mass_lee[286] <- piat summary(lm(log10(data$mass_lee) ~ log10(data$mass_benson))) # ----------------------------------------------------------------------------------------------------------------------------------- # Comparison between phylogenies # Opening phylogeneyic trees library(ape) tree.lee <- read.tree("Lee_Theropod_Phylogeny.txt") tree.benson <- read.tree("Benson_Dinosaur_Phylogeny.txt") tree.benson <- tree.benson[[1]] # Selecting species in common in both trees tree.benson1 <- drop.tip(tree.benson,which(is.na(match(tree.benson$tip.label,tree.lee$tip.label)))) tree.lee1 <- drop.tip(tree.lee,which(is.na(match(tree.lee$tip.label,tree.benson$tip.label)))) # Calculating Baker's Gamma Index library(dendextend) a <- as.hclust(compute.brlen(tree.lee1)) b <- as.hclust(multi2di (compute.brlen(tree.benson1))) gamma <- cor_bakers_gamma(as.dendrogram(a),as.dendrogram(b)) gamma # Null model to test significance of dendogram similarity # Shuffles the tip data of b R <- 10 cor_bakers_gamma_results <- numeric(R) dend_mixed <- as.dendrogram(b) for(i in 1:R) { dend_mixed <- sample.dendrogram(dend_mixed, replace = FALSE) cor_bakers_gamma_results[i] <- cor_bakers_gamma(a, dend_mixed) } P.value <- round(sum(gamma < cor_bakers_gamma_results)/ R, 4) P.value