--- title: "R Code for \"Body size correlates with discrete character morphological proxies of morphology\"" geometry: - scale=0.8 papersize: a4 fontsize: 12pt header-includes: - \usepackage{authblk} - \author{Tom Brougham} - \author{Nicol\'a{s} E. Campione} - \affil{School of Environmental and Rural Science, Universiy of New England, Armidale, NSW, Australia} - \definecolor{lightergray}{gray}{0.9} - \aboverulesep=0mm - \belowrulesep=0mm output: bookdown::pdf_document2: fig_caption: true toc: false linkcolor: blue urlcolor: blue --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # Introduction This document contains the annotated and commented `R` script used to generate the graphical and tabulated results presented in the main manuscript. The code is distributed as an RMarkdown (`.Rmd`) file, together with the associated nexus (`.nex`) files and body mass data in CSV format. A precompiled PDF accompanies the code and data; however, all analyses can be reproduced and compiled into a PDF, HTML, ODT, or Word document with RStudio () after all supplementary files have been placed into the same folder. See the RMarkdown website () for more information about setting up your R environment. The code is divided into two sections, reflecting the two separate analysis presented in the manuscript: 1) tests for correlation between morphology and body size; and 2) evaluation of the effects of size removal from hypothetica; avian stem lineage (ASL) taxa. # Required packages The code blocks below require installation of a selection of packages as prerequisites. The most recent versions of all packages except `Claddis` can be installed from the CRAN; `Claddis` will be installed from the Github development repository. ```{r install-pkgs, eval=FALSE, message=FALSE, warning=FALSE} install.packages(c("ade4", "ape", "devtools", "geomorph", "gtools", "knitr", "kableExtra", "magrittr", "MASS", "nlme", "paleotree", "phylobase", "phangorn", "qpcR", "RColorBrewer")) devtools::install_github("graemetlloyd/Claddis") ``` # Correlation between morphology and body size The first analysis starts by loading and log-transforming the body mass data for all dinosaurian taxa. The first and last appearance dates for each taxon is also loaded, which will be used for time scaling of phylogenetic trees for each matrix. A list of all `.nex` files in the current working directory is also obtained. ```{r read-data, warning=FALSE} mass <- as.data.frame(read.csv("dinosaurs-bm.csv", sep = ",", as.is = T, row.names = 1)) mass$mass <- log10(as.numeric(mass$mass)) taxa_fad_lad <- read.csv("dinosaurs-fad-lad.csv", sep=",", as.is = T, row.names = 1) nex_filenames <- list.files(pattern="*.nex$") ``` The list of file names will be iterated over; each `.nex` file will be read and have a distance matrix calculated for it using the Maximum Observable Rescaled Distance (MORD) metric, the default in `Claddis`. The Principal Coordinate Analysis (PCoA) follows, from which the first axis will be extracted from the results and merged with the body mass data for each taxon. A robust linear regression is then calculated from the merged dataframe. The implementation of robust linear regression (`rlm`) in `MASS` does not calculate an r-squared value, but the weights that are calculated for each taxon can be supplied as an additional argument into the standard linear regression function (`lm`), from which a weighted r-squared value is calculated. In addition, two separate analyses utilising phylogenetic information are performed. The first calculates Pearson correlation coefficients between body mass and the first principle coordinate axis for each character-taxon matrix given a phylogenetic tree. The second uses the phylogenetic generalised least squares (pGLS) method which allows a tree topology to be a covariate along with the body mass and PCo axes. The first such analysis is a regression of body mass with each individual principal coordinate axis from the set that was calculated for each of the nine discrete character-taxon matrices. The second is a regression of body mass with the complete set of principal coordinate axes for each character-taxon matrix. Finally, a Mantel test is performed which assesses the strength of the correlation between the MORD distance matrices of each discrete character-taxon matrix and the log-transformed body mass data for each taxa in each matrix, where possible. The results of the Mantel tests are presented to show the distribution of the correlation coefficients of the randomised permuted distance matrices with respect to that of the original, unpermuted distance matrices. A bivariate plot of the two distance matrices is also presented to assess the linearity of the relationship. The results are presented in Figs. \@ref(fig:pcoa-bm-plots)--\@ref(fig:dists-bivariate) and Table \@ref(tab:pcoa-bm-results). ```{r pcoa-bm-analysis, echo=FALSE, message=FALSE, warning=FALSE, results='hide'} all_results <- c() for (nex_filename in nex_filenames) { # Read in the phylogenetic matrix matrix <- Claddis::ReadMorphNexus(nex_filename) # Read in phylogenetic tree tree <- phylobase::readNexus(nex_filename, type = "tree") if (length(tree) > 1) { tree <- as(tree[[1]], "phylo") } else { tree <- as(tree, "phylo") } # Time-calibrate tree if edge lengths are not present if (!"edge.length" %in% names(tree)) { tree <- paleotree::timePaleoPhy(tree, taxa_fad_lad[rownames(taxa_fad_lad) %in% tree$tip.label, ], "equal", vartime = 1) } # Calculate the distance matrix using MORD distances <- Claddis::MorphDistMatrix(matrix) max_dists <- distances$DistanceMatrix # Trim taxa with NAs if necessary if (any(is.na(max_dists))) { trim <- Claddis::TrimMorphDistMatrix(max_dists) max_dists <- trim$DistMatrix } # Compute PCoA pcoa_result <- ape::pcoa(max_dists, correction = "cailliez") # Merge PCo and body mass data pcoa_bm <- merge(mass, pcoa_result$vectors.cor, by.x = 0, by.y = 0) rownames(pcoa_bm) <- pcoa_bm$Row.names pcoa_bm$Row.names <- NULL pcoa_bm$group <- NULL # Drop taxa that aren't in tree pcoa_bm <- pcoa_bm[rownames(pcoa_bm) %in% tree$tip.label, ] # Drop taxa that have no body mass estimates pcoa_bm <- pcoa_bm[!is.na(pcoa_bm$mass), ] # Drop taxa from tree that have no body mass estimates tree <- ape::drop.tip(tree, tree$tip.label[!tree$tip.label %in% rownames(pcoa_bm)]) # Calculate phylogenetic correlation between body mass and PCo1 corphy_bm_pco1 <- ape::corphylo(pcoa_bm[c(1, 2)], phy = tree) # Calculate robust linear regression and weighted r-squared for PCo1 pco1_rlm <- MASS::rlm(Axis.1 ~ mass, pcoa_bm, method = "MM", maxit = 1000) pco1_rlm <- lm(Axis.1 ~ mass, pco1_rlm$model, weights = pco1_rlm$w) # Perform phylogenetic GLS of body masses and PCo1 scores for all axes pco_n_pgls <- sapply(2:ncol(pcoa_bm), function (idx) { geomorph::procD.pgls(pco ~ mass, tree, data = list(pco = pcoa_bm[, idx], mass = pcoa_bm$mass)) }) # Perform phylogenetic GLS of body masses and ALL PCo scores pco_all_pgls <- geomorph::procD.pgls(pcoa_bm[, -1] ~ pcoa_bm[, 1], tree) # Perform Mantel test on morphological distance and body mass matrixes # Work out which taxa have body mass estimates bm_include <- rownames(pcoa_bm[!is.na(pcoa_bm$mass), ]) # Assemble the distance matrices and perform test trim_dist_matrix <- max_dists[bm_include, bm_include] dist_bm <- dist(mass[bm_include, "mass"]) attr(dist_bm, "Labels") <- bm_include trim_dist_matrix <- dist(trim_dist_matrix) mantel_result <- ade4::mantel.rtest(trim_dist_matrix, dist_bm, nrepet = 9999) dists_lm <- lm(dist_bm ~ trim_dist_matrix) # Accumulate results results <- list(pcoa = pcoa_bm, corphy_bm_pco1 = corphy_bm_pco1, pco1_rlm = pco1_rlm, pco1_pv = pcoa_result$values$Rel_corr_eig[1], pco_n_pgls = pco_n_pgls, pco_all_pgls = pco_all_pgls, mantel = mantel_result, dists_lm = dists_lm, fn = nex_filename) if (!length(all_results)) { all_results <- list(results) } else { all_results[length(all_results) + 1] <- list(results) } } ``` ```{r pcoa-bm-plots, echo=FALSE, fig.align='center', fig.cap="Scatter plots of PCo1 against log body mass (kg) for nine discrete character-taxon matrices of dinosaurs. Unfilled circles represent outgroup taxa. Robust linear regression lines are indicated in red.", fig.height=7, fig.width=7, warning=FALSE} # Define outgroups for each matrix and insert in results list outgroups <- list("brusatte-et-al_2010.nex" = c("Allosaurus_fragilis"), "brusatte-et-al_2014.nex" = c("Allosaurus_fragilis"), "cau_2018.nex" = c("Eodromaeus_murphi", "Guaibasaurus_candelariensis", "Herrerasaurus_ischigualastensis", "Heterodontosaurus_tucki", "Pisanosaurus_mertii", "Plateosaurus_engelhardti", "Saturnalia_tupiniquim", "Staurikosaurus_pricei"), "csiki-et-al_2010.nex" = c("Allosaurus_fragilis"), "currie-rogers_2005.nex" = c("Camarasaurus", "Diplodocus"), "longrich-et-al_2010.nex" = c("Thescelosaurus_neglectus"), "prieto-marquez_2010.nex" = c("Gilmoreosaurus_mongoliensis", "Bactrosaurus_johnsoni"), "sampson-et-al_2010.nex" = c("Leptoceratops_gracilis", "Protoceratops_andrewsi", "Turanoceratops_tardabilis"), "thompson-et-al_2012.nex" = c("Lesothosaurus_diagnosticus", "Scelidosaurus_harrisonii", "Scutellosaurus_lawleri", "Stegosaurus_stenops")) all_results <- lapply(all_results, function (x) { x$outgroups <- outgroups[x$fn][[1]] return(x) }) par(mfrow=c(3,3), mai=c(0.2, 0.2, 0.2, 0.2), oma=c(4, 4, 0, 0)) stat_table <- data.frame(row.names = c("df", "% PCo1 var.", "Phy. Corr.", "Slope.1", "r2.1", "Slope.2", "r2.2", "Statistic", "r2.3", "pco.all.aov r2")) idx <- 1 for (result in all_results) { plot_bg <- ifelse(rownames(result$pcoa) %in% result$outgroups, "white", "black") plot(result$pcoa$Axis.1 ~ result$pcoa$mass, pch=21, cex=0.8, xlab="", ylab="", frame.plot=F, col="black", bg=plot_bg) clip(min(result$pcoa$mass, na.rm = TRUE), max(result$pcoa$mass, na.rm = TRUE), min(result$pcoa$Axis.1), max(result$pcoa$Axis.1)) abline(result$pco1_rlm, col="red") title(chartr("123456789", "abcdefghi", as.character(idx)), adj=0.05, line=-0.5) idx <- idx + 1 pco1_rlm_stars <- gtools::stars.pval(coef(summary(result$pco1_rlm))[8])[1] pco1_r2 <- summary(result$pco1_rlm)$r.squared pco1_pgls_stars <- gtools::stars.pval(result$pco_n_pgls[,1]$aov.table$`Pr(>F)`[1])[1] mantel_stars <- gtools::stars.pval(result$mantel$pvalue)[1] pco_all_pgls_stars <- gtools::stars.pval(result$pco_all_pgls$aov.table$`Pr(>F)`[1])[1] stat_table[result$fn] = c(result$pco1_rlm$df.residual, signif(result$pco1_pv * 100, 2), signif(result$corphy_bm_pco1$cor.matrix[2], 2), paste(signif(result$pco1_rlm$coefficients[2], 2), pco1_rlm_stars), signif(pco1_r2, 2), paste(signif(result$pco_n_pgls[,1]$pgls.coefficients[2], 2), pco1_pgls_stars), signif(result$pco_n_pgls[,1]$aov.table$Rsq[1], 2), paste(signif(result$mantel$obs, 2), mantel_stars), signif(result$mantel$obs**2, 2), paste(signif(result$pco_all_pgls$aov.table$Rsq[1], 2), pco_all_pgls_stars)) } # Print the big axis labels mtext('Log body mass (kg)', side = 1, outer = TRUE, line = 2) mtext('PCo1', side = 2, outer = TRUE, line = 2) ``` ```{r pcoa-bm-results, echo=FALSE} library(magrittr) names(stat_table) <- c("a", "b", "c", "d", "e", "f", "g", "h", "i") stat_table <- t(stat_table) colnames(stat_table)[c(5, 7, 9, 10)] <- "r²" colnames(stat_table)[c(4, 6)] <- "Slope" knitr::kable(stat_table, caption = "The results of the robust regression and multivariate phylogenetic GLS for PCo1, and Mantel tests and multivariate phylogenetic GLS of all PCo axes, for nine discrete character-taxon matrices of dinosaurs.", booktabs = TRUE) %>% kableExtra::kable_styling(latex_options = c("striped", "scale_down"), stripe_color = "lightergray") %>% kableExtra::add_header_above(c(" " = 4, "Robust Linear Regression" = 2, "Phylogenetic GLS" = 2, "Mantel test" = 2, "Phylogenetic GLS" = 1)) ``` ```{r pco-all-pgls, fig.height=7, fig.align='center', fig.cap='Bar plots of the correllation coefficients of body mass against each principal coordinate axis for the nine discrete character-taxon matrices of dinosaurs. The red line represents the two-tailed 95% confidence interval.', fig.width=7, echo=FALSE} par(mfrow=c(3,3), mai=c(0.2, 0.2, 0.2, 0.2), oma=c(4, 4, 0, 0)) idx <- 1 for (result in all_results) { stdv <- qnorm(1 - 0.05/2)/sqrt(ncol(result$pco_n_pgls)) coeffs <- apply(result$pco_n_pgls, 2, function (n) {n$aov.table$Rsq[1]}) barplot(coeffs, col=gray(0.6)) clip(0,ncol(result$pco_n_pgls)*1.2, min(coeffs), max(coeffs)) abline(h=stdv, col="red") title(chartr("123456789", "abcdefghi", as.character(idx)), adj=0.05, line=0.5) idx <- idx + 1 } mtext('PCoA Axis', side = 1, outer = TRUE, line = 2) mtext('pGLS correlation coefficient', side = 2, outer = TRUE, line = 2) ``` ```{r mantel-plots, fig.height=7, fig.align='center', fig.cap="Histograms of the correlations of the permuted body mass and discrete character distance matrices for the nine discrete character-taxon matrices of dinosaurs, with the observed matrix correlation indicated by the black circle marker.", fig.width=7, echo=FALSE} par(mfrow=c(3,3), mai=c(0.2, 0.2, 0.2, 0.2), oma=c(4, 4, 0, 0)) # Plot Mantel tests for all axes idx <- 1 for (result in all_results) { plot(result$mantel$plot$hist, xlim=result$mantel$plot$xlim, main=NULL, col=gray(0.6)) y0 <- max(result$mantel$plot$hist$counts) lines(c(result$mantel$obs, result$mantel$obs), c(y0/2, 0)) points(result$mantel$obs, y0/2, pch = 16, cex = 2) title(chartr("123456789", "abcdefghi", as.character(idx)), adj=0.05, line=-0.5) idx <- idx + 1 } mtext('Distance', side = 1, outer = TRUE, line = 2) mtext('Frequency', side = 2, outer = TRUE, line = 2) ``` ```{r dists-bivariate, fig.height=7, fig.align='center', fig.cap="Kernel density plots of the discrete character distance matrices plotted against body mass for the nine discrete character-taxon dinosaur matrices. The black line represents the linear regression line.", fig.width=7, echo=FALSE} par(mfrow=c(3,3), mai=c(0.2, 0.2, 0.2, 0.2), oma=c(4, 4, 0, 0)) # Distance matrix bivariate plots idx <- 1 palette <- RColorBrewer::brewer.pal(8, "Set1") for (result in all_results) { dens <- MASS::kde2d(result$dists_lm$model$trim_dist_matrix, result$dists_lm$model$dist_bm, n=100) myPal <- colorRampPalette(c("white",palette[c(2,6,5,1)])) plot(result$dists_lm$model$trim_dist_matrix, result$dists_lm$model$dist_bm, type="n", frame.plot=FALSE) image(dens, col=myPal(256), add = TRUE) title(chartr("123456789", "abcdefghi", as.character(idx)), adj=0.05, line=-0.5) clip(min(result$dists_lm$model$trim_dist_matrix), max(result$dists_lm$model$trim_dist_matrix), min(result$dists_lm$model$dist_bm), max(result$dists_lm$model$dist_bm)) abline(result$dists_lm) idx <- idx + 1 } mtext('Discrete character distances', side = 1, outer = TRUE, line = 2) mtext('Body mass distances', side = 2, outer = TRUE, line = 2) ``` # Macroevolution case study The phylogenetic character-taxon matrix and time-calibrated Bayesian interence tree of Cau (2018) will be read. The body mass data is reused from the previous analysis, but truncated to the taxa present in the phylogenetic dataset and sorted in preparation for the analysis. Subsequently, taxa in the tree that do not have a body mass estimate will be pruned. ```{r} # Read in phylogenetic matrix morpho <- Claddis::ReadMorphNexus("cau_2018.nex") # Read in tree files tree <- phylobase::readNexus("cau_2018.nex", type="tree")[[1]] # Drop taxa that have no mass estimates prune_taxa <- rownames(mass[is.na(mass$mass),]) prune_taxa <- prune_taxa[prune_taxa %in% tree@label] tree_pruned <- phylobase::prune(tree, prune_taxa) mass <- mass[!is.na(mass$mass),] # Remove taxa from mass list that aren't in the tree mass <- mass[rownames(mass) %in% tree_pruned@label,"mass",drop=FALSE] # Reorder data to match the order of the tree's tip labels mass <- mass[order(match(rownames(mass),tree_pruned@label)), "mass", drop=FALSE] ``` The ancestral body size of each internal node within the pruned tree will then be reconstructed; the resulting list will be restricted to nodes that compose the ASL. ```{r warning=FALSE} # Define a function to list of all numbers and branch lengths of nodes # ancestral to a given taxon anc_nodes <- function (phy, taxon) { v.anc_nodes <- phylobase::ancestors(phy, taxon, type="ALL") v.anc.lengths <- head(phylobase::edgeLength(phy, v.anc_nodes), -1) return(list(nodes = v.anc_nodes, lengths = v.anc.lengths)) } asl_nodes <- anc_nodes(tree_pruned, "Meleagris_gallopavo") ages <- rev(cumsum(asl_nodes$lengths)) # Perform the ancestral state reconstruction for body mass anc <- phytools::fastAnc(as(tree_pruned, "phylo"), t(mass[1])) # Prune ancestral reconstructions to the ancestral nodes anc <- anc[names(anc) %in% asl_nodes$nodes] ``` Parsimony uninformative characters are removed from the matrix prior to discrete character ancestral state reconstruction. Ancestral states are calculated at each internal node of the original, unpruned tree for each character and compiled into a single ancestral state matrix. This matrix will then be restricted to those taxa that comprise the ASL and converted into a matrix. ```{r warning=FALSE} # Identify and remove uninformative characters uninformative <- c() for (char in 1:ncol(morpho$Matrix_1$Matrix)) { counts <- table(morpho$Matrix_1$Matrix[,char]) if (length(counts) < 2 || any(counts == 1)) { uninformative[length(uninformative) + 1] <- char } } new_matrix <- Claddis::MatrixPruner(morpho, characters2prune = uninformative) # Obtain list of nodes and branches along the avian stem lineage for the # complete tree asl_nodes_full <- anc_nodes(tree, "Meleagris_gallopavo") # Perform discrete character ancestral state reconstruction mtx_phydat <- phangorn::phyDat(new_matrix$Matrix_1$Matrix, type="USER", levels=c("0","1"), ambiguity=NA) anc_matrix <- phangorn::ancestral.pars(as(tree, "phylo"), mtx_phydat) # Limit the ancestral nodes to those that comprise the avian stem lineage anc_matrix <- anc_matrix[names(anc_matrix) %in% asl_nodes_full$nodes[-1],] # Convert the phydat object into an ancestral state matrix anc_matrix <- sapply(anc_matrix, as.character) anc_matrix[anc_matrix == "0.5"] <- "0/1" anc_matrix <- Claddis::MakeMorphMatrix(t(anc_matrix)) ``` As in the first analysis, a distance matrix will be calculated using MORD, followed by PCoA. The first principal coordinate axis (PCo1) is combined into a single dataframe together with the age and the reconstructed body mass of each hypothetical ancestor. Generalised least squares regression is used to assess the strength of the correlation between the PCo1 and the hypothetical ancestral body masses, assuming an AR1 autocorrelation structure using the ages of each hypothetical ancestor. The residuals of this correlation are plotted with the ages of the respective ancestors to provide the size-corrected morphological proxy. The corrected Akaike Information Criterion (AICc) of this model is then compared with a second intercept-only model that uses a constant body mass for each hypothetical ancestor. This allows a determination of whether the inclusion of body mass results in a better fit to the data. The results of the analyses are presented in Table. \@ref(tab:cau-asl-anc-pco1) and Fig. \@ref(fig:cau-asl-resid). ```{r message=FALSE} # Perform the principal coordinates analysis distances <- Claddis::MorphDistMatrix(anc_matrix) max_dists <- distances$DistanceMatrix pcoa_axes <- ape::pcoa(max_dists) # Merge the PCo1 values with the body mass data and plot the results pcoa_ages <- data.frame(ages = rev(cumsum(asl_nodes_full$lengths)), pcoa1 = pcoa_axes$vectors[,1]) pcoa_bm <- merge(pcoa_ages, data.frame(ages, anc), by="ages") # Calculate the GLS models, taking into account autocorrelation in the time # series pcoa_bm_gls <- nlme::gls(pcoa1 ~ anc, pcoa_bm, correlation = nlme::corAR1(form=~ages)) pcoa_bm_gls_null <- nlme::gls(pcoa1 ~ 1, pcoa_bm, correlation = nlme::corAR1(form=~ages)) ``` ```{r cau-asl-anc-pco1, echo=FALSE, result="hide", fig.align='center', fig.height=7, fig.width=6, out.width='\\textwidth'} gls_results <- list(pcoa_bm_gls, pcoa_bm_gls_null) %>% sapply(function (x) { c(signif(x$coefficients[1], 3), ifelse(is.na(x$coefficients[2]), "-", signif(x$coefficients[2], 3)), signif(x$logLik, 3), signif(qpcR::AICc(x), 3)) }) %>% t() %>% as.data.frame(row.names = c("Normal", "Null model")) names(gls_results) <- c("Intercept", "Slope", "Log Likelihood", "AICc") knitr::kable(gls_results, caption = "Results of generalised least squares analyses for hypothetical avian stem lineage taxa of the Cau (2018) matrix. The normal model uses ancestral body mass as the outcome variable of ancestral PCo1 values, whereas the null model uses a constant outcome variable (body mass) of one. Both analyses accommodate for autocorrelation along the time series using an autoregressive function (AR(1)) containing the ages of each hypothetical ancestor.", booktabs = TRUE) %>% kableExtra::kable_styling(latex_options = c("striped"), stripe_color = "lightergray") ``` ```{r cau-asl-resid, echo=FALSE, fig.align='center', fig.height=6, warning=FALSE, message=FALSE, out.width='\\textwidth', fig.cap="Re-analysis of the Cau (2018) ancestral avian morphospace. A line graph showing body mass, PCo1, and the residual of PCo1 vs body mass for each hypothetical taxon on the avian stem lineage, plotted over time."} par(mai=c(0, 0.8, 0.5, 1.6), mfcol=c(2, 1)) branch_cols <- rep("black", length(tree_pruned@edge.length)) branch_pal <- RColorBrewer::brewer.pal(5, "Set2") branch_cols[c(1:4, 164:180)] <- branch_pal[1] branch_cols[c(5:8, 128:163)] <- branch_pal[2] branch_cols[c(9:17, 67:127)] <- branch_pal[3] branch_cols[c(18:23, 51:66)] <- branch_pal[4] branch_cols[c(24:50)] <- branch_pal[5] plot(ape::ladderize(as(tree_pruned, "phylo"), right = T), show.tip.label=F, x.lim=179.90, edge.width=1, edge.color=branch_cols) ape::nodelabels(c("Theropoda", "Tetanurae", "Coelurosauria", "Avialae", "Ornithuromorpha"), c(103,114,123,158,169), frame = "none", adj=c(1.05,1.2), cex=0.5) heights <- c(86, 68.5, 60, 54, 22, 12) lengths <- c(1, 6, 8, 10, 19, 25) invisible(apply(cbind(heights, lengths), 1, function (x) { lines(rep(abs(ages[x[2]] - max(ages)), 2), c(0, x[1]), lty=2, col="dimgray") })) par(mai=c(1.2,0.8,0,1.6), mgp=c(0,0.6,0)) palette = RColorBrewer::brewer.pal(3, "Set2") plot(anc ~ ages, xaxp=c(260,70,19), axes=F, type="n", xlab="", ylab="", xlim=rev(range(ages)), col=palette[1]) rect(185, min(anc), 145, max(anc), border=NA, col=rgb(0.9, 0.9, 0.9)) lines(ages, anc, col=palette[1], lwd=0.8) points(ages, anc, col=palette[1], pch=16, cex=0.6) axis(2, ylim=range(anc), line=0.5, col=palette[1], las=1, cex.axis=0.5, lwd=0.8) mtext(text="Log body mass (kg)", side=2, line=2.0, cex=0.6) par(new=T) plot(pcoa_bm$pcoa1 ~ pcoa_bm$ages, axes=F, type="n", xlab="", ylab="", xlim=rev(range(pcoa_bm$ages)), col=palette[2]) lines(pcoa_bm$ages, pcoa_bm$pcoa1, col=palette[2], lwd=0.8) points(pcoa_bm$ages, pcoa_bm$pcoa1, col=palette[2], pch=16, cex=0.6) axis(4, ylim=range(pcoa_bm$pcoa1), line=0.5, col=palette[2], las=1, cex.axis=0.5, lwd=0.8) mtext(text="PCo1", side=4, line=2.0, cex=0.6) par(new=T) plot(pcoa_bm_gls$residuals ~ pcoa_bm$ages, axes=F, type="n", xlab="", ylab="", xlim=rev(range(pcoa_bm$ages)), ylim=range(c(pcoa_bm_gls$residuals)), pch=21, col=palette[3]) lines(pcoa_bm$ages, pcoa_bm_gls$residuals, col=palette[3], lwd=0.8) points(pcoa_bm$ages, pcoa_bm_gls$residuals, col=palette[3], pch=16, cex=0.6) axis(4, ylim=range(pcoa_bm_gls$residuals), line=3.5, col=palette[3], las=1, cex.axis=0.5, lwd=0.8) mtext(text="PCo1 vs. log body mass (kg) residual", side=4, line=5, cex=0.6) axis(1, pretty(range(ages), 13), line = 0.6, cex.axis=0.5, lwd=0.8) mtext("Time (millions of years ago)", side=1, col="black", line=2, cex=0.6) mtext(paste("Phase I"), side=1, line=-0.5, at=217.5, cex=0.7) mtext(paste("Phase II"), side=1, line=-0.5, at=165, cex=0.7) mtext(paste("Phase III"), side=1, line=-0.5, at=107.58, cex=0.7) ```