Validating the species responses for grampians trait-environment model with linear relationships and location random effect

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

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]]
)

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

responses <- cbind(observed_responses, predicted_responses)

ggplot(responses) +
aes(Prediction, Estimate, color = Region) +
xlim(-7, 7) +
ylim(-7, 7) +
geom_errorbar(
  aes(
    ymin = Estimate - Std..Error * 1.96,
    ymax = Estimate + Std..Error * 1.96
  ),
  alpha = .1
) +
geom_errorbarh(
  aes(
    xmin = Prediction - Predicted..Error * 1.96,
    xmax = Prediction + Predicted..Error * 1.96
  ),
  alpha = .1
) +
geom_point(alpha = .5) +
facet_wrap(vars(Var))

Calculate the correlation between predicted vs observed response coefficients for each covariate among all taxa within a region.

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.

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.

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.

ggplot(correlations) +
aes(GD, cor) +
geom_point() +
geom_smooth(method = "lm") +
xlim(0, 1000) +
ylim(-1, 1) +
facet_wrap("var")