## METADATA =============================================================== ## Hagen et al. Earth history events shaped the evolution of uneven biodiversity across tropical moist forests. PNAS ## Description:Fiting GLS models of species richness and environment relationships ## ## R version: 4.0.0 for Windows ## Date: 2020-07-14 16:20:32 ## License: GPL3 ## Author: Oskar Hagen (oskar@hagen.bio) and Alexander Skeels (alexander.skeels@gmail.com) ##=======================================================================## data_df <- read.csv(file="Dataset_S2.csv", header=T) ###~~~~~~~~~~~~~~~~~ Start Figure S1 # colours color_scale <- c("#FAA22C","#af8dc3","#007E6A") # plot PCA pca_plot <- ggplot(data_df, aes(PC1, PC2, colour = region, fill=region)) + geom_point(alpha=0.6) + stat_ellipse(geom="polygon",level=0.95,alpha=0.2, size=2)+ theme_light()+ scale_fill_manual(values = color_scale)+ scale_colour_manual(values = color_scale) # plot density of env. variables. MAP_p <- ggplot(data_df, aes(MAP, colour = region, fill=region)) + geom_density(alpha=0.6) + theme_light()+ scale_fill_manual(values = color_scale)+ scale_colour_manual(values = color_scale) MAT_p <- ggplot(data_df, aes(MAT, colour = region, fill=region)) + geom_density(alpha=0.6) + theme_light()+ scale_fill_manual(values = color_scale)+ scale_colour_manual(values = color_scale) PET_p <- ggplot(data_df, aes(PET, colour = region, fill=region)) + geom_density(alpha=0.6) + theme_light()+ scale_fill_manual(values = color_scale)+ scale_colour_manual(values = color_scale) pca_plot grid.arrange(PET_p, MAT_p, MAP_p) ###~~~~~~~~~~~~~~~~~ End Figure S1 # subset regional data afro_df <- data_df[which(data_df$region %in% "afro")] neo_df <- data_df[which(data_df$region %in% "neo")] indo_df <- data_df[which(data_df$region %in% "indo")] # scale units so they are comparble across variables afro_df[, c(3:5)] <- scale(afro_df[, c(3:5)]) neo_df[, c(3:5)] <- scale(neo_df[, c(3:5)]) indo_df[, c(3:5)] <- scale(indo_df[, c(3:5)]) # amphibian ~ MAT m_a_b1_i <- gls(log(amph+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_a_b1_n <- gls(log(amph+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_a_b1_a <- gls(log(amph+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # mammal ~ MAT m_m_b1_i <- gls(log(mamm+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_m_b1_n <- gls(log(mamm+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_m_b1_a <- gls(log(mamm+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # bird ~ MAT m_b_b1_i <- gls(log(bird+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_b_b1_n <- gls(log(bird+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_b_b1_a <- gls(log(bird+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # squamate ~ MAT m_s_b1_i <- gls(log(squa+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_s_b1_n <- gls(log(squa+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_s_b1_a <- gls(log(squa+1) ~ MAT, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # amphibian ~ MAP m_a_b12_i <- gls(log(amph+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_a_b12_n <- gls(log(amph+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_a_b12_a <- gls(log(amph+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # mammal ~ MAP m_m_b12_i <- gls(log(mamm+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_m_b12_n <- gls(log(mamm+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_m_b12_a <- gls(log(mamm+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # bird ~ MAP m_b_b12_i <- gls(log(bird+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_b_b12_n <- gls(log(bird+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_b_b12_a <- gls(log(bird+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # squamate ~ MAP m_s_b12_i <- gls(log(squa+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_s_b12_n <- gls(log(squa+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_s_b12_a <- gls(log(squa+1) ~ MAP, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # amphibian ~ PET m_a_p_i <- gls(log(amph+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_a_p_n <- gls(log(amph+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_a_p_a <- gls(log(amph+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # mammal ~ PET m_m_p_i <- gls(log(mamm+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_m_p_n <- gls(log(mamm+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_m_p_a <- gls(log(mamm+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # bird ~ PET m_b_p_i <- gls(log(bird+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_b_p_n <- gls(log(bird+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_b_p_a <- gls(log(bird+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) # squamte ~ PET m_s_p_i <- gls(log(squa+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = indo_df) m_s_p_n <- gls(log(squa+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = neo_df) m_s_p_a <- gls(log(squa+1) ~ PET, correlation = corGaus(form = ~x + y, nugget = TRUE), data = afro_df) gls_model_results <- list( m_a_b1_i = m_a_b1_i, m_a_b1_n = m_a_b1_n, m_a_b1_a = m_a_b1_a, m_m_b1_i = m_m_b1_i, m_m_b1_n = m_m_b1_n, m_m_b1_a = m_m_b1_a, m_b_b1_i = m_b_b1_i, m_b_b1_n = m_b_b1_n, m_b_b1_a = m_b_b1_a, m_s_b1_i = m_s_b1_i, m_s_b1_n = m_s_b1_n, m_s_b1_a = m_s_b1_a, m_a_b12_i = m_a_b12_i, m_a_b12_n = m_a_b12_n, m_a_b12_a = m_a_b12_a, m_m_b12_i = m_m_b12_i, m_m_b12_n = m_m_b12_n, m_m_b12_a = m_m_b12_a, m_b_b12_i = m_b_b12_i, m_b_b12_n = m_b_b12_n, m_b_b12_a = m_b_b12_a, m_s_b12_i = m_s_b12_i, m_s_b12_n = m_s_b12_n, m_s_b12_a = m_s_b12_a, m_s_p_i = m_s_p_i, m_s_p_n = m_s_p_n, m_s_p_a = m_s_p_a, m_a_p_i = m_a_p_i, m_a_p_n = m_a_p_n, m_a_p_a = m_a_p_a, m_m_p_i = m_m_p_i, m_m_p_n = m_m_p_n, m_m_p_a = m_m_p_a, m_b_p_i = m_b_p_i, m_b_p_n = m_b_p_n, m_b_p_a = m_b_p_a) # summarise model outputs gls_model_summary <- data.frame(model_id=names(gls_model_results), slope=NA, se=NA, p=NA, pcor=NA, region=NA, var=NA, clade=NA) gls_model_summary$slope <- round(unlist(lapply(gls_model_results, FUN=function(x){return(x$coefficients[2])})),4) gls_model_summary$se <-round(unlist(lapply(gls_model_results, FUN=function(x){return(summary(x)$tTable[2,2])})), 5) gls_model_summary$p <- round(unlist(lapply(gls_model_results, FUN=function(x){return(summary(x)$tTable[2,4])})), 5) gls_model_summary$pcor <- p.adjust(gls_model_summary$p, method ="holm") gls_model_summary$region <- sapply(names(gls_model_results), FUN=function(x){strsplit(x, "_")[[1]][4]}) gls_model_summary$var <- sapply(names(gls_model_results), FUN=function(x){strsplit(x, "_")[[1]][3]}) gls_model_summary$clade <- sapply(names(gls_model_results), FUN=function(x){strsplit(x, "_")[[1]][2]}) gls_model_summary$clade[which(gls_model_summary$clade == "a")] <- "amphibians" gls_model_summary$clade[which(gls_model_summary$clade == "s")] <- "squamates" gls_model_summary$clade[which(gls_model_summary$clade == "m")] <- "mammals" gls_model_summary$clade[which(gls_model_summary$clade == "b")] <- "birds" gls_model_summary$region[which(gls_model_summary$region == "n")] <- "Neotropics" gls_model_summary$region[which(gls_model_summary$region == "i")] <- "Inodmalaya" gls_model_summary$region[which(gls_model_summary$region == "a")] <- "Afrotropics" gls_model_summary$var[which(gls_model_summary$var == "b12")] <- "MAP" gls_model_summary$var[which(gls_model_summary$var == "b1")] <- "MAT" gls_model_summary$var[which(gls_model_summary$var == "p")] <- "PET" write.csv(gls_model_summary, file="Dataset_S3.csv")