Analysis code for Friedman et al., Evolution of a Multifunctional Trait

Shared effects of foraging ecology and thermoregulation on beak morphology, with consequences for song evolution

Nicholas R. Friedman

Eliot T. Miller

Jason R. Ball

Haruka Kasuga

Vladimír Remeš

Evan P. Economo

19 November 2019

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.

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:

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.

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

plotTangentSpace(avgRostrumLatP)

## 
## PC Summary
## 
## Importance of first k=25 (out of 40) components:
##                            PC1     PC2      PC3      PC4      PC5      PC6
## Standard deviation     0.04756 0.01584 0.007158 0.005358 0.004284 0.002943
## Proportion of Variance 0.85679 0.09506 0.019410 0.010870 0.006950 0.003280
## Cumulative Proportion  0.85679 0.95185 0.971260 0.982130 0.989080 0.992360
##                             PC7      PC8      PC9     PC10     PC11
## Standard deviation     0.002485 0.002042 0.001625 0.001495 0.001095
## Proportion of Variance 0.002340 0.001580 0.001000 0.000850 0.000450
## Cumulative Proportion  0.994700 0.996280 0.997280 0.998130 0.998580
##                             PC12      PC13      PC14      PC15      PC16
## Standard deviation     0.0008419 0.0007961 0.0006216 0.0006129 0.0005465
## Proportion of Variance 0.0002700 0.0002400 0.0001500 0.0001400 0.0001100
## Cumulative Proportion  0.9988500 0.9990900 0.9992300 0.9993800 0.9994900
##                             PC17      PC18      PC19    PC20      PC21
## Standard deviation     0.0004886 0.0004436 0.0004167 0.00039 0.0003254
## Proportion of Variance 0.0000900 0.0000700 0.0000700 0.00006 0.0000400
## Cumulative Proportion  0.9995800 0.9996500 0.9997200 0.99978 0.9998200
##                             PC22      PC23      PC24      PC25
## Standard deviation     0.0003126 0.0002689 0.0002549 0.0002463
## Proportion of Variance 0.0000400 0.0000300 0.0000200 0.0000200
## Cumulative Proportion  0.9998600 0.9998800 0.9999100 0.9999300

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

plotTangentSpace(avgRostrumLatP, axis1 = 2, axis2 = 1)

## 
## PC Summary
## 
## Importance of first k=25 (out of 40) components:
##                            PC1     PC2      PC3      PC4      PC5      PC6
## Standard deviation     0.04756 0.01584 0.007158 0.005358 0.004284 0.002943
## Proportion of Variance 0.85679 0.09506 0.019410 0.010870 0.006950 0.003280
## Cumulative Proportion  0.85679 0.95185 0.971260 0.982130 0.989080 0.992360
##                             PC7      PC8      PC9     PC10     PC11
## Standard deviation     0.002485 0.002042 0.001625 0.001495 0.001095
## Proportion of Variance 0.002340 0.001580 0.001000 0.000850 0.000450
## Cumulative Proportion  0.994700 0.996280 0.997280 0.998130 0.998580
##                             PC12      PC13      PC14      PC15      PC16
## Standard deviation     0.0008419 0.0007961 0.0006216 0.0006129 0.0005465
## Proportion of Variance 0.0002700 0.0002400 0.0001500 0.0001400 0.0001100
## Cumulative Proportion  0.9988500 0.9990900 0.9992300 0.9993800 0.9994900
##                             PC17      PC18      PC19    PC20      PC21
## Standard deviation     0.0004886 0.0004436 0.0004167 0.00039 0.0003254
## Proportion of Variance 0.0000900 0.0000700 0.0000700 0.00006 0.0000400
## Cumulative Proportion  0.9995800 0.9996500 0.9997200 0.99978 0.9998200
##                             PC22      PC23      PC24      PC25
## Standard deviation     0.0003126 0.0002689 0.0002549 0.0002463
## Proportion of Variance 0.0000400 0.0000300 0.0000200 0.0000200
## Cumulative Proportion  0.9998600 0.9998800 0.9999100 0.9999300

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

plotTangentSpace(avgRostrumLatP, axis1 = 3, axis2 = 1)

## 
## PC Summary
## 
## Importance of first k=25 (out of 40) components:
##                            PC1     PC2      PC3      PC4      PC5      PC6
## Standard deviation     0.04756 0.01584 0.007158 0.005358 0.004284 0.002943
## Proportion of Variance 0.85679 0.09506 0.019410 0.010870 0.006950 0.003280
## Cumulative Proportion  0.85679 0.95185 0.971260 0.982130 0.989080 0.992360
##                             PC7      PC8      PC9     PC10     PC11
## Standard deviation     0.002485 0.002042 0.001625 0.001495 0.001095
## Proportion of Variance 0.002340 0.001580 0.001000 0.000850 0.000450
## Cumulative Proportion  0.994700 0.996280 0.997280 0.998130 0.998580
##                             PC12      PC13      PC14      PC15      PC16
## Standard deviation     0.0008419 0.0007961 0.0006216 0.0006129 0.0005465
## Proportion of Variance 0.0002700 0.0002400 0.0001500 0.0001400 0.0001100
## Cumulative Proportion  0.9988500 0.9990900 0.9992300 0.9993800 0.9994900
##                             PC17      PC18      PC19    PC20      PC21
## Standard deviation     0.0004886 0.0004436 0.0004167 0.00039 0.0003254
## Proportion of Variance 0.0000900 0.0000700 0.0000700 0.00006 0.0000400
## Cumulative Proportion  0.9995800 0.9996500 0.9997200 0.99978 0.9998200
##                             PC22      PC23      PC24      PC25
## Standard deviation     0.0003126 0.0002689 0.0002549 0.0002463
## Proportion of Variance 0.0000400 0.0000300 0.0000200 0.0000200
## Cumulative Proportion  0.9998600 0.9998800 0.9999100 0.9999300

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.

FeedPC$L[,1:4]

##                              PC1          PC2          PC3         PC4
## attacks             -0.300868812  0.067473759  0.033046688 -0.20855591
## mean.attack.height   0.273038441  0.458570202 -0.087417029  0.45707504
## mean.percent.canopy -0.216519336  0.714122845  0.235280417  0.11284242
## mean.FD              0.064830556  0.632498181  0.602377329 -0.20220236
## mean.third          -0.269509735 -0.699007255  0.302709086 -0.05388987
## hanging              0.486892774  0.407237183 -0.115717577 -0.04562436
## reaching             0.258278573  0.166015704  0.111585652  0.61823702
## flutter.chasing      0.110685596  0.164816806 -0.036044241  0.39373685
## gaping               0.208746242  0.164365589 -0.736647635 -0.12895108
## gleaning             0.860919306 -0.303537831  0.245981826  0.09350913
## lunging              0.033840633 -0.468942390 -0.060846819  0.09130288
## pecking              0.160863976 -0.622440464 -0.107695747 -0.11860946
## probing             -0.712901820  0.457721972 -0.178032379 -0.23379714
## pulling              0.310809953  0.188791906 -0.645153282 -0.46676840
## sally.gliding       -0.008049029  0.039090115 -0.004873902  0.13363164
## sally.hovering      -0.119516334  0.017868041  0.178996814  0.18405238
## sally.pouncing      -0.001774025 -0.140242486 -0.023011614  0.14447588
## sally.stalling      -0.243979619 -0.101438513 -0.073602168  0.16406321
## sally.striking      -0.723443520 -0.116905551 -0.107692583  0.19109560
## air                 -0.731777632 -0.194176038 -0.158337100  0.28111547
## branches             0.469473076  0.219108578 -0.581592349  0.28631443
## nectarivory         -0.850264287  0.299178537  0.089820071 -0.17770500
## frugivory            0.160609477 -0.027734960  0.229840010 -0.10550339
## ground               0.135447750 -0.929394748 -0.120594279 -0.05996870
## hanging.bark         0.284028729  0.137768003 -0.752354655  0.10866576
## insect.cases         0.366296708  0.127327655  0.112581659 -0.68812984
## leaves               0.669696648  0.163047663  0.446710234  0.03582521
## webs                 0.098889957  0.001466701  0.198964781  0.10233841
## woody.fruits         0.260621692  0.071265592  0.086118698 -0.64879044
## dead                 0.355068176  0.213605418  0.196123587  0.20540258
## mean.flower.length  -0.379308788  0.224955252 -0.247099452 -0.06681988

These plots describe the loadings of those pPC axes. Here’s pPCs 1 and 2.

biplot(FeedPC,cex = 0.4,xlim = c(-40,40),family = "serif")

And here’s pPCs 3 and 4.

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.

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

BothTestBeak1 <- gls(model = BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 + BodyMass.Value, data = FeedShape, correlation = corPagel(value = 1,phy = TreeFP))
summary(BothTestBeak1)

## Generalized least squares fit by REML
##   Model: BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 +      BodyMass.Value 
##   Data: FeedShape 
##         AIC       BIC   logLik
##   -170.3389 -149.0676 95.16947
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.9009919 
## 
## Coefficients:
##                      Value  Std.Error   t-value p-value
## (Intercept)    -0.05770286 0.03455196 -1.670032  0.1000
## FeedPC1         0.00166127 0.00034717  4.785222  0.0000
## FeedPC2        -0.00014885 0.00039231 -0.379418  0.7057
## FeedPC3         0.00063785 0.00035217  1.811175  0.0750
## FeedPC4         0.00057118 0.00043978  1.298775  0.1988
## bio5            0.00019785 0.00010775  1.836241  0.0711
## bio6           -0.00027489 0.00008246 -3.333457  0.0015
## BodyMass.Value  0.00019154 0.00016800  1.140113  0.2586
## 
##  Correlation: 
##                (Intr) FedPC1 FedPC2 FedPC3 FedPC4 bio5   bio6  
## FeedPC1        -0.351                                          
## FeedPC2        -0.101  0.105                                   
## FeedPC3         0.355 -0.094  0.034                            
## FeedPC4        -0.385  0.055  0.030 -0.163                     
## bio5           -0.913  0.464  0.148 -0.341  0.398              
## bio6            0.195 -0.396 -0.232 -0.111  0.022 -0.420       
## BodyMass.Value -0.180 -0.153 -0.001  0.060 -0.147  0.057 -0.016
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.12014590  0.09583374  0.55581680  1.17402987  2.58807907 
## 
## Residual standard error: 0.03961267 
## Degrees of freedom: 70 total; 62 residual

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.

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)

## Generalized least squares fit by REML
##   Model: scale(BeakPC1) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) +      scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value) 
##   Data: FeedShape 
##        AIC      BIC    logLik
##   162.0753 183.3467 -71.03766
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.9009919 
## 
## Coefficients:
##                            Value  Std.Error   t-value p-value
## (Intercept)           -0.4288978 0.23279291 -1.842401  0.0702
## scale(FeedPC1)         0.4850455 0.10136324  4.785222  0.0000
## scale(FeedPC2)        -0.0457750 0.12064517 -0.379418  0.7057
## scale(FeedPC3)         0.1329754 0.07341941  1.811175  0.0750
## scale(FeedPC4)         0.1086868 0.08368412  1.298775  0.1988
## scale(bio5)            0.1687930 0.09192312  1.836241  0.0711
## scale(bio6)           -0.3083274 0.09249477 -3.333457  0.0015
## scale(BodyMass.Value)  0.1101640 0.09662547  1.140113  0.2586
## 
##  Correlation: 
##                       (Intr) s(FPC1 s(FPC2 s(FPC3 s(FPC4 scl(5) scl(6)
## scale(FeedPC1)         0.029                                          
## scale(FeedPC2)        -0.002  0.105                                   
## scale(FeedPC3)        -0.016 -0.094  0.034                            
## scale(FeedPC4)         0.014  0.055  0.030 -0.163                     
## scale(bio5)            0.035  0.464  0.148 -0.341  0.398              
## scale(bio6)           -0.055 -0.396 -0.232 -0.111  0.022 -0.420       
## scale(BodyMass.Value)  0.057 -0.153 -0.001  0.060 -0.147  0.057 -0.016
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.12014590  0.09583374  0.55581680  1.17402987  2.58807907 
## 
## Residual standard error: 0.8110893 
## Degrees of freedom: 70 total; 62 residual

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:

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)

## Generalized least squares fit by REML
##   Model: scale(BeakPC2) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) +      scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value) 
##   Data: FeedShape 
##        AIC      BIC    logLik
##   176.6441 197.9155 -78.32206
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.5755824 
## 
## Coefficients:
##                            Value  Std.Error   t-value p-value
## (Intercept)           -0.0634508 0.18672044 -0.339817  0.7351
## scale(FeedPC1)         0.4077383 0.10587463  3.851143  0.0003
## scale(FeedPC2)        -0.1142615 0.12435060 -0.918865  0.3617
## scale(FeedPC3)        -0.1219028 0.08877174 -1.373216  0.1746
## scale(FeedPC4)        -0.3821025 0.09282679 -4.116295  0.0001
## scale(bio5)           -0.2638018 0.11317308 -2.330959  0.0230
## scale(bio6)           -0.0242679 0.10511918 -0.230861  0.8182
## scale(BodyMass.Value)  0.0428782 0.10242129  0.418646  0.6769
## 
##  Correlation: 
##                       (Intr) s(FPC1 s(FPC2 s(FPC3 s(FPC4 scl(5) scl(6)
## scale(FeedPC1)         0.038                                          
## scale(FeedPC2)         0.010  0.140                                   
## scale(FeedPC3)        -0.015 -0.080 -0.014                            
## scale(FeedPC4)         0.027  0.104  0.076 -0.186                     
## scale(bio5)            0.044  0.405  0.228 -0.352  0.414              
## scale(bio6)           -0.070 -0.361 -0.292 -0.122 -0.070 -0.421       
## scale(BodyMass.Value)  0.065 -0.099 -0.063  0.054 -0.142  0.043  0.026
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.24486839 -0.45007843  0.09352754  0.53625662  2.36218861 
## 
## Residual standard error: 0.7573884 
## Degrees of freedom: 70 total; 62 residual

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.

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)

## Generalized least squares fit by REML
##   Model: scale(BeakPC3) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) +      scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value) 
##   Data: FeedShape 
##        AIC     BIC    logLik
##   192.8096 214.081 -86.40481
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.6825852 
## 
## Coefficients:
##                            Value  Std.Error    t-value p-value
## (Intercept)            0.0460549 0.23584047  0.1952801  0.8458
## scale(FeedPC1)         0.1484494 0.12411849  1.1960298  0.2362
## scale(FeedPC2)         0.4421165 0.14641112  3.0196917  0.0037
## scale(FeedPC3)         0.2384635 0.09933904  2.4005018  0.0194
## scale(FeedPC4)        -0.2061582 0.10640612 -1.9374658  0.0572
## scale(bio5)            0.0589917 0.12701230  0.4644564  0.6439
## scale(bio6)           -0.2651540 0.11992940 -2.2109170  0.0307
## scale(BodyMass.Value) -0.1100857 0.11955667 -0.9207827  0.3607
## 
##  Correlation: 
##                       (Intr) s(FPC1 s(FPC2 s(FPC3 s(FPC4 scl(5) scl(6)
## scale(FeedPC1)         0.036                                          
## scale(FeedPC2)         0.005  0.132                                   
## scale(FeedPC3)        -0.016 -0.085  0.001                            
## scale(FeedPC4)         0.023  0.093  0.057 -0.182                     
## scale(bio5)            0.041  0.422  0.200 -0.354  0.410              
## scale(bio6)           -0.066 -0.371 -0.275 -0.119 -0.045 -0.418       
## scale(BodyMass.Value)  0.064 -0.112 -0.044  0.058 -0.143  0.044  0.012
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -2.260042758 -0.568551620 -0.007110102  0.448293541  2.043195867 
## 
## Residual standard error: 0.9038668 
## Degrees of freedom: 70 total; 62 residual

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:

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)

## Generalized least squares fit by REML
##   Model: scale(Csize) ~ scale(FeedPC1) + scale(FeedPC2) + scale(FeedPC3) +      scale(FeedPC4) + scale(bio5) + scale(bio6) + scale(BodyMass.Value) 
##   Data: FeedShape 
##      AIC      BIC    logLik
##   130.28 151.5514 -55.14002
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.9690184 
## 
## Coefficients:
##                            Value  Std.Error   t-value p-value
## (Intercept)            0.2152897 0.19819639  1.086244  0.2816
## scale(FeedPC1)        -0.3144137 0.07785852 -4.038270  0.0002
## scale(FeedPC2)        -0.0326254 0.09453025 -0.345132  0.7312
## scale(FeedPC3)        -0.0803681 0.05531618 -1.452886  0.1513
## scale(FeedPC4)        -0.0447979 0.06435243 -0.696133  0.4889
## scale(bio5)           -0.0726711 0.06672642 -1.089091  0.2803
## scale(bio6)            0.2887158 0.07045598  4.097819  0.0001
## scale(BodyMass.Value)  0.5666160 0.07483075  7.571967  0.0000
## 
##  Correlation: 
##                       (Intr) s(FPC1 s(FPC2 s(FPC3 s(FPC4 scl(5) scl(6)
## scale(FeedPC1)         0.025                                          
## scale(FeedPC2)        -0.003  0.085                                   
## scale(FeedPC3)        -0.016 -0.105  0.046                            
## scale(FeedPC4)         0.011  0.025  0.034 -0.149                     
## scale(bio5)            0.033  0.486  0.139 -0.324  0.394              
## scale(bio6)           -0.050 -0.404 -0.213 -0.102  0.052 -0.430       
## scale(BodyMass.Value)  0.052 -0.176  0.014  0.059 -0.148  0.063 -0.022
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.89486485 -0.94824863 -0.44697310  0.02109088  3.21306282 
## 
## Residual standard error: 0.6732761 
## Degrees of freedom: 70 total; 62 residual

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.

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).

BeakFreq <- gls(scale(NoteMinFreq) ~ scale(BodyMass.Value) + scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3), data = SongShape,correlation = corPagel(1,TreeSP))
summary(BeakFreq)

## Generalized least squares fit by REML
##   Model: scale(NoteMinFreq) ~ scale(BodyMass.Value) + scale(BeakPC1) +      scale(BeakPC2) + scale(BeakPC3) 
##   Data: SongShape 
##        AIC      BIC    logLik
##   244.1189 261.0507 -115.0594
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.3411759 
## 
## Coefficients:
##                            Value Std.Error    t-value p-value
## (Intercept)           -0.0064275 0.1796048 -0.0357868  0.9715
## scale(BodyMass.Value) -0.3164395 0.1115645 -2.8363809  0.0057
## scale(BeakPC1)        -0.2304052 0.1144886 -2.0124728  0.0474
## scale(BeakPC2)         0.0755945 0.1033735  0.7312752  0.4667
## scale(BeakPC3)         0.1028731 0.1001175  1.0275235  0.3072
## 
##  Correlation: 
##                       (Intr) s(BM.V s(BPC1 s(BPC2
## scale(BodyMass.Value) -0.015                     
## scale(BeakPC1)         0.180 -0.235              
## scale(BeakPC2)         0.000 -0.028  0.022       
## scale(BeakPC3)         0.005  0.203  0.040  0.011
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6492688 -0.7379367 -0.1757267  0.5112410  3.3848071 
## 
## Residual standard error: 0.9013969 
## Degrees of freedom: 88 total; 83 residual

Looks like there’s a marginally significant relationship with beak shape PC1. Birds with long narrow beaks tend to sing higher songs.

BeakMaxFreq <- gls(scale(NoteMaxFreq) ~ scale(BodyMass.Value) + scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3), data = SongShape,correlation = corPagel(1,TreeSP))
summary(BeakMaxFreq)

## Generalized least squares fit by REML
##   Model: scale(NoteMaxFreq) ~ scale(BodyMass.Value) + scale(BeakPC1) +      scale(BeakPC2) + scale(BeakPC3) 
##   Data: SongShape 
##       AIC      BIC    logLik
##   228.841 245.7729 -107.4205
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.6738513 
## 
## Coefficients:
##                            Value  Std.Error    t-value p-value
## (Intercept)           -0.0211034 0.23086408 -0.0914106  0.9274
## scale(BodyMass.Value) -0.2049557 0.11080532 -1.8496920  0.0679
## scale(BeakPC1)        -0.3234271 0.11384653 -2.8409041  0.0057
## scale(BeakPC2)         0.0265772 0.09984132  0.2661943  0.7907
## scale(BeakPC3)         0.1506041 0.09460745  1.5918838  0.1152
## 
##  Correlation: 
##                       (Intr) s(BM.V s(BPC1 s(BPC2
## scale(BodyMass.Value) -0.005                     
## scale(BeakPC1)         0.172 -0.189              
## scale(BeakPC2)         0.020 -0.072  0.087       
## scale(BeakPC3)        -0.005  0.271 -0.014  0.050
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.81010818 -0.60396164 -0.06377713  0.56128702  3.03196395 
## 
## Residual standard error: 0.9218802 
## Degrees of freedom: 88 total; 83 residual

Looks like there’s a significant relationship with beak shape PC1. Birds with long narrow beaks tend to sing higher songs.

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)

## Generalized least squares fit by REML
##   Model: scale(OverallSongPace) ~ scale(Csize) + scale(BeakPC1) + scale(BeakPC2) +      scale(BeakPC3) 
##   Data: SongShape 
##        AIC      BIC    logLik
##   238.5469 255.4787 -112.2734
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##     lambda 
## -0.1662207 
## 
## Coefficients:
##                     Value  Std.Error   t-value p-value
## (Intercept)    -0.0109243 0.03256898 -0.335420  0.7382
## scale(Csize)   -0.4253327 0.07889242 -5.391300  0.0000
## scale(BeakPC1)  0.0934761 0.06517086  1.434324  0.1552
## scale(BeakPC2) -0.1210756 0.07133171 -1.697360  0.0934
## scale(BeakPC3) -0.1065205 0.08155491 -1.306120  0.1951
## 
##  Correlation: 
##                (Intr) scl(C) s(BPC1 s(BPC2
## scale(Csize)    0.033                     
## scale(BeakPC1) -0.756  0.251              
## scale(BeakPC2)  0.412  0.210  0.183       
## scale(BeakPC3)  0.082 -0.267  0.079 -0.188
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6801524 -0.6713135 -0.1204710  0.3218392  4.7179756 
## 
## Residual standard error: 0.7911833 
## Degrees of freedom: 88 total; 83 residual

This is so cool!

Let’s try this one again, but with beak size as body size residuals.

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)

## Generalized least squares fit by REML
##   Model: scale(OverallSongPace) ~ scale(Resid.Beak.Size) + scale(BeakPC1) +      scale(BeakPC2) + scale(BeakPC3) 
##   Data: SongShape 
##        AIC      BIC    logLik
##   259.5755 276.5074 -122.7877
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##       lambda 
## -0.005301327 
## 
## Coefficients:
##                             Value  Std.Error   t-value p-value
## (Intercept)             0.0032331 0.09684081  0.033386  0.9734
## scale(Resid.Beak.Size) -0.5476808 0.16852626 -3.249825  0.0017
## scale(BeakPC1)         -0.1810105 0.16127839 -1.122348  0.2650
## scale(BeakPC2)         -0.1771697 0.10664005 -1.661381  0.1004
## scale(BeakPC3)         -0.0088233 0.10237664 -0.086185  0.9315
## 
##  Correlation: 
##                        (Intr) s(R.B. s(BPC1 s(BPC2
## scale(Resid.Beak.Size)  0.001                     
## scale(BeakPC1)         -0.011  0.784              
## scale(BeakPC2)          0.005  0.351  0.259       
## scale(BeakPC3)          0.000 -0.183 -0.078 -0.107
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6156327 -0.5668975 -0.1936620  0.2685191  4.0574193 
## 
## Residual standard error: 0.9340229 
## Degrees of freedom: 88 total; 83 residual

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.

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)

## Generalized least squares fit by REML
##   Model: scale(OverallSongPace) ~ scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3) +      scale(BodyMass.Value) 
##   Data: SongShape 
##        AIC     BIC   logLik
##   254.8041 271.736 -120.402
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##      lambda 
## -0.06276459 
## 
## Coefficients:
##                            Value  Std.Error   t-value p-value
## (Intercept)            0.0219099 0.04569216  0.479511  0.6328
## scale(BeakPC1)         0.3776303 0.08573848  4.404444  0.0000
## scale(BeakPC2)        -0.0854093 0.08930729 -0.956353  0.3417
## scale(BeakPC3)        -0.0878343 0.09378595 -0.936540  0.3517
## scale(BodyMass.Value) -0.4284832 0.09562462 -4.480888  0.0000
## 
##  Correlation: 
##                       (Intr) s(BPC1 s(BPC2 s(BPC3
## scale(BeakPC1)        -0.830                     
## scale(BeakPC2)         0.274  0.053              
## scale(BeakPC3)         0.084  0.167 -0.080       
## scale(BodyMass.Value)  0.184 -0.271  0.045 -0.037
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6017351 -0.7295044 -0.1745410  0.3402999  3.9136970 
## 
## Residual standard error: 0.8952415 
## Degrees of freedom: 88 total; 83 residual

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)

## Generalized least squares fit by REML
##   Model: scale(OverallSongPace) ~ scale(BeakPC1) + scale(BeakPC2) + scale(BeakPC3) +      scale(Csize) 
##   Data: SongShape 
##        AIC      BIC    logLik
##   238.5469 255.4787 -112.2734
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##     lambda 
## -0.1662207 
## 
## Coefficients:
##                     Value  Std.Error   t-value p-value
## (Intercept)    -0.0109243 0.03256898 -0.335420  0.7382
## scale(BeakPC1)  0.0934761 0.06517086  1.434324  0.1552
## scale(BeakPC2) -0.1210756 0.07133171 -1.697360  0.0934
## scale(BeakPC3) -0.1065205 0.08155491 -1.306120  0.1951
## scale(Csize)   -0.4253327 0.07889242 -5.391300  0.0000
## 
##  Correlation: 
##                (Intr) s(BPC1 s(BPC2 s(BPC3
## scale(BeakPC1) -0.756                     
## scale(BeakPC2)  0.412  0.183              
## scale(BeakPC3)  0.082  0.079 -0.188       
## scale(Csize)    0.033  0.251  0.210 -0.267
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6801524 -0.6713135 -0.1204710  0.3218392  4.7179756 
## 
## Residual standard error: 0.7911833 
## Degrees of freedom: 88 total; 83 residual

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.

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”.

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.

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"))

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.

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

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.

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?

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.

result <- phylo_path(model_set = models,data = FeedPath,tree = TreeFP)

summary(result)

##        model  k  q      C     p   CICc delta_CICc     l     w
## 1 Intershare  5 16 26.366 0.003 68.630      0.000 1.000 0.992
## 2     Shared  7 14 43.849 0.000 79.485     10.855 0.004 0.004
## 3 Intersplit  9 12 50.122 0.000 79.596     10.965 0.004 0.004
## 4      Split 11 10 67.605 0.000 91.334     22.703 0.000 0.000

When we plot this as CICc weights, it’s pretty strongly supportive of the model that has both functions shared and interfering.

CICcplot0<-plot(summary(result))
CICcplot0+theme_classic()

No Allometric Correction

This is what the best model looks like.

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.

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.

otherresult <- phylo_path(model_set = othermodels,data = FeedPath,tree = TreeFP)

summary(otherresult)

##        model  k  q       C p    CICc delta_CICc     l     w
## 1 Intershare  5 16  50.811 0  93.075      0.000 1.000 0.996
## 2     Shared  7 14  68.294 0 103.930     10.855 0.004 0.004
## 3 Intersplit  9 12  87.592 0 117.066     23.991 0.000 0.000
## 4      Split 11 10 105.075 0 128.803     35.728 0.000 0.000

The model selection result is the same.

CICcplot <- plot(summary(otherresult))
CICcplot+theme_classic()

Beak Size Residuals (Figure 2)

This is what the best model looks like.

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

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

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.

otherresult2 <- phylo_path(model_set = othermodels2,data = FeedPath,tree = TreeFP)

summary(otherresult2)

##        model  k  q       C p    CICc delta_CICc     l     w
## 1 Intershare 10 18  48.057 0  97.469      0.000 1.000 0.999
## 2     Shared 12 16  69.114 0 111.378     13.909 0.001 0.001
## 3 Intersplit 14 14  84.531 0 120.167     22.699 0.000 0.000
## 4      Split 16 12 105.588 0 135.061     37.593 0.000 0.000

The model selection result is the same.

CICcplot2 <- plot(summary(otherresult2))
CICcplot2+theme_classic()

Body Size as Covariate

This is what the best model looks like.

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:

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.
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.

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)

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

summary(Rostrum.out)

## 
## Call:
## procD.pgls(f1 = Rostrum ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 +  
##     Winter + Summer + Body.Size + Beak.Size, phy = TreeFP, iter = 5000,  
##     data = Rostrum.gdf)
## 
## Type I (Sequential) Sums of Squares and Cross-products
## Randomized Residual Permutation Procedure Used
## 5001 Permutations
## ANOVA effect sizes and P-values based on empirical F distributions
## 
## 
##           Df        SS        MS     Rsq       F      Z    Pr(>F)    
## FeedPC1    1 0.0033250 0.0033250 0.21839 23.7014 3.5839 0.0002000 ***
## FeedPC2    1 0.0002124 0.0002124 0.01395  1.5141 1.1255 0.1421716    
## FeedPC3    1 0.0004933 0.0004933 0.03240  3.5166 1.8988 0.0231954 *  
## FeedPC4    1 0.0005128 0.0005128 0.03368  3.6551 2.0559 0.0167966 *  
## Winter     1 0.0010058 0.0010058 0.06606  7.1692 2.8832 0.0003999 ***
## Summer     1 0.0007984 0.0007984 0.05244  5.6910 2.7824 0.0005999 ***
## Body.Size  1 0.0000831 0.0000831 0.00546  0.5923 0.4894 0.3217357    
## Beak.Size  1 0.0002368 0.0002368 0.01555  1.6877 1.5579 0.0659868 .  
## Residuals 61 0.0085575 0.0001403 0.56207                             
## Total     69 0.0152250                                               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Okay, so I’ve also included allometry effects with body mass and beak size, and they each explain about 1% of the variation.

sum(Rostrum.out$aov.table$Rsq[1:8], na.rm = T)

## [1] 0.4379301

Re-doing this without song as it’s a response of the response in one of our models.

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.

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.

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.

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)

sum(R2z$R2)

## [1] 0.5128095

Let’s plot them together

Figure 3

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?

c(Size = sum(R2z$R2), Shape = sum(R2s$R2))

##      Size     Shape 
## 0.5128095 0.4379301

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.

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.

ShapeXForag <- two.b.pls(A1 = avgRostrumLatPP, A2 = FeedPCP, iter = 5000,print.progress = F)
ShapeXForag

## 
## Call:
## two.b.pls(A1 = avgRostrumLatPP, A2 = FeedPCP, iter = 5000, print.progress = F) 
## 
## 
## 
## 
## r-PLS: 0.632
## 
## P-value: 0
## 
## Based on 5001 random permutations

Looks like a decent correlation between these variables, but in this case we didn’t correct for phylogeny.

Tree2bP <- treedata(phy = Tree, data = FeedPCP, warnings = F)$phy
ShapeXForagP <- phylo.integration(avgRostrumLatPP, FeedPCP, iter = 5000,print.progress = F, phy = Tree2bP)
ShapeXForagP

## 
## Call:
## phylo.integration(A = avgRostrumLatPP, A2 = FeedPCP, phy = Tree2bP,  
##     iter = 5000, print.progress = F) 
## 
## 
## 
## r-PLS: 0.625
## 
## P-value: 0.003
## 
## Based on 5001 random permutations

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.

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?

summary(lm(data = Measurements, PC1 ~ Csize))

## 
## Call:
## lm(formula = PC1 ~ Csize, data = Measurements)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6977 -1.5051 -0.3248  1.1782  8.0219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -19.79313    0.76735  -25.79   <2e-16 ***
## Csize         1.19633    0.04391   27.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.58 on 107 degrees of freedom
## Multiple R-squared:  0.874,  Adjusted R-squared:  0.8728 
## F-statistic: 742.3 on 1 and 107 DF,  p-value: < 2.2e-16

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.

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.

#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.

anovout <- anova(StraightBeak, CurvyBeak2s, CurvyBeak2w, CurvyBeak2sw, CurvyBeak3s, CurvyBeak3w, CurvyBeak3sw)
dAIC <- anovout$AIC - min(anovout$AIC)
names(dAIC) <- rownames(anovout)
sort(dAIC)

##  CurvyBeak3w CurvyBeak3sw  CurvyBeak2w CurvyBeak2sw StraightBeak 
##     0.000000     2.523677     3.266648     5.046751     9.866792 
##  CurvyBeak2s  CurvyBeak3s 
##    10.692874    12.576101

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.

summary(CurvyBeak3w)

## Generalized least squares fit by maximum likelihood
##   Model: BeakPC1 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + poly(bio6,      3) + BodyMass.Value 
##   Data: FeedShape 
##         AIC       BIC   logLik
##   -295.4435 -270.7101 158.7218
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.8770267 
## 
## Coefficients:
##                      Value  Std.Error   t-value p-value
## (Intercept)    -0.19956025 0.04719567 -4.228359  0.0001
## FeedPC1         0.00149939 0.00032478  4.616695  0.0000
## FeedPC2        -0.00004810 0.00036809 -0.130662  0.8965
## FeedPC3         0.00065815 0.00032941  1.997942  0.0503
## FeedPC4         0.00059021 0.00040607  1.453486  0.1513
## bio5            0.00055522 0.00014138  3.927269  0.0002
## poly(bio6, 3)1 -0.18344785 0.03832921 -4.786111  0.0000
## poly(bio6, 3)2  0.12450136 0.03570155  3.487282  0.0009
## poly(bio6, 3)3 -0.06063942 0.02800536 -2.165279  0.0344
## BodyMass.Value  0.00035576 0.00016079  2.212560  0.0307
## 
##  Correlation: 
##                (Intr) FedPC1 FedPC2 FedPC3 FedPC4 bio5   p(6,3)1 p(6,3)2
## FeedPC1        -0.204                                                   
## FeedPC2        -0.157  0.124                                            
## FeedPC3         0.222 -0.108  0.006                                     
## FeedPC4        -0.275  0.050  0.021 -0.158                              
## bio5           -0.973  0.236  0.163 -0.238  0.299                       
## poly(bio6, 3)1  0.541 -0.304 -0.253 -0.099  0.006 -0.566                
## poly(bio6, 3)2 -0.671 -0.084  0.125 -0.020  0.006  0.682 -0.422         
## poly(bio6, 3)3  0.395  0.187  0.117 -0.104 -0.060 -0.398  0.201  -0.355 
## BodyMass.Value -0.288 -0.187 -0.012  0.076 -0.128  0.209 -0.108   0.206 
##                p(6,3)3
## FeedPC1               
## FeedPC2               
## FeedPC3               
## FeedPC4               
## bio5                  
## poly(bio6, 3)1        
## poly(bio6, 3)2        
## poly(bio6, 3)3        
## BodyMass.Value -0.248 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.5427043 -0.0886201  0.7408374  1.3582160  2.6058291 
## 
## Residual standard error: 0.03311475 
## Degrees of freedom: 70 total; 60 residual

Let’s plot that bivariate relationship using a 3rd order polynomial, as kindly suggested by our reviewer 3.

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.

#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)

## StraightBeak  CurvyBeak2w  CurvyBeak2s  CurvyBeak3w  CurvyBeak3s 
##     0.000000     1.205836     1.356685     2.675259     2.693288 
## CurvyBeak2sw CurvyBeak3sw 
##     4.266925     5.719561

Here the best model has doesn’t have any polynomial fits for winter minimum temperatures or summer maximum temperatures.

summary(StraightBeak)

## Generalized least squares fit by maximum likelihood
##   Model: BeakPC2 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + bio5 + bio6 +      BodyMass.Value 
##   Data: FeedShape 
##         AIC      BIC   logLik
##   -422.7604 -402.524 220.3802
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.7841637 
## 
## Coefficients:
##                      Value   Std.Error   t-value p-value
## (Intercept)     0.03347470 0.013467766  2.485542  0.0156
## FeedPC1         0.00041879 0.000127885  3.274740  0.0017
## FeedPC2        -0.00014727 0.000143602 -1.025537  0.3091
## FeedPC3        -0.00015085 0.000136872 -1.102122  0.2747
## FeedPC4        -0.00066067 0.000164972 -4.004731  0.0002
## bio5           -0.00011238 0.000042680 -2.633019  0.0107
## bio6            0.00000105 0.000031315  0.033542  0.9733
## BodyMass.Value  0.00004643 0.000062209  0.746293  0.4583
## 
##  Correlation: 
##                (Intr) FedPC1 FedPC2 FedPC3 FedPC4 bio5   bio6  
## FeedPC1        -0.343                                          
## FeedPC2        -0.122  0.122                                   
## FeedPC3         0.373 -0.088  0.015                            
## FeedPC4        -0.391  0.080  0.041 -0.177                     
## bio5           -0.931  0.439  0.174 -0.352  0.405              
## bio6            0.203 -0.381 -0.257 -0.116 -0.017 -0.416       
## BodyMass.Value -0.170 -0.127 -0.025  0.060 -0.145  0.048 -0.001
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.24791079 -0.45973542  0.08714408  0.62986989  2.39984863 
## 
## Residual standard error: 0.01287842 
## Degrees of freedom: 70 total; 62 residual

It’s a bit ugly, but presumably the rest of our multivariate regression model

Beak Shape PC3

Doing this again for Beak PC3.

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)

## CurvyBeak3sw CurvyBeak2sw  CurvyBeak2s  CurvyBeak3s  CurvyBeak2w 
##     0.000000     2.520109     3.769897     5.095679     5.805827 
##  CurvyBeak3w StraightBeak 
##     6.427798     6.439875

Here’s the best model. It’s the one with 3rd order polynomials for both summer and winter temperatures.

summary(CurvyBeak3sw)

## Generalized least squares fit by maximum likelihood
##   Model: BeakPC3 ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5,      3) + poly(bio6, 3) + BodyMass.Value 
##   Data: FeedShape 
##         AIC      BIC   logLik
##   -513.0944 -483.864 269.5472
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.7487043 
## 
## Coefficients:
##                       Value   Std.Error    t-value p-value
## (Intercept)     0.003113196 0.002036015  1.5290632  0.1317
## FeedPC1         0.000130191 0.000067747  1.9217404  0.0596
## FeedPC2         0.000197979 0.000074630  2.6527989  0.0103
## FeedPC3         0.000151145 0.000075214  2.0095138  0.0491
## FeedPC4        -0.000140658 0.000085006 -1.6546797  0.1034
## poly(bio5, 3)1  0.001936676 0.011065875  0.1750134  0.8617
## poly(bio5, 3)2 -0.015677179 0.005497360 -2.8517653  0.0060
## poly(bio5, 3)3 -0.006036437 0.005900126 -1.0231031  0.3105
## poly(bio6, 3)1 -0.017747844 0.008696175 -2.0408794  0.0458
## poly(bio6, 3)2 -0.010327528 0.008109615 -1.2734919  0.2079
## poly(bio6, 3)3 -0.012573150 0.006392283 -1.9669263  0.0540
## BodyMass.Value -0.000025583 0.000033044 -0.7742273  0.4419
## 
##  Correlation: 
##                (Intr) FedPC1 FedPC2 FedPC3 FedPC4 p(5,3)1 p(5,3)2 p(5,3)3
## FeedPC1         0.082                                                    
## FeedPC2         0.004  0.135                                             
## FeedPC3        -0.003 -0.154 -0.034                                      
## FeedPC4         0.059  0.086  0.036 -0.162                               
## poly(bio5, 3)1 -0.037  0.217  0.170 -0.218  0.320                        
## poly(bio5, 3)2  0.039 -0.224 -0.040  0.292 -0.038 -0.024                 
## poly(bio5, 3)3 -0.072 -0.034  0.053 -0.116 -0.103 -0.225  -0.023         
## poly(bio6, 3)1 -0.048 -0.273 -0.235 -0.126 -0.070 -0.596  -0.001   0.370 
## poly(bio6, 3)2 -0.048 -0.116  0.102  0.021  0.028  0.693   0.130  -0.184 
## poly(bio6, 3)3  0.031  0.095  0.110 -0.070 -0.075 -0.405   0.212   0.242 
## BodyMass.Value -0.421 -0.128 -0.022  0.019 -0.139  0.134  -0.113   0.142 
##                p(6,3)1 p(6,3)2 p(6,3)3
## FeedPC1                               
## FeedPC2                               
## FeedPC3                               
## FeedPC4                               
## poly(bio5, 3)1                        
## poly(bio5, 3)2                        
## poly(bio5, 3)3                        
## poly(bio6, 3)1                        
## poly(bio6, 3)2 -0.448                 
## poly(bio6, 3)3  0.269  -0.333         
## BodyMass.Value -0.018   0.121  -0.173 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4662326 -0.7222111 -0.1195357  0.3530150  2.3323859 
## 
## Residual standard error: 0.006253303 
## Degrees of freedom: 70 total; 58 residual

Let’s plot those bivariate regression models.

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

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.

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)

## CurvyBeak3sw  CurvyBeak3s CurvyBeak2sw  CurvyBeak2s  CurvyBeak3w 
##     0.000000     4.627534     6.549969     7.537153     9.856765 
## StraightBeak  CurvyBeak2w 
##    13.196291    14.752453

Here’s the best model:

summary(CurvyBeak3sw)

## Generalized least squares fit by maximum likelihood
##   Model: Csize ~ FeedPC1 + FeedPC2 + FeedPC3 + FeedPC4 + poly(bio5, 3) +      poly(bio6, 3) + BodyMass.Value 
##   Data: FeedShape 
##        AIC      BIC    logLik
##   304.0808 333.3113 -139.0404
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
## lambda 
##      1 
## 
## Coefficients:
##                    Value Std.Error   t-value p-value
## (Intercept)    12.885236  0.923004 13.960110  0.0000
## FeedPC1        -0.085135  0.022668 -3.755726  0.0004
## FeedPC2         0.009133  0.027284  0.334741  0.7390
## FeedPC3        -0.055782  0.022976 -2.427802  0.0183
## FeedPC4        -0.038429  0.028978 -1.326127  0.1900
## poly(bio5, 3)1 -8.644991  3.333445 -2.593410  0.0120
## poly(bio5, 3)2 -3.144423  1.375214 -2.286497  0.0259
## poly(bio5, 3)3  4.925861  1.734684  2.839629  0.0062
## poly(bio6, 3)1 14.909836  2.819624  5.287880  0.0000
## poly(bio6, 3)2 -2.963587  2.372469 -1.249157  0.2166
## poly(bio6, 3)3  5.487430  1.999058  2.745008  0.0080
## BodyMass.Value  0.091579  0.011989  7.638277  0.0000
## 
##  Correlation: 
##                (Intr) FedPC1 FedPC2 FedPC3 FedPC4 p(5,3)1 p(5,3)2 p(5,3)3
## FeedPC1         0.096                                                    
## FeedPC2        -0.021  0.097                                             
## FeedPC3        -0.016 -0.148  0.008                                      
## FeedPC4         0.040 -0.034  0.003 -0.106                               
## poly(bio5, 3)1 -0.044  0.261  0.118 -0.145  0.296                        
## poly(bio5, 3)2  0.054 -0.056 -0.090  0.222  0.076  0.133                 
## poly(bio5, 3)3 -0.038 -0.065  0.031 -0.039 -0.076 -0.327  -0.048         
## poly(bio6, 3)1 -0.016 -0.315 -0.192 -0.115  0.035 -0.626  -0.127   0.324 
## poly(bio6, 3)2 -0.060 -0.047  0.118  0.036  0.000  0.699   0.153  -0.201 
## poly(bio6, 3)3  0.069  0.162  0.103 -0.035 -0.099 -0.458   0.230   0.261 
## BodyMass.Value -0.330 -0.223  0.029  0.034 -0.129  0.193  -0.176   0.090 
##                p(6,3)1 p(6,3)2 p(6,3)3
## FeedPC1                               
## FeedPC2                               
## FeedPC3                               
## FeedPC4                               
## poly(bio5, 3)1                        
## poly(bio5, 3)2                        
## poly(bio5, 3)3                        
## poly(bio6, 3)1                        
## poly(bio6, 3)2 -0.462                 
## poly(bio6, 3)3  0.242  -0.364         
## BodyMass.Value -0.081   0.196  -0.301 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.95736343 -0.97143706 -0.57702433  0.02351743  3.29850711 
## 
## Residual standard error: 2.66751 
## Degrees of freedom: 70 total; 58 residual

As above, it has 3rd order polynomials for both summer and winter temperatures. Let’s try to visualize what that looks like.

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

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

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).

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.

SongShape_var <- merge(SongShape, SongVar[,4:6], by.x = "Species", by.y = "row.names")

Doing the analysis.

Song Pace ~ Beak Elongation

PaceXpc1 <- intra_phylm(formula = OverallSongPace ~ BeakPC1, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda")

summary(PaceXpc1)

##                  mean CI_low CI_high
## intercept       5.807  5.694   5.920
## se.intercept    0.457  0.429   0.484
## pval.intercept  0.000  0.000   0.000
## estimate       17.593 15.498  19.687
## se.estimate     9.565  9.218   9.911
## pval.estimate   0.115  0.066   0.165

Song Pace ~ Beak Size

PaceXsize <- intra_phylm(formula = OverallSongPace ~ Csize, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda")

summary(PaceXsize)

##                  mean CI_low CI_high
## intercept      10.316  9.947  10.685
## se.intercept    1.293  1.258   1.328
## pval.intercept  0.000  0.000   0.000
## estimate       -0.277 -0.294  -0.261
## se.estimate     0.075  0.073   0.077
## pval.estimate   0.001  0.001   0.002

Song Pace ~ Beak Curvature

PaceXpc2 <- intra_phylm(formula = OverallSongPace ~ BeakPC2, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda")

summary(PaceXpc2)

##                   mean  CI_low CI_high
## intercept        5.907   5.810   6.005
## se.intercept     0.465   0.436   0.493
## pval.intercept   0.000   0.000   0.000
## estimate       -12.792 -17.578  -8.007
## se.estimate     28.442  27.795  29.088
## pval.estimate    0.600   0.529   0.671

Song Pace ~ Beak Tapering

PaceXpc3 <- intra_phylm(formula = OverallSongPace ~ BeakPC3, data = SongShape_var,phy = TreeSP, "OverallSongPace_var",model = "lambda")

summary(PaceXpc3)

##                   mean  CI_low CI_high
## intercept        5.911   5.805   6.017
## se.intercept     0.466   0.431   0.502
## pval.intercept   0.000   0.000   0.000
## estimate       -46.948 -68.314 -25.583
## se.estimate     57.346  55.642  59.049
## pval.estimate    0.364   0.266   0.462

Max Frequency ~ Beak Elongation

MaxFrXpc1 <- intra_phylm(formula = NoteMaxFreq ~ BeakPC1, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda")

summary(MaxFrXpc1)

##                      mean     CI_low    CI_high
## intercept        3335.965   3308.275   3363.655
## se.intercept      200.818    182.390    219.245
## pval.intercept      0.000      0.000      0.000
## estimate       -13156.906 -13659.384 -12654.427
## se.estimate      3103.768   3050.647   3156.888
## pval.estimate       0.000      0.000      0.000

Max Frequency ~ Beak Curvature

MaxFrXpc2 <- intra_phylm(formula = NoteMaxFreq ~ BeakPC2, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda")

summary(MaxFrXpc2)

##                     mean   CI_low   CI_high
## intercept       3327.440 3296.982  3357.899
## se.intercept     175.700  158.320   193.080
## pval.intercept     0.000    0.000     0.000
## estimate        7817.056 6446.763  9187.349
## se.estimate    10079.611 9877.069 10282.154
## pval.estimate      0.468    0.390     0.545

Max Frequency ~ Beak Tapering

MaxFrXpc3 <- intra_phylm(formula = NoteMaxFreq ~ BeakPC3, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda")

summary(MaxFrXpc3)

##                     mean    CI_low   CI_high
## intercept       3309.306  3276.896  3341.716
## se.intercept     153.938   144.412   163.463
## pval.intercept     0.000     0.000     0.000
## estimate       41888.087 38390.552 45385.622
## se.estimate    20228.311 19795.111 20661.510
## pval.estimate      0.061     0.039     0.083

Max Frequency ~ Beak Size

MaxFrXsize <- intra_phylm(formula = NoteMaxFreq ~ Csize, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda")

summary(MaxFrXsize)

##                    mean   CI_low  CI_high
## intercept      4343.652 4267.257 4420.046
## se.intercept    471.454  461.614  481.293
## pval.intercept    0.000    0.000    0.000
## estimate        -66.246  -70.930  -61.562
## se.estimate      27.694   27.112   28.275
## pval.estimate     0.029    0.019    0.039

Max Frequency ~ Body Size

MaxFrXbody <- intra_phylm(formula = NoteMaxFreq ~ BodyMass.Value, data = SongShape_var,phy = TreeSP, "NoteMaxFreq_var",model = "lambda")

summary(MaxFrXbody)

##                    mean   CI_low  CI_high
## intercept      4052.587 4015.561 4089.614
## se.intercept    218.577  211.678  225.477
## pval.intercept    0.000    0.000    0.000
## estimate        -22.410  -23.142  -21.677
## se.estimate       4.841    4.749    4.932
## pval.estimate     0.000    0.000    0.000

Min Frequency ~ Beak Elongation

MinFrXpc1 <- intra_phylm(formula = NoteMinFreq ~ BeakPC1, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda")

summary(MinFrXpc1)

##                     mean    CI_low   CI_high
## intercept       1989.871  1972.293  2007.450
## se.intercept     143.632   130.963   156.301
## pval.intercept     0.000     0.000     0.000
## estimate       -7398.798 -7851.564 -6946.032
## se.estimate     2299.477  2255.054  2343.901
## pval.estimate      0.006     0.001     0.011

Min Frequency ~ Beak Curvature

MinFrXpc2 <- intra_phylm(formula = NoteMinFreq ~ BeakPC2, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda")

summary(MinFrXpc2)

##                    mean   CI_low  CI_high
## intercept      1969.367 1948.238 1990.497
## se.intercept    129.283  116.791  141.775
## pval.intercept    0.000    0.000    0.000
## estimate       4754.550 3515.735 5993.365
## se.estimate    7218.194 7068.742 7367.647
## pval.estimate     0.511    0.420    0.602

Min Frequency ~ Beak Tapering

MinFrXpc3 <- intra_phylm(formula = NoteMinFreq ~ BeakPC3, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda")

summary(MinFrXpc3)

##                     mean    CI_low   CI_high
## intercept       1943.096  1921.507  1964.684
## se.intercept     125.681   116.280   135.082
## pval.intercept     0.000     0.000     0.000
## estimate       19187.559 16075.869 22299.250
## se.estimate    14899.351 14565.278 15233.424
## pval.estimate      0.258     0.177     0.339

Min Frequency ~ Beak Size

MinFrXsize <- intra_phylm(formula = NoteMinFreq ~ Csize, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda")

summary(MinFrXsize)

##                    mean   CI_low  CI_high
## intercept      2674.778 2638.337 2711.218
## se.intercept    350.790  342.276  359.305
## pval.intercept    0.000    0.000    0.000
## estimate        -44.830  -46.985  -42.675
## se.estimate      20.479   20.036   20.922
## pval.estimate     0.038    0.028    0.049

Min Frequency ~ Body Size

MinFrXbody <- intra_phylm(formula = NoteMinFreq ~ BodyMass.Value, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda")

summary(MinFrXbody)

##                    mean   CI_low  CI_high
## intercept      2448.411 2411.880 2484.942
## se.intercept    155.698  149.760  161.637
## pval.intercept    0.000    0.000    0.000
## estimate        -14.801  -15.513  -14.089
## se.estimate       3.468    3.380    3.557
## pval.estimate     0.000    0.000    0.000

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

MinFrXsize.r <- intra_phylm(formula = NoteMinFreq ~ Resid.Beak.Size, data = SongShape_var,phy = TreeSP, "NoteMinFreq_var",model = "lambda")

summary(MinFrXsize.r)

##                    mean   CI_low  CI_high
## intercept      1975.150 1953.452 1996.849
## se.intercept    147.628  131.631  163.625
## pval.intercept    0.000    0.000    0.000
## estimate       1510.163 1249.113 1771.212
## se.estimate    1112.452 1092.335 1132.570
## pval.estimate     0.254    0.167    0.342

That’s all! Thanks for reading.