############################################################################################################################ ### Supporting Information ### # Title: Identifying factors that boost species discoveries of global reptiles # Authors: Jhonny J. M. Guedes1, Renato N. Feio1, Shai Meiri2 & Mario R. Moura3,4* # 1 Departamento de Biologia Animal, Universidade Federal de Viçosa, 36570-075, Viçosa, MG, Brazil # 2 School of Zoology, Tel Aviv University, 6997801, Tel Aviv, Israel # 3 Yale University, Department of Ecology and Evolutionary Biology, 06511, New Haven, CT, USA # 4 Departamento de Ciências Biológicas, Universidade Federal da Paraíba, 58397-000, Areia, PB, Brazil # * Corresponding author: mariormoura@gmail.com # DOI: 10.1093/zoolinnean/zlaa029. ############################################################################################################################ # STEPS IN THIS SCRIPT # 1. Load and understand the dataset. # 2. Prepare the variables that will be used, check for multicollinearity, non-normality, and standardise. # 3. Identify the best family distribution for the time-to-event models. # 4. Determine which factors (if any) affect the time lag in species discoveries and obtain the pseudo-R2 for the models. # 5. Same step as above, but accounting for lizards and snakes individually. # 6. Generate the Figure 1 - Global map of the holotype collection sites for squamate reptiles. # 7. Perform the Sensitivity Analysis to assess the effect of holotypes collected long ago. # 8. Generate subplot for the results of the Sensitivity Analysis. # 9. Repeat the Sensitivity Analysis for lizards only and create one subplot. # 10. Repeat the Sensitivity Analysis for snakes only and create the subplot. # 11. Create the multi panel plot (Figure 2 in main text) containing all three sensitivity plots above. # 12. Generate the Figure 3 - Proportion of reptiles described based on taxonomic reviews and holotypes collected by authors. # 13. Generate the Figure S1 presented in the supporting information. # 14. Generate the Figure S2 presented in the supporting information. # 15. Create the multi panel plot Figure S2. # 16. Generate the Figure S3 presented in the supporting information. # 17. Generate the Figure S4 presented in the supporting information. # Clean workspace rm(list=ls()); gc() # Install and load R packages needed to run the analysis: needed_packages<-c("survival", "flexsurv", "MuMIn", "usdm", "plyr", "data.table") new.packages<-needed_packages[!(needed_packages %in% installed.packages()[,"Package"])] if(length(new.packages)) install.packages(new.packages) lapply(needed_packages, require, character.only = TRUE) rm(needed_packages,new.packages) # Set the working directory: setwd() ##### ### TIME-TO-EVENT MODELS FOR SPECIES DESCRIPTIONS: # STEP 1 - Load, understand and prepare the dataset for analysis. ########################################################################################################################## # STEP 1 - Load, understand and prepare the dataset for analysis. # Load the response and predictor variables: trait_data <- read.csv("trait_data.csv", h=TRUE, stringsAsFactors = TRUE) ## We have 18 columns in the dataset, each of which is explained below: # Species: Self-explanatory. # Year_of_description: Year in which the species was formally described. # Year_hol: Year in which the holotype was collected. # TimeLag: Time lag in years between collection of the holotype and the species formal description. # Censor: indicates that discovery has occurred for all species in the dataset. # Log10_mass: Maximum known body mass log transformed. # Lat: Latitude, in absolute values, of the collection site of the holotype. # Lat_mz: Latitude, in absolute values, of the scientific collection in which the holotype is housed. # LogN_specimens_TS: Number of specimens within the type series log transformed. # LogN_spp_genus: Number of species in genus log transformed. # LogN_authors: Number of authors (taxonomists) per species descriptions log transformed. # Col_hol_is_author: Is the collector of the holotype an author of the species description? (Yes or No answer). # Taxonomic_review: Was the paper of description based on a Taxonomic review or not (Yes or No answer). # Molecular: Did the authors use molecular as evidence for the species description? (Yes or No answer). # Group: Reptile group to which a species belong (Lizards, Snakes or Turtles). # Biog_realm: Biogeographic regions to which a species belong (Afrotropic, Australasia, IndoMalay, Neartic, Neotropic, Oceania, or Palearctic). # Raw_long: Longitude of the collection site of the holotype in its raw value (needed for plotting figure 1). # Raw_lat: Latitude of the collection site of the holotype in its raw value (needed for plotting figure 1). ##### # STEP 2 - Prepare the variables that will be used, check for multicollinearity, non-normality, and standardise. ########################################################################################################################## # STEP 2 - Prepare the variables that will be used, check for multicollinearity, non-normality, and standardise. # Subset the continuous predictors that will be tested: names(trait_data) continuous_predictors<-trait_data[, c(6:11)] summary(continuous_predictors) # Standardised the predictors to make them comparable (median = 0, SD = 1): continuous_predictors<-scale(continuous_predictors, center=T, scale=T) continuous_predictors<-as.data.frame(continuous_predictors) # Check the Pearson correlation among predictors: corr_matrix<-cor(continuous_predictors[complete.cases(continuous_predictors[,]),], method="pearson") corr_matrix[corr_matrix == 1] <- 0 # replace correlation values of '1' (autocorrelation) by 'NA' range(corr_matrix) # Conclusion: continuous predictors show low correlation # Check for multicollinearity: usdm::vif(continuous_predictors) # variation inflation factor # Conclusion: continuous predictors show low multicollinearity (keep them all) # Subset the categorical predictors that will be tested, and add the group and biogeographic variables for filtering: names(trait_data) categorical_predictors <- trait_data[,c(12:14)] summary(categorical_predictors) # Merge the predictors and response variables in a single dataset: dataset <- cbind(trait_data[,c(1, 4:5, 15:16)], continuous_predictors, categorical_predictors) dataset$TimeLag <- dataset$TimeLag+1 # add +1 to the time lag to avoid erros when running the AFT model dataset <- dataset[, c(1:3,6:14,4:5)] summary(dataset) # Merge amphisbaenians with lizards for the taxon-specific analysis, and exclude turtles due to small sample size (im both datasets) # Removing turtles: dataset <- dataset[ !(dataset$Group %in% "Turtles"), ]; dataset <- droplevels(dataset) trait_data <- trait_data[ !(trait_data$Group %in% "Turtles"), ]; trait_data <- droplevels(trait_data) # Merging lizards and amphisbaenians: levels(trait_data$Group) levels(trait_data$Group)[3] <- "Lizards" levels(dataset$Group)[3] <- "Lizards" # Remove the rows (species) with NA values for further analysis: dataset <- dataset[complete.cases(dataset[,]),] dataset$Species <- droplevels(dataset$Species) summary(dataset) # Clean up the Global Environment rm(categorical_predictors,continuous_predictors,corr_matrix) ##### # STEP 3 - Identify the best family distribution for the time-to-event models. ########################################################################################################################## # STEP 3 - Identify the best family distribution for the time-to-event models. # Identify the best error distribution to be used for the Accelerated Failure Time (AFT) Model: Model_fit<-as.data.frame(matrix(ncol=3, nrow=6)) names(Model_fit)<-(c("Distrib", "npars", "AICc")) Model_fit[,1]<-c("exponential", "weibull", "lognormal", "loglogistic", "gamma", "gompertz") # define the error distributions Model_fit[1,2]<-(flexsurvreg(Surv(TimeLag, Censor)~1, dist="exp", data=dataset))$npars # get the numbers of parameters Model_fit[1,3]<-AICc(flexsurvreg(Surv(TimeLag, Censor)~1, dist="exp", data=dataset)) # get the AIC of the model Model_fit[2,2]<-(flexsurvreg(Surv(TimeLag, Censor)~1, dist="weibull", data=dataset))$npars Model_fit[2,3]<-AICc(flexsurvreg(Surv(TimeLag, Censor)~1, dist="weibull", data=dataset)) Model_fit[3,2]<-(flexsurvreg(Surv(TimeLag, Censor)~1, dist="lnorm", data=dataset))$npars Model_fit[3,3]<-AICc(flexsurvreg(Surv(TimeLag, Censor)~1, dist="lnorm", data=dataset)) Model_fit[4,2]<-(flexsurvreg(Surv(TimeLag, Censor)~1, dist="llogis", data=dataset))$npars Model_fit[4,3]<-AICc(flexsurvreg(Surv(TimeLag, Censor)~1, dist="llogis", data=dataset)) Model_fit[5,2]<-(flexsurvreg(Surv(TimeLag, Censor)~1, dist="gamma", data=dataset))$npars Model_fit[5,3]<-AICc(flexsurvreg(Surv(TimeLag, Censor)~1, dist="gamma", data=dataset)) Model_fit[6,2]<-(flexsurvreg(Surv(TimeLag, Censor)~1, dist="gompertz", data=dataset))$npars Model_fit[6,3]<-AICc(flexsurvreg(Surv(TimeLag, Censor)~1, dist="gompertz", data=dataset)) # Verify which error distribution best fits the data: Model_fit$deltaAICc<-Model_fit[,3]-min(Model_fit$AICc, na.rm=T) Model_fit$wAICc<-Weights(Model_fit$AICc) Model_fit[order(Model_fit$AICc, decreasing=F),] #write.csv(Model_fit[order(Model_fit$AICc, decreasing=F),],"Error distribution by AICc.csv", row.names = F) # Visual inspection of the top-4 best fitted models: par(mfrow=c(2,2)) plot((flexsurvreg(Surv(TimeLag, Censor)~1, dist="lognormal", data=dataset)), ylab="Species to be described", xlab="Time since collection", main="lognormal") plot((flexsurvreg(Surv(TimeLag, Censor)~1, dist="llogis", data=dataset)), ylab="Species to be described", xlab="Time since collection", main="llogis") plot((flexsurvreg(Surv(TimeLag, Censor)~1, dist="gompertz", data=dataset)), ylab="Species to be described", xlab="Time since collection", main="gompertz") plot((flexsurvreg(Surv(TimeLag, Censor)~1, dist="weibull", data=dataset)), ylab="Species to be described", xlab="Time since collection", main="weibull") par(mfrow=c(1,1)) # CONCLUSION: the 'lognormal' distribution shows the best fit to the model residuals. rm(Model_fit) ##### # STEP 4 - Determine which factors affect the time lag in species discoveries of global reptiles. ########################################################################################################################## # STEP 4 - Determine which factors affect the time lag in species discoveries of global reptiles. # Set a global model including the five continuous and three binary categorical predictors. globalMod <- survreg(Surv(TimeLag, Censor)~Log10_mass+Lat+Lat_mz+LogN_specimens_TS+ LogN_spp_genus+LogN_authors+Col_hol_is_author+Taxonomic_review+Molecular, dataset, dist="lognorm", na.action = na.fail) # Performs automated model selection based on subsets of the supplied global model. fits <- dredge(globalMod, trace = TRUE) # Coerce the object 'fits' into a data frame names(fits) sel.table<-as.data.frame(fits)[11:15] # The number of parameters (df) should be K by convention names(sel.table)[1] = "K" head(sel.table) # CONCLUSION: several models holding similar level of support. Therefore, model averaging is necessary. # Get the model averaged coefficients library(dplyr) library(knitr) avgmod <- model.avg(fits, revised.var = TRUE) ma_coefs <- coefTable(avgmod, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(Predictors, Estimate, `Std. Error`, lower95, upper95) kable(ma_coefs, digits = 3) # Save the result as html format; you may need to load the packages "kableExtra" and "magrittr" kable(ma_coefs, digits = 3, "html") %>% cat(., file = "global model avg.html") # Obtain the Pseudo-R2 from the model above: ma_coefs$Predictors row.names(ma_coefs) <- ma_coefs[,1] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","Log10_mass","Lat","Lat_mz","LogN_specimens_TS","LogN_spp_genus", "LogN_authors","Col_hol_is_author","Taxonomic_review","Molecular"),] # Predict the survival probabilities for each species based on the model coefficient: pct<-seq(0, 1, by=0.002) TimeLag_per_spp<-as.data.frame(matrix(nrow=length(pct), ncol=nrow(dataset))) names(TimeLag_per_spp)<-dataset$Species rownames(TimeLag_per_spp)<-(1-pct) for (k in 1:nrow(dataset)){ TimeLag_per_spp[,k]<-qlnorm(pct, meanlog=ma_coefs$Estimate[[1]]+ dataset$Log10_mass[k]*ma_coefs$Estimate[[3]]+ dataset$Lat[k]*ma_coefs$Estimate[[4]]+ dataset$Lat_mz[k]*ma_coefs$Estimate[[5]]+ dataset$LogN_specimens_TS[k]*ma_coefs$Estimate[[6]]+ dataset$LogN_spp_genus[k]*ma_coefs$Estimate[[7]]+ dataset$LogN_authors[k]*ma_coefs$Estimate[[8]]+ dataset$Col_hol_is_author[k]*ma_coefs$Estimate[[9]]+ dataset$Taxonomic_review[k]*ma_coefs$Estimate[[10]]+ dataset$Molecular[k]*ma_coefs$Estimate[[11]], sdlog=exp(ma_coefs$Estimate[[2]])) } # end of k for loop TimeLag_predicted <- as.data.frame(t(TimeLag_per_spp[500,])) colnames(TimeLag_predicted) <- "PredTimeLag" TimeLag_predicted$Species <- row.names(TimeLag_predicted) # Join the dataset and the predicted time lag for each species dataset<-join(dataset, TimeLag_predicted, by="Species", type="left", match="first") # Get the Pseudo-R2: cor(dataset$TimeLag, dataset$PredTimeLag, method = "pearson")^2 # Remove the column with the predicted time lag dataset <- dataset[,-15] # Clean the workspace rm(list=setdiff(ls(),c("dataset", "trait_data"))); gc() ##### # STEP 5 - Same step as above, but accounting for each reptile group individually. ########################################################################################################################## # STEP 5 - Same step as above, but accounting for each reptile group individually. # Prepare subdatasets: levels(dataset$Group) subds_lizards <- subset(dataset, Group=="Lizards", select=Species:Molecular) subds_snakes <- subset(dataset, Group=="Snakes", select=Species:Molecular) # Set a full model formula full_model_formula<-as.formula(Surv(TimeLag, Censor)~Log10_mass + Lat + Lat_mz + LogN_specimens_TS + LogN_spp_genus + LogN_authors + Col_hol_is_author + Taxonomic_review + Molecular) # Lizards # Set a global model globalMod <- survreg(full_model_formula, subds_lizards, dist="lognorm", na.action = na.fail) # Performs automated model selection with subsets of the supplied global model. fits <- dredge(globalMod, trace = TRUE) # Coerce the object 'fits' into a dataframe sel.table<-as.data.frame(fits)[11:15] # The number of parameters (df) should be K by convention names(sel.table)[1] = "K" head(sel.table) # CONCLUSION: several models holding similar level of support. Therefore, model averaging is necessary. # Get the model averaged coefficients: avgmod <- model.avg(fits, revised.var = TRUE) ma_coefs <- coefTable(avgmod, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(Predictors, Estimate, `Std. Error`, lower95, upper95) kable(ma_coefs, digits = 3) # Save the result: kable(ma_coefs, digits = 3, "html") %>% cat(., file = "model avg lizards.html") # Obtain the Pseudo-R2 from the model above: ma_coefs$Predictors row.names(ma_coefs) <- ma_coefs[,1] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","Log10_mass","Lat","Lat_mz","LogN_specimens_TS","LogN_spp_genus", "LogN_authors","Col_hol_is_author","Taxonomic_review","Molecular"),] # Predict the survival probabilities for each species based on the model coefs: pct<-seq(0, 1, by=0.002) TimeLag_per_spp<-as.data.frame(matrix(nrow=length(pct), ncol=nrow(subds_lizards))) names(TimeLag_per_spp)<-subds_lizards$Species rownames(TimeLag_per_spp)<-(1-pct) for (k in 1:nrow(subds_lizards)){ TimeLag_per_spp[,k]<-qlnorm(pct, meanlog=ma_coefs$Estimate[[1]]+ subds_lizards$Log10_mass[k]*ma_coefs$Estimate[[3]]+ subds_lizards$Lat[k]*ma_coefs$Estimate[[4]]+ subds_lizards$Lat_mz[k]*ma_coefs$Estimate[[5]]+ subds_lizards$LogN_specimens_TS[k]*ma_coefs$Estimate[[6]]+ subds_lizards$LogN_spp_genus[k]*ma_coefs$Estimate[[7]]+ subds_lizards$LogN_authors[k]*ma_coefs$Estimate[[8]]+ subds_lizards$Col_hol_is_author[k]*ma_coefs$Estimate[[9]]+ subds_lizards$Taxonomic_review[k]*ma_coefs$Estimate[[10]]+ subds_lizards$Molecular[k]*ma_coefs$Estimate[[11]], sdlog=exp(ma_coefs$Estimate[[2]])) } # end of k for loop TimeLag_predicted <- as.data.frame(t(TimeLag_per_spp[500,])) colnames(TimeLag_predicted) <- "PredTimeLag" TimeLag_predicted$Species <- row.names(TimeLag_predicted) # Join the dataset and the predicted time lag for each species subds_lizards<-join(subds_lizards, TimeLag_predicted, by="Species", type="left", match="first") # Get the Pseudo-R2: cor(subds_lizards$TimeLag, subds_lizards$PredTimeLag, method = "pearson")^2 # Clean the workspace rm(list=setdiff(ls(),c("dataset", "trait_data", "subds_snakes", "full_model_formula"))); gc() # Snakes # Set a global model globalMod <- survreg(full_model_formula, subds_snakes, dist="lognorm", na.action = na.fail) # Performs automated model selection with subsets of the supplied global model. fits <- dredge(globalMod, trace = TRUE) # Coerce the object 'fits' into a dataframe names(fits) sel.table<-as.data.frame(fits)[11:15] # The number of parameters (df) should be K by convention names(sel.table)[1] = "K" head(sel.table) # CONCLUSION: several models holding similar level of support. Therefore, model averaging is necessary. # Get the model averaged coefficients: avgmod <- model.avg(fits, revised.var = TRUE) ma_coefs <- coefTable(avgmod, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(Predictors, Estimate, `Std. Error`, lower95, upper95) kable(ma_coefs, digits = 3) # Save the result: kable(ma_coefs, digits = 3, "html") %>% cat(., file = "model avg snakes.html") # Obtain the Pseudo-R2 from the model above: ma_coefs$Predictors row.names(ma_coefs) <- ma_coefs[,1] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","Log10_mass","Lat","Lat_mz","LogN_specimens_TS","LogN_spp_genus", "LogN_authors","Col_hol_is_author","Taxonomic_review","Molecular"),] # Predict the survival probabilities for each species based on the model coefs: pct<-seq(0, 1, by=0.002) TimeLag_per_spp<-as.data.frame(matrix(nrow=length(pct), ncol=nrow(subds_snakes))) names(TimeLag_per_spp)<-subds_snakes$Species rownames(TimeLag_per_spp)<-(1-pct) for (k in 1:nrow(subds_snakes)){ TimeLag_per_spp[,k]<-qlnorm(pct, meanlog=ma_coefs$Estimate[[1]]+ subds_snakes$Log10_mass[k]*ma_coefs$Estimate[[3]]+ subds_snakes$Lat[k]*ma_coefs$Estimate[[4]]+ subds_snakes$Lat_mz[k]*ma_coefs$Estimate[[5]]+ subds_snakes$LogN_specimens_TS[k]*ma_coefs$Estimate[[6]]+ subds_snakes$LogN_spp_genus[k]*ma_coefs$Estimate[[7]]+ subds_snakes$LogN_authors[k]*ma_coefs$Estimate[[8]]+ subds_snakes$Col_hol_is_author[k]*ma_coefs$Estimate[[9]]+ subds_snakes$Taxonomic_review[k]*ma_coefs$Estimate[[10]]+ subds_snakes$Molecular[k]*ma_coefs$Estimate[[11]], sdlog=exp(ma_coefs$Estimate[[2]])) } # end of k for loop TimeLag_predicted <- as.data.frame(t(TimeLag_per_spp[500,])) colnames(TimeLag_predicted) <- "PredTimeLag" TimeLag_predicted$Species <- row.names(TimeLag_predicted) # Join the dataset and the predicted time lag for each species subds_snakes<-join(subds_snakes, TimeLag_predicted, by="Species", type="left", match="first") # Get the Pseudo-R2: cor(subds_snakes$TimeLag, subds_snakes$PredTimeLag, method = "pearson")^2 # Clean the workspace rm(list=setdiff(ls(),c("dataset", "trait_data"))); gc() ##### # STEP 6 - Generate the Figure 1 - Global map of the holotype collection sites for squamate reptiles. ########################################################################################################################## # STEP 6 - Generate the Figure 1 - Global map of the holotype collection sites for squamate reptiles. # Generate the figures presented in the paper: # FIGURE 1 - Global distribution of holotypes of reptile species described from 1992 to 2017, and their time lags # First, create the world map containing the holotype distribution data. rm(list=setdiff(ls(),c("dataset", "trait_data"))); gc() # Install and load required packages: needed_packages<-c("rgdal", "sp", "raster", "maptools", "maps", "dplyr", "plyr", "ggplot2") new.packages<-needed_packages[!(needed_packages %in% installed.packages()[,"Package"])] if(length(new.packages)) install.packages(new.packages) lapply(needed_packages, require, character.only = TRUE) # Load holotypes distribution data, and prepare it for plotting: names(trait_data) site_descriptors <- trait_data[, c(1,4,15:18)] # this dataset contains raw longitude and latitude (not absolute values) site_descriptors$Raw_long <- as.numeric(as.character(site_descriptors$Raw_long)) site_descriptors$Raw_lat <- as.numeric(as.character(site_descriptors$Raw_lat)) site_descriptors <- site_descriptors[order(site_descriptors$Group, decreasing = F),] # group with less species will be plotted on top site_descriptors <- site_descriptors[complete.cases(site_descriptors[, 5:6]), ] # remove those species lacking longitude and latitude data coords<-site_descriptors[, c(5:6)] # dataframe with only longitude and latitude of holotypes # convert the object into a SpatialPointsDataframe site_descriptors<-SpatialPointsDataFrame(coords=coords, data=site_descriptors, proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) crs(site_descriptors) # Transform native projection of points into CEA (cylindrical equal area) projection cea<-"+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" site_descriptors<-spTransform(site_descriptors, cea) site_descriptors@data$id<-rownames(site_descriptors@data) site_descriptors_df <- data.frame(site_descriptors) extent(site_descriptors) # Download the shapefile of the biogeographic realms at . # Please, unzip it in a folder called 'shapefile' at your working directory. # Read the shapefile; transform them into CEA Projection; prepare them for ggplot2: terr_cover<-readOGR(dsn="shapefile", layer='wwf_terr_ecos') # change directory as needed crs(terr_cover) # Retrieve Coordinate Reference System From Object terr_cover<-spTransform(terr_cover, cea) # Transform native projection terr_cover@data$id<-rownames(terr_cover@data) terr_cover_df<-fortify(terr_cover, region="id") terr_cover_df<-join(terr_cover_df, terr_cover@data, by="id") # joining tables by "id" crs(terr_cover) extent(terr_cover) # Define colors to be used for each reptile group myColors <- c("#fdae61", "#454fff") names(myColors)<-levels(site_descriptors_df$Group) plot <- ggplot(data=site_descriptors_df) + # dataset to plot geom_polygon(data=terr_cover_df, aes(long, lat, group=group), colour=NA, fill="grey90", alpha = .5, size=0.2, show.legend = F) + geom_point(aes(x=Raw_long.1, y=Raw_lat.1, colour=Group, shape=Group), na.rm = T) + scale_shape_manual(values=c(4,2)) + scale_colour_manual(name = "Group", values = myColors) + #geom_polygon(data=terr_cover_df, aes(long, lat, group=group), colour = "transparent", fill="transparent", alpha = .1, size=0.2) + coord_equal() + theme(axis.line=element_blank(), axis.text=element_blank(), axis.ticks=element_blank(), axis.title=element_blank(), panel.grid.minor=element_blank(), panel.grid.major=element_blank(), legend.position=c(.5, 1), legend.direction="horizontal", legend.text=element_text(size=10), legend.key = element_blank(), legend.background=element_blank(), legend.title=element_blank(), legend.justification = "center", panel.background=element_blank(), plot.margin=unit(c(0.1,0.1,0.1,0.1), "lines")) ggsave("test.pdf", width=17.8, height=12.5, units="cm") #ggsave("holotype_map.pdf", width=17.8, height=12.5, units="cm") # Now, create the boxplots based on time lags across biogeographic realms and reptile groups. # Boxplots were added to the map with the program inkscape. rm(list=setdiff(ls(),c("dataset", "trait_data"))); gc() # Load required package library(ggplot2) # Define colors to be used for each reptile group myColors <- c("#fdae61", "#454fff") names(myColors)<-levels(trait_data$Group) # Folder where you want the graphs to be saved: results <- "WRITE HERE/" # Create a list of biogreographic realms in the data to loop over biog_list <- unique(trait_data$Biog_realm) # Create for loop to produce ggplot2 graphs for(i in seq_along(biog_list)){ plot <- ggplot(subset(trait_data, trait_data$Biog_realm==biog_list[[i]]), aes(x=Group, y=TimeLag, fill=Group))+ geom_boxplot(outlier.size=0.5, outlier.alpha = .7, na.rm = T, colour="black", size = .3, fatten = .8)+ labs(x = paste(biog_list[[i]]), y = NULL)+ scale_fill_manual(name = "Group", values = myColors)+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(face = "bold", size = 14), axis.line = element_line(colour="black"), axis.text = element_blank(), axis.text.y = element_text(size = 12, colour = "black"), axis.ticks.x = element_blank(), plot.background=element_rect(fill="transparent", colour=NA), legend.position="none") # save plots as .pdf ggsave(plot, file=paste(results, 'TL_', biog_list[i], ".pdf", sep=''), width=1.5, height = 2, units="in", dpi = "print", bg = "transparent") print(plot) } ##### # STEP 7 - Perform the Sensitivity Analysis to assess the effect of holotypes collected long ago. ########################################################################################################################## # STEP 7 - Perform the Sensitivity Analysis to assess the effect of holotypes collected long ago. rm(list=setdiff(ls(),c("trait_data"))); gc() # Install and load R packages needed to run the analysis: needed_packages<-c("survival","flexsurv","MuMIn","usdm","plyr","data.table","dplyr","knitr") new.packages<-needed_packages[!(needed_packages %in% installed.packages()[,"Package"])] if(length(new.packages)) install.packages(new.packages) lapply(needed_packages, require, character.only = TRUE) # Subset the continuous predictors that will be tested: names(trait_data) continuous_predictors<-trait_data[, c(6:11)] summary(continuous_predictors) # Standardised the predictors to make them comparable (median = 0, SD = 1): continuous_predictors<-scale(continuous_predictors, center=T, scale=T) continuous_predictors<-as.data.frame(continuous_predictors) # Subset the categorical predictors that will be tested, and add the group and biogeographic variables for filtering: names(trait_data) categorical_predictors <- trait_data[,c(12:14)] summary(categorical_predictors) # Merge the predictors and response variables in a single dataset: dataset <- cbind(trait_data[,c(1, 3:5, 15)], continuous_predictors, categorical_predictors) dataset$TimeLag <- dataset$TimeLag+1 # add +1 to the time lag to avoid erros when running the AFT model dataset <- dataset[, c(1:4,6:14,5)] summary(dataset) ## Remove the rows (species) with NA values: dataset <- dataset[complete.cases(dataset[,]),] dataset$Species <- droplevels(dataset$Species) summary(dataset) # Clean up the Global Environment rm(categorical_predictors,continuous_predictors, new.packages, needed_packages) # Generate the three sensitivity plots to be combined and thus, create Figure 2 # 1st - Sensitivity analysis for the global model (all reptiles) # Create an object to 'guide' the for loop: collection_dates<-seq(from=1952, to=1992, by=5) # Set a full model formula full_model_formula<-as.formula(Surv(TimeLag, Censor) ~ Log10_mass + Lat + Lat_mz + LogN_specimens_TS + LogN_spp_genus + LogN_authors + Col_hol_is_author + Taxonomic_review + Molecular) # Create an empty list to store the outputs: my_outputs<-list() # Loop over each time period: for (i in 1:length(collection_dates)){ # i = each one of the time periods # Create subsets based on species collected in a given time period data_subset <- dataset[which(dataset$Year_hol>=collection_dates[i]), ] fits <- dredge(survreg(full_model_formula, data_subset, dist="lognorm", na.action = na.fail)) avgmodel_output <- model.avg(fits, revised.var = T) ma_coefs <- coefTable(avgmodel_output, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, TimePeriod = paste(i), Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(TimePeriod, Predictors, Estimate, `Std. Error`, lower95, upper95) row.names(ma_coefs) <- ma_coefs[,2] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","LogN_authors","LogN_specimens_TS","LogN_spp_genus", "Log10_mass","Lat","Lat_mz","Col_hol_is_author","Taxonomic_review","Molecular"),] print(i) my_outputs[[i]]<-ma_coefs # save the coefficients } # end of the for loop i # Reorder levels of 'Time Period' variable for plotting my_outputs <- bind_rows(my_outputs) my_outputs$TimePeriod <- factor(my_outputs$TimePeriod, labels = c('1952 to 2017','1957 to 2017','1962 to 2017','1967 to 2017','1972 to 2017', '1977 to 2017','1982 to 2017','1987 to 2017','1992 to 2017')) # Get coefficients for the full time period (all reptiles) fits <- dredge(survreg(full_model_formula, dataset, dist="lognorm", na.action = na.fail), trace = TRUE) avgmod <- model.avg(fits, revised.var = TRUE) ma_coefs <- coefTable(avgmod, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, TimePeriod = "Full time period", Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(TimePeriod, Predictors, Estimate, `Std. Error`, lower95, upper95) # Combine 'my_outputs' and 'ma_coefs' SensitAnalysis_data <- rbind(my_outputs, ma_coefs) SensitAnalysis_data <- SensitAnalysis_data[SensitAnalysis_data$Predictors != "(Intercept)", ] SensitAnalysis_data <- SensitAnalysis_data[SensitAnalysis_data$Predictors != "Log(scale)", ] # write.csv(SensitAnalysis_data, "SensitAnalysis_data.csv", row.names = F) str(SensitAnalysis_data) levels(SensitAnalysis_data$Predictors) # Reorder levels to plot and add labels: SensitAnalysis_data$Predictors <- factor(SensitAnalysis_data$Predictors, levels=c("Log10_mass","Lat","Lat_mz","LogN_authors","LogN_specimens_TS", "LogN_spp_genus","Col_hol_is_author","Taxonomic_review","Molecular"), labels=c("Body size","Lat of holotype","Lat of museum","N of authors/spp", "N of type-specimens","N of spp/genus","Collector is author", "Taxonomic review","Molecular analysis")) # Reorder levels to plot and add labels: levels(SensitAnalysis_data$TimePeriod) SensitAnalysis_data$TimePeriod <- factor(SensitAnalysis_data$TimePeriod, levels=c("Full time period","1952 to 2017", "1957 to 2017", "1962 to 2017", "1967 to 2017", "1972 to 2017", "1977 to 2017","1982 to 2017", "1987 to 2017", "1992 to 2017")) ##### # STEP 8 - Generate subplot for the results of the Sensitivity Analysis. ########################################################################################################################## # STEP 8 - Generate subplot for the results of the Sensitivity Analysis. rm(list=setdiff(ls(),c("trait_data"))); gc() library(ggplot2) library(viridis) # Define colors to be used in the plot myColors <- viridis_pal(option="plasma")(10) names(myColors)<-levels(SensitAnalysis_data$TimePeriod) (Sensitivity_plot <- ggplot(SensitAnalysis_data, aes(x = TimePeriod, y = Estimate, ymin = lower95, ymax = upper95))+ geom_pointrange(aes(col = TimePeriod, shape = TimePeriod), size = 0.4, na.rm = T)+ scale_colour_manual(name = "TimePeriod", values = myColors)+ scale_shape_manual(values=c(0,1,2,4,5,6,7,8,9,10))+ geom_errorbar(aes(ymin=lower95, ymax=upper95, col = TimePeriod), width=0.1, na.rm = T)+ geom_hline(yintercept =0, linetype=2)+ labs(x="Predictors", y=NULL)+ scale_x_discrete(limits = rev(levels(SensitAnalysis_data$TimePeriod)))+ facet_wrap(~Predictors, strip.position="left", nrow=9)+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=12, face="bold"), axis.text = element_text(size=10), axis.line = element_line(colour="black"), axis.ticks.y.left = element_blank(), axis.text.x=element_text(face="italic"), axis.text.y = element_blank(), plot.background=element_rect(fill = "white"), strip.background = element_blank(), strip.placement = "outside", strip.text.y = element_text(hjust=1, vjust=0.5, angle=180, size = 10), legend.key = element_blank(), legend.title = element_blank(), legend.position = "none")+ coord_flip()) rm(list=setdiff(ls(),c("dataset", "trait_data", "Sensitivity_plot"))); gc() ##### # STEP 9 - Repeat the Sensitivity Analysis for lizards only and create one subplot. ########################################################################################################################## # STEP 9 - Repeat the Sensitivity Analysis for lizards only and create the plot. # Select lizards only levels(dataset$Group) df_lizards <- subset(dataset, Group == "Lizards", select = Species:Molecular) # Create an object to 'guide' the for loop: collection_dates<-seq(from=1952, to=1992, by=5) # Set a full model formula full_model_formula<-as.formula(Surv(TimeLag, Censor) ~ Log10_mass + Lat + Lat_mz + LogN_specimens_TS + LogN_spp_genus + LogN_authors + Col_hol_is_author + Taxonomic_review + Molecular) # Create an empty list to store the outputs: my_outputs<-list() # Loop over each time period: for (i in 1:length(collection_dates)){ # i = each one of the time periods # Create subsets based on species collected in a given time period data_subset <- df_lizards[which(df_lizards$Year_hol>=collection_dates[i]), ] fits <- dredge(survreg(full_model_formula, data_subset, dist="lognorm", na.action = na.fail)) avgmodel_output <- model.avg(fits, revised.var = T) ma_coefs <- coefTable(avgmodel_output, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, TimePeriod = paste(i), Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(TimePeriod, Predictors, Estimate, `Std. Error`, lower95, upper95) row.names(ma_coefs) <- ma_coefs[,2] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","LogN_authors","LogN_specimens_TS","LogN_spp_genus", "Log10_mass","Lat","Lat_mz","Col_hol_is_author","Taxonomic_review","Molecular"),] print(i) my_outputs[[i]]<-ma_coefs } # end of the for loop i my_outputs <- bind_rows(my_outputs) my_outputs$TimePeriod <- factor(my_outputs$TimePeriod, labels = c('1952 to 2017','1957 to 2017','1962 to 2017','1967 to 2017','1972 to 2017', '1977 to 2017','1982 to 2017','1987 to 2017','1992 to 2017')) # Get coefficients for the full time period (all lizards) fits <- dredge(survreg(full_model_formula, df_lizards, dist="lognorm", na.action = na.fail), trace = TRUE) avgmod <- model.avg(fits, revised.var = TRUE) ma_coefs <- coefTable(avgmod, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, TimePeriod = "Full time period", Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(TimePeriod, Predictors, Estimate, `Std. Error`, lower95, upper95) SensitAnalysis_lizards <- rbind(my_outputs, ma_coefs) SensitAnalysis_lizards <- SensitAnalysis_lizards[SensitAnalysis_lizards$Predictors != "(Intercept)", ] SensitAnalysis_lizards <- SensitAnalysis_lizards[SensitAnalysis_lizards$Predictors != "Log(scale)", ] #write.csv(SensitAnalysis_lizards, "SensitAnalysis_lizards.csv", row.names = F) str(SensitAnalysis_lizards) # Reorder levels to plot and add labels: SensitAnalysis_lizards$Predictors <- factor(SensitAnalysis_lizards$Predictors, levels=c("Log10_mass","Lat","Lat_mz","LogN_authors","LogN_specimens_TS", "LogN_spp_genus","Col_hol_is_author","Taxonomic_review","Molecular"), labels=c("Body size","Lat of holotype","Lat of museum","N of authors/spp", "N of type-specimens","N of spp/genus","Collector is author", "Taxonomic review","Molecular analysis")) # Reorder levels to plot and add labels: levels(SensitAnalysis_lizards$TimePeriod) SensitAnalysis_lizards$TimePeriod <- factor(SensitAnalysis_lizards$TimePeriod, levels=c("Full time period","1952 to 2017", "1957 to 2017", "1962 to 2017", "1967 to 2017", "1972 to 2017", "1977 to 2017","1982 to 2017", "1987 to 2017", "1992 to 2017")) library(ggplot2) library(viridis) # Define colors to be used in the plot myColors <- viridis_pal(option="plasma")(10) names(myColors)<-levels(SensitAnalysis_lizards$TimePeriod) (Sensitivity_plot_lizards <- ggplot(SensitAnalysis_lizards, aes(x = TimePeriod, y = Estimate, ymin = lower95, ymax = upper95))+ geom_pointrange(aes(col = TimePeriod, shape = TimePeriod), size = 0.4)+ scale_colour_manual(name = "TimePeriod", values = myColors)+ scale_shape_manual(values=c(0,1,2,4,5,6,7,8,9,10))+ geom_errorbar(aes(ymin=lower95, ymax=upper95, col = TimePeriod), width=0.1)+ geom_hline(yintercept =0, linetype=2)+ labs(x=NULL, y=NULL)+ scale_x_discrete(limits = rev(levels(SensitAnalysis_lizards$TimePeriod)))+ facet_wrap(~Predictors, strip.position="left", nrow=9)+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=12, face="bold"), axis.text = element_text(size=10), axis.line = element_line(colour="black"), axis.line.y = element_blank(), axis.ticks.y.left = element_blank(), axis.ticks.y=element_blank(), axis.text.x=element_text(face="italic"), axis.text.y = element_blank(), plot.background=element_rect(fill = "white"), strip.background = element_blank(), strip.placement = "outside", strip.text.y = element_blank(), legend.title = element_blank(), legend.position="none")+ coord_flip()) rm(list=setdiff(ls(),c("dataset", "trait_data", "Sensitivity_plot", "Sensitivity_plot_lizards"))); gc() ##### # STEP 10 - Repeat the Sensitivity Analysis for snakes only and create the subplot. ########################################################################################################################## # STEP 10 - Repeat the Sensitivity Analysis for snakes only and create the plot. # Select snakes only levels(dataset$Group) df_snakes <- subset(dataset, Group == "Snakes", select = Species:Molecular) # Create an object to 'guide' the for loop: collection_dates<-seq(from=1952, to=1992, by=5) # Set a full model formula full_model_formula<-as.formula(Surv(TimeLag, Censor) ~ Log10_mass + Lat + Lat_mz + LogN_specimens_TS + LogN_spp_genus + LogN_authors + Col_hol_is_author + Taxonomic_review + Molecular) # Create an empty list to store the outputs: my_outputs<-list() # Loop over each time period: for (i in 1:length(collection_dates)){ # i = each one of the time periods # Create subsets based on species collected in a given time period data_subset <- df_snakes[which(df_snakes$Year_hol>=collection_dates[i]), ] fits <- dredge(survreg(full_model_formula, data_subset, dist="lognorm", na.action = na.fail)) avgmodel_output <- model.avg(fits, revised.var = T) ma_coefs <- coefTable(avgmodel_output, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, TimePeriod = paste(i), Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(TimePeriod, Predictors, Estimate, `Std. Error`, lower95, upper95) row.names(ma_coefs) <- ma_coefs[,2] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","LogN_authors","LogN_specimens_TS","LogN_spp_genus", "Log10_mass","Lat","Lat_mz","Col_hol_is_author","Taxonomic_review","Molecular"),] print(i) my_outputs[[i]]<-ma_coefs } # end of the for loop i my_outputs <- bind_rows(my_outputs) my_outputs$TimePeriod <- factor(my_outputs$TimePeriod, labels = c('1952 to 2017','1957 to 2017','1962 to 2017','1967 to 2017','1972 to 2017', '1977 to 2017','1982 to 2017','1987 to 2017','1992 to 2017')) # Get coefficients for the full time period (all snakes) fits <- dredge(survreg(full_model_formula, df_snakes, dist="lognorm", na.action = na.fail), trace = TRUE) avgmod <- model.avg(fits, revised.var = TRUE) ma_coefs <- coefTable(avgmod, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, TimePeriod = "Full time period", Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(TimePeriod, Predictors, Estimate, `Std. Error`, lower95, upper95) SensitAnalysis_snakes <- rbind(my_outputs, ma_coefs) SensitAnalysis_snakes <- SensitAnalysis_snakes[SensitAnalysis_snakes$Predictors != "(Intercept)", ] SensitAnalysis_snakes <- SensitAnalysis_snakes[SensitAnalysis_snakes$Predictors != "Log(scale)", ] #write.csv(SensitAnalysis_snakes, "SensitAnalysis_snakes.csv", row.names = F) str(SensitAnalysis_snakes) # Reorder levels to plot and add labels: SensitAnalysis_snakes$Predictors <- factor(SensitAnalysis_snakes$Predictors, levels=c("Log10_mass","Lat","Lat_mz","LogN_authors","LogN_specimens_TS", "LogN_spp_genus","Col_hol_is_author","Taxonomic_review","Molecular"), labels=c("Body size","Lat of holotype","Lat of museum","N of authors/spp", "N of type-specimens","N of spp/genus","Collector is author", "Taxonomic review","Molecular analysis")) # Reorder levels to plot and add labels: levels(SensitAnalysis_snakes$TimePeriod) SensitAnalysis_snakes$TimePeriod <- factor(SensitAnalysis_snakes$TimePeriod, levels=c("Full time period","1952 to 2017", "1957 to 2017", "1962 to 2017", "1967 to 2017", "1972 to 2017", "1977 to 2017","1982 to 2017", "1987 to 2017", "1992 to 2017")) library(ggplot2) library(viridis) # Define colors to be used in the plot myColors <- viridis_pal(option="plasma")(10) names(myColors)<-levels(SensitAnalysis_snakes$TimePeriod) (Sensitivity_plot_snakes <- ggplot(SensitAnalysis_snakes, aes(x = TimePeriod, y = Estimate, ymin = lower95, ymax = upper95))+ geom_pointrange(aes(col = TimePeriod, shape = TimePeriod), size = 0.4)+ scale_colour_manual(name = "TimePeriod", values = myColors)+ scale_shape_manual(values=c(0,1,2,4,5,6,7,8,9,10))+ geom_errorbar(aes(ymin=lower95, ymax=upper95, col = TimePeriod), width=0.1)+ geom_hline(yintercept =0, linetype=2)+ labs(x=NULL, y=NULL)+ scale_x_discrete(limits = rev(levels(SensitAnalysis_snakes$TimePeriod)))+ facet_wrap(~Predictors, strip.position="left", nrow=9)+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=12, face="bold"), axis.text = element_text(size=10), axis.line = element_line(colour="black"), axis.line.y = element_blank(), axis.ticks.y.left = element_blank(), axis.ticks.y=element_blank(), axis.text.x=element_text(face="italic"), axis.text.y = element_blank(), plot.background=element_rect(fill = "white"), strip.background = element_blank(), strip.placement = "outside", strip.text.y = element_blank(), legend.key = element_blank(), legend.title = element_blank(), legend.text=element_text(size=8), legend.justification = "center", legend.direction = "vertical", legend.position=c(.85, .5))+ coord_flip()) rm(list=setdiff(ls(),c("dataset", "trait_data", "Sensitivity_plot", "Sensitivity_plot_lizards", "Sensitivity_plot_snakes"))); gc() ##### # STEP 11 - Create the multi panel plot (Figure 2 in main text) containing all three sensitivity plots above. ########################################################################################################################## # STEP 11 - Create the multi panel plot (Figure 2 in main text) containing all three sensitivity plots above. library("ggpubr") (Sensitivity_plot_combined <- ggarrange(Sensitivity_plot,Sensitivity_plot_lizards,Sensitivity_plot_snakes, ncol = 3, nrow = 1, labels = c('All reptiles','Lizards','Snakes'), label.x = c(0.4, 0.1, 0.1), label.y = c(.98,.98,.98), widths = c(1.4,1.0,1.2), font.label = list(size = 10, color = "black"), align = "h")) (Sensitivity_plot_final <- annotate_figure(Sensitivity_plot_combined, bottom = text_grob("Avg. Weighted Coefs. (IC95)", color = "black", face = "bold", size = 12, vjust = 0, hjust = 0))) myfile <- paste("~/WRITE HERE YOUR WORKING DIRECTORY/Sensitivity plot.pdf") myfile.png <- paste("~/WRITE HERE YOUR WORKING DIRECTORY/Sensitivity plot.png") ggsave(myfile, plot=Sensitivity_plot_final, width=10, height = 7, units="in", dpi = "print") ggsave(myfile.png, plot=Sensitivity_plot_final, width=10, height = 7, units="in", dpi = "print") rm(list=setdiff(ls(),c("trait_data"))); gc() ##### # STEP 12 - Generate the Figure 3 - Proportion of reptiles described based on taxonomic reviews and holotypes collected by authors. ########################################################################################################################## # STEP 12 - Generate the Figure 3 - Proportion of reptiles described based on taxonomic reviews and holotypes collected by authors. # Load required packages needed_packages<-c("ggplot2", "egg", "dplyr", "ggpubr", "gridExtra", "plyr", "grid", "forcats") new.packages<-needed_packages[!(needed_packages %in% installed.packages()[,"Package"])] if(length(new.packages)) install.packages(new.packages) lapply(needed_packages, require, character.only = TRUE) # Clean workspace rm(list = ls()); gc() # Load trait data: trait_data <- read.csv("trait_data.csv", h=TRUE, stringsAsFactors = TRUE) # Load the taxonomist_data and extract the family information: taxonomist_data <- read.csv("taxonomist_data.csv", h=T, stringsAsFactors = T, fileEncoding = "Latin1") taxonomist_data <- taxonomist_data[ , -c(3:5)] # Join trait data with taxonomist data matching by species: trait_data <- join(trait_data, taxonomist_data, by = "Species"); rm(taxonomist_data) # Remove species with missing data for who collected the holotype: trait_data <- droplevels(trait_data[!is.na(trait_data$Col_hol_is_author), ] ) # Prepare for plotting: # Order the bars according to the proportion of collectors not authors # Define colors to be used for each reptile group snake.col <- c("#abd9e9", "#454fff") names(snake.col)<-levels(trait_data$col_hol_is_author) lizard.col <- c("#fee090", "#fdae61") names(lizard.col)<-levels(trait_data$col_hol_is_author) fig.lizards <- trait_data[trait_data$Group=="Lizards",] %>% # convert variable to factor, ordered (descending) by the proportion of rows where class == "0" mutate(Family = fct_reorder(.f = Family, .x = as.factor(Col_hol_is_author), .fun = function(.x) mean(.x == "0"), .desc = TRUE)) %>% ggplot(aes(x = Family, fill = as.factor(Col_hol_is_author))) + geom_bar(position = "fill") + geom_hline(yintercept=.5, linetype="dashed", color = "red") + scale_fill_manual(name = "Col_hol_is_author", values = lizard.col) + scale_y_continuous(breaks = seq(0, 1, .25)) + labs(y = "% collectors are authors", x = "Lizard families") + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=8, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=6, colour = "black"), axis.text.x = element_text(angle = 45, hjust=1, vjust=1), # hjust=1, vjust = 0.5 legend.position ="none") fig.lizards fig.snakes <- trait_data[trait_data$Group=="Snakes",] %>% # convert variable to a factor, ordered (descending) by the proportion of rows where the class == "0" mutate(Family = fct_reorder(.f = Family, .x = as.factor(Col_hol_is_author), .fun = function(.x) mean(.x == "0"), .desc = TRUE)) %>% ggplot(aes(x = Family, fill = as.factor(Col_hol_is_author))) + geom_bar(position = "fill") + geom_hline(yintercept=.5, linetype="dashed", color = "red") + scale_y_continuous(breaks = seq(0, 1, .25)) + scale_fill_manual(name = "Col_hol_is_author", values = snake.col) + labs(y = "Snake families", x = NULL) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=8, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=6, colour = "black"), axis.line.y = element_blank(), axis.ticks.y.left = element_blank(), axis.ticks.y=element_blank(), axis.text.y = element_blank(), strip.text.y = element_blank(), axis.text.x = element_text(angle = 45, hjust=1, vjust=1), legend.position ="none") fig.snakes # Order the bars according to the proportion of TR:non-TR # Define colors to be used for each reptile group snake.col <- c("#abd9e9", "#454fff") names(snake.col)<-levels(trait_data$Taxonomic_review) lizard.col <- c("#fee090", "#fdae61") names(lizard.col)<-levels(trait_data$Taxonomic_review) # lizards fig.lizardsTR <- trait_data[trait_data$Group=="Lizards",] %>% # convert variable to factor, ordered (descending) by the proportion of rows where class == "0" mutate(Family = fct_reorder(.f = Family, .x = as.factor(Taxonomic_review), .fun = function(.x) mean(.x == "0"), .desc = TRUE)) %>% ggplot(aes(x = Family, fill = as.factor(Taxonomic_review))) + geom_bar(position = "fill") + geom_hline(yintercept=.5, linetype="dashed", color = "red") + scale_y_continuous(breaks = seq(0, 1, .25)) + scale_fill_manual(name = "Taxonomic_review", values = lizard.col) + labs(y = "% taxonomic reviews", x = NULL) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=8, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=6, colour = "black"), axis.text.x = element_text(angle = 45, hjust=1, vjust=1), legend.position ="none") fig.lizardsTR fig.snakesTR <- trait_data[trait_data$Group=="Snakes",] %>% # convert variable to a factor, ordered (descending) by the proportion of rows where the class == "0" mutate(Family = fct_reorder(.f = Family, .x = as.factor(Taxonomic_review), .fun = function(.x) mean(.x == "0"), .desc = TRUE)) %>% ggplot(aes(x = Family, fill = as.factor(Taxonomic_review))) + geom_bar(position = "fill") + geom_hline(yintercept=.5, linetype="dashed", color = "red") + scale_y_continuous(breaks = seq(0, 1, .25)) + scale_fill_manual(name = "Taxonomic_review", values = snake.col) + labs(y = NULL, x = NULL) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=8, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=6, colour = "black"), axis.line.y = element_blank(), axis.ticks.y.left = element_blank(), axis.ticks.y=element_blank(), axis.text.y = element_blank(), strip.text.y = element_blank(), axis.text.x = element_text(angle = 45, hjust=1, vjust=1), legend.position ="none") fig.snakesTR # Arrange the multiple ggplots of different sizes in a single multipanel plot require(cowplot) multiplot <- align_plots(fig.lizardsTR, fig.snakesTR, fig.lizards, fig.snakes, align = "hv") (p.final <- grid.arrange(grobs = multiplot, layout_matrix = rbind(c(1,1,2), # 1st plot covers first 2 columns (1 means plot 1; 2 = plot 2, ...) c(3,3,4)))) myfile.pdf <- paste("WRITE HERE YOUR WORKING DIRECTORY/Fig. 3 - family guidelines.pdf") myfile.png <- paste("WRITE HERE YOUR WORKING DIRECTORY/Fig. 3 - family guidelines.png") ggsave(myfile.pdf, plot=p.final, width=10, height = 5, units="in", dpi = "print") ggsave(myfile.png, plot=p.final, width=10, height = 5, units="in", dpi = "print") ##### # STEP 13 - Generate the Figure S1 presented in the supporting information. ########################################################################################################################## # STEP 13 - Generate the Figure S1 presented in the supporting information. library(ggplot2) # Supplementary Information S1 - Number of authors per species described between 1992 and 2017 rm(list=setdiff(ls(),c("dataset", "trait_data"))); gc() # Load the taxonomist dataset present in the supplementary material taxonomist_data <- read.csv("taxonomist_data.csv", h=T, stringsAsFactors = T, fileEncoding = "Latin1") # Reorder levels in the variable N_authors for plotting taxonomist_data$N_authors <- as.factor(as.character(taxonomist_data$N_authors)) levels(taxonomist_data$N_authors) taxonomist_data$N_authors <- factor(taxonomist_data$N_authors, levels =c("1","2","3","4","5","6","7","8","9","10","11","13","14","15","16","17","19")) (Fig.S1 <- ggplot(taxonomist_data, aes(N_authors, Description_year))+ geom_boxplot(outlier.size=0.5, na.rm = T, colour="black", size = .5)+ labs(x = "Number of authors/species", y = "Year")+ scale_y_continuous(breaks = c(1992, 1997, 2002, 2007, 2012, 2017))+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt"), axis.title = element_text(size=9, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=10))) # Save the ggplot with the "ggsave" function myfile <- paste(getwd(), "/Authors_per_year.pdf") myfile.png <- paste(getwd(), "/Authors_per_year.png") ggsave(myfile, plot=Fig.S1, width=4, height = 3, units="in", dpi = "print", bg = "transparent") ggsave(myfile.png, plot=Fig.S1, width=4, height = 3, units="in", dpi = "print", bg = "transparent") ##### # STEP 14 - Generate the Figure S2 presented in the supporting information. ########################################################################################################################## # STEP 14 - Generate the Figure S2 presented in the supporting information. # Sensitivity Analysis for the variable taxonomic review in snakes (single predictor models) rm(list=setdiff(ls(),c("trait_data"))); gc() set.seed(123) # Install and load R packages needed to run the analysis: needed_packages<-c("survival","flexsurv","MuMIn","usdm","plyr","data.table","dplyr","knitr") new.packages<-needed_packages[!(needed_packages %in% installed.packages()[,"Package"])] if(length(new.packages)) install.packages(new.packages) lapply(needed_packages, require, character.only = TRUE) # Subset the continuous predictors that will be tested: names(trait_data) continuous_predictors<-trait_data[, c(6:11)] summary(continuous_predictors) # Standardised the predictors to make them comparable (median = 0, SD = 1): continuous_predictors<-scale(continuous_predictors, center=T, scale=T) continuous_predictors<-as.data.frame(continuous_predictors) # Subset the categorical predictors that will be tested, and add the group and biogeographic variables for filtering: names(trait_data) categorical_predictors <- trait_data[,c(12:14)] summary(categorical_predictors) # Merge the predictors and response variables in a single dataset: dataset <- cbind(trait_data[,c(1, 3:5, 15)], continuous_predictors, categorical_predictors) dataset$TimeLag <- dataset$TimeLag+1 # add +1 to the time lag to avoid erros when running the AFT model dataset <- dataset[, c(1:4,6:14,5)] ## Remove the rows (species) with NA values: dataset <- dataset[complete.cases(dataset[,]),] dataset$Species <- droplevels(dataset$Species) # Clean up the Global Environment rm(categorical_predictors,continuous_predictors, new.packages, needed_packages) # Create a new dataset selecting only snakes levels(dataset$Group) df_snakes <- subset(dataset, Group == "Snakes", select = Species:Molecular) # Create an object to 'guide' the for loop: collection_dates<-seq(from=1952, to=1992, by=5) # Set a full model formula full_model_formula<-as.formula(Surv(TimeLag, Censor)~Taxonomic_review) # Create an empty list to store the outputs: my_outputs<-list() # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create an subset based on species collected in a given time period, with equal nº species described on TR and non-TR: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] l <- sum(data_subset$Taxonomic_review==1) out <- lapply(1:100, function(i){ # Create a list with 100 dataframes for each time period data_subset %>% group_by(Taxonomic_review) %>% sample_n(size = l) }) for (i in 1:100){ # each one of the iterations mydata<-as.data.frame(out[[i]]) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the time-period sensitivity_output[i,3] <- fullmodel[["coefficients"]][2] # average coefficient (Estimates) sensitivity_output[i,4] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,1] # lower IC sensitivity_output[i,5] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,2] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_50.50_TRonly.RData") # Re-run the same analysis, but using the 60:40 proportion of species described on TR vs Non-TR. rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes", "Sensitivity_plot", "collection_dates", "full_model_formula"))); gc(); set.seed(123) # Create an empty list to store the outputs: my_outputs<-list() # Object to store subsets of dataframes out <- list() # Function to get subsets based on different proportions of TR:non-TR get_subsets <- function(df, variable){ df1 <- df[sample(which(variable=='0'),round(1.5*length(which(variable=='1')))), ] # 60:40 TR:non-TR df2 <- df[sample(which(variable=='1'),length(which(variable=='1'))), ] df_new <- rbind(df1, df2) } # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create a subset of 100 dataframes selecting 60:40 of species described on TR and non-TR, for a given time period: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] for (i in 1:100){ # each one of the iterations mydata <- get_subsets(data_subset, data_subset$Taxonomic_review) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the time-period sensitivity_output[i,3] <- fullmodel[["coefficients"]][2] # average coefficient (Estimates) sensitivity_output[i,4] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,1] # lower IC sensitivity_output[i,5] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,2] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_60.40_TRonly.RData") # Re-run the same analysis, but using the 70:30 proportion of species described on TR vs Non-TR. rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes", "Sensitivity_plot", "collection_dates", "full_model_formula"))); gc(); set.seed(123) # Create an empty list to store the outputs: my_outputs<-list() # Object to store subsets of dataframes out <- list() # Function to get subsets based on different proportions of TR:non-TR get_subsets <- function(df, variable){ df1 <- df[sample(which(variable=='0'),round(2.33*length(which(variable=='1')))), ] # 70:30 TR:non-TR df2 <- df[sample(which(variable=='1'),length(which(variable=='1'))), ] df_new <- rbind(df1, df2) } # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create a subset of 100 dataframes selecting 60:40 of species described on TR and non-TR, for a given time period: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] for (i in 1:100){ # each one of the iterations mydata <- get_subsets(data_subset, data_subset$Taxonomic_review) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the time-period sensitivity_output[i,3] <- fullmodel[["coefficients"]][2] # average coefficient (Estimates) sensitivity_output[i,4] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,1] # lower IC sensitivity_output[i,5] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,2] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_70.30_TRonly.RData") # Re-run the same analysis, but using the 80:20 proportion of species described on TR vs Non-TR. rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes", "Sensitivity_plot", "collection_dates", "full_model_formula"))); gc(); set.seed(123) # Create an empty list to store the outputs: my_outputs<-list() # Object to store subsets of dataframes out <- list() # Function to get subsets based on different proportions of TR:non-TR get_subsets <- function(df, variable){ df1 <- df[sample(which(variable=='0'),round(4*length(which(variable=='1')))), ] # 80:20 TR:non-TR df2 <- df[sample(which(variable=='1'),length(which(variable=='1'))), ] df_new <- rbind(df1, df2) } # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create a subset of 100 dataframes selecting 60:40 of species described on TR and non-TR, for a given time period: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] for (i in 1:100){ # each one of the iterations mydata <- get_subsets(data_subset, data_subset$Taxonomic_review) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the time-period sensitivity_output[i,3] <- fullmodel[["coefficients"]][2] # average coefficient (Estimates) sensitivity_output[i,4] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,1] # lower IC sensitivity_output[i,5] <- confint(fullmodel, "Taxonomic_review", full = T, level = .95)[1,2] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_80.20_TRonly.RData") # Clean the workspace rm(list=setdiff(ls(),c("trait_data", "df_snakes", "Sensitivity_plot"))); gc() # Create the subplot A from Figure S2. # Load outputs (and rename objects) from previous analysis load("my_outputs_50.50_TRonly.RData"); my_outputs50.50 <- my_outputs load("my_outputs_60.40_TRonly.RData"); my_outputs60.40 <- my_outputs load("my_outputs_70.30_TRonly.RData"); my_outputs70.30 <- my_outputs load("my_outputs_80.20_TRonly.RData"); my_outputs80.20 <- my_outputs; rm(my_outputs) # Merge all dataframes into a single big dataframe outputs50.50 <- bind_rows(my_outputs50.50) outputs60.40 <- bind_rows(my_outputs60.40) outputs70.30 <- bind_rows(my_outputs70.30) outputs80.20 <- bind_rows(my_outputs80.20) # Get mean AvgCoef by TimePeriod outputs50.50$TimePeriod <- as.factor(outputs50.50$TimePeriod) outputs60.40$TimePeriod <- as.factor(outputs60.40$TimePeriod) outputs70.30$TimePeriod <- as.factor(outputs70.30$TimePeriod) outputs80.20$TimePeriod <- as.factor(outputs80.20$TimePeriod) outputs50.50 <- outputs50.50 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs50.50$prop <- '50:50' outputs60.40 <- outputs60.40 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs60.40$prop <- '60:40' outputs70.30 <- outputs70.30 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs70.30$prop <- '70:30' outputs80.20 <- outputs80.20 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs80.20$prop <- '80:20' # Combine all dataframes outputs <- rbind(outputs50.50, outputs60.40, outputs70.30, outputs80.20) outputs$prop <- as.factor(outputs$prop) outputs$TimePeriod <- factor(outputs$TimePeriod, labels = c('1952-2017','1957-2017','1962-2017','1967-2017','1972-2017', '1977-2017','1982-2017','1987-2017','1992-2017')) # Prepare for plotting library(ggplot2) library(viridis) # Define colors to be used in the plot myColors <- viridis_pal(option="plasma")(4) names(myColors)<-levels(outputs$prop) (Sensitivity_plotA <- ggplot(outputs, aes(x = TimePeriod, y = AvgCoef.mean, ymin = lower95.mean, ymax = upper95.mean, group=prop))+ geom_pointrange(aes(col = prop, shape = prop), size = 0.4, position=position_dodge(width=2))+ scale_colour_manual(name = "prop", values = myColors)+ scale_shape_manual(values=c(0,1,2,4))+ geom_errorbar(aes(ymin=lower95.mean, ymax=upper95.mean, col = prop), width=0.1, position=position_dodge(width=2))+ geom_hline(yintercept =0, linetype=2)+ labs(x='Time Period', y='Mean Avg. Coefs. (IC95)')+ ylim(-.1, 1.2) + facet_wrap(~TimePeriod, strip.position="left", nrow=9, scales = "free_y")+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=12, face="bold"), axis.title.x = element_blank(), axis.text = element_text(size=10), axis.line = element_line(colour="black"), axis.ticks.y.left = element_blank(), axis.text.x=element_text(face="italic"), axis.text.y = element_blank(), plot.background=element_rect(fill = "white"), strip.background = element_blank(), strip.placement = "outside", strip.text.y = element_text(hjust=1, vjust=0.5, angle=180, size = 10), legend.key = element_blank(), legend.title = element_blank(), legend.position = "none")+ coord_flip()) # Same as above, but with taxonomic review included in multi predictor models. rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes", "Sensitivity_plotA", "collection_dates"))); gc() set.seed(123) # Set a full model formula full_model_formula<-as.formula(Surv(TimeLag, Censor)~Log10_mass + Lat + Lat_mz + LogN_specimens_TS + LogN_spp_genus + LogN_authors + Col_hol_is_author + Taxonomic_review + Molecular) # Create an empty list to store the outputs: my_outputs<-list() # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create a subset based on species collected in a given time period, with equal nº of species described based on TR and non-TR: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] l <- sum(data_subset$Taxonomic_review==1) out <- lapply(1:100, function(i){ # Create a list with 100 dataframes for each time period data_subset %>% group_by(Taxonomic_review) %>% sample_n(size = l) }) for (i in 1:100){ # each one of the iterations mydata<-as.data.frame(out[[i]]) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) fits <- dredge(fullmodel) # Performs automated model selection with subsets of the supplied global model. avgmodel_output <- model.avg(fits, revised.var = T) # Get the model averaged coefficients ma_coefs <- coefTable(avgmodel_output, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(Predictors, Estimate, lower95, upper95) row.names(ma_coefs) <- ma_coefs[,1] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","Log10_mass","Lat","Lat_mz","LogN_specimens_TS","LogN_spp_genus", "LogN_authors","Col_hol_is_author","Taxonomic_review","Molecular"),] sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the timeperiod sensitivity_output[i,3] <- ma_coefs$Estimate[[10]] # average coefficient (Estimates) sensitivity_output[i,4] <- ma_coefs$lower95[[10]] # lower IC sensitivity_output[i,5] <- ma_coefs$upper95[[10]] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_50.50.RData") # Re-run the same analysis, but using the 60:40 proportion of species described on TR vs Non-TR. rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes", "collection_dates", "full_model_formula"))); gc() set.seed(123) # Create an empty list to store the outputs: my_outputs<-list() # Object to store subsets of dataframes out <- list() # Function to get subsets based on different proportions of TR:non-TR get_subsets <- function(df, variable){ df1 <- df[sample(which(variable=='0'),round(1.5*length(which(variable=='1')))), ] # 60:40 TR:non-TR df2 <- df[sample(which(variable=='1'),length(which(variable=='1'))), ] df_new <- rbind(df1, df2) } # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create a subset of 100 dataframes selecting 60:40 of species described based on TR and non-TR, for a given time period: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] for (i in 1:100){ # each one of the iterations mydata <- get_subsets(data_subset, data_subset$Taxonomic_review) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) fits <- dredge(fullmodel) # Performs automated model selection with subsets of the supplied global model avgmodel_output <- model.avg(fits, revised.var = T) # Get the model averaged coefficients ma_coefs <- coefTable(avgmodel_output, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(Predictors, Estimate, lower95, upper95) row.names(ma_coefs) <- ma_coefs[,1] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","Log10_mass","Lat","Lat_mz","LogN_specimens_TS", "LogN_spp_genus","LogN_authors","Col_hol_is_author","Taxonomic_review","Molecular"),] sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the timeperiod sensitivity_output[i,3] <- ma_coefs$Estimate[[10]] # average coefficient (Estimates) sensitivity_output[i,4] <- ma_coefs$lower95[[10]] # lower IC sensitivity_output[i,5] <- ma_coefs$upper95[[10]] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_60.40.RData") # Re-run the same analysis, but using the 70:30 proportion of species described on TR vs Non-TR. rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes", "collection_dates", "full_model_formula"))); gc() set.seed(123) # Create an empty list to store the outputs: my_outputs<-list() # Object to store subsets of dataframes out <- list() # Function to get subsets based on different proportions of TR:non-TR get_subsets <- function(df, variable){ df1 <- df[sample(which(variable=='0'),round(2.33*length(which(variable=='1')))), ] # 70:30 TR:non-TR df2 <- df[sample(which(variable=='1'),length(which(variable=='1'))), ] df_new <- rbind(df1, df2) } # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create a subset of 100 dataframes selecting 60:40 of species described based on TR and non-TR, for a given time period: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] for (i in 1:100){ # each one of the iterations mydata <- get_subsets(data_subset, data_subset$Taxonomic_review) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) fits <- dredge(fullmodel) # Performs automated model selection with subsets of the supplied global model avgmodel_output <- model.avg(fits, revised.var = T) # Get the model averaged coefficients ma_coefs <- coefTable(avgmodel_output, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(Predictors, Estimate, lower95, upper95) row.names(ma_coefs) <- ma_coefs[,1] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","Log10_mass","Lat","Lat_mz","LogN_specimens_TS", "LogN_spp_genus","LogN_authors","Col_hol_is_author","Taxonomic_review","Molecular"),] sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the timeperiod sensitivity_output[i,3] <- ma_coefs$Estimate[[10]] # average coefficient (Estimates) sensitivity_output[i,4] <- ma_coefs$lower95[[10]] # lower IC sensitivity_output[i,5] <- ma_coefs$upper95[[10]] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_70.30.RData") # Re-run the same analysis, but using the 80:20 proportion of species described on TR vs Non-TR. rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes", "collection_dates", "full_model_formula"))); gc() set.seed(123) # Create an empty list to store the outputs: my_outputs<-list() # Object to store subsets of dataframes out <- list() # Function to get subsets based on different proportions of TR:non-TR get_subsets <- function(df, variable){ df1 <- df[sample(which(variable=='0'),round(4*length(which(variable=='1')))), ] # 80:20 TR:non-TR df2 <- df[sample(which(variable=='1'),length(which(variable=='1'))), ] df_new <- rbind(df1, df2) } # Create the loop: for (j in 1:length(collection_dates)){ # j = each one of the time periods # Create an empty dataframe sensitivity_output<-as.data.frame(matrix(nrow=100, ncol=5)) names(sensitivity_output)<-c("Iter", "TimePeriod", "AvgCoef", "LowerCoef", "UpperCoef") # Create a subset of 100 dataframes selecting 60:40 of species described based on TR and non-TR, for a given time period: data_subset<-df_snakes[which(df_snakes$Year_hol>=collection_dates[j]),] for (i in 1:100){ # each one of the iterations mydata <- get_subsets(data_subset, data_subset$Taxonomic_review) fullmodel <- survreg(full_model_formula, mydata, dist="lognorm", na.action = na.fail) fits <- dredge(fullmodel) # Performs automated model selection with subsets of the supplied global model avgmodel_output <- model.avg(fits, revised.var = T) # Get the model averaged coefficients ma_coefs <- coefTable(avgmodel_output, full=TRUE) coefnames <- row.names(ma_coefs) ma_coefs <- as.data.frame(ma_coefs) ma_coefs <- mutate(ma_coefs, Predictors = coefnames, lower95 = Estimate - 1.96 * `Std. Error`, upper95 = Estimate + 1.96 * `Std. Error`) %>% dplyr::select(Predictors, Estimate, lower95, upper95) row.names(ma_coefs) <- ma_coefs[,1] ma_coefs <- ma_coefs[c("(Intercept)","Log(scale)","Log10_mass","Lat","Lat_mz","LogN_specimens_TS", "LogN_spp_genus","LogN_authors","Col_hol_is_author","Taxonomic_review","Molecular"),] sensitivity_output[i,1] <- paste(i) # store the number of the iteration sensitivity_output[i,2] <- paste(j) # store the timeperiod sensitivity_output[i,3] <- ma_coefs$Estimate[[10]] # average coefficient (Estimates) sensitivity_output[i,4] <- ma_coefs$lower95[[10]] # lower IC sensitivity_output[i,5] <- ma_coefs$upper95[[10]] # upper IC print(i) } # end of the for loop i # store the output of a given time period my_outputs[[j]]<-sensitivity_output } # end of the for loop j save(my_outputs, file = "my_outputs_80.20.RData") # Create the subplot B from Figure S2. # Load outputs (and rename objects) from previous analysis load("my_outputs_50.50.RData"); my_outputs50.50 <- my_outputs load("my_outputs_60.40.RData"); my_outputs60.40 <- my_outputs load("my_outputs_70.30.RData"); my_outputs70.30 <- my_outputs load("my_outputs_80.20.RData"); my_outputs80.20 <- my_outputs; rm(my_outputs) # Merge all dataframes into a single big dataframe outputs50.50 <- bind_rows(my_outputs50.50) outputs60.40 <- bind_rows(my_outputs60.40) outputs70.30 <- bind_rows(my_outputs70.30) outputs80.20 <- bind_rows(my_outputs80.20) # Get mean AvgCoef by TimePeriod outputs50.50$TimePeriod <- as.factor(outputs50.50$TimePeriod) outputs60.40$TimePeriod <- as.factor(outputs60.40$TimePeriod) outputs70.30$TimePeriod <- as.factor(outputs70.30$TimePeriod) outputs80.20$TimePeriod <- as.factor(outputs80.20$TimePeriod) outputs50.50 <- outputs50.50 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs50.50$prop <- '50:50' outputs60.40 <- outputs60.40 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs60.40$prop <- '60:40' outputs70.30 <- outputs70.30 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs70.30$prop <- '70:30' outputs80.20 <- outputs80.20 %>% group_by(TimePeriod) %>% summarise(AvgCoef.mean = mean(AvgCoef), lower95.mean = mean(LowerCoef), upper95.mean = mean(UpperCoef)); outputs80.20$prop <- '80:20' # Combine all dataframes outputs <- rbind(outputs50.50, outputs60.40, outputs70.30, outputs80.20) outputs$prop <- as.factor(outputs$prop) # Prepare for plotting library(ggplot2) library(viridis) myColors <- viridis_pal(option="plasma")(4) names(myColors)<-levels(outputs$prop) (Sensitivity_plotB <- ggplot(outputs, aes(x = TimePeriod, y = AvgCoef.mean, ymin = lower95.mean, ymax = upper95.mean, group=prop))+ geom_pointrange(aes(col = prop, shape = prop), size = 0.4, position=position_dodge(width=2))+ scale_colour_manual(name = "prop", values = myColors)+ scale_shape_manual(values=c(0,1,2,4))+ geom_errorbar(aes(ymin=lower95.mean, ymax=upper95.mean, col = prop), width=0.1, position=position_dodge(width=2))+ geom_hline(yintercept =0, linetype=2)+ labs(x=NULL, y=NULL)+ ylim(-.15, 1) + facet_wrap(~TimePeriod, strip.position="left", nrow=9, scales = "free_y")+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(size=12, face="bold"), axis.text = element_text(size=10), axis.line = element_line(colour="black"), axis.line.y = element_blank(), axis.ticks.y.left = element_blank(), axis.ticks.y=element_blank(), axis.text.x=element_text(face="italic"), axis.text.y = element_blank(), plot.background=element_rect(fill = "white"), strip.background = element_blank(), strip.placement = "outside", strip.text.y = element_blank(), legend.key = element_blank(), legend.title = element_blank(), legend.position = c(0.9,0.2))+ coord_flip()) ##### # STEP 15 - Create the multi panel plot Figure S2. ########################################################################################################################## # STEP 15 - Create the multi panel plot Figure S2. library("ggpubr") (Sensitivity_plot_comb <- ggarrange(Sensitivity_plotA,Sensitivity_plotB, ncol = 2, nrow = 1, labels = c('A', 'B'), label.x = c(0.25, 0.1), label.y = c(.98,.98), widths = c(1,1), font.label = list(size = 10, color = "black"), align = "h")) (Sensitivity_plot_comb <- annotate_figure(Sensitivity_plot_comb, bottom = text_grob("Mean Avg. Coefs. (IC95)", color = "black", face = "bold", size = 12, vjust = 0, hjust = 0))) # Save the plot myfile <- paste(getwd(), "/Sensitivity plot.pdf") myfile.png <- paste(getwd(), "/Sensitivity plot.png") ggsave(myfile, plot=Sensitivity_plot_final, width=10, height = 7, units="in", dpi = "print") ggsave(myfile.png, plot=Sensitivity_plot_final, width=10, height = 7, units="in", dpi = "print") # This figure was later edited in the program Inkscape for minor aesthetics adjustments, # and for adding the sample size (obtained below) to the Figure S2 # Get the number of species (sample size) per time period and proportion rm(list=setdiff(ls(),c("dataset", "trait_data", "df_snakes"))); gc() set.seed(123) # Time period to 'guide' the for loop below collection_dates<-seq(from=1952, to=1992, by=5) # Create an empty datafrase to store the outputs N_spp<-as.data.frame(matrix(nrow=9, ncol=6)) names(N_spp)<-c("TimePeriod","n",'50.50','60.40','70.30','80.20') # Get the sample size for (i in 1:length(collection_dates)){ # i = each one of the time periods data_subset <- df_snakes[which(df_snakes$Year_hol>=collection_dates[i]), ] N_spp[i,1] <- paste(i) # TimePeriod N_spp[i,2] <- summarise(data_subset, count = n()) # Sample size N_spp[i,3] <- sum(data_subset$Taxonomic_review==1)*2 # proportion of 50.50 for TR:non-TR, respectively N_spp[i,4] <- round(sum(data_subset$Taxonomic_review==1)*1.5)+sum(data_subset$Taxonomic_review==1) # proportion of 60.40 for TR:non-TR, respectively N_spp[i,5] <- round(sum(data_subset$Taxonomic_review==1)*2.33)+sum(data_subset$Taxonomic_review==1) # proportion of 70.30 for TR:non-TR, respectively N_spp[i,6] <- (sum(data_subset$Taxonomic_review==1)*4)+sum(data_subset$Taxonomic_review==1) # proportion of 80.20 for TR:non-TR, respectively } # Save the dataframe write.csv(N_spp, "N_spp.csv", row.names = F) ##### # STEP 16 - Generate the Figure S3 presented in the supporting information. ########################################################################################################################## # STEP 16 - Generate the Figure S3 presented in the supporting information. # Variation in collection dates by size of type series rm(list=setdiff(ls(),c("dataset", "trait_data"))); gc() # Load the lizard collection dates dataset lizard_collection_dates <- read.csv("lizard_type_series.csv", h=T, stringsAsFactors = T) # Remove unique rows (species with 50+ type-specimens - not included in the compilation) lizard_collection_dates <- subset(lizard_collection_dates, duplicated(Species) | duplicated(Species, fromLast=TRUE)) lizard_collection_dates <- droplevels(lizard_collection_dates) # Prepare for plotting library(dplyr) library(raster) library(stats) collection_dates <- lizard_collection_dates %>% group_by(Species) %>% summarise(range = max(Collection_year, na.rm = TRUE) - min(Collection_year, na.rm = TRUE), # get the range in collection dates for each species mad = mad(Collection_year, na.rm = T), # get the median absolute deviation n = sum(!is.na(Collection_year))) # Select only species with 2 or more collection dates available in the type series (for plotting) coll_dates_2orMore <- subset(collection_dates, n >= 2) # for the range plot coll_dates_3orMore <- subset(collection_dates, n >= 3) # for the MAD plot # Install needed packages to create the density plot in Fig. S2 install.packages("remotes") remotes::install_github("caleblareau/BuenColors") install.packages("remotes") remotes::install_github("sckott/rphylopic") library(rphylopic) library(BuenColors) library(tidyverse) library(data.table) library(cowplot) library(ggpubr) library(RColorBrewer) library(ggforce) library(ggplot2) library(viridis) library(scales) coll_dates_2orMore$density<-get_density(x=coll_dates_2orMore$n, y=coll_dates_2orMore$range, n=500) # Fig. A (fig.A <- ggplot(coll_dates_2orMore) + geom_point(aes(x=n, y=range, color = density), size =.9) + scale_color_viridis(option="cividis", direction=1) + stat_smooth(method="loess", aes(x=n, y=range), size=.8, col="red") + labs(x=NULL, y="Range of collection dates") + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(face = "bold", size = 9), axis.text = element_text(size = 7), axis.line = element_line(colour="black"), plot.margin = rep(unit(0,"null"),4), panel.spacing = unit(0,"null"), legend.position="none", plot.background=element_rect(fill="transparent", colour=NA)) + scale_y_continuous(limits = c(0, 170))) coll_dates_3orMore$density<-get_density(x=coll_dates_3orMore$n, y=coll_dates_3orMore$range, n=500) # Fig. B (fig.B <- ggplot(coll_dates_3orMore) + geom_point(aes(x=n, y=mad, color = density), size = .9) + scale_color_viridis(option="cividis", direction=1) + stat_smooth(method="loess", aes(x=n, y=mad), size=.8, col="red") + labs(x="N of type-specimens", y="M.A.D. of collection dates") + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), axis.title = element_text(face = "bold", size = 9), axis.text = element_text(size = 7), axis.line = element_line(colour="black"), plot.margin = rep(unit(0,"null"),4), panel.spacing = unit(0,"null"), legend.position="none", plot.background=element_rect(fill="transparent", colour=NA)) + scale_y_continuous(limits = c(0, 60))) ## Create a Multipanel plot (fig <- ggarrange(fig.A, fig.B, ncol = 1, nrow = 2, labels = c('A', 'B'), label.x = c(0.13, 0.13), label.y = c(1,1), font.label = list(size = 8, color = "black"), align = "hv")) # Save the ggplot myfile <- paste(getwd(), "/Collection_dates_variation.pdf") myfile.png <- paste(getwd(), "/Collection_dates_variation.png") ggsave(myfile, plot=fig, width=4, height = 5, units="in", dpi = "print") ggsave(myfile.png, plot=fig, width=4, height = 5, units="in", dpi = "print") ##### # STEP 17 - Generate the Figure S4 presented in the supporting information. ########################################################################################################################## # STEP 17 - Generate the Figure S4 presented in the supporting information. # Correlation plot rm(list = ls()); gc() library(ggplot2) library(plyr) # Load trait data trait_data <- read.csv("trait_data.csv", h=T, stringsAsFactors=T) # Fig S4.A - Number of species described per year from 1992 to 2017 (Fig.A <- ggplot(trait_data, aes(Year_of_description))+ geom_histogram(binwidth = 1, na.rm = T, colour="black", fill="grey50")+ labs(x =NULL, y = "Species described/year")+ scale_x_continuous(breaks = c(1992, 1997, 2002, 2007, 2012, 2017))+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt"), axis.title = element_text(size=13, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=10))) # Fig S4.B - Frequency plot with the number of holotypes by a given time lag value # Histogram of observed time lag (Fig.B <- ggplot(trait_data, aes(TimeLag))+ geom_histogram(binwidth = 1, na.rm = T, colour="black", fill="grey50")+ labs(x = NULL, y = "Number of holotypes")+ xlim(0, 40)+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt"), axis.title = element_text(size=14, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=10, colour = "black"))) # Create another histogram with the full range to plot inside the previous one (Fig.B_inset <- ggplot(trait_data, aes(TimeLag))+ geom_histogram(binwidth = 1, na.rm = T, colour="black", fill="grey50")+ labs(x = NULL, y = NULL)+ scale_x_continuous(breaks = c(0 , 75, 150))+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt"), axis.title = element_text(size=13, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=6, colour = "black"))) # Place the last plot within Fig.B (Fig.B_final <- Fig.B + annotation_custom(ggplotGrob(Fig.B_inset), xmin = 25, xmax = 40, ymin = 200, ymax = 330)); rm(Fig.B, Fig.B_inset) ## Fig S4.C - frequency plot with the year of collection of holotypes (Fig.C <- ggplot(trait_data, aes(Year_hol))+ geom_histogram(binwidth = 1, na.rm = T, colour="black", fill="grey50")+ labs(x ="Year", y = "Holotypes collected/year")+ scale_x_continuous(breaks = seq(1952, 2017, 10), lim = c(1952, 2018))+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt"), axis.title = element_text(size=13, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=10))) # Create an inset histogram with the full range (Fig.C_inset <- ggplot(trait_data, aes(Year_hol))+ geom_histogram(binwidth = 1, na.rm = T, colour="black", fill="grey50")+ labs(x =NULL, y = NULL)+ scale_x_continuous(breaks = c(1852, 1934, 2017))+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt"), axis.title = element_text(size=13, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=6, colour = "black"))) # Place the last plot within Fig.C (Fig.C_final <- Fig.C + annotation_custom(ggplotGrob(Fig.C_inset), xmin = 1952, xmax = 1974, ymin = 72, ymax = 110)); rm(Fig.C, Fig.C_inset) # Fig S4.D - Get the pearson correlation between the number of holotypes collected/year and the number of descriptions/year by time lag values # Load correlation data df <- read.csv("correlation_data.csv", h=T) df$Timelag <- as.factor(df$Timelag) # Create a correlation function corr <- function(x, y) { corr=(cor.test(x, y, method="pearson")) } # Get the correlation for each time lag value correlation <- ddply(df, .(Timelag), summarise, COR=corr(n_description,n_holotype)$estimate, lower=corr(n_description,n_holotype)$conf.int[[1]], upper=corr(n_description,n_holotype)$conf.int[[2]]) df <- join(df, correlation, by = "Timelag"); rm(correlation) (Fig.D <- ggplot(df, aes(as.numeric(as.character(Timelag)), COR))+ geom_smooth(method="loess", se=F, colour="black", fill="grey50", span = .25)+ geom_ribbon(aes(ymin=lower, ymax=upper), alpha=.1, colour = "gray50") + labs(x = "Time lag (years)", y= "Pearson correlation")+ theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background = element_blank(), plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt"), axis.title = element_text(size=13, face="bold"), axis.line = element_line(colour="black"), axis.text = element_text(size=10, colour = "black"))) # Create a multi panel plot, which is the Fig S1 from the paper library(ggpubr) (Fig.S4 <- ggarrange(Fig.A, Fig.B_final, Fig.C_final, Fig.D, ncol = 2, nrow = 2, labels = c('A', 'B', 'C', 'D'), label.x = c(0, 0, 0, 0), label.y = c(1,1,1,1), font.label = list(size = 12, color = "black"), align = "hv")) # Save the ggplot myfile <- paste(getwd(), "/correlation_plot.pdf") myfile.png <- paste(getwd(), "/correlation_plot.png") ggsave(myfile, plot=Fig.S4, width=7, height = 6, units="in", dpi = "print") ggsave(myfile.png, plot=Fig.S4, width=7, height = 6, units="in", dpi = "print") rm(list=setdiff(ls(),c("trait_data"))); gc() #####