## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 10.5,
  fig.height = 6.5,
  out.width = "100%",
  fig.align = "center"
)

species_example <- requireNamespace("disdat", quietly = TRUE) &&
  requireNamespace("randomForest", quietly = TRUE)

ensemble_example <- species_example &&
  requireNamespace("mgcv", quietly = TRUE)

## ----eval = species_example, message = FALSE, warning = FALSE-----------------
spID <- "nsw14"
pr <- disdat::disPo("NSW")
bg <- disdat::disBg("NSW")

pr <- pr[pr$spid == spID, ]
training <- rbind(pr, bg)

covars <- c(
  "disturb", "mi", "rainann", "raindq",
  "rugged", "soildepth",
  "tempann", "topo", "vegsys"
)

training <- training[, c("occ", covars)]
training$vegsys <- as.factor(training$vegsys)

numeric_covars <- setdiff(covars, "vegsys")
training[numeric_covars] <- lapply(
  training[numeric_covars],
  function(x) as.numeric(scale(x))
)

head(training[, 1:5])

## ----eval = species_example, message = FALSE, warning = FALSE-----------------
# get the number of presences (before converting to factors);
pr_num <- sum(training$occ)
smpsize <- c("0" = pr_num, "1" = pr_num)

training$occ <- as.factor(training$occ)

set.seed(20260501)
rf_model <- randomForest::randomForest(
  occ ~ .,
  data = training,
  ntree = 500,
  sampsize = smpsize,
  replace = TRUE
)

## ----eval = species_example, fig.height = 8, message = FALSE, warning = FALSE----
predictors <- training[, covars]

curves::univariate(
  rf_model,
  predictors,
  method = "profile",
  type = "prob",
  response = "1",
  ylab = "Relative likelihood"
)

## ----eval = species_example, fig.height = 8, message = FALSE, warning = FALSE----
set.seed(20260501)
curves::univariate(
  rf_model,
  predictors,
  method = "pdp",
  n = 100,
  background_n = 200,
  interval = "quantile",
  interval_level = 0.8,
  type = "prob",
  response = "1",
  ylab = "Relative likelihood"
)

## ----eval = species_example, fig.height = 8, message = FALSE, warning = FALSE----
set.seed(20260502)
curves::univariate(
  rf_model,
  predictors,
  method = "ice+pdp",
  n = 100,
  background_n = 200,
  type = "prob",
  response = "1",
  ylab = "Relative likelihood"
)

## ----eval = species_example, fig.height = 8, message = FALSE, warning = FALSE----
curves::univariate(
  rf_model,
  predictors,
  method = "ale",
  n = 40,
  type = "prob",
  response = "1",
  ylab = "Relative likelihood"
)

## ----eval = species_example, fig.width = 5.6, fig.height = 4, message = FALSE, warning = FALSE----
set.seed(20260503)
curves::bivariate(
  rf_model,
  predictors,
  pairs = c(2, 3),
  method = "pdp",
  background_n = 200,
  rug = TRUE,
  plot_type = "heatmap",
  type = "prob",
  response = "1",
  zlab = "Relative likelihood"
)

## ----eval = species_example, fig.width = 5.6, fig.height = 4, message = FALSE, warning = FALSE----
curves::bivariate(
  rf_model,
  predictors,
  pairs = c(2, 3),
  method = "ale",
  n = 40,
  plot_type = "heatmap",
  type = "prob",
  response = "1",
  zlab = "ALE"
)

## ----eval = ensemble_example, fig.width = 7, fig.height = 4.5, message = FALSE, warning = FALSE----
ensemble_covars <- c("mi", "rainann")
ensemble_training <- training[, c("occ", ensemble_covars)]

set.seed(20260504)
rf_member <- randomForest::randomForest(
  occ ~ .,
  data = ensemble_training,
  ntree = 300,
  sampsize = smpsize,
  replace = TRUE
)

gam_training <- transform(
  ensemble_training,
  occ_num = as.integer(as.character(occ))
)

gam_member <- mgcv::gam(
  occ_num ~ s(mi, k = 6) + s(rainann, k = 6),
  data = gam_training,
  family = binomial(),
  method = "REML"
)


rf_likelihood <- function(model, newdata) {
  predict(model, newdata, type = "prob")[, "1"]
}

gam_likelihood <- function(model, newdata) {
  as.numeric(stats::predict(model, newdata, type = "response"))
}

set.seed(20260505)
curves::multimodel(
  list(rf_member, gam_member),
  ensemble_training[, ensemble_covars],
  fun = list(rf_likelihood, gam_likelihood),
  method = "pdp",
  n = 60,
  background_n = 200,
  show_models = TRUE,
  ylab = "Relative likelihood"
)

## ----eval = species_example && requireNamespace("plotly", quietly = TRUE), message = FALSE, warning = FALSE----
rf_likelihood <- function(model, newdata) {
  predict(model, newdata, type = "prob")[, "1"]
}

curves::bivariate(
  rf_model,
  predictors,
  pairs = c(2, 3),
  plot_type = "surface",
  fun = rf_likelihood,
  zlab = "Relative likelihood"
)