Speciation Rate

To obtain speciation rate estimates, we used Bayesian Analysis of Macroevolutionary Mixtures (BAMM). Additionally, we compared our BAMM output files (event data) to an existing rate dataset (Harvey et al. 2017).

Rather than obtaining one estimate of speciation rate from one tree, we ran BAMM 100 times on 100 trees plus an MCC tree to obtain uncertainty estimates of the speciation rate generated from the phylogeny produced by Jetz et al. (2012). We also estimated extinction rate from BAMM and speciation rates using the diversification rate and node density statistic (\(\lambda_{DR}\) and (\(\lambda_{ND}\)). See Title and Rabosky (2018) for a comparison of methods.

Sexual dichromatism and Male-biased Sexual Selection

The first proxy for sexual selection that we used was based on measures of sexual dichromatism. We used two existing datasets; a complete set of male and female plumage scores in Passerines from Dale et al. (2015) and reflectance data from 1000 birds from Armenta et al. (2008).

Briefly, Dale et al. (2015) obtained male and female plumage score, the ratio of which provides an index of sexual dichromatism (SDi). The plumage score was extracted using the Handbook of the Birds of the World Del Hoyo et al. (2011). By scanning images of males and females in this book across multiple patches, Dale et al. (2015) obtained mean plumage scores from RGB values. Dale et al. (2015) compared their RGB scores to a dataset of reflectance measurements of Australian birds. The two datasets were correlated but the relationship was expectedly noisy (R² = 0.67). Here we compare the Dale et al. (2015) data against another reflectance dataset from Armenta et al. (2008), which estimated sexual dichromatism as the mean number of just noticeable differences (JNDs) between male and female plumage colour (discriminability), based on a model of bird colour vision.

Alongside sexual dichromatism, Dale et al. (2015) compiled available data for traits such as: Body size, tropical life history, Sexual selection, Cooperative breeding and Migration. As a second measure of sexual selection, we used a phylogenetic PCA (PC1) that estimates male-biased sexual selection from a species level dataset of sexual size dimorphism, social polygyny and paternal care. Higher values indicate higher size dimorphism (larger males), higher polygyny and less paternal care.

Environmental predictors

We obtained species range distribution maps from Birdlife International (BirdLife International and Handbook of the Birds of the World 2017). These species range maps cover nearly all species of birds. From these spatial information files we obtained estimates of the following (note that not all the data extracted was subsequently used in the analysis):

• Range Size
• Average and standard deviation in 19 bioclimatic variables (each range randomly sampled 1000 times)
• Average and standard deviation in 19 bioclimatic variables from the last-glacial maximum (LGM) and the last-inter-glacial (LIG)
• Net primary productivity (NPP) estimates and variability
• Average and standard deviation of human population density in the species ranges.

Below we outline the code used for extracting spatial data from bird ranges, however no raw data is provided alongside this file due to the large file size of the shapefile and raster information. We have put the extracted environmental variables in a csv filed called complete.dataframe.csv for convenience. Briefly for each of the ~ 6,000 species, we:

1) Randomly sampled the range polygon 1000 times:

#Increse the iterations (defaukt is 4) so we can obtain complete samples of each range
bird.points <- lapply(bird.ranges, 
                       function(x) {spsample(x, n=1000, type="random", iter = 30)})

2) Extracted the bioclim or other (e.g. NPP) data from each of the points:

bird.values <- lapply(bird.points, function(x) {raster::extract(bioclim_data, x)})
bird.values <- as.list(data.frame((bird.values)))

3) Summarised the extracted data across the 1,000 points into a summary value of interest (means and sd) for each species and exported that data.

#Obtain means per variable per species
bird.frame <- lapply(bird.values, function(x) {as.data.frame(x)})
bird.summary <- lapply(bird.frame, function(x) {
(as.data.frame(apply(x,2,mean, na.rm =T)))})

#transpose
bird.means <- t(as.data.frame(bird.summary))
bird.means <- (split(bird.means, 1:19))

#Now add column of species name: Same order carried through
bird.means <- cbind.data.frame(bird.names, bird.means)

rownames(bird.means) <- NULL
colnames(bird.means) <- c("binomial","bioclim1", "bioclim2", "bioclim3", "bioclim4", "bioclim5", "bioclim6", "bioclim7", "bioclim8", "bioclim9", "bioclim10", "bioclim11", "bioclim12", "bioclim13", "bioclim14", "bioclim15", "bioclim16", "bioclim17", "bioclim18", "bioclim19")

#Write csv
write.csv(bird.means, 'data/bird.means.csv')

4) The process above was repeated when we used gridded data for the LGM, LIG, NPP and human population density, while range size was obtained from the following code:

bird.range.size <- sapply(bird.ranges,
                       function(x) {(area(x))})
bird.range.size <- sapply(bird.range.size, function(x) {sum(x)})
bird.range.size <- as.data.frame(bird.range.size)
bird.range.size <- cbind.data.frame(shps.jetz, bird.range.size)
rownames(bird.range.size) <- NULL
colnames(bird.range.size) <- c("binomial", "range.size.m2")

#write.csv
write.csv(bird.range.size, 'data/bird.range.size.csv')

Environmental Spatial Variability PCAs

To obtain estimates of spatial variability in environment, we performed PCAs of the standard errors of the bioclimatic variables across the 1000 random points extracted per species (excluding seasonality in temperature and precipitation). We expect variability to increase with increased range size. To correct for this association, and to be able to include both variables in the model, we took the residuals from GAM models of the relationship between PCA components and range size.

Table S1: Loadings for the first two PCs, from a PCA on the variation (standard error, where n ~ 1,000) of each bioclim variable except the standard deviation of temperature seasonality (sd.bioclim4) and the standard error of precipitation seasonality (sd.bioclim15), as these are already measures of variability and based on other bioclimatic variables.

restricted.data <- complete.dataframe %>% drop_na(bioclim1) #Drop species without environmental data
PCA.bioclim <- prcomp(restricted.data[c('sd.bioclim1', 'sd.bioclim2', 'sd.bioclim3', 'sd.bioclim5', 'sd.bioclim6', 'sd.bioclim7', 'sd.bioclim8', 'sd.bioclim9', 'sd.bioclim10', 'sd.bioclim11', 'sd.bioclim12', 'sd.bioclim13', 'sd.bioclim14', 'sd.bioclim16', 'sd.bioclim17', 'sd.bioclim18', 'sd.bioclim19')], #Select sd.bioclim
                      scale = TRUE, center = TRUE)
PCA.predictions <- predict(PCA.bioclim)
restricted.data <- cbind(restricted.data, PCA.predictions)

#Create table with PC1 and PC2
as.data.frame(PCA.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1", "PC2")) %>% pander()

#run GAM of association between range size and spatial climatic variation
PC1.gam <- mgcv::gam(PC1*(-1) ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian") #Inverse as the PC1 loads to the negative (COUNTER-INTUITIVE)

#take residuals
restricted.data$residuals.PC1 <- residuals.gam(PC1.gam) 

#We can do the same for PC2
PC2.gam <- mgcv::gam(PC2 ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian")

#Take residuals
restricted.data$residuals.PC2 <- residuals.gam(PC2.gam)

#plot
par(mfrow = c(1,2))
plot.gam(PC1.gam, residuals = T, main = 'PC1 residuls (temp)')
plot.gam(PC2.gam, residuals = T, main = 'PC2 residuals (precip)')

Figure S1: The relationship between spatial variability in the temperature components (PC1) and log-range size is relatively strong. The relationship between spatial variability in precipitation (PC2) and log-range size is weaker. In the analysis we only used the residuals from PC1. PC1 accounts for 48.1 % of the variation.

Long term climate variability (LIG anomaly)

Climate stability through time can potentially affect diversification dynamics. To gain estimates for change in climate over the past ~130,000 years we used the difference in bioclim variables between the LIG and present values. The plots show the two PCs, broadly representing temperature and precipitation. We used the difference between the present and LIG climates as these represent a longer (evolutionarily meaningful) time-scale than the difference between the LGM and present climates.

Table S2: Loadings for the first two PCs of each of the PCAs for the difference in bioclimatic variables between today and the LIG. Here, PC1 is more heavily loaded for absolute temperature difference, while PC2 is more heavily loaded for absolute difference in precipitation.

historical.variation.data <- as.data.frame(restricted.data[1])

#FOR LIG

historical.variation.data$bio1.LIG.diff <- abs(restricted.data$bioclim1 - restricted.data$LIG.bi1)
historical.variation.data$bio2.LIG.diff <- abs(restricted.data$bioclim2 - restricted.data$LIG.bi2)
historical.variation.data$bio3.LIG.diff <- abs(restricted.data$bioclim3 - restricted.data$LIG.bi3)
historical.variation.data$bio4.LIG.diff <- abs(restricted.data$bioclim4 - restricted.data$LIG.bi4)
historical.variation.data$bio5.LIG.diff <- abs(restricted.data$bioclim5 - restricted.data$LIG.bi5)
historical.variation.data$bio6.LIG.diff <- abs(restricted.data$bioclim6 - restricted.data$LIG.bi6)
historical.variation.data$bio7.LIG.diff <- abs(restricted.data$bioclim7 - restricted.data$LIG.bi7)
historical.variation.data$bio8.LIG.diff <- abs(restricted.data$bioclim8 - restricted.data$LIG.bi8)
historical.variation.data$bio9.LIG.diff <- abs(restricted.data$bioclim9 - restricted.data$LIG.bi9)
historical.variation.data$bio10.LIG.diff <- abs(restricted.data$bioclim10 - restricted.data$LIG.bi10)
historical.variation.data$bio11.LIG.diff <- abs(restricted.data$bioclim11 - restricted.data$LIG.bi11)
historical.variation.data$bio12.LIG.diff <- abs(restricted.data$bioclim12 - restricted.data$LIG.bi12)
historical.variation.data$bio13.LIG.diff <- abs(restricted.data$bioclim13 - restricted.data$LIG.bi13)
historical.variation.data$bio14.LIG.diff <- abs(restricted.data$bioclim14 - restricted.data$LIG.bi14)
historical.variation.data$bio15.LIG.diff <- abs(restricted.data$bioclim15 - restricted.data$LIG.bi15)
historical.variation.data$bio16.LIG.diff <- abs(restricted.data$bioclim16 - restricted.data$LIG.bi16)
historical.variation.data$bio17.LIG.diff <- abs(restricted.data$bioclim17 - restricted.data$LIG.bi17)
historical.variation.data$bio18.LIG.diff <- abs(restricted.data$bioclim18 - restricted.data$LIG.bi18)
historical.variation.data$bio19.LIG.diff <- abs(restricted.data$bioclim19 - restricted.data$LIG.bi19)

historical.variation.data <- historical.variation.data %>% drop_na(bio1.LIG.diff)

#Run a PCA of the difference removing the variation variables (4 and 15)

#For LIG
PCA.LIG.bioclim <- prcomp(historical.variation.data[c(
'bio1.LIG.diff',
'bio2.LIG.diff',
'bio3.LIG.diff',
'bio5.LIG.diff',
'bio6.LIG.diff',
'bio7.LIG.diff',
'bio8.LIG.diff',
'bio9.LIG.diff',
'bio10.LIG.diff',
'bio11.LIG.diff',
'bio12.LIG.diff',
'bio13.LIG.diff',
'bio14.LIG.diff',
'bio16.LIG.diff',
'bio17.LIG.diff',
'bio18.LIG.diff',
'bio19.LIG.diff')],
                      scale = TRUE, center = TRUE)

#Create table with PC1 and PC2
as.data.frame(PCA.LIG.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1.LIG", "PC2.LIG")) %>% pander()

#Now we can predict the PCA results: 
PCA.LIG.predictions <- as.data.frame(predict(PCA.LIG.bioclim)*(-1)) #So that higher numbers mean more variation
PCA.LIG.predictions <- rename(PCA.LIG.predictions, PC1.LIG = PC1,
       PC2.LIG = PC2)

#Bind them to dataframe, taking only the first two PCAs
historical.variation.data <- cbind(historical.variation.data, PCA.LIG.predictions[1:2])

#Now to the restricted dataframe
restricted.data <- right_join(restricted.data, historical.variation.data %>% dplyr::select(binomial, PC1.LIG, PC2.LIG), by = 'binomial')

PCA.plot <- historical.variation.data %>% ggplot(aes(x = PC1.LIG, y= PC2.LIG))+
  geom_point(shape = 21)+
  theme_minimal()

PCA.plot.m <- ggExtra::ggMarginal(
  p = PCA.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

grid.arrange(PCA.plot.m)

Figure S2: The distribution of the two PCA components for the absolute difference in bioclimatic variables between today and the LIG.

Correlations between environmental predictors

We checked the correlations between environmental predictors used in subsequent phylogenetic comparative analyses (PGLS models). Specifically we tested whether the following environmental variables are correlated:

• Range size
  
• Short temporal variability of temperature (mean BIOCLIM4)
  
• Spatial variability of temperature (PCA1) [residual.PCA1]
  
• Long-term variability of temperature (LIG)
  
• NPP

To check the correlations between the environmental predictors to be used in the PGLS models, we inspected the following correlation plots with the correlation value plotted:

#Draw plot
restricted.data$log.range.size <- log(restricted.data$range.size.m2)
pairs(restricted.data[,c("log.range.size", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP")], lower.panel=panel.smooth, upper.panel=panel.cor)

Figure S3: The highest correlation is between Long term temperature variation and seasonal variation. However, this correlation is moderate (0.53) so we can be confident that collinearity with not be an issue for our PGLS models.

Sexual Dimorphism

The sexual dimorphism dataset looks to have an over-dispersed distribution (see phylogenic plot in manuscript). Unfortunately, transformations do not greatly improve the distribution. Our model fit is expected to be reduced by this but further transformations are not obvious and would risk problems of interpretation. In any case, the phylogenetically corrected models we employ do not assume normality in the response variable.

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

We obtained two tip-rate metrics of speciation using statistics derived from the properties of the nodes and branches along root-to-tip paths of the phylogeny. Node density (ND) is a simple statistic calculating the density of nodes from the phylogenetic root-to-tip, while the diversification rate (DR) (e.g., Jetz et al. 2012; Quintero and Jetz 2018; Rabosky et al. 2018), is derived from the sum of edge lengths branching from a node, with each more basal node having the sum of lengths down-weighted.

The functions that generate these tip rates estimates are: calculate.log.es, calculate.nd and calculate.tb. Here they are run on 1000 trees plus a full MCC tree of passerine birds based on 2500 sample trees of the posterior (Jetz et al. 2012). This was done using the phanghorn package (Schliep 2010).

#Obtain DR (log(es)) estimates, by calling the calculate.dr function
passerine.trees <- read.nexus('data/trees/Passerine_Trees_Full.nex')

#DR, ND and TB can be calculated on the 100 trees by: 
  
# es.list <- sapply(passerine.trees, calculate.log.es)
# nd.list <- sapply(passerine.trees, calculate.nd)
# tb.list <- sapply(passerine.trees, calculate.tb)

#To avoid re-running the time-consuming rates, we can save them and load them: 
es.list <- readRDS("data/es.list.rds")
nd.list <- readRDS("data/nd.list.rds")
tb.list <- readRDS("data/tb.list.rds")

#Also load in our MCC tree
MCC.passerine <- read.tree('data/MCC.passerine.tre') #Tree based on 2500 samples of posterior

#For ND
#transpose the list so that the elements are the species
t.nd.list <- nd.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
nd.values <- as.data.frame(t(sapply(t.nd.list, function(x) sapply(x, function(x) (x)))))
nd.values <- nd.values %>% tibble::rownames_to_column("binomial")
nd.summary <- melt(nd.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.nd = mean(value), #use normal mean, not a rate
  var.nd = var(value), #Don't really need this if we use CV
  CV.nd = (sd(value)/mean(value)) #Given normal distribution of ND, standard CV formula is appropriate
)

#For ES
#transpose the list so that the elements are the species
t.es.list <- es.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
es.values <- as.data.frame(t(sapply(t.es.list, function(x) sapply(x, function(x) (x)))))
es.values <- es.values %>% tibble::rownames_to_column("binomial")
es.summary <- melt(es.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.loges = log(1/mean(1/exp(value))), #use harmonic mean for rates and log-transform post calculation of HM
  var.loges = var(value), #Don't really need this if we use CV
  CV.loges = sqrt(exp(var(value)) - 1) #CV for log-normal distribution of ES
)

#For TB
#transpose the list so that the elements are the species
t.tb.list <- tb.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
tb.values <- as.data.frame(t(sapply(t.tb.list, function(x) sapply(x, function(x) (x)))))
tb.values <- tb.values %>% tibble::rownames_to_column("binomial")
tb.summary <- melt(tb.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.tb = mean(value), #use normal mean; not a rate
  var.tb = var(value), #Don't really need this if we use CV
  CV.tb = (sd(value)/mean(value)) #Given normal distribution of tb, standard CV formula is appropriate
)

MCC.nd.df <- as.data.frame(calculate.nd(MCC.passerine)) %>% tibble::rownames_to_column("binomial")
MCC.dr.df <- as.data.frame(calculate.log.es(MCC.passerine)) %>% tibble::rownames_to_column("binomial")

#Join the dataframes
dr.summary <- full_join(nd.summary, es.summary, by = "binomial")
dr.summary <- full_join(dr.summary, tb.summary, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.nd.df, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.dr.df, by = "binomial")

dr.summary <- dr.summary %>% rename(TipLabel = binomial, MCC.ND = nd, MCC.DR = es)

saveRDS(dr.summary, 'data/dr.summary.rds')

dr.summary <- readRDS('data/dr.summary.rds')
#Plot a summary of logES and TB with the DRs with their weightings
dr.plot <- dr.summary %>% ggplot(aes(x = MCC.DR, y = MCC.ND, 
                          fill = MCC.DR, 
                      size = MCC.ND))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  theme(legend.position="bottom")+
  labs(size = 'Node Density [ND]', y='Node Density [ND]', x= 'Diversification Rate [DR]\n (log-equal splits)', fill = 'Diversification Rate [DR]')+
  scale_fill_distiller(guide = "colorbar", palette = "RdPu", direction = 1)
marginal.dr <-ggMarginal(
  p = dr.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = '#0000009C',
  fill = '#D12E6769'
)

grid.newpage()
grid.draw(marginal.dr)

Figure S4: \(\lambda_{DR}\) and \(\lambda_{ND}\) estimates of speciation rate are correlated for the MCC tree used in this study. Both of these response variables are normally distributed.

restricted.data <- left_join(restricted.data, dr.summary, by = 'TipLabel')

es.plot <- restricted.data %>% ggplot(aes(y = mean.loges, x = SDi, fill = mean.loges))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "YlOrBr", direction = 1, guide = FALSE)+
ylab("mean.DR")+
theme_minimal()

nd.plot <- restricted.data %>% ggplot(aes(y = mean.nd, x = SDi, fill = mean.nd))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "Greens", direction = 1, guide = FALSE)+
theme_minimal()

tb.plot <- restricted.data %>% ggplot(aes(y = mean.tb, x = SDi, fill = mean.tb))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "PuBu", direction = 1, guide = FALSE)+
theme_minimal()

grid.arrange(es.plot, nd.plot, tb.plot, nrow = 3)

Figure S5: Scatter plots showing the raw relationship between sexual dichromatism (SDi) and speciation rates. Across all three measures of speciation the pattern and spread is similar, with no obvious relationship. Note that terminal branch length (mean.tb) was not used in the analysis.

BAMM measures of speciation and extinction

name.passerine.tree <- names(passerine.trees)

priors <- sapply(name.passerine.tree, function(x) {
  setBAMMpriors(passerine.trees[[x]], outfile = NULL)
})

sapply(name.passerine.tree, function(x) {
  write.tree(passerine.trees[[x]], paste("data/bamm_files/", x, ".tre", sep=""))
})

# Here is a block of parameters for the control file. We can make a control file for each tree:
params <- list()
for (x in name.passerine.tree) {

# GENERAL SETUP AND DATA INPUT

params[[x]] <- list(modeltype = 'speciationextinction',
# Specify "speciationextinction" or "trait" analysis
                                  
treefile = paste(x, ".tre", sep=""),
# File name of the phylogenetic tree to be analyzed

runInfoFilename = 'run_info.txt',
# File name to output general information about this run

sampleFromPriorOnly = 0,
# Whether to perform analysis sampling from prior only (no likelihoods computed)

runMCMC = 1,
# Whether to perform the MCMC simulation. If runMCMC = 0, the program will only
# check whether the data file can be read and the initial likelihood computed

loadEventData = 0,                       
# Whether to load a previous event data file

eventDataInfile = 'event_data_in.txt',
# File name of the event data file to load, used only if loadEventData = 1

initializeModel = 1,
# Whether to initialize (but not run) the MCMC. If initializeModel = 0, the
# program will only ensure that the data files (e.g., treefile) can be read

useGlobalSamplingProbability = 1,
# Whether to use a "global" sampling probability. If False (0), expects a file
# name for species-specific sampling probabilities (see sampleProbsFilename)
                                        
globalSamplingFraction = 1,
# The sampling probability. If useGlobalSamplingProbability = 0, this is ignored
# and BAMM looks for a file name with species-specific sampling fractions

sampleProbsFilename = 'sample_probs.txt',
# File name containing species-specific sampling fractions

seed = as.numeric(gsub("tree_", "", x, perl = TRUE)),
# Seed for the random number generator. Set for reproducibility to the number of the treefile

overwrite = 1,
# If True (1), the program will overwrite any output files in the current
# directory (if present)


# PRIORS

expectedNumberOfShifts = 100,
# prior on the number of shifts in diversification
# Suggested values: 
#     expectedNumberOfShifts = 1.0 for small trees (< 500 tips)
#  expectedNumberOfShifts = 10 or even 50 for large trees (> 5000 tips) 
 
lambdaInitPrior = as.numeric(priors['lambdaInitPrior', x]),
# Prior (rate parameter of exponential) on the initial lambda value for rate
# regimes

lambdaShiftPrior = 0.05,
# Prior (std dev of normal) on lambda shift parameter for rate regimes
# You cannot adjust the mean of this distribution (fixed at zero, which is
# equal to a constant rate diversification process)

muInitPrior = as.numeric(priors['muInitPrior', x]),
# Prior (rate parameter of exponential) on extinction rates  

lambdaIsTimeVariablePrior = 1,
# Prior (probability) of the time mode being time-variable (vs. time-constant)
            

# MCMC SIMULATION SETTINGS & OUTPUT OPTIONS

numberOfGenerations = '100000000',
# Number of generations to perform MCMC simulation

mcmcOutfile = 'mcmc_out.txt',
# File name for the MCMC output, which only includes summary information about
# MCMC simulation (e.g., log-likelihoods, log-prior, number of processes)

mcmcWriteFreq = 1000,
# Frequency in which to write the MCMC output to a file

eventDataOutfile = 'event_data.txt',
# The raw event data (these are the main results). ALL of the results are
# contained in this file, and all branch-specific speciation rates, shift
# positions, marginal distributions etc can be reconstructed from this output.
# See R package BAMMtools for working with this output

eventDataWriteFreq = 1000,
# Frequency in which to write the event data to a file

printFreq = 10000,
# Frequency in which to print MCMC status to the screen

acceptanceResetFreq = 1000,
# Frequency in which to reset the acceptance rate calculation
# The acceptance rate is output to both the MCMC data file and the screen

outName = x,
# Optional name that will be prefixed on all output files (separated with "_")
# If commented out, no prefix will be used


# OPERATORS: MCMC SCALING OPERATORS

updateLambdaInitScale = 2,
# Scale parameter for updating the initial speciation rate for each process

updateLambdaShiftScale = 0.1,
# Scale parameter for the exponential change parameter for speciation

updateMuInitScale = 2,
# Scale parameter for updating initial extinction rate for each process

updateEventLocationScale = 0.1,
# Scale parameter for updating LOCAL moves of events on the tree
# This defines the width of the sliding window proposal
 
updateEventRateScale = 4,
# Scale parameter (proportional shrinking/expanding) for updating
# the rate parameter of the Poisson process

# OPERATORS: MCMC MOVE FREQUENCIES

updateRateEventNumber = 1,
# Relative frequency of MCMC moves that change the number of events

updateRateEventPosition = 0.25,
# Relative frequency of MCMC moves that change the location of an event on the
# tree

updateRateEventRate = 1,
# Relative frequency of MCMC moves that change the rate at which events occur 

updateRateLambda0 = 1,
# Relative frequency of MCMC moves that change the initial speciation rate
# associated with an event

updateRateLambdaShift = 1,
# Relative frequency of MCMC moves that change the exponential shift parameter
# of the speciation rate associated with an event

updateRateMu0 = 1,
# Relative frequency of MCMC moves that change the extinction rate for a given
# event

updateRateLambdaTimeMode = 0,
# Relative frequency of MCMC moves that flip the time mode
# (time-constant <=> time-variable)

localGlobalMoveRatio = 10,
# Ratio of local to global moves of events 


# INITIAL PARAMETER VALUES

lambdaInit0 = 0.032,
# Initial speciation rate (at the root of the tree)

lambdaShift0 = 0,
# Initial shift parameter for the root process

muInit0 = 0.005,
# Initial value of extinction (at the root)

initialNumberEvents = 0,
# Initial number of non-root processes


# METROPOLIS COUPLED MCMC

numberOfChains = 1,
# Number of Markov chains to run

deltaT = 0.01,
# Temperature increment parameter. This value should be > 0
# The temperature for the i-th chain is computed as 1 / [1 + deltaT * (i - 1)]

swapPeriod = 1000,
# Number of generations in which to propose a chain swap

chainSwapFileName = 'chain_swap.txt',
# File name in which to output data about each chain swap proposal.
# The format of each line is [generation],[rank_1],[rank_2],[swap_accepted]
# where [generation] is the generation in which the swap proposal was made,
# [rank_1] and [rank_2] are the chains that were chosen, and [swap_accepted] is
# whether the swap was made. The cold chain has a rank of 1.


# NUMERICAL AND OTHER PARAMETERS

minCladeSizeForShift = 3,
# Allows you to constrain location of possible rate-change events to occur
# only on branches with at least this many descendant tips. A value of 1
# allows shifts to occur on all branches. 

segLength = 0.025,
# Controls the "grain" of the likelihood calculations. Approximates the
# continuous-time change in diversification rates by breaking each branch into
# a constant-rate diversification segments, with each segment given a length
# determined by segLength. segLength is in units of the root-to-tip distance of
# the tree. So, if the segLength parameter is 0.01, and the crown age of your
# tree is 50, the "step size" of the constant rate approximation will be 0.5.
# If the value is greater than the branch length (e.g., you have a branch of
# length < 0.5 in the preceding example) BAMM will not break the branch into
# segments but use the mean rate across the entire branch.

outName = x)
  }

bammcontrolfile <- list()
for (x in name.passerine.tree) {
  bammcontrolfile[x] <- paste("data/bamm_files/control_", x, ".txt", sep="")
}

# Now writing control parameters to file

for (x in name.passerine.tree) {generateControlFile(file = bammcontrolfile[[x]], type = "diversification", params = params[[x]])}

bamm.commands <- list()
for (x in name.passerine.tree) {
  bamm.commands[x] <- paste("bamm -c control_", x, ".txt", sep="")
}

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.mcmc <- read.csv( "data/BAMM_MCC/tree_MCC_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.mcmc))
postburn <- MCC.BAMM.mcmc[burnstart:nrow(MCC.BAMM.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()

ESS <- readRDS('data/ESS.rds')
summary(ESS) %>% pander()

# Read in Event Data
MCC.BAMM.ED <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/tree_MCC_event_data.txt", burnin=0.1, nsamples=1000)

## Reading event datafile:  data/BAMM_MCC/tree_MCC_event_data.txt 
##      ...........
## Read a total of 100000 samples from posterior
## 
## Discarded as burnin: GENERATIONS <  9999000
## Analyzing  1000  samples from posterior
## 
## Setting recursive sequence on tree...
## 
## Done with recursive sequence

saveRDS(MCC.BAMM.ED, 'data/MCC.BAMM.ED.rds')
#From the Event Data we can extract 

library(purrr)
#BAMM.EventData <- readRDS('data/BAMM.EventData.rds')

#Big lapply over each tree in BAMM event data
BAMM.extraction.function <- function(x) {
######Get mean and var for lambda
#Transpose list so each element in the list is a species
transposed.lambda <- lapply(purrr::transpose(x$tipLambda), unlist)

#Now turn it into a df with mean and variance
lambda <- sapply(transposed.lambda, function(x) {
  mean.lambda = mean(log(x))
  var.lambda = var(log(x))
  return(c(mean.lambda, var.lambda))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.lambda", "var.lambda"))

lambda$TipLabel <- x[["tip.label"]]

#####NOW FOR Extinction

#Transpose list so each element in the list is a species
transposed.mu <- lapply(purrr::transpose(x$tipMu), unlist)

#Now turn it into a df with mean and variance
mu <- sapply(transposed.mu, function(x) {
  mean.mu = mean(log(x))
  var.mu = var(log(x))
  return(c(mean.mu, var.mu))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.mu", "var.mu"))

mu$TipLabel <- x[["tip.label"]]

left_join(lambda, mu, by = "TipLabel")
}

MCC.BAMM.df <- BAMM.extraction.function(MCC.BAMM.ED)

#Save df for later use
saveRDS(MCC.BAMM.df, 'data/MCC.BAMM.df.rds')

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
MCC.BAMM.df$CV.lambda <- sqrt(exp(MCC.BAMM.df$var.lambda) - 1)*100
#Plot a summary of logES and TB with the DRs with their weightings
BAMM.MCC.plot <- MCC.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/var.lambda, 
                      size = 1/var.mu))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)),
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)),
                size = 0.0025)+
  # scale_y_continuous(trans = "log")+
  # scale_x_continuous(trans = "log")+
  # xlim(-8, 3)+
  # ylim(-10,2)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Inverse log(var) / Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Inverse log(var) / Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

BAMM.variance <- ggExtra::ggMarginal(
  p = BAMM.MCC.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = '#BA3B1C91'
)

grid.newpage()
grid.draw(BAMM.variance)

css <- credibleShiftSet(MCC.BAMM.ED, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

summary(css)

## 
##  95 % credible set of rate shift configurations sampled with BAMM
## 
## Distinct shift configurations in credible set:  950
## 
## Frequency of 9 shift configurations with highest posterior probability:
## 
## 
##    rank     probability cumulative  Core_shifts
##          1      0.001      0.001         54
##          2      0.001      0.002         44
##          3      0.001      0.003         46
##          4      0.001      0.004         47
##          5      0.001      0.005         49
##          6      0.001      0.006         45
##          7      0.001      0.007         47
##          8      0.001      0.008         43
##          9      0.001      0.009         45
## 
## ...omitted 941 additional distinct shift configurations
## from the credible set. You can access the full set from your 
## credibleshiftset object

round(computeBayesFactors(MCC.BAMM.mcmc, expectedNumberOfShifts=100, burnin=0.1)[,1], digits = 2)

##      42      43      44      45      46      47      48      49      50 
##    1.00   17.17   24.48   49.45  120.71  283.77  454.33  784.80 1291.85 
##      51      52      53      54      55      56      57      58      59 
## 1813.33 2561.62 3608.07 4574.91 5492.44 6322.11 7141.12 7720.26 7886.28 
##      60      61      62      63      64      65      66      67      68 
## 7846.73 7661.83 7116.15 6466.36 5749.34 4811.16 4136.80 3338.17 2674.70 
##      69      70      71      72      73      74      75      76      77 
## 2137.61 1593.48 1218.40  868.01  643.91  452.36  311.07  199.17  135.99 
##      78      79      80      81      82      83      84      85      87 
##   85.85   65.03   40.87   14.74   13.40    9.02   10.63    1.53    1.56

plotPrior(MCC.BAMM.mcmc, expectedNumberOfShifts=100)

load('data/Hackett_split_eventsample.rda')
css3 <- credibleShiftSet(ed, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

Harvey.BAMM.df <- BAMM.extraction.function(ed)

Harvey.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/(var.lambda), 
                      size = 1/(var.mu)))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)), 
                size = 0.0025)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

Harvey.vs.us <- full_join(MCC.BAMM.df %>% select(TipLabel, mean.lambda), Harvey.BAMM.df %>% select(TipLabel, Harvey.lambda = mean.lambda), by = "TipLabel")
saveRDS(Harvey.vs.us, "data/Harvey.vs.us.rds")
# BAMM.variance <- ggExtra::ggMarginal(
#   p = BAMM.MCC.plot,
#   type = 'density',
#   margins = 'both',
#   size = 5,
#   colour = 'black',
#   fill = '#BA3B1C91'
# )

PGLS Models

The method behind running the PGLS models is as follows:

1: Estimate the phylogenetic signal (\(\lambda\)) by running a model without interactions and all six predictor variables. The value of \(\lambda\) obtained here will be fixed in all successive models. The value is fixed for subsequent models as independently estimating it in each case becomes computationally intensive.

2: Create a global model of six predictor variables plus the five interactions between sexual dichromatism and the other variables.

3: Dredge the global model but fix the six independent predictors, hence conducting model selection on the interaction terms.

4: Take the top model in the MCC model and run it on the 100 phylogenetic trees.

5: Repeat this step for DR, ND, BAMM-speciation and BAMM-extinction.

#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(restricted.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(restricted.data) <- restricted.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.corPagel <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.DR.corPagel, 'data/MCC.DR.corPagel.rds')

#Run a corPagel model to estimate lambda for ND
MCC.ND.corPagel <- gls(MCC.ND ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.ND.corPagel, 'data/MCC.ND.corPagel.rds')

MCC.DR.corPagel <- readRDS('data/MCC.DR.corPagel.rds')
MCC.ND.corPagel <- readRDS('data/MCC.ND.corPagel.rds')
MCC.DR.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()

MCC.ND.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("ND lambda") %>% pander()

#Run model for DR
MCC.DR <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Run model for ND
MCC.ND <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP
                      + SDi*log(range.size.m2)
                      + SDi*bioclim4
                      + SDi*residuals.PC1
                      + SDi*PC1.LIG
                      + SDi*NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("restricted.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.ND.model <- pdredge(MCC.ND, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.ND.model, "data/dredged.ND.model.rds")

dredged.DR.model <- pdredge(MCC.DR, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.DR.model, "data/dredged.DR.model.rds")

dredged.DR.model <- readRDS("data/dredged.DR.model.rds")
dredged.ND.model <- readRDS("data/dredged.ND.model.rds")
kable(dredged.DR.model, "html", caption = "DR Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

kable(dredged.ND.model, "html", caption = "ND Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
MCC.DR.top <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.DR.top, 'data/MCC.DR.top.rds')

#Run model for ND
MCC.ND.top <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.ND.top, 'data/MCC.ND.top.rds')
#Run the 100 models for DR and ND using the best model:


# No longer used now that we use 1000 trees

# #Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
# DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
#   left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
#             x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ES")), 
#             by = "TipLabel")
# })
# 
# #PGLS needs tiplabel as rowname
# DR.model.data <- lapply(DR.model.data, function(x) {
#   tibble::column_to_rownames(x, "TipLabel")})
# 
# #Prune the trees
# pruned.trees<-lapply(passerine.trees, function(x) {
#   drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
# })
# 
# #Use mapply to create a list of PGLS global models
# DR.pgls.models <- mcmapply(function(x,y) {
#   gls(ES ~ SDi 
#          + log(range.size.m2)
#          + bioclim4 #Seasonal variation
#          + residuals.PC1 #Spatial variation
#          + PC1.LIG #Long-term climate variation
#          + NPP,
#     corPagel(0.985, phy = y, fixed = TRUE), 
#     data = x, 
#     method = "REML")
# }, x = DR.model.data, y = pruned.trees,
# SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
# mc.cores = 8) #Specify core number 
# 
# saveRDS(DR.pgls.models, "data/DR.pgls.models.rds")
# 
# #Now for Node Density:
# ND.model.data <- lapply(nd.list, function(x) {
#   left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
#             x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
#             by = "TipLabel")
# })
# 
# #PGLS needs tiplabel as rowname
# ND.model.data <- lapply(ND.model.data, function(x) {
#   tibble::column_to_rownames(x, "TipLabel")})
# 
# #Use mapply to create a list of PGLS global models
# ND.pgls.models <- mcmapply(function(x,y) {
# gls(ND ~ SDi 
#          + log(range.size.m2)
#          + bioclim4 #Seasonal variation
#          + residuals.PC1 #Spatial variation
#          + PC1.LIG #Long-term climate variation
#          + NPP,
#     corPagel(0.9996, phy = y, fixed = TRUE), 
#     data = x, 
#     method = "REML")
# }, x = ND.model.data, y = pruned.trees, 
# SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
# mc.cores = 8) #Specify core number 
# 
# saveRDS(ND.pgls.models, "data/ND.pgls.models.rds")

#Read in DR models and extract estimates:

#We ran DR and ND models on 1000 trees
files.DR.SD <- list.files(path = "/Users/justincally/Dropbox/Runs Spartan/DR_SD/", pattern = "\\.rds$", full.names = TRUE) #1000 models
df <- list()
lapply(files.DR.SD, function(x) {
  tryCatch({
  model.x <- readRDS(x)
  name.x <- str_sub(x, -19, -5)
  df[[name.x]] <<- data.frame(model.x$coefficients,
                              confint(model.x),
                              coef(summary(model.x))[,2], #Std.Error
                              coef(summary(model.x))[,3], #t-val
                              coef(summary(model.x))[,4], #pval
                              model.x[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = name.x) %>% tibble::rownames_to_column()
  rm(model.x)
  gc(verbose = FALSE)
  },
  error = function(e) NULL
  )
})
DR.pgls.summary <- bind_rows(df) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))
#saveRDS(DR.pgls.summary, "data/DR.pgls.summary.rds") #Save a simple df with the model and coeff


# Read in ND models and extract estimates

files.ND.SD <- list.files(path = "/Users/justincally/Dropbox/Runs Spartan/ND_SD/", pattern = "\\.rds$", full.names = TRUE) #1000 models
df <- list()
lapply(files.ND.SD, function(x) {
  tryCatch({
  model.x <- readRDS(x) 
  name.x <- str_sub(x, -19, -5)
  df[[name.x]] <<- data.frame(model.x$coefficients, 
                              confint(model.x), 
                              coef(summary(model.x))[,2], #Std.Error
                              coef(summary(model.x))[,3], #t-val
                              coef(summary(model.x))[,4], #pval
                              model.x[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = name.x) %>% tibble::rownames_to_column()
  rm(model.x)
  gc(verbose = FALSE)
  },
  error = function(e) NULL
  )
})
ND.pgls.summary <- bind_rows(df) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))
#saveRDS(ND.pgls.summary, "data/ND.pgls.summary.rds") #Save a simple df with the model and coeff

#Read in the 1000 tree summary df and the MCC tree model
DR.pgls.summary <- readRDS("data/DR.pgls.summary.rds")
MCC.DR.top <- readRDS('data/MCC.DR.top.rds')
MCC.DR.summary <-  data.frame(MCC.DR.top$coefficients,
                              confint(MCC.DR.top),
                              coef(summary(MCC.DR.top))[,2], #Std.Error
                              coef(summary(MCC.DR.top))[,3], #t-val
                              coef(summary(MCC.DR.top))[,4], #pval
                              MCC.DR.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))

#Now plot this for Node Density (ND):  

ND.pgls.summary <- readRDS("data/ND.pgls.summary.rds")
MCC.ND.top <- readRDS('data/MCC.ND.top.rds')
MCC.ND.summary <- data.frame(MCC.ND.top$coefficients,
                              confint(MCC.ND.top),
                              coef(summary(MCC.ND.top))[,2], #Std.Error
                              coef(summary(MCC.ND.top))[,3], #t-val
                              coef(summary(MCC.ND.top))[,4], #pval
                              MCC.ND.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )



DR.plot <-DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  #geom_errorbarh(aes(xmin = LCI, xmax = UCI, colour = Parameter), position = position_jitter(seed = 1), height = 0)+
  geom_point(shape = 21, alpha = 0.5, size = 0.75, position = position_jitter(seed = 1))+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

DR.plot <- DR.plot + geom_errorbarh(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("DR Models")


ND.plot <-ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

ND.plot <- ND.plot + geom_errorbarh(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("ND Models")

PGLS Models on BAMM estimates

Using the previously generated BAMM dataframe we undertook model selection for speciation and extinction. The process is similar to that used on DR and ND, except given the use of a Bayesian approach in BAMM we can make use of varying levels of uncertainty between tips (species) by constructing a weighted model, where the weight is the inverse of the variance, such that more precise estimates of speciation or extinction at a given tip (species) holds higher weight in the model.

Based on preliminary findings we found that Pagel’s lambda was = 1 and running corPagel lead to problems of convergence. Therefore we ran the following models assuming Brownian Motion with corBrownian.

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')

#Create model dataframe for use in models
MCC.BAMM.model.data <- left_join(restricted.data %>% 
                                   dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
                                 MCC.BAMM.df %>% as.data.frame(), 
                                 by = "TipLabel")

saveRDS(MCC.BAMM.model.data, 'data/MCC.BAMM.model.data.rds')
#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.BAMM.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.BAMM.model.data) <- MCC.BAMM.model.data$TipLabel

#Run model for DR
MCC.BAMM.lambda <- gls(mean.lambda ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.lambda), #sqrt to account for overdispersedskewed variance distribution
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

#Run model for DR
MCC.BAMM.mu <- gls(mean.mu ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")


#Dredge the global MCC models

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.BAMM.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.MCC.lambda <- pdredge(MCC.BAMM.lambda, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.lambda, "data/dredged.MCC.lambda.rds")

dredged.MCC.mu <- pdredge(MCC.BAMM.mu, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.mu, "data/dredged.MCC.mu.rds")

Table S6: The dredged models for BAMM speciation and BAMM extinction both show the top model is one with no interactions, with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

dredged.MCC.lambda <- readRDS("data/dredged.MCC.lambda.rds")
dredged.MCC.mu <- readRDS("data/dredged.MCC.mu.rds")
kable(dredged.MCC.lambda, "html", caption = "BAMM-Speciation Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

kable(dredged.MCC.mu, "html", caption = "BAMM-Extinction Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

Given that we ran models weighted according to the variance of the response (\(\lambda_{BAMM}\) and \(\mu_{BAMM}\)), we checked to see whether a weighted model is favourable to an unweighted model using an anova to compare AIC values.

#Run the top model for the MCC
#Run model for ND
MCC.Lambda.top <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top, 'data/MCC.Lambda.top.rds')

MCC.Mu.top <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

saveRDS(MCC.Mu.top, 'data/MCC.Mu.top.rds')


#We can also see how the models look without the weightings: 
MCC.Lambda.top.unweighted <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top.unweighted, 'data/MCC.Lambda.top.unweighted.rds')


MCC.Mu.top.unweighted <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Mu.top.unweighted, 'data/MCC.Mu.top.unweighted.rds')

Table S7: Models with a set of weights for the response based on the inverse of the variance of the posterior distribution in \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) have much lower AIC values, indicating that the weighting scheme improves the fit of the model. In any case, we found that an unweighted model did not change qualitative results.

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.Lambda.top.unweighted <- readRDS('data/MCC.Lambda.top.unweighted.rds')
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.Mu.top.unweighted <- readRDS('data/MCC.Mu.top.unweighted.rds')

anova(MCC.Lambda.top, MCC.Lambda.top.unweighted) %>% pander(split.table = Inf)

anova(MCC.Mu.top, MCC.Mu.top.unweighted) %>% pander(split.table = Inf)

We ran the top (no interactions) model on the 100 trees, each with unique estimates of speciation and extinction from the BAMM runs.

#Read in the BAMM data for the 100 trees
BAMM.df <- readRDS('data/BAMM.df.rds')

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the BAMM for the 100 trees

BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
BAMM.model.data <- lapply(BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.lambda.pgls.models, "data/BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.mu.pgls.models, "data/BAMM.mu.pgls.models.rds")

BAMM.lambda.pgls.models <- readRDS("data/BAMM.lambda.pgls.models.rds")

BAMM.lambda.pgls.summary <- lapply(BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

#Now plot this for Lambda

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.lambda.summary <- data.frame(MCC.Lambda.top$coefficients,
                              confint(MCC.Lambda.top),
                              coef(summary(MCC.Lambda.top))[,2], #Std.Error
                              coef(summary(MCC.Lambda.top))[,3], #t-val
                              coef(summary(MCC.Lambda.top))[,4], #pval
                              MCC.Lambda.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))


parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

BAMM.lambda.pgls.summary <- bind_rows(BAMM.lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

  # for (x in unique(BAMM.lambda.pgls.summary$Parameter)[2:7]){
  #   filter(Parameter == x & between(Estimate, left = as.numeric(hpd.Lambda.top[1,x]), right = as.numeric(hpd.Lambda.top[2,x])))
  #   } 


remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

BAMM.lambda.pgls.summary.RO <- dcast(BAMM.lambda.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
BAMM.lambda.pgls.summary.RO$Estimate <- NULL
BAMM.lambda.pgls.summary.RO <- sapply(BAMM.lambda.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.lambda.pgls.summary.RO <-melt(BAMM.lambda.pgls.summary.RO) %>% na.omit()
BAMM.lambda.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.lambda.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.lambda.plot <- BAMM.lambda.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.lambda.plot <- BAMM.lambda.plot + geom_errorbarh(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("BAMM Speciation")


#Figure for extinction

BAMM.mu.pgls.models <- readRDS("data/BAMM.mu.pgls.models.rds")

BAMM.mu.pgls.summary <- lapply(BAMM.mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

BAMM.mu.pgls.summary <- bind_rows(BAMM.mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

#Now plot this for Lambda

MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.mu.summary <- data.frame(MCC.Mu.top$coefficients,
                              confint(MCC.Mu.top),
                              coef(summary(MCC.Mu.top))[,2], #Std.Error
                              coef(summary(MCC.Mu.top))[,3], #t-val
                              coef(summary(MCC.Mu.top))[,4], #pval
                              MCC.Mu.top[["modelStruct"]][["corStruct"]][1], #lambda value
                              model = "MCC_model") %>% 
  tibble::rownames_to_column() %>% 
  `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI", "SE", "tval", "pval" , "lambda" ,"model_name"))

#Get rid of outliers
BAMM.mu.pgls.summary.RO <- dcast(BAMM.mu.pgls.summary %>% select(Parameter,Estimate) %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter, value.var = "Estimate")
BAMM.mu.pgls.summary.RO$Estimate <- NULL
BAMM.mu.pgls.summary.RO <- sapply(BAMM.mu.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.mu.pgls.summary.RO <-melt(BAMM.mu.pgls.summary.RO) %>% na.omit()
BAMM.mu.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.mu.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.mu.plot <-BAMM.mu.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.mu.plot <- BAMM.mu.plot + geom_errorbarh(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,0.5,.5),"cm"))+
  ggtitle("BAMM Extinction")

grid.arrange(
  symmetrise_scale(DR.plot, "x"),
  symmetrise_scale(ND.plot, "x"),
  symmetrise_scale(BAMM.lambda.plot, "x"),
  symmetrise_scale(BAMM.mu.plot, "x"),
  nrow = 1
) 

MCC.DR.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates", split.table = Inf)

MCC.ND.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates", split.table = Inf)

MCC.lambda.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates", split.table = Inf)

MCC.mu.summary %>% select(-model_name) %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates", split.table = Inf)

hpd.DR.top <- list()
for(x in unique(DR.pgls.summary$Parameter)) {
hpd.DR.top[[x]] = hdi(DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select("Estimate"))
}
hpd.DR.top <- bind_rows(hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.DR.top, 'data/hpd.DR.top.rds')

#For ND
hpd.ND.top <- list()
for(x in unique(ND.pgls.summary$Parameter)) {
hpd.ND.top[[x]] = hdi(ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.ND.top <- bind_rows(hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 

saveRDS(hpd.ND.top, 'data/hpd.ND.top.rds')

hpd.Lambda.top <- list()
for(x in unique(BAMM.lambda.pgls.summary$Parameter)) {
hpd.Lambda.top[[x]] = hdi(BAMM.lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Lambda.top <- bind_rows(hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 


saveRDS(hpd.Lambda.top, 'data/hpd.Lambda.top.rds')

hpd.Mu.top <- list()
for(x in unique(BAMM.mu.pgls.summary$Parameter)) {
hpd.Mu.top[[x]] = hdi(BAMM.mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Mu.top <- bind_rows(hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.Mu.top, 'data/hpd.Mu.top.rds')

hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "DR HPD Interval")

hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "ND HPD Interval")

hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Speciation HPD Interval")

hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Extinction HPD Interval")

hpd.DR.top2 <- hpd.DR.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.DR.top2$HPD.Interval <- hpd.DR.top2$Upper - hpd.DR.top2$Lower
DR.intervals <- left_join(hpd.DR.top2, MCC.DR.summary, by = "Parameter")
DR.intervals$MCC.Interval <- DR.intervals$UCI - DR.intervals$LCI
DR.intervals$DR.IntervalRatio <- hpd.DR.top2$HPD.Interval/DR.intervals$MCC.Interval

hpd.ND.top2 <- hpd.ND.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.ND.top2$HPD.Interval <- hpd.ND.top2$Upper - hpd.ND.top2$Lower
ND.intervals <- left_join(hpd.ND.top2, MCC.ND.summary, by = "Parameter")
ND.intervals$MCC.Interval <- ND.intervals$UCI - ND.intervals$LCI
ND.intervals$ND.IntervalRatio <- hpd.ND.top2$HPD.Interval/ND.intervals$MCC.Interval


hpd.Lambda.top2 <- hpd.Lambda.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Lambda.top2$HPD.Interval <- hpd.Lambda.top2$Upper - hpd.Lambda.top2$Lower
Lambda.intervals <- left_join(hpd.Lambda.top2, MCC.lambda.summary, by = "Parameter")
Lambda.intervals$MCC.Interval <- Lambda.intervals$UCI - Lambda.intervals$LCI
Lambda.intervals$Lambda.IntervalRatio <- hpd.Lambda.top2$HPD.Interval/Lambda.intervals$MCC.Interval


hpd.Mu.top2 <- hpd.Mu.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Mu.top2$HPD.Interval <- hpd.Mu.top2$Upper - hpd.Mu.top2$Lower
Mu.intervals <- left_join(hpd.Mu.top2, MCC.mu.summary, by = "Parameter")
Mu.intervals$MCC.Interval <- Mu.intervals$UCI - Mu.intervals$LCI
Mu.intervals$Mu.IntervalRatio <- hpd.Mu.top2$HPD.Interval/Mu.intervals$MCC.Interval


plyr::join_all(list(DR.intervals %>% select(Parameter, DR.IntervalRatio),
          ND.intervals %>% select(Parameter, ND.IntervalRatio),
          Lambda.intervals %>% select(Parameter, Lambda.IntervalRatio),
          Mu.intervals %>% select(Parameter, Mu.IntervalRatio)), by = "Parameter", type = "left") %>% `colnames<-`(c("Parameter", "λDR.Interval.Ratio", "λND.Interval.Ratio", "λBAMM.Interval.Ratio", "μBAMM.Interval.Ratio")) %>%
  pander(digits = 3, split.table = Inf) 

Subsetted analysis with spectrophotometry data

Using the dataset from Armenta et al. (2008) we conducted the analysis on a subset of species for which dichromatism values were derived from spectrophotometry data. The measure of dichromatism is a difference measure between the sexes, based on a model of bird colour vision. To make this dataset comparable with the RGB measure of dichromatism, we use the absolute difference between the sexes; thereby making the scale from monochromatism to dichromatism rather than female colouration to male colouration.

#Read in Armenta data
Armenta.data <- read.csv('data/Armenta_2008.csv', stringsAsFactors = F)
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#A couple of rows have "-", remove from dataset
Armenta.data <- Armenta.data %>% dplyr::select(binomial, Colour.discriminability) %>% mutate(Colour.discriminability = replace(Colour.discriminability, Colour.discriminability == "-", "NA")) %>% filter(Colour.discriminability != "NA")
Armenta.data$Colour.discriminability <- Armenta.data$Colour.discriminability %>% as.numeric()

MCC.Armenta.model.data <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), Armenta.data, by = "binomial")

MCC.Armenta.model.data <- inner_join(MCC.Armenta.model.data, MCC.BAMM.df, by = "TipLabel") %>% filter(Colour.discriminability != "NA")

#Plot the correlation
p1 <- MCC.Armenta.model.data %>% ggplot(aes(x = SDi, y = Colour.discriminability))+
  geom_point()+
  geom_smooth(method = 'loess')+
  theme_minimal()

grid.newpage()
ggExtra::ggMarginal(
  p = p1,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

Figure S10: There is a correlation between the RGB measures and the Spectophotometry measures. The RGB measures seem to be more noisy around the lower values.

R-squared and correlation for relationship between RGB and Spectrophotometry Data:

data_frame(
R2 = summary(lm(SDi ~ Colour.discriminability,
   data = MCC.Armenta.model.data))$r.squared,
r = cor(MCC.Armenta.model.data$SDi, MCC.Armenta.model.data$Colour.discriminability)) %>% pander()

Run PGLS model for the data

#Prune tree
pruned.MCC.Armenta.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.Armenta.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.Armenta.model.data) <- MCC.Armenta.model.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Armenta.global <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Armenta.global <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Armenta.global <- gls(mean.lambda ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Armenta.global <- gls(mean.mu ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.Armenta.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Armenta.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Armenta.dredged.ND.model <- pdredge(MCC.ND.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.ND.model, 'data/Armenta.dredged.ND.model.rds')

Armenta.dredged.DR.model <- pdredge(MCC.DR.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.DR.model, 'data/Armenta.dredged.DR.model.rds')

Armenta.dredged.spec.model <- pdredge(MCC.Lambda.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.spec.model, 'data/Armenta.dredged.spec.model.rds')

Armenta.dredged.extinct.model <- pdredge(MCC.Extinction.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.extinct.model, 'data/Armenta.dredged.extinct.model.rds')

Table S11: The dredged models using a subset of species for which sexual dichromatism was measured using spectrophotometry. All four show the top model is one with no interactions, with \(\delta AICc > 4\).

Armenta.dredged.ND.model <- readRDS("data/Armenta.dredged.ND.model.rds")
Armenta.dredged.DR.model <- readRDS("data/Armenta.dredged.DR.model.rds")
Armenta.dredged.spec.model <- readRDS("data/Armenta.dredged.spec.model.rds")
Armenta.dredged.extinct.model <- readRDS("data/Armenta.dredged.extinct.model.rds")


kable(Armenta.dredged.ND.model, "html", caption = "ND Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

kable(Armenta.dredged.DR.model, "html", caption = "DR Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

kable(Armenta.dredged.spec.model, "html", caption = "BAMM Speciation Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

kable(Armenta.dredged.extinct.model, "html", caption = "BAMM Extinction Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

Run the top model and 100 models for each one:

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
Armenta.MCC.DR.top <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.DR.top, 'data/Armenta.MCC.DR.top.rds')

#Run model for ND
Armenta.MCC.ND.top <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.ND.top, 'data/Armenta.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Armenta.data.noMCC <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), Armenta.data, by = "binomial") %>% filter(Colour.discriminability != "NA")
#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Armenta.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.DR.model.data <- lapply(Armenta.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Armenta.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(MCC.Armenta.model.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Armenta.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(1, phy = y, fixed = FALSE), 
    data = x, 
    method = "REML")
}, x = Armenta.DR.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.DR.pgls.models, "data/Armenta.DR.pgls.models.rds")

#Now for Node Density:
Armenta.ND.model.data <- lapply(nd.list, function(x) {
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.ND.model.data <- lapply(Armenta.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Armenta.ND.model.data, y = Armenta.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.ND.pgls.models, "data/Armenta.ND.pgls.models.rds")

#Run the BAMM models

Armenta.MCC.Lambda.top <- gls(mean.lambda ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Lambda.top, 'data/Armenta.MCC.Lambda.top.rds')

Armenta.MCC.Mu.top <- gls(mean.mu ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Mu.top, 'data/Armenta.MCC.Mu.top.rds')

Armenta.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.BAMM.model.data <- lapply(Armenta.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.lambda.pgls.models, "data/Armenta.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Armenta.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ Colour.discriminability
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.mu.pgls.models, "data/Armenta.BAMM.mu.pgls.models.rds")

Using the MCC models as well as the 100 PGLS models we can generate the plots similar to Figure S8:

Armenta.DR.pgls.models <- readRDS("data/Armenta.DR.pgls.models.rds")
Armenta.MCC.DR.top <- readRDS('data/Armenta.MCC.DR.top.rds')

Armenta.DR.pgls.summary <- lapply(Armenta.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.DR.summary <- data.frame(Armenta.MCC.DR.top$coefficients, confint(Armenta.MCC.DR.top)) %>% tibble::rownames_to_column()

Armenta.DR.pgls.summary <- bind_rows(Armenta.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Colour.discriminability` = "Sexual Dichromatism"
                    )

Armenta.DR.pgls.summary$Parameter = factor(Armenta.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.DR.summary$Parameter = factor(Armenta.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.DR.plot <-Armenta.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.DR.plot <- Armenta.DR.plot + geom_errorbarh(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta DR")

#For ND

Armenta.ND.pgls.models <- readRDS("data/Armenta.ND.pgls.models.rds")
Armenta.MCC.ND.top <- readRDS('data/Armenta.MCC.ND.top.rds')

Armenta.ND.pgls.summary <- lapply(Armenta.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.ND.summary <- data.frame(Armenta.MCC.ND.top$coefficients, confint(Armenta.MCC.ND.top)) %>% tibble::rownames_to_column()

Armenta.ND.pgls.summary <- bind_rows(Armenta.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.ND.pgls.summary$Parameter = factor(Armenta.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.ND.summary$Parameter = factor(Armenta.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.ND.plot <-Armenta.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.ND.plot <- Armenta.ND.plot + geom_errorbarh(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta ND")

#For Lambda
Armenta.BAMM.lambda.pgls.models <- readRDS("data/Armenta.BAMM.lambda.pgls.models.rds")
Armenta.MCC.Lambda.top <- readRDS('data/Armenta.MCC.Lambda.top.rds')

Armenta.Lambda.pgls.summary <- lapply(Armenta.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Lambda.summary <- data.frame(Armenta.MCC.Lambda.top$coefficients, confint(Armenta.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Armenta.Lambda.pgls.summary <- bind_rows(Armenta.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Lambda.pgls.summary$Parameter = factor(Armenta.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Lambda.summary$Parameter = factor(Armenta.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Lambda.plot <-Armenta.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Lambda.plot <- Armenta.Lambda.plot + geom_errorbarh(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Speciaton")

#For Mu

Armenta.BAMM.Mu.pgls.models <- readRDS("data/Armenta.BAMM.mu.pgls.models.rds")
Armenta.MCC.Mu.top <- readRDS('data/Armenta.MCC.Mu.top.rds')

Armenta.Mu.pgls.summary <- lapply(Armenta.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Mu.summary <- data.frame(Armenta.MCC.Mu.top$coefficients, confint(Armenta.MCC.Mu.top)) %>% tibble::rownames_to_column()

Armenta.Mu.pgls.summary <- bind_rows(Armenta.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Mu.pgls.summary$Parameter = factor(Armenta.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Mu.summary$Parameter = factor(Armenta.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Mu.plot <-Armenta.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Mu.plot <- Armenta.Mu.plot + geom_errorbarh(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Extinction")

Plot the results

grid.arrange(symmetrise_scale(Armenta.DR.plot, "x"),
             symmetrise_scale(Armenta.ND.plot, "x"),
             symmetrise_scale(Armenta.Lambda.plot, "x"), 
             symmetrise_scale(Armenta.Mu.plot, "x"), 
             nrow = 1)

Figure S11: Measures of sexual dichromatism using spectrophotometry do not change the main patterns. This dataset is a subset of the complete dataset (n = 581), thus drawing conclusions for the other predictors (e.g. borderline long term temperature variation and spatial temperature variation) potentially risks type I error. Additionally, these environmental predictors were not measured differently in this analysis compared to the analysis using RGB measures of dichromatism for all species so we refrain from drawing conclusions on the effects of the environmental predictors based on this smaller dataset.

Table S12: MCC estimates from the above plot are presented as numerical values below.

Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates from Spec measures")

Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates from Spec measures")

Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates from Spec measures")

Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates from Spec measures")

Table S13: HPD intervals calculated for the above figure (i.e. models using Spectrophotometry values of sexual dichromatism). These intervals do not take into account the variance associated with each estimate and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that are in the same direction (+ or -). These intervals are calculated in the same way as in Table S9.

Armenta.hpd.DR.top <- list()
Armenta.DR.pgls.summary <- na.omit(Armenta.DR.pgls.summary)
for(x in unique(Armenta.DR.pgls.summary$Parameter)) {
Armenta.hpd.DR.top[[x]] = hdi(Armenta.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.DR.top <- bind_rows(Armenta.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Armenta.hpd.DR.top, 'data/Armenta.hpd.DR.top.rds')

#For ND
Armenta.hpd.ND.top <- list()
Armenta.ND.pgls.summary <- na.omit(Armenta.ND.pgls.summary)
for(x in unique(Armenta.ND.pgls.summary$Parameter)) {
Armenta.hpd.ND.top[[x]] = hdi(Armenta.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.ND.top <- bind_rows(Armenta.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.ND.top, 'data/Armenta.hpd.ND.top.rds')

Armenta.hpd.Lambda.top <- list()
Armenta.Lambda.pgls.summary <- na.omit(Armenta.Lambda.pgls.summary)
for(x in unique(Armenta.Lambda.pgls.summary$Parameter)) {
Armenta.hpd.Lambda.top[[x]] = hdi(Armenta.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Lambda.top <- bind_rows(Armenta.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Armenta.hpd.Lambda.top, 'data/Armenta.hpd.Lambda.top.rds')

Armenta.hpd.Mu.top <- list()
Armenta.Mu.pgls.summary <- na.omit(Armenta.Mu.pgls.summary)
for(x in unique(Armenta.Mu.pgls.summary$Parameter)) {
Armenta.hpd.Mu.top[[x]] = hdi(Armenta.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Mu.top <- bind_rows(Armenta.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.Mu.top, 'data/Armenta.hpd.Mu.top.rds')

Armenta.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta DR HPD Interval")

Armenta.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta ND HPD Interval")

Armenta.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Speciation HPD Interval")

Armenta.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Extinction HPD Interval")

Analysis using male-biased measure of sexual selection

Sexual dichromatism is expected to be a good measure of sexual selection strength in birds. However, the relationship between sexual dichromatism and male-biased measures of sexual selection (social mating system, sexual size dimorphism and paternal care) is expected to be relatively noisy given the precision of the measurement used (Dale et al. 2015). Here we use a dataset of sexual selection for 2,465 species of the 5,812 species with sexual dichromatism scores from RGB measures. This male-biased sexual selection score is based on three components, combined in a phylogenetic PCA (ppca) with the following loadings:

• sexual size dimorphism: 0.37
• social polygyny: 0.57
• paternal care: -0.57

As such, species with high values for this score are expected to have high male-biased sexual selection (high dimorphism, high polygyny and low paternal care). PGLS models using this dataset were run using the same process as above.

SS.subset <- plumage.scores %>% filter(Sexual_selection_ppca != "NA")
SS.subset <- inner_join(SS.subset, restricted.data %>% dplyr::select(TipLabel, binomial, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), by = "TipLabel")
SS.subset <- inner_join(SS.subset, MCC.BAMM.df, by = "TipLabel")

SS.subset %>% ggplot(aes(x = SDi, y = Sexual_selection_ppca))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_minimal()

#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel

Figure S12: The relationship between absolute sexual dichromatism and male-baised sexual selection. Although positively correlated, the distribution is very noisy with some species having high sexual dichromatism and low male-biased sexual selection and vice versa.

data_frame(
R2 = summary(lm(SDi ~ Sexual_selection_ppca,
   data = SS.subset))$r.squared,
r = cor(SS.subset$SDi, SS.subset$Sexual_selection_ppca)) %>% pander(digits = 2)

#Prune tree
pruned.MCC.Subset.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(SS.subset$TipLabel, MCC.passerine$tip.label)])

saveRDS(pruned.MCC.Subset.tree, 'data/pruned.MCC.Subset.tree.rds')

#Run a corPagel model to estimate lambda for DR
MCC.DR.Subset.global <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Subset.global <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Subset.global <- gls(mean.lambda ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variatioColour.discriminability
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Subset.global <- gls(mean.mu ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("SS.subset"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Subset.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Subset.dredged.ND.model <- pdredge(MCC.ND.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.ND.model, 'data/Subset.dredged.ND.model.rds')

Subset.dredged.DR.model <- pdredge(MCC.DR.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.DR.model, 'data/Subset.dredged.DR.model.rds')

Subset.dredged.spec.model <- pdredge(MCC.Lambda.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.spec.model.rds')

Subset.dredged.extinct.model <- pdredge(MCC.Extinction.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.extinct.model.rds')

Table S14: All top models (\(\delta AICc > 4\)) using male-biased sexual selection measures do not contain interactions between sexual selection and the environmental variables. Thus interaction terms were not included in further analysis.

Subset.dredged.ND.model <- readRDS("data/Subset.dredged.ND.model.rds")
Subset.dredged.DR.model <- readRDS("data/Subset.dredged.DR.model.rds")
Subset.dredged.spec.model <- readRDS("data/Subset.dredged.spec.model.rds")
Subset.dredged.extinct.model <- readRDS("data/Subset.dredged.extinct.model.rds")


kable(Subset.dredged.ND.model, "html", caption = "ND Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")