#Set working directory setwd("/Users/Macrobird/Google Drive (1)/PhD/Chapter2/MANUSCRIPT/Dryad/") ### INDEX: # 1: Results i) Avian Morphospace. Code to plot Figure 1. # Dataset required: "AllTraits-PCs1-8.csv" # 2: Results: ii) Morphological variance and density maps. Code to plot Figure 2. # Datasets required: "ALLcombinedmetrics.csv", "BirdRealms_Behrman100km.rds" # 3:Results: ii) Morpho Variance and Density vs Species Richness. Code to plot Figure 3. # Datasets required: "ALLcombinedmetrics.csv" # !!!!!!!!!!!!! SPECIES RICHNESS is required and not included in the above dataset as it is derived # from BirdLife Range maps. See below for more info # 4: Results: iii) MMorpho Variance SES vs Density SES maps # Datasets required: "ALLcombinedmetrics.csv", "BirdRealms_Behrman100km.rds" # 5: Results: iv) Environmental and evolutionary drivers. Code for GLS analysis and Figure 5 # Datasets required: "ALLcombinedmetrics.csv" # !!!!!!!!!!!!! SPECIES RICHNESS, ASSEMBLAGE EVOLUTIONARY DISTINCTIVENESS, GROSS PRIMARY PRODUCTIVITIY, # HABITAT HETEROGENEITY, MAIN HABITAT TYPE and ALTITUDINAL RANGE # are required and not included in the above dataset as they are derived # from external data sources. See below for more info ### ------------------------------------- ### ##### 1: Results: i) Avian Morphospace ##### ### ------------------------------------- ### # A, B, C, D and silhouttes were added in powerpoint. Silhouttes were Open Source via phylopic.org # Load libraries. Install if not already. library(ggpointdensity) library(ggplot2) library(viridis) library(ggthemes) # Read files # Contains principal components scores 1 - 8 for each species traitmatrix <- read.csv("AllTraits-PCs1-8.csv", row.names = 1) #### PLOT PC1 and PC2 pdf("PCs1-2.pdf", width = 6, height = 5) ggplot(data = traitmatrix, mapping = aes(x = PC1, y = PC2)) + # Colour points according to the number of neighbouring points within one standard deviation #of the Euclidean distance of each species to all other species. geom_pointdensity(adjust = sd(dist((traitmatrix[,1:2]), method = "euclidean"))) + xlab("PC1 (35.1%)") + ylab("PC2 (10.0%)") + scale_y_continuous(breaks=c(-3,-2,-1,0,1,2,3), labels=c("-3.0","-2.0","-1.0","0.0","1.0","2.0","3.0")) + scale_color_viridis() + theme_classic(base_size = 14) + theme(legend.title = element_blank()) dev.off() #### PLOT PC3 and PC4 pdf("PCs3-4.pdf", width = 6, height = 5) ggplot(data = traitmatrix, mapping = aes(x = PC3, y = PC4)) + geom_pointdensity(adjust = sd(dist((traitmatrix[,3:4]), method = "euclidean"))) + xlab("PC3 (9.3%)") + ylab("PC4 (9.1%)") + scale_y_continuous(breaks=c(-5,-2.5,0,2.5,5.0,7.5), labels=c("-5.0","-2.5","0.0","2.5","5.0","7.5")) + scale_color_viridis() + theme_classic(base_size = 14) + theme(legend.title = element_blank()) dev.off() #### PLOT PC5 and PC6 pdf("PCs5-6.pdf", width = 6, height = 5) ggplot(data = traitmatrix, mapping = aes(x = PC5, y = PC6)) + geom_pointdensity(adjust = sd(dist((traitmatrix[,5:6]), method = "euclidean"))) + xlab("PC5 (9.1%)") + ylab("PC6 (9.1%)") + scale_x_continuous(breaks=c(-8,-4,0,4,8), labels=c("-8.0","-4.0","0.0","4.0","8.0")) + scale_y_continuous(breaks=c(-4,-0,4), labels=c("-4.0","0.0","4.0")) + scale_color_viridis() + theme_classic(base_size = 14) + theme(legend.title = element_blank()) dev.off() #### PLOT PC7 and PC8 pdf("PCs7-8.pdf", width = 6, height = 5) ggplot(data = traitmatrix, mapping = aes(x = PC7, y = PC8)) + geom_pointdensity(adjust = sd(dist((traitmatrix[,7:8]), method = "euclidean"))) + xlab("PC7 (9.1%)") + ylab("PC8 (5.4%)") + scale_y_continuous(breaks=c(-2,0,2), labels=c("-2.0","0.0","2.0")) + scale_x_continuous(breaks=c(-5,0,5,10), labels=c("-5.0","0.0","5.0","10.0")) + scale_color_viridis() + theme_classic(base_size = 14) + theme(legend.title = element_blank()) dev.off() ### ------------------------------------------------------------ ### ##### 2: Results: ii) Morphological variance and density maps ##### ### ------------------------------------------------------------ ### # Annotations were added in powerpoint. library(ggspatial) library(ggplot2) library(viridis) library(ggthemes) library(raster) # Read files # All Morphological Variance and Density values for each grid cell: Morpho_metrics <- read.csv("ALLcombinedmetrics_27OCT20.csv", stringsAsFactors = FALSE) # Raster of realms. Use a raster in the same format your data were calculated from PAM_richness_raster <- readRDS("BirdRealms_Behrman100km.rds") # PLOT FIGURE 2A - MORPHOLOGICAL VARIANCE # Make sure the crs is set to the correct projection crs(PAM_richness_raster) <- "+proj=cea +lat_ts=30 +units=km" # Assign the values of the metric of your choice. e.g. Sum of Variance (Morphological Variance): values(PAM_richness_raster) <- as.numeric(Morpho_metrics$SUM_VAR) pdf("MAP_MORPHO_VARIANCE.pdf", width = 9, height = 7) ggplot() + layer_spatial(PAM_richness_raster) + scale_fill_continuous(name="", type = "viridis", na.value = NA) + theme(legend.position = "bottom", legend.text = element_text(size = 18)) + guides(fill = guide_colourbar(barwidth = 25)) dev.off() # PLOT FIGURE 2B - MORPHOLOGICAL DENSITY values(PAM_richness_raster) <- as.numeric(Morpho_metrics$NEIGH) # Morphological Density pdf("MAP_MORPHO_DENSITY.pdf", width = 9, height = 7) ggplot() + layer_spatial(PAM_richness_raster) + scale_fill_continuous(name="", type = "viridis", na.value = NA) + theme(legend.position = "bottom", legend.text = element_text(size = 18)) + guides(fill = guide_colourbar(barwidth = 25)) dev.off() # PLOT FIGURE 2C - MORPHOLOGICAL VARIANCE SES GLOBAL values(PAM_richness_raster) <- as.numeric(Morpho_metrics$SES_SUM_VAR_G) # Morphological Variance SES (Global) pdf("MORPHO_VARIANCE_SES_GLOBAL.pdf", width = 9, height = 7) ggplot() + layer_spatial(PAM_richness_raster) + scale_fill_continuous(name="", type = "viridis", na.value = NA) + theme(legend.position = "bottom", legend.text = element_text(size = 18)) + guides(fill = guide_colourbar(barwidth = 25)) dev.off() # PLOT FIGURE 2D - MORPHOLOGICAL DENSITY SES GLOBAL values(PAM_richness_raster) <- as.numeric(Morpho_metrics$SES_NEIGH_G) # Morphological Density SES (Global) pdf("MORPHO_DENSITY_SES_GLOBAL.pdf", width = 9, height = 7) ggplot() + layer_spatial(PAM_richness_raster) + scale_fill_continuous(name="", type = "viridis", na.value = NA) + theme(legend.position = "bottom", legend.text = element_text(size = 18)) + guides(fill = guide_colourbar(barwidth = 25)) dev.off() # PLOT FIGURE 2E - MORPHOLOGICAL VARIANCE SES PHYLO values(PAM_richness_raster) <- as.numeric(Morpho_metrics$SES_SUM_VAR_R) # Morphological VARIANCE SES (Phyloregion) pdf("MORPHO_VARIANCE_SES_PHYLO.pdf", width = 9, height = 7) ggplot() + layer_spatial(PAM_richness_raster) + scale_fill_continuous(name="", type = "viridis", na.value = NA) + theme(legend.position = "bottom", legend.text = element_text(size = 18)) + guides(fill = guide_colourbar(barwidth = 25)) dev.off() # PLOT FIGURE 2F - MORPHOLOGICAL DENSITY SES PHYLOREGION values(PAM_richness_raster) <- as.numeric(Morpho_metrics$SES_NEIGH_R) # Morphological Density SES (Phyloregion) pdf("MORPHO_DENSITY_SES_PHYLO.pdf", width = 9, height = 7) ggplot() + layer_spatial(PAM_richness_raster) + scale_fill_continuous(name="", type = "viridis", na.value = NA) + theme(legend.position = "bottom", legend.text = element_text(size = 18)) + guides(fill = guide_colourbar(barwidth = 25)) dev.off() ### -------------------------------------------------------------------- ### ##### 3: Results: ii) Morpho Variance and Density vs Species Richness ##### ### -------------------------------------------------------------------- ### # Load libraries. Install if not already. library(ggpointdensity) library(ggplot2) library(ggpubr) library(ggthemes) # Dataset: # All Morphological Variance and Density values for each grid cell: Morpho_metrics <- read.csv("ALLcombinedmetrics.csv", stringsAsFactors = FALSE) # !!!!!!!!!!!!!!!! # SPECIES RICHNESS is required and not included in the above dataset as it is derived # from BirdLife Range maps. # BirdLife Range maps can be requested from the BirdLife datazone: http://www.birdlife.org/datazone/home # See methods of Hughes et al. Global biogeographic patterns of avian morphological diversity. Ecology Letters. # for detailed info. # In brief, the range maps were used to generate a presence absence matrix containing each species for # each 100km grid cell, under Behrman equal area projection. # We extracted the species richness for each grid cell, and added it to our master dataset: # "ALLcombinedmetrics.csv" ### Species richness (SR) and Morphological Variance (SUM_VAR) VarSR <- ggplot(data = Morpho_metrics, aes(x = SR, y = SUM_VAR)) + geom_pointdensity(size = 0.05) + scale_color_viridis() + # Plot upper quartiles from null expected values (generated from global pool of species): geom_line(data = Morpho_metrics, aes(x = SR , y = SES_SUM_VAR97.5_G), colour = "black") + # Plot lower quartiles from null expected values (generated from global pool of species): geom_line(data = Morpho_metrics, aes(x = SR , y = SES_SUM_VAR2.5_G), colour = "black") + ylab("Morphological Variance") + xlab("Species Richness") + theme_base(base_size=10) + theme(legend.position = "none") ### Species richness (SR) and Morphological Density (NEIGH) DenSR <- ggplot(data = Morpho_metrics, aes(x = SR, y = NEIGH)) + geom_pointdensity(size = 0.05) + scale_color_viridis() + # Plot upper quartiles from null expected values (generated from global pool of species): geom_line(data = Morpho_metrics, aes(x = SR , y = SES_NEIGH97.5_G), colour = "black") + # Plot lower quartiles from null expected values (generated from global pool of species) geom_line(data = Morpho_metrics, aes(x = SR , y = SES_NEIGH2.5_G), colour = "black") + ylab("Morphological Density") + xlab("Species Richness") + theme_base(base_size=10) + theme(legend.position = "none") pdf("SR_Variance_Density.pdf", width = 6, height = 3) ggarrange(VarSR,DenSR, ncol = 2) dev.off() ### -------------------------------------------------------------------- ### ##### 4: Results: iii) Morpho Variance SES vs Density SES maps ##### ### -------------------------------------------------------------------- ### # Load libraries. Install if not already. library(ggplot2) library(viridis) library(ggthemes) library(raster) # Datasets: # All Morphological Variance and Density SES values for each grid cell: Morpho_metrics <- read.csv("MANUSCRIPT/Dryad/ALLcombinedmetrics.csv", stringsAsFactors = FALSE) # Raster of realms. Use a raster in the same format your data were calculated from PAM_richness_raster <- readRDS("MANUSCRIPT/Dryad/BirdRealms_Behrman100km.rds") ## FIG 4 - collated and annoted in powerpoint. The 4 figures that make up Fig 4. are outlined below: ### PLOT global map of variance SES vs density SES calculated from GLOBAL SPECIES POOL crs(PAM_richness_raster) <- "+proj=cea +lat_ts=30 +units=km" values(PAM_richness_raster) <- as.numeric(Morpho_metrics$Fig4_sumvar_neigh_G) PAM_richness_raster@data@values[is.na(PAM_richness_raster@data@values)] <- NA PAM_richness_raster@data@values[(PAM_richness_raster@data@values=='NA')] <- NA # Define colour scheme. Here I use viridis plasma viridis(8, direction = -1, option = "C") cols <- c("#F0F921FF", "#FEBC2AFF", "#F48849FF", "#DB5C68FF", # starting at 1 (> -2 variance SES,> -2 density SES) "#6D6D64", # grey "#B93289FF", "#8B0AA5FF", "#5402A3FF", "#0D0887FF") # ending at 9 (>2 variance SES, > 2 density SES) pdf("Variance_vs_Density_GlobalSES_Map.pdf", width = 18, height = 12) plot(PAM_richness_raster, col = cols, axes = FALSE, box = FALSE, legend = FALSE) dev.off() #### X-Y PLOT of variance SES vs density SES calculated from GLOBAL SPECIES POOL Morpho_metrics <- Morpho_metrics[(Morpho_metrics$Terrestrial_SR9==TRUE),] # Reduce to calculations # Define colour scheme. Here I use viridis plasma table(Morpho_metrics$Fig4_sumvar_neigh_G) # Need to only select colours for those values present # 6 isn't present, so need to remove that from the 9 colour options: viridis(8, direction = -1, option = "C") # colour options 1-4, 6-9 cols2 <- c("#F0F921FF", "#FEBC2AFF", "#F48849FF", "#DB5C68FF", # 1, 2, 3, 4 "#6D6D64", # grey (5) "#8B0AA5FF", "#5402A3FF", "#0D0887FF") #7, 8, 9 pdf("Variance_vs_Density_GlobalSES_CORR.pdf", width = 6, height = 6) ggplot(data = Morpho_metrics, aes(x = SES_SUM_VAR_G, y = SES_NEIGH_G, colour = as.factor(Fig4_sumvar_neigh_G))) + geom_point(size = 0.1) + scale_color_manual(values=cols2) + geom_vline(xintercept=c(-2,2), linetype="dotted") + geom_hline(yintercept=c(-2,2), linetype="dotted") + scale_y_continuous(limits = c(-6.7,2.5), breaks=c(-6,-4,-2,0,2), labels=c("-6","-4","-2","0","2")) + ylab("Morphological Density SES") + xlab("Morphological Variance SES") + theme_foundation() + theme(line = element_line(colour = "black", lineend = "round", linetype = "solid"), rect = element_rect(fill = "white", colour = "black", linetype = "solid"), text = element_text(colour = "black", face = "plain", size = 25, vjust = 0.5, hjust = 0.5, lineheight = 1), panel.grid = element_blank(), strip.background = element_rect(colour = NA), legend.key = element_rect(colour = NA), title = element_text(size = rel(1)), plot.title = element_text(size = rel(1.2), face = "bold"), strip.text = element_text(), legend.title = element_blank(), legend.position = "none", axis.ticks.length = unit(0.5, "lines")) dev.off() ### PLOT global map of variance SES vs density SES calculated from PHYLOREGIONAL SPECIES POOL crs(PAM_richness_raster) <- "+proj=cea +lat_ts=30 +units=km" values(PAM_richness_raster) <- as.numeric(Morpho_metrics$Fig4_sumvar_neigh_R) PAM_richness_raster@data@values[is.na(PAM_richness_raster@data@values)] <- NA PAM_richness_raster@data@values[(PAM_richness_raster@data@values=='NA')] <- NA # Define colour scheme. Here I use viridis plasma viridis(8, direction = -1, option = "C") cols <- c("#F0F921FF", "#FEBC2AFF", "#F48849FF", "#DB5C68FF", # starting at 1 (> -2 variance SES,> -2 density SES) "#6D6D64", # grey "#B93289FF", "#8B0AA5FF", "#5402A3FF", "#0D0887FF") # ending at 9 (>2 variance SES, > 2 density SES) pdf("Variance_vs_Density_PHYLOREGION_SES_Map.pdf", width = 18, height = 12) plot(PAM_richness_raster, col = cols, axes = FALSE, box = FALSE, legend = FALSE) dev.off() #### X-Y PLOT of variance SES vs density SES calculated from GLOBAL SPECIES POOL # Define colour scheme. Here I use viridis plasma table(Morpho_metrics$Fig4_sumvar_neigh_R) # Need to only select colours for those values present # 6 isn't present, so need to remove that from the 9 colour options: viridis(8, direction = -1, option = "C") # colour options 1-4, 6-9 cols2 <- c("#F0F921FF", "#FEBC2AFF", "#F48849FF", "#DB5C68FF", # 1, 2, 3, 4 "#6D6D64", # grey (5) "#8B0AA5FF", "#5402A3FF", "#0D0887FF") #7, 8, 9 pdf("Variance_vs_Density_PHYLOREGION_SES_CORR.pdf", width = 6, height = 6) ggplot(data = Morpho_metrics, aes(x = SES_SUM_VAR_R, y = SES_NEIGH_R, colour = as.factor(Fig4_sumvar_neigh_R))) + geom_point(size = 0.1) + #facet_wrap(vars(RealmName)) + scale_color_manual(values=cols2) + scale_y_continuous(breaks=c(-6,-4,-2,0,2,4), labels=c("-6","-4","-2","0","2","4")) + geom_vline(xintercept=c(-2,2), linetype="dotted") + geom_hline(yintercept=c(-2,2), linetype="dotted") + ylab("Morphological Density SES") + xlab("Morphological Variance SES") + theme_foundation() + theme(line = element_line(colour = "black", lineend = "round", linetype = "solid"), rect = element_rect(fill = "white", colour = "black", linetype = "solid"), text = element_text(colour = "black", face = "plain", size = 25, vjust = 0.5, hjust = 0.5, lineheight = 1), panel.grid = element_blank(), strip.background = element_rect(colour = NA), legend.key = element_rect(colour = NA), title = element_text(size = rel(1)), plot.title = element_text(size = rel(1.2), face = "bold"), strip.text = element_text(), legend.title = element_blank(), legend.position = "none", axis.ticks.length = unit(0.5, "lines")) dev.off() ### --------------------------------------------------------------------------------- ### ##### 5: Results: iv) Environmental and evolutionary drivers.Analysis and Figure 5 ##### ### --------------------------------------------------------------------------------- ### # Load libraries. Install if not already. library(ggplot2) library(viridis) library(ggthemes) library(raster) library(effects) library(ggpubr) library(RColorBrewer) library(sf) library(raster) library(gdalUtilities) library(vegan) library(rgdal) library(nlme) # Datasets: # All Morphological Variance and Density SES values for each grid cell: Morpho_metrics <- read.csv("MANUSCRIPT/Dryad/ALLcombinedmetrics.csv", stringsAsFactors = FALSE) # Figure is annotated in powerpoint # !!!!!!!!!!!!!!!! # SPECIES RICHNESS is required and not included in the above dataset as it is derived # from BirdLife Range maps. # BirdLife Range maps can be requested from the BirdLife datazone: http://www.birdlife.org/datazone/home # See methods of Hughes et al. Global biogeographic patterns of avian morphological diversity. Ecology Letters. # for detailed info. # In brief, the range maps were used to generate a presence absence matrix containing each species for # each 100km grid cell, under Behrman equal area projection. # We extracted the species richness for each grid cell, and added it to our master dataset: # "ALLcombinedmetrics.csv" # !!!!!!!!!!!!!!!! # ASSEMBLAGE EVOLUTIONARY DISTINCTIVENESS SES is required and not included in the above dataset as # it is derived from the BirdTree (Jetz et al. 2012). ED scores for each species were calculated using the # 'equal splits' derivation. SES scores were caclulated. # Trees with the Hackett backbone can be downloaded from http://birdtree.org/ # See methods of Hughes et al. Global biogeographic patterns of avian morphological diversity. Ecology Letters. # for detailed info. # We extracted the Assemblage ED SES for each grid cell, and added it to our master dataset: # "ALLcombinedmetrics.csv" # !!!!!!!!!!!!!!!! # GROSS PRIMARY PRODUCTIVITY (GPP) is required and not included in the above dataset as # it is an external dataset available: https://doi.org/10.6084/m9.figshare.c.3789814 ACCESSED 23rd September 2020 # GPP was calculated for each 100 km grid cell using annual data across the years 2000 – 2016 # For each year, GPP data was transferred to a Behrman projection at 1km resolution, and then the # average GPP value calculated for each 100 km grid cell. The mean value of GPP across years was # then calculated. # See methods of Hughes et al. Global biogeographic patterns of avian morphological diversity. Ecology Letters. # for detailed info. # We extracted the GPP for each grid cell, and added it to our master dataset: # "ALLcombinedmetrics.csv" # !!!!!!!!!!!!!!!! # LAND COVER TYPE and HABITAT HETEROGENEITY is required and not included in the above dataset # as it is an external dataset available: http://doi.org/10.5281/zenodo.3243509 ACCESSED 23rd September 2020 # We used a 2019 dataset of global habitat types from the Copernicus Land Monitoring Services “Global Land Cover Map” # This map was projected under a Behrmann projection, and the modal habitat type from 19 habitat classifications # was calculated for each 1 km grid cell. We then calculated the aggregated modal habitat type (major habitat type), # as well as Shannon’s index (as a score of habitat heterogeneity) for each 100 km grid cell. # See methods of Hughes et al. Global biogeographic patterns of avian morphological diversity. Ecology Letters. # for detailed info. # We extracted the LAND COVER TYPE and HABITAT HETEROGENEITY score for each grid cell, and added it to our master dataset: # "ALLcombinedmetrics.csv" # !!!!!!!!!!!!!!!! # ALTITUDINAL RANGE is required and not included in the above dataset # as it is an external dataset available from the WorldClim database. # Altitudinal data at a 30 arc-second resolution was accessed via the WorldClim database # (Fick & Hijmans 2017) using the getData function in the R package raster (version 3.3.13: Hijmans, 2020). # This data was re-projected to a Behrmann projection, and the average altitude calculated at a 1 km # grid cell resolution. For each 100 km grid cell, the altitudinal range was then calculated from the # minimum and maximum altitudes of the 1 km grid cells contained within. # See methods of Hughes et al. Global biogeographic patterns of avian morphological diversity. Ecology Letters. # for detailed info. # We extracted the ALTITUDINAL RANGE for each grid cell, and added it to our master dataset: # "ALLcombinedmetrics.csv" #### GLS MODELS # Run for each dependent variable: # Morphological Variance SES (sum of variance) GLOBAL # Morphological Variance SES (sum of variance) PHYLOREGIONAL # Morphological Density SES (mean nearest neighbour distance) GLOBAL # Morphological Density SES (mean nearest neighbour distance) GLOBAL ## Reduce to dataset containing no missing values LAND_data <- Morpho_metrics[Morpho_metrics$Terrestrial_SR9==TRUE,] # REMOVE GRID CELLS WITH <9 SPECIES LC_data <- LAND_data[!is.na(LAND_data$LC_SHANNON),] GPP_data <- LC_data[!is.na(LC_data$GPP),] TRIM_data <- GPP_data[!is.na(GPP_data$ALT_RANGE),] ### LOG TRANSFORM NUMERIC VARIABLES TRIM_data$SR_LOG <- log10(TRIM_data$SR) TRIM_data$GPP_LOG <- log10(TRIM_data$GPP+1) TRIM_data$LC_SHANNON_LOG <- log10(TRIM_data$LC_SHANNON+1) TRIM_data$ALT_RANGE_LOG <- log10(TRIM_data$ALT_RANGE+1) ### REDUCE THE SIZE OF THE DATASET - FULL DATASET IS 15461 observations 15277 obs #### RUN BELOW MODELS WITH EACH SUBSET OF DATA #### EXTRACT AIC SCORE FOR EACH MODEL AND SAVE BEST MODEL A <- seq(from = 1, to = 15277, by = 4) DATA_A <- TRIM_data[c(A),] B <- seq(from = 2, to = 15277, by = 4) DATA_B <- TRIM_data[c(B),] C <- seq(from = 3, to = 15277, by = 4) DATA_C <- TRIM_data[c(C),] D <- seq(from = 4, to = 15277, by = 4) DATA_D <- TRIM_data[c(D),] #### SELECT THE DATASET TO RUN MODELS ON DATA_TEST <- DATA_A ### FULL MODEL ########## ### FULL MODEL #### # 1. Fit gls model (nlme library) without accounting for spatial autocorrelation m1 <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE) + LC_DOM, data=DATA_TEST) # 3. If there is evidence for spatial autocorrelation (i.e. the semivariogram is not flat) then fit gls models with different correlation structures (this will be slow) m2 <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE) + LC_DOM, correlation = corExp(form = ~x + y, nugget = TRUE), data = DATA_TEST) m3 <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE) + LC_DOM, correlation = corGaus(form = ~x + y, nugget = TRUE), data = DATA_TEST) m4 <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE) + LC_DOM, correlation = corSpher(form = ~x + y, nugget = TRUE), data = DATA_TEST) # 4. Compare model fits AIC(m1, m2, m3, m4) # Run best model without a nugget - best model is exponential m2_N <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE) + LC_DOM, correlation = corExp(form = ~x + y), data = DATA_TEST) AIC(m1, m2, m3, m4, m2_N) # Best model is exponential without nugget ####### FULL MODEL MINUS LAND-USE TYPE #### Keep nugget, remove Land-use type m1_NoLand <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE), data=DATA_TEST) m2_NoLand <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE), correlation = corExp(form = ~x + y, nugget = TRUE), data = DATA_TEST) m3_NoLand <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE), correlation = corGaus(form = ~x + y, nugget = TRUE), data = DATA_TEST) m4_NoLand <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE), correlation = corSpher(form = ~x + y, nugget = TRUE), data = DATA_TEST) AIC(m1_NoLand, m2_NoLand, m3_NoLand, m4_NoLand) # Best model is exponential #NO NUGGET m2_NoLand_N <- gls(SES_SUM_VAR_G ~ poly(SR_LOG,2, raw = TRUE) + poly(AssemblageED_SES_R,2, raw = TRUE) + poly(GPP_LOG,2, raw = TRUE) + poly(LC_SHANNON_LOG,2, raw = TRUE) + poly(ALT_RANGE_LOG,2, raw = TRUE), correlation = corExp(form = ~x + y), data = DATA_TEST) AIC(m1_NoLand, m2_NoLand, m3_NoLand, m4_NoLand, m2_NoLand_N) # Best model is exponential without nugget ####### LAND-USE TYPE #### Keep nugget, keep land-use only m1_Land <- gls(SES_SUM_VAR_G ~ LC_DOM, data=DATA_TEST) m2_Land <- gls(SES_SUM_VAR_G ~ LC_DOM, correlation = corExp(form = ~x + y, nugget = TRUE), data = DATA_TEST) m3_Land <- gls(SES_SUM_VAR_G ~ LC_DOM, correlation = corGaus(form = ~x + y, nugget = TRUE), data = DATA_TEST) m4_Land <- gls(SES_SUM_VAR_G ~ LC_DOM, correlation = corSpher(form = ~x + y, nugget = TRUE), data = DATA_TEST) #AIC(m1_Land, m2_Land, m3_Land, m4_B_Land) # Best is exponential #NO NUGGET m2_Land_N <- gls(SES_SUM_VAR_R ~ LC_DOM, correlation = corExp(form = ~x + y), data = DATA_TEST) #AIC(m1_Land, m2_Land, m3_Land, m4_Land, m2_Land_N) # Best is exponential without nugget ######## END - VARIOGRAMS ### PUT IN BEST MODELS #### MOST SUPPORT FOR MODEL WITHOUT LAND_TYPE (LC_DOM), EXPONENTIAL CORR STRUCTURE AND NO NUGGET # 2. Fit semivariogram and plot to assess deree of spatial autocorrelation (looking to see if semivariogram (y axis) increases with distance) vario1 <- Variogram(m2_NoLand_N, form = ~x + y, resType = "pearson", data=DATA_TEST) plot(vario1, smooth = TRUE) # Check variogram sit for the best model - estimated line should fit data well vario2 <- Variogram(m2_NoLand_N, form = ~x + y, resType = "pearson") plot(vario2, smooth = FALSE, ylim = c(0, 1.2)) # Check variogram sit for the best model - data should have no trend vario3 <- Variogram(m2_NoLand_N, form = ~x + y, resType = "normalized") plot(vario3, smooth = FALSE, ylim = c(0, 1.2)) #### OUTPUT STATS summary(m2_NoLand_N) anova(m2_NoLand_N) ##### SAVE BEST MODELS ### ALL MINUS LAND USE (BEST MODEL) m2_B_NoLand_N no nugget saveRDS(m2_NoLand_N, "A_EcoVarG_MinusLandUse_Exp_no-nug.RDS") #### PLOTTING FIGURE 5 DATA_TEST <- TRIM_data #### ECOMORPHOLOGICAL VARIANCE SES (GLOBAL) ##### ### DATASET A ##### BestA <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/A_EcoVolG_MinusLandUse_Exp_no-nug.RDS") summary(BestA) # significant = sum_var_log, neigh_log, gpp_log, shannon_log, alt_range_log # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestA) SR_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_A) <- "x" SR_LOG_A$y <- predictX2$fit[,1] SR_LOG_A$upper <- predictX2$upper[,1] SR_LOG_A$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestA) ED_SES_A <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_A) <- "x" ED_SES_A$y <- predictX2$fit[,1] ED_SES_A$upper <- predictX2$upper[,1] ED_SES_A$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestA) GPP_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_A) <- "x" GPP_LOG_A$y <- predictX2$fit[,1] GPP_LOG_A$upper <- predictX2$upper[,1] GPP_LOG_A$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestA) SHANNON_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_A) <- "x" SHANNON_LOG_A$y <- predictX2$fit[,1] SHANNON_LOG_A$upper <- predictX2$upper[,1] SHANNON_LOG_A$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestA) ALT_RANGE_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_A) <- "x" ALT_RANGE_LOG_A$y <- predictX2$fit[,1] ALT_RANGE_LOG_A$upper <- predictX2$upper[,1] ALT_RANGE_LOG_A$lower <- predictX2$lower[,1] ##### ### DATASET B ###### BestB <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/B_EcoVolG_MinusLandUse_Exp_nug.RDS") summary(BestB) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestB) SR_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_B) <- "x" SR_LOG_B$y <- predictX2$fit[,1] SR_LOG_B$upper <- predictX2$upper[,1] SR_LOG_B$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestB) ED_SES_B <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_B) <- "x" ED_SES_B$y <- predictX2$fit[,1] ED_SES_B$upper <- predictX2$upper[,1] ED_SES_B$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestB) GPP_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_B) <- "x" GPP_LOG_B$y <- predictX2$fit[,1] GPP_LOG_B$upper <- predictX2$upper[,1] GPP_LOG_B$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestB) SHANNON_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_B) <- "x" SHANNON_LOG_B$y <- predictX2$fit[,1] SHANNON_LOG_B$upper <- predictX2$upper[,1] SHANNON_LOG_B$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestB) ALT_RANGE_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_B) <- "x" ALT_RANGE_LOG_B$y <- predictX2$fit[,1] ALT_RANGE_LOG_B$upper <- predictX2$upper[,1] ALT_RANGE_LOG_B$lower <- predictX2$lower[,1] ##### ### DATASET C ###### BestC <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/C_EcoVolG_MinusLandUse_Exp_no-nug.RDS") summary(BestC) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestC) SR_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_C) <- "x" SR_LOG_C$y <- predictX2$fit[,1] SR_LOG_C$upper <- predictX2$upper[,1] SR_LOG_C$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestC) ED_SES_C <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_C) <- "x" ED_SES_C$y <- predictX2$fit[,1] ED_SES_C$upper <- predictX2$upper[,1] ED_SES_C$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestC) GPP_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_C) <- "x" GPP_LOG_C$y <- predictX2$fit[,1] GPP_LOG_C$upper <- predictX2$upper[,1] GPP_LOG_C$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestB) SHANNON_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_C) <- "x" SHANNON_LOG_C$y <- predictX2$fit[,1] SHANNON_LOG_C$upper <- predictX2$upper[,1] SHANNON_LOG_C$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestC) ALT_RANGE_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_C) <- "x" ALT_RANGE_LOG_C$y <- predictX2$fit[,1] ALT_RANGE_LOG_C$upper <- predictX2$upper[,1] ALT_RANGE_LOG_C$lower <- predictX2$lower[,1] ##### ### DATASET D ###### BestD <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/D_EcoVolG_MinusLandUse_Exp_nug.RDS") summary(BestD) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestD) SR_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_D) <- "x" SR_LOG_D$y <- predictX2$fit[,1] SR_LOG_D$upper <- predictX2$upper[,1] SR_LOG_D$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestD) ED_SES_D <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_D) <- "x" ED_SES_D$y <- predictX2$fit[,1] ED_SES_D$upper <- predictX2$upper[,1] ED_SES_D$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestD) GPP_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_D) <- "x" GPP_LOG_D$y <- predictX2$fit[,1] GPP_LOG_D$upper <- predictX2$upper[,1] GPP_LOG_D$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestD) SHANNON_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_D) <- "x" SHANNON_LOG_D$y <- predictX2$fit[,1] SHANNON_LOG_D$upper <- predictX2$upper[,1] SHANNON_LOG_D$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestD) ALT_RANGE_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_D) <- "x" ALT_RANGE_LOG_D$y <- predictX2$fit[,1] ALT_RANGE_LOG_D$upper <- predictX2$upper[,1] ALT_RANGE_LOG_D$lower <- predictX2$lower[,1] ##### ##### SR PLOTS ####### SumG_SR <- ggplot(data = TRIM_data, aes(x = SR_LOG, y = SES_SUM_VAR_G)) + #xlab("log10(Species Richness)") + ylab("Variance SES") + scale_y_continuous(limits = c(-1,7), breaks=c(-1,3,7), labels=c("-1","3","7")) + scale_x_continuous(breaks=c(1,2,3), labels=c("1","2","3")) + geom_ribbon(data = SR_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SR_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SR_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SR_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SR_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = SR_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = SR_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = SR_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ED PLOTS SumG_EDG <- ggplot(data = TRIM_data, aes(x = AssemblageED_SES_G, y = SES_SUM_VAR_G)) + #xlab("Assemblage ED SES") + ylab("Variance SES") + scale_y_continuous(limits = c(-1,7), breaks=c(-1,3,7), labels=c("-1","3","7")) + scale_x_continuous(breaks=c(0,6), labels=c("0","6")) + geom_ribbon(data = ED_SES_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ED_SES_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ED_SES_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ED_SES_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ED_SES_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ED_SES_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ED_SES_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ED_SES_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### GPP PLOTS SumG_GPP <- ggplot(data = TRIM_data, aes(x = GPP_LOG, y = SES_SUM_VAR_G)) + #xlab("log10(GPP)") + ylab("Variance SES") + scale_y_continuous(limits = c(-1,7), breaks=c(-1,3,7), labels=c("-1","3","7")) + geom_ribbon(data = GPP_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = GPP_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = GPP_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = GPP_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = GPP_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1], linetype = 2) + geom_line(data = GPP_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2], linetype = 2) + geom_line(data = GPP_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3], linetype = 2) + geom_line(data = GPP_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4], linetype = 2) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ALTITUDINAL RANGE SumG_ALT <- ggplot(data = TRIM_data, aes(x = ALT_RANGE_LOG, y = SES_SUM_VAR_G)) + #xlab("log10(Altitudinal Range)") + ylab("Variance SES") + scale_y_continuous(limits = c(-1,7), breaks=c(-1,3,7), labels=c("-1","3","7")) + scale_x_continuous(breaks=c(0,2,4), labels=c("0","2","4")) + geom_ribbon(data = ALT_RANGE_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ALT_RANGE_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ALT_RANGE_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ALT_RANGE_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ALT_RANGE_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### HABITAT HETEROGENEITY SumG_SHAN <- ggplot(data = TRIM_data, aes(x = LC_SHANNON_LOG, y = SES_SUM_VAR_G)) + #xlab("log10(Habitat Heterogeneity)") + ylab("Variance SES") + scale_y_continuous(limits = c(-1,7), breaks=c(-1,3,7), labels=c("-1","3","7")) + scale_x_continuous(breaks=c(0,0.2,0.4), labels=c("0.0","0.2","0.4")) + geom_ribbon(data = SHANNON_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SHANNON_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1], linetype = 2) + geom_line(data = SHANNON_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2], linetype = 2) + geom_line(data = SHANNON_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3], linetype = 2) + geom_line(data = SHANNON_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") #### ECOMORPHOLOGICAL VARIANCE SES (PHYLOREGION) ##### ### DATASET A ##### BestA <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/A_EcoVolR_MinusLandUse_Exp_no-nug.RDS") summary(BestA) # significant = sum_var_log, neigh_log, gpp_log, shannon_log, alt_range_log # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestA) SR_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_A) <- "x" SR_LOG_A$y <- predictX2$fit[,1] SR_LOG_A$upper <- predictX2$upper[,1] SR_LOG_A$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestA) ED_SES_A <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_A) <- "x" ED_SES_A$y <- predictX2$fit[,1] ED_SES_A$upper <- predictX2$upper[,1] ED_SES_A$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestA) GPP_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_A) <- "x" GPP_LOG_A$y <- predictX2$fit[,1] GPP_LOG_A$upper <- predictX2$upper[,1] GPP_LOG_A$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestA) SHANNON_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_A) <- "x" SHANNON_LOG_A$y <- predictX2$fit[,1] SHANNON_LOG_A$upper <- predictX2$upper[,1] SHANNON_LOG_A$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestA) ALT_RANGE_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_A) <- "x" ALT_RANGE_LOG_A$y <- predictX2$fit[,1] ALT_RANGE_LOG_A$upper <- predictX2$upper[,1] ALT_RANGE_LOG_A$lower <- predictX2$lower[,1] ##### ### DATASET B ###### BestB <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/B_EcoVolR_MinusLandUse_Exp_nug.RDS") summary(BestB) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestB) SR_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_B) <- "x" SR_LOG_B$y <- predictX2$fit[,1] SR_LOG_B$upper <- predictX2$upper[,1] SR_LOG_B$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestB) ED_SES_B <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_B) <- "x" ED_SES_B$y <- predictX2$fit[,1] ED_SES_B$upper <- predictX2$upper[,1] ED_SES_B$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestB) GPP_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_B) <- "x" GPP_LOG_B$y <- predictX2$fit[,1] GPP_LOG_B$upper <- predictX2$upper[,1] GPP_LOG_B$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestB) SHANNON_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_B) <- "x" SHANNON_LOG_B$y <- predictX2$fit[,1] SHANNON_LOG_B$upper <- predictX2$upper[,1] SHANNON_LOG_B$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestB) ALT_RANGE_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_B) <- "x" ALT_RANGE_LOG_B$y <- predictX2$fit[,1] ALT_RANGE_LOG_B$upper <- predictX2$upper[,1] ALT_RANGE_LOG_B$lower <- predictX2$lower[,1] ##### ### DATASET C ##### BestC <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/C_EcoVolR_MinusLandUse_Exp_no-nug.RDS") summary(BestC) # significant = sum_var_log, neigh_log, gpp_log, shannon_log, alt_range_log # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestC) SR_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_C) <- "x" SR_LOG_C$y <- predictX2$fit[,1] SR_LOG_C$upper <- predictX2$upper[,1] SR_LOG_C$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestC) ED_SES_C <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_C) <- "x" ED_SES_C$y <- predictX2$fit[,1] ED_SES_C$upper <- predictX2$upper[,1] ED_SES_C$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestC) GPP_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_C) <- "x" GPP_LOG_C$y <- predictX2$fit[,1] GPP_LOG_C$upper <- predictX2$upper[,1] GPP_LOG_C$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestC) SHANNON_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_C) <- "x" SHANNON_LOG_C$y <- predictX2$fit[,1] SHANNON_LOG_C$upper <- predictX2$upper[,1] SHANNON_LOG_C$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestC) ALT_RANGE_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_C) <- "x" ALT_RANGE_LOG_C$y <- predictX2$fit[,1] ALT_RANGE_LOG_C$upper <- predictX2$upper[,1] ALT_RANGE_LOG_C$lower <- predictX2$lower[,1] ##### ### DATASET D ##### BestD <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoVol/D_EcoVolR_MinusLandUse_Exp_nug.RDS") summary(BestD) # significant = sum_var_log, neigh_log, gpp_log, shannon_log, alt_range_log # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestD) SR_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_D) <- "x" SR_LOG_D$y <- predictX2$fit[,1] SR_LOG_D$upper <- predictX2$upper[,1] SR_LOG_D$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestD) ED_SES_D <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_D) <- "x" ED_SES_D$y <- predictX2$fit[,1] ED_SES_D$upper <- predictX2$upper[,1] ED_SES_D$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestD) GPP_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_D) <- "x" GPP_LOG_D$y <- predictX2$fit[,1] GPP_LOG_D$upper <- predictX2$upper[,1] GPP_LOG_D$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestD) SHANNON_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_D) <- "x" SHANNON_LOG_D$y <- predictX2$fit[,1] SHANNON_LOG_D$upper <- predictX2$upper[,1] SHANNON_LOG_D$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestD) ALT_RANGE_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_D) <- "x" ALT_RANGE_LOG_D$y <- predictX2$fit[,1] ALT_RANGE_LOG_D$upper <- predictX2$upper[,1] ALT_RANGE_LOG_D$lower <- predictX2$lower[,1] ##### ##### SR PLOTS SumR_SR <- ggplot(data = TRIM_data, aes(x = SR_LOG, y = SES_SUM_VAR_R)) + scale_y_continuous(limits = c(-2,7.5), breaks=c(-1,3,7), labels=c("-1","3","7")) + scale_x_continuous(breaks=c(1,2,3), labels=c("1","2","3")) + #xlab("log10(Species Richness)") + ylab("Variance SES") + geom_ribbon(data = SR_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SR_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SR_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SR_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SR_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = SR_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = SR_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = SR_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ED PLOTS SumR_EDR <- ggplot(data = TRIM_data, aes(x = AssemblageED_SES_R, y = SES_SUM_VAR_R)) + #xlab("Assemblage ED SES") + ylab("Variance SES") + scale_y_continuous(limits = c(-2,7.5), breaks=c(-1,3,7), labels=c("-1","3","7")) + geom_ribbon(data = ED_SES_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ED_SES_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ED_SES_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ED_SES_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ED_SES_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ED_SES_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ED_SES_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ED_SES_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### GPP PLOTS SumR_GPP <- ggplot(data = TRIM_data, aes(x = GPP_LOG, y = SES_SUM_VAR_R)) + #xlab("log10(GPP)") + ylab("Variance SES") + scale_y_continuous(limits = c(-2,7.5), breaks=c(-1,3,7), labels=c("-1","3","7")) + geom_ribbon(data = GPP_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = GPP_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = GPP_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = GPP_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = GPP_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1], linetype = 2) + geom_line(data = GPP_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2], linetype = 2) + geom_line(data = GPP_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3], linetype = 2) + geom_line(data = GPP_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4], linetype = 2) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ALTITUDINAL RANGE SumR_ALT <- ggplot(data = TRIM_data, aes(x = ALT_RANGE_LOG, y = SES_SUM_VAR_R)) + #xlab("log10(Altitudinal Range)") + ylab("Variance SES") + scale_y_continuous(limits = c(-2,7.5), breaks=c(-1,3,7), labels=c("-1","3","7")) + scale_x_continuous(breaks=c(0,2,4), labels=c("0","2","4")) + geom_ribbon(data = ALT_RANGE_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ALT_RANGE_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ALT_RANGE_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ALT_RANGE_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ALT_RANGE_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### HABITAT HETEROGENEITY SumR_SHAN <- ggplot(data = TRIM_data, aes(x = LC_SHANNON_LOG, y = SES_SUM_VAR_R)) + #xlab("log10(Habitat Heterogeneity)") + ylab("Variance SES") + scale_y_continuous(limits = c(-2,7.5), breaks=c(-1,3,7), labels=c("-1","3","7")) + scale_x_continuous(breaks = c(0,0.2,0.4), labels=c("0.0","0.2","0.4")) + geom_ribbon(data = SHANNON_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SHANNON_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1], linetype = 2) + geom_line(data = SHANNON_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2], linetype = 2) + geom_line(data = SHANNON_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3], linetype = 2) + geom_line(data = SHANNON_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4], linetype = 2) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") #### ECOMORPHOLOGICAL DENSITY SES (GLOBAL) ##### ### DATASET A ##### BestA <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/A_EcoDisG_MinusLandUse_Exp_nug.RDS") summary(BestA) # significant = sum_var_log, neigh_log, gpp_log, shannon_log, alt_range_log # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestA) SR_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_A) <- "x" SR_LOG_A$y <- predictX2$fit[,1] SR_LOG_A$upper <- predictX2$upper[,1] SR_LOG_A$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestA) ED_SES_A <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_A) <- "x" ED_SES_A$y <- predictX2$fit[,1] ED_SES_A$upper <- predictX2$upper[,1] ED_SES_A$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestA) GPP_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_A) <- "x" GPP_LOG_A$y <- predictX2$fit[,1] GPP_LOG_A$upper <- predictX2$upper[,1] GPP_LOG_A$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestA) SHANNON_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_A) <- "x" SHANNON_LOG_A$y <- predictX2$fit[,1] SHANNON_LOG_A$upper <- predictX2$upper[,1] SHANNON_LOG_A$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestA) ALT_RANGE_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_A) <- "x" ALT_RANGE_LOG_A$y <- predictX2$fit[,1] ALT_RANGE_LOG_A$upper <- predictX2$upper[,1] ALT_RANGE_LOG_A$lower <- predictX2$lower[,1] ##### ### DATASET B ###### BestB <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/B_EcoDisG_MinusLandUse_Exp_nug.RDS") summary(BestB) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestB) SR_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_B) <- "x" SR_LOG_B$y <- predictX2$fit[,1] SR_LOG_B$upper <- predictX2$upper[,1] SR_LOG_B$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestB) ED_SES_B <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_B) <- "x" ED_SES_B$y <- predictX2$fit[,1] ED_SES_B$upper <- predictX2$upper[,1] ED_SES_B$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestB) GPP_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_B) <- "x" GPP_LOG_B$y <- predictX2$fit[,1] GPP_LOG_B$upper <- predictX2$upper[,1] GPP_LOG_B$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestB) SHANNON_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_B) <- "x" SHANNON_LOG_B$y <- predictX2$fit[,1] SHANNON_LOG_B$upper <- predictX2$upper[,1] SHANNON_LOG_B$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestB) ALT_RANGE_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_B) <- "x" ALT_RANGE_LOG_B$y <- predictX2$fit[,1] ALT_RANGE_LOG_B$upper <- predictX2$upper[,1] ALT_RANGE_LOG_B$lower <- predictX2$lower[,1] ##### ### DATASET C ##### BestC <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/C_EcoDisG_MinusLandUse_Exp_nug.RDS") summary(BestC) # significant = sum_var_log, neigh_log, gpp_log, shannon_log, alt_range_log # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestC) SR_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_C) <- "x" SR_LOG_C$y <- predictX2$fit[,1] SR_LOG_C$upper <- predictX2$upper[,1] SR_LOG_C$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestC) ED_SES_C <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_C) <- "x" ED_SES_C$y <- predictX2$fit[,1] ED_SES_C$upper <- predictX2$upper[,1] ED_SES_C$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestC) GPP_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_C) <- "x" GPP_LOG_C$y <- predictX2$fit[,1] GPP_LOG_C$upper <- predictX2$upper[,1] GPP_LOG_C$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestC) SHANNON_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_C) <- "x" SHANNON_LOG_C$y <- predictX2$fit[,1] SHANNON_LOG_C$upper <- predictX2$upper[,1] SHANNON_LOG_C$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestC) ALT_RANGE_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_C) <- "x" ALT_RANGE_LOG_C$y <- predictX2$fit[,1] ALT_RANGE_LOG_C$upper <- predictX2$upper[,1] ALT_RANGE_LOG_C$lower <- predictX2$lower[,1] ##### ### DATASET D ##### BestD <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/D_EcoDisG_MinusLandUse_Exp_nug.RDS") summary(BestD) # significant = sum_var_log, neigh_log, gpp_log, shannon_log, alt_range_log # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestD) SR_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_D) <- "x" SR_LOG_D$y <- predictX2$fit[,1] SR_LOG_D$upper <- predictX2$upper[,1] SR_LOG_D$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_G", mod=BestD) ED_SES_D <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_D) <- "x" ED_SES_D$y <- predictX2$fit[,1] ED_SES_D$upper <- predictX2$upper[,1] ED_SES_D$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestD) GPP_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_D) <- "x" GPP_LOG_D$y <- predictX2$fit[,1] GPP_LOG_D$upper <- predictX2$upper[,1] GPP_LOG_D$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestD) SHANNON_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_D) <- "x" SHANNON_LOG_D$y <- predictX2$fit[,1] SHANNON_LOG_D$upper <- predictX2$upper[,1] SHANNON_LOG_D$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestD) ALT_RANGE_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_D) <- "x" ALT_RANGE_LOG_D$y <- predictX2$fit[,1] ALT_RANGE_LOG_D$upper <- predictX2$upper[,1] ALT_RANGE_LOG_D$lower <- predictX2$lower[,1] ##### ##### SR PLOTS NeighG_SR <- ggplot(data = TRIM_data, aes(x = SR_LOG, y = SES_NEIGH_G)) + #xlab("log10 (Species Richness)") + ylab("Density SES (Global)") + scale_y_continuous(limits = c(-5.5,1), breaks=c(-4,-2,0), labels=c("-4","-2","0")) + scale_x_continuous(breaks=c(1,2,3), labels=c("1","2","3")) + geom_ribbon(data = SR_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SR_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SR_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SR_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SR_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = SR_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = SR_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = SR_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ED PLOTS NeighG_EDG <- ggplot(data = TRIM_data, aes(x = AssemblageED_SES_G, y = SES_NEIGH_G)) + #xlab("Assemblage ED SES (Global)") + ylab("Density SES (Global)") + scale_y_continuous(limits = c(-5.5,1), breaks=c(-4,-2,0), labels=c("-4","-2","0")) + scale_x_continuous(breaks=c(0,6), labels=c("0","6")) + geom_ribbon(data = ED_SES_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ED_SES_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ED_SES_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ED_SES_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ED_SES_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ED_SES_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ED_SES_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ED_SES_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### GPP PLOTS NeighG_GPP <- ggplot(data = TRIM_data, aes(x = GPP_LOG, y = SES_NEIGH_G)) + #geom_point(size = 1) + #xlab("log10 (GPP)") + ylab("Density SES (Global)") + scale_y_continuous(limits = c(-5.5,1), breaks=c(-4,-2,0), labels=c("-4","-2","0")) + scale_x_continuous(breaks=c(0,1,2,3), labels=c("0","1","2","3")) + geom_ribbon(data = GPP_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = GPP_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = GPP_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = GPP_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = GPP_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = GPP_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = GPP_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = GPP_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ALTITUDINAL RANGE NeighG_ALT <- ggplot(data = TRIM_data, aes(x = ALT_RANGE_LOG, y = SES_NEIGH_G)) + #xlab("log10 (Altitudinal Range)") + ylab("Density SES (Global)") + scale_y_continuous(limits = c(-5.5,1), breaks=c(-4,-2,0), labels=c("-4","-2","0")) + scale_x_continuous(breaks=c(0,2,4), labels=c("0","2","4")) + geom_ribbon(data = ALT_RANGE_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ALT_RANGE_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ALT_RANGE_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ALT_RANGE_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ALT_RANGE_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### HABITAT HETEROGENEITY NeighG_SHAN <- ggplot(data = TRIM_data, aes(x = LC_SHANNON_LOG, y = SES_NEIGH_G)) + #xlab("log10 (Habitat Heterogeneity)") + ylab("Density SES (Global)") + scale_y_continuous(limits = c(-5.5,1), breaks=c(-4,-2,0), labels=c("-4","-2","0")) + scale_x_continuous(breaks=c(0.0,0.2,0.4), labels=c("0.0","0.2","0.4")) + geom_ribbon(data = SHANNON_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SHANNON_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1], linetype = 2) + geom_line(data = SHANNON_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2], linetype = 2) + geom_line(data = SHANNON_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3], linetype = 2) + geom_line(data = SHANNON_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4], linetype = 2) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") #### ECOMORPHOLOGICAL DENSITY SES (PHYLOREGION) ##### ### DATASET A ##### BestA <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/A_EcoDisR_MinusLandUse_Exp_nug.RDS") summary(BestA) # significant = SR,ED,GPP,ALT # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestA) SR_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_A) <- "x" SR_LOG_A$y <- predictX2$fit[,1] SR_LOG_A$upper <- predictX2$upper[,1] SR_LOG_A$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestA) ED_SES_A <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_A) <- "x" ED_SES_A$y <- predictX2$fit[,1] ED_SES_A$upper <- predictX2$upper[,1] ED_SES_A$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestA) GPP_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_A) <- "x" GPP_LOG_A$y <- predictX2$fit[,1] GPP_LOG_A$upper <- predictX2$upper[,1] GPP_LOG_A$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestA) SHANNON_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_A) <- "x" SHANNON_LOG_A$y <- predictX2$fit[,1] SHANNON_LOG_A$upper <- predictX2$upper[,1] SHANNON_LOG_A$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestA) ALT_RANGE_LOG_A <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_A) <- "x" ALT_RANGE_LOG_A$y <- predictX2$fit[,1] ALT_RANGE_LOG_A$upper <- predictX2$upper[,1] ALT_RANGE_LOG_A$lower <- predictX2$lower[,1] ##### ### DATASET B ###### BestB <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/B_EcoDisR_MinusLandUse_Exp_nug.RDS") summary(BestB) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestB) SR_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_B) <- "x" SR_LOG_B$y <- predictX2$fit[,1] SR_LOG_B$upper <- predictX2$upper[,1] SR_LOG_B$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestB) ED_SES_B <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_B) <- "x" ED_SES_B$y <- predictX2$fit[,1] ED_SES_B$upper <- predictX2$upper[,1] ED_SES_B$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestB) GPP_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_B) <- "x" GPP_LOG_B$y <- predictX2$fit[,1] GPP_LOG_B$upper <- predictX2$upper[,1] GPP_LOG_B$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestB) SHANNON_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_B) <- "x" SHANNON_LOG_B$y <- predictX2$fit[,1] SHANNON_LOG_B$upper <- predictX2$upper[,1] SHANNON_LOG_B$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestB) ALT_RANGE_LOG_B <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_B) <- "x" ALT_RANGE_LOG_B$y <- predictX2$fit[,1] ALT_RANGE_LOG_B$upper <- predictX2$upper[,1] ALT_RANGE_LOG_B$lower <- predictX2$lower[,1] ##### ### DATASET C ###### BestC <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/C_EcoDisR_MinusLandUse_Exp_nug.RDS") summary(BestC) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestC) SR_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_C) <- "x" SR_LOG_C$y <- predictX2$fit[,1] SR_LOG_C$upper <- predictX2$upper[,1] SR_LOG_C$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestC) ED_SES_C <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_C) <- "x" ED_SES_C$y <- predictX2$fit[,1] ED_SES_C$upper <- predictX2$upper[,1] ED_SES_C$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestC) GPP_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_C) <- "x" GPP_LOG_C$y <- predictX2$fit[,1] GPP_LOG_C$upper <- predictX2$upper[,1] GPP_LOG_C$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestC) SHANNON_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_C) <- "x" SHANNON_LOG_C$y <- predictX2$fit[,1] SHANNON_LOG_C$upper <- predictX2$upper[,1] SHANNON_LOG_C$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestC) ALT_RANGE_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_C) <- "x" ALT_RANGE_LOG_C$y <- predictX2$fit[,1] ALT_RANGE_LOG_C$upper <- predictX2$upper[,1] ALT_RANGE_LOG_C$lower <- predictX2$lower[,1] ##### ### DATASET D ###### BestD <- readRDS("DATA/OUTPUTS/GLS OUTPUTS/EcoDis/D_EcoDisR_MinusLandUse_Exp_nug.RDS") summary(BestD) # significant predictors: sum var, neigh, gpp, alt_range # SR predictX2 <- predictorEffect(predictor="SR_LOG", mod=BestD) SR_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(SR_LOG_D) <- "x" SR_LOG_D$y <- predictX2$fit[,1] SR_LOG_D$upper <- predictX2$upper[,1] SR_LOG_D$lower <- predictX2$lower[,1] # ED predictX2 <- predictorEffect(predictor="AssemblageED_SES_R", mod=BestD) ED_SES_D <- as.data.frame(predictX2$x[,1]) colnames(ED_SES_D) <- "x" ED_SES_D$y <- predictX2$fit[,1] ED_SES_D$upper <- predictX2$upper[,1] ED_SES_D$lower <- predictX2$lower[,1] # GPP predictX2 <- predictorEffect(predictor="GPP_LOG", mod=BestD) GPP_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(GPP_LOG_D) <- "x" GPP_LOG_D$y <- predictX2$fit[,1] GPP_LOG_D$upper <- predictX2$upper[,1] GPP_LOG_D$lower <- predictX2$lower[,1] # SHANNON predictX2 <- predictorEffect(predictor="LC_SHANNON_LOG", mod=BestD) SHANNON_LOG_C <- as.data.frame(predictX2$x[,1]) colnames(SHANNON_LOG_C) <- "x" SHANNON_LOG_C$y <- predictX2$fit[,1] SHANNON_LOG_C$upper <- predictX2$upper[,1] SHANNON_LOG_C$lower <- predictX2$lower[,1] # ALTITUDE predictX2 <- predictorEffect(predictor="ALT_RANGE_LOG", mod=BestD) ALT_RANGE_LOG_D <- as.data.frame(predictX2$x[,1]) colnames(ALT_RANGE_LOG_D) <- "x" ALT_RANGE_LOG_D$y <- predictX2$fit[,1] ALT_RANGE_LOG_D$upper <- predictX2$upper[,1] ALT_RANGE_LOG_D$lower <- predictX2$lower[,1] ##### ##### SR PLOTS NeighR_SR <- ggplot(data = TRIM_data, aes(x = SR_LOG, y = SES_NEIGH_R)) + #xlab("log10 (Species Richness)") + ylab("Density SES (Phyloregion)") + scale_y_continuous(limits = c(-2.5,2.5), breaks=c(-2,0,2), labels=c("-2","0","2")) + scale_x_continuous(breaks=c(1,2,3), labels=c("1","2","3")) + geom_ribbon(data = SR_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SR_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SR_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SR_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SR_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = SR_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = SR_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3], linetype = 2) + geom_line(data = SR_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ED PLOTS NeighR_EDR <- ggplot(data = TRIM_data, aes(x = AssemblageED_SES_R, y = SES_NEIGH_R)) + #xlab("Assemblage ED SES (Phyloregion)") + ylab("Density SES (Phyloregion)") + scale_y_continuous(limits = c(-2.5,2.5), breaks=c(-2,0,2), labels=c("-2","0","2")) + geom_ribbon(data = ED_SES_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ED_SES_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ED_SES_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ED_SES_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ED_SES_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ED_SES_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ED_SES_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ED_SES_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### GPP PLOTS NeighR_GPP <- ggplot(data = TRIM_data, aes(x = GPP_LOG, y = SES_NEIGH_R)) + #geom_point(size = 1) + #xlab("log10 (GPP)") + ylab("Density SES (Phyloregion)") + scale_y_continuous(limits = c(-2.5,2.5), breaks=c(-2,0,2), labels=c("-2","0","2")) + geom_ribbon(data = GPP_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = GPP_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = GPP_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = GPP_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = GPP_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = GPP_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = GPP_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = GPP_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### ALTITUDINAL RANGE NeighR_ALT <- ggplot(data = TRIM_data, aes(x = ALT_RANGE_LOG, y = SES_NEIGH_R)) + #xlab("log10 (Altitudinal Range)") + ylab("Density SES (Phyloregion)") + scale_y_continuous(limits = c(-2.5,2.5), breaks=c(-2,0,2), labels=c("-2","0","2")) + scale_x_continuous(breaks=c(0,2,4), labels=c("0","2","4")) + geom_ribbon(data = ALT_RANGE_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = ALT_RANGE_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = ALT_RANGE_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1]) + geom_line(data = ALT_RANGE_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2]) + geom_line(data = ALT_RANGE_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3]) + geom_line(data = ALT_RANGE_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4]) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") ##### HABITAT HETEROGENEITY NeighR_SHAN <- ggplot(data = TRIM_data, aes(x = LC_SHANNON_LOG, y = SES_NEIGH_R)) + xlab("log10 (Habitat Heterogeneity)") + ylab("Density SES (Phyloregion)") + scale_y_continuous(limits = c(-2.5,2.5), breaks=c(-2,0,2), labels=c("-2","0","2")) + scale_x_continuous(breaks=c(0.0,0.2,0.4), labels=c("0.0","0.2","0.4")) + geom_ribbon(data = SHANNON_LOG_A, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[1], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_B, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[2], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_C, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[3], alpha = 0.05) + geom_ribbon(data = SHANNON_LOG_D, aes(x = x, y = y, ymin = lower, ymax = upper), fill = brewer.pal(4, name = "Dark2")[4], alpha = 0.05) + geom_line(data = SHANNON_LOG_A, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[1], linetype = 2) + geom_line(data = SHANNON_LOG_B, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[2], linetype = 2) + geom_line(data = SHANNON_LOG_C, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[3], linetype = 2) + geom_line(data = SHANNON_LOG_D, aes(x = x, y = y), colour = brewer.pal(4, name = "Dark2")[4], linetype = 2) + theme_base(base_size=14) + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), legend.position = "none") pdf("SESVAR_SESDIS_ENVTRAITS.pdf", width = 12, height = 10) ggarrange(SumG_SR, SumG_EDG, SumG_GPP, SumG_SHAN, SumG_ALT, SumR_SR, SumR_EDR, SumR_GPP, SumR_SHAN, SumR_ALT, NeighG_SR, NeighG_EDG, NeighG_GPP, NeighG_SHAN, NeighG_ALT, NeighR_SR,NeighR_EDR,NeighR_GPP,NeighR_SHAN,NeighR_ALT, ncol = 5, nrow = 4) dev.off()