--- title: "Analysis code for Friedman et al., Evolution of a Multifunctional Trait" subtitle: "Shared effects of foraging ecology and thermoregulation on beak morphology, with consequences for song evolution" author: - Nicholas R. Friedman - Eliot T. Miller - Jason R. Ball - Haruka Kasuga - Vladimír Remeš - Evan P. Economo date: "`r format(Sys.time(), '%d %B %Y')`" abstract: While morphological traits are often associated with multiple functions, it remains unclear how evolution balances the selective effects of different functions. Birds' beaks function in foraging, but also in thermoregulating and singing, among other behaviours. Studies of beak evolution abound, however most focus on a single function. Thus, we quantified relative contributions of different functions over an evolutionary time scale. We measured beak shape using geometric morphometrics and compared this trait to foraging behaviour, climatic variables, and song characteristics in a phylogenetic comparative study of an Australasian radiation of songbirds (Meliphagidae). We found that both climate and foraging behaviour were significantly correlated with beak shape and size. However, foraging ecology had a greater effect on shape, and climate had a greater effect on size. We also found that evolutionary changes in beak morphology had significant consequences for vocal performance; species with elongate-shaped beaks sang at higher frequencies, while species with large beaks sang at a slower pace. The evolution of the avian beak exemplifies how morphological traits can represent an evolutionary compromise among functions, and suggest that specialization along any functional axis may increase ecological divergence or reproductive isolation along others. output: html_document: toc: true toc_depth: 2 toc_float: true --- # Introduction This R Markdown notebook describes analyses comparing honeyeater beak shape and size to their various functions: feeding, thermoregulation, and singing. The goal of this document is to describe the methods behind this study, and reproduce its key findings in total. The analyses in this notebook are not presented in any particular order with respect to the paper, but are organized instead to preserve the flow of various R objects defined above to code below. In particular: - The first section of the document describes the dataset, especially the loadings of Principal component axes. - The second section of the document describes a set of of multivariate PGLS analyses conducted to test for correlations between foraging behavior, themoregulation, and beak morphology. They also test for correlations between beak morphology and song characteristics. - The third section describes a set of phylogenetically-corrected path analyses conducted to test whether beak characteristics are map to individual functions, or whether they are shared by multiple functions. This analysis is repeated using several different approaches to dealing with allometry in a comparative framework. - The fourth section describes a set of analyses intended to estimate how much evolutionary variation in beak morphology each function can explain. - The fifth section describes a two-block Partial Least Squares analysis aimed at estimating the overall correlation between two highly dimensional sets of traits: foraging behavior and beak morphology. - The sixth section compares estimates of beak size using centroid size to more conventional measurements made using calipers. - The seventh section explores the effect of temperature on beak size and shape using polynomial regression terms, as in Greenberg and Danner (2012). - The eighth section repeats phylogenetic regressions of beak size and song characteristics while accounting for intraspecific variation in song characteristics. Here are all the libraries that were used in these analyses: ```{r, message=FALSE, warning=FALSE} library(geomorph, quietly = T) library(extrafont,quietly = T) library(phytools, quietly = T) library(geiger, quietly = T) library(nlme,quietly = T) library(ggplot2,quietly = T) library(phylopath,quietly = T) library(wesanderson, quietly = T) #fun but not really necessary! library(sensiPhy, quietly = T) ``` # Data First I'm going to load the data. These datasets are organized to include all species relevant to the paper. ```{r} load("Friedman_Data.RData") ``` ## Description of the Dataset Features of this dataset, described in more detail in the paper, include the following: - *avgRostrumLatP* - data on beak shape, which is organized in an array of matrices, each one describing the x and y coordinates landmarks for a specimen I digitized in TPSdig. These data were aligned using a Generalized Procrustes analysis and what we're looking at are the species averages. So we can procede to visualize these in _geomorph_. - *Measurements* - linear measurements taken from museum specimens at Tring. Table includes the first principal component of these measurements, as well as the centroid size of each species' beak landmarks for comparison. - *FeedShape* - A data frame for running comparative analyses, including information on feeding behavior PCs 1-4, beak shape PCs 1-4, average temperature for each species' range, body mass from the literature, and beak size both as centroid size and residuals of its correlation with body size. - *SongShape* - A data frame for running comparative analyses, including information on beak shape PCs 1-4, beak size and its residuals from body size, body size from the literature, and three song characteristics of each species: note minimum frequency, note maximum frequency, and overall song pace. - *Tree* - A dated phylogenetic tree from Marki et al. 2017 MPE. _This was not produced by our data but is required for our analysis; please be courteous and contact Petter Marki before using this tree in your project_. - *MeliSong_AllP* - A data frame including song data from individuals, for use in Ives & Garland style analyses that account for intraspecific variation. ## PC1 - Figure 1 ```{r} plotTangentSpace(avgRostrumLatP) ``` PC1 explains 86% of beak variation. It describes variation between narrow, elongate beaks on the left and short, deep beaks on the right. We refer to this as "depth" in the paper; "elongation" also works but isn't in line with the polarity of the axis as shown here. ## PC2 ```{r, warning=FALSE, message=FALSE} plotTangentSpace(avgRostrumLatP, axis1 = 2, axis2 = 1) ``` The morphospace here makes the same figure, rotated 90º to the left. What we're most interested in here are the warp grids explaining the variation on this axis. Beaks vary on PC2 between curved on the left, to straight on the right. This was confusing to explain in the manuscript as "straightness", so for many of the results the sign is flipped to give "curvature" instead. This axis explains 10% of variation. ## PC3 - Figure S3 ```{r, warning=FALSE, message=FALSE} plotTangentSpace(avgRostrumLatP, axis1 = 3, axis2 = 1) ``` PC3 explains only 2% of variation, and is a subtle difference: the beaks on the left begin to taper distally along the rostrum, while the beaks on the right taper more proximally, coming to a more acute angle at the tip. ## Foraging pPCs (Figure S1) We ran a phylogenetically corrected PCA on Eliot Miller and Sarah Wagner's amazing field data. ```{r} FeedPC$L[,1:4] ``` These plots describe the loadings of those pPC axes. Here's pPCs 1 and 2. ```{r} biplot(FeedPC,cex = 0.4,xlim = c(-40,40),family = "serif") ``` And here's pPCs 3 and 4. ```{r} biplot(FeedPC,cex = 0.4,xlim = c(-40,40),family = "serif",choices= 3:4) ``` # Multivariate PGLS With all the datasets loaded, what's next is to do the analyses. For each analysis we're going to do PGLS using predictor and response variables scaled by their standard deviation. That way the output model will have the slope and its confidence interval equivalent to standardized effect size, which we can then compare using forest plots. ## Feeding behavior and climate as predictors of beak evolution First we need to prune the phylogeny to include only the relevant taxa. ```{r} drops <- name.check(phy = Tree,data = FeedShape)$tree_not_data TreeFP <- drop.tip(phy = Tree,tip = drops) ``` For these analyses, I will be using body size as the allometry covariate. This data comes from the literature (Wilman et al. 2014). Doing our first analysis, which compares the relative role of beak shape (PC1) to feeding behavior and winter/summer temperatures ```{r} BothTestBeak1 <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 + BodyMass.Value, data = FeedShape, correlation = corPagel(value = 1,phy = TreeFP)) summary(BothTestBeak1) ``` Notice how I'm scaling by SD and centering values below and throughout, this is to give standardized Beta regression coefficients that are useful as effect size. Significance tests are unaffected. ```{r} BothTestBeak1 <- gls(model = scale(BeakPC1) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) + scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value), data = FeedShape, correlation = corPagel(value = 1, phy = TreeFP)) summary(BothTestBeak1) ``` There's a significant effect of climate on beak depth, but a much greater effect of feeding PC1 Let's do another for beak curvature: ```{r} BothTestBeak2 <- gls(model = scale(BeakPC2) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) + scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value),data = FeedShape,correlation = corPagel(value = 1,phy = TreeFP)) summary(BothTestBeak2) ``` Again it looks like there's a really strong effect of feeding behavior and diet on beak shape - curvature this time. Let's do this again for beak tapering. ```{r} BothTestBeak3 <- gls(model = scale(BeakPC3) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) + scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value),data = FeedShape,correlation = corPagel(value = 1,phy = TreeFP)) summary(BothTestBeak3) ``` Again there are significant relationships between feeding and our third axis of beak shape, as well as a significant positive correlation with winter minimum temperatures. Let's do this one more time for beak size overall: ```{r} BothTestBeak4 <- gls(model = scale(Csize) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) + scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value),data = FeedShape,correlation = corPagel(value = 1,phy = TreeFP)) summary(BothTestBeak4) ``` Obviously there is a very tight correlation between body size and beak size. Even so beak size is correlated with both winter temperature and feeding behavior PC 1 (nectarivory, etc). There's a significant relationship between beak size and body size (obviously), but also a very strong effect of winter minimum temperatures on beak sizes. There is a marginally weaker effect of feeding behavior and diet on beak size. ## Beak Shape and Size as predictors of Song evolution As before, we have to make a pruned tree for the analysis first. ```{r} drops <- name.check(phy = Tree,data = SongShape)$tree_not_data TreeSP <- drop.tip(phy = Tree,tip = drops) ``` Note that body size was used for the frequency correlations while beak size was used for the song pace correlation. The relationship between body size and frequency is well-documented, and likewise we have a priori reasons to believe that a bigger beak should be more difficult to move rapidly. Using both together creates a statistical issue, because they are highly collinear. This is solved by using only the relevant one, here and in previous studies of beak evolution (e.g., Grant & Grant 2002). ```{r} BeakFreq <- gls(scale(NoteMinFreq) ~ scale(BodyMass.Value) + scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3), data = SongShape,correlation = corPagel(1,TreeSP)) summary(BeakFreq) ``` Looks like there's a marginally significant relationship with beak shape PC1. Birds with long narrow beaks tend to sing higher songs. ```{r} BeakMaxFreq <- gls(scale(NoteMaxFreq) ~ scale(BodyMass.Value) + scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3), data = SongShape,correlation = corPagel(1,TreeSP)) summary(BeakMaxFreq) ``` Looks like there's a significant relationship with beak shape PC1. Birds with long narrow beaks tend to sing higher songs. ```{r} BeakPace <- gls(scale(OverallSongPace) ~ scale(Csize) + scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3), data = SongShape,correlation = corPagel(0.5,TreeSP)) #tweaked init lambda to avoid false convergence summary(BeakPace) ``` This is so cool! Let's try this one again, but with beak size as body size residuals. ```{r} BeakPace2 <- gls(scale(OverallSongPace) ~ scale(Resid.Beak.Size) + scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3), data = SongShape,correlation = corPagel(0.5,TreeSP)) #tweaked init lambda to avoid false convergence summary(BeakPace2) ``` Gee what do you know, the results are the same. Let's try this again again, throwing beak size and body size into a models of their own to resolve the collinearity. ```{r} BeakPace3a <- gls(scale(OverallSongPace) ~ scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3) + scale(BodyMass.Value) , data = SongShape,correlation = corPagel(0.5,TreeSP)) #tweaked init lambda to avoid false convergence summary(BeakPace3a) ``` ```{r} BeakPace3b <- gls(scale(OverallSongPace) ~ scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3) + scale(Csize) , data = SongShape,correlation = corPagel(0.5,TreeSP)) #tweaked init lambda to avoid false convergence summary(BeakPace3b) ``` ### Forest plot of song effect sizes This code organizes the output of the code above into data frames that behave well in ggplot. Apologies if this is ugly. ```{r} BeakPaceOut <- data.frame(Effect = BeakPace$coefficients,Error = summary(BeakPace)$tTable[,2],Upper = BeakPace$coefficients + summary(BeakPace)$tTable[,2]*1.96, Lower = BeakPace$coefficients - summary(BeakPace)$tTable[,2]*1.96) BeakPaceOut <- BeakPaceOut[-1,] BeakPaceOut <- cbind(BeakPaceOut,Response = "Pace", Predictor = c("Beak Size","Elongation","Curvature","Taper")) BeakMaxFreqOut <- data.frame(Effect = BeakMaxFreq$coefficients,Error = summary(BeakMaxFreq)$tTable[,2],Upper = BeakMaxFreq$coefficients + summary(BeakMaxFreq)$tTable[,2]*1.96, Lower = BeakMaxFreq$coefficients - summary(BeakMaxFreq)$tTable[,2]*1.96) BeakMaxFreqOut <- BeakMaxFreqOut[-1,] BeakMaxFreqOut <- cbind(BeakMaxFreqOut,Response = "MaxFreq", Predictor = c("Body Size","Elongation","Curvature","Taper")) BeakFreqOut <- data.frame(Effect = BeakFreq$coefficients,Error = summary(BeakFreq)$tTable[,2],Upper = BeakFreq$coefficients + summary(BeakFreq)$tTable[,2]*1.96, Lower = BeakFreq$coefficients - summary(BeakFreq)$tTable[,2]*1.96) BeakFreqOut <- BeakFreqOut[-1,] BeakFreqOut <- cbind(BeakFreqOut,Response = "MinFreq", Predictor = c("Body Size","Elongation","Curvature","Taper")) ``` ## Figure 4 Here's the code for the forest plot. *Note* that the direction of the morphology PC axes may be reversed from the figure - this was to match textual descriptions of "curvature", as well as "elongation" vs "depth". ```{r,} SongOutput <- rbind(BeakPaceOut,BeakMaxFreqOut,BeakFreqOut) ggplot(data = SongOutput,aes(x = Predictor,y = Effect))+ facet_grid( ~ Response)+ geom_pointrange(aes(ymin = Lower, ymax = Upper),fatten = 0.5)+ coord_flip()+ theme_classic()+ theme(text = element_text(size = 9,family = "serif"))+ geom_hline(yintercept = 0,linetype = "dashed")+ scale_y_continuous(limits = c(-0.75,0.75)) ``` And here's a second one using the residuals. ```{r} BeakPaceOut <- data.frame(Effect = BeakPace2$coefficients,Error = summary(BeakPace2)$tTable[,2],Upper = BeakPace2$coefficients + summary(BeakPace2)$tTable[,2]*1.96, Lower = BeakPace2$coefficients - summary(BeakPace2)$tTable[,2]*1.96) BeakPaceOut <- BeakPaceOut[-1,] BeakPaceOut <- cbind(BeakPaceOut,Response = "Pace", Predictor = c("Beak Size","Elongation","Curvature","Taper")) BeakMaxFreqOut <- data.frame(Effect = BeakMaxFreq$coefficients,Error = summary(BeakMaxFreq)$tTable[,2],Upper = BeakMaxFreq$coefficients + summary(BeakMaxFreq)$tTable[,2]*1.96, Lower = BeakMaxFreq$coefficients - summary(BeakMaxFreq)$tTable[,2]*1.96) BeakMaxFreqOut <- BeakMaxFreqOut[-1,] BeakMaxFreqOut <- cbind(BeakMaxFreqOut,Response = "MaxFreq", Predictor = c("Body Size","Elongation","Curvature","Taper")) BeakFreqOut <- data.frame(Effect = BeakFreq$coefficients,Error = summary(BeakFreq)$tTable[,2],Upper = BeakFreq$coefficients + summary(BeakFreq)$tTable[,2]*1.96, Lower = BeakFreq$coefficients - summary(BeakFreq)$tTable[,2]*1.96) BeakFreqOut <- BeakFreqOut[-1,] BeakFreqOut <- cbind(BeakFreqOut,Response = "MinFreq", Predictor = c("Body Size","Elongation","Curvature","Taper")) ``` ```{r,} SongOutput <- rbind(BeakPaceOut,BeakMaxFreqOut,BeakFreqOut) ggplot(data = SongOutput,aes(x = Predictor,y = Effect))+ facet_grid( ~ Response)+ geom_pointrange(aes(ymin = Lower, ymax = Upper),fatten = 0.5)+ coord_flip()+ theme_classic()+ theme(text = element_text(size = 9,family = "serif"))+ geom_hline(yintercept = 0,linetype = "dashed")+ scale_y_continuous(limits = c(-0.88,0.88)) ``` And here's a third one using body size instead of beak size as correction for pace. ```{r} BeakPaceOut2 <- data.frame(Effect = BeakPace3a$coefficients,Error = summary(BeakPace3a)$tTable[,2],Upper = BeakPace3a$coefficients + summary(BeakPace3a)$tTable[,2]*1.96, Lower = BeakPace3a$coefficients - summary(BeakPace3a)$tTable[,2]*1.96) BeakPaceOut2 <- BeakPaceOut2[-1,] BeakPaceOut2 <- cbind(BeakPaceOut2,Response = "Pace", Predictor = c("Depth","Curvature","Taper","Body Size")) ``` ## Figure S7 ```{r,} ggplot(data = BeakPaceOut2,aes(x = Predictor,y = Effect))+ geom_pointrange(aes(ymin = Lower, ymax = Upper),fatten = 0.5)+ coord_flip()+ theme_classic()+ theme(text = element_text(size = 9,family = "serif"))+ geom_hline(yintercept = 0,linetype = "dashed")+ scale_y_continuous(limits = c(-0.88,0.88)) ``` # Path Analysis So how are all these different traits related? Many of them ought to be correlated in ways we're interested, and some in ways we aren't. Instead of plugging everything into one big model we need to define an explicit hypothesis. Here I'm testing whether function is split (size vs. shape) or shared (size and shape). This ties into a set of alternative hypotheses defined in the paper: are different characteristics of the beak specialized for different functions, or are there trade-offs between the functions? #### Note: I'm switching the direction of two axes here: Beak PC2 (so it becomes curvature) and Feed PC1 (so it becomes nectarivory). This makes things easier to discuss in the text. ```{r} FeedPath <- FeedShape colnames(FeedPath) <- c("Species","Elong","Curve","Taper","Nect","Gape","Lerp","FeedPC4","Size","SmaxT","WminT","Body", "Beak.R") FeedPath$Curve <- -FeedPath$Curve FeedPath$Nect <- -FeedPath$Nect ``` I've decided to organize the path diagram into some distinct model sets as a hypothesis test. These vary in the degree to which functions are shared between trait axes of variation, or to which they are split between these. They also vary in the degree to which functional explanations for trait evolution "interfere", which is to say that functional predictor variables have correlations with each other that create indirect effects. - *Split*: Functions are split between different aspects of trait variation, based on a priori expectations - *Shared*: Functions are shared between different axes of trait variation - *Intersplit*: The same as *Split*, but with interference between predictor variables - *Intershare*: The same as *Shared*, but with interference between predictor variables Here is a slight re-organization. We've seen that diet and climate are both factors. But are they shared or split? And are they interfering with one another? ```{r} models <- define_model_set( Split = c(Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Size ~ SmaxT), Intersplit = c(Nect ~ SmaxT,Nect ~ WminT,Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Size ~ SmaxT), Shared = c(Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Elong ~ WminT, Size ~ Nect, Size ~ SmaxT, Curve ~ SmaxT, Elong ~ SmaxT), Intershare = c(Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Elong ~ WminT, Size ~ Nect, Nect ~ WminT, Size ~ SmaxT, Elong ~ SmaxT, Curve ~ SmaxT, Nect ~ SmaxT)) plot_model_set(models,algorithm = "circle") ``` We are going to do a path analysis here using a Pagel's lambda model. ```{r,warning=FALSE,message=FALSE,results='hide'} result <- phylo_path(model_set = models,data = FeedPath,tree = TreeFP) ``` ```{r} summary(result) ``` When we plot this as CICc weights, it's pretty strongly supportive of the model that has both functions shared and interfering. ```{r} CICcplot0<-plot(summary(result)) CICcplot0+theme_classic() ``` ## No Allometric Correction This is what the best model looks like. ```{r} plot(average(result,cut_off = 2),curvature = 0.05, arrow = grid::arrow(type = "open"),algorithm = "circle") ``` *Note* - one some computers the red arrow between Nect and Elong doesn't print! The variables may also be in different corners from figure to figure. Please not also that Elongation is not a proper description of the variable, and it has been re-named to the proper "Depth" in Figure 2. ### Allometry Tweak Here I'm repeating the path analysis above using a tweak, which is to look at the body-size residuals of beak size instead of just the raw beak size. This should correct for changes that are explained by body size, but without adding heaps of new parameters to the model. ```{r} othermodels <- define_model_set( Split = c(Elong ~ Nect, Curve ~ Nect, Beak.R ~ WminT, Beak.R ~ SmaxT), Intersplit = c(Nect ~ SmaxT,Nect ~ WminT,Elong ~ Nect, Curve ~ Nect, Beak.R ~ WminT, Beak.R ~ SmaxT), Shared = c(Elong ~ Nect, Curve ~ Nect, Beak.R ~ WminT, Elong ~ WminT, Beak.R ~ Nect, Beak.R ~ SmaxT, Curve ~ SmaxT, Elong ~ SmaxT), Intershare = c(Elong ~ Nect, Curve ~ Nect, Beak.R ~ WminT, Elong ~ WminT, Beak.R ~ Nect, Nect ~ WminT, Beak.R ~ SmaxT, Elong ~ SmaxT, Curve ~ SmaxT, Nect ~ SmaxT)) plot_model_set(othermodels,algorithm = "circle") ``` So we'll run the path analysis here. ```{r,warning=FALSE,message=FALSE,results='hide'} otherresult <- phylo_path(model_set = othermodels,data = FeedPath,tree = TreeFP) ``` ```{r} summary(otherresult) ``` The model selection result is the same. ```{r} CICcplot <- plot(summary(otherresult)) CICcplot+theme_classic() ``` ## Beak Size Residuals (Figure 2) This is what the best model looks like. ```{r} plot(average(otherresult,cut_off = 2),curvature = 0.05, arrow = grid::arrow(type = "open"),algorithm = "circle", text_size = 3) ``` Here we are reporting the indirect effects in the model. ## Figure S4 ```{r,warning=FALSE} x <- coef_plot(average(otherresult,cut_off = 2)) y <- x+ theme_classic()+ coord_flip()+ theme(axis.text.y = element_text(size = 9)) y ``` It's about the same, except that the relationships between beak size (now residuals) are stronger for winter temperature and nectarivory. ### Another Allometry Tweak - Including Body Size as a Covariate Here I'm repeating the path analysis above using a tweak, which is to look at the body-size residuals of beak size instead of just the raw beak size. This should correct for changes that are explained by body size, but without adding heaps of new parameters to the model. ```{r} othermodels2 <- define_model_set( Split = c(Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Size ~ SmaxT), Intersplit = c(Nect ~ SmaxT,Nect ~ WminT,Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Size ~ SmaxT), Shared = c(Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Elong ~ WminT, Size ~ Nect, Size ~ SmaxT, Curve ~ SmaxT, Elong ~ SmaxT), Intershare = c(Elong ~ Nect, Curve ~ Nect, Size ~ WminT, Elong ~ WminT, Size ~ Nect, Nect ~ WminT, Size ~ SmaxT, Elong ~ SmaxT, Curve ~ SmaxT, Nect ~ SmaxT),.common = c(Size ~ Body)) plot_model_set(othermodels2,algorithm = "circle") ``` So we'll run the path analysis here. ```{r,warning=FALSE,message=FALSE,results='hide'} otherresult2 <- phylo_path(model_set = othermodels2,data = FeedPath,tree = TreeFP) ``` ```{r} summary(otherresult2) ``` The model selection result is the same. ```{r} CICcplot2 <- plot(summary(otherresult2)) CICcplot2+theme_classic() ``` ## Body Size as Covariate This is what the best model looks like. ```{r} plot(average(otherresult2,cut_off = 2),curvature = 0.05, arrow = grid::arrow(type = "open"),algorithm = "circle", text_size = 3) ``` # Multivariate Contribution to Variation We want to be able to measure how much evolutionary variation each function of the beak is able to explain. So I'm going to do two analyses: 1) One analysis for beak shape as a multivariate response variable, using the methods described by Dean Adams and his research group and implemented in the R package _geomorph_. 2) Another analysis for beak size as a univariate response variable. In this case, I will be using _phylopath_ and building a simplified model. The goal of each of these is to get an R-squared value for each function that describes its effect on the overall trait. ## Beak Shape So first lets save this as a geomorph.data.frame. ```{r} RostrumP <- avgRostrumLatP[,,dimnames(avgRostrumLatP)[[3]] %in% FeedShape$Species] Rostrum.gdf <- geomorph.data.frame(Rostrum = RostrumP, FeedPC1 = FeedShape$FeedPC1, FeedPC2 = FeedShape$FeedPC2, FeedPC3 = FeedShape$FeedPC3, FeedPC4 = FeedShape$FeedPC4, Summer = FeedShape$bio5, Winter = FeedShape$bio6, Beak.Size = FeedShape$Csize, Body.Size = FeedShape$BodyMass.Value) ``` ```{r, message=FALSE, warning=FALSE, results='hide'} Rostrum.out <- procD.pgls(Rostrum ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + Winter + Summer + Body.Size + Beak.Size, data = Rostrum.gdf, phy = TreeFP, iter = 5000) ``` Here are the results ```{r} summary(Rostrum.out) ``` Okay, so I've also included allometry effects with body mass and beak size, and they each explain about 1% of the variation. ```{r} sum(Rostrum.out$aov.table$Rsq[1:8], na.rm = T) ``` Re-doing this without song as it's a response of the response in one of our models. ```{r} cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") R2s <- Rostrum.out$R2 names(R2s) <- c("FeedPC1","FeedPC2","FeedPC3","FeedPC4","Summer","Winter","Body.M","Beak.S") R2s <- data.frame(Predictor = names(R2s),R2 = R2s,Category = c("Foraging","Foraging","Foraging","Foraging","Thermoregulation","Thermoregulation","Allometry","Allometry")) shapebar <- ggplot(data = R2s,mapping = aes(x = Predictor,y = R2, fill = Category))+ geom_col(width = 0.85)+ scale_fill_manual(values=cbbPalette)+ scale_y_continuous(limits = c(0,0.25))+ ggtitle("Beak Shape")+ theme_minimal() plot(shapebar) ``` So foraging PC1 explains most of the variation here. ## Beak Size I'm using a phylogenetic path model that is re-organized into just a normal multivariate regression. This is to produce comparable coefficients of determination. ```{r} models <- define_model_set( All = c(Csize ~ FeedPC1, Csize ~ FeedPC2, Csize ~ FeedPC3, Csize ~ FeedPC4, Csize ~ bio5, Csize ~ bio6, Csize ~ BodyMass.Value)) ``` Here's the analysis. ```{r,warning=FALSE,message=FALSE,results='hide'} resultz <- phylo_path(model_set = models,data = FeedShape,tree = TreeFP) resultz.rz <- best(resultz)$coef[,"Csize"] ``` And here I'm just organizing it and plotting it. ```{r} require(ggplot2,quietly = T) names(resultz.rz) <- c("Body.M","Winter","Summer","FeedPC4","FeedPC3","FeedPC2","FeedPC1", "Csize") resultz.rz <- resultz.rz[-length(resultz.rz)] resultz.rz <- resultz.rz[sort(names(resultz.rz))] R2z <- data.frame(Predictor = names(resultz.rz),R2 = resultz.rz^2,Category = c("Allometry","Foraging","Foraging","Foraging","Foraging","Thermoregulation","Thermoregulation")) sizebar <- ggplot(data = R2z,mapping = aes(x = Predictor,y = R2, fill = Category))+ geom_col()+ scale_fill_manual(values=cbbPalette)+ ggtitle("Beak Size")+ theme_minimal() plot(sizebar) ``` ```{r} sum(R2z$R2) ``` Let's plot them together ## Figure 3 ```{r} R2z <- cbind(R2z,Type = "Size") R2s <- cbind(R2s,Type = "Shape") allR2s <- rbind(R2z,R2s) positionz <- c("Body.M","Beak.S","Winter","Summer","FeedPC4","FeedPC3","FeedPC2","FeedPC1") lifeaquatic <- wes_palette("Zissou1")[c(1,3,5)] ggplot(data = allR2s,mapping = aes(x = Predictor,y = R2, fill = Category))+ facet_grid(Type~.)+ geom_col()+ scale_fill_manual(values=lifeaquatic)+ scale_x_discrete(limits = positionz)+ coord_flip()+ theme_minimal() ``` How much variation explained by each model? ```{r} c(Size = sum(R2z$R2), Shape = sum(R2s$R2)) ``` I'm still only explaining half of the variation. That's pretty low, but it makes sense to me in light of Jen Bright's results. Much of the variation should be explained by integration with the skull. # Two-block PLS This analysis was kindly suggested by one of our reviewers as a way to test patterns of covariation between two multivariate datasets, one describing feeding behavior and the other describing beak shape. In the future someone should do this with song as well, but to include it in this paper would probably explode our methods section. First I have to arrange the Foraging PCs with the landmark data, make sure we're pruned to a matching dataset. ```{r} FeedPCP <- FeedPC$S[rownames(FeedPC$S) %in% dimnames(avgRostrumLatP)[[3]],] avgRostrumLatPP <- avgRostrumLatP[,,dimnames(avgRostrumLatP)[[3]] %in% rownames(FeedPC$S)] ``` This code will run a two-block PLS analysis. ```{r} ShapeXForag <- two.b.pls(A1 = avgRostrumLatPP, A2 = FeedPCP, iter = 5000,print.progress = F) ShapeXForag ``` Looks like a decent correlation between these variables, but in this case we didn't correct for phylogeny. ```{r} Tree2bP <- treedata(phy = Tree, data = FeedPCP, warnings = F)$phy ShapeXForagP <- phylo.integration(avgRostrumLatPP, FeedPCP, iter = 5000,print.progress = F, phy = Tree2bP) ShapeXForagP ``` # Size Measurement Comparisons Here in response to one of the reviewers' comments, I'm comparing the caliper-measured beak size versus the centroid-size-derived beak size. It's probably worth showing that they're similar, though mathematically the centroid-size should be a less biased approach. Just the same, since we didn't include some portions of the beak this may be a necessary confirmation. I'm pulling in a file of measurements I made at the museum. ```{r} ggplot(data = Measurements, mapping = aes(x = Csize, y = PC1))+ geom_point()+ geom_smooth(method = "lm")+ ylab("PC1 of Beak Measurements")+ xlab("Centroid Size of Landmarks")+ theme_classic() ``` Is this a tight correlation? Yes, but how tight? ```{r} summary(lm(data = Measurements, PC1 ~ Csize)) ``` The R-squared value is 0.874. It could be higher, but there are a few species off from the trendline. # Polynomial Temperature PGLS ## Figure S5 One of the reviewers asked if we could test a curvilinear relationship as predicting beak size and shape. This should follow Russ Greenberg and Ray Danner's work on Sparrow bills. Apologies from the author are due because the code here is functional but quite ugly. Let's explore the data first to get an idea which model might be most appropriate. ```{r} loess.gg <- ggplot(data = FeedShape, mapping = aes(x = bio6, y = Resid.Beak.Size))+ geom_point()+ geom_smooth(method = "loess")+ theme_classic()+ ggtitle("Loess Smoothed Beak Size")+ xlab("Winter minimum temperature ºC*10")+ theme_classic()+ ylab("Beak Size (Body Size Residuals)") loess.gg ``` ## Beak Shape PC1 Here I am fitting polynomial regression models. To interpret the results, I will be using a model test at the end. ```{r} #optimization issue - setting fixed lambda beakPC1 <- FeedShape$BeakPC1 names(beakPC1) <- rownames(FeedShape) beakPC1.lam <- fitContinuous(TreeFP,dat = beakPC1, model = "lambda") StraightBeak <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC1.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak2w <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6,2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC1.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak2s <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 2) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC1.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak2sw <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,2) + poly(bio6,2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC1.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak3w <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6,3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC1.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak3s <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 3) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC1.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak3sw <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,3) + poly(bio6,3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC1.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") ``` Is it a better fit than the linear model? Let's make some more and compare their fit. ```{r} anovout <- anova(StraightBeak, CurvyBeak2s, CurvyBeak2w, CurvyBeak2sw, CurvyBeak3s, CurvyBeak3w, CurvyBeak3sw) dAIC <- anovout$AIC - min(anovout$AIC) names(dAIC) <- rownames(anovout) sort(dAIC) ``` The formula with a third order polynomial for winter temperatures was the best fit to the data. Here's what the best model looks like. ```{r} summary(CurvyBeak3w) ``` Let's plot that bivariate relationship using a 3rd order polynomial, as kindly suggested by our reviewer 3. ```{r} beak1win.gg <- ggplot(data = FeedShape, mapping = aes(x = bio6, y = BeakPC1))+ geom_point()+ stat_smooth(method = "lm", se = T, formula=y ~ poly(x, 3, raw=TRUE))+ ggtitle("Polynomial Regression Model for Beak Shape PC1")+ xlab("Winter minimum temperature ºC*10")+ theme_classic() beak1win.gg ``` ## Beak Shape PC2 Doing this again for Beak PC2. ```{r,warning=FALSE} #optimization issue - setting fixed lambda beakPC2 <- FeedShape$BeakPC2 names(beakPC2) <- rownames(FeedShape) beakPC2.lam <- fitContinuous(TreeFP,dat = beakPC2, model = "lambda") StraightBeak <- gls(model = BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC2.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak2w <- gls(model = BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6, 2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC2.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak2s <- gls(model = BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 2) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC2.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak2sw <- gls(model = BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,3) + poly(bio6,2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC2.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak3w <- gls(model = BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6, 3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC2.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak3s <- gls(model = BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 3) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC2.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") CurvyBeak3sw <- gls(model = BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,3) + poly(bio6,3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC2.lam$opt$lambda,fixed = T,phy = TreeFP), method = "ML") anovout <- anova(StraightBeak, CurvyBeak2s, CurvyBeak2w, CurvyBeak2sw, CurvyBeak3s, CurvyBeak3w, CurvyBeak3sw) dAIC <- anovout$AIC - min(anovout$AIC) names(dAIC) <- rownames(anovout) sort(dAIC) ``` Here the best model has doesn't have any polynomial fits for winter minimum temperatures or summer maximum temperatures. ```{r} summary(StraightBeak) ``` It's a bit ugly, but presumably the rest of our multivariate regression model ## Beak Shape PC3 Doing this again for Beak PC3. ```{r,warning=FALSE} beakPC3 <- FeedShape$BeakPC3 names(beakPC3) <- rownames(FeedShape) beakPC3.lam <- fitContinuous(TreeFP,dat = beakPC3, model = "lambda") StraightBeak <- gls(model = BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC3.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak2w <- gls(model = BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6, 2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC3.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak2s <- gls(model = BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 2) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC3.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak2sw <- gls(model = BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,3) + poly(bio6,2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC3.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak3w <- gls(model = BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6, 3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC3.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak3s <- gls(model = BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 3) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC3.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak3sw <- gls(model = BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,3) + poly(bio6,3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = beakPC3.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") anovout <- anova(StraightBeak, CurvyBeak2s, CurvyBeak2w, CurvyBeak2sw, CurvyBeak3s, CurvyBeak3w, CurvyBeak3sw) dAIC <- anovout$AIC - min(anovout$AIC) names(dAIC) <- rownames(anovout) sort(dAIC) ``` Here's the best model. It's the one with 3rd order polynomials for both summer and winter temperatures. ```{r} summary(CurvyBeak3sw) ``` Let's plot those bivariate regression models. ```{r} beak3win.gg <- ggplot(data = FeedShape, mapping = aes(x = bio6, y = BeakPC3))+ geom_point()+ stat_smooth(method = "lm", se = T, formula=y ~ poly(x, 3, raw=TRUE))+ ggtitle("Polynomial Regression Model for Beak Shape PC3")+ xlab("Winter minimum temperature ºC*10")+ theme_classic() beak3win.gg ``` ```{r} beak3sum.gg <- ggplot(data = FeedShape, mapping = aes(x = bio5, y = BeakPC3))+ geom_point()+ stat_smooth(method = "lm", se = T, formula=y ~ poly(x, 3, raw=TRUE))+ ggtitle("Polynomial Regression Model for Beak Shape PC3")+ xlab("Summer maximum temperature ºC*10")+ theme_classic() beak3sum.gg ``` Now for beak size. ```{r,warning=FALSE} csize <- FeedShape$Csize names(csize) <- rownames(FeedShape) csize.lam <- fitContinuous(TreeFP,dat = csize, model = "lambda") StraightBeak <- gls(model = Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = csize.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak2w <- gls(model = Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6, 2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = csize.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak2s <- gls(model = Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 2) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = csize.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak2sw <- gls(model = Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,3) + poly(bio6,2) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = csize.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak3w <- gls(model = Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6, 3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = csize.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak3s <- gls(model = Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 3) + bio6 + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = csize.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") CurvyBeak3sw <- gls(model = Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,3) + poly(bio6,3) + BodyMass.Value,data = FeedShape ,correlation = corPagel(value = csize.lam$opt$lambda, fixed = T,phy = TreeFP), method = "ML") anovout <- anova(StraightBeak, CurvyBeak2s, CurvyBeak2w, CurvyBeak2sw, CurvyBeak3s, CurvyBeak3w, CurvyBeak3sw) dAIC <- anovout$AIC - min(anovout$AIC) names(dAIC) <- rownames(anovout) sort(dAIC) ``` Here's the best model: ```{r} summary(CurvyBeak3sw) ``` As above, it has 3rd order polynomials for both summer and winter temperatures. Let's try to visualize what that looks like. ```{r} sizewin.gg <- ggplot(data = FeedShape, mapping = aes(x = bio6, y = Csize))+ geom_point()+ stat_smooth(method = "lm", se = T, formula=y ~ poly(x, 3, raw=TRUE))+ ggtitle("Polynomial Regression Model for Beak Size")+ xlab("Winter minimum temperature ºC*10")+ ylab("Beak Size (landmark centroid size)")+ theme_classic() sizewin.gg ``` ```{r} sizesum.gg <- ggplot(data = FeedShape, mapping = aes(x = bio5, y = Csize))+ geom_point()+ stat_smooth(method = "lm", se = T, formula=y ~ poly(x, 3, raw=TRUE))+ ggtitle("Polynomial Regression Model for Beak Size")+ xlab("Summer maximum temperature ºC*10")+ ylab("Beak Size (landmark centroid size)")+ theme_classic() sizesum.gg ``` # Song PGLS w/ Intraspecific Variation ```{r} load("Processed_Data.RData") ``` Instead of using raw values, which cause ML optimization errors, first I'm going to first estimate some variances for my response variables (song characteristics). ```{r} MeliSong_AllPP <- MeliSong_AllP[MeliSong_AllP$species %in% SongShape$Species,] SongVar <- data.frame( NoteMaxFreq = tapply(X = MeliSong_AllPP$NoteMaxFreq,INDEX = MeliSong_AllPP$species, FUN = mean), NoteMinFreq = tapply(X = MeliSong_AllPP$NoteMinFreq,INDEX = MeliSong_AllPP$species, FUN = mean), OverallSongPace = tapply(X = MeliSong_AllPP$OverallSongPace,INDEX = MeliSong_AllPP$species, FUN = mean), NoteMaxFreq_var = tapply(X = MeliSong_AllPP$NoteMaxFreq,INDEX = MeliSong_AllPP$species, FUN = sd), NoteMinFreq_var = tapply(X = MeliSong_AllPP$NoteMinFreq,INDEX = MeliSong_AllPP$species, FUN = sd), OverallSongPace_var = tapply(X = MeliSong_AllPP$OverallSongPace,INDEX = MeliSong_AllPP$species, FUN = sd) ) ``` Putting this together. ```{r} SongShape_var <- merge(SongShape, SongVar[,4:6], by.x = "Species", by.y = "row.names") ``` Doing the analysis. ### Song Pace ~ Beak Elongation ```{r, message=FALSE, warning=FALSE, results='hide'} PaceXpc1 <- intra_phylm(formula = OverallSongPace ~ BeakPC1, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda") ``` ```{r} summary(PaceXpc1) ``` ### Song Pace ~ Beak Size ```{r, message=FALSE, warning=FALSE, results='hide'} PaceXsize <- intra_phylm(formula = OverallSongPace ~ Csize, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda") ``` ```{r} summary(PaceXsize) ``` ### Song Pace ~ Beak Curvature ```{r, message=FALSE, warning=FALSE, results='hide'} PaceXpc2 <- intra_phylm(formula = OverallSongPace ~ BeakPC2, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda") ``` ```{r} summary(PaceXpc2) ``` ### Song Pace ~ Beak Tapering ```{r, message=FALSE, warning=FALSE, results='hide'} PaceXpc3 <- intra_phylm(formula = OverallSongPace ~ BeakPC3, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda") ``` ```{r} summary(PaceXpc3) ``` ### Max Frequency ~ Beak Elongation ```{r, message=FALSE, warning=FALSE, results='hide'} MaxFrXpc1 <- intra_phylm(formula = NoteMaxFreq ~ BeakPC1, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda") ``` ```{r} summary(MaxFrXpc1) ``` ### Max Frequency ~ Beak Curvature ```{r, message=FALSE, warning=FALSE, results='hide'} MaxFrXpc2 <- intra_phylm(formula = NoteMaxFreq ~ BeakPC2, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda") ``` ```{r} summary(MaxFrXpc2) ``` ### Max Frequency ~ Beak Tapering ```{r, message=FALSE, warning=FALSE, results='hide'} MaxFrXpc3 <- intra_phylm(formula = NoteMaxFreq ~ BeakPC3, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda") ``` ```{r} summary(MaxFrXpc3) ``` ### Max Frequency ~ Beak Size ```{r, message=FALSE, warning=FALSE, results='hide'} MaxFrXsize <- intra_phylm(formula = NoteMaxFreq ~ Csize, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda") ``` ```{r} summary(MaxFrXsize) ``` ### Max Frequency ~ Body Size ```{r, message=FALSE, warning=FALSE, results='hide'} MaxFrXbody <- intra_phylm(formula = NoteMaxFreq ~ BodyMass.Value, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda") ``` ```{r} summary(MaxFrXbody) ``` ### Min Frequency ~ Beak Elongation ```{r, message=FALSE, warning=FALSE, results='hide'} MinFrXpc1 <- intra_phylm(formula = NoteMinFreq ~ BeakPC1, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda") ``` ```{r} summary(MinFrXpc1) ``` ### Min Frequency ~ Beak Curvature ```{r, message=FALSE, warning=FALSE, results='hide'} MinFrXpc2 <- intra_phylm(formula = NoteMinFreq ~ BeakPC2, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda") ``` ```{r} summary(MinFrXpc2) ``` ### Min Frequency ~ Beak Tapering ```{r, message=FALSE, warning=FALSE, results='hide'} MinFrXpc3 <- intra_phylm(formula = NoteMinFreq ~ BeakPC3, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda") ``` ```{r} summary(MinFrXpc3) ``` ### Min Frequency ~ Beak Size ```{r, message=FALSE, warning=FALSE, results='hide'} MinFrXsize <- intra_phylm(formula = NoteMinFreq ~ Csize, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda") ``` ```{r} summary(MinFrXsize) ``` ### Min Frequency ~ Body Size ```{r, message=FALSE, warning=FALSE, results='hide'} MinFrXbody <- intra_phylm(formula = NoteMinFreq ~ BodyMass.Value, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda") ``` ```{r} summary(MinFrXbody) ``` It's significant. But it's difficult to distinguish between this and beak size because they are correalted with each other AND we are restricted to a bivariate regression. But we understand the proximate mechanism by which body size influences minimum frequency, and indeed minimum frequency is an honest index of body size in many organisms because of this. We can try to tease them apart a little bit by doing a regression with the body size residuals of beak size. Interpret with caution because this approach is not recommended for reasons detailed in Freckleton (2009). ### Min Frequency ~ Beak Size Residuals ```{r, message=FALSE, warning=FALSE, results='hide'} MinFrXsize.r <- intra_phylm(formula = NoteMinFreq ~ Resid.Beak.Size, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda") ``` ```{r} summary(MinFrXsize.r) ``` That's all! Thanks for reading.