Internal validation of the quadratic grampians model

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

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000644 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.44 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 24.2374 seconds (Warm-up)
Chain 1:                8.41223 seconds (Sampling)
Chain 1:                32.6496 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000509 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 5.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 21.9177 seconds (Warm-up)
Chain 2:                8.39957 seconds (Sampling)
Chain 2:                30.3172 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000476 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.76 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 24.1555 seconds (Warm-up)
Chain 3:                8.38095 seconds (Sampling)
Chain 3:                32.5365 seconds (Total)
Chain 3: 

library(arm)
library(forcats)
library(ggrepel)
library(tidyr)

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

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

ggplot(spread(gramps_auc, random_effect, AUC)) +
aes(`FALSE`, `TRUE`, label = taxon) +
geom_abline(slope = 1, intercept = 0, color = "grey") +
geom_point() +
geom_text_repel(fontface = "italic", size = 3) +
xlab("AUC (no random effect)") +
ylab("AUC (random effect included)") +
xlim(.4, 1) +
ylim(.4, 1) +
theme_bw()

vars <- c("mlq", "mlq2", "twi", "twi2","r1k", "r1k2", "tn", "tn2")
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(
  grampians_responses[order(grampians_responses$AUC), ]
) +
xlim(-3, 3) +
ylim(-3, 3) +
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, 4) +
theme_bw()