Plots for model of Grampians Occupancy and Trait Data with linear relationships and location random effect

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

Coefficient Plot

coef_grampians <- as.data.frame.table(coef_sims(model_grampians))

ggplot(coef_grampians) +
aes(Var1, Freq) +
geom_hline(yintercept = 0) +
geom_violin(fill = "grey35", scale = "width", width = .5, trim = TRUE) +
coord_flip() +
facet_grid(
  ~Var2,
  labeller = as_labeller(
    c(
      `(Intercept)` = "Intercept",
      mlq  = "Mean Moisture\nIndex Lowest\nQuarter",
      twi  = "Topographic\nWetness\nIndex",
      r1k  = "Relief 1000m\nRadius",
      tn  = "Total Nitrogen"
    )
  )
) +
scale_x_discrete(name = "", limits = rev(levels(coef_grampians$Var1))) +
ylab("Standardised coefficients") +
theme_bw()

Trait-Environmental Response Relationships

trait_env_plot(
  model_grampians,
  modeldata_grampians,
  gvars = c(
    sla = "sla_mm2_per_mg", sm = "seed_mass_mg", mh = "max_height_m"
  ),
  xlabs = list(
    expression(SLA~(mm^2~mg^-1)), "Seed Mass (mg)", "Maximum Height (m)"
  ),
  ylabs = c(
    "Mean Moisture Index\nLowest Quarter", "Topographic\nWetness Index",
    "Relief 1000m Radius", "Total Nitrogen"
  )
)

Traits-Environmental Response 2D Surfaces

panel_data <- mapply(
  function(sla_raw, sm_raw) {
    sla    <- log(mdg$sla_mm2_per_mg)
    sla_m  <- mean(sla)
    sla_sd <- sd(sla)
    sm     <- log(mdg$seed_mass_mg)
    sm_m   <- mean(sm)
    sm_sd  <- sd(sm)
    x <- rep(seq_rng(mdg$mlq, na.rm = TRUE, length.out = 50), times = 50)
    y <- rep(seq_rng(mdg$twi, na.rm = TRUE, length.out = 50), each = 50)
    data.frame(
      x = x,
      y = y,
      z = predict_prob_new_taxon(
        model_grampians,
        data.frame(
          mlq = x,
          twi = y,
          r1k = 0,
          tn = 0,
          sla = (log(sla_raw) - sla_m) / sla_sd,
          sm  = (log(sm_raw) - sm_m) / sm_sd,
          mh = 0
        )
      ),
      sla = sla_raw,
      sm  = sm_raw
    )
  },
  sla_raw = c(3, 5, 3, 5),
  sm_raw = c(.5, .5, 1.5, 1.5),
  SIMPLIFY = FALSE
)

panel_data <- do.call(rbind, panel_data)
colnames(panel_data)[4:5] <- c("SLA", "Seed Mass")

ggplot(panel_data) +
aes(x, y, z = z, fill = z) +
geom_raster(alpha = .7) +
geom_contour(col = "grey") +
scale_fill_gradientn(
  name = "Pr(Y=1)", colours = grey.colors(50, 1, .1, gamma = .3)
) +
xlab("Mean Moisture Index Lowest Quarter") +
ylab("Topographic Wetness Index") +
facet_grid(`Seed Mass` ~ SLA, labeller = label_both) +
theme_bw()

3D Response Surface

persp(
  list(
    x = seq_rng(mdg$mlq, na.rm = TRUE, length.out = 50),
    y = seq_rng(mdg$twi, na.rm = TRUE, length.out = 50)
  ),
  z = matrix(
    predict_prob_taxon(
      model_grampians,
      "Eucalyptus radiata",
      data.frame(
        mlq = rep(seq_rng(mdg$mlq, na.rm = TRUE, length.out = 50), times = 50),
        twi = rep(seq_rng(mdg$twi, na.rm = TRUE, length.out = 50), each = 50),
        r1k = 0,
        tn = 0
      )
    ),
    50
  ),
  main = "Eucalyptus radiata",
  zlab = "\nPr(Y=1)",
  xlab = "\nMean Moisture Index Lowest Quarter",
  ylab = "\nTopographic Wetness Index",
  col = "cornflowerblue", border = "salmon",
  theta = 25, phi = 45,
  font.main = 3
)

2D Surface Predictions

for (i in levels(model_grampians@flist$taxon)) {
  panel_data <- data.frame(
    x = rep(seq_rng(mdg$mlq, na.rm = TRUE, length.out = 50), times = 50),
    y = rep(seq_rng(mdg$twi, na.rm = TRUE, length.out = 50), each = 50),
    z = predict_prob_taxon(
      model_grampians, i,
      cbind(
        mlq = rep(seq_rng(mdg$mlq, na.rm = TRUE, length.out = 50), times = 50),
        twi = rep(seq_rng(mdg$twi, na.rm = TRUE, length.out = 50), each = 50),
        r1k = 0,
        tn = 0
      )
    )
  )

  rug_data <- lapply(
    c(absent = FALSE, present = TRUE),
    function(x) {
      ans <- subset(
        mdg,
        taxon == i & occupancy == x,
        c("mlq", "twi")
      )
      names(ans) <- c("x", "y")
      ans
    }
  )

  densityplot <-
    ggplot(panel_data) +
    aes(x, y, z = z, fill = z) +
    geom_raster(alpha = .7) +
    geom_rug(
      data = rug_data$absent, aes(x, y), inherit.aes = FALSE, col = "grey"
    ) +
    geom_rug(
      data = rug_data$present, aes(x, y), inherit.aes = FALSE, col = "black"
    ) +
    scale_fill_gradientn(name = "Pr(Y=1)", colours = grey.colors(50, 1, 0)) +
    xlab("Mean Moisture Index Lowest Quarter") +
    ylab("Topographic Wetness Index") +
    ggtitle(i) +
    theme_minimal() +
    theme(plot.title = element_text(face = "italic"))

  print(densityplot)
}

Partial Dependence Plots

x <- list(
  mean_moisture_index_of_low_qtr =
    seq_rng(mdg$mean_moisture_index_of_low_qtr, na.rm = TRUE),
  topographic_wetness_index_3sec =
    seq_rng(mdg$topographic_wetness_index_3sec, na.rm = TRUE),
  relief_1000m_radius =
    seq_rng(mdg$relief_1000m_radius, na.rm = TRUE),
  total_nitrogen =
    seq_rng(mdg$total_nitrogen, na.rm = TRUE)
)
np <- list(
  "Mean Moisture Index Lowest Quarter" = data.frame(
    mlq = seq_rng(mdg$mlq, na.rm = TRUE), twi = 0, r1k = 0, tn = 0
  ),
  "Topographic Wetness Index" = data.frame(
    twi = seq_rng(mdg$mlq, na.rm = TRUE), mlq = 0, r1k = 0, tn = 0
  ),
  "Relief 1000m Radius" =data.frame(
    r1k = seq_rng(mdg$mlq, na.rm = TRUE), twi = 0, mlq = 0, tn = 0
  ),
  "Total Nitrogen" = data.frame(
    tn = seq_rng(mdg$mlq, na.rm = TRUE), twi = 0, r1k = 0, mlq = 0
  )
)

par(mfrow = c(1, 4), oma = c(0, 4, 4, 0), font.main = 3, las = 1)
for(i in levels(model_grampians@flist$taxon)) {
  for (j in seq_along(np)) {
    par(mar = c(5, 0, 0, 1))
    plot(x[[j]], predict_prob_taxon(model_grampians, i, np[[j]]),
      ylim = 0:1, type = "l", lwd = 3,
      ylab = ifelse(j == 1, "Pr(Y=1)", ""),
      xlab = names(np)[j], yaxt = ifelse(j == 1, "s", "n"),
      panel.first = grid(), xpd = NA
    )
    with(
      subset(mdg, taxon == i),
      points(
        get(names(x[j])), ifelse(occupancy, 1.02, -.02), pch = 15,
        col = rgb(.1, .1, .1, .1)
      )
    )
  }
  title(i, outer = TRUE)
}