Validation of linear model with location random effect

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

The following is a validation of the predictions/inference from a trait-environment model of the Grampians Bioregion ~20 taxa dataset.

The following steps are undertaken:

• Fit a model to Grampians data.
• Calculate Performance Metrics:
  • Predict probability of occurrence
    • Predict within Grampians training set:
      • Predict for each species including species random effects.
      • Predict for each species without species random effects.
    • Predict to the South-east validation set:
      • Predict to each IBRA subregion within the South-east:
        • Predict for each species common to Grampians and South-east subregion:
          • Predict including species random effects.
          • Predict without species random effects.
        • Predict for each species not found in the Grampians training set.
  • For each prediction type (within-sample/out-of-sample, with-re/without-re, and modelled-species/unmodelled-species), bioregion, and species (and their combinations) calculate performance metrics

mds <- local({
  nms1 <-   c(
    "mean_moisture_index_of_low_qtr", "topographic_wetness_index_3sec",
    "relief_1000m_radius", "total_nitrogen"
  )
  nms2 <- c("sla_mm2_per_mg", "seed_mass_mg", "max_height_m")
  vars <- lapply(nms1, function(x) log(grampians[[x]]))
  vars <- do.call(cbind, vars)
  vars <- cbind(
    vars, do.call(cbind, lapply(nms2, function(x) log(grampians[[x]])))
  )
  vars <- scale(vars)
  cntrs <- attr(vars, "scaled:center")
  scales <- attr(vars, "scaled:scale")

  vars_new_raw <- na.omit(southeast)
  vars_new <- lapply(seq_along(nms1), function(x) log(vars_new_raw[[nms1[x]]]))
  vars_new <- do.call(cbind, vars_new)
  vars_new <- cbind(
    vars_new, do.call(cbind, lapply(nms2, function(x) log(vars_new_raw[[x]])))
  )
  vars_new <- scale(vars_new, center = cntrs, scale = scales)

  setNames(
    data.frame(
      vars_new_raw[, c("occupancy", "taxon", "ibra_subregion")], vars_new
    ),
    c(
      "y", "taxon", "ibra_subregion", "mlq", "twi", "r1k", "tn", "sla", "sm",
      "mh"
    )
  )
})

perf <- local({
  perf <- mapply(
    function(f, re.form) {
      f(
        model_grampians, y, list(taxon, ibra_subregion),
        rbind(mds, mdg[names(mds)]), re.form
      )
    },
    f = rep(list(auroc, pred_inf, dev_explained), each = 2),
    re.form = list(~(mlq + twi + r1k + tn | taxon),  NA)
  )
  taxon <- sub("(.*)\\..*", "\\1", rownames(perf))
  ibra_subregion <- sub(".*\\.", "", rownames(perf))
  dim(perf) <- dim(perf) * c(2, .5)
  colnames(perf) <- c("AUC", "mean_predictive_info", "deviance_explained")
  cbind(
    taxon = rep(taxon, 2),
    ibra_subregion = rep(ibra_subregion, 2),
    random_effect = rep(c(TRUE, FALSE), each = nrow(perf) / 2),
    as.data.frame(perf)
  )
})

ggplot(
  transform(
    subset(perf, ibra_subregion != "Greater Grampians" & !random_effect),
    ibra_subregion = factor(
      ibra_subregion, names(sort(tapply(AUC, ibra_subregion, median)))
    )
  )
) +
aes(AUC) + geom_histogram() + facet_wrap(~ibra_subregion, 3) +
xlab("AUC") +
ylab("No. of taxa") +
ggtitle(
  paste(
    "Area Under Receiver Operator Curve for taxa within IBRA subregions of",
    "Southeast Australia"
  )
)

ggplot(
  transform(
    subset(perf, ibra_subregion != "Greater Grampians" & !random_effect),
    ibra_subregion = factor(
      ibra_subregion,
      names(sort(tapply(mean_predictive_info, ibra_subregion, median)))
    )
  )
) +
aes(mean_predictive_info) + geom_histogram() + facet_wrap(~ibra_subregion, 3) +
xlab("Mean Predictive Information") +
ylab("No. of taxa") +
ggtitle(
  paste(
    "Mean Predictive Information for taxa within IBRA subregions of",
    "Southeast Australia"
  )
)

ggplot(
  transform(
    subset(perf, taxon %in% mdg$taxon & random_effect),
    taxon = factor(
      taxon,
      names(
        sort(
          tapply(
            AUC[ibra_subregion == "Greater Grampians"],
            taxon[ibra_subregion == "Greater Grampians"],
            median
          )
        )
      )
    )
  )
) +
aes(
  AUC,
  mean_predictive_info,
  col = ifelse(
    ibra_subregion == "Greater Grampians",
    "Grampians",
    "Other Region"
  )
) +
geom_point() +
facet_wrap(~taxon, labeller = label_wrap_gen(40)) +
labs(col = NULL) +
xlab("AUC") +
ylab("Mean Predictive Information") +
ggtitle(
  paste(
    ""
  )
) +
theme(legend.position = "bottom")