Plots for model of Grampians Occupancy and Trait Data with quadratic relationships

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

Coefficient Plot

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

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

coef_grampians$Var1 <- fct_recode(
  coef_grampians$Var1,
  `Eucalyptus goniocalyx` = "Eucalyptus (goniocalyx|nortonii)",
  `Eucalyptus viminalis subsp. viminalis` =
    "Eucalyptus viminalis subsp. (pryoriana|viminalis)"
)

coef_grampians$Var1 <- as.character(coef_grampians$Var1)
coef_grampians$Var1 <- gsub("Eucalyptus", "E.", coef_grampians$Var1)
coef_grampians$Var1 <- gsub("subsp.", "ssp.", coef_grampians$Var1)

ggplot(coef_grampians) +
aes(factor(Var1, levels = rev(levels(factor(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", covariate_nms))
) +
scale_x_discrete(name = "") +
ylab("Standardised coefficients") +
theme_bw() +
theme(
  axis.text.y = element_text(face = "italic"),
  strip.text.x = element_text(size = 6)
)

Trait-Environmental Response Relationships

trait_env_plot(
  model_grampians,
  mdg,
  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 = covariate_nms,
  plot.log = c("", "", "x")
)

Traits-Environmental Response 2D Surfaces

unscaled_covariates <- mdg[
  c(
    "mean_moisture_index_of_low_qtr", "topographic_wetness_index_3sec",
    "relief_1000m_radius", "total_nitrogen"
  )
]

covariate_polys <- new_poly(unscaled_covariates, 1:2, length.out = 50)
covariate_polys_expd <- covariate_polys[rep(seq(50), times = 50), ]
covariate_polys_expd[, 3:4] <- covariate_polys[rep(seq(50), each = 50), 3:4]
colnames(covariate_polys_expd) <- colnames(covariate_polys) <-
  names(covariate_nms)

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)
    data.frame(
      x = covariate_polys_expd[, "mlq"],
      y = covariate_polys_expd[, "twi"],
      x_raw = rep(seq_rng(mdg$mean_moisture_index_of_low_qtr, na.rm = TRUE, length.out = 50), times = 50),
      y_raw = rep(seq_rng(mdg$topographic_wetness_index_3sec, na.rm = TRUE, length.out = 50), each = 50),
      z = predict_prob_new_taxon(
        model_grampians,
        data.frame(
          covariate_polys_expd,
          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)[6:7] <- c("SLA", "Seed Mass")

ggplot(panel_data) +
aes(x_raw, y_raw, 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 = covariate_polys[, "mlq"],
    y = covariate_polys[, "twi"]
  ),
  z = matrix(
    predict_prob_taxon(
      model_grampians,
      "Eucalyptus radiata",
      covariate_polys_expd
    ),
    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 = covariate_polys_expd[, "mlq"],
    y = covariate_polys_expd[, "twi"],
    z = predict_prob_taxon(model_grampians, i, covariate_polys_expd)
  )

  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" =
    new_poly(unscaled_covariates, 1, colnames = names(covariate_nms)
  ),
  "Topographic Wetness Index" =
    new_poly(unscaled_covariates, 2, colnames = names(covariate_nms)
  ),
  "Relief 1000m Radius" =
    new_poly(unscaled_covariates, 3, colnames = names(covariate_nms)
  ),
  "Total Nitrogen" =
    new_poly(unscaled_covariates, 4, colnames = names(covariate_nms)
  )
)

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

par(mfrow = c(1, 4), oma = c(0, 4, 1, 0), las = 1)
for (j in seq_along(np)) {
  par(mar = c(5, 0, 0, 1))
  plot(
    x[[j]],
    predict_prob_taxon(model_grampians, "Eucalyptus melliodora", np[[j]]),
    ylim = 0:1, type = "l", lwd = 3, col = "black",
    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 == "Eucalyptus melliodora"),
    points(
      get(names(x[j])), ifelse(occupancy, 1.025, -.01), pch = 15,
      col = rgb(.2, .2, .2, .1)
    )
  )
  points(
    x[[j]],
    predict_prob_taxon(model_grampians, "Eucalyptus obliqua", np[[j]]),
    type = "l", lwd = 3, col = "grey"
  )
  with(
    subset(mdg, taxon == "Eucalyptus obliqua"),
    points(
      get(names(x[j])), ifelse(occupancy, 1.01, -.025), pch = 15,
      col =  rgb(.8, .8, .8, .1)
    )
  )
}
legend(
  "topright",
  legend = c("Eucalyptus melliodora", "Eucalyptus obliqua"),
  fill = c("black", "grey"), ncol = 1, bty = "n", text.font = 3,
  border = NA, y.intersp = 2, inset = .05
)