Validation of quadratic model

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

The following is a validation of the predictions/inference from a trait-environment model of the Grampians Bioregion ~20 taxa dataset.

The following steps are undertaken:

• Fit a model to Grampians data.
• Calculate Performance Metrics:
  • Predict probability of occurrence
    • Predict within Grampians training set:
      • Predict for each species including species random effects.
      • Predict for each species without species random effects.
    • Predict to the South-east validation set:
      • Predict to each IBRA subregion within the South-east:
        • Predict for each species common to Grampians and South-east subregion:
          • Predict including species random effects.
          • Predict without species random effects.
        • Predict for each species not found in the Grampians training set.
  • For each prediction type (within-sample/out-of-sample, with-re/without-re, and modelled-species/unmodelled-species), bioregion, and species (and their combinations) calculate performance metrics

mds <- local({
  nms1 <-   c(
    "mean_moisture_index_of_low_qtr", "topographic_wetness_index_3sec",
    "relief_1000m_radius", "total_nitrogen"
  )
  nms2 <- c("sla_mm2_per_mg", "seed_mass_mg", "max_height_m")
  vars <- lapply(nms1, function(x) poly(log(grampians[[x]]), 2))
  coefs <- lapply(vars, function(x) attr(x, "coefs"))
  vars <- do.call(cbind, vars)
  vars <- cbind(
    vars, do.call(cbind, lapply(nms2, function(x) log(grampians[[x]])))
  )
  vars <- scale(vars)
  cntrs <- attr(vars, "scaled:center")
  scales <- attr(vars, "scaled:scale")

  vars_new_raw <- na.omit(southeast)
  vars_new <- lapply(
    seq_along(nms1), function(x) {
      poly(log(vars_new_raw[[nms1[x]]]), degree = 2, coefs = coefs[[x]])
    }
  )
  vars_new <- do.call(cbind, vars_new)
  vars_new <- cbind(
    vars_new, do.call(cbind, lapply(nms2, function(x) log(vars_new_raw[[x]])))
  )
  vars_new <- scale(vars_new, center = cntrs, scale = scales)


  setNames(
    data.frame(
      vars_new_raw[, c("occupancy", "taxon", "ibra_subregion")], vars_new
    ),
    c(
      "y", "taxon", "ibra_subregion", "mlq", "mlq2","twi", "twi2", "r1k",
      "r1k2", "tn", "tn2", "sla", "sm", "mh"
    )
  )
})

perf <- local({
  perf <- mapply(
    function(f, re.form) {
      f(
        model_grampians, y, list(taxon, ibra_subregion),
        rbind(mds, mdg[names(mds)]), re.form
      )
    },
    f = rep(list(auroc, pred_inf, dev_explained), each = 2),
    re.form = list(NULL,  NA)
  )
  taxon <- sub("(.*)\\..*", "\\1", rownames(perf))
  ibra_subregion <- sub(".*\\.", "", rownames(perf))
  dim(perf) <- dim(perf) * c(2, .5)
  colnames(perf) <- c("AUC", "mean_predictive_info", "deviance_explained")
  cbind(
    taxon = rep(taxon, 2),
    ibra_subregion = rep(ibra_subregion, 2),
    random_effect = rep(c(TRUE, FALSE), each = nrow(perf) / 2),
    as.data.frame(perf)
  )
})

aucdist <- data.frame(
  perf[c("taxon", "ibra_subregion", "AUC")],
  `Geographic (km)`= dists[as.character(perf$ibra_subregion)],
  Environmental = kldists[as.character(perf$ibra_subregion)],
  Community = community_dist[as.character(perf$ibra_subregion)],
  "Dummy" = "Dummy",
  check.names = FALSE
)

aucdist <- gather(
  aucdist, "distance", "value", `Geographic (km)`, Environmental, Community
)

aucdist <- na.omit(aucdist)

aucdistmean <- summarise(
  group_by(aucdist, ibra_subregion, distance),
  meanauc = mean(AUC), value = mean(value)
)

aucdistplot <- ggplot(aucdist) +
  aes(AUC, value) +
  geom_vline(xintercept = .5, color = "grey", size = 1.5) +
  geom_point(alpha=.1) +
  geom_point(
    aes(meanauc), data = aucdistmean, size = 2, shape = 21, fill = "white"
  ) +
  xlim(0, 1) +
  ylab("Distance") +
  facet_grid(distance ~ ., scales = "free_y") +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    strip.background = element_blank()
  )

aucbw <- ggplot(aucdist) +
  aes(x = Dummy, y = AUC) +
  geom_boxplot(fill = "grey") +
  coord_flip() +
  facet_grid(Dummy ~ .) +
  ylim(0, 1) +
  scale_x_discrete(labels = "0.95") +
  xlab("Dummy") +
  theme_bw() +
  theme(
    axis.title.y = element_text(color = "transparent"),
    axis.text.y  = element_text(color = "transparent"),
    axis.ticks.y = element_line(color = "transparent"),
    strip.background = element_blank(),
    strip.text.y = element_text(color = "transparent"),
    plot.margin = margin(0, 5.5, 5.5, 5.5)
  )

grid.arrange(aucdistplot, aucbw, ncol = 1, heights = c(.86, .14))

ggplot(
  transform(
    subset(perf, ibra_subregion != "Greater Grampians" & !random_effect),
    ibra_subregion = factor(
      ibra_subregion, names(sort(tapply(AUC, ibra_subregion, median)))
    )
  )
) +
aes(AUC) + geom_histogram() + facet_wrap(~ibra_subregion, 3) +
xlab("AUC") +
ylab("No. of taxa") +
ggtitle(
  paste(
    "Area Under Receiver Operator Curve for taxa within IBRA subregions of",
    "Southeast Australia"
  )
)

ggplot(
  transform(
    subset(perf, ibra_subregion != "Greater Grampians" & !random_effect),
    ibra_subregion = factor(
      ibra_subregion,
      names(sort(tapply(mean_predictive_info, ibra_subregion, median)))
    )
  )
) +
aes(mean_predictive_info) + geom_histogram() + facet_wrap(~ibra_subregion, 3) +
xlab("Mean Predictive Information") +
ylab("No. of taxa") +
ggtitle(
  paste(
    "Mean Predictive Information for taxa within IBRA subregions of",
    "Southeast Australia"
  )
)

ggplot(
  transform(
    subset(perf, taxon %in% mdg$taxon & random_effect),
    taxon = factor(
      taxon,
      names(
        sort(
          tapply(
            AUC[ibra_subregion == "Greater Grampians"],
            taxon[ibra_subregion == "Greater Grampians"],
            median
          )
        )
      )
    )
  )
) +
aes(
  AUC,
  mean_predictive_info,
  col = ifelse(
    ibra_subregion == "Greater Grampians",
    "Grampians",
    "Other Region"
  )
) +
geom_point() +
facet_wrap(~taxon, labeller = label_wrap_gen(40)) +
labs(col = NULL) +
xlab("AUC") +
ylab("Mean Predictive Information") +
ggtitle(
  paste(
    ""
  )
) +
theme(legend.position = "bottom")

model_auc <- stan_glmer(
  (AUC - .001) ~ value + (1 | taxon) + (1 | ibra_subregion),
  family = mgcv::betar,
  iter = 1000, chains = 3,
  data = aucdist,
  subset = ibra_subregion != "Greater Grampians" & distance == "Geographic (km)"
)

summary(model_auc, regex_pars = "^[^b]", probs = c(.025, .975))

Model Info:
 function:     stan_glmer
 family:       beta [logit, link.phi=log]
 formula:      (AUC - 0.001) ~ value + (1 | taxon) + (1 | ibra_subregion)
 algorithm:    sampling
 sample:       1500 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 1120
 groups:       taxon (82), ibra_subregion (18)

Estimates:
                                                mean   sd   2.5%   97.5%
(Intercept)                                    0.9    0.2  0.5    1.4   
value                                          0.0    0.0  0.0    0.0   
(phi)                                          9.9    0.4  9.1   10.8   
Sigma[taxon:(Intercept),(Intercept)]           0.2    0.0  0.1    0.2   
Sigma[ibra_subregion:(Intercept),(Intercept)]  0.1    0.0  0.0    0.1   

Fit Diagnostics:
           mean   sd   2.5%   97.5%
mean_PPD 0.7    0.0  0.7    0.7    

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
                                              mcse Rhat n_eff
(Intercept)                                   0.0  1.0   751 
value                                         0.0  1.0   723 
(phi)                                         0.0  1.0  2248 
Sigma[taxon:(Intercept),(Intercept)]          0.0  1.0   429 
Sigma[ibra_subregion:(Intercept),(Intercept)] 0.0  1.0   723 
mean_PPD                                      0.0  1.0  1485 
log-posterior                                 0.5  1.0   375 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).