Validating the species responses for grampians trait-environment model with quadratic relationships

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

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000645 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.45 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: 21.7004 seconds (Warm-up)
Chain 1:                8.58763 seconds (Sampling)
Chain 1:                30.2881 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000918 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.18 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: 26.6919 seconds (Warm-up)
Chain 2:                8.89649 seconds (Sampling)
Chain 2:                35.5884 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000529 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 5.29 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: 29.9949 seconds (Warm-up)
Chain 3:                9.73212 seconds (Sampling)
Chain 3:                39.7271 seconds (Total)
Chain 3: 

library(arm)
library(forcats)
library(ggplot2)
library(gridExtra)
library(tidyr)

vars <- c("mlq", "mlq2", "twi", "twi2", "r1k", "r1k2", "tn", "tn2")
taxa <- unique(mds$taxon)
taxa_x_region <- unique(mds[, c("taxon", "ibra_subregion")])
traits <- sapply(
  taxa,
  function(x) unlist(unique(mds[mds$taxon == x, c("sla", "sm", "mh")]))
)

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

observed_responses <- mapply(
  function(taxoni, region) {
    df <- subset(mds, taxon == taxoni & ibra_subregion == region)
    model <- bayesglm(
      formula(sprintf("y ~ %s", paste(vars, collapse = " + "))),
      binomial, df, prior.scale = 1, prior.df = Inf
    )
    data.frame(
      Taxon = taxoni, Region = region, Var = vars,
      summary(model)$coefficients[vars, 1:2]
    )
  },
  taxa_x_region[, 1], taxa_x_region[, 2],
  SIMPLIFY = FALSE
)

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)[taxa_x_region[, 1]]
)

responses <- cbind(observed_responses, predicted_responses)

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

plot_regions <- c(
  "Snowy Mountains", "Jervis", "Victorian Alps", "Strzelecki Ranges",
  "Greater Grampians"
)

responses <- rbind(responses, grampians_responses)

responses$Region <- forcats::fct_relevel(responses$Region, plot_regions)

scatters <- ggplot(
    subset(
      responses[order(responses$AUC),],
      Var %in% c("mlq", "twi") & Region %in% plot_regions
    )
  ) +
  xlim(-2, 6) +
  ylim(-7, 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_grid(
    Region ~ Var,
    labeller = labeller(
      Var = as_labeller(
        c(
          mlq  = "Moisture Lowest Quarter",
          twi  = "Topographic Wetness"
        )
      )
    )
  ) +
  theme_bw() +
  theme(
    legend.position = c(0.35, .89),
    legend.key.width = unit(.07, "inches"),
    legend.key.height = unit(.09, "inches"),
    legend.background = element_rect(fill = "transparent"),
    legend.title = element_text(size = 8),
    legend.text = element_text(size = 6),
    strip.background = element_blank(),
    strip.text.y = element_blank()
  )

perf$ibra_subregion <- forcats::fct_relevel(perf$ibra_subregion, plot_regions)

hist_data <- subset(perf, ibra_subregion %in% plot_regions & !random_effect)
med_data <- with(hist_data, tapply(AUC, ibra_subregion, median))
med_data <- as.data.frame(med_data)
med_data <- na.omit(med_data)
med_data$ibra_subregion <- row.names(med_data)
med_data$ibra_subregion <- fct_relevel(med_data$ibra_subregion, plot_regions)

hists <- ggplot(hist_data) +
  aes(AUC) +
  geom_histogram() +
  facet_grid(ibra_subregion ~ random_effect) +
  geom_vline(aes(xintercept = med_data), med_data, col = "grey", size = 2) +
  xlab("AUC") +
  ylab("No. of taxa") +
  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text.x = element_text(color = "transparent")
  )

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

grid.arrange(scatters, hists, ncol = 2, widths = c(.45, .55))

Plot of the observed vs. predicted responses to TWI grey-coded by trait values.

ggplot(
  merge(
    subset(
      responses,
      Var == "twi" & Region %in% plot_regions[c(1, 3, 5)]
    ),
    gather(
      unique(
        rbind(
          mdg[c("taxon", "ibra_subregion", "sla", "sm")],
          mds[c("taxon", "ibra_subregion", "sla", "sm")]
        )
      ),
      "trait",
      "value",
      sla, sm
    ),
    by.x = c("Taxon", "Region"), by.y = c("taxon", "ibra_subregion")
  )
) +
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 = value), size = 2) +
scale_color_gradient(
  name = "SD", low = "white", high = "black"
) +
facet_grid(
  Region ~ trait,
  labeller = labeller(
    trait = as_labeller(
      c(
        sla = "Specific Leaf Area",
        sm  = "Seed Mass"
      )
    )
  )
) +
theme_bw() +
theme(
  legend.position = c(0.07, .91),
  legend.key.width = unit(.07, "inches"),
  legend.key.height = unit(.09, "inches"),
  legend.background = element_rect(fill = "transparent"),
  legend.title = element_text(size = 8),
  legend.text = element_text(size = 6),
  strip.background = element_blank()
)

# Calculate the correlation between predicted vs observed response coefficients
# for each covariate among all taxa within a region.
# calccor
# correlations <- lapply(
#   vars,
#   function(var) {
#     ans <- sapply(
#       rgn_names,
#       function(region) {
#         ss <- subset(responses, Region == region & Var == var)
#         cor(ss$Estimate, ss$Prediction)
#       }
#     )
#     data.frame(
#       var2 = var, region = names(ans), cor = ans, stringsAsFactors = FALSE,
#       row.names = NULL
#     )
#   }
# )
#
# correlations <- do.call(rbind, correlations)
#
# correlations <- data.frame(
#   correlations,
#   var = gsub("2", "", correlations$var),
#   KLD = kldists[correlations$region],
#   GD = dists[correlations$region],
#   stringsAsFactors = FALSE
# )
#
# correlations <- merge(correlations, kldists_univariate)
#
# Plot the correlations vs. the KL distance of the regions environment from the
# Grampians.
# plotcors1, fig.width=11
# ggplot(correlations) +
# aes(KLD, cor) +
# geom_point() +
# geom_smooth(method = "lm") +
# xlim(0, 35) +
# ylim(-1, 1) +
# facet_wrap("var2", 2, dir = "v")
#
# Plot the correlations vs. the univariate KL distance of the regions
# environment from the Grampians.
# plotcors2, fig.width=11
# ggplot(correlations) +
# aes(kldist_var, cor) +
# geom_point() +
# geom_smooth(method = "lm") +
# ylim(-1, 1) +
# facet_wrap("var2", 2, scales = "free_x", dir = "v")
#
# Plot the correlations vs. the geographic (centroid) distance of the region
# to the Grampians.
# plotcors3, fig.width=11
# ggplot(correlations) +
# aes(GD, cor) +
# geom_point() +
# geom_smooth(method = "lm") +
# xlim(0, 1000) +
# ylim(-1, 1) +
# facet_wrap("var2", 2, dir = "v")