# reading in packages required #### line 34-78 # functions #### line 82-1430 # making master occurrence files- pres = present holo = holocene #### acer pres master line 1433 #### acer holo master line 1631 #### apal pres master line 1631 #### apal holo master line 1858 #### cnat pres master line 1861 #### cnat holo master line 2057 #### past pres master line 2061 #### past holo master line 2280 # ecospat line #-# #### acer ecospat line 2293 #### apal ecospat line 2550 #### cnat ecospat line 2805 #### past ecospat line 3079 #ecological niche models pres = present holo = holocene #### acer pres ENM line 3390 #### acer holo ENM line 3613 #### apal pres ENM line 3868 #### apal holo ENM line 4105 #### cnat pres ENM line 4338 #### cnat holo ENM line 4593 #### past pres ENM line 4849 #### past holo ENM line 5113 ############ plotting by area and latitude ### latitude line 5365 #################&&&&&&&&&&&&&&&&&&&&&##############Reading in packages # loading the requisite libraries for generating the master.background and master.occurrence files library(raster) library(dplyr) #ecospat packages library(ENMTools) library(dismo) library(raster) library(rgdal) # for spatial data analysis library(dplyr) library(rgeos) # for spatial data analysis library(scales) library(tidyr) library(colorRamps) library(ENMeval) # for a few new tools in ENM/SDM library(ggplot2) library(maptools) library(spThin) # error here library(sdm) library(ecospat) library(ade4) library(dichromat) library(terra) # ENM model packages packs <- c('ade4', 'adehabitatMA', 'adehabitatHR', 'alphahull', 'dismo', 'dplyr', 'ecospat', 'ggplot2', 'jsonlite', 'kuenm', 'matrixStats', 'raster', 'rgdal', 'rgeos', 'rJava', 'sf', 'sp', 'splitstackshape', 'tidyr', 'utils', 'wesanderson', 'PerformanceAnalytics', 'SDMTools') # loading in the packages lapply(packs, library, character.only=T) # giving the package versions packs.df <- as.data.frame(matrix(NA, nrow=length(packs), ncol=2)) colnames(packs.df) <- c('pkg.name', 'pkg.version') for(i in 1:length(packs)){ packs.df[i,1] <- packs[i] packs.df[i,2] <- as.character(packageVersion(packs[i])) } packs.df library(raster) library(rasterize) ########################### FUNCTIONS ######### ##### Prep Parameters ------------------------------------------------------------------------------------------------------ # Prep Parameters for Maxent models in R with the dismo package # sourced from https://github.com/shandongfx/workshop_maxent_R/blob/master/code/Appendix2_prepPara.R # https://github.com/shandongfx/workshop_maxent_R/blob/master/code # should be used in concert with Appendix 3 from Feng et al. 2017 (PeerJ Preprints) # https://peerj.com/preprints/3346.pdf (manuscript still unpublished as of August 2020) # Appendix 3 # https://github.com/shandongfx/workshop_maxent_R/blob/master/code/Appendix3_maxentParameters_v2.pdf # some arguments may change, or may/may not be used depending on if you're using raster data vs "samples-with-data" (SWD; column data) # A function that implements Maxent parameters using the general R manner # leave "doclamp" as default - later in the code (internal fxns run), "doclamp" is set to FALSE prepPara <- function(userfeatures=NULL, # 41 NULL=autofeature, could be any combination of # c("L", "Q", "H", "P", "T") # MUST be specified as a single string (e.g., "LQ", "LQP", "LQHPT", etc.) responsecurves=TRUE, # 1 jackknife=TRUE, # 3 outputformat="logistic", # 4 outputfiletype="asc", # 5 projectionlayers=NULL, # 7 randomseed=FALSE, # 10 removeduplicates=TRUE, # 16 betamultiplier=NULL, # 20, 53-56 biasfile=NULL, # 22 testsamplesfile=NULL, # 23 replicates=1, # 24-25 replicatetype="crossvalidate", # 24-25 writeplotdata=TRUE, # 37 extrapolate=TRUE, # 39 doclamp=TRUE, # 42 beta_threshold=NULL, # 20, 53-56 beta_categorical=NULL, # 20, 53-56 beta_lqp=NULL, # 20, 53-56 beta_hinge=NULL, # 20, 53-56 applythresholdrule=NULL # 60 ){ #20, 29-33, & 41 features, default is autofeature if(is.null(userfeatures)){ args_out <- c("autofeature") } else { args_out <- c("noautofeature") if(grepl("L",userfeatures)) args_out <- c(args_out,"linear") else args_out <- c(args_out,"nolinear") if(grepl("Q",userfeatures)) args_out <- c(args_out,"quadratic") else args_out <- c(args_out,"noquadratic") if(grepl("H",userfeatures)) args_out <- c(args_out,"hinge") else args_out <- c(args_out,"nohinge") if(grepl("P",userfeatures)) args_out <- c(args_out,"product") else args_out <- c(args_out,"noproduct") if(grepl("T",userfeatures)) args_out <- c(args_out,"threshold") else args_out <- c(args_out,"nothreshold") } # 1 - generate response curves for each variable if(responsecurves) args_out <- c(args_out,"responsecurves") else args_out <- c(args_out,"noresponsecurves") # 2 #if(picture) args_out <- c(args_out,"pictures") else args_out <- c(args_out,"nopictures") # 3 - apply variable jackknife to see how the model changes if that variable is omitted, then if it's the ONLY variable used if(jackknife) args_out <- c(args_out,"jackknife") else args_out <- c(args_out,"nojackknife") # 4 - output format type. choose from c("logistic", "cumulative", "raw") args_out <- c(args_out,paste0("outputformat=",outputformat)) # 5 - output file type. choose from c("asc", "mxe", "grd", "bil") args_out <- c(args_out,paste0("outputfiletype=",outputfiletype)) # 7 - pathway to projection layers. # it seems that the projection layers should be the only files in that folder, just like with the MaxEnt .jar file if(!is.null(projectionlayers)) args_out <- c(args_out,paste0("projectionlayers=",projectionlayers)) # 10 - will use different random number generators for selecting training vs testing data and background points (if applicable) if(randomseed) args_out <- c(args_out,"randomseed") else args_out <- c(args_out,"norandomseed") # 16 - remove duplicate coordinates that are in the same grid - ONLY for raster data, not SWD if(removeduplicates) args_out <- c(args_out,"removeduplicates") else args_out <- c(args_out,"noremoveduplicates") # 20 & 53-56 - various beta (regularization) multipliers to be applied. default = 1 # 20 applies all parameters by this regularization multiplier. # 53-56 can apply uniquely to different feature types # check if negative betas <- c( betamultiplier,beta_threshold,beta_categorical,beta_lqp,beta_hinge) if(! is.null(betas) ){ for(i in 1:length(betas)){ if(betas[i] <0) stop("betamultiplier has to be positive") } } if ( !is.null(betamultiplier) ){ args_out <- c(args_out,paste0("betamultiplier=",betamultiplier)) } else { if(!is.null(beta_threshold)) args_out <- c(args_out,paste0("beta_threshold=",beta_threshold)) if(!is.null(beta_categorical)) args_out <- c(args_out,paste0("beta_categorical=",beta_categorical)) if(!is.null(beta_lqp)) args_out <- c(args_out,paste0("beta_lqp=",beta_lqp)) if(!is.null(beta_hinge)) args_out <- c(args_out,paste0("beta_hinge=",beta_hinge)) } # 22 - pathway to a bias file for selecting background points - ONLY for raster data, not SWD if(!is.null(biasfile)) args_out <- c(args_out,paste0("biasfile=",biasfile)) # 23 - pathway to a test data file - can be in csv format if(!is.null(testsamplesfile)) args_out <- c(args_out,paste0("testsamplesfile=",testsamplesfile)) # 24 - replicates = number of replicates to run (integer) # 25 - replicatetype = what type of replicates to run. choose from c('crossvalidate', 'bootstrap', 'subsample') replicates <- as.integer(replicates) if(replicates>1 ){ args_out <- c(args_out, paste0("replicates=",replicates), paste0("replicatetype=",replicatetype) ) } # 37 - write output files containing the data used to make response curves if(writeplotdata) args_out <- c(args_out,"writeplotdata") else args_out <- c(args_out,"nowriteplotdata") # 39 - allow extrapolation beyond the limits of the training data if(extrapolate) args_out <- c(args_out,"extrapolate") else args_out <- c(args_out,"noextrapolate") # 42 - apply clamping when projecting if(doclamp) args_out <- c(args_out,"doclamp") else args_out <- c(args_out,"nodoclamp") # 60 - threshold your model to binary 1/0 # options are: c('Fixed cumulative value 1', 'Fixed cumulative value 5', 'Fixed cumulative value 10', 'Minimum training presence', # '10 percentile training presence', 'Equal training sensitivity and specificity', 'Maximum training sensitivity plus specificity'). if(!is.null(applythresholdrule)) args_out <- c(args_out,paste0("applythresholdrule=",applythresholdrule)) return(args_out) } # prepPara() # [1] "autofeature" "responsecurves" "jackknife" "outputformat=logistic" "outputfiletype=asc" # [6] "norandomseed" "removeduplicates" "writeplotdata" "extrapolate" "doclamp" ##### Create Null Object, Summary Object, Eval Object, and Nested File Structure ------------------------------------------- ### functions that create the null data.frame, summary data.frame, and the eval data.frame ### ARGUMENTS ### ### arguments that exist to keep track of everything and do not change how the functions run ## taxon.name <- name of the entity that will get assigned to any null/summary/eval objects. Does not change how function runs. ## time.bin <- name of the time bin that will get assigned to any null/summary/eval objects. Does not change how function runs. ## extent <- name of the extent that will get assigned to any null/summary/eval objects. Does not change how function runs. ### arguments that set the model hyper-parameters and will change how the functions run ## cv.runs <- name(s) of the cross-validation types. folders will be generated with these names in create.folders.for.maxent. # recommended framework is to treat each cross-validation fold as a letter. # In the case of 5-fold cross validation, for running all the data (using the null.aic and optimize.maxent.likelihood . . . # . . . functions), you should specify 'abcde'. For running cross validation models with maxent.crossval.error, you . . . # . . . should specify c('abcd', 'abce', 'abde', 'acde', 'bcde'). ## f.class <- name(s) of the feature classes used in analysis. ## beta.values <- regularization multipliers used in analysis # create.summary.df and create.eval.df will create every possible combination of cv.runs, f.class, and beta.values # function to create the null data.frame create.null.df <- function(taxon.name, time.bin, extent){ null.object <- as.data.frame(matrix(NA, nrow=1, ncol=17)) colnames(null.object) <- c( "Taxa", "Time.Bin", "Extent", "CrossVal", "cv.num", "Features", "Betas", "n", "k", "ln.L", "AIC", "AICc", "delta.i", "delta.i.c", "w.i", "w.i.c", "lambdas" ) null.object$Taxa <- taxon.name null.object$Time.Bin <- time.bin null.object$Extent <- extent null.object$CrossVal <- 'abcde' null.object$cv.num <- 1 null.object$lambdas <- 'Incercept' return(null.object) } # function to create the summary data.frame create.summary.df <- function(taxon.name, time.bin, extent, cv.runs = 'abcde', f.class = c('LQP', 'Q'), beta.values = c(0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5) ){ summary.object <- expand.grid(CrossVal=cv.runs, Features = f.class, Betas = beta.values, stringsAsFactors = T) summary.object$cv.num <- as.numeric(summary.object$CrossVal) # ifelse functions defining what to do if taxon.name, time.bin, and extent are not specified if( !is.null(taxon.name) ){ summary.object$Taxa <- taxon.name } else { summary.object$Taxa <- 'Taxon1' } if( !is.null(time.bin) ){ summary.object$Time.Bin <- time.bin } else { summary.object$Time.Bin <- 'TimeBin1' } if( !is.null(extent) ){ summary.object$Extent <- extent } else { summary.object$Extent <- 'Extent1' } # re-ordering columns summary.object <- summary.object[,c(5:7, 1, 4, 2:3)] # adding in all the other parameters summary.object$n <- NA # sample size summary.object$k <- NA # number of non-zero lambdas summary.object$ln.L <- NA # log-likelihood summary.object$AIC <- NA # AIC summary.object$AICc <- NA # AICc corrected for small sample size summary.object$delta.i <- NA # delta.i for AIC - wont calculate everything until all models have been run summary.object$delta.i.c <- NA # delta.i for AICc - wont calculate everything until all models have been run summary.object$w.i <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run summary.object$w.i.c <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run summary.object$lambdas <- NA # list of all the non-zero lambdas summary.object[,4] <- as.character(summary.object[,4]) summary.object[,6] <- as.character(summary.object[,6]) return(summary.object) } # function to create the eval data.frame create.eval.df <- function(taxon.name, time.bin, extent, cv.runs = c('abcd', 'abce', 'abde', 'acde', 'bcde'), f.class = 'LQP', beta.values = 1 ){ summary.object <- expand.grid(CrossVal=cv.runs, Features = f.class, Betas = beta.values, stringsAsFactors = T) summary.object$cv.num <- 1:length(cv.runs) # ifelse functions defining what to do if taxon.name, time.bin, and extent are not specified if( !is.null(taxon.name) ){ summary.object$Taxa <- taxon.name } else { summary.object$Taxa <- 'Taxon1' } if( !is.null(time.bin) ){ summary.object$Time.Bin <- time.bin } else { summary.object$Time.Bin <- 'TimeBin1' } if( !is.null(extent) ){ summary.object$Extent <- extent } else { summary.object$Extent <- 'Extent1' } # re-ordering columns summary.object <- summary.object[,c(5:7, 1, 4, 2:3)] # adding in all the other parameters summary.object$n <- NA # sample size summary.object$k <- NA # number of non-zero lambdas summary.object$ln.L <- NA # log-likelihood summary.object$AIC <- NA # AIC summary.object$AICc <- NA # AICc corrected for small sample size summary.object$delta.i <- NA # delta.i for AIC - wont calculate everything until all models have been run summary.object$delta.i.c <- NA # delta.i for AICc - wont calculate everything until all models have been run summary.object$w.i <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run summary.object$w.i.c <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run summary.object$lambdas <- NA # list of all the non-zero lambdas summary.object[,4] <- as.character(summary.object[,4]) summary.object[,6] <- as.character(summary.object[,6]) # renaming the columns colnames(summary.object) <- c('Taxa', 'Time.Bin', 'Extent', 'CrossVal', 'cv.num', 'Features', 'Betas', 'n', 'thresh', 'test.sens', 'pROC_0.1', 'pval_0.1', 'pROC_1', 'pval_1', 'pROC_5', 'pval_5', 'lambdas') return(summary.object) } # function to create the nested fie structure that it will use to store the output (including html files) of maxent models # summary.eval.df <- the summary.df data.frame or the eval.df data.frame. # the function will use information from the cv.runs, f.class, and beta.values columns to create this nested file structure # wd <- working directory that the nested file structure is going to be generated in. if running multiple species, it is recommended . . . # . . . that you make a folder for each species, then run this function in each species folder create.folders.for.maxent <- function(summary.eval.df, wd = getwd() ){ # setting the working directory setwd( getwd() ) # for loop that goes through the summary/evaluation data.frame and creates the nested file structure for( i in 1:nrow(summary.eval.df) ){ if( !dir.exists( paste(getwd(), summary.eval.df$CrossVal[i], sep='/' ) ) ){ # if the cross-validation folder exists writeLines( c('Creating folder:', paste(getwd(), summary.eval.df$CrossVal[i], sep='/'), sep='') ) dir.create( paste(getwd(), summary.eval.df$CrossVal[i], sep='/') ) setwd( paste(getwd(), summary.eval.df$CrossVal[i], sep='/') ) } else { setwd( paste(getwd(), summary.eval.df$CrossVal[i], sep='/') ) } if( !dir.exists( paste(getwd(), summary.eval.df$Features[i], sep='/' ) ) ){ # if the feature class folder exists writeLines( c('Creating folder:', paste(getwd(), summary.eval.df$Features[i], sep='/'), sep='') ) dir.create( paste(getwd(), summary.eval.df$Features[i], sep='/') ) setwd( paste(getwd(), summary.eval.df$Features[i], sep='/') ) } else { setwd( paste(getwd(), summary.eval.df$Features[i], sep='/') ) } if( !dir.exists( paste(getwd(), summary.eval.df$Features[i], sep='/' ) ) ){ # if the regularization folder exists writeLines( c('Creating folder:', paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/'), sep='') ) dir.create( paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/') ) setwd( paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/') ) } else { setwd( paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/') ) } setwd('../../../') # go three directories back } # finishing main for loop } # setwd("~/Dropbox/SVP_Models/ModelOutput/Tyrano") # # # examp <- create.summary.df('Tyrano', 'K', 'Laurimidia') # examp2 <- create.eval.df('Tyrano', 'K', 'Laurimidia') # # # create.folders.for.maxent(examp) # create.folders.for.maxent(examp2) ##### Optimize MaxEnt Likelihood ------------------------------------------------------------------------------------------- # function that runs all the maxent models. # function will return a filled out summary model optimize.maxent.likelihood <- function(summary.df, # summary object that will keep all the model output # the nested file structure created from create.folders.for.maxent MUST exist occs, # species occurrence object background, # sampled Background object (COLUMNS MUST BE IDENTICAL TO occs) predic, # column numbers of the predictor variables first.occ.col, # number of the first cross-validation column in occs/background home=getwd(), # directory where all the models will be ran. all.models = TRUE # do you want to keep all versions of all the models? # helpful for quickly checking some models, but may consume loads (e.g., >1GB) . . . # . . . of hard disk space. Recommended to set to FALSE for exploratory analyses. # if FALSE, function will create a folder called "RunOver" and will . . . # . . . continuously write-over it for all models, and the only model you see . . . # . . . at the end will be the last model that was ran ){ # Calculating the total number of models nmodels <- nrow(summary.df) # prompting the user if they want to store models in the RAM print(paste0('The time is ', Sys.time(), '. You are running ', nmodels, ' total MaxEnt models.')) # setting up the progress bar prog <- txtProgressBar(min=0, max=nrow(summary.df), style=3, char='+') for(i in 1:nrow(summary.df)){ # setting the working directory for each folder if(all.models == TRUE){ setwd( paste(home, summary.df[i,4], summary.df[i,6], paste('beta', summary.df[i,7], sep='_'), sep='/') ) } else { # create the RunOver folder if( !dir.exists( paste(home, 'RunOver', sep='/') ) ){ dir.create( paste(home, 'RunOver', sep='/') ) setwd( paste(home, 'RunOver', sep='/') ) } else { setwd( paste(home, 'RunOver', sep='/') ) } } ########################################################################### ### MaxEnt things happen here ### preparing the data for the maxent model # filtering the occ object by it's respective cross-validation identity cv.number <- summary.df$cv.num[i] # assigning the column number to be sent through maxent col.number <- first.occ.col + cv.number - 1 # filtering the species dataset by col.number and assigning to summary.df sp <- occs %>% dplyr::filter(occs[,col.number] == 1) n <- nrow(sp) summary.df[i,8] <- n # bending the occ and background data together mx.data <- rbind(sp, background) ### running the actual maxent model mx.model <- dismo::maxent(p = mx.data[,col.number], x = mx.data[,predic], path = paste0(getwd()), args = prepPara(userfeatures = summary.df[i,6], betamultiplier = summary.df[i,7], doclamp = FALSE) ) ### calculating k and the names of the lambdas lambda.file <- as.data.frame(mx.model@lambdas) %>% `colnames<-`('lambdas') # lambdas data frame lambdas.df <- splitstackshape::cSplit(lambda.file, sep=',', splitCols='lambdas') colnames(lambdas.df) <- c('feature', 'lambda', 'min', 'max') # finding the non-zero lambdas non.zero.lambdas <- lambdas.df %>% dplyr::filter(!is.na(max)) %>% dplyr::filter(lambda != 0) class(non.zero.lambdas) <- 'data.frame' # assigning the number of parameters k <- nrow(non.zero.lambdas) # if beta is too high, the model gets over-regularized to the point that all the lambda coefficients get set to 0 # this effectively becomes an intercept only model, which is effectively the global mean # in this sense, k should get set to 0 if(k == 0){ k <- 1 } summary.df[i,9] <- k # giving noting the variables/features/hyperparameters with non-zero lambdas summary.df[i,17] <- toString(non.zero.lambdas[,1]) ### calculating the log-likelihood, then AIC and AICc # if statement calculating if there is an appropriate AIC value # e.g., can't fit 4 observations (occurrence points) with 5 variables if(n - k < 2){ summary.df[i,10] <- NA summary.df[i,11] <- NA summary.df[i,12] <- NA } else { # if it is possible to calculate AIC and AICc # logistic model output mx.back <- dismo::predict(mx.model, background[,predic]) # sum of all background point values - mx.back / back.sum should = 1 back.sum <- sum(mx.back) # logistic values of the (k-1)/k occurrences mx.occs <- dismo::predict(mx.model, sp[,predic]) # scaling to make compatible for calculating AIC occs.raw <- mx.occs / back.sum # log(likelihood) log.like <- sum(log(occs.raw)) summary.df[i,10]<- log.like ### calculating AIC and AICc # AIC AIC <- 2*k - 2*log.like summary.df[i,11] <- AIC # AICc summary.df[i,12] <- AIC + 2*((k^2 + k) / (n - k - 1)) } # updating the prograss bar for each run to get an idea of how long things will take setTxtProgressBar(prog, i) ########################################################################### # returning to the home directory setwd(home) } # closes the for loop return(summary.df) } ##### Null AICc ------------------------------------------------------------------------------------------------------------ # function that calculates AIC and AICc values for a null intercept-only model # the arguments are the same as the optimize.maxent.likelihood model, except that null.df should be a data.frame created from the . . . # . . . create.null.df object # function will return null.aic <- function(null.df, occs, background, first.occ.col){ # number of parameters null.df$k <- 1 # assigning the column number to be sent through maxent col.number <- first.occ.col # filtering the species dataset by col.number and assigning to null.df sp <- occs %>% dplyr::filter(occs[,col.number] == 1) n <- nrow(sp) null.df[1,8] <- n # giving noting the variables/features/hyperparameters with non-zero lambdas null.df[1,17] <- 'Intercept' out.scale <- rep(1, times=n ) / nrow(background) log.like <- sum(log(out.scale)) null.df[1,10] <- log.like ### calculating AIC and AICc # AIC AIC <- 2 - 2*log.like null.df[1,11] <- AIC # AICc for one term null model null.df[1,12] <- AIC + 4/(n - 2) return(null.df) } ##### MaxEnt Cross-Validation Error ----------------------------------------------------------------------------------------------- # function that runs five-fold cross-validation once you've found the optimum hyper-parameters for your maxent model(s) # functions returns a list containing 3 objects: 1.) a filled out eval object; 2.) an object containing all the maxent models; and . . . # . . . 3.) a data.frame with dimensions [ 1:nrow(background), 1:nrow(eval.df) ] containing the projections of all the models in maxent.crossval.error <- function(eval.df, # object generated from create.eval.df that will keep all the model output occs, # species occurrence object background, # sampled Background object (COLUMNS MUST BE IDENTICAL TO occs) predic, # column numbers of the predictor variables first.occ.col, # number of the first training in occs/background (assumes k = 5) # NOTE: this is not the presence column! it's the column to the right of it first.test.col, # number of the first testing column in occs/background (assumes k = 5) home=getwd(), # directory where all the models will be ran all.models = TRUE, # do you want to keep all versions of all the models? # helpful for quickly checking some models, but may consume loads (e.g., >1GB) . . . # . . . of hard disk space. Recommended to set to FALSE for exploratory analyses. # if FALSE, function will create a folder called "RunOver" and will . . . # . . . continuously write-over it for all models, and the only model you see . . . # . . . at the end will be the last model that was ran omission.rate=0, # user specified omission rate to be calculated for the threshold # express as a proportion from 0-1 all.background, # background data from the extents the species exists in # if using all the potential background points, all.background should == background pROC.reps = 500 # number of iterations each partialROC test will go through ){ # Calculating the total number of models nmodels <- nrow(eval.df) # stop if an incompatible omission rate is specified if(omission.rate > 1 || omission.rate < 0){ stop('Specify an omission rate between 0-1.') } # prompting the user if they want to store models in the RAM print(paste0('The time is ', Sys.time(), '. You are running ', nmodels, ' total MaxEnt models.')) # setting up the progress bar prog <- txtProgressBar(min=0, max=nrow(eval.df), style=3, char='+') # setting up various list objects to send output to # list for the maxent models maxent.list <- list() # list for the testing background data test.back.list <- list() # list for all the occs all.occs.list <- list() for(i in 1:nrow(eval.df)){ # setting the working directory for each folder if(all.models == TRUE){ setwd( paste(home, eval.df[i,4], eval.df[i,6], paste('beta', eval.df[i,7], sep='_'), sep='/') ) } else { # create the RunOver folder if( !dir.exists( paste(home, 'RunOver', sep='/') ) ){ dir.create( paste(home, 'RunOver', sep='/') ) setwd( paste(home, 'RunOver', sep='/') ) } else { setwd( paste(home, 'RunOver', sep='/') ) } } ########################################################################### ### MaxEnt things happen here ### preparing the data for the maxent model # filtering the occ object by it's respective cross-validation identity cv.number <- eval.df$cv.num[i] # assigning the column number to be sent through maxent col.number <- first.occ.col + cv.number - 1 # filtering the species dataset by col.number and assigning to eval.df sp <- occs %>% dplyr::filter(occs[,col.number] == 1) n <- nrow(sp) eval.df[i,8] <- n # generating the testing data test.col.number <- first.test.col + cv.number - 1 sp.test <- occs %>% dplyr::filter(occs[,test.col.number] == 1) # bending the occ and background data together mx.data <- rbind(sp, background) ### running the actual maxent model mx.model <- dismo::maxent(p = mx.data[,col.number], x = mx.data[,predic], path = paste0(getwd()), args = prepPara(userfeatures = eval.df[i,6], betamultiplier = eval.df[i,7], doclamp = FALSE) ) ### calculating k and the names of the lambdas lambda.file <- as.data.frame(mx.model@lambdas) %>% `colnames<-`('lambdas') # lambdas data frame lambdas.df <- splitstackshape::cSplit(lambda.file, sep=',', splitCols='lambdas') colnames(lambdas.df) <- c('feature', 'lambda', 'min', 'max') # finding the non-zero lambdas non.zero.lambdas <- lambdas.df %>% dplyr::filter(!is.na(max)) %>% dplyr::filter(lambda != 0) class(non.zero.lambdas) <- 'data.frame' # giving noting the variables/features/hyperparameters with non-zero lambdas eval.df[i,17] <- toString( non.zero.lambdas[,1] ) ## evaluating the models # predicting the training data mx.train.occ <- dismo::predict(mx.model, sp[,predic]) # predicting the training background data mx.train.back <- dismo::predict(mx.model, background[,predic]) # predicting the testing data mx.test.occ <- dismo::predict(mx.model, sp.test[,predic]) # projecting to all extents mx.test.back <- dismo::predict(mx.model, all.background[,predic]) # projecting to all occ points mx.all.occs <- dismo::predict(mx.model, occs[,predic]) # calculate the threshold if(omission.rate == 0){ # if using the LTP threshold, dismo calculates slightly too high of a threshold thresh <- min(mx.train.occ) } else { # make a model evaluation object mx.eval <- dismo::evaluate(p=mx.train.occ, a=mx.train.back) thresh <- dismo::threshold(mx.eval, stat='sensitivity', sensitivity= (1 - omission.rate) ) } # assigning the threshold value to the output file eval.df[i,9] <- thresh # calculating the test sensitivity sens <- length(which(mx.test.occ >= thresh)) / length(mx.test.occ) eval.df[i,10] <- sens # making a raster of the testing background data # for some reason, kuenm calculates wonky AUC_ratio values (i.e., > 2, which is impossible) unless you specify a raster r <- raster(nrows=1, ncols=length(mx.test.back) ) r[r] <- mx.test.back ## calculating AUC_ratios from a partialROC test # error = 0.1% pROC_0.1 <- kuenm::kuenm_proc(occ.test = mx.test.occ, # numeric vector of the predicted suitability values on the testing data model = r, # raster model of the predicted suitability values for the background threshold = 0.1, # potential error threshold (expressed as a percent) rand.percent = 50, # percentage of data to be used in each bootstrap rep iterations = pROC.reps # number of repititions ) # assigning the average AUC_ratio from pROC.reps iterations to eval.df eval.df[i,11] <- as.numeric(pROC_0.1$pROC_summary[1]) # assigning the partialROC p-value to eval.df eval.df[i,12] <- as.numeric(pROC_0.1$pROC_summary[2]) # error = 1% pROC_1 <- kuenm::kuenm_proc(occ.test = mx.test.occ, # numeric vector of the predicted suitability values on the testing data model = r, # raster model of the predicted suitability values for the background threshold = 1, # potential error threshold (expressed as a percent) rand.percent = 50, # percentage of data to be used in each bootstrap rep iterations = pROC.reps # number of repititions ) # assigning the average AUC_ratio from pROC.reps iterations to eval.df eval.df[i,13] <- as.numeric(pROC_1$pROC_summary[1]) # assigning the partialROC p-value to eval.df eval.df[i,14] <- as.numeric(pROC_1$pROC_summary[2]) # error = 5% pROC_5 <- kuenm::kuenm_proc(occ.test = mx.test.occ, # numeric vector of the predicted suitability values on the testing data model = r, # raster model of the predicted suitability values for the background threshold = 5, # potential error threshold (expressed as a percent) rand.percent = 50, # percentage of data to be used in each bootstrap rep iterations = pROC.reps # number of repititions ) # assigning the average AUC_ratio from pROC.reps iterations to eval.df eval.df[i,15] <- as.numeric(pROC_5$pROC_summary[1]) # assigning the partialROC p-value to eval.df eval.df[i,16] <- as.numeric(pROC_5$pROC_summary[2]) mx.test.back.df <- as.data.frame( as.matrix(mx.test.back, ncol=1) ) # assigning the objects to the various lists maxent.list[[i]] <- mx.model test.back.list[[i]] <- mx.test.back all.occs.list[[i]] <- mx.all.occs # updating the prograss bar for each run to get an idea of how long things will take setTxtProgressBar(prog, i) ########################################################################### # returning to the home directory setwd(home) } # closes the for loop # bind the testing background data into a single data back.projections <- as.data.frame( do.call('cbind', test.back.list ) ) colnames(back.projections) <- eval.df$CrossVal # bind all occs together occ.projections <- as.data.frame(do.call('cbind', all.occs.list)) colnames(occ.projections) <- eval.df$CrossVal # make a list of the output out.list <- list() out.list$maxent.models <- maxent.list out.list$back.projection <- back.projections out.list$occ.projection <- occ.projections out.list$summary <- eval.df return(out.list) } ##### MaxEnt Evaluation ---------------------------------------------------------------------------------------------------- # function that takes the eval object generated from maxent.crossval.error and calculates the weighted mean and standard deviation # the weighted mean is calculating my testing sensitivity (1 - omission rate) multiplied by the partial ROC/AUC value. # this ensures that if a model does not discriminate between presences/non-presences well, it will receive comparatively lower weight # later package versions will include Boyce index as an additional calibration technique and will offer the user the ability to . . . # . . . choose which metrics to use for assessing model reliability # the weighted mean and standard deviation can be plotted to infer model variability/uncertainty ### ARGUMENTS ### ## eval <- eval object generated from maxent.crossval.error that will project the model to every grid cell used in the training region ## pROC.error <- the user-specified omission rate. # (choose from 0.1, 1, 5 - however they almost always end up super correlated with each other) maxent.eval <- function(eval, pROC.error=1){ # extracting model sensitivity sens <- eval$summary$test.sens # which pROC error amount to use? if(pROC.error == 0.1){ AUC_ratio <- eval$summary$pROC_0.1 } else if(pROC.error == 1){ AUC_ratio <- eval$summary$pROC_1 } else if(pROC.error == 5){ AUC_ratio <- eval$summary$pROC_5 } else { stop('Select an appropriate partialROC error amount.') } # weights = sensitivity*AUC_ratio weights <- sens * AUC_ratio weights[is.nan(weights)] <- 0 # making a matrix of the background points models.mat <- as.matrix(eval$back.projection) # weighted means and standard deviations w.mean <- matrixStats::rowWeightedMeans(x=models.mat, w=weights) w.sd <- matrixStats::rowWeightedSds(x=models.mat, w=weights) # merging the weighted means/sd's of all background points back.df <- data.frame(w.mean, w.sd) # making a matrix of the background points occs.mat <- as.matrix(eval$occ.projection) # weighted means and standard deviations w.mean <- matrixStats::rowWeightedMeans(x=occs.mat, w=weights) w.sd <- matrixStats::rowWeightedSds(x=occs.mat, w=weights) # merging the weighted means/sd's of all occ points occ.df <- data.frame(w.mean, w.sd) # making the output list out.list <- list() out.list$back <- back.df out.list$occ <- occ.df out.list$weights <- weights return(out.list) } ##### Calculating the MaxEnt Threshold ------------------------------------------------------------------------------------- # function that calculates the threshold value for all the training data ### ARGUMENTS ### ## eval <- eval object generated from maxent.crossval.error that will project the model to every grid cell used in the training region ## occs <- the occurrence data.frame ## predic <- column numbers of the predictor variables ## pROC.error <- the user-specified omission rate. # (choose from 0.1, 1, 5 - however they almost always end up super correlated with each other) # function will return the lowest training preference threshold. # future package versions will give the opportunity to select different thresholds maxent.thresh <- function(eval, occs, predic, pROC.error=1){ n.rows <- nrow(occs) n.cols <- nrow(eval$summary) # making a data.frame of the occurrences occ.values <- as.data.frame( matrix(NA, nrow=n.rows, ncol=n.cols ) ) colnames(occ.values) <- eval$summary$CrossVal # for loop predicting all the cross validation maxent models for(i in 1:nrow(eval$summary) ){ mx.occs <- dismo::predict(object=eval$maxent.models[[i]], x=occs[,predic]) occ.values[,i] <- mx.occs } # # extracting model sensitivity sens <- eval$summary$test.sens # which pROC error amount to use? if(pROC.error == 0.1){ AUC_ratio <- eval$summary$pROC_0.1 } else if(pROC.error == 1){ AUC_ratio <- eval$summary$pROC_1 } else if(pROC.error == 5){ AUC_ratio <- eval$summary$pROC_5 } else { stop('Select an appropriate partialROC error amount.') } # weights = sensitivity*AUC_ratio weights <- sens * AUC_ratio weights[is.nan(weights)]<- 0 # making a matrix of the background points models.mat <- as.matrix(occ.values) # weighted means w.mean <- matrixStats::rowWeightedMeans(x=models.mat, w=weights) # setting the threshold as the lowest training presence ltp <- min(w.mean) return(ltp) } ##### Projecting the Best Models to All the Extents ------------------------------------------------------------------------ # function to project models to any and all extents you want to # effectively making new master occurrence and master background files that can easily be saved and re-loaded again maxent.everything <- function(eval, means, everything, thresh, pROC.error=1, predic, name='pc'){ # making a data.frame of the occurrences every.value <- as.data.frame( matrix(NA, nrow=nrow(everything), ncol=nrow(eval$summary)) ) colnames(every.value) <- eval$summary$CrossVal # for loop predicting all the cross validation maxent models for(i in 1:nrow(eval$summary) ){ mx.everything <- dismo::predict(object=eval$maxent.models[[i]], x=everything[,predic]) every.value[,i] <- mx.everything } # extracting model sensitivity sens <- eval$summary$test.sens # which pROC error amount to use? if(pROC.error == 0.1){ AUC_ratio <- eval$summary$pROC_0.1 } else if(pROC.error == 1){ AUC_ratio <- eval$summary$pROC_1 } else if(pROC.error == 5){ AUC_ratio <- eval$summary$pROC_5 } else { stop('Select an appropriate partialROC error amount.') } # weights = sensitivity*AUC_ratio weights <- sens * AUC_ratio weights[is.nan(weights)] <- 0 # making a matrix of the background points models.mat <- as.matrix(every.value) # weighted means w.mean <- matrixStats::rowWeightedMeans(x=models.mat, w=weights) # making the output data frame out.df <- as.data.frame( matrix(NA, nrow=length(w.mean), ncol=2) ) colnames(out.df) <- c(paste0(eval$summary[1,1], '_contin'), paste0(eval$summary[1,1], '_thresh') ) # assigning the weighted means to the output data.frame out.df[,1] <- w.mean # replicating the weighted means so they can be thresolded out.df[,2] <- w.mean # calculating the threshold thresh <- min(means$occ$w.mean) # thresolding the values out.df[,2][out.df[,2] >= thresh] <- 1 out.df[,2][out.df[,2] < thresh] <- 0 return(out.df) } ##### Informed Analysis ----------------------------------------------------------------------------------------------- # function that calculates Multivariate Environmental Suitability Surface (MESS), but allows for some extrapolation . . . # . . . based on how response curves look. informed.mess <- function(ref.extent, # reference data.frame - must contain same column names as 'data' mess.extent, # data.frame containing the extent to be projected to - ideally this should be all extents coord.cols=2:3, # column numbers of the coordinates of the data/ref.data (e.g., long/lat) predic=6:11, # column numbers of ONLY the final predictor variables tolerance=NULL, # tolerance vector to extrapolate beyond the limits of the data. see example below # if tolerance is not specified (default), function will calculate basic MESS keep.layers=FALSE # decision to retain the MESS values for all variables instead of just the final MESS data # changing to TRUE can help easily identify which variables are causing the MESS value . . . # . . . to be so low (i.e., it may be only 1 variable that is responsible) # can easily re-check this later ){ # isolating the coordinates and predictor variables new.vars <- mess.extent[, predic] ref.vars <- ref.extent[, predic] # if( !is.null(tolerance) ){ if(length(predic) != length(tolerance)/2 ){ stop('Ensure that tolerance is 2x as long as predic') } # changing tolerance to a data.frame tolerance <- as.data.frame( matrix(tolerance, ncol=2, byrow=TRUE) ) nvars <- length(predic) # pre-changing values outside the range of the input data # extra rows extra.rows <- matrix(NA, nrow=2, ncol=nvars) colnames(extra.rows) <- colnames(ref.vars) ref.vars <- rbind(ref.vars, extra.rows ) for(i in 1:nvars){ # if we can extrapolate beyond the lower end of variable i, allow mess to not if(tolerance[i,1] == 1){ min.var <- min(ref.vars[,i], na.rm=TRUE) - 0.0001 new.vars[ new.vars[,i] < min.var, i ] <- min.var ref.vars[nrow(ref.vars)-1, i] <- min.var } else { ref.vars[nrow(ref.vars)-1, i] <- min(ref.vars[,i], na.rm=TRUE) } if(tolerance[i,2] == 1){ max.var <- max(ref.vars[,i], na.rm=TRUE) + 0.0001 new.vars[ new.vars[,i] > max.var, i ] <- max.var ref.vars[nrow(ref.vars), i] <- max.var } else { ref.vars[nrow(ref.vars), i] <- max(ref.vars[,i], na.rm=TRUE) } } # end for(i in 1:nvars) loop } # running the mess analysis mess.vars <- as.data.frame( sapply(1:ncol(new.vars), function(i) .messi3(new.vars[, i], ref.vars[, i])) ) # re-asigning of the extrapolating points to have a mess value of 0. nref <- nrow(ref.vars) fix.vars <- as.data.frame( sapply(1:ncol(new.vars), function(i) .messi3(ref.vars[ (nref-1):nref , i], ref.vars[, i])) ) for(i in 1:nvars){ mess.vars[mess.vars[,i] == fix.vars[1,i],i] <- 0 mess.vars[mess.vars[,i] == fix.vars[2,i],i] <- 0 } # making a simple thresholded version of the mess analysis for easy plotting final.mess <- as.data.frame( apply(mess.vars, 1, min) ) colnames(final.mess) <- 'mess.raw' mess.thresh <- final.mess$mess.raw mess.thresh[mess.thresh > 0] <- 1 mess.thresh[mess.thresh < 0] <- -1 final.mess <- cbind(final.mess, mess.thresh) # deciding whether to keep individual mess layers if(keep.layers==TRUE){ colnames(ref.vars) <- paste0('mess_', colnames(ref.vars)) final.mess <- cbind(final.mess, ref.vars) return(final.mess) } else { return(final.mess) } } # internal function originally from dismo that runs the actual mess analysis in data.frame format .messi3 <- function(p,v) { # p=new.vars v=ref.vars # seems 2-3 times faster than messi2 v <- stats::na.omit(v) f <- 100*findInterval(p, sort(v)) / length(v) minv <- min(v) maxv <- max(v) res <- 2*f f[is.na(f)] <- -99 i <- f>50 & f<100 res[i] <- 200-res[i] i <- f==0 res[i] <- 100*(p[i]-minv)/(maxv-minv) i <- f==100 res[i] <- 100*(maxv-p[i])/(maxv-minv) res } ##### uncert.suit.plot ----------------------------------------------------------------------------------------------------- suit.uncert.plot <- function(means){ back <- means$back occs <- means$occ ltp <- min(occs$w.mean) output <- ggplot(data=back, aes(x=w.mean, y=w.sd)) + geom_point(colour='black', size=0.75) + geom_point(data=occs, aes(x=w.mean, y=w.sd), size=2, colour='red', shape=18 ) + ylim(0, 0.55) + xlim(0,1) + theme_classic() + coord_fixed(1/0.55) + annotate(geom='text', x=0.05, y=0.45, label=round(ltp, 4), hjust=0 ) return(output) } # function that adds k-folds to the dataset ### Arguments ### # x <- vector or data.frame of your occurrences # # k <- number of folds you want to make # # seed <- value for set.seed in case you want to reproduce your exact k-fold samples in the future # returns a vector that has the same length as nrow(x) that contains the k-fold bin that the occurrences got assigned to make.kfolds <- function(x, k=5, seed=NULL){ if(class(x) != "data.frame"){ class(x) <- "data.frame" } n <- nrow(x) rep.times <- n %/% k # number of full reps for the rep function k.bins <- rep.times*k # number of occs minus any remainder when dividing by k remainder <- n - k.bins # find the remainder output <- rep(1:k, rep.times) # make k bins of equal size if(!is.null(seed)){ set.seed(seed) } if(remainder != 0){ # if there is a remainder, . . . extra <- sample(x=1:k, size=remainder) # randomly sample it. . . output <- c(output, extra) # . . . and add it to the output } output <- sample(output, size=n) # randomize the order of the output return(output) } # function that makes a background file for each extent (time.bin and/or region) in the analysis # merging the background files for every extent will create the master background file ### ARGUMENTS ### # r.stack <- a raster stack containing all the predictor variables for analysis for a given extent # if using multiple extents, ensure that all the predictor variables have the exact same name and are in the same order # if testing multiple types of variables (i.e., GCM-based vs sedimentology-based), it is ideal to seperate them into . . . # . . . those two categories before reading them into R # # x.col <- the column name of the x-coordinate to be put in the background/occurrence objects # this MUST be consistent for all background/occurrence objects # # y.col <- the column name of the y-coordinate to be put in the background/occurrence objects # this MUST be consistent for all background/occurrence objects # # time.bin <- the name of the time.bin to be put in the background/occurrence objects # if only using a single time.bin, keep this consistent for all background/occurrence objects # # region <- the name of the region to be put in the background/occurrence objects # if only using a single region, keep this consistent for all background/occurrence objects # # k.folds <- the number of folds you want to run for cross-validation # k stops at 26 because you don't want to run 27+ fold cross-validation. The resulting data.frame becomes unwieldy. make.background.df <- function(r.stack, x.col='long', y.col='lat', time.bin='time1', region=1, k.folds=5){ # ensuring that k.folds is a positive integer between [1,26] if( !is.null(k.folds) ){ if( k.folds < 1){ k.folds <- 1 print('NOTE: Making background data.frame without any k.folds.') } else if( k.folds > 26 ){ k.folds <- 26 warning('k.folds has been set to 26.') } else if( k.folds%%1 != 0 ){ warning('k.folds has been rounded to the nearest whole number.') k.folds <- round(k.folds) } } else { k.folds <- 1 print('NOTE: Making background data.frame without any k.folds.') } # extracting the env-variables to a data.frame coords.df <- raster::sampleRandom(r.stack, ncell(r.stack), xy=TRUE, sp=FALSE, na.rm=FALSE) colnames(coords.df)[1:2] <- c(x.col, y.col) coords.df <- coords.df[complete.cases(coords.df), ] n.points <- nrow(coords.df) # making the data.frame of the name, time.bin, region, and presence columns d.cols <- as.data.frame( matrix(0, nrow=n.points, ncol=4) ) d.cols[,1] <- 'Background' d.cols[,2] <- time.bin d.cols[,3] <- region colnames(d.cols) <- c('Name', 'time.bin', 'region', 'presence') # merging everything together out <- do.call('cbind', list( d.cols[,1], coords.df[,1:2], d.cols[,2:3], coords.df[,3:ncol(coords.df)], d.cols[,4] )) colnames(out)[1] <- 'Name' colnames(out)[ ncol(out) ] <- 'presence' # making the cross-validation columns if( k.folds > 1 ){ ltrs <- letters[1:k.folds] # training columns train <- as.data.frame( matrix(0, nrow=n.points, ncol=k.folds) ) for(i in 1:k.folds){ colnames(train)[i] <- paste0('tr.', paste(ltrs[-(k.folds+1-i)], collapse='') ) } # testing columns test <- as.data.frame( matrix(0, nrow=n.points, ncol=k.folds) ) for(i in 1:k.folds){ colnames(test)[i] <- paste0('test.', paste(ltrs[(k.folds+1-i)], collapse='') ) } out <- do.call('cbind', list( out, train, test ) ) } # output return(out) } # function that makes an occurrence file for each extent (time.bin and/or region) in the analysis # merging the occurrence files for every extent will create the master occurrence file # MUST have the make.kfolds function also loaded ### ARGUMENTS ### # r.stack <- a raster stack containing all the predictor variables for analysis for a given extent # if using multiple extents, ensure that all the predictor variables have the exact same name and are in the same order # if testing multiple types of variables (i.e., GCM-based vs sedimentology-based), it is ideal to seperate them into . . . # . . . those two categories before reading them into R # # taxa.df <- a data.frame of occurrence with three columns that have: # 1.) the names of the taxa you are modeling # 2.) the x-coordinates of the occurrences # 3.) the y-coordinates of the occurrences # THESE COLUMNS MUST BE IN THIS ORDER!!! # this is the same format as the SWD (species with data) format for the regular maxent.jar file # # x.col <- the column name of the x-coordinate to be put in the background/occurrence objects # this MUST be consistent for all background/occurrence objects # # y.col <- the column name of the y-coordinate to be put in the background/occurrence objects # this MUST be consistent for all background/occurrence objects # # time.bin <- the name of the time.bin to be put in the background/occurrence objects # if only using a single time.bin, keep this consistent for all background/occurrence objects # # region <- the name of the region to be put in the background/occurrence objects # if only using a single region, keep this consistent for all background/occurrence objects # # k.folds <- the number of folds you want to run for cross-validation # k stops at 26 because you don't want to run 27+ fold cross-validation. The resulting data.frame becomes unwieldy. # # k.seed <- value for set.seed in case you want to reproduce your exact k-fold samples in the future # NOTE that this set.seed will apply exactly the same to every species you are modeling # make.occurrence.df <- function(r.stack, taxa.df, x.col='long', y.col='lat', time.bin='time1', region=1, k.folds=5, k.seed=NULL){ # ensuring that k.folds is a positive integer between [1,26] if( !is.null(k.folds) ){ if( k.folds < 1){ k.folds <- 1 print('NOTE: Making occurrence data.frame without any k.folds.') } else if( k.folds > 26 ){ k.folds <- 26 warning('k.folds has been set to 26.') } else if( k.folds%%1 != 0 ){ warning('k.folds has been rounded to the nearest whole number.') k.folds <- round(k.folds) } } else { k.folds <- 1 print('NOTE: Making occurrence data.frame without any k.folds.') } # extracting the env-variables to a data.frame occs.df <- raster::extract(x=r.stack, y=taxa.df[,2:3]) occs.df <- cbind(taxa.df, occs.df) colnames(occs.df)[2:3] <- c(x.col, y.col) occs.df <- occs.df[complete.cases(occs.df), ] colnames(occs.df)[1] <-'Name' n.occs <- nrow(occs.df) # making the data.frame of the name, time.bin, region, and presence columns d.cols <- as.data.frame( matrix(0, nrow=n.occs, ncol=3) ) d.cols[,1] <- time.bin d.cols[,2] <- region colnames(d.cols) <- c('time.bin', 'region', 'presence') # merging everything together occs.df <- do.call('cbind', list(occs.df[,1:3], d.cols[,1:2], occs.df[,4:ncol(occs.df)], d.cols[,3] ) ) colnames(occs.df)[ ncol(occs.df) ] <- 'presence' # making the cross-validation columns if( k.folds > 1 ){ ltrs <- letters[1:k.folds] # training columns train <- as.data.frame( matrix(1, nrow=n.occs, ncol=k.folds) ) for(i in 1:k.folds){ colnames(train)[i] <- paste0('tr.', paste(ltrs[-(k.folds+1-i)], collapse='') ) } # testing columns test <- as.data.frame( matrix(0, nrow=n.occs, ncol=k.folds) ) for(i in 1:k.folds){ colnames(test)[i] <- paste0('test.', paste(ltrs[(k.folds+1-i)], collapse='') ) } # splitting up the occurrence data.frame by taxon, then assigning the k.folds occs.list <- base::split(occs.df, occs.df[,1]) n.taxa <- length(occs.list) for(i in 1:n.taxa){ # setting all taxa with < k.folds occurrences to 0's. they will later be converted to exist in all the training and testing subsets if( nrow( occs.list[[i]] ) < k.folds ){ occs.list[[i]]$presence <- 0 } else { # filling in cross-validation folds for taxa who have > k.folds in terms of occurrences occs.list[[i]]$presence <- make.kfolds(occs.list[[i]], k=k.folds, seed=k.seed) } } # closing for(i in 1:n.taxa) # merging all the different taxa back together occs.df <- do.call('rbind', occs.list) colnames(occs.df)[ ncol(occs.df) ] <- 'presence' row.names(occs.df) <- 1:nrow(occs.df) # assigning the k.fold parameters in occs.df for(i in 1:nrow(occs.df) ){ # for taxa with fewer occs than k.folds if( occs.df$presence[i] == 0 ){ # train[i,] <- 1 # train is filled with 1's by default test[i,] <- 1 } else { # for taxa with greater occs than k.folds cv <- occs.df$presence[i] train[i, (k.folds+1-cv) ] <- 0 test[i, (k.folds+1-cv) ] <- 1 } } # closing for(i in 1:nrow(occs.df) ) out <- do.call('cbind', list( occs.df, train, test ) ) } # closing if( k.folds > 1 ) # setting all the presence rows to 1 out$presence <- 1 # output return(out) } # function that reduces occurrences to one per grid cell ### ARGUMENTS ### ## occ.data = data.frame of the occurrence file ## rast = raster of the entire background training extent # there will be issues if any occs are in grid cells with NA values (i.e., outside the extent) ## name = name of the column that has the taxa names ## long= column name that has longitude ## lat = column name that has latitude ## max.dist = argument from seegSDM. distance is in map units (e.g., degrees) if the raster is projected, otherwise, it is in meters. # if any occurrences lie JUST BARELY outside the extent, this function will assign them to the nearest gril cell . . . # . . . inside the extent if the distance from that occurrence to the nearest cell is <= max.dist. # otherwise, that point is ignored. # recommended that you remove points outside the extent first. it's easier to clear up that way ## round.to is how many decimal places to round the coordinates to # fewer decimal plaes = faster run times # nearestLand function extracted from seegSDM as that package doesn't appear to exist for more versions (> 3.6.2) of R # this version of nearestLand is from seegSDM version 0.1-9 seegSDM_nearestLand <- function (points, raster, max_distance) { nearest <- function(lis, raster) { neighbours <- matrix(lis[[1]], ncol = 2) point <- lis[[2]] land <- !is.na(neighbours[, 2]) if (!any(land)) { return(c(NA, NA)) } else { coords <- xyFromCell(raster, neighbours[land, 1]) if (nrow(coords) == 1) { return(coords[1, ]) } dists <- sqrt((coords[, 1] - point[1])^2 + (coords[, 2] - point[2])^2) return(coords[which.min(dists), ]) } # ending else } # ending nearest neighbour_list <- extract(raster, points, buffer = max_distance, cellnumbers = TRUE) neighbour_list <- lapply(1:nrow(points), function(i) { list(neighbours = neighbour_list[[i]], point = as.numeric(points[i, ])) }) return(t(sapply(neighbour_list, nearest, raster))) } # one.occ.per.grid.cell.no.SEEG <- function(occ.data, rast, name = 'species', long = 'long', lat = 'lat', max.dist = .0833, round.to = 5) {# I CHANGED MAX DIST TO .083 because 5 arc minutes? require(raster) occ.mat <- as.matrix(occ.data[,c(long, lat)]) occ.mat <- round(occ.mat, digits = round.to) moved <- seegSDM_nearestLand(occ.mat, raster=rast, max_distance=max.dist) # centers occ. points within the grid cell they occur in moved <- as.data.frame(moved) moved <- cbind(occ.data[,name], moved) # bind names and the paleo-longitudes colnames(moved) <- c(name, 'long.thin', 'lat.thin') moved <- unique(moved) # returns only one occurrence per entity per grid numbs <- as.numeric(row.names(moved)) # row numbers of the thinned data out <- occ.data[numbs,] # thinning the original data with the rows of the thinned data return(out) } #################&&&&&&&&&&&&&&&&&&&&&############MAKE MASTER################# #################### ACER PRES MASTER # reading in the occurrence file Pres.occs <- read.csv('AcroporaCervicornisPresent.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Holo.occs <- read.csv('AcroporaCervicornisHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Pres.occs <- Pres.occs[,-1] Holo.occs <- Holo.occs[,-1] head(Pres.occs) head(Holo.occs) identical(colnames(Pres.occs), colnames(Holo.occs)) nrow(Pres.occs) nrow(Holo.occs) Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3) Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3) #assigning wgs84 projection if not already assigned in layers wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present") Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$') Prescomp.files # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas) #Pres.files <- Pres.files[ c(1,2,3,4,5) ] Prescomp.rasters <- raster::stack( Prescomp.files) # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED") Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$') Holo.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size Holo.rasters <- raster::stack( Holo.files) Holo.rasters <- Holo.rasters/100 # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) # re-naming the names of rasters if you want to names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Prescomp.rasters), names(Holo.rasters) ) ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050") RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$') RCP452050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452050.rasters <- raster::stack( RCP452050.files) # optionally assigning a crs if one isn't provided crs(RCP452050.rasters) <- wgs1984 names(RCP452050.rasters) # re-naming the names of turo.rasters if you want to names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452050.rasters) ) # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050") RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$') RCP852050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852050.rasters <- raster::stack( RCP852050.files) # optionally assigning a crs if one isn't provided crs(RCP852050.rasters) <- wgs1984 names(RCP852050.rasters) # re-naming the names of turo.rasters if you want to names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852050.rasters) ) # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100") RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$') RCP452100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452100.rasters <- raster::stack( RCP452100.files) # optionally assigning a crs if one isn't provided crs(RCP452100.rasters) <- wgs1984 names(RCP452100.rasters) names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452100.rasters) ) # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100") RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$') RCP852100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852100.rasters <- raster::stack( RCP852100.files) # optionally assigning a crs if one isn't provided crs(RCP852100.rasters) <- wgs1984 names(RCP852100.rasters) names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852100.rasters) ) # making the background files Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean') RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean') RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean') RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean') #spatial thinning Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat") Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat") #how many occurrences are there - compare to original before spatial thinning to see if less nrow(Pres.occs.thin) nrow(Holo.occs.thin) #look at new spatially thinned dataframe View(Pres.occs.thin) # making the occurrence files # c(8,2,1) = (GENUS, Longitude, Latitude) Pres.occs.df <- make.occurrence.df(r.stack=Prescomp.rasters, taxa.df=Pres.occs.thin[,c(3,1,2)], #I think this bit is reorganizing the columns and I made it 3 not 8 and switched 2 1 for 1 2 x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)], x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') # checking to see if everything has identical column names for your different extents identical( colnames(Prescomp.back), colnames(Holo.back) ) identical( colnames(Pres.occs.df), colnames(Holo.occs.df) ) setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/Master_files_acerv_add_mask") masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back, RCP452050.back,RCP852050.back, RCP452100.back, RCP852100.back ) ) write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE) masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) ) View(masterall.occs) write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE) masterpres.occs <- Pres.occs.df write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE) masterholo.occs <- Holo.occs.df write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE) #################################### ACER HOLO MASTER same code but with area masked for cuba jamaica and hispanola #################################### APAL PRES MASTER # reading in the occurrence file Pres.occs <- read.csv('AcroporaPalmataPresent2.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Holo.occs <- read.csv('AcroporaPalmataHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Pres.occs <- Pres.occs[,-1] Holo.occs <- Holo.occs[,-1] head(Pres.occs) head(Holo.occs) identical(colnames(Pres.occs), colnames(Holo.occs)) nrow(Pres.occs) nrow(Holo.occs) Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3) Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3) ################################################################################################################################################ #assignin wgs84 projection if not already assigned in layers wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$') Prescomp.files # re-ordering any files # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas) #Pres.files <- Pres.files[ c(1,2,3,4,5) ] Prescomp.rasters <- raster::stack( Prescomp.files) # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam', pattern='\\.tif$') Holo.files #reorder Holo.files <- Holo.files[ c(1,3,2,4) ] Holo.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size Holo.rasters <- raster::stack( Holo.files) Holo.rasters <- Holo.rasters/100 # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED") Holo2.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$') Holo2.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size Holo2.rasters <- raster::stack( Holo2.files) Holo2.rasters <- Holo2.rasters/100 # optionally assigning a crs if one isn't provided crs(Holo2.rasters) <- wgs1984 names(Holo2.rasters) # re-naming the names of rasters if you want to names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) names(Holo2.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Prescomp.rasters), names(Holo.rasters) ) # ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$') RCP452050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452050.rasters <- raster::stack( RCP452050.files) # optionally assigning a crs if one isn't provided crs(RCP452050.rasters) <- wgs1984 names(RCP452050.rasters) # re-naming the names of turo.rasters if you want to names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452050.rasters) ) # # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050") RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$') RCP852050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852050.rasters <- raster::stack( RCP852050.files) # optionally assigning a crs if one isn't provided crs(RCP852050.rasters) <- wgs1984 names(RCP852050.rasters) # re-naming the names of turo.rasters if you want to names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852050.rasters) ) # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100") RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$') RCP452100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452100.rasters <- raster::stack( RCP452100.files) # optionally assigning a crs if one isn't provided crs(RCP452100.rasters) <- wgs1984 names(RCP452100.rasters) names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452100.rasters) ) # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100") RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$') RCP852100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852100.rasters <- raster::stack( RCP852100.files) # optionally assigning a crs if one isn't provided crs(RCP852100.rasters) <- wgs1984 names(RCP852100.rasters) names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852100.rasters) ) # # # making the background files Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') Holo2.back <- make.background.df(r.stack=Holo2.rasters, x.col='long', y.col='lat', time.bin='Holocene2', region='Caribbean') # RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean') # RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean') # RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean') # RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean') #spatial thinning Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat") Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat") #how many occurrences are there - compare to original before spatial thinning to see if less nrow(Pres.occs.thin) nrow(Holo.occs.thin) #look at new spatially thinned dataframe View(Pres.occs.thin) # making the occurrence files # c(8,2,1) = (GENUS, Longitude, Latitude) Pres.occs.df <- make.occurrence.df(r.stack=Prescomp.rasters, taxa.df=Pres.occs.thin[,c(3,1,2)], #I think this bit is reorganizing the columns and I made it 3 not 8 and switched 2 1 for 1 2 x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)], x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') # checking to see if everything has identical column names for your different extents identical( colnames(Prescomp.back), colnames(Holo.back) ) identical( colnames(Pres.occs.df), colnames(Holo.occs.df) ) setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/Master_files_apal_add_mask") # merging together the files masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back, Holo2.back, RCP452050.back,RCP852050.back, RCP452100.back, RCP852100.back ) ) write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE) masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) ) View(masterall.occs) write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE) masterpres.occs <- Pres.occs.df write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE) masterholo.occs <- Holo.occs.df write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE) ########################### APAL HOLO MASTER same but with masked background #################### CNAT PRES MASTER setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/Cnatans") # reading in the occurrence file Pres.occs <- read.csv('CnatansPresent.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Holo.occs <- read.csv('CnatansHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Pres.occs <- Pres.occs[,-1] Holo.occs <- Holo.occs[,-1] head(Pres.occs) head(Holo.occs) identical(colnames(Pres.occs), colnames(Holo.occs)) nrow(Pres.occs) nrow(Holo.occs) Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3) Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3) ################################################################################################################################################ #assignin wgs84 projection if not already assigned in layers wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present") Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$') Prescomp.files Prescomp.rasters <- raster::stack( Prescomp.files) # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED") Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$') Holo.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size Holo.rasters <- raster::stack( Holo.files) #scaling the holocene raster Holo.rasters <- Holo.rasters/100 # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) # re-naming the names of rasters if you want to names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Prescomp.rasters), names(Holo.rasters) ) ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050") RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$') RCP452050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452050.rasters <- raster::stack( RCP452050.files) # optionally assigning a crs if one isn't provided crs(RCP452050.rasters) <- wgs1984 names(RCP452050.rasters) # re-naming the names of turo.rasters if you want to names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452050.rasters) ) # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050") RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$') RCP852050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852050.rasters <- raster::stack( RCP852050.files) # optionally assigning a crs if one isn't provided crs(RCP852050.rasters) <- wgs1984 names(RCP852050.rasters) # re-naming the names of turo.rasters if you want to names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852050.rasters) ) # # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100") RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$') RCP452100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452100.rasters <- raster::stack( RCP452100.files) # optionally assigning a crs if one isn't provided crs(RCP452100.rasters) <- wgs1984 names(RCP452100.rasters) names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452100.rasters) ) # # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100") RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$') RCP852100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852100.rasters <- raster::stack( RCP852100.files) # optionally assigning a crs if one isn't provided crs(RCP852100.rasters) <- wgs1984 names(RCP852100.rasters) names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852100.rasters) ) # # # making the background files Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean') RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean') RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean') RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean') #spatial thinning Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat") Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat") #how many occurrences are there - compare to original before spatial thinning to see if less nrow(Pres.occs.thin) nrow(Holo.occs.thin) #look at new spatially thinned dataframe View(Pres.occs.thin) # making the occurrence files # c(8,2,1) = (GENUS, Longitude, Latitude) Pres.occs.df <- make.occurrence.df(r.stack=Prescomp.rasters, taxa.df=Pres.occs.thin[,c(3,1,2)], #I think this bit is reorganizing the columns and I made it 3 not 8 and switched 2 1 for 1 2 x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)], x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') # checking to see if everything has identical column names for your different extents identical( colnames(Prescomp.back), colnames(Holo.back) ) identical( colnames(Pres.occs.df), colnames(Holo.occs.df) ) identical( colnames(Prescomp.back), colnames(Holo.occs.df) ) setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/Master_files_cnat_add_mask") # merging together the files masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back, RCP452050.back,RCP852050.back, RCP452100.back, RCP852100.back) ) write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE) masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) )#LGM.occs View(masterall.occs) write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE) masterpres.occs <- Pres.occs.df write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE) masterholo.occs <- Holo.occs.df write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE) #################### CNAT HOLO MASTER same with masked extent #################### PAST PRES MASTER setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/Pasteroides") # reading in the occurrence file Pres.occs <- read.csv('PoritesAsteroidesPresent.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Holo.occs <- read.csv('PoritesAsteroidesHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1") Pres.occs <- Pres.occs[,-1] Holo.occs <- Holo.occs[,-1] head(Pres.occs) head(Holo.occs) identical(colnames(Pres.occs), colnames(Holo.occs)) nrow(Pres.occs) nrow(Holo.occs) Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3) Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3) ################################################################################################################################################ #assignin wgs84 projection if not already assigned in layers wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present") Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$') Prescomp.files # re-ordering any files # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas) #Pres.files <- Pres.files[ c(1,2,3,4,5) ] Prescomp.rasters <- raster::stack( Prescomp.files) # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED") Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$') Holo.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size Holo.rasters <- raster::stack( Holo.files) Holo.rasters <- Holo.rasters/100 # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) # re-naming the names of rasters if you want to names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Prescomp.rasters), names(Holo.rasters) ) ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050") RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$') RCP452050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452050.rasters <- raster::stack( RCP452050.files) # optionally assigning a crs if one isn't provided crs(RCP452050.rasters) <- wgs1984 names(RCP452050.rasters) # re-naming the names of turo.rasters if you want to names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452050.rasters) ) # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++= setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050") RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$') RCP852050.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852050.rasters <- raster::stack( RCP852050.files) # optionally assigning a crs if one isn't provided crs(RCP852050.rasters) <- wgs1984 names(RCP852050.rasters) # re-naming the names of turo.rasters if you want to names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" ) #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852050.rasters) ) # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100") RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$') RCP452100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP452100.rasters <- raster::stack( RCP452100.files) # optionally assigning a crs if one isn't provided crs(RCP452100.rasters) <- wgs1984 names(RCP452100.rasters) names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP452100.rasters) ) # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100") RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$') RCP852100.files #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size RCP852100.rasters <- raster::stack( RCP852100.files) # optionally assigning a crs if one isn't provided crs(RCP852100.rasters) <- wgs1984 names(RCP852100.rasters) names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp") #make sure e-layer names of each extent (time bin or region) is the same name and same order; if it is, will say TRUE identical( names(Holo.rasters), names(RCP852100.rasters) ) # # # making the background files Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean') RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean') RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean') RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean') #spatial thinning Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat") Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat") #how many occurrences are there - compare to original before spatial thinning to see if less nrow(Pres.occs.thin) nrow(Holo.occs.thin) #look at new spatially thinned dataframe View(Pres.occs.thin) # making the occurrence files # c(8,2,1) = (GENUS, Longitude, Latitude) Pres.occs.df <- make.occurrence.df(r.stack=Prescomp.rasters, taxa.df=Pres.occs.thin[,c(3,1,2)], #I think this bit is reorganizing the columns and I made it 3 not 8 and switched 2 1 for 1 2 x.col='long', y.col='lat', time.bin='Present', region='Caribbean') Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)], x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean') # checking to see if everything has identical column names for your different extents identical( colnames(Prescomp.back), colnames(Holo.back) ) identical( colnames(Pres.occs.df), colnames(Holo.occs.df) ) setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/Master_files_past_add_mask") # merging together the files masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back, RCP452050.back,RCP852050.back, RCP452100.back, RCP852100.back) ) write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE) masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) ) View(masterall.occs) write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE) masterpres.occs <- Pres.occs.df write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE) masterholo.occs <- Holo.occs.df write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE) #################### PAST HOLO MASTER same with masked extent #################&&&&&&&&&&&&&&&&&&&&&############ECOSPAT################# ####################################### ACER ECOSPAT ########################### #### MODEL COMPARISON #### #need cleaned occurrence data and environmental data rasters ## Read in shapefiles of cleaned occ data ## ########################occurrence data###################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/Master_files_acerv_add_mask") AcervPres = read.csv("masterpres.occs.csv") AcervPres = as(AcervPres,'data.frame') AcervHolo = read.csv("masterholo.occs.csv") AcervHolo = as(AcervHolo,'data.frame') #convert to dataframe for use in ecospat ###################environmental################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present") # Upload basic rasters for First interval (Present) Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$') Prescomp.files maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif") maxsalinitypres maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif") maxtemppres rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif") rangesalinitypres rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif") rangetemppres # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas) Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres)) wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') #Upload basic rasters for Second interval (Holocene) ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam") Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$') #scaling holocene layers maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif") maxsalinityholo2 <- maxsalinityholo/100 maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif") maxtempholo2 <- maxtempholo/100 rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif") rangesalinityholo2 <- rangesalinityholo/100 rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif") rangetempholo2 <- rangetempholo/100 Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2) # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') ### Ecospat Analysis of Niche Equivalency and Similary #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Set working directory for saving plots and test results setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/Ecospat_acerv_10_2_24") combine.lat = c(AcervPres$lat, AcervHolo$lat) combine.lon = c(AcervPres$long, AcervHolo$long) #I need to turn the points into a spatvector ( library(terra) ptspres <- vect(cbind(AcervPres$long, AcervPres$lat)) ptsholo <- vect(cbind(AcervHolo$long, AcervHolo$lat)) # get random sample - background point locations - rasterobject <- rast(Prescomp.rasters) bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) rasterobject <- rast(Holo.rasters) bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) # get environmental data from occ points extractpres = na.omit(cbind(AcervPres[2:3], raster::extract(Prescomp.rasters, AcervPres[2:3]), rep(1, nrow(AcervPres)))) extractholo = na.omit(cbind(AcervHolo[2:3], raster::extract(Holo.rasters, AcervHolo[2:3]), rep(1, nrow(AcervHolo)))) #change name of last column to occ (represents presence = 1) colnames(extractpres)[ncol(extractpres)] = 'occ' colnames(extractholo)[ncol(extractholo)] = 'occ' # get environmental data from bg points extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres))))) extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo))))) # change name of last column to occ (represents absence = 0) colnames(extbgpres)[ncol(extbgpres)] = 'occ' colnames(extbgholo)[ncol(extbgholo)] = 'occ' colnames(extbgpres)[1] = 'long' colnames(extbgpres)[2] = 'lat' colnames(extbgholo)[1] = 'long' colnames(extbgholo)[2] = 'lat' # merge data from occ and bg datPresAcerv = rbind(extractpres, extbgpres) datHoloAcerv = rbind(extractholo, extbgholo) #principle component analysis pca.env <- dudi.pca( rbind(datHoloAcerv, datPresAcerv)[,3:6], scannf=FALSE, nf=2 ) # look at variable contribution ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) # Save pdf of PCA plot pdf('ecospat_Acervpresholo.pdf') ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) scores.globclim<-pca.env$li # PCA scores for the whole study area scores.globclim<-pca.env$li # PCA scores for the whole study area (all points) scores.sppres <- suprow(pca.env, extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution scores.spholo <- suprow(pca.env, extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution scores.climholo <- suprow(pca.env,datHoloAcerv[,3:6])$li # PCA scores for the first interval study area scores.climpres <- suprow(pca.env,datPresAcerv[,3:6])$li # PCA scores for the second interval study area # density distribtuion for first interval grid.Acervholo <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climholo, sp = scores.spholo, R = 100, th.sp = 0 ) # density distribution for second interval grid.Acervpres <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climpres, sp = scores.sppres, R = 100, th.sp = 0 ) histpres <- hist(grid.Acervpres$glob1) histholo <- hist(grid.Acervholo$glob1) plot( histpres, col=rgb(0,0,1,1/4)) # first histogram plot( histholo, col=rgb(1,0,0,1/4), add=T) # Schoener's D metric and I metric D.overlap <- ecospat.niche.overlap (grid.Acervholo, grid.Acervpres, cor=T) D.overlap write.csv(D.overlap,"presAcerv_holoAcerv_I_D.csv",row.names=FALSE) ## Niche Equivalency Test ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D ## p.I: pvalue of the test on I ## Test for greater or lower equivalency eq.testgr <- ecospat.niche.equivalency.test(grid.Acervholo, grid.Acervpres, rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs eq.testlw <- ecospat.niche.equivalency.test(grid.Acervholo, grid.Acervpres, rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs nichdynindex <- ecospat.niche.dyn.index(grid.Acervholo, grid.Acervpres, intersection = NA) #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w nichdynindex # write p values of equivalency test p_EQ_DI = cbind("p.D_GR"=eq.testgr$p.D,"p.I_GR"=eq.testgr$p.I,"p.D_LW"=eq.testlw$p.D,"p.I_LW"=eq.testlw$p.D) write.csv(p_EQ_DI,"glycim_EQ_TestAcervpresholo.csv",row.names=FALSE) p_EQ_DI ## Test for greater (niche conservatism) or lower (niche divergence) similarity sim.testgr <- ecospat.niche.similarity.test(grid.Acervholo, grid.Acervpres, rep=1000, overlap.alternative = "higher", rand.type=2,ncores=4) sim.testlw <- ecospat.niche.similarity.test(grid.Acervholo, grid.Acervpres, rep=1000, overlap.alternative = "lower", rand.type=2,ncores=4) # write p values of similarity test p_SIM_DI = cbind("p.D_GR"=sim.testgr$p.D,"p.I_GR"=sim.testgr$p.I,"p.D_LW"=sim.testlw$p.D,"p.I_LW"=sim.testlw$p.D) write.csv(p_SIM_DI,"presholoAcerv_SIM_Test.csv",row.names=FALSE) p_SIM_DI # Plot test distributions ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency") ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency") ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity") ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity") ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals ecospat.plot.niche (grid.Acervpres, title='Acervpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE) ecospat.plot.niche (grid.Acervholo, title='Acervholo', name.axis1='PC1', name.axis2='PC2') ## Plot niche overlap ecospat.plot.niche.dyn(z1 = grid.Acervholo, z2 = grid.Acervpres, quant=0.25, interest=2, title= "A. cerv Holo to Pres", name.axis1="PC1", name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = '#74add1', colZ1 = "#313695", colZ2 = "#abd9d9", transparency = 40) ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white') rstudioapi::savePlotAsImage("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/Ecospat_acerv_10_2_24/figures/A. cerv Holo to Pres.png",width=1000,height=750) # Save pdf of plots ### Look at niche expantion, stability, and unfilling # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Acervholo, z2 = grid.Acervpres, intersection=NA) dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Acervholo, z2 = grid.Acervpres, intersection=0) # write csv of niche dynmaic percentages Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w, "Overlapping"=dynam_overlap$dynamic.index.w) Niche_Ex_St write.csv(Niche_Ex_St,"GlycymDynamicsAcervholopres.csv",row.names=TRUE) ############################################ APAL ECOSPAT ################### ## Read in shapefiles of cleaned occ data ## ########################occurrence data###################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/Master_files_apal_add_mask") ApalPres = read.csv("masterpres.occs.csv") ApalPres = as(ApalPres,'data.frame') ApalHolo = read.csv("masterholo.occs.csv") ApalHolo = as(ApalHolo,'data.frame') #convert to dataframe for use in ecospat ###################environmental################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present") # Upload basic rasters for First interval (Present) Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$') Prescomp.files maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif") maxsalinitypres maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif") maxtemppres rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif") rangesalinitypres rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif") rangetemppres # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas) Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres)) wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') #Upload basic rasters for Second interval (Holocene) ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam") Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$') #Holo.files # scaling holocene layers maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif") maxsalinityholo2 <- maxsalinityholo/100 maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif") maxtempholo2 <- maxtempholo/100 rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif") rangesalinityholo2 <- rangesalinityholo/100 rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif") rangetempholo2 <- rangetempholo/100 Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2) # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') ### Ecospat Analysis of Niche Equivalency and Similary #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Set working directory for saving plots and test results setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/Ecospat_apal_10_2_24") combine.lat = c(ApalPres$lat, ApalHolo$lat) #, ApalLGM$lat) combine.lon = c(ApalPres$long, ApalHolo$long) #, ApalLGM$long) #I need to turn the points into a spatvector ( ptspres <- vect(cbind(ApalPres$long, ApalPres$lat)) ptsholo <- vect(cbind(ApalHolo$long, ApalHolo$lat)) # get random sample - background point locations - rasterobject <- rast(Prescomp.rasters) bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) rasterobject <- rast(Holo.rasters) bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) # get environmental data from occ points extractpres = na.omit(cbind(ApalPres[2:3], raster::extract(Prescomp.rasters, ApalPres[2:3]), rep(1, nrow(ApalPres)))) extractholo = na.omit(cbind(ApalHolo[2:3], raster::extract(Holo.rasters, ApalHolo[2:3]), rep(1, nrow(ApalHolo)))) #change name of last column to occ (represents presence = 1) colnames(extractpres)[ncol(extractpres)] = 'occ' colnames(extractholo)[ncol(extractholo)] = 'occ' # get environmental data from bg points #cant do stack-- need to do all??? extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres))))) extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo))))) # change name of last column to occ (represents absence = 0) colnames(extbgpres)[ncol(extbgpres)] = 'occ' colnames(extbgholo)[ncol(extbgholo)] = 'occ' colnames(extbgpres)[1] = 'long' colnames(extbgpres)[2] = 'lat' colnames(extbgholo)[1] = 'long' colnames(extbgholo)[2] = 'lat' # merge data from occ and bg datPresApal = rbind(extractpres, extbgpres) datHoloApal = rbind(extractholo, extbgholo) #principle component analysis pca.env <- dudi.pca( rbind(datHoloApal, datPresApal)[,3:6], scannf=FALSE, nf=2 ) # look at variable contribution ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) # Save pdf of PCA plot pdf('ecospat_Apalpresholo.pdf') ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) scores.globclim<-pca.env$li # PCA scores for the whole study area scores.globclim<-pca.env$li # PCA scores for the whole study area (all points) scores.sppres <- suprow(pca.env, extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution scores.spholo <- suprow(pca.env, extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution scores.climholo <- suprow(pca.env,datHoloApal[,3:6])$li # PCA scores for the first interval study area scores.climpres <- suprow(pca.env,datPresApal[,3:6])$li # PCA scores for the second interval study area # density distribtuion for first interval grid.Apalholo <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climholo, sp = scores.spholo, R = 100, th.sp = 0 ) # density distribution for second interval grid.Apalpres <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climpres, sp = scores.sppres, R = 100, th.sp = 0 ) histpres <- hist(grid.Apalpres$glob1) histholo <- hist(grid.Apalholo$glob1) plot( histpres, col=rgb(0,0,1,1/4)) # first histogram plot( histholo, col=rgb(1,0,0,1/4), add=T) # Schoener's D metric and I metric D.overlap <- ecospat.niche.overlap (grid.Apalholo, grid.Apalpres, cor=T) D.overlap write.csv(D.overlap,"presApal_holoApal_I_D.csv",row.names=FALSE) ## Niche Equivalency Test ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D ## p.I: pvalue of the test on I ## Test for greater or lower equivalency eq.testgr <- ecospat.niche.equivalency.test(grid.Apalholo, grid.Apalpres, rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs eq.testlw <- ecospat.niche.equivalency.test(grid.Apalholo, grid.Apalpres, rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs nichdynindex <- ecospat.niche.dyn.index(grid.Apalholo, grid.Apalpres, intersection = NA) #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w nichdynindex # write p values of equivalency test p_EQ_DI = cbind("p.D_GR"=eq.testgr$p.D,"p.I_GR"=eq.testgr$p.I,"p.D_LW"=eq.testlw$p.D,"p.I_LW"=eq.testlw$p.D) write.csv(p_EQ_DI,"glycim_EQ_TestApalpresholo.csv",row.names=FALSE) p_EQ_DI ## Test for greater (niche conservatism) or lower (niche divergence) similarity sim.testgr <- ecospat.niche.similarity.test(grid.Apalholo, grid.Apalpres, rep=1000, overlap.alternative = "higher", rand.type=2,ncores=4) sim.testlw <- ecospat.niche.similarity.test(grid.Apalholo, grid.Apalpres, rep=1000, overlap.alternative = "lower", rand.type=2,ncores=4) # write p values of similarity test p_SIM_DI = cbind("p.D_GR"=sim.testgr$p.D,"p.I_GR"=sim.testgr$p.I,"p.D_LW"=sim.testlw$p.D,"p.I_LW"=sim.testlw$p.D) write.csv(p_SIM_DI,"presholoApal_SIM_Test.csv",row.names=FALSE) p_SIM_DI #p.D_GR p.I_GR p.D_LW p.I_LW #[1,] 0.04595405 0.03996004 0.9480519 0.9480519 #dev.off() # Plot test distributions ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency") ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency") ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity") ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity") ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals ecospat.plot.niche (grid.Apalpres, title='Apalpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE) ecospat.plot.niche (grid.Apalholo, title='Apalholo', name.axis1='PC1', name.axis2='PC2') ## Plot niche overlap ecospat.plot.niche.dyn(z1 = grid.Apalholo, z2 = grid.Apalpres, quant=0.25, interest=2, title= "A. pal Holo to Pres", name.axis1="PC1", name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = 'black', colZ1 = "#313695", colZ2 = "#abd9d9", transparency = 40) ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white') rstudioapi::savePlotAsImage("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/Ecospat_apal_10_2_24/figures/A. pal Holo to Pres.png",width=1000,height=750) #___________________________ # Save pdf of plots ### Look at niche expantion, stability, and unfilling # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Apalholo, z2 = grid.Apalpres, intersection=NA) dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Apalholo, z2 = grid.Apalpres, intersection=0) # write csv of niche dynmaic percentages Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w, "Overlapping"=dynam_overlap$dynamic.index.w) Niche_Ex_St write.csv(Niche_Ex_St,"GlycymDynamicsApalholopres.csv",row.names=TRUE) ########################## CNAT ECOSPAT ############# #need cleaned occurrence data and environmental data rasters ## Read in shapefiles of cleaned occ data ## ########################occurrence data###################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/Master_files_cnat_add_mask") CnatPres = read.csv("masterpres.occs.csv") CnatPres = as(CnatPres,'data.frame') CnatHolo = read.csv("masterholo.occs.csv") CnatHolo = as(CnatHolo,'data.frame') #convert to dataframe for use in ecospat ##################environmental################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present") # Upload basic rasters for First interval (Present) Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$') Prescomp.files maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif") maxsalinitypres maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif") maxtemppres rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif") rangesalinitypres rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif") rangetemppres # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas) Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres)) wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') #Upload basic rasters for Second interval (Holocene) ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam") Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$') #Holo.files maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif") maxsalinityholo2 <- maxsalinityholo/100 maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif") maxtempholo2 <- maxtempholo/100 rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif") rangesalinityholo2 <- rangesalinityholo/100 rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif") rangetempholo2 <- rangetempholo/100 Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2) # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') ### Ecospat Analysis of Niche Equivalency and Similary #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Set working directory for saving plots and test results setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/Ecospat_cnat_10_2_24") combine.lat = c(CnatPres$lat, CnatHolo$lat) #, CnatLGM$lat) combine.lon = c(CnatPres$long, CnatHolo$long) #, CnatLGM$long) #combine.lat = c(glyMaa$coords.x1, glyDan$coords.x1) #combine.lon = c(glyMaa$coords.x2, glyDan$coords.x2) #extt=extent(c(min(combine.lon)-5, max(combine.lon)+5, min(combine.lat)-5, max(combine.lat)+5)) # useful if mapping pts and env #I need to turn the points into a spatvector ( library(terra) ptspres <- vect(cbind(CnatPres$long, CnatPres$lat)) ptsholo <- vect(cbind(CnatHolo$long, CnatHolo$lat)) # get random sample - background point locations - rasterobject <- rast(Prescomp.rasters) bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) rasterobject <- rast(Holo.rasters) bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) # get environmental data from occ points extractpres = na.omit(cbind(CnatPres[2:3], raster::extract(Prescomp.rasters, CnatPres[2:3]), rep(1, nrow(CnatPres)))) extractholo = na.omit(cbind(CnatHolo[2:3], raster::extract(Holo.rasters, CnatHolo[2:3]), rep(1, nrow(CnatHolo)))) #change name of last column to occ (represents presence = 1) colnames(extractpres)[ncol(extractpres)] = 'occ' colnames(extractholo)[ncol(extractholo)] = 'occ' # get environmental data from bg points #cant do stack-- need to do all??? extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres))))) extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo))))) # change name of last column to occ (represents absence = 0) colnames(extbgpres)[ncol(extbgpres)] = 'occ' colnames(extbgholo)[ncol(extbgholo)] = 'occ' colnames(extbgpres)[1] = 'long' colnames(extbgpres)[2] = 'lat' colnames(extbgholo)[1] = 'long' colnames(extbgholo)[2] = 'lat' #colnames(extbgholo)[ncol(extbgholo)] = 'occ' # merge data from occ and bg datPresCnat = rbind(extractpres, extbgpres) datHoloCnat = rbind(extractholo, extbgholo) #principle component analysis pca.env <- dudi.pca( rbind(datHoloCnat, datPresCnat)[,3:6], scannf=FALSE, nf=2 ) # look at variable contribution ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) #dev.off() # Save pdf of PCA plot pdf('ecospat_Cnatpresholo.pdf') ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) scores.globclim<-pca.env$li # PCA scores for the whole study area scores.globclim<-pca.env$li # PCA scores for the whole study area (all points) scores.sppres <- suprow(pca.env, extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution scores.spholo <- suprow(pca.env, extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution scores.climholo <- suprow(pca.env,datHoloCnat[,3:6])$li # PCA scores for the first interval study area scores.climpres <- suprow(pca.env,datPresCnat[,3:6])$li # PCA scores for the second interval study area #dev.on() # density distribtuion for first interval grid.Cnatholo <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climholo, sp = scores.spholo, R = 100, th.sp = 0 ) # density distribution for second interval grid.Cnatpres <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climpres, sp = scores.sppres, R = 100, th.sp = 0 ) #dev.off() histpres <- hist(grid.Cnatpres$glob1) histholo <- hist(grid.Cnatholo$glob1) plot( histpres, col=rgb(0,0,1,1/4)) # first histogram plot( histholo, col=rgb(1,0,0,1/4), add=T) # Schoener's D metric and I metric D.overlap <- ecospat.niche.overlap (grid.Cnatholo, grid.Cnatpres, cor=T) D.overlap write.csv(D.overlap,"presCnat_holoCnat_I_D.csv",row.names=FALSE) ## Niche Equivalency Test ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D ## p.I: pvalue of the test on I ## Test for greater or lower equivalency eq.testgr <- ecospat.niche.equivalency.test(grid.Cnatholo, grid.Cnatpres, rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs eq.testlw <- ecospat.niche.equivalency.test(grid.Cnatholo, grid.Cnatpres, rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs nichdynindex <- ecospat.niche.dyn.index(grid.Cnatholo, grid.Cnatpres, intersection = NA) #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w nichdynindex # write p values of equivalency test p_EQ_DI = cbind("p.D_GR"=eq.testgr$p.D,"p.I_GR"=eq.testgr$p.I,"p.D_LW"=eq.testlw$p.D,"p.I_LW"=eq.testlw$p.D) write.csv(p_EQ_DI,"glycim_EQ_TestCnatpresholo.csv",row.names=FALSE) p_EQ_DI ## Test for greater (niche conservatism) or lower (niche divergence) similarity sim.testgr <- ecospat.niche.similarity.test(grid.Cnatholo, grid.Cnatpres, rep=1000, overlap.alternative = "higher", rand.type=2,ncores=4) sim.testlw <- ecospat.niche.similarity.test(grid.Cnatholo, grid.Cnatpres, rep=1000, overlap.alternative = "lower", rand.type=2,ncores=4) # write p values of similarity test p_SIM_DI = cbind("p.D_GR"=sim.testgr$p.D,"p.I_GR"=sim.testgr$p.I,"p.D_LW"=sim.testlw$p.D,"p.I_LW"=sim.testlw$p.D) write.csv(p_SIM_DI,"presholoCnat_SIM_Test.csv",row.names=FALSE) p_SIM_DI #p.D_GR p.I_GR p.D_LW p.I_LW #[1,] 0.04595405 0.03996004 0.9480519 0.9480519 #dev.off() # Plot test distributions ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency") ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency") ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity") ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity") ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals ecospat.plot.niche (grid.Cnatpres, title='Cnatpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE) ecospat.plot.niche (grid.Cnatholo, title='Cnatholo', name.axis1='PC1', name.axis2='PC2') ## Plot niche overlap ecospat.plot.niche.dyn(z1 = grid.Cnatholo, z2 = grid.Cnatpres, quant=0.25, interest=2, title= "C. nat Holo to Pres", name.axis1="PC1", name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = 'black', colZ1 = "#313695", colZ2 = "#abd9d9", transparency = 40) ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white') rstudioapi::savePlotAsImage("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/Ecospat_cnat_10_2_24/figures/C. nat Holo to Pres.pdf",width=1000,height=750) #___________________________ ### Look at niche expantion, stability, and unfilling # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Cnatholo, z2 = grid.Cnatpres, intersection=NA) dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Cnatholo, z2 = grid.Cnatpres, intersection=0) # write csv of niche dynmaic percentages Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w, "Overlapping"=dynam_overlap$dynamic.index.w) Niche_Ex_St write.csv(Niche_Ex_St,"GlycymDynamicsCnatholopres.csv",row.names=TRUE) ############################ PAST ECOSPAT ###################### #### MODEL COMPARISON #### #need cleaned occurrence data and environmental data rasters ## Read in shapefiles of cleaned occ data ## ########################occurrence data###################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/Master_files_past_add_mask") PastPres = read.csv("masterpres.occs.csv") PastPres = as(PastPres,'data.frame') PastHolo = read.csv("masterholo.occs.csv") PastHolo = as(PastHolo,'data.frame') #convert to dataframe for use in ecospat ###################environmental################### setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present") # Upload basic rasters for First interval (Present) Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$') Prescomp.files maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif") maxsalinitypres maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif") maxtemppres rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif") rangesalinitypres rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif") rangetemppres # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas) Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres)) wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # optionally assigning a crs if one isn't provided crs(Prescomp.rasters) <- wgs1984 names(Prescomp.rasters) names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') #Upload basic rasters for Second interval (Holocene) ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam") Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$') #Holo.files maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif") maxsalinityholo2 <- maxsalinityholo/100 maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif") maxtempholo2 <- maxtempholo/100 rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif") rangesalinityholo2 <- rangesalinityholo/100 rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif") rangetempholo2 <- rangetempholo/100 Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2) # optionally assigning a crs if one isn't provided crs(Holo.rasters) <- wgs1984 names(Holo.rasters) names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp') ### Ecospat Analysis of Niche Equivalency and Similary #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Set working directory for saving plots and test results setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/Ecospat_past_10_2_24") combine.lat = c(PastPres$lat, PastHolo$lat) #, PastLGM$lat) combine.lon = c(PastPres$long, PastHolo$long) #, PastLGM$long) #I need to turn the points into a spatvector ( ptspres <- vect(cbind(PastPres$long, PastPres$lat)) ptsholo <- vect(cbind(PastHolo$long, PastHolo$lat)) # get random sample - background point locations - rasterobject <- rast(Prescomp.rasters) bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) rasterobject <- rast(Holo.rasters) bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject, return.type = "points", n = 500) # get environmental data from occ points extractpres = na.omit(cbind(PastPres[2:3], raster::extract(Prescomp.rasters, PastPres[2:3]), rep(1, nrow(PastPres)))) extractholo = na.omit(cbind(PastHolo[2:3], raster::extract(Holo.rasters, PastHolo[2:3]), rep(1, nrow(PastHolo)))) #change name of last column to occ (represents presence = 1) colnames(extractpres)[ncol(extractpres)] = 'occ' colnames(extractholo)[ncol(extractholo)] = 'occ' # get environmental data from bg points #cant do stack-- need to do all??? extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres))))) extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo))))) # change name of last column to occ (represents absence = 0) colnames(extbgpres)[ncol(extbgpres)] = 'occ' colnames(extbgholo)[ncol(extbgholo)] = 'occ' colnames(extbgpres)[1] = 'long' colnames(extbgpres)[2] = 'lat' colnames(extbgholo)[1] = 'long' colnames(extbgholo)[2] = 'lat' #colnames(extbgholo)[ncol(extbgholo)] = 'occ' # merge data from occ and bg datPresPast = rbind(extractpres, extbgpres) datHoloPast = rbind(extractholo, extbgholo) #principle component analysis pca.env <- dudi.pca( rbind(datHoloPast, datPresPast)[,3:6], scannf=FALSE, nf=2 ) # look at variable contribution ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) #dev.off() # Save pdf of PCA plot pdf('ecospat_Pastpresholo.pdf') ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig) scores.globclim<-pca.env$li # PCA scores for the whole study area scores.globclim<-pca.env$li # PCA scores for the whole study area (all points) scores.sppres <- suprow(pca.env, extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution scores.spholo <- suprow(pca.env, extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution scores.climholo <- suprow(pca.env,datHoloPast[,3:6])$li # PCA scores for the first interval study area scores.climpres <- suprow(pca.env,datPresPast[,3:6])$li # PCA scores for the second interval study area #dev.on() # density distribtuion for first interval grid.Pastholo <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climholo, sp = scores.spholo, R = 100, th.sp = 0 ) # density distribution for second interval grid.Pastpres <- ecospat.grid.clim.dyn( glob = scores.globclim, glob1 = scores.climpres, sp = scores.sppres, R = 100, th.sp = 0 ) histpres <- hist(grid.Pastpres$glob1) histholo <- hist(grid.Pastholo$glob1) plot( histpres, col=rgb(0,0,1,1/4)) # first histogram plot( histholo, col=rgb(1,0,0,1/4), add=T) # Schoener's D metric and I metric D.overlap <- ecospat.niche.overlap (grid.Pastholo, grid.Pastpres, cor=T) D.overlap write.csv(D.overlap,"presPast_holoPast_I_D.csv",row.names=FALSE) ## Niche Equivalency Test ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D ## p.I: pvalue of the test on I ## Test for greater or lower equivalency eq.testgr <- ecospat.niche.equivalency.test(grid.Pastholo, grid.Pastpres, rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs eq.testlw <- ecospat.niche.equivalency.test(grid.Pastholo, grid.Pastpres, rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs #this is if they are the same niche I think nichdynindex <- ecospat.niche.dyn.index(grid.Pastholo, grid.Pastpres, intersection = NA) #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w nichdynindex # write p values of equivalency test p_EQ_DI = cbind("p.D_GR"=eq.testgr$p.D,"p.I_GR"=eq.testgr$p.I,"p.D_LW"=eq.testlw$p.D,"p.I_LW"=eq.testlw$p.D) write.csv(p_EQ_DI,"glycim_EQ_TestPastpresholo.csv",row.names=FALSE) p_EQ_DI ## Test for greater (niche conservatism) or lower (niche divergence) similarity sim.testgr <- ecospat.niche.similarity.test(grid.Pastholo, grid.Pastpres, rep=1000, overlap.alternative = "higher", rand.type=2,ncores=4) sim.testlw <- ecospat.niche.similarity.test(grid.Pastholo, grid.Pastpres, rep=1000, overlap.alternative = "lower", rand.type=2,ncores=4) # write p values of similarity test p_SIM_DI = cbind("p.D_GR"=sim.testgr$p.D,"p.I_GR"=sim.testgr$p.I,"p.D_LW"=sim.testlw$p.D,"p.I_LW"=sim.testlw$p.D) write.csv(p_SIM_DI,"presholoPast_SIM_Test.csv",row.names=FALSE) p_SIM_DI #p.D_GR p.I_GR p.D_LW p.I_LW # Plot test distributions ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency") ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency") ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity") ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity") ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals ecospat.plot.niche (grid.Pastpres, title='Pastpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE) ecospat.plot.niche (grid.Pastholo, title='Pastholo', name.axis1='PC1', name.axis2='PC2') ## Plot niche overlap ecospat.plot.niche.dyn(z1 = grid.Pastholo, z2 = grid.Pastpres, quant=0.25, interest=2, title= "P. ast Holo to Pres", name.axis1="PC1", name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = '#74add1', colZ1 = "#313695", colZ2 = "#abd9d9", transparency = 40) ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white') rstudioapi::savePlotAsImage("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/Ecospat_past_10_2_24/figures/P. ast Holo to Pres.png",width=1000,height=750) #___________________________ ### Look at niche expantion, stability, and unfilling # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Pastholo, z2 = grid.Pastpres, intersection=NA) dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Pastholo, z2 = grid.Pastpres, intersection=0) # write csv of niche dynmaic percentages Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w, "Overlapping"=dynam_overlap$dynamic.index.w) Niche_Ex_St write.csv(Niche_Ex_St,"GlycymDynamicsPastholopres.csv",row.names=TRUE) #################&&&&&&&&&&&&&&&&&&&&&############ ENM ################# ########################################## ACER PRES ENM ################################ ##### Loading in the Data -------------------------------------------------------------------------------------------------- # master background file; setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/Master_files_acerv_add_mask") All.back <- read.csv('masterall.background.csv', header=T) #filter to time period background pres.back <- All.back %>% filter(time.bin== c('Present')) Holo.back <- All.back %>% filter(time.bin=='Holocene') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(pres.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/models_Acerv_pres_add_mask/pres_back.csv", row.names=FALSE) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files pres.occs <- read.csv('masterpres.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(pres.back), colnames(pres.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # create null df - empty dataframe that null data will go into later pres.null <- create.null.df('Acervicornis', 'Present', 'Caribbean') #function(taxon.name, time.bin, extent){ # create summary df pres.summary <- create.summary.df('Acervicornis', 'Present', 'Caribbean') create.folders.for.maxent(pres.summary) # run null model Acerv.null <- null.aic(null.df = pres.null, occs = pres.occs, background = pres.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Acerv.optim <- optimize.maxent.likelihood(pres.summary, # name of the output summary file occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Acerv.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Acerv.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/models_Acerv_pres_add_mask/Acerv_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### LQP 0.10 # create eval object # default f.class = LQP, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP Acervpres.eval_LQP0.10 <- create.eval.df('Acerv', 'Present', 'Caribbean', beta.values = 0.10, f.class = 'LQP') Acervpres.eval_LQP0.10 # create folders for Baculites in the Cenomanian create.folders.for.maxent(Acervpres.eval_LQP0.10) # run eval object Acervpres.eval_LQP0.10 <- maxent.crossval.error(eval.df = Acervpres.eval_LQP0.10, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0.025, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) # response curves good ### LQP 0.25 Acervpres.eval_LQP0.25 <- create.eval.df('Acerv', 'Present', 'Caribbean', beta.values = 0.25, f.class = 'LQP') Acervpres.eval_LQP0.25 # create folders create.folders.for.maxent(Acervpres.eval_LQP0.25) # run eval object Acervpres.eval_LQP0.25 <- maxent.crossval.error(eval.df = Acervpres.eval_LQP0.25, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #response curves similar- not as good, donmt go down as far ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps Acervpres.eval_LQP0.10$summary # pROC 1.672626 Acervpres.eval_LQP0.25$summary # pROC 1.662092 # going with LQP 0.10 ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### LQP 0.10 # evaluate model to calculate weighted suitibility and stdev Acervpres.means_LQP0.10 <- maxent.eval(eval = Acervpres.eval_LQP0.10) thresh = min(Acervpres.means_LQP0.10$occ$w.mean) print(thresh) # this is the threshold # save this value!! ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(Acervpres.eval_LQP0.10) ggsave('Acervpres.eval_LQP0.10.pdf') ##### Projecting Model to all extents -------------------------------------------------------------------------------------- Acervpres.thresh_LQP0.10 <- 0.001166019 Acervpres.everything_LQP0.10 <- maxent.everything(eval = Acervpres.eval_LQP0.10, thresh = Acervpres.thresh_LQP0.10, means = Acervpres.means_LQP0.10, everything= All.back, predic = 6:9) View(Acervpres.everything_LQP0.10) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Acerv.tolerance <- c(1,1, # max salinity 1,1, # max temp 1,1, # range salinity 1,1) # range temp ################### running mess holocene Acerv.mess.Holo <- informed.mess(ref.extent = pres.back, mess.extent = Holo.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 45 2050 Acerv.mess.452050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 45 2100 Acerv.mess.452100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 85 2050 Acerv.mess.852050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 85 2100 Acerv.mess.852100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) #### Acerv mess all Acerv.mess.all <- informed.mess(ref.extent = pres.back, mess.extent = All.back, coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- pres.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') pres.output <- cbind(pres.output,Acervpres.everything_LQP0.10) pres.output$Name <- 'Acervicornis' pres.with.mess <- cbind(pres.output, Acerv.mess.all) write.csv(pres.output, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/models_Acerv_pres_add_mask/pres.output.all.csv", row.names = F) # null model write.csv(Acerv.null, 'Acerv.null.csv', row.names = F) # model optimized summary object write.csv(Acerv.optim, 'Acerv.optim.summary.csv', row.names = F) # model projected to all extents write.csv(Acervpres.everything_LQP0.10, 'Acerv.LQP0.10.csv', row.names = F) # model mess analysis write.csv(Acerv.mess.Holo, 'Acerv.mess.Holo.csv', row.names = F) write.csv(Acerv.mess.452050, 'Acerv.mess.452050.csv', row.names = F) write.csv(Acerv.mess.452100, 'Acerv.mess.452100.csv', row.names = F) write.csv(Acerv.mess.852050, 'Acerv.mess.852050.csv', row.names = F) write.csv(Acerv.mess.852100, 'Acerv.mess.852100.csv', row.names = F) ########################################## ACER HOLO ENM ################################ # master background file; All.back <- read.csv('masterall.background.csv', header=T) pres.back <- All.back %>% filter(time.bin=='Present') Holo.back <- All.back %>% filter(time.bin=='Holocene') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(Holo.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/models_Acerv_Holo_add_mask/Holo_back.csv", row.names=FALSE) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files Pres.occs <- read.csv('masterpres.occs.csv', header=T) Holo.occs <- read.csv('masterholo.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(Holo.back), colnames(Holo.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # create null df - empty dataframe that null data will go into later Holo.null <- create.null.df('Acervicornis', 'Holocene', 'Caribbean') #function(taxon.name, time.bin, extent) # create summary df Holo.summary <- create.summary.df('Acervicornis', 'Holocene', 'Caribbean') create.folders.for.maxent(Holo.summary) # run null model Acerv.null <- null.aic(null.df = Holo.null, occs = Holo.occs, background = Holo.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Acerv.optim <- optimize.maxent.likelihood(Holo.summary, # name of the output summary file occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Acerv.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Acerv.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/models_Acerv_Holo_add_mask/Acerv_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### Q 0.05 # create eval object # default f.class = Q, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP AcervHolo.eval_Q0.05 <- create.eval.df('Acerv', 'Holocene', 'Caribbean', beta.values = 0.05, f.class = 'Q') AcervHolo.eval_Q0.05 # create folders for Baculites in the Cenomanian create.folders.for.maxent(AcervHolo.eval_Q0.05) # run eval object AcervHolo.eval_Q0.05 <- maxent.crossval.error(eval.df = AcervHolo.eval_Q0.05, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) # WORKED!!!! curves are curves!!!!!!! or at least temp is curves ### Q0.025 AcervHolo.eval_Q0.025 <- create.eval.df('Acerv', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'Q') AcervHolo.eval_Q0.025 # create folders create.folders.for.maxent(AcervHolo.eval_Q0.025) # run eval object AcervHolo.eval_Q0.025 <- maxent.crossval.error(eval.df = AcervHolo.eval_Q0.025, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #curves are not as good going with other #still a good model compared to null and has good curves ### Q 0.10 AcervHolo.eval_LQP0.025 <- create.eval.df('Acerv', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'LQP') AcervHolo.eval_LQP0.025 # create folders create.folders.for.maxent(AcervHolo.eval_LQP0.025) # run eval object AcervHolo.eval_LQP0.025 <- maxent.crossval.error(eval.df = AcervHolo.eval_LQP0.025, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps AcervHolo.eval_Q0.025$summary # 1.823454 AcervHolo.eval_Q0.05$summary # 1.821950 AcervHolo.eval_LQP0.025$summary #1.887521 # neither mode seems to be systematically better basically identical, i will do LQP0.025 ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### Q 0.10 # evaluate model to calculate weighted suitibility and stdev AcervHolo.means_LQP0.025 <- maxent.eval(eval = AcervHolo.eval_LQP0.025) thresh = min(AcervHolo.means_LQP0.025$occ$w.mean) print(thresh) # 0.08827532 # this is the threshold # save this value!! ### Q 0.5 ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(AcervHolo.eval_LQP0.025) ggsave('AcervHolo.eval_LQP0.025.pdf') ### deciding to go with Q 0.10 since it has slightly more realistic response curves ##### Projecting Model to all extents -------------------------------------------------------------------------------------- AcervHolo.thresh_LQP0.025 <- 0.08827532 AcervHolo.everything_LQP0.025 <- maxent.everything(eval = AcervHolo.eval_LQP0.025, thresh = AcervHolo.thresh_LQP0.025, means = AcervHolo.means_LQP0.025, everything= All.back, predic = 6:9) View(AcervHolo.everything_LQP0.025) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Acerv.tolerance <- c(1,1, # max salinity 1,1, # max temp 1,1, # range salinity 1,0) # range temp ################### running mess present Acerv.mess.pres <- informed.mess(ref.extent = Holo.back, mess.extent = pres.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 45 2050 Acerv.mess.452050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 45 2100 Acerv.mess.452100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 85 2050 Acerv.mess.852050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ################### running mess 85 2100 Acerv.mess.852100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) #mess all Acerv.mess.all <- informed.mess(ref.extent = Holo.back, mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Acerv.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- pres.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') pres.output <- cbind(pres.output, AcervHolo.everything_LQP0.025) pres.output$Name <- 'Acervicornis' pres.with.mess <- cbind(pres.output, Acerv.mess.all) write.csv(pres.output, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/models_Acerv_Holo_add_mask/pres.output.all.csv", row.names = F) # null model write.csv(Acerv.null, 'Acerv.null.csv', row.names = F) # model optimized summary object write.csv(Acerv.optim, 'Acerv.optim.summary.csv', row.names = F) # model projected to all extents write.csv(AcervHolo.everything_LQP0.025, 'Acerv.LQP0.025.csv', row.names = F) # model mess analysis write.csv(Acerv.mess.pres, 'Acerv.mess.pres.csv', row.names = F) write.csv(Acerv.mess.452050, 'Acerv.mess.452050.csv', row.names = F) write.csv(Acerv.mess.452100, 'Acerv.mess.452100.csv', row.names = F) write.csv(Acerv.mess.852050, 'Acerv.mess.852050.csv', row.names = F) write.csv(Acerv.mess.852100, 'Acerv.mess.852100.csv', row.names = F) ######################################### APAL PRES ENM ############ ##### Loading in the Data -------------------------------------------------------------------------------------------------- # master background file; setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/Master_files_apal_add_mask") All.back <- read.csv('masterall.background.csv', header=T) #filter to time period background pres.back <- All.back %>% filter(time.bin== c('Present')) Holo.back <- All.back %>% filter(time.bin=='Holocene') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(pres.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/models_Apal_pres_add_mask/pres_back.csv", row.names=FALSE) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files pres.occs <- read.csv('masterpres.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(pres.back), colnames(pres.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # working directory for models setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/models_Apal_pres_add_mask") # create null df - empty dataframe that null data will go into later pres.null <- create.null.df('Apalmata', 'Present', 'Caribbean') #function(taxon.name, time.bin, extent){ # create summary df pres.summary <- create.summary.df('Apalmata', 'Present', 'Caribbean') create.folders.for.maxent(pres.summary) # run null model Apal.null <- null.aic(null.df = pres.null, occs = pres.occs, background = pres.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Apal.optim <- optimize.maxent.likelihood(pres.summary, # name of the output summary file occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Apal.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Apal.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/models_Apal_pres_add_mask/Apal_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### LQP .05 # create eval object # default f.class = LQP, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP Apalpres.eval_LQP0.05 <- create.eval.df('Apal', 'Present', 'Caribbean', beta.values = 0.05, f.class = 'LQP') Apalpres.eval_LQP0.05 # create folders for Baculites in the Cenomanian create.folders.for.maxent(Apalpres.eval_LQP0.05) # run eval object Apalpres.eval_LQP0.05 <- maxent.crossval.error(eval.df = Apalpres.eval_LQP0.05, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) # response curves good max temp has a true hump that goes all the way down, only range tamp increases to lower ranges but that makes perfect sense right? ### LQP 0.025 Apalpres.eval_LQP0.025 <- create.eval.df('Apal', 'Present', 'Caribbean', beta.values = 0.025, f.class = 'LQP') Apalpres.eval_LQP0.025 # create folders create.folders.for.maxent(Apalpres.eval_LQP0.025) # run eval object Apalpres.eval_LQP0.025 <- maxent.crossval.error(eval.df = Apalpres.eval_LQP0.025, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #response curves similar- only thing is max temp starts going back down at high temps but not all the way. ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps Apalpres.eval_LQP0.05$summary # pROC 1.690565 Apalpres.eval_LQP0.025$summary # pROC 1.685479 # going with LQP 0.025 since slightly higher pROC?? im not sure which is better? ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### LQP 0.025 # evaluate model to calculate weighted suitibility and stdev Apalpres.means_LQP0.025 <- maxent.eval(eval = Apalpres.eval_LQP0.025) thresh = min(Apalpres.means_LQP0.025$occ$w.mean) print(thresh) # 0.2613874 # this is the threshold # save this value!! ### Q 0.5 ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(Apalpres.eval_LQP0.025) ggsave('Apalpres.eval_LQP0.025.pdf') ### deciding to go with Q 2.0 since it has slightly more realistic response curves ##### Projecting Model to all extents -------------------------------------------------------------------------------------- Apalpres.thresh_LQP0.025 <- 0.2613874 Apalpres.everything_LQP0.025 <- maxent.everything(eval = Apalpres.eval_LQP0.025, thresh = Apalpres.thresh_LQP0.025, means = Apalpres.means_LQP0.025, everything= All.back, predic = 6:9) View(Apalpres.everything_LQP0.025) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Apal.tolerance <- c(1,1, # max salinity 1,1, # max temp 1,1, # range salinity 1,1) # range temp ################## running mess holocene Apal.mess.Holo <- informed.mess(ref.extent = pres.back, mess.extent = Holo.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 45 2050 Apal.mess.452050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 45 2100 Apal.mess.452100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 85 2050 Apal.mess.852050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 85 2100 Apal.mess.852100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) #### apal mess all Apal.mess.all <- informed.mess(ref.extent = pres.back, mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- apal.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') apal.output <- cbind(apal.output, Apalpres.everything_LQP0.025) apal.output$Name <- 'Apalmata' apal.with.mess <- cbind(apal.output, Apal.mess.all) write.csv(apal.with.mess, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/models_Apal_pres_add_mask/apal.output.all.csv", row.names = F) # null model write.csv(Apal.null, 'Apal.null.csv', row.names = F) # model optimized summary object write.csv(Apal.optim, 'Apal.optim.summary.csv', row.names = F) # model projected to all extents write.csv(Apalpres.everything_LQP0.025, 'Apal.LQP0.025.csv', row.names = F) # model mess analysis write.csv(Apal.mess.Holo, 'Apal.mess.Holo.csv', row.names = F) write.csv(Apal.mess.452050, 'Apal.mess.452050.csv', row.names = F) write.csv(Apal.mess.452100, 'Apal.mess.452100.csv', row.names = F) write.csv(Apal.mess.852050, 'Apal.mess.852050.csv', row.names = F) write.csv(Apal.mess.852100, 'Apal.mess.852100.csv', row.names = F) ####################################### APAL HOLO ENM ############# ##### Loading in the Data -------------------------------------------------------------------------------------------------- # master background file; setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/Master_files_apal_add_mask") All.back <- read.csv('masterall.background.csv', header=T) #filter to time period background LGM.back <- All.back %>% filter(time.bin== c('LGM')) pres.back <- All.back %>% filter(time.bin=='Present') Holo.back <- All.back %>% filter(time.bin=='Holocene') Holo2.back <- All.back %>% filter(time.bin=='Holocene2') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(Holo.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/models_Apal_Holo_add_mask/Holo_back.csv", row.names=FALSE) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files Pres.occs <- read.csv('masterpres.occs.csv', header=T) Holo.occs <- read.csv('masterholo.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(Holo.back), colnames(Holo.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # create null df - empty dataframe that null data will go into later Holo.null <- create.null.df('Apalmata', 'Holocene', 'Caribbean') #function(taxon.name, time.bin, extent) # create summary df Holo.summary <- create.summary.df('Apalmata', 'Holocene', 'Caribbean') create.folders.for.maxent(Holo.summary) # run null model Apal.null <- null.aic(null.df = Holo.null, occs = Holo.occs, background = Holo.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Apal.optim <- optimize.maxent.likelihood(Holo.summary, # name of the output summary file occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Apal.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Apal.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/models_Apal_Holo_add_mask/Apal_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### LQP 0.025 # create eval object # default f.class = Q, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP ApalHolo.eval_LQP0.025 <- create.eval.df('Apal', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'LQP') ApalHolo.eval_LQP0.025 # create folders for Baculites in the Cenomanian create.folders.for.maxent(ApalHolo.eval_LQP0.025) # run eval object ApalHolo.eval_LQP0.025 <- maxent.crossval.error(eval.df = ApalHolo.eval_LQP0.025, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) # WORKED!!!! curves are curves!!!!!!! ### LQP 0.10 ApalHolo.eval_LQP0.10 <- create.eval.df('Apal', 'Holocene', 'Caribbean', beta.values = 0.10, f.class = 'LQP') ApalHolo.eval_LQP0.10 # create folders create.folders.for.maxent(ApalHolo.eval_LQP0.10) # run eval object ApalHolo.eval_LQP0.10 <- maxent.crossval.error(eval.df = ApalHolo.eval_LQP0.10, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #curves are not as good going with other ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps ApalHolo.eval_LQP0.10$summary # 1.958130 ApalHolo.eval_LQP0.025$summary #1.971328 # neither mode seems to be systematically better basically identical, i will do LQP0.025 ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### Q 0.10 # evaluate model to calculate weighted suitibility and stdev ApalHolo.means_LQP0.025 <- maxent.eval(eval = ApalHolo.eval_LQP0.025) thresh = min(ApalHolo.means_LQP0.025$occ$w.mean) print(thresh) # 0.01926973 # this is the threshold # save this value!! ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(ApalHolo.eval_LQP0.025) ggsave('ApalHolo.eval_LQP0.025.pdf') ### deciding to go with Q 0.10 since it has slightly more realistic response curves ##### Projecting Model to all extents -------------------------------------------------------------------------------------- ApalHolo.thresh_LQP0.025 <- 0.01926973 ApalHolo.everything_LQP0.025 <- maxent.everything(eval = ApalHolo.eval_LQP0.025, thresh = ApalHolo.thresh_LQP0.025, means = ApalHolo.means_LQP0.025, everything= All.back, predic = 6:9) View(ApalHolo.everything_LQP0.025) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Apal.tolerance <- c(1,1, # max salinity 1,1, # max temp 0,1, # range salinity 1,1) # range temp ################### running mess present Apal.mess.pres <- informed.mess(ref.extent = Holo.back, mess.extent = pres.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 45 2050 Apal.mess.452050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 45 2100 Apal.mess.452100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 85 2050 Apal.mess.852050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ################### running mess 85 2100 Apal.mess.852100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) #### apal mess all Apal.mess.all <- informed.mess(ref.extent = Holo.back, mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Apal.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- apal.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') apal.output <- cbind(apal.output, ApalHolo.everything_LQP0.025) apal.output$Name <- 'Apalmata' apal.with.mess <- cbind(apal.output, Apal.mess.all) write.csv(apal.output, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/models_Apal_holo_add_mask/apal.output.all.csv", row.names = F) # null model write.csv(Apal.null, 'Apal.null.csv', row.names = F) # model optimized summary object write.csv(Apal.optim, 'Apal.optim.summary.csv', row.names = F) # model projected to all extents write.csv(ApalHolo.everything_LQP0.025, 'Apal.LQP0.025.csv', row.names = F) # model mess analysis write.csv(Apal.mess.pres, 'Apal.mess.pres.csv', row.names = F) write.csv(Apal.mess.452050, 'Apal.mess.452050.csv', row.names = F) write.csv(Apal.mess.452100, 'Apal.mess.452100.csv', row.names = F) write.csv(Apal.mess.852050, 'Apal.mess.852050.csv', row.names = F) write.csv(Apal.mess.852100, 'Apal.mess.852100.csv', row.names = F) ############################## CNAT PRES ENM ########### ##### Loading in the Data -------------------------------------------------------------------------------------------------- # master background file; setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/Master_files_cnat_add_mask") All.back <- read.csv('masterall.background.csv', header=T) #filter to time period background pres.back <- All.back %>% filter(time.bin== c('Present')) Holo.back <- All.back %>% filter(time.bin=='Holocene') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(pres.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/models_Cnat_pres_add_mask/pres_back.csv", row.names=FALSE) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files pres.occs <- read.csv('masterpres.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(pres.back), colnames(pres.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # working directory for models setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/models_Cnat_pres_add_mask") # create null df - empty dataframe that null data will go into later pres.null <- create.null.df('Cnatans', 'Present', 'Caribbean') #function(taxon.name, time.bin, extent){ # create summary df pres.summary <- create.summary.df('Cnatans', 'Present', 'Caribbean') create.folders.for.maxent(pres.summary) # run null model Cnat.null <- null.aic(null.df = pres.null, occs = pres.occs, background = pres.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Cnat.optim <- optimize.maxent.likelihood(pres.summary, # name of the output summary file occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Cnat.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Cnat.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/models_Cnat_pres_add_mask/Cnat_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### LQP 0.1 # create eval object # default f.class = LQP, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP Cnatpres.eval_LQP0.1 <- create.eval.df('Cnat', 'Present', 'Caribbean', beta.values = 0.1, f.class = 'LQP') Cnatpres.eval_LQP0.1 # create folders for Baculites in the Cenomanian create.folders.for.maxent(Cnatpres.eval_LQP0.1) # run eval object Cnatpres.eval_LQP0.1 <- maxent.crossval.error(eval.df = Cnatpres.eval_LQP0.1, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) # response curves good, not all are complete humps ### LQP 0.025 Cnatpres.eval_LQP0.025 <- create.eval.df('Cnat', 'Present', 'Caribbean', beta.values = 0.025, f.class = 'LQP') Cnatpres.eval_LQP0.025 # create folders create.folders.for.maxent(Cnatpres.eval_LQP0.025) # run eval object Cnatpres.eval_LQP0.025 <- maxent.crossval.error(eval.df = Cnatpres.eval_LQP0.025, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #response curves better but quite similar, more clear humps ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps Cnatpres.eval_LQP0.1$summary # pROC 1.477102 Cnatpres.eval_LQP0.025$summary # pROC 1.446531 # going with LQP0.025 because response curves are better ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### LQP 0.025 # evaluate model to calculate weighted suitibility and stdev Cnatpres.means_LQP0.025 <- maxent.eval(eval = Cnatpres.eval_LQP0.025) thresh = min(Cnatpres.means_LQP0.025$occ$w.mean) print(thresh) # 0.03913349 # this is the threshold # save this value!! ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(Cnatpres.eval_LQP0.025) ggsave('Cnatpres.eval_LQP0.025.pdf') ### deciding to go with Q 2.0 since it has slightly more realistic response curves ##### Projecting Model to all extents -------------------------------------------------------------------------------------- Cnatpres.thresh_LQP0.025 <- 0.008100324 Cnatpres.everything_LQP0.025 <- maxent.everything(eval = Cnatpres.eval_LQP0.025, thresh = Cnatpres.thresh_LQP0.025, means = Cnatpres.means_LQP0.025, everything= All.back, predic = 6:9) View(Cnatpres.everything_LQP0.025) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Cnat.tolerance <- c(1,1, # max salinity 1,1, # max temp 1,1, # range salinity 1,1) # range temp ################### running mess holocene Cnat.mess.Holo <- informed.mess(ref.extent = pres.back, mess.extent = Holo.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 45 2050 Cnat.mess.452050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 45 2100 Cnat.mess.452100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 85 2050 Cnat.mess.852050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 85 2100 Cnat.mess.852100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) #### apal mess all Cnat.mess.all <- informed.mess(ref.extent = pres.back, mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- pres.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') pres.output <- cbind(pres.output, Cnatpres.everything_LQP0.025) pres.output$Name <- 'Cnatans' pres.with.mess <- cbind(pres.output, Cnat.mess.all) write.csv(pres.output, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/models_Cnat_pres_add_mask/pres.output.all.csv", row.names = F) #Holo______________________ pres.output.Holocene <- pres.with.mess %>% filter(time.bin == 'Holocene') write.csv(pres.output.Holocene, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/models_Cnat_pres_add_mask/pres.output.Holocene.csv", row.names = F) # null model write.csv(Cnat.null, 'Cnat.null.csv', row.names = F) # model optimized summary object write.csv(Cnat.optim, 'Cnat.optim.summary.csv', row.names = F) # model projected to all extents write.csv(Cnatpres.everything_LQP0.025, 'Cnat.LQP0.025.csv', row.names = F) # model mess analysis write.csv(Cnat.mess.Holo, 'Cnat.mess.Holo.csv', row.names = F) write.csv(Cnat.mess.452050, 'Cnat.mess.452050.csv', row.names = F) write.csv(Cnat.mess.452100, 'Cnat.mess.452100.csv', row.names = F) write.csv(Cnat.mess.852050, 'Cnat.mess.852050.csv', row.names = F) write.csv(Cnat.mess.852100, 'Cnat.mess.852100.csv', row.names = F) ############################## CNAT HOLO ENM ########### ##### Loading in the Data -------------------------------------------------------------------------------------------------- # master background file; setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/Master_files_cnat_add_mask") All.back <- read.csv('masterall.background.csv', header=T) #filter to time period background #LGM.back <- All.back %>% filter(time.bin== c('LGM')) pres.back <- All.back %>% filter(time.bin=='Present') Holo.back <- All.back %>% filter(time.bin=='Holocene') Holo2.back <- All.back %>% filter(time.bin=='Holocene2') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(Holo.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/models_Cnat_Holo_add_mask/Holo_back.csv", row.names=FALSE) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files Pres.occs <- read.csv('masterpres.occs.csv', header=T) Holo.occs <- read.csv('masterholo.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(Holo.back), colnames(Holo.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # sourcing the R code to run all the functions setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code") source("KPG_ENM_SVP2021_Functions-copy.R") #this is the fixed functions file # working directory for models setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/models_Cnat_Holo_add_mask") # create null df - empty dataframe that null data will go into later Holo.null <- create.null.df('Cnatans', 'Holocene', 'Caribbean') #function(taxon.name, time.bin, extent) # create summary df Holo.summary <- create.summary.df('Cnatans', 'Holocene', 'Caribbean') create.folders.for.maxent(Holo.summary) # run null model Cnat.null <- null.aic(null.df = Holo.null, occs = Holo.occs, background = Holo.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Cnat.optim <- optimize.maxent.likelihood(Holo.summary, # name of the output summary file occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Cnat.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Cnat.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/models_Cnat_Holo_add_mask/Cnat_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### LQQ 0.025 # create eval object # default f.class = Q, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP CnatHolo.eval_LQP0.025 <- create.eval.df('Cnat', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'LQP') CnatHolo.eval_LQP0.025 # create folders for Baculites in the Cenomanian create.folders.for.maxent(CnatHolo.eval_LQP0.025) # run eval object CnatHolo.eval_LQP0.025 <- maxent.crossval.error(eval.df = CnatHolo.eval_LQP0.025, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) ### LQP 0.05 CnatHolo.eval_LQP0.05 <- create.eval.df('Cnat', 'Holocene', 'Caribbean', beta.values = 0.05, f.class = 'LQP') CnatHolo.eval_LQP0.05 # create folders create.folders.for.maxent(CnatHolo.eval_LQP0.05) # run eval object CnatHolo.eval_LQP0.05 <- maxent.crossval.error(eval.df = CnatHolo.eval_LQP0.05, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #curves are not as good going with other ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps CnatHolo.eval_LQP0.05$summary # 1.820528 CnatHolo.eval_LQP0.025$summary #1.825881 # curves of 0.025 bas not flat line for max sal so i will do LLQPP0.025 ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### LQP 0.10 # evaluate model to calculate weighted suitibility and stdev CnatHolo.means_LQP0.025 <- maxent.eval(eval = CnatHolo.eval_LQP0.025) thresh = min(CnatHolo.means_LQP0.025$occ$w.mean) print(thresh) # 0.005275519 # this is the threshold # save this value!! ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(CnatHolo.eval_LQP0.025) ggsave('CnatHolo.eval_LQP0.025.pdf') ##### Projecting Model to all extents -------------------------------------------------------------------------------------- CnatHolo.thresh_LQP0.025 <- 0.005275519 CnatHolo.everything_LQP0.025 <- maxent.everything(eval = CnatHolo.eval_LQP0.025, thresh = CnatHolo.thresh_LQP0.025, means = CnatHolo.means_LQP0.025, everything= All.back, predic = 6:9) View(CnatHolo.everything_LQP0.025) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Cnat.tolerance <- c(1,1, # max salinity 1,0, # max temp 1,1, # range salinity 1,1) # range temp ################### running mess present Cnat.mess.pres <- informed.mess(ref.extent = Holo.back, mess.extent = pres.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 45 2050 Cnat.mess.452050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 45 2100 Cnat.mess.452100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 85 2050 Cnat.mess.852050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ################### running mess 85 2100 Cnat.mess.852100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) #### cnat mess all Cnat.mess.all <- informed.mess(ref.extent = Holo.back, mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Cnat.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- pres.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') pres.output <- cbind(pres.output, CnatHolo.everything_LQP0.025) pres.output$Name <- 'Cnatans' pres.with.mess <- cbind(pres.output, Cnat.mess.all) write.csv(pres.output, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/models_Cnat_Holo_add_mask/pres.output.all.csv", row.names = F) # null model write.csv(Cnat.null, 'Cnat.null.csv', row.names = F) # model optimized summary object write.csv(Cnat.optim, 'Cnat.optim.summary.csv', row.names = F) # model projected to all extents write.csv(CnatHolo.everything_LQP0.025, 'Cnat.LQP0.025.csv', row.names = F) # model mess analysis write.csv(Cnat.mess.pres, 'Cnat.mess.pres.csv', row.names = F) write.csv(Cnat.mess.452050, 'Cnat.mess.452050.csv', row.names = F) write.csv(Cnat.mess.452100, 'Cnat.mess.452100.csv', row.names = F) write.csv(Cnat.mess.852050, 'Cnat.mess.852050.csv', row.names = F) write.csv(Cnat.mess.852100, 'Cnat.mess.852100.csv', row.names = F) ############################## PAST PRES ENM ########### # master background file; setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/Master_files_past_add_mask") All.back <- read.csv('masterall.background.csv', header=T) #filter to time period background pres.back <- All.back %>% filter(time.bin== c('Present')) Holo.back <- All.back %>% filter(time.bin=='Holocene') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(pres.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/models_Past_pres_add_mask/pres_back.csv", row.names=FALSE) #plotting histogram to see distribution (can do for other background layers if want) hist(pres.back$maxtemp) rstudioapi::savePlotAsImage("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/models_Past_pres_add_mask/figures/presmaxtemphist.png",width=1000,height=750) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files pres.occs <- read.csv('masterpres.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(pres.back), colnames(pres.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # working directory for models setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/models_Past_pres_add_mask") # create null df - empty dataframe that null data will go into later pres.null <- create.null.df('Pasteroides', 'Present', 'Caribbean') #function(taxon.name, time.bin, extent){ # create summary df pres.summary <- create.summary.df('Pasteroides', 'Present', 'Caribbean') create.folders.for.maxent(pres.summary) # run null model Past.null <- null.aic(null.df = pres.null, occs = pres.occs, background = pres.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Past.optim <- optimize.maxent.likelihood(pres.summary, # name of the output summary file occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Past.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Past.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/models_Past_pres_add_mask/Past_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### LQP .05 # create eval object # default f.class = LQP, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP Pastpres.eval_LQP0.25 <- create.eval.df('Past', 'Present', 'Caribbean', beta.values = 0.25, f.class = 'LQP') Pastpres.eval_LQP0.25 # create folders for Baculites in the Cenomanian create.folders.for.maxent(Pastpres.eval_LQP0.25) # run eval object Pastpres.eval_LQP0.25 <- maxent.crossval.error(eval.df = Pastpres.eval_LQP0.25, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) # response curves good. all are humps ### LQP 0.10 Pastpres.eval_LQP0.10 <- create.eval.df('Past', 'Present', 'Caribbean', beta.values = 0.10, f.class = 'LQP') Pastpres.eval_LQP0.10 # create folders create.folders.for.maxent(Pastpres.eval_LQP0.10) # run eval object Pastpres.eval_LQP0.10 <- maxent.crossval.error(eval.df = Pastpres.eval_LQP0.10, occs = pres.occs, # species occurrences background = pres.back, # background for the pres predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0.1, #express as proportion all.background = All.back) #response curves similar- the humps are more complete and cleaner, will go with this one ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps Pastpres.eval_LQP0.25$summary # pROC 1.273761 Pastpres.eval_LQP0.10$summary # pROC 1.302769 # going with LQP 0.10 since slightly higher pROC and its curves are better ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### LQP 0.10 # evaluate model to calculate weighted suitibility and stdev Pastpres.means_LQP0.10 <- maxent.eval(eval = Pastpres.eval_LQP0.10) thresh = min(Pastpres.means_LQP0.10$occ$w.mean) print(thresh) # 8.434429e-05 # this is the threshold # save this value!! ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(Pastpres.eval_LQP0.10) ggsave('Pastpres.eval_LQP0.10.pdf') ### deciding to go with Q 2.0 since it has slightly more realistic response curves ##### Projecting Model to all extents -------------------------------------------------------------------------------------- Pastpres.thresh_LQP0.10 <- 8.928984e-05 Pastpres.everything_LQP0.10 <- maxent.everything(eval = Pastpres.eval_LQP0.10, thresh = Pastpres.thresh_LQP0.10, means = Pastpres.means_LQP0.10, everything= All.back, predic = 6:9) View(Pastpres.everything_LQP0.10) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Past.tolerance <- c(1,1, # max salinity 1,1, # max temp 1,1, # range salinity 1,1) # range temp ################### running mess holocene Past.mess.Holo <- informed.mess(ref.extent = pres.back, mess.extent = Holo.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess LGM Past.mess.LGM <- informed.mess(ref.extent = pres.back, mess.extent = LGM.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 45 2050 Past.mess.452050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 45 2100 Past.mess.452100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 85 2050 Past.mess.852050 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 85 2100 Past.mess.852100 <- informed.mess(ref.extent = pres.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) #mess all for export #### past mess all Past.mess.all <- informed.mess(ref.extent = pres.back, mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- past.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') past.output <- cbind(past.output, Pastpres.everything_LQP0.10) past.output$Name <- 'Pasteroides' past.with.mess <- cbind(past.output, Past.mess.all) write.csv(past.output, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/models_Past_pres_add_mask/past.output.all.csv", row.names = F) # null model write.csv(Past.null, 'Past.null.csv', row.names = F) # model optimized summary object write.csv(Past.optim, 'Past.optim.summary.csv', row.names = F) # model projected to all extents write.csv(Pastpres.everything_LQP0.10, 'Past.LQP0.10.csv', row.names = F) # model mess analysis write.csv(Past.mess.Holo, 'Past.mess.Holo.csv', row.names = F) write.csv(Past.mess.452050, 'Past.mess.452050.csv', row.names = F) write.csv(Past.mess.452100, 'Past.mess.452100.csv', row.names = F) write.csv(Past.mess.852050, 'Past.mess.852050.csv', row.names = F) write.csv(Past.mess.852100, 'Past.mess.852100.csv', row.names = F) ############################## PAST HOLO ENM ########### ##### Loading in the Data -------------------------------------------------------------------------------------------------- # master background file; setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/Master_files_past_add_mask") All.back <- read.csv('masterall.background.csv', header=T) #filter to time period background pres.back <- All.back %>% filter(time.bin=='Present') Holo.back <- All.back %>% filter(time.bin=='Holocene') Holo2.back <- All.back %>% filter(time.bin=='Holocene2') RCP452050.back <- All.back %>% filter(time.bin=='RCP452050') RCP852050.back <- All.back %>% filter(time.bin=='RCP852050') RCP452100.back <- All.back %>% filter(time.bin=='RCP452100') RCP852100.back <- All.back %>% filter(time.bin=='RCP852100') #save background I am using write.csv(Holo.back, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/models_Past_Holo_add_mask/Holo_back.csv", row.names=FALSE) # master occurrence file; all.occ <- read.csv('masterall.occs.csv', header=T) #reading in specific time period files Pres.occs <- read.csv('masterpres.occs.csv', header=T) Holo.occs <- read.csv('masterholo.occs.csv', header=T) #check if background and occurance have same column names identical(colnames(Holo.back), colnames(Holo.occs)) ##### Setting up the models ------------------------------------------------------------------------------------------------ # working directory for models setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/models_Past_Holo_add_mask") # create null df - empty dataframe that null data will go into later Holo.null <- create.null.df('Pasteroides', 'Holocene', 'Caribbean') #function(taxon.name, time.bin, extent) # create summary df Holo.summary <- create.summary.df('Pasteroides', 'Holocene', 'Caribbean') create.folders.for.maxent(Holo.summary) # run null model Past.null <- null.aic(null.df = Holo.null, occs = Holo.occs, background = Holo.back, first.occ.col = 10) ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------ # find most optimized model Past.optim <- optimize.maxent.likelihood(Holo.summary, # name of the output summary file occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10) View(Past.optim) #want model to be at least more than 2 AICc lower than the null #also plan to check pROC and response curves to pick best model write.csv(Past.optim, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/models_Past_Holo_add_mask/Past_optim.csv", row.names=FALSE) ### CHECK RESPONSE CURVES FOR REALISM!!!!! ##### Running Cross-Validation --------------------------------------------------------------------------------------------- ### LQP2.0 # create eval object # default f.class = Q, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP PastHolo.eval_LQP2.0 <- create.eval.df('Past', 'Holocene', 'Caribbean', beta.values = 2.0, f.class = 'LQP') PastHolo.eval_LQP2.0 # create folders for Baculites in the Cenomanian create.folders.for.maxent(PastHolo.eval_LQP2.0) # run eval object PastHolo.eval_LQP2.0 <- maxent.crossval.error(eval.df = PastHolo.eval_LQP2.0, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.) #still a good model compared to null and has good curves ### LQP 0.025 PastHolo.eval_LQP0.025 <- create.eval.df('Past', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'LQP') PastHolo.eval_LQP0.025 # create folders create.folders.for.maxent(PastHolo.eval_LQP0.025) # run eval object PastHolo.eval_LQP0.025 <- maxent.crossval.error(eval.df = PastHolo.eval_LQP0.025, occs = Holo.occs, # species occurrences background = Holo.back, # background for the Holo predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp) first.occ.col = 10, first.test.col = 16, omission.rate = 0, #express as proportion all.background = All.back) ### analyzing summary files ## if running multiple eval models, check for: # 1.) does one model setting have systematically higher omission rates? # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type? # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps PastHolo.eval_LQP2.0$summary PastHolo.eval_LQP0.025$summary # neither mode seems to be systematically better basically identical, i will do LQP0.025 ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------ ### Q 0.10 # evaluate model to calculate weighted suitibility and stdev PastHolo.means_LQP0.025 <- maxent.eval(eval = PastHolo.eval_LQP0.025) thresh = min(PastHolo.means_LQP0.025$occ$w.mean) print(thresh) # 0.1492019 # this is the threshold # save this value!! ### plot model suitability/uncertainty plots # the only input you need is the object made from the maxent.eval function suit.uncert.plot(PastHolo.eval_LQP0.025) ggsave('PastHolo.eval_LQP0.025.pdf') ### deciding to go with Q 0.10 since it has slightly more realistic response curves ##### Projecting Model to all extents -------------------------------------------------------------------------------------- PastHolo.thresh_LQP0.025 <- 0.1492019 PastHolo.everything_LQP0.025 <- maxent.everything(eval = PastHolo.eval_LQP0.025, thresh = PastHolo.thresh_LQP0.025, means = PastHolo.means_LQP0.025, everything= All.back, predic = 6:9) View(PastHolo.everything_LQP0.025) ##### Running Informed MESS Analysis --------------------------------------------------------------------------------------- # defining the tolerance vector # 1 = can extrapolate to non-analog conditions # 0 = cannot extrapolate to non-analog conditions Past.tolerance <- c(1,1, # max salinity 1,0, # max temp 1,1, # range salinity 1,1) # range temp ################### running mess LGM Past.mess.LGM <- informed.mess(ref.extent = Holo.back, mess.extent = LGM.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess present Past.mess.pres <- informed.mess(ref.extent = Holo.back, mess.extent = pres.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 45 2050 Past.mess.452050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 45 2100 Past.mess.452100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 85 2050 Past.mess.852050 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ################### running mess 85 2100 Past.mess.852100 <- informed.mess(ref.extent = Holo.back, mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) #mess all #### apal mess all Past.mess.all <- informed.mess(ref.extent = Holo.back, mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff coord.cols = 2:3, predic = 6:9, tolerance = Past.tolerance) ##### Saving Everything ---------------------------------------------------------------------------------------------------- holo.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin') holo.output <- cbind(holo.output, PastHolo.everything_LQP0.025) holo.output$Name <- 'Pasteroides' pres.with.mess <- cbind(holo.output, Past.mess.all) write.csv(holo.output, "/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/models_Past_Holo_add_mask/holo.output.all.csv", row.names = F) # null model write.csv(Past.null, 'Past.null.csv', row.names = F) # model optimized summary object write.csv(Past.optim, 'Past.optim.summary.csv', row.names = F) # model projected to all extents write.csv(PastHolo.everything_LQP0.025, 'Past.LQP0.025.csv', row.names = F) # model mess analysis write.csv(Past.mess.pres, 'Past.mess.pres.csv', row.names = F) write.csv(Past.mess.452050, 'Past.mess.452050.csv', row.names = F) write.csv(Past.mess.452100, 'Past.mess.452100.csv', row.names = F) write.csv(Past.mess.852050, 'Past.mess.852050.csv', row.names = F) write.csv(Past.mess.852100, 'Past.mess.852100.csv', row.names = F) #######latitude setwd("/Users/clairewilliams/Documents/ENM_thresh_plot/data_illustrator") df <- read.csv('combined_data.csv', stringsAsFactors=FALSE,header = T) df_binned <- df %>% mutate(lat_bin = floor(lat / 5) * 5) %>% # Binning latitude by 5 degrees group_by(Name, lat_bin, time.bin) %>% # Grouping by species, latitude bin, and time bin summarize(count_1 = sum(maxSSS_thresholded == 1), .groups = 'drop') # Counting occurrences of 1 in maxSSS_thresholded # CONVERTING INTO AREA # Constants cell_lat_deg <- 5 / 60 # 5 arc-minutes = 1/12 degree earth_radius_km <- 6371 # Average Earth radius # Approximate latitudinal and longitudinal length of a cell in km cell_height_km <- cell_lat_deg * 111.32 # Roughly constant df_area <- df %>% mutate(lat_bin = floor(lat / 5) * 5) %>% group_by(Name, lat_bin, time.bin) %>% summarize( count_1 = sum(maxSSS_thresholded == 1), .groups = 'drop' ) %>% mutate( lat_rad = lat_bin * pi / 180, # Convert latitude to radians cell_width_km = cell_lat_deg * 111.32 * cos(lat_rad), # Width of each cell cell_area_km2 = cell_height_km * cell_width_km, # Area of one cell total_area_km2 = count_1 * cell_area_km2 # Total area for species/bin )