Explanatory power of traits in the trait-environment model with linear responses

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

covariate_nms <- c(
  mlq  = "Mean\nMoisture\nIndex\nLowest\nQuarter",
  twi  = "Topographic\nWetness\nIndex",
  r1k  = "Relief\n1000m\nRadius",
  tn   = "Total\nNitrogen"
)

residual_variance <- VarCorr(model_grampians)$taxon

trait_values <- unique(mdg[c("taxon", "sla", "sm", "mh")])

full_sd <- sqrt(diag(cov(apply(coef_sims(model_grampians), 1:2, mean))))

resid_sd <- attr(residual_variance, "stddev")

par(mar = c(4.1, 5, 1, 1), lend = 1)
barplot(
  full_sd[-1], las = 1, col = "darkgrey", border = NA,
  mgp = c(3, 3.1, 0), ylim = c(0, 1), yaxt = "n", ylab = "",
  names.arg = covariate_nms,
  panel.last = for(i in seq(.2, .8, .2)) {
    abline(h=i, lty = 3, xpd = FALSE, col = "grey")
  }
)
barplot(
  pmax(0, full_sd[-1] - resid_sd[-1]), col = "lightgrey", add = TRUE,
  axes = FALSE, names.arg = NA, border = NA
)
legend(
  "topright",
  legend = c("Variation\nexplained\nby traits", "Variation\nUnexplained"),
  col = c("lightgrey", "darkgrey"), ncol = 2, bty = "n",
  border = NA, lty = 1, lwd = 20
)
box()
axis(2, las = 1)
mtext("Standard Deviations", 2, 3)