Niche optima validation of linear 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 for inhulln()
library(ks) # use own fork with bugfix in kde
library(ggplot2)

vars <- c("mlq", "twi", "r1k", "tn")
regions <- unique(mds$ibra_subregion)
taxa <- unique(mds$taxon)
taxa_x_region <- unique(mds[, c("taxon", "ibra_subregion")])

traits <- sapply(
  taxa,
  function(x) unique(mds[mds$taxon == x, c("sla", "sm", "mh")]),
  simplify = FALSE
)

environs <- sapply(
  unique(mds$ibra_subregion),
  function(x) {
    as.matrix(
      mds[
        mds$ibra_subregion == x &
        mds$taxon ==
          taxa_x_region[taxa_x_region$ibra_subregion == x, "taxon"][1],
        c("y", vars)
      ]
    )
  }
)

environ_extents <- sapply(
  environs, function(x) apply(x[, vars], 2, range), simplify = FALSE
)

chulls <- sapply(environs, function(x) convhulln(x[, vars]))

empirical_niche_optima <- function(taxon, region) {
  environ <- as.matrix(
    mds[mds$ibra_subregion == region & mds$taxon == taxon, c("y", vars)]
  )
  environ_extent <- environ_extents[[region]]
  gridsize <- rep(30, length(vars))

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

  if (class(d[[2]]) == "try-error") return(NA)

  chull <- chulls[[region]]

  empirical_niche_optima <- arrayInd(
    which.max(
      d[[2]]$estimate / d[[1]]$estimate *
        inhulln(chull, as.matrix(do.call(expand.grid, d[[1]]$eval.points)))
    ),
    gridsize
  )

  empirical_niche_optima <- sapply(
    1:4, function(i) d[[1]]$eval.points[[i]][empirical_niche_optima[, i]]
  )

  names(empirical_niche_optima) <- vars
  empirical_niche_optima
}

file <- paste0(eucs_cache_path(), "/empirical_nich_optima.rds")

if (!file_test("-f", file)) {
  empirical_niche_optima <- parallel::mcmapply(
    empirical_niche_optima, taxa_x_region[, 1], taxa_x_region[, 2],
    mc.cores = 3
  )
  saveRDS(empirical_niche_optima, file, compress = "xz")
}  else {
  empirical_niche_optima <- readRDS(file)
}

chull_vertices <- mapply(
  function(x, y) unique(t(sapply(x, function(z) y[z, vars]))),
  chulls,
  environs
)

predicted_niche_optima <- function(taxon, region) {
  traits <- traits[[taxon]]
  chull_verts <- chull_vertices[[region]]

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

  names(predicted_niche_optima) <- vars
  predicted_niche_optima
}

predicted_niche_optima <- mapply(
  predicted_niche_optima, taxa_x_region[, 1], taxa_x_region[, 2],
  SIMPLIFY = FALSE
)

rand_dist <- mapply(
  function(x, y) sqrt(sum((x - y)^2)),
  lapply(
    taxa_x_region$ibra_subregion,
    function(x) {
      ans <- rep(-999, 4)
      while (!inhulln(chulls[[x]], rbind(ans))) {
        ans <- runif(
          4,
          as.vector(environ_extents[[x]][1, ]),
          as.vector(environ_extents[[x]][2, ])
        )
      }
      ans
    }
  ),
  empirical_niche_optima
)

max_dist <- sapply(chull_vertices, function(x) max(dist(x)))
names(max_dist) <- regions

distance <- data.frame(
  taxa_x_region,
  relative_distance = mapply(
    function(x, y) sqrt(sum((x - y)^2)),
    predicted_niche_optima,
    empirical_niche_optima
  ) / max_dist[taxa_x_region[[2]]],
  random_distance = rand_dist / max_dist[taxa_x_region[[2]]]
)

ggplot(distance) + aes(relative_distance) + geom_histogram() +
  facet_wrap(~ibra_subregion)

ggplot(distance) + aes(random_distance) + geom_histogram() +
  facet_wrap(~ibra_subregion)