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

knitr::opts_chunk$set(cache = TRUE, dev = c("png", "pdf"))
library(here)
source(here("src", "model_linear_grampians_internal_validation.R"))
library(arm)
library(forcats)
library(ggplot2)
library(gridExtra)
library(tidyr)

vars <- c("mlq", "twi", "r1k", "tn")
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)

taxon_region_prevalence <- summarise(
  group_by(rbind(mds, mdg[names(mds)]), taxon, ibra_subregion),
  prevalence = mean(y)
)

responses <- left_join(
  responses, taxon_region_prevalence,
  by = c("Taxon" = "taxon",  "Region" = "ibra_subregion")
)

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

perf <- left_join(perf, taxon_region_prevalence)

Plot of the observed vs. predicted responses for each covariate grey-coded by region (AUPRC/Prevalence).

scatters <- ggplot(
    subset(
      responses[order(responses$`AUC-PR` / responses$prevalence), ],
      Var %in% c("mlq", "twi") & Region %in% plot_regions
    )
  ) +
  xlim(-3, 3) +
  ylim(-3, 3) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_abline(slope = 1, intercept = 0, lty = 2) +
  aes(Prediction, Estimate) +
  geom_point(aes(color = `AUC-PR` / prevalence), size = 2) +
  scale_color_gradient(
    name = "AUPRC /\nPrevalence", low = "white", high = "black",
    trans = "log2"
  ) +
  facet_grid(
    Region ~ Var,
    labeller = labeller(
      Var = as_labeller(
        c(
          mlq  = "Moisture Lowest Quarter",
          twi  = "Topographic Wetness"
        )
      )
    )
  ) +
  xlab("Trait-SDM prediction") +
  ylab("Taxon regression estimate") +
  theme_bw() +
  theme(
    legend.position = c(0.11, .94),
    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()
  )

hist_data <- subset(perf, ibra_subregion %in% plot_regions & !random_effect)
med_data <- with(
  hist_data, tapply(`AUC-PR` / prevalence, 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-PR` / prevalence) +
  geom_histogram() +
  scale_x_continuous(trans = "log2") +
  facet_grid(ibra_subregion ~ random_effect) +
  geom_vline(aes(xintercept = med_data), med_data, col = "grey", size = 2) +
  xlab("AUPRC / Prevalence") +
  ylab("No. of taxa") +
  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text.x = element_text(color = "transparent")
  )

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

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

scatters_roc <- ggplot(
    subset(
      responses[order(responses$AUC), ],
      Var %in% c("mlq", "twi") & Region %in% plot_regions
    )
  ) +
  xlim(-3, 3) +
  ylim(-3, 3) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_abline(slope = 1, intercept = 0, lty = 2) +
  aes(Prediction, Estimate) +
  geom_point(aes(color = AUC), size = 2) +
  scale_color_gradient(
    name = "AUROC", low = "white", high = "black"
  ) +
  facet_grid(
    Region ~ Var,
    labeller = labeller(
      Var = as_labeller(
        c(
          mlq  = "Moisture Lowest Quarter",
          twi  = "Topographic Wetness"
        )
      )
    )
  ) +
  xlab("Trait-SDM prediction") +
  ylab("Taxon regression estimate") +
  theme_bw() +
  theme(
    legend.position = c(0.09, .95),
    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()
  )

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

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

grid.arrange(scatters_roc, hists_roc, 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(Prediction, Estimate) +
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()
)

Plot of the observed vs. predicted responses for each covariate

ggplot(responses) +
aes(
  Prediction, Estimate,
  fill = Region == "Greater Grampians",
  alpha =  Region == "Greater Grampians"
) +
geom_abline(slope = 1, intercept = 0, col = "grey") +
geom_point(col = "black", pch = 21, size = 2) +
scale_fill_grey(start = 0, end = 1) +
scale_alpha_discrete(range = c(0.3, 1)) +
facet_wrap(
  ~Var, 2,
  labeller = labeller(
    Var = as_labeller(
      c(
        mlq = "Moisture Lowest Quarter",
        twi = "Topographic Wetness",
        r1k = "Relief 1000m Radius",
        tn  = "Total Nitrogen"
      )
    )
  )
) +
xlab("Trait-SDM prediction") +
ylab("Taxon regression estimate") +
lims(x = c(-4, 4), y = c(-4, 4)) +
geom_text(
  data    = cbind(
    summarise(group_by(responses, Var), r = cor(Prediction, Estimate)),
    Region = ""
  ),
  mapping = aes(
    x = -Inf, y = Inf,
    label = sprintf("italic('r') == %s", round(r, digits = 2))
  ),
  hjust   = -0.1,
  vjust   = 1.5,
  parse = TRUE
) +
theme_bw() +
theme(
  legend.position = "none",
  strip.background = element_blank()
)

Plot of the miscalibration for each covariate vs. performance

ggplot(responses) +
aes(
  Prediction - Estimate, `AUC-PR`/prevalence,
  fill = Region == "Greater Grampians",
  alpha =  Region == "Greater Grampians"
) +
geom_point(col = "black", pch = 21, size = 2) +
scale_y_log10() +
scale_fill_grey(start = 0, end = 1) +
scale_alpha_discrete(range = c(0.3, 1)) +
facet_wrap(
  ~Var, 2,
  labeller = labeller(
    Var = as_labeller(
      c(
        mlq = "Moisture Lowest Quarter",
        twi = "Topographic Wetness",
        r1k = "Relief 1000m Radius",
        tn  = "Total Nitrogen"
      )
    )
  )
) +
xlim(-5, 5) +
ylab("AUPRC / Prevalence") +
xlab("Miscalibration (Prediction - Estimate)") +
theme_bw() +
theme(
  legend.position = "none",
  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(
#      var = var, region = names(ans), cor = ans, stringsAsFactors = FALSE,
#      row.names = NULL
#    )
#  }
#)

#correlations <- do.call(rbind, correlations)

#correlations <- data.frame(
#  correlations,
#  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("var")

# 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("var", scales = "free_x")

# 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("var")