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

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

covariate_nms <- c(
  mlq  = "Mean Moisture\nIndex Lowest\nQuarter",
  mlq2 = "Mean Moisture\nIndex Lowest\nQuarter\u00b2",
  twi  = "Topographic\nWetness\nIndex",
  twi2 = "Topographic\nWetness\nIndex\u00b2",
  r1k  = "Relief 1000m\nRadius",
  r1k2 = "Relief 1000m\nRadius\u00b2",
  tn   = "Total Nitrogen",
  tn2  = "Total Nitrogen\u00b2"
)

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, 5, 1, 1), lend = 1)
barplot(
  full_sd[-1], las = 1, col = "darkgrey", border = NA, cex.names = .8,
  mgp = c(3, 2, 0), ylim = c(0, 1.5), 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 = 30, inset = .05
)
box()
axis(2, las = 1)
mtext("Standard Deviations", 2, 3)