Integrating fossil data in ecological niche models to improve predictions of future habitat of Caribbean corals.

1. # reading in packages required
2. #### line 34-78
4. # functions
5. #### line 82-1430
7. # making master occurrence files- pres = present holo = holocene
8. #### acer pres master line 1433
9. #### acer holo master line 1631
10. #### apal pres master line 1631
11. #### apal holo master line 1858
12. #### cnat pres master line 1861
13. #### cnat holo master line 2057
14. #### past pres master line 2061
15. #### past holo master line 2280
17. # ecospat line #-#
18. #### acer ecospat line 2293
19. #### apal ecospat line 2550
20. #### cnat ecospat line 2805
21. #### past ecospat line 3079
23. #ecological niche models pres = present holo = holocene
24. #### acer pres ENM line 3390
25. #### acer holo ENM line 3613
26. #### apal pres ENM line 3868
27. #### apal holo ENM line 4105
28. #### cnat pres ENM line 4338
29. #### cnat holo ENM line 4593
30. #### past pres ENM line 4849
31. #### past holo ENM line 5113
33. ############ plotting by area and latitude
34. ### latitude line 5365
36. #################&&&&&&&&&&&&&&&&&&&&&##############Reading in packages
37. # loading the requisite libraries for generating the master.background and master.occurrence files
38. library(raster)
39. library(dplyr)
40. #ecospat packages
41. library(ENMTools)
42. library(dismo)
43. library(raster)
44. library(rgdal) # for spatial data analysis
45. library(dplyr)
46. library(rgeos) # for spatial data analysis
47. library(scales)
48. library(tidyr)
49. library(colorRamps)
50. library(ENMeval) # for a few new tools in ENM/SDM
51. library(ggplot2)
52. library(maptools)
53. library(spThin) # error here
54. library(sdm)
55. library(ecospat)
56. library(ade4)
57. library(dichromat)
58. library(terra)
60. # ENM model packages
61. packs <- c('ade4', 'adehabitatMA', 'adehabitatHR', 'alphahull', 'dismo', 'dplyr', 'ecospat', 'ggplot2', 'jsonlite',
62. 'kuenm', 'matrixStats', 'raster', 'rgdal', 'rgeos', 'rJava',
63. 'sf', 'sp', 'splitstackshape', 'tidyr', 'utils', 'wesanderson', 'PerformanceAnalytics', 'SDMTools')
65. # loading in the packages
66. lapply(packs, library, character.only=T)
68. # giving the package versions
69. packs.df <- as.data.frame(matrix(NA, nrow=length(packs), ncol=2))
70. colnames(packs.df) <- c('pkg.name', 'pkg.version')
71. for(i in 1:length(packs)){
72. packs.df[i,1] <- packs[i]
73. packs.df[i,2] <- as.character(packageVersion(packs[i]))
74. }
75. packs.df
77. library(raster)
78. library(rasterize)
82. ########################### FUNCTIONS #########
83. ##### Prep Parameters ------------------------------------------------------------------------------------------------------
85. # Prep Parameters for Maxent models in R with the dismo package
87. # sourced from https://github.com/shandongfx/workshop_maxent_R/blob/master/code/Appendix2_prepPara.R
88. # https://github.com/shandongfx/workshop_maxent_R/blob/master/code
90. # should be used in concert with Appendix 3 from Feng et al. 2017 (PeerJ Preprints)
91. # https://peerj.com/preprints/3346.pdf (manuscript still unpublished as of August 2020)
92. # Appendix 3
93. # https://github.com/shandongfx/workshop_maxent_R/blob/master/code/Appendix3_maxentParameters_v2.pdf
95. # 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)
97. # A function that implements Maxent parameters using the general R manner
98. # leave "doclamp" as default - later in the code (internal fxns run), "doclamp" is set to FALSE
99. prepPara <- function(userfeatures=NULL, # 41 NULL=autofeature, could be any combination of # c("L", "Q", "H", "P", "T")
100. # MUST be specified as a single string (e.g., "LQ", "LQP", "LQHPT", etc.)
101. responsecurves=TRUE, # 1
102. jackknife=TRUE, # 3
103. outputformat="logistic", # 4
104. outputfiletype="asc", # 5
105. projectionlayers=NULL, # 7
106. randomseed=FALSE, # 10
107. removeduplicates=TRUE, # 16
108. betamultiplier=NULL, # 20, 53-56
109. biasfile=NULL, # 22
110. testsamplesfile=NULL, # 23
111. replicates=1, # 24-25
112. replicatetype="crossvalidate", # 24-25
113. writeplotdata=TRUE, # 37
114. extrapolate=TRUE, # 39
115. doclamp=TRUE, # 42
116. beta_threshold=NULL, # 20, 53-56
117. beta_categorical=NULL, # 20, 53-56
118. beta_lqp=NULL, # 20, 53-56
119. beta_hinge=NULL, # 20, 53-56
120. applythresholdrule=NULL # 60
121. ){
122. #20, 29-33, & 41 features, default is autofeature
123. if(is.null(userfeatures)){
124. args_out <- c("autofeature")
125. } else {
126. args_out <- c("noautofeature")
127. if(grepl("L",userfeatures)) args_out <- c(args_out,"linear") else args_out <- c(args_out,"nolinear")
128. if(grepl("Q",userfeatures)) args_out <- c(args_out,"quadratic") else args_out <- c(args_out,"noquadratic")
129. if(grepl("H",userfeatures)) args_out <- c(args_out,"hinge") else args_out <- c(args_out,"nohinge")
130. if(grepl("P",userfeatures)) args_out <- c(args_out,"product") else args_out <- c(args_out,"noproduct")
131. if(grepl("T",userfeatures)) args_out <- c(args_out,"threshold") else args_out <- c(args_out,"nothreshold")
132. }
134. # 1 - generate response curves for each variable
135. if(responsecurves) args_out <- c(args_out,"responsecurves") else args_out <- c(args_out,"noresponsecurves")
137. # 2
138. #if(picture) args_out <- c(args_out,"pictures") else args_out <- c(args_out,"nopictures")
140. # 3 - apply variable jackknife to see how the model changes if that variable is omitted, then if it's the ONLY variable used
141. if(jackknife) args_out <- c(args_out,"jackknife") else args_out <- c(args_out,"nojackknife")
143. # 4 - output format type. choose from c("logistic", "cumulative", "raw")
144. args_out <- c(args_out,paste0("outputformat=",outputformat))
146. # 5 - output file type. choose from c("asc", "mxe", "grd", "bil")
147. args_out <- c(args_out,paste0("outputfiletype=",outputfiletype))
149. # 7 - pathway to projection layers.
150. # it seems that the projection layers should be the only files in that folder, just like with the MaxEnt .jar file
151. if(!is.null(projectionlayers)) args_out <- c(args_out,paste0("projectionlayers=",projectionlayers))
153. # 10 - will use different random number generators for selecting training vs testing data and background points (if applicable)
154. if(randomseed) args_out <- c(args_out,"randomseed") else args_out <- c(args_out,"norandomseed")
156. # 16 - remove duplicate coordinates that are in the same grid - ONLY for raster data, not SWD
157. if(removeduplicates) args_out <- c(args_out,"removeduplicates") else args_out <- c(args_out,"noremoveduplicates")
159. # 20 & 53-56 - various beta (regularization) multipliers to be applied. default = 1
160. # 20 applies all parameters by this regularization multiplier.
161. # 53-56 can apply uniquely to different feature types
162. # check if negative
163. betas <- c( betamultiplier,beta_threshold,beta_categorical,beta_lqp,beta_hinge)
164. if(! is.null(betas) ){
165. for(i in 1:length(betas)){
166. if(betas[i] <0) stop("betamultiplier has to be positive")
167. }
168. }
169. if ( !is.null(betamultiplier) ){
170. args_out <- c(args_out,paste0("betamultiplier=",betamultiplier))
171. } else {
172. if(!is.null(beta_threshold)) args_out <- c(args_out,paste0("beta_threshold=",beta_threshold))
173. if(!is.null(beta_categorical)) args_out <- c(args_out,paste0("beta_categorical=",beta_categorical))
174. if(!is.null(beta_lqp)) args_out <- c(args_out,paste0("beta_lqp=",beta_lqp))
175. if(!is.null(beta_hinge)) args_out <- c(args_out,paste0("beta_hinge=",beta_hinge))
176. }
178. # 22 - pathway to a bias file for selecting background points - ONLY for raster data, not SWD
179. if(!is.null(biasfile)) args_out <- c(args_out,paste0("biasfile=",biasfile))
181. # 23 - pathway to a test data file - can be in csv format
182. if(!is.null(testsamplesfile)) args_out <- c(args_out,paste0("testsamplesfile=",testsamplesfile))
184. # 24 - replicates = number of replicates to run (integer)
185. # 25 - replicatetype = what type of replicates to run. choose from c('crossvalidate', 'bootstrap', 'subsample')
186. replicates <- as.integer(replicates)
187. if(replicates>1 ){
188. args_out <- c(args_out,
189. paste0("replicates=",replicates),
190. paste0("replicatetype=",replicatetype) )
191. }
193. # 37 - write output files containing the data used to make response curves
194. if(writeplotdata) args_out <- c(args_out,"writeplotdata") else args_out <- c(args_out,"nowriteplotdata")
196. # 39 - allow extrapolation beyond the limits of the training data
197. if(extrapolate) args_out <- c(args_out,"extrapolate") else args_out <- c(args_out,"noextrapolate")
199. # 42 - apply clamping when projecting
200. if(doclamp) args_out <- c(args_out,"doclamp") else args_out <- c(args_out,"nodoclamp")
202. # 60 - threshold your model to binary 1/0
203. # options are: c('Fixed cumulative value 1', 'Fixed cumulative value 5', 'Fixed cumulative value 10', 'Minimum training presence',
204. # '10 percentile training presence', 'Equal training sensitivity and specificity', 'Maximum training sensitivity plus specificity').
205. if(!is.null(applythresholdrule)) args_out <- c(args_out,paste0("applythresholdrule=",applythresholdrule))
207. return(args_out)
208. }
211. # prepPara()
212. # [1] "autofeature" "responsecurves" "jackknife" "outputformat=logistic" "outputfiletype=asc"
213. # [6] "norandomseed" "removeduplicates" "writeplotdata" "extrapolate" "doclamp"
218. ##### Create Null Object, Summary Object, Eval Object, and Nested File Structure -------------------------------------------
220. ### functions that create the null data.frame, summary data.frame, and the eval data.frame
223. ### ARGUMENTS ###
225. ### arguments that exist to keep track of everything and do not change how the functions run
226. ## taxon.name <- name of the entity that will get assigned to any null/summary/eval objects. Does not change how function runs.
227. ## time.bin <- name of the time bin that will get assigned to any null/summary/eval objects. Does not change how function runs.
228. ## extent <- name of the extent that will get assigned to any null/summary/eval objects. Does not change how function runs.
230. ### arguments that set the model hyper-parameters and will change how the functions run
231. ## cv.runs <- name(s) of the cross-validation types. folders will be generated with these names in create.folders.for.maxent.
232. # recommended framework is to treat each cross-validation fold as a letter.
233. # In the case of 5-fold cross validation, for running all the data (using the null.aic and optimize.maxent.likelihood . . .
234. # . . . functions), you should specify 'abcde'. For running cross validation models with maxent.crossval.error, you . . .
235. # . . . should specify c('abcd', 'abce', 'abde', 'acde', 'bcde').
236. ## f.class <- name(s) of the feature classes used in analysis.
237. ## beta.values <- regularization multipliers used in analysis
239. # create.summary.df and create.eval.df will create every possible combination of cv.runs, f.class, and beta.values
242. # function to create the null data.frame
243. create.null.df <- function(taxon.name, time.bin, extent){
244. null.object <- as.data.frame(matrix(NA, nrow=1, ncol=17))
245. colnames(null.object) <- c( "Taxa", "Time.Bin", "Extent", "CrossVal", "cv.num", "Features", "Betas", "n", "k",
246. "ln.L", "AIC", "AICc", "delta.i", "delta.i.c", "w.i", "w.i.c", "lambdas" )
247. null.object$Taxa <- taxon.name
248. null.object$Time.Bin <- time.bin
249. null.object$Extent <- extent
250. null.object$CrossVal <- 'abcde'
251. null.object$cv.num <- 1
252. null.object$lambdas <- 'Incercept'
254. return(null.object)
255. }
258. # function to create the summary data.frame
259. create.summary.df <- function(taxon.name,
260. time.bin,
261. extent,
262. cv.runs = 'abcde',
263. f.class = c('LQP', 'Q'),
264. 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) ){
266. summary.object <- expand.grid(CrossVal=cv.runs, Features = f.class, Betas = beta.values, stringsAsFactors = T)
267. summary.object$cv.num <- as.numeric(summary.object$CrossVal)
269. # ifelse functions defining what to do if taxon.name, time.bin, and extent are not specified
270. if( !is.null(taxon.name) ){
271. summary.object$Taxa <- taxon.name
272. } else {
273. summary.object$Taxa <- 'Taxon1'
274. }
275. if( !is.null(time.bin) ){
276. summary.object$Time.Bin <- time.bin
277. } else {
278. summary.object$Time.Bin <- 'TimeBin1'
279. }
280. if( !is.null(extent) ){
281. summary.object$Extent <- extent
282. } else {
283. summary.object$Extent <- 'Extent1'
284. }
286. # re-ordering columns
287. summary.object <- summary.object[,c(5:7, 1, 4, 2:3)]
288. # adding in all the other parameters
289. summary.object$n <- NA # sample size
290. summary.object$k <- NA # number of non-zero lambdas
291. summary.object$ln.L <- NA # log-likelihood
292. summary.object$AIC <- NA # AIC
293. summary.object$AICc <- NA # AICc corrected for small sample size
294. summary.object$delta.i <- NA # delta.i for AIC - wont calculate everything until all models have been run
295. summary.object$delta.i.c <- NA # delta.i for AICc - wont calculate everything until all models have been run
296. summary.object$w.i <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run
297. summary.object$w.i.c <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run
298. summary.object$lambdas <- NA # list of all the non-zero lambdas
299. summary.object[,4] <- as.character(summary.object[,4])
300. summary.object[,6] <- as.character(summary.object[,6])
302. return(summary.object)
303. }
306. # function to create the eval data.frame
307. create.eval.df <- function(taxon.name,
308. time.bin,
309. extent,
310. cv.runs = c('abcd', 'abce', 'abde', 'acde', 'bcde'),
311. f.class = 'LQP',
312. beta.values = 1 ){
314. summary.object <- expand.grid(CrossVal=cv.runs, Features = f.class, Betas = beta.values, stringsAsFactors = T)
315. summary.object$cv.num <- 1:length(cv.runs)
317. # ifelse functions defining what to do if taxon.name, time.bin, and extent are not specified
318. if( !is.null(taxon.name) ){
319. summary.object$Taxa <- taxon.name
320. } else {
321. summary.object$Taxa <- 'Taxon1'
322. }
323. if( !is.null(time.bin) ){
324. summary.object$Time.Bin <- time.bin
325. } else {
326. summary.object$Time.Bin <- 'TimeBin1'
327. }
328. if( !is.null(extent) ){
329. summary.object$Extent <- extent
330. } else {
331. summary.object$Extent <- 'Extent1'
332. }
334. # re-ordering columns
335. summary.object <- summary.object[,c(5:7, 1, 4, 2:3)]
336. # adding in all the other parameters
337. summary.object$n <- NA # sample size
338. summary.object$k <- NA # number of non-zero lambdas
339. summary.object$ln.L <- NA # log-likelihood
340. summary.object$AIC <- NA # AIC
341. summary.object$AICc <- NA # AICc corrected for small sample size
342. summary.object$delta.i <- NA # delta.i for AIC - wont calculate everything until all models have been run
343. summary.object$delta.i.c <- NA # delta.i for AICc - wont calculate everything until all models have been run
344. summary.object$w.i <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run
345. summary.object$w.i.c <- NA # Akaike Weights for AIC - wont calculate everything until all models have been run
346. summary.object$lambdas <- NA # list of all the non-zero lambdas
347. summary.object[,4] <- as.character(summary.object[,4])
348. summary.object[,6] <- as.character(summary.object[,6])
350. # renaming the columns
351. colnames(summary.object) <- c('Taxa', 'Time.Bin', 'Extent', 'CrossVal', 'cv.num', 'Features', 'Betas', 'n', 'thresh', 'test.sens',
352. 'pROC_0.1', 'pval_0.1', 'pROC_1', 'pval_1', 'pROC_5', 'pval_5', 'lambdas')
354. return(summary.object)
355. }
358. # function to create the nested fie structure that it will use to store the output (including html files) of maxent models
359. # summary.eval.df <- the summary.df data.frame or the eval.df data.frame.
360. # the function will use information from the cv.runs, f.class, and beta.values columns to create this nested file structure
361. # wd <- working directory that the nested file structure is going to be generated in. if running multiple species, it is recommended . . .
362. # . . . that you make a folder for each species, then run this function in each species folder
363. create.folders.for.maxent <- function(summary.eval.df, wd = getwd() ){
365. # setting the working directory
366. setwd( getwd() )
368. # for loop that goes through the summary/evaluation data.frame and creates the nested file structure
369. for( i in 1:nrow(summary.eval.df) ){
371. if( !dir.exists( paste(getwd(), summary.eval.df$CrossVal[i], sep='/' ) ) ){ # if the cross-validation folder exists
372. writeLines( c('Creating folder:', paste(getwd(), summary.eval.df$CrossVal[i], sep='/'), sep='') )
373. dir.create( paste(getwd(), summary.eval.df$CrossVal[i], sep='/') )
374. setwd( paste(getwd(), summary.eval.df$CrossVal[i], sep='/') )
375. } else {
376. setwd( paste(getwd(), summary.eval.df$CrossVal[i], sep='/') )
377. }
379. if( !dir.exists( paste(getwd(), summary.eval.df$Features[i], sep='/' ) ) ){ # if the feature class folder exists
380. writeLines( c('Creating folder:', paste(getwd(), summary.eval.df$Features[i], sep='/'), sep='') )
381. dir.create( paste(getwd(), summary.eval.df$Features[i], sep='/') )
382. setwd( paste(getwd(), summary.eval.df$Features[i], sep='/') )
383. } else {
384. setwd( paste(getwd(), summary.eval.df$Features[i], sep='/') )
385. }
387. if( !dir.exists( paste(getwd(), summary.eval.df$Features[i], sep='/' ) ) ){ # if the regularization folder exists
388. writeLines( c('Creating folder:', paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/'), sep='') )
389. dir.create( paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/') )
390. setwd( paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/') )
391. } else {
392. setwd( paste(getwd(), paste('beta', summary.eval.df$Betas[i], sep='_'), sep='/') )
393. }
395. setwd('../../../') # go three directories back
397. } # finishing main for loop
402. }
405. # setwd("~/Dropbox/SVP_Models/ModelOutput/Tyrano")
406. #
407. #
408. # examp <- create.summary.df('Tyrano', 'K', 'Laurimidia')
409. # examp2 <- create.eval.df('Tyrano', 'K', 'Laurimidia')
410. #
411. #
412. # create.folders.for.maxent(examp)
413. # create.folders.for.maxent(examp2)
418. ##### Optimize MaxEnt Likelihood -------------------------------------------------------------------------------------------
420. # function that runs all the maxent models.
421. # function will return a filled out summary model
423. optimize.maxent.likelihood <- function(summary.df, # summary object that will keep all the model output
424. # the nested file structure created from create.folders.for.maxent MUST exist
425. occs, # species occurrence object
426. background, # sampled Background object (COLUMNS MUST BE IDENTICAL TO occs)
427. predic, # column numbers of the predictor variables
428. first.occ.col, # number of the first cross-validation column in occs/background
429. home=getwd(), # directory where all the models will be ran.
430. all.models = TRUE # do you want to keep all versions of all the models?
431. # helpful for quickly checking some models, but may consume loads (e.g., >1GB) . . .
432. # . . . of hard disk space. Recommended to set to FALSE for exploratory analyses.
433. # if FALSE, function will create a folder called "RunOver" and will . . .
434. # . . . continuously write-over it for all models, and the only model you see . . .
435. # . . . at the end will be the last model that was ran
436. ){
437. # Calculating the total number of models
438. nmodels <- nrow(summary.df)
440. # prompting the user if they want to store models in the RAM
441. print(paste0('The time is ', Sys.time(), '. You are running ', nmodels, ' total MaxEnt models.'))
443. # setting up the progress bar
444. prog <- txtProgressBar(min=0, max=nrow(summary.df), style=3, char='+')
446. for(i in 1:nrow(summary.df)){
448. # setting the working directory for each folder
449. if(all.models == TRUE){
450. setwd( paste(home, summary.df[i,4], summary.df[i,6], paste('beta', summary.df[i,7], sep='_'), sep='/') )
451. } else {
452. # create the RunOver folder
453. if( !dir.exists( paste(home, 'RunOver', sep='/') ) ){
454. dir.create( paste(home, 'RunOver', sep='/') )
455. setwd( paste(home, 'RunOver', sep='/') )
456. } else {
457. setwd( paste(home, 'RunOver', sep='/') )
458. }
459. }
461. ###########################################################################
463. ### MaxEnt things happen here
465. ### preparing the data for the maxent model
466. # filtering the occ object by it's respective cross-validation identity
467. cv.number <- summary.df$cv.num[i]
468. # assigning the column number to be sent through maxent
469. col.number <- first.occ.col + cv.number - 1
470. # filtering the species dataset by col.number and assigning to summary.df
471. sp <- occs %>% dplyr::filter(occs[,col.number] == 1)
472. n <- nrow(sp)
473. summary.df[i,8] <- n
474. # bending the occ and background data together
475. mx.data <- rbind(sp, background)
477. ### running the actual maxent model
478. mx.model <- dismo::maxent(p = mx.data[,col.number],
479. x = mx.data[,predic],
480. path = paste0(getwd()),
481. args = prepPara(userfeatures = summary.df[i,6], betamultiplier = summary.df[i,7], doclamp = FALSE)
482. )
484. ### calculating k and the names of the lambdas
485. lambda.file <- as.data.frame(mx.model@lambdas) %>% `colnames<-`('lambdas')
486. # lambdas data frame
487. lambdas.df <- splitstackshape::cSplit(lambda.file, sep=',', splitCols='lambdas')
488. colnames(lambdas.df) <- c('feature', 'lambda', 'min', 'max')
489. # finding the non-zero lambdas
490. non.zero.lambdas <- lambdas.df %>% dplyr::filter(!is.na(max)) %>% dplyr::filter(lambda != 0)
491. class(non.zero.lambdas) <- 'data.frame'
492. # assigning the number of parameters
493. k <- nrow(non.zero.lambdas)
494. # if beta is too high, the model gets over-regularized to the point that all the lambda coefficients get set to 0
495. # this effectively becomes an intercept only model, which is effectively the global mean
496. # in this sense, k should get set to 0
497. if(k == 0){
498. k <- 1
499. }
500. summary.df[i,9] <- k
501. # giving noting the variables/features/hyperparameters with non-zero lambdas
502. summary.df[i,17] <- toString(non.zero.lambdas[,1])
504. ### calculating the log-likelihood, then AIC and AICc
505. # if statement calculating if there is an appropriate AIC value
506. # e.g., can't fit 4 observations (occurrence points) with 5 variables
507. if(n - k < 2){
508. summary.df[i,10] <- NA
509. summary.df[i,11] <- NA
510. summary.df[i,12] <- NA
511. } else { # if it is possible to calculate AIC and AICc
512. # logistic model output
513. mx.back <- dismo::predict(mx.model, background[,predic])
514. # sum of all background point values - mx.back / back.sum should = 1
515. back.sum <- sum(mx.back)
516. # logistic values of the (k-1)/k occurrences
517. mx.occs <- dismo::predict(mx.model, sp[,predic])
518. # scaling to make compatible for calculating AIC
519. occs.raw <- mx.occs / back.sum
520. # log(likelihood)
521. log.like <- sum(log(occs.raw))
522. summary.df[i,10]<- log.like
524. ### calculating AIC and AICc
525. # AIC
526. AIC <- 2*k - 2*log.like
527. summary.df[i,11] <- AIC
528. # AICc
529. summary.df[i,12] <- AIC + 2*((k^2 + k) / (n - k - 1))
531. }
532. # updating the prograss bar for each run to get an idea of how long things will take
533. setTxtProgressBar(prog, i)
536. ###########################################################################
537. # returning to the home directory
538. setwd(home)
540. } # closes the for loop
542. return(summary.df)
543. }
548. ##### Null AICc ------------------------------------------------------------------------------------------------------------
550. # function that calculates AIC and AICc values for a null intercept-only model
551. # the arguments are the same as the optimize.maxent.likelihood model, except that null.df should be a data.frame created from the . . .
552. # . . . create.null.df object
553. # function will return
555. null.aic <- function(null.df, occs, background, first.occ.col){
557. # number of parameters
558. null.df$k <- 1
560. # assigning the column number to be sent through maxent
561. col.number <- first.occ.col
562. # filtering the species dataset by col.number and assigning to null.df
563. sp <- occs %>% dplyr::filter(occs[,col.number] == 1)
564. n <- nrow(sp)
565. null.df[1,8] <- n
567. # giving noting the variables/features/hyperparameters with non-zero lambdas
568. null.df[1,17] <- 'Intercept'
570. out.scale <- rep(1, times=n ) / nrow(background)
572. log.like <- sum(log(out.scale))
574. null.df[1,10] <- log.like
576. ### calculating AIC and AICc
577. # AIC
578. AIC <- 2 - 2*log.like
579. null.df[1,11] <- AIC
580. # AICc for one term null model
581. null.df[1,12] <- AIC + 4/(n - 2)
583. return(null.df)
584. }
589. ##### MaxEnt Cross-Validation Error -----------------------------------------------------------------------------------------------
591. # function that runs five-fold cross-validation once you've found the optimum hyper-parameters for your maxent model(s)
592. # functions returns a list containing 3 objects: 1.) a filled out eval object; 2.) an object containing all the maxent models; and . . .
593. # . . . 3.) a data.frame with dimensions [ 1:nrow(background), 1:nrow(eval.df) ] containing the projections of all the models in
595. maxent.crossval.error <- function(eval.df, # object generated from create.eval.df that will keep all the model output
596. occs, # species occurrence object
597. background, # sampled Background object (COLUMNS MUST BE IDENTICAL TO occs)
598. predic, # column numbers of the predictor variables
599. first.occ.col, # number of the first training in occs/background (assumes k = 5)
600. # NOTE: this is not the presence column! it's the column to the right of it
601. first.test.col, # number of the first testing column in occs/background (assumes k = 5)
602. home=getwd(), # directory where all the models will be ran
603. all.models = TRUE, # do you want to keep all versions of all the models?
604. # helpful for quickly checking some models, but may consume loads (e.g., >1GB) . . .
605. # . . . of hard disk space. Recommended to set to FALSE for exploratory analyses.
606. # if FALSE, function will create a folder called "RunOver" and will . . .
607. # . . . continuously write-over it for all models, and the only model you see . . .
608. # . . . at the end will be the last model that was ran
609. omission.rate=0, # user specified omission rate to be calculated for the threshold
610. # express as a proportion from 0-1
611. all.background, # background data from the extents the species exists in
612. # if using all the potential background points, all.background should == background
613. pROC.reps = 500 # number of iterations each partialROC test will go through
614. ){
615. # Calculating the total number of models
616. nmodels <- nrow(eval.df)
618. # stop if an incompatible omission rate is specified
619. if(omission.rate > 1 || omission.rate < 0){
620. stop('Specify an omission rate between 0-1.')
621. }
623. # prompting the user if they want to store models in the RAM
624. print(paste0('The time is ', Sys.time(), '. You are running ', nmodels, ' total MaxEnt models.'))
626. # setting up the progress bar
627. prog <- txtProgressBar(min=0, max=nrow(eval.df), style=3, char='+')
629. # setting up various list objects to send output to
630. # list for the maxent models
631. maxent.list <- list()
632. # list for the testing background data
633. test.back.list <- list()
634. # list for all the occs
635. all.occs.list <- list()
637. for(i in 1:nrow(eval.df)){
639. # setting the working directory for each folder
640. if(all.models == TRUE){
641. setwd( paste(home, eval.df[i,4], eval.df[i,6], paste('beta', eval.df[i,7], sep='_'), sep='/') )
642. } else {
643. # create the RunOver folder
644. if( !dir.exists( paste(home, 'RunOver', sep='/') ) ){
645. dir.create( paste(home, 'RunOver', sep='/') )
646. setwd( paste(home, 'RunOver', sep='/') )
647. } else {
648. setwd( paste(home, 'RunOver', sep='/') )
649. }
650. }
652. ###########################################################################
654. ### MaxEnt things happen here
656. ### preparing the data for the maxent model
657. # filtering the occ object by it's respective cross-validation identity
658. cv.number <- eval.df$cv.num[i]
659. # assigning the column number to be sent through maxent
660. col.number <- first.occ.col + cv.number - 1
661. # filtering the species dataset by col.number and assigning to eval.df
662. sp <- occs %>% dplyr::filter(occs[,col.number] == 1)
663. n <- nrow(sp)
664. eval.df[i,8] <- n
666. # generating the testing data
668. test.col.number <- first.test.col + cv.number - 1
669. sp.test <- occs %>% dplyr::filter(occs[,test.col.number] == 1)
672. # bending the occ and background data together
673. mx.data <- rbind(sp, background)
675. ### running the actual maxent model
676. mx.model <- dismo::maxent(p = mx.data[,col.number],
677. x = mx.data[,predic],
678. path = paste0(getwd()),
679. args = prepPara(userfeatures = eval.df[i,6], betamultiplier = eval.df[i,7], doclamp = FALSE)
680. )
682. ### calculating k and the names of the lambdas
683. lambda.file <- as.data.frame(mx.model@lambdas) %>% `colnames<-`('lambdas')
684. # lambdas data frame
685. lambdas.df <- splitstackshape::cSplit(lambda.file, sep=',', splitCols='lambdas')
686. colnames(lambdas.df) <- c('feature', 'lambda', 'min', 'max')
687. # finding the non-zero lambdas
688. non.zero.lambdas <- lambdas.df %>% dplyr::filter(!is.na(max)) %>% dplyr::filter(lambda != 0)
689. class(non.zero.lambdas) <- 'data.frame'
690. # giving noting the variables/features/hyperparameters with non-zero lambdas
691. eval.df[i,17] <- toString( non.zero.lambdas[,1] )
694. ## evaluating the models
695. # predicting the training data
696. mx.train.occ <- dismo::predict(mx.model, sp[,predic])
697. # predicting the training background data
698. mx.train.back <- dismo::predict(mx.model, background[,predic])
699. # predicting the testing data
700. mx.test.occ <- dismo::predict(mx.model, sp.test[,predic])
701. # projecting to all extents
702. mx.test.back <- dismo::predict(mx.model, all.background[,predic])
703. # projecting to all occ points
704. mx.all.occs <- dismo::predict(mx.model, occs[,predic])
706. # calculate the threshold
707. if(omission.rate == 0){ # if using the LTP threshold, dismo calculates slightly too high of a threshold
708. thresh <- min(mx.train.occ)
709. } else {
710. # make a model evaluation object
711. mx.eval <- dismo::evaluate(p=mx.train.occ, a=mx.train.back)
712. thresh <- dismo::threshold(mx.eval, stat='sensitivity', sensitivity= (1 - omission.rate) )
713. }
714. # assigning the threshold value to the output file
715. eval.df[i,9] <- thresh
716. # calculating the test sensitivity
717. sens <- length(which(mx.test.occ >= thresh)) / length(mx.test.occ)
718. eval.df[i,10] <- sens
720. # making a raster of the testing background data
721. # for some reason, kuenm calculates wonky AUC_ratio values (i.e., > 2, which is impossible) unless you specify a raster
722. r <- raster(nrows=1, ncols=length(mx.test.back) )
723. r[r] <- mx.test.back
725. ## calculating AUC_ratios from a partialROC test
726. # error = 0.1%
727. pROC_0.1 <- kuenm::kuenm_proc(occ.test = mx.test.occ, # numeric vector of the predicted suitability values on the testing data
728. model = r, # raster model of the predicted suitability values for the background
729. threshold = 0.1, # potential error threshold (expressed as a percent)
730. rand.percent = 50, # percentage of data to be used in each bootstrap rep
731. iterations = pROC.reps # number of repititions
732. )
733. # assigning the average AUC_ratio from pROC.reps iterations to eval.df
734. eval.df[i,11] <- as.numeric(pROC_0.1$pROC_summary[1])
735. # assigning the partialROC p-value to eval.df
736. eval.df[i,12] <- as.numeric(pROC_0.1$pROC_summary[2])
738. # error = 1%
739. pROC_1 <- kuenm::kuenm_proc(occ.test = mx.test.occ, # numeric vector of the predicted suitability values on the testing data
740. model = r, # raster model of the predicted suitability values for the background
741. threshold = 1, # potential error threshold (expressed as a percent)
742. rand.percent = 50, # percentage of data to be used in each bootstrap rep
743. iterations = pROC.reps # number of repititions
744. )
745. # assigning the average AUC_ratio from pROC.reps iterations to eval.df
746. eval.df[i,13] <- as.numeric(pROC_1$pROC_summary[1])
747. # assigning the partialROC p-value to eval.df
748. eval.df[i,14] <- as.numeric(pROC_1$pROC_summary[2])
750. # error = 5%
751. pROC_5 <- kuenm::kuenm_proc(occ.test = mx.test.occ, # numeric vector of the predicted suitability values on the testing data
752. model = r, # raster model of the predicted suitability values for the background
753. threshold = 5, # potential error threshold (expressed as a percent)
754. rand.percent = 50, # percentage of data to be used in each bootstrap rep
755. iterations = pROC.reps # number of repititions
756. )
757. # assigning the average AUC_ratio from pROC.reps iterations to eval.df
758. eval.df[i,15] <- as.numeric(pROC_5$pROC_summary[1])
759. # assigning the partialROC p-value to eval.df
760. eval.df[i,16] <- as.numeric(pROC_5$pROC_summary[2])
763. mx.test.back.df <- as.data.frame( as.matrix(mx.test.back, ncol=1) )
765. # assigning the objects to the various lists
766. maxent.list[[i]] <- mx.model
767. test.back.list[[i]] <- mx.test.back
768. all.occs.list[[i]] <- mx.all.occs
770. # updating the prograss bar for each run to get an idea of how long things will take
771. setTxtProgressBar(prog, i)
774. ###########################################################################
775. # returning to the home directory
776. setwd(home)
778. } # closes the for loop
780. # bind the testing background data into a single data
781. back.projections <- as.data.frame( do.call('cbind', test.back.list ) )
782. colnames(back.projections) <- eval.df$CrossVal
784. # bind all occs together
785. occ.projections <- as.data.frame(do.call('cbind', all.occs.list))
786. colnames(occ.projections) <- eval.df$CrossVal
788. # make a list of the output
789. out.list <- list()
790. out.list$maxent.models <- maxent.list
791. out.list$back.projection <- back.projections
792. out.list$occ.projection <- occ.projections
793. out.list$summary <- eval.df
795. return(out.list)
796. }
801. ##### MaxEnt Evaluation ----------------------------------------------------------------------------------------------------
803. # function that takes the eval object generated from maxent.crossval.error and calculates the weighted mean and standard deviation
804. # the weighted mean is calculating my testing sensitivity (1 - omission rate) multiplied by the partial ROC/AUC value.
805. # this ensures that if a model does not discriminate between presences/non-presences well, it will receive comparatively lower weight
806. # later package versions will include Boyce index as an additional calibration technique and will offer the user the ability to . . .
807. # . . . choose which metrics to use for assessing model reliability
808. # the weighted mean and standard deviation can be plotted to infer model variability/uncertainty
810. ### ARGUMENTS ###
811. ## eval <- eval object generated from maxent.crossval.error that will project the model to every grid cell used in the training region
812. ## pROC.error <- the user-specified omission rate.
813. # (choose from 0.1, 1, 5 - however they almost always end up super correlated with each other)
814. maxent.eval <- function(eval, pROC.error=1){
816. # extracting model sensitivity
817. sens <- eval$summary$test.sens
819. # which pROC error amount to use?
820. if(pROC.error == 0.1){
821. AUC_ratio <- eval$summary$pROC_0.1
822. } else if(pROC.error == 1){
823. AUC_ratio <- eval$summary$pROC_1
824. } else if(pROC.error == 5){
825. AUC_ratio <- eval$summary$pROC_5
826. } else {
827. stop('Select an appropriate partialROC error amount.')
828. }
830. # weights = sensitivity*AUC_ratio
831. weights <- sens * AUC_ratio
832. weights[is.nan(weights)] <- 0
834. # making a matrix of the background points
835. models.mat <- as.matrix(eval$back.projection)
836. # weighted means and standard deviations
837. w.mean <- matrixStats::rowWeightedMeans(x=models.mat, w=weights)
838. w.sd <- matrixStats::rowWeightedSds(x=models.mat, w=weights)
839. # merging the weighted means/sd's of all background points
840. back.df <- data.frame(w.mean, w.sd)
842. # making a matrix of the background points
843. occs.mat <- as.matrix(eval$occ.projection)
844. # weighted means and standard deviations
845. w.mean <- matrixStats::rowWeightedMeans(x=occs.mat, w=weights)
846. w.sd <- matrixStats::rowWeightedSds(x=occs.mat, w=weights)
847. # merging the weighted means/sd's of all occ points
848. occ.df <- data.frame(w.mean, w.sd)
851. # making the output list
852. out.list <- list()
853. out.list$back <- back.df
854. out.list$occ <- occ.df
855. out.list$weights <- weights
857. return(out.list)
859. }
863. ##### Calculating the MaxEnt Threshold -------------------------------------------------------------------------------------
865. # function that calculates the threshold value for all the training data
867. ### ARGUMENTS ###
868. ## eval <- eval object generated from maxent.crossval.error that will project the model to every grid cell used in the training region
869. ## occs <- the occurrence data.frame
870. ## predic <- column numbers of the predictor variables
871. ## pROC.error <- the user-specified omission rate.
872. # (choose from 0.1, 1, 5 - however they almost always end up super correlated with each other)
874. # function will return the lowest training preference threshold.
875. # future package versions will give the opportunity to select different thresholds
877. maxent.thresh <- function(eval, occs, predic, pROC.error=1){
879. n.rows <- nrow(occs)
880. n.cols <- nrow(eval$summary)
882. # making a data.frame of the occurrences
883. occ.values <- as.data.frame( matrix(NA, nrow=n.rows, ncol=n.cols ) )
884. colnames(occ.values) <- eval$summary$CrossVal
886. # for loop predicting all the cross validation maxent models
887. for(i in 1:nrow(eval$summary) ){
888. mx.occs <- dismo::predict(object=eval$maxent.models[[i]], x=occs[,predic])
889. occ.values[,i] <- mx.occs
890. }
892. #
893. # extracting model sensitivity
894. sens <- eval$summary$test.sens
896. # which pROC error amount to use?
897. if(pROC.error == 0.1){
898. AUC_ratio <- eval$summary$pROC_0.1
899. } else if(pROC.error == 1){
900. AUC_ratio <- eval$summary$pROC_1
901. } else if(pROC.error == 5){
902. AUC_ratio <- eval$summary$pROC_5
903. } else {
904. stop('Select an appropriate partialROC error amount.')
905. }
907. # weights = sensitivity*AUC_ratio
908. weights <- sens * AUC_ratio
909. weights[is.nan(weights)]<- 0
912. # making a matrix of the background points
913. models.mat <- as.matrix(occ.values)
914. # weighted means
915. w.mean <- matrixStats::rowWeightedMeans(x=models.mat, w=weights)
917. # setting the threshold as the lowest training presence
918. ltp <- min(w.mean)
920. return(ltp)
921. }
925. ##### Projecting the Best Models to All the Extents ------------------------------------------------------------------------
927. # function to project models to any and all extents you want to
928. # effectively making new master occurrence and master background files that can easily be saved and re-loaded again
930. maxent.everything <- function(eval, means, everything, thresh, pROC.error=1, predic, name='pc'){
932. # making a data.frame of the occurrences
933. every.value <- as.data.frame( matrix(NA, nrow=nrow(everything), ncol=nrow(eval$summary)) )
934. colnames(every.value) <- eval$summary$CrossVal
936. # for loop predicting all the cross validation maxent models
937. for(i in 1:nrow(eval$summary) ){
938. mx.everything <- dismo::predict(object=eval$maxent.models[[i]], x=everything[,predic])
939. every.value[,i] <- mx.everything
940. }
942. # extracting model sensitivity
943. sens <- eval$summary$test.sens
945. # which pROC error amount to use?
946. if(pROC.error == 0.1){
947. AUC_ratio <- eval$summary$pROC_0.1
948. } else if(pROC.error == 1){
949. AUC_ratio <- eval$summary$pROC_1
950. } else if(pROC.error == 5){
951. AUC_ratio <- eval$summary$pROC_5
952. } else {
953. stop('Select an appropriate partialROC error amount.')
954. }
956. # weights = sensitivity*AUC_ratio
957. weights <- sens * AUC_ratio
958. weights[is.nan(weights)] <- 0
960. # making a matrix of the background points
961. models.mat <- as.matrix(every.value)
962. # weighted means
963. w.mean <- matrixStats::rowWeightedMeans(x=models.mat, w=weights)
965. # making the output data frame
966. out.df <- as.data.frame( matrix(NA, nrow=length(w.mean), ncol=2) )
967. colnames(out.df) <- c(paste0(eval$summary[1,1], '_contin'), paste0(eval$summary[1,1], '_thresh') )
968. # assigning the weighted means to the output data.frame
969. out.df[,1] <- w.mean
970. # replicating the weighted means so they can be thresolded
971. out.df[,2] <- w.mean
972. # calculating the threshold
973. thresh <- min(means$occ$w.mean)
974. # thresolding the values
975. out.df[,2][out.df[,2] >= thresh] <- 1
976. out.df[,2][out.df[,2] < thresh] <- 0
980. return(out.df)
982. }
987. ##### Informed Analysis -----------------------------------------------------------------------------------------------
989. # function that calculates Multivariate Environmental Suitability Surface (MESS), but allows for some extrapolation . . .
990. # . . . based on how response curves look.
992. informed.mess <- function(ref.extent, # reference data.frame - must contain same column names as 'data'
993. mess.extent, # data.frame containing the extent to be projected to - ideally this should be all extents
994. coord.cols=2:3, # column numbers of the coordinates of the data/ref.data (e.g., long/lat)
995. predic=6:11, # column numbers of ONLY the final predictor variables
996. tolerance=NULL, # tolerance vector to extrapolate beyond the limits of the data. see example below
997. # if tolerance is not specified (default), function will calculate basic MESS
998. keep.layers=FALSE # decision to retain the MESS values for all variables instead of just the final MESS data
999. # changing to TRUE can help easily identify which variables are causing the MESS value . . .
1000. # . . . to be so low (i.e., it may be only 1 variable that is responsible)
1001. # can easily re-check this later
1002. ){
1003. # isolating the coordinates and predictor variables
1004. new.vars <- mess.extent[, predic]
1005. ref.vars <- ref.extent[, predic]
1006. #
1007. if( !is.null(tolerance) ){
1008. if(length(predic) != length(tolerance)/2 ){
1009. stop('Ensure that tolerance is 2x as long as predic')
1010. }
1011. # changing tolerance to a data.frame
1012. tolerance <- as.data.frame( matrix(tolerance, ncol=2, byrow=TRUE) )
1013. nvars <- length(predic)
1014. # pre-changing values outside the range of the input data
1015. # extra rows
1016. extra.rows <- matrix(NA, nrow=2, ncol=nvars)
1017. colnames(extra.rows) <- colnames(ref.vars)
1018. ref.vars <- rbind(ref.vars, extra.rows )
1019. for(i in 1:nvars){
1020. # if we can extrapolate beyond the lower end of variable i, allow mess to not
1021. if(tolerance[i,1] == 1){
1022. min.var <- min(ref.vars[,i], na.rm=TRUE) - 0.0001
1023. new.vars[ new.vars[,i] < min.var, i ] <- min.var
1024. ref.vars[nrow(ref.vars)-1, i] <- min.var
1025. } else {
1026. ref.vars[nrow(ref.vars)-1, i] <- min(ref.vars[,i], na.rm=TRUE)
1027. }
1029. if(tolerance[i,2] == 1){
1030. max.var <- max(ref.vars[,i], na.rm=TRUE) + 0.0001
1031. new.vars[ new.vars[,i] > max.var, i ] <- max.var
1032. ref.vars[nrow(ref.vars), i] <- max.var
1033. } else {
1034. ref.vars[nrow(ref.vars), i] <- max(ref.vars[,i], na.rm=TRUE)
1035. }
1037. } # end for(i in 1:nvars) loop
1039. }
1041. # running the mess analysis
1042. mess.vars <- as.data.frame( sapply(1:ncol(new.vars), function(i) .messi3(new.vars[, i], ref.vars[, i])) )
1043. # re-asigning of the extrapolating points to have a mess value of 0.
1044. nref <- nrow(ref.vars)
1045. fix.vars <- as.data.frame( sapply(1:ncol(new.vars), function(i) .messi3(ref.vars[ (nref-1):nref , i], ref.vars[, i])) )
1046. for(i in 1:nvars){
1047. mess.vars[mess.vars[,i] == fix.vars[1,i],i] <- 0
1048. mess.vars[mess.vars[,i] == fix.vars[2,i],i] <- 0
1050. }
1052. # making a simple thresholded version of the mess analysis for easy plotting
1053. final.mess <- as.data.frame( apply(mess.vars, 1, min) )
1054. colnames(final.mess) <- 'mess.raw'
1055. mess.thresh <- final.mess$mess.raw
1056. mess.thresh[mess.thresh > 0] <- 1
1057. mess.thresh[mess.thresh < 0] <- -1
1058. final.mess <- cbind(final.mess, mess.thresh)
1060. # deciding whether to keep individual mess layers
1061. if(keep.layers==TRUE){
1062. colnames(ref.vars) <- paste0('mess_', colnames(ref.vars))
1063. final.mess <- cbind(final.mess, ref.vars)
1064. return(final.mess)
1065. } else {
1066. return(final.mess)
1067. }
1069. }
1073. # internal function originally from dismo that runs the actual mess analysis in data.frame format
1074. .messi3 <- function(p,v) { # p=new.vars v=ref.vars
1075. # seems 2-3 times faster than messi2
1076. v <- stats::na.omit(v)
1077. f <- 100*findInterval(p, sort(v)) / length(v)
1078. minv <- min(v)
1079. maxv <- max(v)
1080. res <- 2*f
1081. f[is.na(f)] <- -99
1082. i <- f>50 & f<100
1083. res[i] <- 200-res[i]
1085. i <- f==0
1086. res[i] <- 100*(p[i]-minv)/(maxv-minv)
1087. i <- f==100
1088. res[i] <- 100*(maxv-p[i])/(maxv-minv)
1089. res
1090. }
1093. ##### uncert.suit.plot -----------------------------------------------------------------------------------------------------
1095. suit.uncert.plot <- function(means){
1097. back <- means$back
1098. occs <- means$occ
1100. ltp <- min(occs$w.mean)
1102. output <- ggplot(data=back, aes(x=w.mean, y=w.sd)) + geom_point(colour='black', size=0.75) +
1103. geom_point(data=occs, aes(x=w.mean, y=w.sd), size=2, colour='red', shape=18 ) +
1104. ylim(0, 0.55) + xlim(0,1) + theme_classic() + coord_fixed(1/0.55) +
1105. annotate(geom='text', x=0.05, y=0.45, label=round(ltp, 4), hjust=0 )
1107. return(output)
1109. }
1112. # function that adds k-folds to the dataset
1113. ### Arguments ###
1114. # x <- vector or data.frame of your occurrences
1115. #
1116. # k <- number of folds you want to make
1117. #
1118. # seed <- value for set.seed in case you want to reproduce your exact k-fold samples in the future
1120. # returns a vector that has the same length as nrow(x) that contains the k-fold bin that the occurrences got assigned to
1122. make.kfolds <- function(x, k=5, seed=NULL){
1123. if(class(x) != "data.frame"){
1124. class(x) <- "data.frame"
1125. }
1126. n <- nrow(x)
1127. rep.times <- n %/% k # number of full reps for the rep function
1128. k.bins <- rep.times*k # number of occs minus any remainder when dividing by k
1129. remainder <- n - k.bins # find the remainder
1130. output <- rep(1:k, rep.times) # make k bins of equal size
1131. if(!is.null(seed)){
1132. set.seed(seed)
1133. }
1134. if(remainder != 0){ # if there is a remainder, . . .
1135. extra <- sample(x=1:k, size=remainder) # randomly sample it. . .
1136. output <- c(output, extra) # . . . and add it to the output
1137. }
1138. output <- sample(output, size=n) # randomize the order of the output
1139. return(output)
1140. }
1145. # function that makes a background file for each extent (time.bin and/or region) in the analysis
1146. # merging the background files for every extent will create the master background file
1148. ### ARGUMENTS ###
1149. # r.stack <- a raster stack containing all the predictor variables for analysis for a given extent
1150. # if using multiple extents, ensure that all the predictor variables have the exact same name and are in the same order
1151. # if testing multiple types of variables (i.e., GCM-based vs sedimentology-based), it is ideal to seperate them into . . .
1152. # . . . those two categories before reading them into R
1153. #
1154. # x.col <- the column name of the x-coordinate to be put in the background/occurrence objects
1155. # this MUST be consistent for all background/occurrence objects
1156. #
1157. # y.col <- the column name of the y-coordinate to be put in the background/occurrence objects
1158. # this MUST be consistent for all background/occurrence objects
1159. #
1160. # time.bin <- the name of the time.bin to be put in the background/occurrence objects
1161. # if only using a single time.bin, keep this consistent for all background/occurrence objects
1162. #
1163. # region <- the name of the region to be put in the background/occurrence objects
1164. # if only using a single region, keep this consistent for all background/occurrence objects
1165. #
1166. # k.folds <- the number of folds you want to run for cross-validation
1167. # k stops at 26 because you don't want to run 27+ fold cross-validation. The resulting data.frame becomes unwieldy.
1170. make.background.df <- function(r.stack, x.col='long', y.col='lat', time.bin='time1', region=1, k.folds=5){
1172. # ensuring that k.folds is a positive integer between [1,26]
1173. if( !is.null(k.folds) ){
1174. if( k.folds < 1){
1175. k.folds <- 1
1176. print('NOTE: Making background data.frame without any k.folds.')
1177. } else if( k.folds > 26 ){
1178. k.folds <- 26
1179. warning('k.folds has been set to 26.')
1180. } else if( k.folds%%1 != 0 ){
1181. warning('k.folds has been rounded to the nearest whole number.')
1182. k.folds <- round(k.folds)
1183. }
1184. } else {
1185. k.folds <- 1
1186. print('NOTE: Making background data.frame without any k.folds.')
1187. }
1189. # extracting the env-variables to a data.frame
1190. coords.df <- raster::sampleRandom(r.stack, ncell(r.stack), xy=TRUE, sp=FALSE, na.rm=FALSE)
1191. colnames(coords.df)[1:2] <- c(x.col, y.col)
1192. coords.df <- coords.df[complete.cases(coords.df), ]
1193. n.points <- nrow(coords.df)
1195. # making the data.frame of the name, time.bin, region, and presence columns
1196. d.cols <- as.data.frame( matrix(0, nrow=n.points, ncol=4) )
1197. d.cols[,1] <- 'Background'
1198. d.cols[,2] <- time.bin
1199. d.cols[,3] <- region
1200. colnames(d.cols) <- c('Name', 'time.bin', 'region', 'presence')
1202. # merging everything together
1203. 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] ))
1204. colnames(out)[1] <- 'Name'
1205. colnames(out)[ ncol(out) ] <- 'presence'
1207. # making the cross-validation columns
1208. if( k.folds > 1 ){
1209. ltrs <- letters[1:k.folds]
1210. # training columns
1211. train <- as.data.frame( matrix(0, nrow=n.points, ncol=k.folds) )
1212. for(i in 1:k.folds){
1213. colnames(train)[i] <- paste0('tr.', paste(ltrs[-(k.folds+1-i)], collapse='') )
1214. }
1215. # testing columns
1216. test <- as.data.frame( matrix(0, nrow=n.points, ncol=k.folds) )
1217. for(i in 1:k.folds){
1218. colnames(test)[i] <- paste0('test.', paste(ltrs[(k.folds+1-i)], collapse='') )
1219. }
1221. out <- do.call('cbind', list( out, train, test ) )
1222. }
1224. # output
1225. return(out)
1226. }
1230. # function that makes an occurrence file for each extent (time.bin and/or region) in the analysis
1231. # merging the occurrence files for every extent will create the master occurrence file
1232. # MUST have the make.kfolds function also loaded
1234. ### ARGUMENTS ###
1235. # r.stack <- a raster stack containing all the predictor variables for analysis for a given extent
1236. # if using multiple extents, ensure that all the predictor variables have the exact same name and are in the same order
1237. # if testing multiple types of variables (i.e., GCM-based vs sedimentology-based), it is ideal to seperate them into . . .
1238. # . . . those two categories before reading them into R
1239. #
1240. # taxa.df <- a data.frame of occurrence with three columns that have:
1241. # 1.) the names of the taxa you are modeling
1242. # 2.) the x-coordinates of the occurrences
1243. # 3.) the y-coordinates of the occurrences
1244. # THESE COLUMNS MUST BE IN THIS ORDER!!!
1245. # this is the same format as the SWD (species with data) format for the regular maxent.jar file
1246. #
1247. # x.col <- the column name of the x-coordinate to be put in the background/occurrence objects
1248. # this MUST be consistent for all background/occurrence objects
1249. #
1250. # y.col <- the column name of the y-coordinate to be put in the background/occurrence objects
1251. # this MUST be consistent for all background/occurrence objects
1252. #
1253. # time.bin <- the name of the time.bin to be put in the background/occurrence objects
1254. # if only using a single time.bin, keep this consistent for all background/occurrence objects
1255. #
1256. # region <- the name of the region to be put in the background/occurrence objects
1257. # if only using a single region, keep this consistent for all background/occurrence objects
1258. #
1259. # k.folds <- the number of folds you want to run for cross-validation
1260. # k stops at 26 because you don't want to run 27+ fold cross-validation. The resulting data.frame becomes unwieldy.
1261. #
1262. # k.seed <- value for set.seed in case you want to reproduce your exact k-fold samples in the future
1263. # NOTE that this set.seed will apply exactly the same to every species you are modeling
1264. #
1266. 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){
1269. # ensuring that k.folds is a positive integer between [1,26]
1270. if( !is.null(k.folds) ){
1271. if( k.folds < 1){
1272. k.folds <- 1
1273. print('NOTE: Making occurrence data.frame without any k.folds.')
1274. } else if( k.folds > 26 ){
1275. k.folds <- 26
1276. warning('k.folds has been set to 26.')
1277. } else if( k.folds%%1 != 0 ){
1278. warning('k.folds has been rounded to the nearest whole number.')
1279. k.folds <- round(k.folds)
1280. }
1281. } else {
1282. k.folds <- 1
1283. print('NOTE: Making occurrence data.frame without any k.folds.')
1284. }
1286. # extracting the env-variables to a data.frame
1287. occs.df <- raster::extract(x=r.stack, y=taxa.df[,2:3])
1288. occs.df <- cbind(taxa.df, occs.df)
1289. colnames(occs.df)[2:3] <- c(x.col, y.col)
1290. occs.df <- occs.df[complete.cases(occs.df), ]
1291. colnames(occs.df)[1] <-'Name'
1292. n.occs <- nrow(occs.df)
1294. # making the data.frame of the name, time.bin, region, and presence columns
1295. d.cols <- as.data.frame( matrix(0, nrow=n.occs, ncol=3) )
1296. d.cols[,1] <- time.bin
1297. d.cols[,2] <- region
1298. colnames(d.cols) <- c('time.bin', 'region', 'presence')
1300. # merging everything together
1301. occs.df <- do.call('cbind', list(occs.df[,1:3], d.cols[,1:2], occs.df[,4:ncol(occs.df)], d.cols[,3] ) )
1302. colnames(occs.df)[ ncol(occs.df) ] <- 'presence'
1304. # making the cross-validation columns
1305. if( k.folds > 1 ){
1307. ltrs <- letters[1:k.folds]
1308. # training columns
1309. train <- as.data.frame( matrix(1, nrow=n.occs, ncol=k.folds) )
1310. for(i in 1:k.folds){
1311. colnames(train)[i] <- paste0('tr.', paste(ltrs[-(k.folds+1-i)], collapse='') )
1312. }
1313. # testing columns
1314. test <- as.data.frame( matrix(0, nrow=n.occs, ncol=k.folds) )
1315. for(i in 1:k.folds){
1316. colnames(test)[i] <- paste0('test.', paste(ltrs[(k.folds+1-i)], collapse='') )
1317. }
1319. # splitting up the occurrence data.frame by taxon, then assigning the k.folds
1320. occs.list <- base::split(occs.df, occs.df[,1])
1321. n.taxa <- length(occs.list)
1322. for(i in 1:n.taxa){
1324. # setting all taxa with < k.folds occurrences to 0's. they will later be converted to exist in all the training and testing subsets
1325. if( nrow( occs.list[[i]] ) < k.folds ){
1326. occs.list[[i]]$presence <- 0
1327. } else {
1328. # filling in cross-validation folds for taxa who have > k.folds in terms of occurrences
1329. occs.list[[i]]$presence <- make.kfolds(occs.list[[i]], k=k.folds, seed=k.seed)
1330. }
1333. } # closing for(i in 1:n.taxa)
1335. # merging all the different taxa back together
1336. occs.df <- do.call('rbind', occs.list)
1337. colnames(occs.df)[ ncol(occs.df) ] <- 'presence'
1338. row.names(occs.df) <- 1:nrow(occs.df)
1340. # assigning the k.fold parameters in occs.df
1341. for(i in 1:nrow(occs.df) ){
1343. # for taxa with fewer occs than k.folds
1344. if( occs.df$presence[i] == 0 ){
1345. # train[i,] <- 1 # train is filled with 1's by default
1346. test[i,] <- 1
1347. } else {
1348. # for taxa with greater occs than k.folds
1349. cv <- occs.df$presence[i]
1350. train[i, (k.folds+1-cv) ] <- 0
1351. test[i, (k.folds+1-cv) ] <- 1
1352. }
1353. } # closing for(i in 1:nrow(occs.df) )
1357. out <- do.call('cbind', list( occs.df, train, test ) )
1359. } # closing if( k.folds > 1 )
1361. # setting all the presence rows to 1
1362. out$presence <- 1
1365. # output
1366. return(out)
1367. }
1370. # function that reduces occurrences to one per grid cell
1372. ### ARGUMENTS ###
1373. ## occ.data = data.frame of the occurrence file
1374. ## rast = raster of the entire background training extent
1375. # there will be issues if any occs are in grid cells with NA values (i.e., outside the extent)
1376. ## name = name of the column that has the taxa names
1377. ## long= column name that has longitude
1378. ## lat = column name that has latitude
1379. ## max.dist = argument from seegSDM. distance is in map units (e.g., degrees) if the raster is projected, otherwise, it is in meters.
1380. # if any occurrences lie JUST BARELY outside the extent, this function will assign them to the nearest gril cell . . .
1381. # . . . inside the extent if the distance from that occurrence to the nearest cell is <= max.dist.
1382. # otherwise, that point is ignored.
1383. # recommended that you remove points outside the extent first. it's easier to clear up that way
1384. ## round.to is how many decimal places to round the coordinates to
1385. # fewer decimal plaes = faster run times
1388. # nearestLand function extracted from seegSDM as that package doesn't appear to exist for more versions (> 3.6.2) of R
1389. # this version of nearestLand is from seegSDM version 0.1-9
1390. seegSDM_nearestLand <- function (points, raster, max_distance) {
1391. nearest <- function(lis, raster) {
1392. neighbours <- matrix(lis[[1]], ncol = 2)
1393. point <- lis[[2]]
1394. land <- !is.na(neighbours[, 2])
1395. if (!any(land)) {
1396. return(c(NA, NA))
1397. } else {
1398. coords <- xyFromCell(raster, neighbours[land, 1])
1399. if (nrow(coords) == 1) {
1400. return(coords[1, ])
1401. }
1402. dists <- sqrt((coords[, 1] - point[1])^2 + (coords[, 2] - point[2])^2)
1403. return(coords[which.min(dists), ])
1404. } # ending else
1405. } # ending nearest
1407. neighbour_list <- extract(raster, points, buffer = max_distance, cellnumbers = TRUE)
1408. neighbour_list <- lapply(1:nrow(points), function(i) {
1409. list(neighbours = neighbour_list[[i]], point = as.numeric(points[i, ]))
1410. })
1411. return(t(sapply(neighbour_list, nearest, raster)))
1412. }
1415. #
1416. 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?
1417. require(raster)
1418. occ.mat <- as.matrix(occ.data[,c(long, lat)])
1419. occ.mat <- round(occ.mat, digits = round.to)
1420. moved <- seegSDM_nearestLand(occ.mat, raster=rast, max_distance=max.dist) # centers occ. points within the grid cell they occur in
1421. moved <- as.data.frame(moved)
1422. moved <- cbind(occ.data[,name], moved) # bind names and the paleo-longitudes
1423. colnames(moved) <- c(name, 'long.thin', 'lat.thin')
1424. moved <- unique(moved) # returns only one occurrence per entity per grid
1425. numbs <- as.numeric(row.names(moved)) # row numbers of the thinned data
1426. out <- occ.data[numbs,] # thinning the original data with the rows of the thinned data
1427. return(out)
1428. }
1431. #################&&&&&&&&&&&&&&&&&&&&&############MAKE MASTER#################
1433. #################### ACER PRES MASTER
1435. # reading in the occurrence file
1436. Pres.occs <- read.csv('AcroporaCervicornisPresent.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
1438. Holo.occs <- read.csv('AcroporaCervicornisHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
1440. Pres.occs <- Pres.occs[,-1]
1441. Holo.occs <- Holo.occs[,-1]
1443. head(Pres.occs)
1444. head(Holo.occs)
1445. identical(colnames(Pres.occs), colnames(Holo.occs))
1447. nrow(Pres.occs)
1448. nrow(Holo.occs)
1450. Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3)
1451. Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3)
1453. #assigning wgs84 projection if not already assigned in layers
1454. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
1456. ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type
1458. ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1459. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present")
1461. Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$')
1462. Prescomp.files
1464. # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas)
1465. #Pres.files <- Pres.files[ c(1,2,3,4,5) ]
1466. Prescomp.rasters <- raster::stack( Prescomp.files)
1468. # optionally assigning a crs if one isn't provided
1469. crs(Prescomp.rasters) <- wgs1984
1470. names(Prescomp.rasters)
1472. # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1473. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED")
1474. Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$')
1475. Holo.files
1477. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1478. Holo.rasters <- raster::stack( Holo.files)
1479. Holo.rasters <- Holo.rasters/100
1481. # optionally assigning a crs if one isn't provided
1482. crs(Holo.rasters) <- wgs1984
1483. names(Holo.rasters)
1485. # re-naming the names of rasters if you want to
1486. names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1487. names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1489. #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
1490. identical( names(Prescomp.rasters), names(Holo.rasters) )
1493. ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
1494. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050")
1495. RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$')
1496. RCP452050.files
1498. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1499. RCP452050.rasters <- raster::stack( RCP452050.files)
1501. # optionally assigning a crs if one isn't provided
1502. crs(RCP452050.rasters) <- wgs1984
1503. names(RCP452050.rasters)
1505. # re-naming the names of turo.rasters if you want to
1506. names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1508. #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
1509. identical( names(Holo.rasters), names(RCP452050.rasters) )
1511. # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
1512. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050")
1513. RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$')
1514. RCP852050.files
1516. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1517. RCP852050.rasters <- raster::stack( RCP852050.files)
1519. # optionally assigning a crs if one isn't provided
1520. crs(RCP852050.rasters) <- wgs1984
1521. names(RCP852050.rasters)
1523. # re-naming the names of turo.rasters if you want to
1524. names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1526. #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
1527. identical( names(Holo.rasters), names(RCP852050.rasters) )
1529. # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1530. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100")
1531. RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$')
1532. RCP452100.files
1534. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1535. RCP452100.rasters <- raster::stack( RCP452100.files)
1537. # optionally assigning a crs if one isn't provided
1538. crs(RCP452100.rasters) <- wgs1984
1539. names(RCP452100.rasters)
1541. names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
1543. #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
1544. identical( names(Holo.rasters), names(RCP452100.rasters) )
1547. # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1548. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100")
1549. RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$')
1550. RCP852100.files
1552. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1553. RCP852100.rasters <- raster::stack( RCP852100.files)
1555. # optionally assigning a crs if one isn't provided
1556. crs(RCP852100.rasters) <- wgs1984
1557. names(RCP852100.rasters)
1559. names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
1561. #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
1562. identical( names(Holo.rasters), names(RCP852100.rasters) )
1564. # making the background files
1565. Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
1567. Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
1569. RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean')
1571. RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean')
1573. RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean')
1575. RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean')
1578. #spatial thinning
1579. Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat")
1581. Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat")
1583. #how many occurrences are there - compare to original before spatial thinning to see if less
1584. nrow(Pres.occs.thin)
1585. nrow(Holo.occs.thin)
1587. #look at new spatially thinned dataframe
1588. View(Pres.occs.thin)
1590. # making the occurrence files
1591. # c(8,2,1) = (GENUS, Longitude, Latitude)
1592. 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
1593. x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
1595. Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)],
1596. x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
1599. # checking to see if everything has identical column names for your different extents
1600. identical( colnames(Prescomp.back), colnames(Holo.back) )
1602. identical( colnames(Pres.occs.df), colnames(Holo.occs.df) )
1605. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/Master_files_acerv_add_mask")
1608. masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back,
1609. RCP452050.back,RCP852050.back,
1610. RCP452100.back, RCP852100.back
1611. ) )
1613. write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE)
1615. masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) )
1617. View(masterall.occs)
1618. write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE)
1620. masterpres.occs <- Pres.occs.df
1621. write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE)
1623. masterholo.occs <- Holo.occs.df
1624. write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE)
1627. #################################### ACER HOLO MASTER same code but with area masked for cuba jamaica and hispanola
1631. #################################### APAL PRES MASTER
1634. # reading in the occurrence file
1635. Pres.occs <- read.csv('AcroporaPalmataPresent2.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
1637. Holo.occs <- read.csv('AcroporaPalmataHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
1639. Pres.occs <- Pres.occs[,-1]
1640. Holo.occs <- Holo.occs[,-1]
1642. head(Pres.occs)
1643. head(Holo.occs)
1644. identical(colnames(Pres.occs), colnames(Holo.occs))
1646. nrow(Pres.occs)
1647. nrow(Holo.occs)
1649. Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3)
1650. Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3)
1653. ################################################################################################################################################
1655. #assignin wgs84 projection if not already assigned in layers
1656. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
1658. ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type
1660. ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1662. Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$')
1663. Prescomp.files
1664. # re-ordering any files
1666. # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas)
1667. #Pres.files <- Pres.files[ c(1,2,3,4,5) ]
1668. Prescomp.rasters <- raster::stack( Prescomp.files)
1670. # optionally assigning a crs if one isn't provided
1671. crs(Prescomp.rasters) <- wgs1984
1672. names(Prescomp.rasters)
1674. # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1675. Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam', pattern='\\.tif$')
1676. Holo.files
1678. #reorder
1679. Holo.files <- Holo.files[ c(1,3,2,4) ]
1680. Holo.files
1682. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1683. Holo.rasters <- raster::stack( Holo.files)
1684. Holo.rasters <- Holo.rasters/100
1686. # optionally assigning a crs if one isn't provided
1687. crs(Holo.rasters) <- wgs1984
1688. names(Holo.rasters)
1690. names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1691. names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1693. # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1694. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED")
1695. Holo2.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$')
1696. Holo2.files
1700. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1701. Holo2.rasters <- raster::stack( Holo2.files)
1702. Holo2.rasters <- Holo2.rasters/100
1704. # optionally assigning a crs if one isn't provided
1705. crs(Holo2.rasters) <- wgs1984
1706. names(Holo2.rasters)
1709. # re-naming the names of rasters if you want to
1710. names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1711. names(Holo2.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1714. #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
1715. identical( names(Prescomp.rasters), names(Holo.rasters) )
1717. #
1718. ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
1719. RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$')
1720. RCP452050.files
1722. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1723. RCP452050.rasters <- raster::stack( RCP452050.files)
1725. # optionally assigning a crs if one isn't provided
1726. crs(RCP452050.rasters) <- wgs1984
1727. names(RCP452050.rasters)
1729. # re-naming the names of turo.rasters if you want to
1730. names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1732. #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
1733. identical( names(Holo.rasters), names(RCP452050.rasters) )
1735. #
1736. # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
1737. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050")
1738. RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$')
1739. RCP852050.files
1741. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1742. RCP852050.rasters <- raster::stack( RCP852050.files)
1744. # optionally assigning a crs if one isn't provided
1745. crs(RCP852050.rasters) <- wgs1984
1746. names(RCP852050.rasters)
1748. # re-naming the names of turo.rasters if you want to
1749. names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1751. #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
1752. identical( names(Holo.rasters), names(RCP852050.rasters) )
1754. # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1755. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100")
1756. RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$')
1757. RCP452100.files
1759. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1760. RCP452100.rasters <- raster::stack( RCP452100.files)
1762. # optionally assigning a crs if one isn't provided
1763. crs(RCP452100.rasters) <- wgs1984
1764. names(RCP452100.rasters)
1766. names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
1768. #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
1769. identical( names(Holo.rasters), names(RCP452100.rasters) )
1771. # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1772. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100")
1773. RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$')
1774. RCP852100.files
1776. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1777. RCP852100.rasters <- raster::stack( RCP852100.files)
1779. # optionally assigning a crs if one isn't provided
1780. crs(RCP852100.rasters) <- wgs1984
1781. names(RCP852100.rasters)
1783. names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
1785. #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
1786. identical( names(Holo.rasters), names(RCP852100.rasters) )
1787. #
1788. #
1789. # making the background files
1790. Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
1792. Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
1794. Holo2.back <- make.background.df(r.stack=Holo2.rasters, x.col='long', y.col='lat', time.bin='Holocene2', region='Caribbean')
1795. #
1796. RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean')
1797. #
1798. RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean')
1799. #
1800. RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean')
1801. #
1802. RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean')
1805. #spatial thinning
1806. Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat")
1808. Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat")
1811. #how many occurrences are there - compare to original before spatial thinning to see if less
1812. nrow(Pres.occs.thin)
1813. nrow(Holo.occs.thin)
1816. #look at new spatially thinned dataframe
1817. View(Pres.occs.thin)
1819. # making the occurrence files
1820. # c(8,2,1) = (GENUS, Longitude, Latitude)
1821. 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
1822. x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
1824. Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)],
1825. x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
1828. # checking to see if everything has identical column names for your different extents
1829. identical( colnames(Prescomp.back), colnames(Holo.back) )
1832. identical( colnames(Pres.occs.df), colnames(Holo.occs.df) )
1835. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/Master_files_apal_add_mask")
1837. # merging together the files
1839. masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back, Holo2.back,
1840. RCP452050.back,RCP852050.back,
1841. RCP452100.back, RCP852100.back
1842. ) )
1844. write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE)
1847. masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) )
1849. View(masterall.occs)
1850. write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE)
1852. masterpres.occs <- Pres.occs.df
1853. write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE)
1855. masterholo.occs <- Holo.occs.df
1856. write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE)
1858. ########################### APAL HOLO MASTER same but with masked background
1861. #################### CNAT PRES MASTER
1862. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/Cnatans")
1864. # reading in the occurrence file
1865. Pres.occs <- read.csv('CnatansPresent.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
1867. Holo.occs <- read.csv('CnatansHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
1869. Pres.occs <- Pres.occs[,-1]
1870. Holo.occs <- Holo.occs[,-1]
1872. head(Pres.occs)
1873. head(Holo.occs)
1874. identical(colnames(Pres.occs), colnames(Holo.occs))
1876. nrow(Pres.occs)
1877. nrow(Holo.occs)
1879. Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3)
1880. Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3)
1883. ################################################################################################################################################
1885. #assignin wgs84 projection if not already assigned in layers
1886. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
1888. ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type
1890. ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1891. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present")
1893. Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$')
1894. Prescomp.files
1895. Prescomp.rasters <- raster::stack( Prescomp.files)
1898. # optionally assigning a crs if one isn't provided
1899. crs(Prescomp.rasters) <- wgs1984
1900. names(Prescomp.rasters)
1902. # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1903. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED")
1904. Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$')
1905. Holo.files
1907. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1908. Holo.rasters <- raster::stack( Holo.files)
1909. #scaling the holocene raster
1910. Holo.rasters <- Holo.rasters/100
1912. # optionally assigning a crs if one isn't provided
1913. crs(Holo.rasters) <- wgs1984
1914. names(Holo.rasters)
1916. # re-naming the names of rasters if you want to
1917. names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1918. names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1921. #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
1922. identical( names(Prescomp.rasters), names(Holo.rasters) )
1924. ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
1925. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050")
1926. RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$')
1927. RCP452050.files
1929. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1930. RCP452050.rasters <- raster::stack( RCP452050.files)
1932. # optionally assigning a crs if one isn't provided
1933. crs(RCP452050.rasters) <- wgs1984
1934. names(RCP452050.rasters)
1936. # re-naming the names of turo.rasters if you want to
1937. names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1939. #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
1940. identical( names(Holo.rasters), names(RCP452050.rasters) )
1941. # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
1942. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050")
1943. RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$')
1944. RCP852050.files
1946. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1947. RCP852050.rasters <- raster::stack( RCP852050.files)
1949. # optionally assigning a crs if one isn't provided
1950. crs(RCP852050.rasters) <- wgs1984
1951. names(RCP852050.rasters)
1953. # re-naming the names of turo.rasters if you want to
1954. names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
1956. #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
1957. identical( names(Holo.rasters), names(RCP852050.rasters) )
1959. #
1960. # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1961. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100")
1962. RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$')
1963. RCP452100.files
1965. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1966. RCP452100.rasters <- raster::stack( RCP452100.files)
1968. # optionally assigning a crs if one isn't provided
1969. crs(RCP452100.rasters) <- wgs1984
1970. names(RCP452100.rasters)
1972. names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
1974. #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
1975. identical( names(Holo.rasters), names(RCP452100.rasters) )
1977. #
1978. # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1979. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100")
1980. RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$')
1981. RCP852100.files
1983. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
1984. RCP852100.rasters <- raster::stack( RCP852100.files)
1986. # optionally assigning a crs if one isn't provided
1987. crs(RCP852100.rasters) <- wgs1984
1988. names(RCP852100.rasters)
1990. names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
1992. #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
1993. identical( names(Holo.rasters), names(RCP852100.rasters) )
1994. #
1995. #
1996. # making the background files
1997. Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
1999. Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
2000. RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean')
2001. RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean')
2002. RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean')
2003. RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean')
2004. #spatial thinning
2005. Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat")
2007. Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat")
2009. #how many occurrences are there - compare to original before spatial thinning to see if less
2010. nrow(Pres.occs.thin)
2011. nrow(Holo.occs.thin)
2013. #look at new spatially thinned dataframe
2014. View(Pres.occs.thin)
2016. # making the occurrence files
2017. # c(8,2,1) = (GENUS, Longitude, Latitude)
2018. 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
2019. x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
2021. Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)],
2022. x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
2026. # checking to see if everything has identical column names for your different extents
2027. identical( colnames(Prescomp.back), colnames(Holo.back) )
2029. identical( colnames(Pres.occs.df), colnames(Holo.occs.df) )
2031. identical( colnames(Prescomp.back), colnames(Holo.occs.df) )
2034. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/Master_files_cnat_add_mask")
2036. # merging together the files
2037. masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back,
2038. RCP452050.back,RCP852050.back,
2039. RCP452100.back, RCP852100.back) )
2041. write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE)
2044. masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) )#LGM.occs
2046. View(masterall.occs)
2047. write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE)
2049. masterpres.occs <- Pres.occs.df
2050. write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE)
2052. masterholo.occs <- Holo.occs.df
2053. write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE)
2057. #################### CNAT HOLO MASTER same with masked extent
2061. #################### PAST PRES MASTER
2063. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/Pasteroides")
2065. # reading in the occurrence file
2066. Pres.occs <- read.csv('PoritesAsteroidesPresent.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
2068. Holo.occs <- read.csv('PoritesAsteroidesHolocene.csv', header = TRUE, stringsAsFactors=FALSE, fileEncoding="latin1")
2070. Pres.occs <- Pres.occs[,-1]
2071. Holo.occs <- Holo.occs[,-1]
2074. head(Pres.occs)
2075. head(Holo.occs)
2076. identical(colnames(Pres.occs), colnames(Holo.occs))
2078. nrow(Pres.occs)
2079. nrow(Holo.occs)
2081. Pres.occs[,1:2] <- round(Pres.occs[,1:2], 3)
2082. Holo.occs[,1:2] <- round(Holo.occs[,1:2], 3)
2085. ################################################################################################################################################
2087. #assignin wgs84 projection if not already assigned in layers
2088. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
2090. ### reading in the raster files (if they are geotiff, then change below to '\\.tif$'); all layers must be the same file type
2092. ## present+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2093. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present")
2095. Prescomp.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Bio-Oracle Future/MASKED/present', pattern='\\.tif$')
2096. Prescomp.files
2097. # re-ordering any files
2100. # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas)
2101. #Pres.files <- Pres.files[ c(1,2,3,4,5) ]
2102. Prescomp.rasters <- raster::stack( Prescomp.files)
2105. # optionally assigning a crs if one isn't provided
2106. crs(Prescomp.rasters) <- wgs1984
2107. names(Prescomp.rasters)
2109. # ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2110. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED")
2111. Holo.files <- list.files(path='/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/HOLOMASKED', pattern='\\.tif$')
2112. Holo.files
2114. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
2115. Holo.rasters <- raster::stack( Holo.files)
2116. Holo.rasters <- Holo.rasters/100
2118. # optionally assigning a crs if one isn't provided
2119. crs(Holo.rasters) <- wgs1984
2120. names(Holo.rasters)
2122. # re-naming the names of rasters if you want to
2123. names(Prescomp.rasters) <- c( "maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
2124. names(Holo.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
2127. #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
2128. identical( names(Prescomp.rasters), names(Holo.rasters) )
2131. ######################future 2050 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
2132. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050")
2133. RCP452050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452050", pattern='\\.tif$')
2134. RCP452050.files
2136. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
2137. RCP452050.rasters <- raster::stack( RCP452050.files)
2139. # optionally assigning a crs if one isn't provided
2140. crs(RCP452050.rasters) <- wgs1984
2141. names(RCP452050.rasters)
2143. # re-naming the names of turo.rasters if you want to
2144. names(RCP452050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
2146. #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
2147. identical( names(Holo.rasters), names(RCP452050.rasters) )
2150. # #############################future 2050 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
2151. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050")
2152. RCP852050.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852050", pattern='\\.tif$')
2153. RCP852050.files
2155. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
2156. RCP852050.rasters <- raster::stack( RCP852050.files)
2158. # optionally assigning a crs if one isn't provided
2159. crs(RCP852050.rasters) <- wgs1984
2160. names(RCP852050.rasters)
2162. # re-naming the names of turo.rasters if you want to
2163. names(RCP852050.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp" )
2165. #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
2166. identical( names(Holo.rasters), names(RCP852050.rasters) )
2169. # ###############future 2100 RCP45+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2170. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100")
2171. RCP452100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP452100", pattern='\\.tif$')
2172. RCP452100.files
2174. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
2175. RCP452100.rasters <- raster::stack( RCP452100.files)
2177. # optionally assigning a crs if one isn't provided
2178. crs(RCP452100.rasters) <- wgs1984
2179. names(RCP452100.rasters)
2181. names(RCP452100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
2183. #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
2184. identical( names(Holo.rasters), names(RCP452100.rasters) )
2187. # #################future 2100 RCP85+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2188. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100")
2189. RCP852100.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle future/MASKED/RCP852100", pattern='\\.tif$')
2190. RCP852100.files
2192. #when loading things as a raster::stack, must have the same file extent - nrows, ncol, ncell, nlayers, pixel size
2193. RCP852100.rasters <- raster::stack( RCP852100.files)
2195. # optionally assigning a crs if one isn't provided
2196. crs(RCP852100.rasters) <- wgs1984
2197. names(RCP852100.rasters)
2199. names(RCP852100.rasters) <- c("maxsalinity", "maxtemp", "rangesalinity", "rangetemp")
2201. #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
2202. identical( names(Holo.rasters), names(RCP852100.rasters) )
2203. #
2204. #
2205. # making the background files
2206. Prescomp.back <- make.background.df(r.stack=Prescomp.rasters, x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
2208. Holo.back <- make.background.df(r.stack=Holo.rasters, x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
2210. RCP452050.back <- make.background.df(r.stack=RCP452050.rasters, x.col='long', y.col='lat', time.bin='RCP452050', region='Caribbean')
2212. RCP852050.back <- make.background.df(r.stack=RCP852050.rasters, x.col='long', y.col='lat', time.bin='RCP852050', region='Caribbean')
2214. RCP452100.back <- make.background.df(r.stack=RCP452100.rasters, x.col='long', y.col='lat', time.bin='RCP452100', region='Caribbean')
2216. RCP852100.back <- make.background.df(r.stack=RCP852100.rasters, x.col='long', y.col='lat', time.bin='RCP852100', region='Caribbean')
2219. #spatial thinning
2220. Pres.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Pres.occs, rast = Prescomp.rasters[[2]], name = "name", long = "long", lat = "lat")
2222. Holo.occs.thin <- one.occ.per.grid.cell.no.SEEG(occ.data = Holo.occs, rast = Holo.rasters[[2]], name = "name", long = "long", lat = "lat")
2226. #how many occurrences are there - compare to original before spatial thinning to see if less
2227. nrow(Pres.occs.thin)
2228. nrow(Holo.occs.thin)
2231. #look at new spatially thinned dataframe
2232. View(Pres.occs.thin)
2234. # making the occurrence files
2235. # c(8,2,1) = (GENUS, Longitude, Latitude)
2236. 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
2237. x.col='long', y.col='lat', time.bin='Present', region='Caribbean')
2239. Holo.occs.df <- make.occurrence.df(r.stack=Holo.rasters, taxa.df=Holo.occs.thin[,c(3,1,2)],
2240. x.col='long', y.col='lat', time.bin='Holocene', region='Caribbean')
2243. # checking to see if everything has identical column names for your different extents
2244. identical( colnames(Prescomp.back), colnames(Holo.back) )
2246. identical( colnames(Pres.occs.df), colnames(Holo.occs.df) )
2250. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/Master_files_past_add_mask")
2252. # merging together the files
2254. masterall.background <- do.call('rbind', list(Prescomp.back, Holo.back,
2255. RCP452050.back,RCP852050.back,
2256. RCP452100.back, RCP852100.back) )
2258. write.csv(masterall.background, 'masterall.background.csv', row.names = FALSE)
2260. masterall.occs <- do.call('rbind', list(Pres.occs, Holo.occs) )
2262. View(masterall.occs)
2263. write.csv(masterall.occs, 'masterall.occs.csv', row.names = FALSE)
2265. masterpres.occs <- Pres.occs.df
2266. write.csv(masterpres.occs, 'masterpres.occs.csv', row.names = FALSE)
2268. masterholo.occs <- Holo.occs.df
2269. write.csv(masterholo.occs, 'masterholo.occs.csv', row.names = FALSE)
2280. #################### PAST HOLO MASTER same with masked extent
2290. #################&&&&&&&&&&&&&&&&&&&&&############ECOSPAT#################
2293. ####################################### ACER ECOSPAT ###########################
2294. #### MODEL COMPARISON ####
2296. #need cleaned occurrence data and environmental data rasters
2298. ## Read in shapefiles of cleaned occ data ##
2299. ########################occurrence data######################
2300. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/Master_files_acerv_add_mask")
2301. AcervPres = read.csv("masterpres.occs.csv")
2302. AcervPres = as(AcervPres,'data.frame')
2304. AcervHolo = read.csv("masterholo.occs.csv")
2305. AcervHolo = as(AcervHolo,'data.frame')
2307. #convert to dataframe for use in ecospat
2309. ###################environmental###################
2310. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present")
2312. # Upload basic rasters for First interval (Present)
2313. Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$')
2314. Prescomp.files
2316. maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif")
2317. maxsalinitypres
2319. maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif")
2320. maxtemppres
2322. rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif")
2323. rangesalinitypres
2325. rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif")
2326. rangetemppres
2329. # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas)
2330. Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres))
2333. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
2334. # optionally assigning a crs if one isn't provided
2335. crs(Prescomp.rasters) <- wgs1984
2336. names(Prescomp.rasters)
2337. names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
2341. #Upload basic rasters for Second interval (Holocene)
2342. ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2343. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam")
2344. Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$')
2347. #scaling holocene layers
2348. maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif")
2349. maxsalinityholo2 <- maxsalinityholo/100
2351. maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif")
2352. maxtempholo2 <- maxtempholo/100
2354. rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif")
2355. rangesalinityholo2 <- rangesalinityholo/100
2357. rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif")
2358. rangetempholo2 <- rangetempholo/100
2360. Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2)
2362. # optionally assigning a crs if one isn't provided
2363. crs(Holo.rasters) <- wgs1984
2364. names(Holo.rasters)
2366. names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
2369. ### Ecospat Analysis of Niche Equivalency and Similary
2371. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2373. # Set working directory for saving plots and test results
2374. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/maskedtwice/Ecospat_acerv_10_2_24")
2376. combine.lat = c(AcervPres$lat, AcervHolo$lat)
2377. combine.lon = c(AcervPres$long, AcervHolo$long)
2379. #I need to turn the points into a spatvector (
2380. library(terra)
2381. ptspres <- vect(cbind(AcervPres$long, AcervPres$lat))
2382. ptsholo <- vect(cbind(AcervHolo$long, AcervHolo$lat))
2385. # get random sample - background point locations -
2386. rasterobject <- rast(Prescomp.rasters)
2387. bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
2388. return.type = "points", n = 500)
2390. rasterobject <- rast(Holo.rasters)
2391. bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
2392. return.type = "points", n = 500)
2394. # get environmental data from occ points
2395. extractpres = na.omit(cbind(AcervPres[2:3], raster::extract(Prescomp.rasters, AcervPres[2:3]), rep(1, nrow(AcervPres))))
2397. extractholo = na.omit(cbind(AcervHolo[2:3], raster::extract(Holo.rasters, AcervHolo[2:3]), rep(1, nrow(AcervHolo))))
2399. #change name of last column to occ (represents presence = 1)
2400. colnames(extractpres)[ncol(extractpres)] = 'occ'
2401. colnames(extractholo)[ncol(extractholo)] = 'occ'
2403. # get environmental data from bg points
2404. extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres)))))
2406. extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo)))))
2409. # change name of last column to occ (represents absence = 0)
2410. colnames(extbgpres)[ncol(extbgpres)] = 'occ'
2411. colnames(extbgholo)[ncol(extbgholo)] = 'occ'
2414. colnames(extbgpres)[1] = 'long'
2415. colnames(extbgpres)[2] = 'lat'
2417. colnames(extbgholo)[1] = 'long'
2418. colnames(extbgholo)[2] = 'lat'
2420. # merge data from occ and bg
2421. datPresAcerv = rbind(extractpres, extbgpres)
2423. datHoloAcerv = rbind(extractholo, extbgholo)
2426. #principle component analysis
2427. pca.env <- dudi.pca(
2428. rbind(datHoloAcerv, datPresAcerv)[,3:6],
2429. scannf=FALSE,
2430. nf=2
2431. )
2433. # look at variable contribution
2434. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
2436. # Save pdf of PCA plot
2437. pdf('ecospat_Acervpresholo.pdf')
2438. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
2440. scores.globclim<-pca.env$li # PCA scores for the whole study area
2441. scores.globclim<-pca.env$li # PCA scores for the whole study area (all points)
2443. scores.sppres <- suprow(pca.env,
2444. extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution
2446. scores.spholo <- suprow(pca.env,
2447. extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution
2449. scores.climholo <- suprow(pca.env,datHoloAcerv[,3:6])$li # PCA scores for the first interval study area
2450. scores.climpres <- suprow(pca.env,datPresAcerv[,3:6])$li # PCA scores for the second interval study area
2453. # density distribtuion for first interval
2454. grid.Acervholo <- ecospat.grid.clim.dyn(
2455. glob = scores.globclim,
2456. glob1 = scores.climholo,
2457. sp = scores.spholo,
2458. R = 100,
2459. th.sp = 0
2460. )
2462. # density distribution for second interval
2463. grid.Acervpres <- ecospat.grid.clim.dyn(
2464. glob = scores.globclim,
2465. glob1 = scores.climpres,
2466. sp = scores.sppres,
2467. R = 100,
2468. th.sp = 0
2469. )
2471. histpres <- hist(grid.Acervpres$glob1)
2472. histholo <- hist(grid.Acervholo$glob1)
2473. plot( histpres, col=rgb(0,0,1,1/4)) # first histogram
2474. plot( histholo, col=rgb(1,0,0,1/4), add=T)
2476. # Schoener's D metric and I metric
2477. D.overlap <- ecospat.niche.overlap (grid.Acervholo, grid.Acervpres, cor=T)
2478. D.overlap
2479. write.csv(D.overlap,"presAcerv_holoAcerv_I_D.csv",row.names=FALSE)
2481. ## Niche Equivalency Test
2482. ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D
2483. ## p.I: pvalue of the test on I
2484. ## Test for greater or lower equivalency
2485. eq.testgr <- ecospat.niche.equivalency.test(grid.Acervholo, grid.Acervpres,
2486. rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs
2487. eq.testlw <- ecospat.niche.equivalency.test(grid.Acervholo, grid.Acervpres,
2488. rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs
2491. nichdynindex <- ecospat.niche.dyn.index(grid.Acervholo, grid.Acervpres, intersection = NA)
2492. #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w
2493. nichdynindex
2496. # write p values of equivalency test
2497. 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)
2498. write.csv(p_EQ_DI,"glycim_EQ_TestAcervpresholo.csv",row.names=FALSE)
2500. p_EQ_DI
2503. ## Test for greater (niche conservatism) or lower (niche divergence) similarity
2504. sim.testgr <- ecospat.niche.similarity.test(grid.Acervholo, grid.Acervpres,
2505. rep=1000, overlap.alternative = "higher",
2506. rand.type=2,ncores=4)
2507. sim.testlw <- ecospat.niche.similarity.test(grid.Acervholo, grid.Acervpres,
2508. rep=1000, overlap.alternative = "lower",
2509. rand.type=2,ncores=4)
2511. # write p values of similarity test
2512. 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)
2513. write.csv(p_SIM_DI,"presholoAcerv_SIM_Test.csv",row.names=FALSE)
2515. p_SIM_DI
2517. # Plot test distributions
2518. ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency")
2519. ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency")
2520. ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity")
2521. ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity")
2524. ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals
2525. ecospat.plot.niche (grid.Acervpres, title='Acervpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE)
2526. ecospat.plot.niche (grid.Acervholo, title='Acervholo', name.axis1='PC1', name.axis2='PC2')
2529. ## Plot niche overlap
2530. ecospat.plot.niche.dyn(z1 = grid.Acervholo, z2 = grid.Acervpres, quant=0.25, interest=2,
2531. title= "A. cerv Holo to Pres", name.axis1="PC1",
2532. name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = '#74add1', colZ1 =
2533. "#313695", colZ2 = "#abd9d9", transparency = 40)
2534. ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white')
2535. 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)
2537. # Save pdf of plots
2538. ### Look at niche expantion, stability, and unfilling
2539. # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection
2540. dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Acervholo, z2 = grid.Acervpres, intersection=NA)
2541. dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Acervholo, z2 = grid.Acervpres, intersection=0)
2543. # write csv of niche dynmaic percentages
2544. Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w,
2545. "Overlapping"=dynam_overlap$dynamic.index.w)
2547. Niche_Ex_St
2548. write.csv(Niche_Ex_St,"GlycymDynamicsAcervholopres.csv",row.names=TRUE)
2550. ############################################ APAL ECOSPAT ###################
2552. ## Read in shapefiles of cleaned occ data ##
2553. ########################occurrence data######################
2554. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/Master_files_apal_add_mask")
2555. ApalPres = read.csv("masterpres.occs.csv")
2556. ApalPres = as(ApalPres,'data.frame')
2558. ApalHolo = read.csv("masterholo.occs.csv")
2559. ApalHolo = as(ApalHolo,'data.frame')
2560. #convert to dataframe for use in ecospat
2563. ###################environmental###################
2565. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present")
2567. # Upload basic rasters for First interval (Present)
2568. Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$')
2569. Prescomp.files
2571. maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif")
2572. maxsalinitypres
2574. maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif")
2575. maxtemppres
2577. rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif")
2578. rangesalinitypres
2580. rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif")
2581. rangetemppres
2584. # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas)
2585. Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres))
2588. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
2589. # optionally assigning a crs if one isn't provided
2590. crs(Prescomp.rasters) <- wgs1984
2591. names(Prescomp.rasters)
2592. names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
2596. #Upload basic rasters for Second interval (Holocene)
2597. ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2598. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam")
2599. Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$')
2600. #Holo.files
2602. # scaling holocene layers
2603. maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif")
2604. maxsalinityholo2 <- maxsalinityholo/100
2606. maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif")
2607. maxtempholo2 <- maxtempholo/100
2609. rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif")
2610. rangesalinityholo2 <- rangesalinityholo/100
2612. rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif")
2613. rangetempholo2 <- rangetempholo/100
2615. Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2)
2617. # optionally assigning a crs if one isn't provided
2618. crs(Holo.rasters) <- wgs1984
2619. names(Holo.rasters)
2621. names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
2624. ### Ecospat Analysis of Niche Equivalency and Similary
2626. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2628. # Set working directory for saving plots and test results
2629. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/Ecospat_apal_10_2_24")
2631. combine.lat = c(ApalPres$lat, ApalHolo$lat) #, ApalLGM$lat)
2632. combine.lon = c(ApalPres$long, ApalHolo$long) #, ApalLGM$long)
2635. #I need to turn the points into a spatvector (
2636. ptspres <- vect(cbind(ApalPres$long, ApalPres$lat))
2637. ptsholo <- vect(cbind(ApalHolo$long, ApalHolo$lat))
2640. # get random sample - background point locations -
2641. rasterobject <- rast(Prescomp.rasters)
2642. bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
2643. return.type = "points", n = 500)
2645. rasterobject <- rast(Holo.rasters)
2646. bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
2647. return.type = "points", n = 500)
2650. # get environmental data from occ points
2651. extractpres = na.omit(cbind(ApalPres[2:3], raster::extract(Prescomp.rasters, ApalPres[2:3]), rep(1, nrow(ApalPres))))
2653. extractholo = na.omit(cbind(ApalHolo[2:3], raster::extract(Holo.rasters, ApalHolo[2:3]), rep(1, nrow(ApalHolo))))
2656. #change name of last column to occ (represents presence = 1)
2657. colnames(extractpres)[ncol(extractpres)] = 'occ'
2658. colnames(extractholo)[ncol(extractholo)] = 'occ'
2660. # get environmental data from bg points
2661. #cant do stack-- need to do all???
2662. extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres)))))
2664. extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo)))))
2666. # change name of last column to occ (represents absence = 0)
2667. colnames(extbgpres)[ncol(extbgpres)] = 'occ'
2668. colnames(extbgholo)[ncol(extbgholo)] = 'occ'
2670. colnames(extbgpres)[1] = 'long'
2671. colnames(extbgpres)[2] = 'lat'
2673. colnames(extbgholo)[1] = 'long'
2674. colnames(extbgholo)[2] = 'lat'
2676. # merge data from occ and bg
2677. datPresApal = rbind(extractpres, extbgpres)
2679. datHoloApal = rbind(extractholo, extbgholo)
2682. #principle component analysis
2683. pca.env <- dudi.pca(
2684. rbind(datHoloApal, datPresApal)[,3:6],
2685. scannf=FALSE,
2686. nf=2
2687. )
2689. # look at variable contribution
2690. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
2692. # Save pdf of PCA plot
2693. pdf('ecospat_Apalpresholo.pdf')
2694. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
2696. scores.globclim<-pca.env$li # PCA scores for the whole study area
2697. scores.globclim<-pca.env$li # PCA scores for the whole study area (all points)
2699. scores.sppres <- suprow(pca.env,
2700. extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution
2702. scores.spholo <- suprow(pca.env,
2703. extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution
2705. scores.climholo <- suprow(pca.env,datHoloApal[,3:6])$li # PCA scores for the first interval study area
2706. scores.climpres <- suprow(pca.env,datPresApal[,3:6])$li # PCA scores for the second interval study area
2708. # density distribtuion for first interval
2709. grid.Apalholo <- ecospat.grid.clim.dyn(
2710. glob = scores.globclim,
2711. glob1 = scores.climholo,
2712. sp = scores.spholo,
2713. R = 100,
2714. th.sp = 0
2715. )
2717. # density distribution for second interval
2718. grid.Apalpres <- ecospat.grid.clim.dyn(
2719. glob = scores.globclim,
2720. glob1 = scores.climpres,
2721. sp = scores.sppres,
2722. R = 100,
2723. th.sp = 0
2724. )
2726. histpres <- hist(grid.Apalpres$glob1)
2727. histholo <- hist(grid.Apalholo$glob1)
2728. plot( histpres, col=rgb(0,0,1,1/4)) # first histogram
2729. plot( histholo, col=rgb(1,0,0,1/4), add=T)
2731. # Schoener's D metric and I metric
2732. D.overlap <- ecospat.niche.overlap (grid.Apalholo, grid.Apalpres, cor=T)
2733. D.overlap
2734. write.csv(D.overlap,"presApal_holoApal_I_D.csv",row.names=FALSE)
2736. ## Niche Equivalency Test
2737. ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D
2738. ## p.I: pvalue of the test on I
2739. ## Test for greater or lower equivalency
2740. eq.testgr <- ecospat.niche.equivalency.test(grid.Apalholo, grid.Apalpres,
2741. rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs
2742. eq.testlw <- ecospat.niche.equivalency.test(grid.Apalholo, grid.Apalpres,
2743. rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs
2746. nichdynindex <- ecospat.niche.dyn.index(grid.Apalholo, grid.Apalpres, intersection = NA)
2747. #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w
2748. nichdynindex
2750. # write p values of equivalency test
2751. 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)
2752. write.csv(p_EQ_DI,"glycim_EQ_TestApalpresholo.csv",row.names=FALSE)
2754. p_EQ_DI
2756. ## Test for greater (niche conservatism) or lower (niche divergence) similarity
2757. sim.testgr <- ecospat.niche.similarity.test(grid.Apalholo, grid.Apalpres,
2758. rep=1000, overlap.alternative = "higher",
2759. rand.type=2,ncores=4)
2760. sim.testlw <- ecospat.niche.similarity.test(grid.Apalholo, grid.Apalpres,
2761. rep=1000, overlap.alternative = "lower",
2762. rand.type=2,ncores=4)
2764. # write p values of similarity test
2765. 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)
2766. write.csv(p_SIM_DI,"presholoApal_SIM_Test.csv",row.names=FALSE)
2768. p_SIM_DI
2769. #p.D_GR p.I_GR p.D_LW p.I_LW
2770. #[1,] 0.04595405 0.03996004 0.9480519 0.9480519
2771. #dev.off()
2772. # Plot test distributions
2773. ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency")
2774. ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency")
2775. ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity")
2776. ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity")
2778. ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals
2779. ecospat.plot.niche (grid.Apalpres, title='Apalpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE)
2780. ecospat.plot.niche (grid.Apalholo, title='Apalholo', name.axis1='PC1', name.axis2='PC2')
2783. ## Plot niche overlap
2784. ecospat.plot.niche.dyn(z1 = grid.Apalholo, z2 = grid.Apalpres, quant=0.25, interest=2,
2785. title= "A. pal Holo to Pres", name.axis1="PC1",
2786. name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = 'black', colZ1 =
2787. "#313695", colZ2 = "#abd9d9", transparency = 40)
2788. ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white')
2789. 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)
2791. #___________________________
2793. # Save pdf of plots
2794. ### Look at niche expantion, stability, and unfilling
2795. # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection
2796. dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Apalholo, z2 = grid.Apalpres, intersection=NA)
2797. dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Apalholo, z2 = grid.Apalpres, intersection=0)
2799. # write csv of niche dynmaic percentages
2800. Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w,
2801. "Overlapping"=dynam_overlap$dynamic.index.w)
2803. Niche_Ex_St
2804. write.csv(Niche_Ex_St,"GlycymDynamicsApalholopres.csv",row.names=TRUE)
2806. ########################## CNAT ECOSPAT #############
2808. #need cleaned occurrence data and environmental data rasters
2810. ## Read in shapefiles of cleaned occ data ##
2811. ########################occurrence data######################
2812. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/Master_files_cnat_add_mask")
2813. CnatPres = read.csv("masterpres.occs.csv")
2814. CnatPres = as(CnatPres,'data.frame')
2816. CnatHolo = read.csv("masterholo.occs.csv")
2817. CnatHolo = as(CnatHolo,'data.frame')
2818. #convert to dataframe for use in ecospat
2820. ##################environmental###################
2821. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present")
2823. # Upload basic rasters for First interval (Present)
2824. Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$')
2825. Prescomp.files
2828. maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif")
2829. maxsalinitypres
2831. maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif")
2832. maxtemppres
2834. rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif")
2835. rangesalinitypres
2837. rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif")
2838. rangetemppres
2841. # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas)
2842. Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres))
2845. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
2846. # optionally assigning a crs if one isn't provided
2847. crs(Prescomp.rasters) <- wgs1984
2848. names(Prescomp.rasters)
2849. names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
2852. #Upload basic rasters for Second interval (Holocene)
2853. ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2854. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam")
2855. Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$')
2856. #Holo.files
2859. maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif")
2860. maxsalinityholo2 <- maxsalinityholo/100
2862. maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif")
2863. maxtempholo2 <- maxtempholo/100
2865. rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif")
2866. rangesalinityholo2 <- rangesalinityholo/100
2868. rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif")
2869. rangetempholo2 <- rangetempholo/100
2871. Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2)
2873. # optionally assigning a crs if one isn't provided
2874. crs(Holo.rasters) <- wgs1984
2875. names(Holo.rasters)
2877. names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
2881. ### Ecospat Analysis of Niche Equivalency and Similary
2883. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2885. # Set working directory for saving plots and test results
2886. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/maskedtwice/Ecospat_cnat_10_2_24")
2888. combine.lat = c(CnatPres$lat, CnatHolo$lat) #, CnatLGM$lat)
2889. combine.lon = c(CnatPres$long, CnatHolo$long) #, CnatLGM$long)
2891. #combine.lat = c(glyMaa$coords.x1, glyDan$coords.x1)
2892. #combine.lon = c(glyMaa$coords.x2, glyDan$coords.x2)
2893. #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
2895. #I need to turn the points into a spatvector (
2896. library(terra)
2897. ptspres <- vect(cbind(CnatPres$long, CnatPres$lat))
2898. ptsholo <- vect(cbind(CnatHolo$long, CnatHolo$lat))
2901. # get random sample - background point locations -
2902. rasterobject <- rast(Prescomp.rasters)
2903. bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
2904. return.type = "points", n = 500)
2906. rasterobject <- rast(Holo.rasters)
2907. bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
2908. return.type = "points", n = 500)
2911. # get environmental data from occ points
2912. extractpres = na.omit(cbind(CnatPres[2:3], raster::extract(Prescomp.rasters, CnatPres[2:3]), rep(1, nrow(CnatPres))))
2914. extractholo = na.omit(cbind(CnatHolo[2:3], raster::extract(Holo.rasters, CnatHolo[2:3]), rep(1, nrow(CnatHolo))))
2917. #change name of last column to occ (represents presence = 1)
2918. colnames(extractpres)[ncol(extractpres)] = 'occ'
2919. colnames(extractholo)[ncol(extractholo)] = 'occ'
2922. # get environmental data from bg points
2923. #cant do stack-- need to do all???
2924. extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres)))))
2926. extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo)))))
2928. # change name of last column to occ (represents absence = 0)
2929. colnames(extbgpres)[ncol(extbgpres)] = 'occ'
2930. colnames(extbgholo)[ncol(extbgholo)] = 'occ'
2933. colnames(extbgpres)[1] = 'long'
2934. colnames(extbgpres)[2] = 'lat'
2936. colnames(extbgholo)[1] = 'long'
2937. colnames(extbgholo)[2] = 'lat'
2940. #colnames(extbgholo)[ncol(extbgholo)] = 'occ'
2942. # merge data from occ and bg
2943. datPresCnat = rbind(extractpres, extbgpres)
2945. datHoloCnat = rbind(extractholo, extbgholo)
2949. #principle component analysis
2951. pca.env <- dudi.pca(
2952. rbind(datHoloCnat, datPresCnat)[,3:6],
2953. scannf=FALSE,
2954. nf=2
2955. )
2957. # look at variable contribution
2958. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
2959. #dev.off()
2960. # Save pdf of PCA plot
2961. pdf('ecospat_Cnatpresholo.pdf')
2962. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
2964. scores.globclim<-pca.env$li # PCA scores for the whole study area
2965. scores.globclim<-pca.env$li # PCA scores for the whole study area (all points)
2967. scores.sppres <- suprow(pca.env,
2968. extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution
2970. scores.spholo <- suprow(pca.env,
2971. extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution
2973. scores.climholo <- suprow(pca.env,datHoloCnat[,3:6])$li # PCA scores for the first interval study area
2974. scores.climpres <- suprow(pca.env,datPresCnat[,3:6])$li # PCA scores for the second interval study area
2976. #dev.on()
2977. # density distribtuion for first interval
2978. grid.Cnatholo <- ecospat.grid.clim.dyn(
2979. glob = scores.globclim,
2980. glob1 = scores.climholo,
2981. sp = scores.spholo,
2982. R = 100,
2983. th.sp = 0
2984. )
2987. # density distribution for second interval
2988. grid.Cnatpres <- ecospat.grid.clim.dyn(
2989. glob = scores.globclim,
2990. glob1 = scores.climpres,
2991. sp = scores.sppres,
2992. R = 100,
2993. th.sp = 0
2994. )
2996. #dev.off()
2997. histpres <- hist(grid.Cnatpres$glob1)
2998. histholo <- hist(grid.Cnatholo$glob1)
2999. plot( histpres, col=rgb(0,0,1,1/4)) # first histogram
3000. plot( histholo, col=rgb(1,0,0,1/4), add=T)
3002. # Schoener's D metric and I metric
3003. D.overlap <- ecospat.niche.overlap (grid.Cnatholo, grid.Cnatpres, cor=T)
3004. D.overlap
3005. write.csv(D.overlap,"presCnat_holoCnat_I_D.csv",row.names=FALSE)
3007. ## Niche Equivalency Test
3008. ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D
3009. ## p.I: pvalue of the test on I
3010. ## Test for greater or lower equivalency
3011. eq.testgr <- ecospat.niche.equivalency.test(grid.Cnatholo, grid.Cnatpres,
3012. rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs
3013. eq.testlw <- ecospat.niche.equivalency.test(grid.Cnatholo, grid.Cnatpres,
3014. rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs
3017. nichdynindex <- ecospat.niche.dyn.index(grid.Cnatholo, grid.Cnatpres, intersection = NA)
3018. #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w
3019. nichdynindex
3022. # write p values of equivalency test
3023. 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)
3024. write.csv(p_EQ_DI,"glycim_EQ_TestCnatpresholo.csv",row.names=FALSE)
3026. p_EQ_DI
3029. ## Test for greater (niche conservatism) or lower (niche divergence) similarity
3030. sim.testgr <- ecospat.niche.similarity.test(grid.Cnatholo, grid.Cnatpres,
3031. rep=1000, overlap.alternative = "higher",
3032. rand.type=2,ncores=4)
3033. sim.testlw <- ecospat.niche.similarity.test(grid.Cnatholo, grid.Cnatpres,
3034. rep=1000, overlap.alternative = "lower",
3035. rand.type=2,ncores=4)
3037. # write p values of similarity test
3038. 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)
3039. write.csv(p_SIM_DI,"presholoCnat_SIM_Test.csv",row.names=FALSE)
3041. p_SIM_DI
3042. #p.D_GR p.I_GR p.D_LW p.I_LW
3043. #[1,] 0.04595405 0.03996004 0.9480519 0.9480519
3044. #dev.off()
3045. # Plot test distributions
3046. ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency")
3047. ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency")
3048. ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity")
3049. ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity")
3052. ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals
3053. ecospat.plot.niche (grid.Cnatpres, title='Cnatpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE)
3054. ecospat.plot.niche (grid.Cnatholo, title='Cnatholo', name.axis1='PC1', name.axis2='PC2')
3057. ## Plot niche overlap
3058. ecospat.plot.niche.dyn(z1 = grid.Cnatholo, z2 = grid.Cnatpres, quant=0.25, interest=2,
3059. title= "C. nat Holo to Pres", name.axis1="PC1",
3060. name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = 'black', colZ1 =
3061. "#313695", colZ2 = "#abd9d9", transparency = 40)
3062. ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white')
3063. 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)
3065. #___________________________
3067. ### Look at niche expantion, stability, and unfilling
3068. # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection
3069. dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Cnatholo, z2 = grid.Cnatpres, intersection=NA)
3070. dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Cnatholo, z2 = grid.Cnatpres, intersection=0)
3072. # write csv of niche dynmaic percentages
3073. Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w,
3074. "Overlapping"=dynam_overlap$dynamic.index.w)
3076. Niche_Ex_St
3077. write.csv(Niche_Ex_St,"GlycymDynamicsCnatholopres.csv",row.names=TRUE)
3079. ############################ PAST ECOSPAT ######################
3081. #### MODEL COMPARISON ####
3083. #need cleaned occurrence data and environmental data rasters
3085. ## Read in shapefiles of cleaned occ data ##
3086. ########################occurrence data######################
3087. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/Master_files_past_add_mask")
3088. PastPres = read.csv("masterpres.occs.csv")
3089. PastPres = as(PastPres,'data.frame')
3091. PastHolo = read.csv("masterholo.occs.csv")
3092. PastHolo = as(PastHolo,'data.frame')
3093. #convert to dataframe for use in ecospat
3095. ###################environmental###################
3098. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present")
3100. # Upload basic rasters for First interval (Present)
3101. Prescomp.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present", pattern='\\.tif$')
3102. Prescomp.files
3104. maxsalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxsalinity.tif")
3105. maxsalinitypres
3107. maxtemppres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/maxtemp.tif")
3108. maxtemppres
3110. rangesalinitypres <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangesalinity.tif")
3111. rangesalinitypres
3113. rangetemppres <- raster ("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/BIO-Oracle Future/MASKED/present/rangetemp.tif")
3114. rangetemppres
3117. # make sure order of how e-layers are read in is the same for each analysis extent (i.e., time bins or spatial areas)
3118. Prescomp.rasters <- raster::stack(c(maxsalinitypres, maxtemppres, rangesalinitypres, rangetemppres))
3121. wgs1984 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
3122. # optionally assigning a crs if one isn't provided
3123. crs(Prescomp.rasters) <- wgs1984
3124. names(Prescomp.rasters)
3125. names(Prescomp.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
3129. #Upload basic rasters for Second interval (Holocene)
3130. ## holocene+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3131. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam")
3132. Holo.files <- list.files(path="/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam", pattern='\\.tif$')
3133. #Holo.files
3136. maxsalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salmaxmask2.tif")
3137. maxsalinityholo2 <- maxsalinityholo/100
3139. maxtempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/tempmaxmask2.tif")
3140. maxtempholo2 <- maxtempholo/100
3142. rangesalinityholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/salrangemask2.tif")
3143. rangesalinityholo2 <- rangesalinityholo/100
3145. rangetempholo <- raster("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/marspec_env/holomaskcubadrjam/temprangemask2.tif")
3146. rangetempholo2 <- rangetempholo/100
3148. Holo.rasters <- raster::stack(maxsalinityholo2, maxtempholo2, rangesalinityholo2, rangetempholo2)
3150. # optionally assigning a crs if one isn't provided
3151. crs(Holo.rasters) <- wgs1984
3152. names(Holo.rasters)
3154. names(Holo.rasters) <- c('maxsalinity', 'maxtemp', 'rangesalinity', 'rangetemp')
3157. ### Ecospat Analysis of Niche Equivalency and Similary
3159. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3161. # Set working directory for saving plots and test results
3162. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/PAST/masked/maskedtwice/Ecospat_past_10_2_24")
3164. combine.lat = c(PastPres$lat, PastHolo$lat) #, PastLGM$lat)
3165. combine.lon = c(PastPres$long, PastHolo$long) #, PastLGM$long)
3168. #I need to turn the points into a spatvector (
3169. ptspres <- vect(cbind(PastPres$long, PastPres$lat))
3170. ptsholo <- vect(cbind(PastHolo$long, PastHolo$lat))
3172. # get random sample - background point locations -
3173. rasterobject <- rast(Prescomp.rasters)
3174. bgpres = ENMTools::background.buffer(points = ptspres , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
3175. return.type = "points", n = 500)
3177. rasterobject <- rast(Holo.rasters)
3178. bgholo = ENMTools::background.buffer(points = ptsholo , buffer.width = 200, buffer.type = "circles", mask = rasterobject,
3179. return.type = "points", n = 500)
3181. # get environmental data from occ points
3182. extractpres = na.omit(cbind(PastPres[2:3], raster::extract(Prescomp.rasters, PastPres[2:3]), rep(1, nrow(PastPres))))
3184. extractholo = na.omit(cbind(PastHolo[2:3], raster::extract(Holo.rasters, PastHolo[2:3]), rep(1, nrow(PastHolo))))
3188. #change name of last column to occ (represents presence = 1)
3189. colnames(extractpres)[ncol(extractpres)] = 'occ'
3190. colnames(extractholo)[ncol(extractholo)] = 'occ'
3193. # get environmental data from bg points
3194. #cant do stack-- need to do all???
3195. extbgpres = na.omit(cbind(geom(bgpres)[,3:4], raster::extract(Prescomp.rasters, geom(bgpres)[,3:4]), rep(0, nrow(geom(bgpres)))))
3197. extbgholo = na.omit(cbind(geom(bgholo)[,3:4], raster::extract(Holo.rasters, geom(bgholo)[,3:4]), rep(0, nrow(geom(bgholo)))))
3199. # change name of last column to occ (represents absence = 0)
3200. colnames(extbgpres)[ncol(extbgpres)] = 'occ'
3201. colnames(extbgholo)[ncol(extbgholo)] = 'occ'
3204. colnames(extbgpres)[1] = 'long'
3205. colnames(extbgpres)[2] = 'lat'
3207. colnames(extbgholo)[1] = 'long'
3208. colnames(extbgholo)[2] = 'lat'
3211. #colnames(extbgholo)[ncol(extbgholo)] = 'occ'
3213. # merge data from occ and bg
3214. datPresPast = rbind(extractpres, extbgpres)
3216. datHoloPast = rbind(extractholo, extbgholo)
3218. #principle component analysis
3220. pca.env <- dudi.pca(
3221. rbind(datHoloPast, datPresPast)[,3:6],
3222. scannf=FALSE,
3223. nf=2
3224. )
3226. # look at variable contribution
3227. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
3228. #dev.off()
3229. # Save pdf of PCA plot
3230. pdf('ecospat_Pastpresholo.pdf')
3231. ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
3233. scores.globclim<-pca.env$li # PCA scores for the whole study area
3234. scores.globclim<-pca.env$li # PCA scores for the whole study area (all points)
3236. scores.sppres <- suprow(pca.env,
3237. extractpres[which(extractpres[, 7]==1),3:6])$li # PCA scores for the first interval occ pts distribution
3239. scores.spholo <- suprow(pca.env,
3240. extractholo[which(extractholo[,7]==1),3:6])$li # PCA scores for the second interal occ pts distribution
3242. scores.climholo <- suprow(pca.env,datHoloPast[,3:6])$li # PCA scores for the first interval study area
3243. scores.climpres <- suprow(pca.env,datPresPast[,3:6])$li # PCA scores for the second interval study area
3245. #dev.on()
3246. # density distribtuion for first interval
3247. grid.Pastholo <- ecospat.grid.clim.dyn(
3248. glob = scores.globclim,
3249. glob1 = scores.climholo,
3250. sp = scores.spholo,
3251. R = 100,
3252. th.sp = 0
3253. )
3256. # density distribution for second interval
3257. grid.Pastpres <- ecospat.grid.clim.dyn(
3258. glob = scores.globclim,
3259. glob1 = scores.climpres,
3260. sp = scores.sppres,
3261. R = 100,
3262. th.sp = 0
3263. )
3265. histpres <- hist(grid.Pastpres$glob1)
3266. histholo <- hist(grid.Pastholo$glob1)
3267. plot( histpres, col=rgb(0,0,1,1/4)) # first histogram
3268. plot( histholo, col=rgb(1,0,0,1/4), add=T)
3270. # Schoener's D metric and I metric
3271. D.overlap <- ecospat.niche.overlap (grid.Pastholo, grid.Pastpres, cor=T)
3272. D.overlap
3273. write.csv(D.overlap,"presPast_holoPast_I_D.csv",row.names=FALSE)
3275. ## Niche Equivalency Test
3276. ## Obs: observed overlaps, sim: similulated overlaps, p.D: pvalue of the test on D
3277. ## p.I: pvalue of the test on I
3278. ## Test for greater or lower equivalency
3279. eq.testgr <- ecospat.niche.equivalency.test(grid.Pastholo, grid.Pastpres,
3280. rep=1000, overlap.alternative = "higher",ncores=4) ##rep = 1000 recommended for operational runs
3281. eq.testlw <- ecospat.niche.equivalency.test(grid.Pastholo, grid.Pastpres,
3282. rep=1000, overlap.alternative = "lower",ncores=4) ##rep = 1000 recommended for operational runs
3283. #this is if they are the same niche I think
3285. nichdynindex <- ecospat.niche.dyn.index(grid.Pastholo, grid.Pastpres, intersection = NA)
3286. #ecospat.niche.dyn.index(nativeGrid, invasiveGrid, intersection = 0.1)$dynamic.index.w
3287. nichdynindex
3290. # write p values of equivalency test
3291. 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)
3292. write.csv(p_EQ_DI,"glycim_EQ_TestPastpresholo.csv",row.names=FALSE)
3294. p_EQ_DI
3297. ## Test for greater (niche conservatism) or lower (niche divergence) similarity
3298. sim.testgr <- ecospat.niche.similarity.test(grid.Pastholo, grid.Pastpres,
3299. rep=1000, overlap.alternative = "higher",
3300. rand.type=2,ncores=4)
3301. sim.testlw <- ecospat.niche.similarity.test(grid.Pastholo, grid.Pastpres,
3302. rep=1000, overlap.alternative = "lower",
3303. rand.type=2,ncores=4)
3305. # write p values of similarity test
3306. 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)
3307. write.csv(p_SIM_DI,"presholoPast_SIM_Test.csv",row.names=FALSE)
3309. p_SIM_DI
3310. #p.D_GR p.I_GR p.D_LW p.I_LW
3312. # Plot test distributions
3313. ecospat.plot.overlap.test(eq.testgr, "D", "Greater Equivalency")
3314. ecospat.plot.overlap.test(eq.testlw, "D", "Lower Equivalency")
3315. ecospat.plot.overlap.test(sim.testgr, "D", "Greater Similarity")
3316. ecospat.plot.overlap.test(sim.testlw, "D", "Lower Similarity")
3320. ## Plot niche space for first (grid.climholo) and second (grid.climpres) intervals
3321. ecospat.plot.niche (grid.Pastpres, title='Pastpres', name.axis1='PC1', name.axis2='PC2', cor=FALSE)
3322. ecospat.plot.niche (grid.Pastholo, title='Pastholo', name.axis1='PC1', name.axis2='PC2')
3325. ## Plot niche overlap
3326. ecospat.plot.niche.dyn(z1 = grid.Pastholo, z2 = grid.Pastpres, quant=0.25, interest=2,
3327. title= "P. ast Holo to Pres", name.axis1="PC1",
3328. name.axis2="PC2", col.unf = "#313695", col.exp = "#abd9d9", col.stab = '#74add1', colZ1 =
3329. "#313695", colZ2 = "#abd9d9", transparency = 40)
3330. ecospat.shift.centroids(sp1 = scores.spholo, sp2 = scores.sppres, scores.climholo, scores.climpres, col = 'white')
3331. 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)
3333. #___________________________
3336. ### Look at niche expantion, stability, and unfilling
3337. # NA=analysis on entire area, 0=analysis on only overlapping, 0.05=analysis on 5th quantile intersection
3338. dynam_allAreas = ecospat.niche.dyn.index (z1 = grid.Pastholo, z2 = grid.Pastpres, intersection=NA)
3339. dynam_overlap = ecospat.niche.dyn.index (z1 = grid.Pastholo, z2 = grid.Pastpres, intersection=0)
3341. # write csv of niche dynmaic percentages
3342. Niche_Ex_St = cbind("AllAreas"=dynam_allAreas$dynamic.index.w,
3343. "Overlapping"=dynam_overlap$dynamic.index.w)
3345. Niche_Ex_St
3346. write.csv(Niche_Ex_St,"GlycymDynamicsPastholopres.csv",row.names=TRUE)
3388. #################&&&&&&&&&&&&&&&&&&&&&############ ENM #################
3390. ########################################## ACER PRES ENM ################################
3391. ##### Loading in the Data --------------------------------------------------------------------------------------------------
3392. # master background file;
3393. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/ACERV/masked/Master_files_acerv_add_mask")
3395. All.back <- read.csv('masterall.background.csv', header=T)
3397. #filter to time period background
3398. pres.back <- All.back %>% filter(time.bin== c('Present'))
3400. Holo.back <- All.back %>% filter(time.bin=='Holocene')
3401. RCP452050.back <- All.back %>% filter(time.bin=='RCP452050')
3402. RCP852050.back <- All.back %>% filter(time.bin=='RCP852050')
3403. RCP452100.back <- All.back %>% filter(time.bin=='RCP452100')
3404. RCP852100.back <- All.back %>% filter(time.bin=='RCP852100')
3406. #save background I am using
3407. 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)
3409. # master occurrence file;
3410. all.occ <- read.csv('masterall.occs.csv', header=T)
3412. #reading in specific time period files
3413. pres.occs <- read.csv('masterpres.occs.csv', header=T)
3415. #check if background and occurance have same column names
3416. identical(colnames(pres.back), colnames(pres.occs))
3418. ##### Setting up the models ------------------------------------------------------------------------------------------------
3420. # create null df - empty dataframe that null data will go into later
3421. pres.null <- create.null.df('Acervicornis', 'Present', 'Caribbean') #function(taxon.name, time.bin, extent){
3423. # create summary df
3424. pres.summary <- create.summary.df('Acervicornis', 'Present', 'Caribbean')
3426. create.folders.for.maxent(pres.summary)
3428. # run null model
3429. Acerv.null <- null.aic(null.df = pres.null,
3430. occs = pres.occs,
3431. background = pres.back,
3432. first.occ.col = 10)
3434. ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------
3436. # find most optimized model
3437. Acerv.optim <- optimize.maxent.likelihood(pres.summary, # name of the output summary file
3438. occs = pres.occs, # species occurrences
3439. background = pres.back, # background for the pres
3440. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3441. first.occ.col = 10)
3443. View(Acerv.optim)
3444. #want model to be at least more than 2 AICc lower than the null
3445. #also plan to check pROC and response curves to pick best model
3447. 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)
3448. ### CHECK RESPONSE CURVES FOR REALISM!!!!!
3450. ##### Running Cross-Validation ---------------------------------------------------------------------------------------------
3452. ### LQP 0.10
3453. # create eval object
3454. # default f.class = LQP, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP
3455. Acervpres.eval_LQP0.10 <- create.eval.df('Acerv', 'Present', 'Caribbean', beta.values = 0.10, f.class = 'LQP')
3456. Acervpres.eval_LQP0.10
3458. # create folders for Baculites in the Cenomanian
3459. create.folders.for.maxent(Acervpres.eval_LQP0.10)
3461. # run eval object
3462. Acervpres.eval_LQP0.10 <- maxent.crossval.error(eval.df = Acervpres.eval_LQP0.10,
3463. occs = pres.occs, # species occurrences
3464. background = pres.back, # background for the pres
3465. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3466. first.occ.col = 10,
3467. first.test.col = 16,
3468. omission.rate = 0.025, #express as proportion
3469. all.background = All.back)
3471. #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.)
3472. # response curves good
3474. ### LQP 0.25
3475. Acervpres.eval_LQP0.25 <- create.eval.df('Acerv', 'Present', 'Caribbean', beta.values = 0.25, f.class = 'LQP')
3476. Acervpres.eval_LQP0.25
3478. # create folders
3479. create.folders.for.maxent(Acervpres.eval_LQP0.25)
3481. # run eval object
3482. Acervpres.eval_LQP0.25 <- maxent.crossval.error(eval.df = Acervpres.eval_LQP0.25,
3483. occs = pres.occs, # species occurrences
3484. background = pres.back, # background for the pres
3485. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3486. first.occ.col = 10,
3487. first.test.col = 16,
3488. omission.rate = 0, #express as proportion
3489. all.background = All.back)
3491. #response curves similar- not as good, donmt go down as far
3493. ### analyzing summary files
3494. ## if running multiple eval models, check for:
3495. # 1.) does one model setting have systematically higher omission rates?
3496. # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type?
3497. # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps
3499. Acervpres.eval_LQP0.10$summary # pROC 1.672626
3500. Acervpres.eval_LQP0.25$summary # pROC 1.662092
3503. # going with LQP 0.10
3506. ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------
3508. ### LQP 0.10
3509. # evaluate model to calculate weighted suitibility and stdev
3510. Acervpres.means_LQP0.10 <- maxent.eval(eval = Acervpres.eval_LQP0.10)
3511. thresh = min(Acervpres.means_LQP0.10$occ$w.mean)
3512. print(thresh)
3513. # this is the threshold
3514. # save this value!!
3517. ### plot model suitability/uncertainty plots
3518. # the only input you need is the object made from the maxent.eval function
3519. suit.uncert.plot(Acervpres.eval_LQP0.10)
3520. ggsave('Acervpres.eval_LQP0.10.pdf')
3525. ##### Projecting Model to all extents --------------------------------------------------------------------------------------
3526. Acervpres.thresh_LQP0.10 <- 0.001166019
3528. Acervpres.everything_LQP0.10 <- maxent.everything(eval = Acervpres.eval_LQP0.10,
3529. thresh = Acervpres.thresh_LQP0.10,
3530. means = Acervpres.means_LQP0.10,
3531. everything= All.back,
3532. predic = 6:9)
3533. View(Acervpres.everything_LQP0.10)
3535. ##### Running Informed MESS Analysis ---------------------------------------------------------------------------------------
3537. # defining the tolerance vector
3538. # 1 = can extrapolate to non-analog conditions
3539. # 0 = cannot extrapolate to non-analog conditions
3540. Acerv.tolerance <- c(1,1, # max salinity
3541. 1,1, # max temp
3542. 1,1, # range salinity
3543. 1,1) # range temp
3545. ################### running mess holocene
3546. Acerv.mess.Holo <- informed.mess(ref.extent = pres.back,
3547. mess.extent = Holo.back, #maybe they are crazy high because project to same, lets try diff
3548. coord.cols = 2:3,
3549. predic = 6:9,
3550. tolerance = Acerv.tolerance)
3552. ################### running mess 45 2050
3553. Acerv.mess.452050 <- informed.mess(ref.extent = pres.back,
3554. mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff
3555. coord.cols = 2:3,
3556. predic = 6:9,
3557. tolerance = Acerv.tolerance)
3559. ################### running mess 45 2100
3560. Acerv.mess.452100 <- informed.mess(ref.extent = pres.back,
3561. mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff
3562. coord.cols = 2:3,
3563. predic = 6:9,
3564. tolerance = Acerv.tolerance)
3566. ################### running mess 85 2050
3567. Acerv.mess.852050 <- informed.mess(ref.extent = pres.back,
3568. mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff
3569. coord.cols = 2:3,
3570. predic = 6:9,
3571. tolerance = Acerv.tolerance)
3573. ################### running mess 85 2100
3574. Acerv.mess.852100 <- informed.mess(ref.extent = pres.back,
3575. mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff
3576. coord.cols = 2:3,
3577. predic = 6:9,
3578. tolerance = Acerv.tolerance)
3580. #### Acerv mess all
3581. Acerv.mess.all <- informed.mess(ref.extent = pres.back,
3582. mess.extent = All.back,
3583. coord.cols = 2:3,
3584. predic = 6:9,
3585. tolerance = Acerv.tolerance)
3588. ##### Saving Everything ----------------------------------------------------------------------------------------------------
3589. pres.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin')
3590. pres.output <- cbind(pres.output,Acervpres.everything_LQP0.10)
3591. pres.output$Name <- 'Acervicornis'
3592. pres.with.mess <- cbind(pres.output, Acerv.mess.all)
3593. 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)
3595. # null model
3596. write.csv(Acerv.null, 'Acerv.null.csv', row.names = F)
3597. # model optimized summary object
3598. write.csv(Acerv.optim, 'Acerv.optim.summary.csv', row.names = F)
3599. # model projected to all extents
3600. write.csv(Acervpres.everything_LQP0.10, 'Acerv.LQP0.10.csv', row.names = F)
3601. # model mess analysis
3602. write.csv(Acerv.mess.Holo, 'Acerv.mess.Holo.csv', row.names = F)
3603. write.csv(Acerv.mess.452050, 'Acerv.mess.452050.csv', row.names = F)
3604. write.csv(Acerv.mess.452100, 'Acerv.mess.452100.csv', row.names = F)
3605. write.csv(Acerv.mess.852050, 'Acerv.mess.852050.csv', row.names = F)
3606. write.csv(Acerv.mess.852100, 'Acerv.mess.852100.csv', row.names = F)
3613. ########################################## ACER HOLO ENM ################################
3614. # master background file;
3615. All.back <- read.csv('masterall.background.csv', header=T)
3617. pres.back <- All.back %>% filter(time.bin=='Present')
3618. Holo.back <- All.back %>% filter(time.bin=='Holocene')
3619. RCP452050.back <- All.back %>% filter(time.bin=='RCP452050')
3620. RCP852050.back <- All.back %>% filter(time.bin=='RCP852050')
3621. RCP452100.back <- All.back %>% filter(time.bin=='RCP452100')
3622. RCP852100.back <- All.back %>% filter(time.bin=='RCP852100')
3624. #save background I am using
3625. 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)
3628. # master occurrence file;
3629. all.occ <- read.csv('masterall.occs.csv', header=T)
3631. #reading in specific time period files
3632. Pres.occs <- read.csv('masterpres.occs.csv', header=T)
3633. Holo.occs <- read.csv('masterholo.occs.csv', header=T)
3635. #check if background and occurance have same column names
3636. identical(colnames(Holo.back), colnames(Holo.occs))
3639. ##### Setting up the models ------------------------------------------------------------------------------------------------
3641. # create null df - empty dataframe that null data will go into later
3642. Holo.null <- create.null.df('Acervicornis', 'Holocene', 'Caribbean') #function(taxon.name, time.bin, extent)
3644. # create summary df
3645. Holo.summary <- create.summary.df('Acervicornis', 'Holocene', 'Caribbean')
3647. create.folders.for.maxent(Holo.summary)
3649. # run null model
3650. Acerv.null <- null.aic(null.df = Holo.null,
3651. occs = Holo.occs,
3652. background = Holo.back,
3653. first.occ.col = 10)
3657. ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------
3659. # find most optimized model
3660. Acerv.optim <- optimize.maxent.likelihood(Holo.summary, # name of the output summary file
3661. occs = Holo.occs, # species occurrences
3662. background = Holo.back, # background for the Holo
3663. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3664. first.occ.col = 10)
3666. View(Acerv.optim)
3667. #want model to be at least more than 2 AICc lower than the null
3668. #also plan to check pROC and response curves to pick best model
3670. 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)
3671. ### CHECK RESPONSE CURVES FOR REALISM!!!!!
3673. ##### Running Cross-Validation ---------------------------------------------------------------------------------------------
3675. ### Q 0.05
3676. # create eval object
3677. # default f.class = Q, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP
3678. AcervHolo.eval_Q0.05 <- create.eval.df('Acerv', 'Holocene', 'Caribbean', beta.values = 0.05, f.class = 'Q')
3679. AcervHolo.eval_Q0.05
3681. # create folders for Baculites in the Cenomanian
3682. create.folders.for.maxent(AcervHolo.eval_Q0.05)
3684. # run eval object
3685. AcervHolo.eval_Q0.05 <- maxent.crossval.error(eval.df = AcervHolo.eval_Q0.05,
3686. occs = Holo.occs, # species occurrences
3687. background = Holo.back, # background for the Holo
3688. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3689. first.occ.col = 10,
3690. first.test.col = 16,
3691. omission.rate = 0, #express as proportion
3692. all.background = All.back)
3695. #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.)
3696. # WORKED!!!! curves are curves!!!!!!! or at least temp is curves
3698. ### Q0.025
3699. AcervHolo.eval_Q0.025 <- create.eval.df('Acerv', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'Q')
3700. AcervHolo.eval_Q0.025
3702. # create folders
3703. create.folders.for.maxent(AcervHolo.eval_Q0.025)
3705. # run eval object
3706. AcervHolo.eval_Q0.025 <- maxent.crossval.error(eval.df = AcervHolo.eval_Q0.025,
3707. occs = Holo.occs, # species occurrences
3708. background = Holo.back, # background for the Holo
3709. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3710. first.occ.col = 10,
3711. first.test.col = 16,
3712. omission.rate = 0, #express as proportion
3713. all.background = All.back)
3715. #curves are not as good going with other
3719. #still a good model compared to null and has good curves
3720. ### Q 0.10
3721. AcervHolo.eval_LQP0.025 <- create.eval.df('Acerv', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'LQP')
3722. AcervHolo.eval_LQP0.025
3724. # create folders
3725. create.folders.for.maxent(AcervHolo.eval_LQP0.025)
3727. # run eval object
3728. AcervHolo.eval_LQP0.025 <- maxent.crossval.error(eval.df = AcervHolo.eval_LQP0.025,
3729. occs = Holo.occs, # species occurrences
3730. background = Holo.back, # background for the Holo
3731. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3732. first.occ.col = 10,
3733. first.test.col = 16,
3734. omission.rate = 0, #express as proportion
3735. all.background = All.back)
3737. ### analyzing summary files
3738. ## if running multiple eval models, check for:
3739. # 1.) does one model setting have systematically higher omission rates?
3740. # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type?
3741. # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps
3743. AcervHolo.eval_Q0.025$summary # 1.823454
3744. AcervHolo.eval_Q0.05$summary # 1.821950
3745. AcervHolo.eval_LQP0.025$summary #1.887521
3747. # neither mode seems to be systematically better basically identical, i will do LQP0.025
3750. ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------
3752. ### Q 0.10
3753. # evaluate model to calculate weighted suitibility and stdev
3754. AcervHolo.means_LQP0.025 <- maxent.eval(eval = AcervHolo.eval_LQP0.025)
3755. thresh = min(AcervHolo.means_LQP0.025$occ$w.mean)
3756. print(thresh)
3757. # 0.08827532
3758. # this is the threshold
3759. # save this value!!
3762. ### Q 0.5
3766. ### plot model suitability/uncertainty plots
3767. # the only input you need is the object made from the maxent.eval function
3768. suit.uncert.plot(AcervHolo.eval_LQP0.025)
3769. ggsave('AcervHolo.eval_LQP0.025.pdf')
3771. ### deciding to go with Q 0.10 since it has slightly more realistic response curves
3774. ##### Projecting Model to all extents --------------------------------------------------------------------------------------
3775. AcervHolo.thresh_LQP0.025 <- 0.08827532
3777. AcervHolo.everything_LQP0.025 <- maxent.everything(eval = AcervHolo.eval_LQP0.025,
3778. thresh = AcervHolo.thresh_LQP0.025,
3779. means = AcervHolo.means_LQP0.025,
3780. everything= All.back,
3781. predic = 6:9)
3782. View(AcervHolo.everything_LQP0.025)
3784. ##### Running Informed MESS Analysis ---------------------------------------------------------------------------------------
3786. # defining the tolerance vector
3787. # 1 = can extrapolate to non-analog conditions
3788. # 0 = cannot extrapolate to non-analog conditions
3789. Acerv.tolerance <- c(1,1, # max salinity
3790. 1,1, # max temp
3791. 1,1, # range salinity
3792. 1,0) # range temp
3796. ################### running mess present
3797. Acerv.mess.pres <- informed.mess(ref.extent = Holo.back,
3798. mess.extent = pres.back, #maybe they are crazy high because project to same, lets try diff
3799. coord.cols = 2:3,
3800. predic = 6:9,
3801. tolerance = Acerv.tolerance)
3804. ################### running mess 45 2050
3805. Acerv.mess.452050 <- informed.mess(ref.extent = Holo.back,
3806. mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff
3807. coord.cols = 2:3,
3808. predic = 6:9,
3809. tolerance = Acerv.tolerance)
3812. ################### running mess 45 2100
3813. Acerv.mess.452100 <- informed.mess(ref.extent = Holo.back,
3814. mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff
3815. coord.cols = 2:3,
3816. predic = 6:9,
3817. tolerance = Acerv.tolerance)
3820. ################### running mess 85 2050
3821. Acerv.mess.852050 <- informed.mess(ref.extent = Holo.back,
3822. mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff
3823. coord.cols = 2:3,
3824. predic = 6:9,
3825. tolerance = Acerv.tolerance)
3828. ################### running mess 85 2100
3829. Acerv.mess.852100 <- informed.mess(ref.extent = Holo.back,
3830. mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff
3831. coord.cols = 2:3,
3832. predic = 6:9,
3833. tolerance = Acerv.tolerance)
3835. #mess all
3836. Acerv.mess.all <- informed.mess(ref.extent = Holo.back,
3837. mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff
3838. coord.cols = 2:3,
3839. predic = 6:9,
3840. tolerance = Acerv.tolerance)
3844. ##### Saving Everything ----------------------------------------------------------------------------------------------------
3846. pres.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin')
3847. pres.output <- cbind(pres.output, AcervHolo.everything_LQP0.025)
3848. pres.output$Name <- 'Acervicornis'
3849. pres.with.mess <- cbind(pres.output, Acerv.mess.all)
3850. 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)
3854. # null model
3855. write.csv(Acerv.null, 'Acerv.null.csv', row.names = F)
3856. # model optimized summary object
3857. write.csv(Acerv.optim, 'Acerv.optim.summary.csv', row.names = F)
3858. # model projected to all extents
3859. write.csv(AcervHolo.everything_LQP0.025, 'Acerv.LQP0.025.csv', row.names = F)
3860. # model mess analysis
3861. write.csv(Acerv.mess.pres, 'Acerv.mess.pres.csv', row.names = F)
3862. write.csv(Acerv.mess.452050, 'Acerv.mess.452050.csv', row.names = F)
3863. write.csv(Acerv.mess.452100, 'Acerv.mess.452100.csv', row.names = F)
3864. write.csv(Acerv.mess.852050, 'Acerv.mess.852050.csv', row.names = F)
3865. write.csv(Acerv.mess.852100, 'Acerv.mess.852100.csv', row.names = F)
3868. ######################################### APAL PRES ENM ############
3870. ##### Loading in the Data --------------------------------------------------------------------------------------------------
3872. # master background file;
3873. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/Master_files_apal_add_mask")
3875. All.back <- read.csv('masterall.background.csv', header=T)
3877. #filter to time period background
3878. pres.back <- All.back %>% filter(time.bin== c('Present'))
3880. Holo.back <- All.back %>% filter(time.bin=='Holocene')
3881. RCP452050.back <- All.back %>% filter(time.bin=='RCP452050')
3882. RCP852050.back <- All.back %>% filter(time.bin=='RCP852050')
3883. RCP452100.back <- All.back %>% filter(time.bin=='RCP452100')
3884. RCP852100.back <- All.back %>% filter(time.bin=='RCP852100')
3886. #save background I am using
3887. 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)
3889. # master occurrence file;
3890. all.occ <- read.csv('masterall.occs.csv', header=T)
3892. #reading in specific time period files
3893. pres.occs <- read.csv('masterpres.occs.csv', header=T)
3895. #check if background and occurance have same column names
3896. identical(colnames(pres.back), colnames(pres.occs))
3899. ##### Setting up the models ------------------------------------------------------------------------------------------------
3901. # working directory for models
3902. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/models_Apal_pres_add_mask")
3904. # create null df - empty dataframe that null data will go into later
3905. pres.null <- create.null.df('Apalmata', 'Present', 'Caribbean') #function(taxon.name, time.bin, extent){
3907. # create summary df
3908. pres.summary <- create.summary.df('Apalmata', 'Present', 'Caribbean')
3910. create.folders.for.maxent(pres.summary)
3912. # run null model
3913. Apal.null <- null.aic(null.df = pres.null,
3914. occs = pres.occs,
3915. background = pres.back,
3916. first.occ.col = 10)
3920. ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------
3922. # find most optimized model
3923. Apal.optim <- optimize.maxent.likelihood(pres.summary, # name of the output summary file
3924. occs = pres.occs, # species occurrences
3925. background = pres.back, # background for the pres
3926. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3927. first.occ.col = 10)
3929. View(Apal.optim)
3930. #want model to be at least more than 2 AICc lower than the null
3931. #also plan to check pROC and response curves to pick best model
3933. 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)
3934. ### CHECK RESPONSE CURVES FOR REALISM!!!!!
3936. ##### Running Cross-Validation ---------------------------------------------------------------------------------------------
3938. ### LQP .05
3939. # create eval object
3940. # default f.class = LQP, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP
3941. Apalpres.eval_LQP0.05 <- create.eval.df('Apal', 'Present', 'Caribbean', beta.values = 0.05, f.class = 'LQP')
3942. Apalpres.eval_LQP0.05
3944. # create folders for Baculites in the Cenomanian
3945. create.folders.for.maxent(Apalpres.eval_LQP0.05)
3947. # run eval object
3948. Apalpres.eval_LQP0.05 <- maxent.crossval.error(eval.df = Apalpres.eval_LQP0.05,
3949. occs = pres.occs, # species occurrences
3950. background = pres.back, # background for the pres
3951. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3952. first.occ.col = 10,
3953. first.test.col = 16,
3954. omission.rate = 0, #express as proportion
3955. all.background = All.back)
3958. #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.)
3959. # 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?
3961. ### LQP 0.025
3962. Apalpres.eval_LQP0.025 <- create.eval.df('Apal', 'Present', 'Caribbean', beta.values = 0.025, f.class = 'LQP')
3963. Apalpres.eval_LQP0.025
3965. # create folders
3966. create.folders.for.maxent(Apalpres.eval_LQP0.025)
3968. # run eval object
3969. Apalpres.eval_LQP0.025 <- maxent.crossval.error(eval.df = Apalpres.eval_LQP0.025,
3970. occs = pres.occs, # species occurrences
3971. background = pres.back, # background for the pres
3972. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
3973. first.occ.col = 10,
3974. first.test.col = 16,
3975. omission.rate = 0, #express as proportion
3976. all.background = All.back)
3978. #response curves similar- only thing is max temp starts going back down at high temps but not all the way.
3980. ### analyzing summary files
3981. ## if running multiple eval models, check for:
3982. # 1.) does one model setting have systematically higher omission rates?
3983. # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type?
3984. # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps
3986. Apalpres.eval_LQP0.05$summary # pROC 1.690565
3987. Apalpres.eval_LQP0.025$summary # pROC 1.685479
3990. # going with LQP 0.025 since slightly higher pROC?? im not sure which is better?
3993. ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------
3995. ### LQP 0.025
3996. # evaluate model to calculate weighted suitibility and stdev
3997. Apalpres.means_LQP0.025 <- maxent.eval(eval = Apalpres.eval_LQP0.025)
3998. thresh = min(Apalpres.means_LQP0.025$occ$w.mean)
3999. print(thresh)
4000. # 0.2613874
4001. # this is the threshold
4002. # save this value!!
4005. ### Q 0.5
4009. ### plot model suitability/uncertainty plots
4010. # the only input you need is the object made from the maxent.eval function
4011. suit.uncert.plot(Apalpres.eval_LQP0.025)
4012. ggsave('Apalpres.eval_LQP0.025.pdf')
4014. ### deciding to go with Q 2.0 since it has slightly more realistic response curves
4017. ##### Projecting Model to all extents --------------------------------------------------------------------------------------
4018. Apalpres.thresh_LQP0.025 <- 0.2613874
4020. Apalpres.everything_LQP0.025 <- maxent.everything(eval = Apalpres.eval_LQP0.025,
4021. thresh = Apalpres.thresh_LQP0.025,
4022. means = Apalpres.means_LQP0.025,
4023. everything= All.back,
4024. predic = 6:9)
4025. View(Apalpres.everything_LQP0.025)
4027. ##### Running Informed MESS Analysis ---------------------------------------------------------------------------------------
4029. # defining the tolerance vector
4030. # 1 = can extrapolate to non-analog conditions
4031. # 0 = cannot extrapolate to non-analog conditions
4032. Apal.tolerance <- c(1,1, # max salinity
4033. 1,1, # max temp
4034. 1,1, # range salinity
4035. 1,1) # range temp
4037. ################## running mess holocene
4038. Apal.mess.Holo <- informed.mess(ref.extent = pres.back,
4039. mess.extent = Holo.back, #maybe they are crazy high because project to same, lets try diff
4040. coord.cols = 2:3,
4041. predic = 6:9,
4042. tolerance = Apal.tolerance)
4045. ################### running mess 45 2050
4046. Apal.mess.452050 <- informed.mess(ref.extent = pres.back,
4047. mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff
4048. coord.cols = 2:3,
4049. predic = 6:9,
4050. tolerance = Apal.tolerance)
4052. ################### running mess 45 2100
4053. Apal.mess.452100 <- informed.mess(ref.extent = pres.back,
4054. mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff
4055. coord.cols = 2:3,
4056. predic = 6:9,
4057. tolerance = Apal.tolerance)
4059. ################### running mess 85 2050
4060. Apal.mess.852050 <- informed.mess(ref.extent = pres.back,
4061. mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff
4062. coord.cols = 2:3,
4063. predic = 6:9,
4064. tolerance = Apal.tolerance)
4066. ################### running mess 85 2100
4067. Apal.mess.852100 <- informed.mess(ref.extent = pres.back,
4068. mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff
4069. coord.cols = 2:3,
4070. predic = 6:9,
4071. tolerance = Apal.tolerance)
4073. #### apal mess all
4074. Apal.mess.all <- informed.mess(ref.extent = pres.back,
4075. mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff
4076. coord.cols = 2:3,
4077. predic = 6:9,
4078. tolerance = Apal.tolerance)
4081. ##### Saving Everything ----------------------------------------------------------------------------------------------------
4083. apal.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin')
4084. apal.output <- cbind(apal.output, Apalpres.everything_LQP0.025)
4085. apal.output$Name <- 'Apalmata'
4086. apal.with.mess <- cbind(apal.output, Apal.mess.all)
4087. 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)
4089. # null model
4090. write.csv(Apal.null, 'Apal.null.csv', row.names = F)
4091. # model optimized summary object
4092. write.csv(Apal.optim, 'Apal.optim.summary.csv', row.names = F)
4093. # model projected to all extents
4094. write.csv(Apalpres.everything_LQP0.025, 'Apal.LQP0.025.csv', row.names = F)
4095. # model mess analysis
4096. write.csv(Apal.mess.Holo, 'Apal.mess.Holo.csv', row.names = F)
4097. write.csv(Apal.mess.452050, 'Apal.mess.452050.csv', row.names = F)
4098. write.csv(Apal.mess.452100, 'Apal.mess.452100.csv', row.names = F)
4099. write.csv(Apal.mess.852050, 'Apal.mess.852050.csv', row.names = F)
4100. write.csv(Apal.mess.852100, 'Apal.mess.852100.csv', row.names = F)
4105. ####################################### APAL HOLO ENM #############
4107. ##### Loading in the Data --------------------------------------------------------------------------------------------------
4109. # master background file;
4110. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/APAL/masked/maskedtwice/Master_files_apal_add_mask")
4112. All.back <- read.csv('masterall.background.csv', header=T)
4114. #filter to time period background
4115. LGM.back <- All.back %>% filter(time.bin== c('LGM'))
4117. pres.back <- All.back %>% filter(time.bin=='Present')
4118. Holo.back <- All.back %>% filter(time.bin=='Holocene')
4119. Holo2.back <- All.back %>% filter(time.bin=='Holocene2')
4120. RCP452050.back <- All.back %>% filter(time.bin=='RCP452050')
4121. RCP852050.back <- All.back %>% filter(time.bin=='RCP852050')
4122. RCP452100.back <- All.back %>% filter(time.bin=='RCP452100')
4123. RCP852100.back <- All.back %>% filter(time.bin=='RCP852100')
4125. #save background I am using
4126. 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)
4128. # master occurrence file;
4129. all.occ <- read.csv('masterall.occs.csv', header=T)
4131. #reading in specific time period files
4132. Pres.occs <- read.csv('masterpres.occs.csv', header=T)
4133. Holo.occs <- read.csv('masterholo.occs.csv', header=T)
4135. #check if background and occurance have same column names
4136. identical(colnames(Holo.back), colnames(Holo.occs))
4139. ##### Setting up the models ------------------------------------------------------------------------------------------------
4142. # create null df - empty dataframe that null data will go into later
4143. Holo.null <- create.null.df('Apalmata', 'Holocene', 'Caribbean') #function(taxon.name, time.bin, extent)
4145. # create summary df
4146. Holo.summary <- create.summary.df('Apalmata', 'Holocene', 'Caribbean')
4148. create.folders.for.maxent(Holo.summary)
4150. # run null model
4151. Apal.null <- null.aic(null.df = Holo.null,
4152. occs = Holo.occs,
4153. background = Holo.back,
4154. first.occ.col = 10)
4158. ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------
4160. # find most optimized model
4161. Apal.optim <- optimize.maxent.likelihood(Holo.summary, # name of the output summary file
4162. occs = Holo.occs, # species occurrences
4163. background = Holo.back, # background for the Holo
4164. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
4165. first.occ.col = 10)
4167. View(Apal.optim)
4168. #want model to be at least more than 2 AICc lower than the null
4169. #also plan to check pROC and response curves to pick best model
4171. 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)
4172. ### CHECK RESPONSE CURVES FOR REALISM!!!!!
4174. ##### Running Cross-Validation ---------------------------------------------------------------------------------------------
4176. ### LQP 0.025
4177. # create eval object
4178. # default f.class = Q, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP
4179. ApalHolo.eval_LQP0.025 <- create.eval.df('Apal', 'Holocene', 'Caribbean', beta.values = 0.025, f.class = 'LQP')
4180. ApalHolo.eval_LQP0.025
4182. # create folders for Baculites in the Cenomanian
4183. create.folders.for.maxent(ApalHolo.eval_LQP0.025)
4185. # run eval object
4186. ApalHolo.eval_LQP0.025 <- maxent.crossval.error(eval.df = ApalHolo.eval_LQP0.025,
4187. occs = Holo.occs, # species occurrences
4188. background = Holo.back, # background for the Holo
4189. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
4190. first.occ.col = 10,
4191. first.test.col = 16,
4192. omission.rate = 0, #express as proportion
4193. all.background = All.back)
4196. #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.)
4197. # WORKED!!!! curves are curves!!!!!!!
4199. ### LQP 0.10
4200. ApalHolo.eval_LQP0.10 <- create.eval.df('Apal', 'Holocene', 'Caribbean', beta.values = 0.10, f.class = 'LQP')
4201. ApalHolo.eval_LQP0.10
4203. # create folders
4204. create.folders.for.maxent(ApalHolo.eval_LQP0.10)
4206. # run eval object
4207. ApalHolo.eval_LQP0.10 <- maxent.crossval.error(eval.df = ApalHolo.eval_LQP0.10,
4208. occs = Holo.occs, # species occurrences
4209. background = Holo.back, # background for the Holo
4210. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
4211. first.occ.col = 10,
4212. first.test.col = 16,
4213. omission.rate = 0, #express as proportion
4214. all.background = All.back)
4216. #curves are not as good going with other
4218. ### analyzing summary files
4219. ## if running multiple eval models, check for:
4220. # 1.) does one model setting have systematically higher omission rates?
4221. # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type?
4222. # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps
4224. ApalHolo.eval_LQP0.10$summary # 1.958130
4225. ApalHolo.eval_LQP0.025$summary #1.971328
4228. # neither mode seems to be systematically better basically identical, i will do LQP0.025
4231. ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------
4233. ### Q 0.10
4234. # evaluate model to calculate weighted suitibility and stdev
4235. ApalHolo.means_LQP0.025 <- maxent.eval(eval = ApalHolo.eval_LQP0.025)
4236. thresh = min(ApalHolo.means_LQP0.025$occ$w.mean)
4237. print(thresh)
4238. # 0.01926973
4239. # this is the threshold
4240. # save this value!!
4244. ### plot model suitability/uncertainty plots
4245. # the only input you need is the object made from the maxent.eval function
4246. suit.uncert.plot(ApalHolo.eval_LQP0.025)
4247. ggsave('ApalHolo.eval_LQP0.025.pdf')
4249. ### deciding to go with Q 0.10 since it has slightly more realistic response curves
4252. ##### Projecting Model to all extents --------------------------------------------------------------------------------------
4253. ApalHolo.thresh_LQP0.025 <- 0.01926973
4255. ApalHolo.everything_LQP0.025 <- maxent.everything(eval = ApalHolo.eval_LQP0.025,
4256. thresh = ApalHolo.thresh_LQP0.025,
4257. means = ApalHolo.means_LQP0.025,
4258. everything= All.back,
4259. predic = 6:9)
4260. View(ApalHolo.everything_LQP0.025)
4262. ##### Running Informed MESS Analysis ---------------------------------------------------------------------------------------
4264. # defining the tolerance vector
4265. # 1 = can extrapolate to non-analog conditions
4266. # 0 = cannot extrapolate to non-analog conditions
4267. Apal.tolerance <- c(1,1, # max salinity
4268. 1,1, # max temp
4269. 0,1, # range salinity
4270. 1,1) # range temp
4272. ################### running mess present
4273. Apal.mess.pres <- informed.mess(ref.extent = Holo.back,
4274. mess.extent = pres.back, #maybe they are crazy high because project to same, lets try diff
4275. coord.cols = 2:3,
4276. predic = 6:9,
4277. tolerance = Apal.tolerance)
4279. ################### running mess 45 2050
4280. Apal.mess.452050 <- informed.mess(ref.extent = Holo.back,
4281. mess.extent = RCP452050.back, #maybe they are crazy high because project to same, lets try diff
4282. coord.cols = 2:3,
4283. predic = 6:9,
4284. tolerance = Apal.tolerance)
4286. ################### running mess 45 2100
4287. Apal.mess.452100 <- informed.mess(ref.extent = Holo.back,
4288. mess.extent = RCP452100.back, #maybe they are crazy high because project to same, lets try diff
4289. coord.cols = 2:3,
4290. predic = 6:9,
4291. tolerance = Apal.tolerance)
4293. ################### running mess 85 2050
4294. Apal.mess.852050 <- informed.mess(ref.extent = Holo.back,
4295. mess.extent = RCP852050.back, #maybe they are crazy high because project to same, lets try diff
4296. coord.cols = 2:3,
4297. predic = 6:9,
4298. tolerance = Apal.tolerance)
4300. ################### running mess 85 2100
4301. Apal.mess.852100 <- informed.mess(ref.extent = Holo.back,
4302. mess.extent = RCP852100.back, #maybe they are crazy high because project to same, lets try diff
4303. coord.cols = 2:3,
4304. predic = 6:9,
4305. tolerance = Apal.tolerance)
4307. #### apal mess all
4308. Apal.mess.all <- informed.mess(ref.extent = Holo.back,
4309. mess.extent = All.back, #maybe they are crazy high because project to same, lets try diff
4310. coord.cols = 2:3,
4311. predic = 6:9,
4312. tolerance = Apal.tolerance)
4316. ##### Saving Everything ----------------------------------------------------------------------------------------------------
4317. apal.output <- All.back %>% dplyr::select('Name', 'long', 'lat', 'time.bin')
4318. apal.output <- cbind(apal.output, ApalHolo.everything_LQP0.025)
4319. apal.output$Name <- 'Apalmata'
4320. apal.with.mess <- cbind(apal.output, Apal.mess.all)
4321. 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)
4323. # null model
4324. write.csv(Apal.null, 'Apal.null.csv', row.names = F)
4325. # model optimized summary object
4326. write.csv(Apal.optim, 'Apal.optim.summary.csv', row.names = F)
4327. # model projected to all extents
4328. write.csv(ApalHolo.everything_LQP0.025, 'Apal.LQP0.025.csv', row.names = F)
4329. # model mess analysis
4330. write.csv(Apal.mess.pres, 'Apal.mess.pres.csv', row.names = F)
4331. write.csv(Apal.mess.452050, 'Apal.mess.452050.csv', row.names = F)
4332. write.csv(Apal.mess.452100, 'Apal.mess.452100.csv', row.names = F)
4333. write.csv(Apal.mess.852050, 'Apal.mess.852050.csv', row.names = F)
4334. write.csv(Apal.mess.852100, 'Apal.mess.852100.csv', row.names = F)
4338. ############################## CNAT PRES ENM ###########
4340. ##### Loading in the Data --------------------------------------------------------------------------------------------------
4342. # master background file;
4343. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/Master_files_cnat_add_mask")
4345. All.back <- read.csv('masterall.background.csv', header=T)
4347. #filter to time period background
4348. pres.back <- All.back %>% filter(time.bin== c('Present'))
4350. Holo.back <- All.back %>% filter(time.bin=='Holocene')
4351. RCP452050.back <- All.back %>% filter(time.bin=='RCP452050')
4352. RCP852050.back <- All.back %>% filter(time.bin=='RCP852050')
4353. RCP452100.back <- All.back %>% filter(time.bin=='RCP452100')
4354. RCP852100.back <- All.back %>% filter(time.bin=='RCP852100')
4356. #save background I am using
4357. 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)
4360. # master occurrence file;
4361. all.occ <- read.csv('masterall.occs.csv', header=T)
4363. #reading in specific time period files
4364. pres.occs <- read.csv('masterpres.occs.csv', header=T)
4366. #check if background and occurance have same column names
4367. identical(colnames(pres.back), colnames(pres.occs))
4370. ##### Setting up the models ------------------------------------------------------------------------------------------------
4372. # working directory for models
4373. setwd("/Users/clairewilliams/Library/CloudStorage/Box-Box/Grad_school/ENM/Jone_code/CNAT/masked/models_Cnat_pres_add_mask")
4376. # create null df - empty dataframe that null data will go into later
4377. pres.null <- create.null.df('Cnatans', 'Present', 'Caribbean') #function(taxon.name, time.bin, extent){
4379. # create summary df
4380. pres.summary <- create.summary.df('Cnatans', 'Present', 'Caribbean')
4382. create.folders.for.maxent(pres.summary)
4384. # run null model
4385. Cnat.null <- null.aic(null.df = pres.null,
4386. occs = pres.occs,
4387. background = pres.back,
4388. first.occ.col = 10)
4391. ##### Optimizing Model Parameters ------------------------------------------------------------------------------------------
4393. # find most optimized model
4394. Cnat.optim <- optimize.maxent.likelihood(pres.summary, # name of the output summary file
4395. occs = pres.occs, # species occurrences
4396. background = pres.back, # background for the pres
4397. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
4398. first.occ.col = 10)
4400. View(Cnat.optim)
4401. #want model to be at least more than 2 AICc lower than the null
4402. #also plan to check pROC and response curves to pick best model
4404. 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)
4405. ### CHECK RESPONSE CURVES FOR REALISM!!!!!
4407. ##### Running Cross-Validation ---------------------------------------------------------------------------------------------
4409. ### LQP 0.1
4410. # create eval object
4411. # default f.class = LQP, so if f.class argument is missing, will assume LQP; only need it there in case NOT LQP
4412. Cnatpres.eval_LQP0.1 <- create.eval.df('Cnat', 'Present', 'Caribbean', beta.values = 0.1, f.class = 'LQP')
4413. Cnatpres.eval_LQP0.1
4415. # create folders for Baculites in the Cenomanian
4416. create.folders.for.maxent(Cnatpres.eval_LQP0.1)
4418. # run eval object
4419. Cnatpres.eval_LQP0.1 <- maxent.crossval.error(eval.df = Cnatpres.eval_LQP0.1,
4420. occs = pres.occs, # species occurrences
4421. background = pres.back, # background for the pres
4422. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
4423. first.occ.col = 10,
4424. first.test.col = 16,
4425. omission.rate = 0, #express as proportion
4426. all.background = All.back)
4429. #these results go into the folders that are only 4 of the 5 abcde (ex abcd, acde etc.)
4430. # response curves good, not all are complete humps
4431. ### LQP 0.025
4432. Cnatpres.eval_LQP0.025 <- create.eval.df('Cnat', 'Present', 'Caribbean', beta.values = 0.025, f.class = 'LQP')
4433. Cnatpres.eval_LQP0.025
4435. # create folders
4436. create.folders.for.maxent(Cnatpres.eval_LQP0.025)
4438. # run eval object
4439. Cnatpres.eval_LQP0.025 <- maxent.crossval.error(eval.df = Cnatpres.eval_LQP0.025,
4440. occs = pres.occs, # species occurrences
4441. background = pres.back, # background for the pres
4442. predic = 6:9, # column numbers of the predictor variables (maxsalinity, maxtemp, rangesalinity, rangetemp)
4443. first.occ.col = 10,
4444. first.test.col = 16,
4445. omission.rate = 0, #express as proportion
4446. all.background = All.back)
4448. #response curves better but quite similar, more clear humps
4450. ### analyzing summary files
4451. ## if running multiple eval models, check for:
4452. # 1.) does one model setting have systematically higher omission rates?
4453. # 2.) do the lambdas of one model change a lot more between cross-validation reps than another model type?
4454. # 3.) do the response curves of one model look non-sensical? check all 5 cross-validation reps
4456. Cnatpres.eval_LQP0.1$summary # pROC 1.477102
4457. Cnatpres.eval_LQP0.025$summary # pROC 1.446531
4460. # going with LQP0.025 because response curves are better
4463. ##### Evaluate Model, Calculate Threshold and Model Variability ------------------------------------------------------------
4465. ### LQP 0.025
4466. # evaluate model to calculate weighted suitibility and stdev
4467. Cnatpres.means_LQP0.025 <- maxent.eval(eval = Cnatpres.eval_LQ

COBI-40-e70323-s014.R, no license · at the source