Predicting niche optima from trait-environment models

This noteboook tests some ideas for finding the niche optima (a vector) for a given set of traits from a trait environment model.

knitr::opts_chunk$set(cache = TRUE, dev = c("png", "pdf"))
library(here)
source(here("src", "model_linear_validation_grampians.R"))

Warning: Column `taxon` joining factor and character vector, coercing into
character vector

Warning: Column `ibra_subregion` joining factor and character vector, coercing
into character vector

Warning: Removed 60 rows containing missing values (geom_point).

Warning: Removed 3 rows containing missing values (geom_point).

Warning: Removed 60 rows containing missing values (geom_point).

Warning: Removed 3 rows containing missing values (geom_point).

Warning: Removed 30 rows containing missing values (geom_point).

Warning: Removed 3 rows containing missing values (geom_point).

Warning: Removed 30 rows containing missing values (geom_point).

Warning: Removed 3 rows containing missing values (geom_point).

Warning: Column `taxon` joining factor and character vector, coercing into
character vector

Warning: Column `ibra_subregion` joining factor and character vector, coercing
into character vector

library(geometry) # use dev version from github
library(ks) # use own fork with bugfix in kde

---

Pick a region and a taxon

region <- "South East Coastal Ranges"
taxon <- "Angophora costata subsp. costata"
traits <- unique(mds[mds$taxon == taxon, c("sla", "sm", "mh")])
vars <- c("mlq", "twi", "r1k", "tn")
environ <- as.matrix(
  mds[mds$ibra_subregion == region & mds$taxon == taxon, c("y", vars)]
)
environ_extent <- apply(environ[, vars], 2, range)
gridsize <- rep(30, length(vars))

---

Calculate the kernel density of the absence and presence points separately.

d <- lapply(
  0:1,
  function(x) {
    kde(
      environ[environ[, "y"] == x, vars],
      xmin = as.vector(environ_extent[1, ]),
      xmax = as.vector(environ_extent[2, ]),
      gridsize = gridsize
    )
  }
)

---

Calculate a convex hull for the environmental space.

chull <- convhulln(environ[, vars])

---

Calculate the empirical niche optima from the presence absence data. The optima is the point inside the hull that maximises the ratio of kernel densities of presence and absence.

empirical_niche_optima <- arrayInd( # get the indices
  which.max(
    d[[2]]$estimate / d[[1]]$estimate * # ratio of kernel densities
      # check if inside hull
      inhulln(chull, as.matrix(do.call(expand.grid, d[[1]]$eval.points)))
  ),
  gridsize
)

# extract values from the grid
empirical_niche_optima <- sapply(
  1:4, function(i) d[[1]]$eval.points[[i]][empirical_niche_optima[, i]]
)

names(empirical_niche_optima) <- vars

---

Convert convex hull object to a matrix of hull vertices

chull <- unique(t(sapply(chull, function(x) environ[x, vars])))

---

Calculate the predicted niche optima based on the model by finding which vertex of the environment convex hull has the highest predicted probability of occurrance.

predicted_niche_optima <- as.vector(
  chull[
    which.max(
      predict(
        model_grampians,
        data.frame(traits[rep(1, nrow(chull)), ], chull, row.names = NULL),
        re.form = ~0
      )
    ),
  ]
)

names(predicted_niche_optima) <- vars

environ_extent

          mlq       twi       r1k        tn
[1,] 1.698242 -2.356902 -1.747737 -1.608506
[2,] 5.519156  2.507559  1.587113  3.089763

predicted_niche_optima

     mlq      twi      r1k       tn 
4.611756 1.052513 1.384686 1.193711 

empirical_niche_optima

        mlq         twi         r1k          tn 
 3.54282107 -1.18272135 -0.82777814  0.01158644