1 Preparing packages and loading data

As always, we need to load all packages we will be working with. In addition to the packages from last session, we will also need the rJava package to be able to interact with the Maxent (that works on java).

# load packages
library(sf)
library(terra)
library(mapview)
library(predicts)
library(tidyverse)
library(rJava)

Load the data (presence points and predictors rasters) with the same functions as in the sloth exercise.

python_presence <- read.table("data_python/python_presence.txt", header = TRUE)

current_asia <- terra::rast("data_python/current_asia.grd")
current_america <- terra::rast("data_python/current_america.grd")
future45_america <- terra::rast("data_python/future45_america.grd")
future85_america <- terra::rast("data_python/future85_america.grd")

2 Running Maxent

To create the “model_df” data.frame needed for building Maxent models, you have to (1) sample background points and (2) extract environmental predictor values at presence and background locations, as demonstrated in the .html file for the sloth exercise. Adapt the code for the new python dataset before continuing with building the maxent models.

We will use the predicts package to build Maxent models. To understand how the MaxEnt() function works, first check out its help file (?predicts::MaxEnt).

Similar to building a GLM, we need to define specify whether and observation is a presence or background point and select which variables to include in the model. Instead of using the formula interface, we will give the function (1) the modeling data frame containing only the predictor variables to be used (x argument) and (2) a vector corresponding to the occurrence column from your modeling data frame (i.e., 1s and 0s; p argument).

We can additionally specify parameter settings via the args argument to select feature types to be used and set the betamultiplier to control the amount of regularization (higher values will result in more restricted and less complex models). Check out the hidden code section below for an example.

Once you have fit the model, you can plot response curves with the response() function (no custom prediction function needed for Maxent models!). In addition, you can check out a summary file generated by the Maxent software by applying the show() function to the Maxent model object, which gives you an overview of the model performance, potential thresholds for creating binary habitat maps and the relative importance of predictor variables in the model.

Test how different settings of the betamultiplier and allowed feature classess change the response curves of the model!

To predict a Maxent model to the study area using the stack of environmental predictors, we can use the predict() method from the raster package (same approach as for GLMs, no extra argument needed to specify the prediction scale).

predictors <- c("annual_mean_temp", "temp_min_cold")

maxent_args <-c('linear=true', 
          'quadratic=true',
          'product=true',
          'hinge=true', 
          'threshold=false',
          'betamultiplier=1')

maxent_model <- predicts::MaxEnt(x = model_df[,predictors], 
                          p = model_df$occ,
                         args = maxent_args)

predicts::partialResponse(maxent_model)

pred_maxent_asia <- terra::predict(current_asia, maxent_model)
## |---------|---------|---------|---------|=========================================                                          
plot(pred_maxent_asia)

train_ids <- sample(1:nrow(model_df), size = 0.8 * nrow(model_df)) 


training_data <- model_df[train_ids, ]   
testing_data  <- model_df[-train_ids, ]  

e <- predicts::pa_evaluate(
  p = testing_data[testing_data$occ == 1, ],  
  a = testing_data[testing_data$occ == 0, ],  
  maxent_model                                
)

plot(e, "ROC")