Internal validation of the linear grampians 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(arm)
library(forcats)
library(ggrepel)
library(tidyr)

gramps_auc <- subset(
  perf, ibra_subregion == "Greater Grampians",
  select = c(taxon, random_effect, AUC, `AUC-PR`)
)

gramps_auc <- inner_join(
  mutate(gramps_auc, ibra_subregion = "Greater Grampians"),
  taxon_region_prevalence
)

gramps_auc <- mutate(gramps_auc, `AUPRC / prevalence` = `AUC-PR` / prevalence)

gramps_auc <- dplyr::select(gramps_auc, -`AUC-PR`)

gramps_auc <- rename(gramps_auc, AUROC = AUC)

gramps_auc$taxon <- fct_recode(
  gramps_auc$taxon,
  `Eucalyptus goniocalyx` = "Eucalyptus (goniocalyx|nortonii)",
  `Eucalyptus viminalis subsp. viminalis` =
    "Eucalyptus viminalis subsp. (pryoriana|viminalis)"
)

gramps_auc$taxon <- gsub("Eucalyptus ",  "", gramps_auc$taxon)

gramps_auc <- pivot_longer(
  gramps_auc, cols = starts_with("AU"), names_to = "type", values_to = "value"
)

gramps_auc <- pivot_wider(
  gramps_auc, names_from = random_effect, values_from = value
)

groc <-
  ggplot(filter(gramps_auc, type == "AUROC")) +
  aes(`FALSE`, `TRUE`, label = taxon) +
  geom_abline(slope = 1, intercept = 0, color = "grey") +
  geom_point() +
  geom_text_repel(fontface = "italic", size = 2) +
  xlab("AUROC (no random effect)") +
  ylab("AUROC (random effect included)") +
  xlim(.4, 1) +
  ylim(.4, 1) +
  theme_bw()

gprc <-
  ggplot(filter(gramps_auc, type != "AUROC")) +
  aes(`FALSE`, `TRUE`, label = taxon) +
  scale_y_continuous(trans = "log2") + scale_x_continuous(trans = "log2") +
  geom_abline(slope = 1, intercept = 0, color = "grey") +
  geom_point() +
  geom_text_repel(fontface = "italic", size = 2) +
  xlab("AUPRC / Prevalence (no random effect)") +
  ylab("AUPRC / Prevalence (random effect included)") +
  theme_bw()

grid.arrange(groc, gprc, nrow = 1, widths = c(.5, .5))

vars <- c("mlq", "twi", "r1k", "tn")
taxa <- unique(mdg$taxon)
traits <- sapply(
  taxa,
  function(x) unlist(unique(mdg[mdg$taxon == x, c("sla", "sm", "mh")]))
)

Calculate the responses of each species in each region using single-species/single-region glms.

observed_responses <- lapply(
  taxa,
  function(taxoni) {
    df <- subset(mdg, taxon == taxoni)
    model <- bayesglm(
      formula(sprintf("y ~ %s", paste(vars, collapse = " + "))),
      binomial, df, prior.scale = 1, prior.df = Inf
    )
    data.frame(
      Taxon = taxoni, Region = "Greater Grampians", Var = vars,
      summary(model)$coefficients[vars, 1:2]
    )
  }
)

observed_responses <- do.call(rbind, observed_responses)

Calculate the predicted responses of each species in each region based on the Grampians model and their trait values.

nsims <- 1000
sims <- fixef(sim(model_grampians, nsims))

predicted_responses <- function(taxon) {
  ans <- vapply(
    vars, function(x) {
      rowSums(
        sweep(sims[, grep(x, colnames(sims))], 2, c(1, traits[, taxon]), "*")
      )
    },
    numeric(nsims)
  )
  t(
    apply(
      ans, 2, function(x) c(Prediction = mean(x), "Predicted..Error" = sd(x))
    )
  )
}

predicted_responses <- do.call(
  rbind, sapply(taxa, predicted_responses, simplify = FALSE)
)

Plot of the observed vs. predicted responses for each covariate color-coded by region.

responses <- cbind(observed_responses, predicted_responses)

grampians_responses <- merge(
  responses,
  subset(perf, !random_effect & ibra_subregion == "Greater Grampians"),
  by.x = c("Taxon", "Region"), by.y = c("taxon", "ibra_subregion")
)

ggplot(
  subset(
    grampians_responses[order(grampians_responses$AUC),],
    Var %in% c("mlq", "twi")
  )
) +
xlim(-3, 3) +
ylim(-1, 2) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline(slope = 1, intercept = 0, lty = 2) +
aes(Estimate, Prediction) +
geom_point(aes(color = AUC), size = 2) +
scale_color_gradient(low = "white", high = "black") +
facet_wrap(
  ~ Var,
    labeller = labeller(
      Var = as_labeller(
        c(
          mlq  = "Mean Moisture\nIndex Lowest\nQuarter",
          twi  = "Topographic\nWetness Index"
        )
      )
    )
) +
theme_bw()