#===========================================================================# # # Allometry reveals trade-offs between Bergmann’s and Allen’s rules, # and different avian adaptive strategies for thermoregulation # # Arkadiusz Frohlich frohlich@iop.krakow.pl # #============================================================== # data preparation and choice of evolutionary models # working directory setwd("D:/") # load database avonet <- read_excel("AVONET Supplementary dataset 1.xlsx", sheet = "AVONET3_BirdTree") %>% # flie available at https://doi.org/10.1111/ele.13898 with(data.frame(species = Species3 %>% loopsub(" ", "_"), bodyMass = Mass, beakLengthCulmen = Beak.Length_Culmen, tarsusLength = Tarsus.Length, migration. = Migration %>% loopsub(1:3, c("Resident", "PartialMigrant", "FullMigrant")))) rownames(avonet) <- avonet$species avonet %<>% within(rm(species)) birdTherm <- read_excel("birdTherm.xlsx", sheet = "data") %>% # flie available at https://doi.org/10.5061/dryad.9ghx3ffn7 as.data.frame rownames(birdTherm) <- birdTherm$species dat <- mergeByRownames(birdTherm, avonet) # load phylogenies phylo <- read.nexus("AllBirdsEricson1Summary.tre") # this is consensus phylogenetic tree made from phylos (see below) phylo %<>% prune(dat$species) phylos <- read.tree("AllBirdsEricson100.tre") # this is random set of 100 phylogenies randomply selected from birdtree.org phylos %<>% prune(dat$species) names(phylos) <- paste0("phy", 1:length(phylos)) # define clusters of variables and abbreviations Temp <- c("absLatitude", "minTmin", "minTavg", "avgTavg", "maxTavg", "maxTmax") TempFullNames <- c("absolute latitude (°)", "min temp all months (°C)", "avg temp coldest month (°C)", "avg temp all months (°C)", "avg temp hottest month (°C)", "max temp all months (°C)") Morph <- Morph0 <- c("BodyMass", "BeakLengthCulmen", "TarsusLength") MorphFullNames <- c("body mass", "beak length", "tarsus length") logMorph <- paste0("log", Morph) logMorphFullNames <- paste("log", MorphFullNames) relMorph <- c("relBeakLength", "relTarsusLength") relMorphFullNames <- paste("relative", MorphFullNames[2:3]) Morph <- c(logMorph[1], relMorph) MorphFullNames <- c(logMorphFullNames[1], relMorphFullNames) # migration factor dat$migration. %<>% factor(c("Resident", "PartialMigrant", "FullMigrant")) # normalize the phenotype variables dat$logBodyMass <- dat$bodyMass %>% log dat$logBeakLengthCulmen <- dat$beakLengthCulmen %>% log dat$logTarsusLength <- dat$tarsusLength %>% log # tgransformations for negative skew in temperature variables: # https://www.datanovia.com/en/lessons/transform-data-to-normal-distribution-in-r/ tran0 <- function(x) sqrt(x) tran1 <- function(x) sqrt(max(x+1) - x) * -1 tran2 <- function(x) log10(max(x+1) - x) * -1 tran3 <- function(x) 1/(max(x+1) - x) * -1 tranList <- list(donothing = donothing, tran0 = tran0, tran1 = tran1, tran2 = tran2, tran3 = tran3) tranList %>% lapply2(dat$absLatitude) %>% multiHist dat$sqrtabsLatitude <- dat$absLatitude %>% tran0 %>% scale %>% as.numeric tranList %>% lapply2(dat$minTmin) %>% multiHist dat$tr2minTmin <- dat$minTmin %>% tran2 %>% scale %>% as.numeric tranList %>% lapply2(dat$minTavg) %>% multiHist dat$tr2minTavg <- dat$minTavg %>% tran2 %>% scale %>% as.numeric tranList %>% lapply2(dat$avgTavg) %>% multiHist dat$tr2avgTavg <- dat$avgTavg %>% tran2 %>% scale %>% as.numeric tranList %>% lapply2(dat$maxTavg) %>% multiHist dat$tr1maxTavg <- dat$maxTavg %>% tran1 %>% scale %>% as.numeric tranList %>% lapply2(dat$maxTmax) %>% multiHist dat$tr1maxTmax <- dat$maxTmax %>% tran1 %>% scale %>% as.numeric trTemp <- c("sqrtabsLatitude", "tr2minTmin", "tr2minTavg", "tr2avgTavg", "tr1maxTavg", "tr1maxTmax") trTempFullNames <- c("absolute latitude", "min temp all months", "avg temp coldest month", "avg temp all months", "avg temp hottest month", "max temp all months") # dictionary to translate variable names dict <- data.frame(n1 = c(trTemp, Temp, relMorph, logMorph), n2 = c(trTempFullNames, TempFullNames, relMorphFullNames, logMorphFullNames)) dict$n3 <- dict$n2 %>% loopsub(c(" \\(°C)", " \\(°)"), "") %>% loopsub(c(" mass", " length", "log ", "relative "), "") dict %<>% rbind(data.frame(n1 = c("migration.", "Resident", "PartialMigrant", "FullMigrant"), n2 = c("", "residents", "partial migrants", "full migrants"), n3 = c("", "residents", "partial migrants", "full migrants"))) #===================================================== # calculate phenotype derivatives # residuals mpackAllometry <- data.frame(response = c("logBeakLengthCulmen", "logTarsusLength")) %>% merge(data.frame(fixed = "logBodyMass")) %>% nameModels("Allometry") %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(ncores = 2, model = "lambda") mpackAllometry %>% lapply(plotSimpleTrend, pointColor = "grey") %>% ggCollage(ncol = 1) relBeakLength <- mpackAllometry[["Allometry_resp[logBeakLengthCulmen]_fix[logBodyMass]_"]] %>% residuals dat$relBeakLength <- relBeakLength[dat$species] relTarsusLength <- mpackAllometry[["Allometry_resp[logTarsusLength]_fix[logBodyMass]_"]] %>% residuals dat$relTarsusLength <- relTarsusLength[dat$species] # ratios ratios <- c("logBodyMass", "logRatioBeakBody", "logRatioTarsusBody") dat$ratioBeakBody <- dat$beakLengthCulmen/dat$bodyMass dat$ratioBeakBody %>% ggHist dat$logRatioBeakBody <- dat$ratioBeakBody %>% log dat$logRatioBeakBody %>% ggHist ggplot() + geom_point(data = dat, aes(x = logRatioBeakBody, y = logBodyMass)) + theme_bw() cor(dat$logBodyMass, dat$logRatioBeakBody) ggplot() + geom_point(data = dat, aes(x = logRatioBeakBody, y = logBeakLengthCulmen)) cor(dat$logBeakLengthCulmen, dat$logRatioBeakBody) dat$ratioTarsusBody <- dat$tarsusLength/dat$bodyMass dat$ratioTarsusBody %>% ggHist dat$logRatioTarsusBody <- dat$ratioTarsusBody %>% log dat$logRatioTarsusBody %>% ggHist ggplot() + geom_point(data = dat, aes(x = logRatioTarsusBody, y = logBodyMass)) cor(dat$logBodyMass, dat$logRatioTarsusBody) ggplot() + geom_point(data = dat, aes(x = logRatioTarsusBody, y = logTarsusLength)) cor(dat$logTarsusLength, dat$logRatioTarsusBody) ggarrange(ggplot() + geom_point(data = dat, aes(x = log(beakLengthCulmen/bodyMass), y = logBodyMass), alpha = .5, shape = 16) + theme_bw(), ggplot() + geom_point(data = dat, aes(x = log(tarsusLength/bodyMass), y = logBodyMass), alpha = .5, shape = 16) + theme_bw()) # PCA PCs <- paste0("PC", 1:3) PhenotypePCA <- phyl.pca(phylo, dat[logMorph]) dat %<>% mergeByRownames(PhenotypePCA$S) # define functions as a checkpoints to investigate interactions between two continuous variables #===================================================== q1 <- function(x) x %>% quantile(probs = .25) %>% as.numeric q3 <- function(x) x %>% quantile(probs = .75) %>% as.numeric contFuns <- list(min = min, q1 = q1, median = median, q3 = q3, max = max) modFuns1 <- rep(list(contFuns), 15) names(modFuns1) <- c(Temp, trTemp, Morph) #===================================================== # define simple color gradients applied in entire article myProbs <- c(0,.25,.5,.75,1) meanValue <- myProbs[3] palBody <- ggarrange(rgb(myProbs, meanValue, meanValue) %>% showColors, rgb(myProbs, meanValue, meanValue) %>% lighten(c(1.2,1.1,1,.9,.8)) %>% showColors, rgb2(myProbs, meanValue, meanValue, lighting = NULL) %>% showColors, rgb2(myProbs, meanValue, meanValue, lighting = .25) %>% showColors, rgb2(myProbs, meanValue, meanValue, lighting = .5) %>% showColors, rgb2(myProbs, meanValue, meanValue, lighting = .75) %>% showColors, rgb2(myProbs, meanValue, meanValue, lighting = 1) %>% showColors, nrow = 1) palBeak <- ggarrange(rgb(meanValue, myProbs,meanValue) %>% showColors, rgb(meanValue, myProbs,meanValue) %>% lighten(c(1.2,1.1,1,.9,.8)) %>% showColors, rgb2(meanValue, myProbs,meanValue, lighting = NULL) %>% showColors, rgb2(meanValue, myProbs,meanValue, lighting = .25) %>% showColors, rgb2(meanValue, myProbs,meanValue, lighting = .5) %>% showColors, rgb2(meanValue, myProbs,meanValue, lighting = .75) %>% showColors, rgb2(meanValue, myProbs,meanValue, lighting = 1) %>% showColors, nrow = 1) palTarsus <- ggarrange(rgb(meanValue,meanValue, myProbs) %>% showColors, rgb(meanValue,meanValue, myProbs) %>% lighten(c(1.2,1.1,1,.9,.8)) %>% showColors, rgb2(meanValue,meanValue, myProbs, lighting = NULL) %>% showColors, rgb2(meanValue,meanValue, myProbs, lighting = .25) %>% showColors, rgb2(meanValue,meanValue, myProbs, lighting = .5) %>% showColors, rgb2(meanValue,meanValue, myProbs, lighting = .75) %>% showColors, rgb2(meanValue,meanValue, myProbs, lighting = 1) %>% showColors, nrow = 1) ggarrange(palBody, palBeak, palTarsus, ncol = 1) bodyGradient <- rgb(myProbs, meanValue, meanValue) %>% lighten(c(1.2,1.1,1,.9,.8)) %>% as.list beakGradient <- rgb(meanValue, myProbs, meanValue) %>% lighten(c(1.2,1.1,1,.9,.8)) %>% as.list tarsusGradient <- rgb(meanValue, meanValue, myProbs) %>% lighten(c(1.2,1.1,1,.9,.8)) %>% as.list temperatureGradient <- colorGradient(x = 1:5, low = "#377cb8ff", high = "#e4181cff") names(bodyGradient) <- names(beakGradient) <- names(tarsusGradient) <- names(temperatureGradient) <- c("min", "q1", "median", "q3", "max") gradients <- rep(list(temperatureGradient), 12) %>% append(list(bodyGradient, beakGradient, tarsusGradient)) names(gradients) <- c(Temp, trTemp, Morph) dat$color <- names(gradients) %>% autoNames %>% lapply(function(i){ if (i %in% relMorph) f <- function(x) 0 else f = median colorGradient2(x = dat[,i], low = gradients[[i]][["min"]], high = gradients[[i]][["max"]], middle = gradients[[i]][["median"]], middleFun = f) } ) %>% as.data.frame # show legend for color gradients phenotypeGradient3d <- data.frame(body = myProbs) %>% merge(data.frame(beak = myProbs)) %>% merge(data.frame(tarsus = myProbs)) %>% within(color <- rgb(body, beak, tarsus)) #library(gg3D) #ggplot(data = phenotypeGradient3d, aes(x = body, y = beak, z = tarsus, color = color)) + # axes_3D(theta = 180, phi = 180) + # stat_3D(theta = 180, phi = 180) + # theme(legend.position = "none", panel.background = element_blank()) #===================================================== # scale all variables dat.scaled <- dat.nonscaled <- dat dat.scaled[,c(logMorph, relMorph, Temp, trTemp)] %<>% scale #===================================================== # explore variation in species traits closestTo <- function(x, to) (abs(to-x)) == min(abs(to-x)) exploreData <- function(trait, fun, trait2 = NULL) dat[closestTo(dat[,trait], fun(dat[,trait])), c(trait, trait2, "species")] exploreData("bodyMass", min) exploreData("bodyMass", max) exploreData("bodyMass", median, "rangeArea") %>% ordby("rangeArea", rev) exploreData("bodyMass", q1, "rangeArea") %>% ordby("rangeArea", rev) exploreData("bodyMass", q3, "rangeArea") %>% ordby("rangeArea", rev) exploreData("relBeakLength", min) exploreData("relBeakLength", max) exploreData("relBeakLength", median, "rangeArea") %>% ordby("rangeArea", rev) exploreData("relTarsusLength", min) exploreData("relTarsusLength", max) exploreData("relTarsusLength", median, "rangeArea") %>% ordby("rangeArea", rev) exploreData("rangeArea", max) exploreData("maxTmax", min) exploreData("maxTmax", median, "rangeArea") %>% ordby("rangeArea", rev) exploreData("maxTmax", max) exploreData("minTmin", min) exploreData("minTmin", median, "rangeArea") %>% ordby("rangeArea", rev) exploreData("minTmin", max) #============================================================== # define models to select phylogny settings in phylolm models BergmannModelTab <- data.frame(response = logMorph[1]) %>% merge(data.frame(fixed = c("1", Temp))) %>% nameModels("Bergmann") AllenModelTab <- data.frame(response = c("logBeakLengthCulmen", "logTarsusLength")) %>% merge(data.frame(fixed = c("1", "logBodyMass", paste0(logMorph[1], " + ", Temp), paste0(logMorph[1], " * ", Temp)))) %>% nameModels("Allen") TemperatureModelTab <- data.frame(response = trTemp) %>% merge(data.frame(fixed = c("1", "logBodyMass", relMorph, paste0(logMorph[1], " + ", relMorph), paste0(logMorph[1], " * ", relMorph), paste0(logMorph[1], " + ", relMorph[1]), paste0(logMorph[1], " * ", relMorph[1]), paste0(logMorph[1], " + ", relMorph[2]), paste0(logMorph[1], " * ", relMorph[2]), paste0(relMorph[1], " + ", relMorph[2]), paste0(relMorph[1], " * ", relMorph[2]), paste0(logMorph[1], " + ", relMorph[1], " + ", relMorph[2]), paste0(logMorph[1], " * (", relMorph[1], " + ", relMorph[2], ")"), paste0(relMorph[1]," * (logBodyMass + ", relMorph[2], ")"), paste0(relMorph[2]," * (logBodyMass + ", relMorph[1], ")"), paste0(logMorph[1], " + ", relMorph[1], " * ", relMorph[2]), paste0(logMorph[1], " * ", relMorph[1], " + ", relMorph[2]), paste0(logMorph[1], " * ", relMorph[2], " + ", relMorph[1]), paste0(logMorph[1], " * ", relMorph[1], " * ", relMorph[2])))) %>% nameModels("Temperature") # compute brownian motion models mpackBM <- rbind(BergmannModelTab, AllenModelTab, TemperatureModelTab) %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(model = "BM", full.matrix = T) save(mpackBM, file = "mpackBM.rda") # compute OU fixed root models mpackOUfixroot <- rbind(BergmannModelTab, AllenModelTab, TemperatureModelTab) %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(model = "OUfixedRoot", full.matrix = T) save(mpackOUfixroot, file = "mpackOUfixroot.rda") # compute OU random root models mpackOUrandroot <- rbind(BergmannModelTab, AllenModelTab, TemperatureModelTab) %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(model = "OUrandomRoot", full.matrix = T) save(mpackOUrandroot, file = "mpackOUrandroot.rda") # compute lambda models mpackLambda <- rbind(BergmannModelTab, AllenModelTab, TemperatureModelTab) %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(model = "lambda", full.matrix = T) save(mpackLambda, file = "mpackLambda.rda") # compute OU random root models mpackKappa <- rbind(BergmannModelTab, AllenModelTab, TemperatureModelTab) %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(model = "kappa", full.matrix = T) save(mpackKappa, file = "mpackKappa.rda") # compute OU random root models mpackDelta <- rbind(BergmannModelTab, AllenModelTab, TemperatureModelTab) %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(model = "delta", full.matrix = T) save(mpackDelta, file = "mpackDelta.rda") # compute OU random root models mpackEB <- rbind(BergmannModelTab, AllenModelTab, TemperatureModelTab) %>% sliceTab %>% lapply(attachDataPhylo) %>% parPhyLm(model = "EB", full.matrix = T) save(mpackEB, file = "mpackEB.rda") # check consistency of different phylogenetic models compPlot <- function(resp) { df <- rbind(mpackBM %>% selectPhyloLmm() %>% spreadResponses(resp), mpackOUrandroot %>% selectPhyloLmm() %>% spreadResponses(resp), mpackOUfixroot %>% selectPhyloLmm() %>% spreadResponses(resp), mpackLambda %>% selectPhyloLmm() %>% spreadResponses(resp), mpackDelta %>% selectPhyloLmm() %>% spreadResponses(resp), mpackKappa %>% selectPhyloLmm() %>% spreadResponses(resp), mpackEB %>% selectPhyloLmm() %>% spreadResponses(resp)) df <- df[!(grepl("Length]", rownames(df)) & grepl("\\[1]", rownames(df))),] df$model %<>% factor(levels = unique(df$model)) p <- ggplot(df, aes(x = delta, y = model, group = fixed, color = fixed, shape = fixed)) + geom_point() + geom_line() return(p)} logMorph %>% rev %>% autoNames %>% lapply(compPlot) %>% entitle %>% ggCollage(common.legend = T)