13 Pseudo-RCS & Distance to fireworks

confidence.intervals <- c(0.95, 0.8, 0.5)
colors <- c("urban" = "#94346E", 
            "agricultural" = "#73AF48", 
            "semiopen" = "#EDAD08", 
            "forests" = "#0F8554", 
            "wetlands" = "#38A6A5", 
            "waterbodies" = "#1D6996",
            "dist_urban" = "#94346E", 
            "total_rcs" = "#CC503E", 
            "dist_radar" = "#666666")
x_labels <- c("urban" = "Urban", 
              "agricultural" = "Agricultural", 
              "semiopen" = "Semi-open", 
              "forests" = "Forests", 
              "wetlands" = "Wetlands", 
              "waterbodies" = "Water bodies", 
              "dist_urban" = "Distance to fireworks (m)", 
              "total_rcs" = "Total RCS (g^(2/3))", 
              "dist_radar" = "Distance from radar (m)")
limit_quantiles <- c("dist_urban", "total_rcs")
manual_limits <- list("dist_urban" = c(min(data_cleaned$dist_urban), quantile(data_cleaned$dist_urban, probs = 0.975)),
                      "total_rcs" = c(0, quantile(data_cleaned$total_rcs, probs = 0.975)))
# ylims <- c(-1.8, 1.5) + modelci$model$offset
ylims <- c(-1, 1)

plot_mboost_pdp <- function(modelci, data, which, confints, colors = NULL, ylims = ylims, varimp = NULL, offset = FALSE) {
  # Extract bootstrapped quantiles
  quants <- t(mboost_bootstrapped_quantiles(modelci, confidence.intervals, which = which))
  if (offset) quants <- quants + modelci$model$offset
  bootstrapped_quantiles <- as.data.frame(quants)
  # bootstrapped_quantiles <- as.data.frame(t(mboost_bootstrapped_quantiles(modelci, confidence.intervals, which = which)))
  bootstrapped_quantiles$x <- modelci$data[modelci$model$which(which)][[1]][, 1]
  bootstrapped_quantiles$y <- plot.mboost_adjusted(modelci$model, which = which, newdata = modelci$data[[modelci$model$which(which)]])[[2]]
  if (offset) {
    bootstrapped_quantiles$y <- bootstrapped_quantiles$y + modelci$model$offset
    ylims <- ylims + modelci$model$offset
  }
  
  p <- ggplot(bootstrapped_quantiles)
  i <- 1
  sorted_confints <- sort(confints, decreasing = TRUE)
  
  # Add variable importance
  if (is.null(varimp)) {
    variable_importance <- function(model, which, exclude = c("dist_radar"), round = 0) {
      vi <- as.data.frame(varimp(model)) %>%
        filter(!variable %in% exclude)
      vi$vi <- vi$reduction / sum(vi$reduction) * 100
      vi %>%
        filter(variable == which) %>%
        dplyr::select(vi) %>%
        as.numeric() %>%
        round(round)
    }
    vi <- variable_importance(modelci$model, which)
  } else {
    vi <- varimp %>%
      filter(variable == which) %>%
      dplyr::select(mean_round)
  }
  
  offset_right <- 0.95
  offset_top <- 0.35
  if (which == "dist_urban") {
    p <- p +
      annotate("text", x = manual_limits$dist_urban[[2]] * offset_right, y = ylims[2] - offset_top, label = paste0(vi, "%"), hjust = 1, color = colors[which],
               fontface = "bold", size = 12, alpha = 0.5)
  } else if(which == "total_rcs") {
    p <- p +
      annotate("text", x = manual_limits$total_rcs[[2]] * offset_right, y = ylims[2] - offset_top, label = paste0(vi, "%"), hjust = 1, color = colors[which],
               fontface = "bold", size = 12, alpha = 0.5)
  } else {
    p <- p +
      annotate("text", x = 1 * offset_right, y = ylims[2] - offset_top, label = paste0(x_labels[which], ": ", vi, "%"), hjust = 1, color = colors[which], 
               fontface = "bold", size = 12, alpha = 0.5)
  }
  
  # Add density
  if (which == "dist_urban") {
    dens <- density(data[, which], bw = 225, from = manual_limits$dist_urban[[1]], to = manual_limits$dist_urban[[2]])
  } else if (which == "total_rcs") {
    dens <- density(data[, which], bw = 2000, from = manual_limits$total_rcs[[1]], to = manual_limits$total_rcs[[2]])
  } else {
    dens <- density(data[, which], bw = 0.03556, from = 0, to = 1)
  }
  scale_lims <- c(ylims[1] + 0.05, ylims[2] - 0.05)
  dens$y <- scales::rescale(dens$y, to = scale_lims)
  dens <- as.data.frame(list(as.matrix(dens$x), as.matrix(dens$y)))
  colnames(dens) <- c("x", "y")
  p <- p +
    geom_line(aes(x = x, y = y), data =  dens, color = "grey50", lineend = "round", linetype = 3) +
    scale_y_continuous(sec.axis = sec_axis(~ ., name = "Predictor frequency", breaks = scale_lims, labels = c("min", "max")))
  
  ## Add horizontal line
  if (offset) {
    yintercept <- modelci$model$offset
  } else {
    yintercept <- 0
  }
  p <- p +
    geom_hline(yintercept = yintercept, col = "darkgrey")
  
  for (ci in sorted_confints) {
    alphas <- rev(factor(sorted_confints))
    ymax <- as.name(paste((1 - (1 - ci) / 2) * 100, "%", sep = ""))
    ymax <- enquo(ymax)
    ymin <- as.name(paste(((1 - ci) / 2) * 100, "%", sep = ""))
    ymin <- enquo(ymin)
    
    p <- p +
      geom_ribbon(aes(x = x, ymin = !!ymin, ymax = !!ymax, alpha = !!alphas[i]), fill = colors[which]) +
      geom_line(aes(x = x, y = y), color = colors[which], size = 1.5, lineend = "round")
    
    if (offset) {
      p <- p + scale_y_continuous(labels = function(x) {10^x})
    }

    i <- i + 1
  }
  
  p <- p +
    scale_alpha_discrete(range = c(0.2, 0.5)) +
    coord_cartesian(ylim = ylims, expand = FALSE) +
    guides(alpha = FALSE)
  
  if (which == "dist_urban") {
    xlabel <- x_labels[which]
    p <- p +
      coord_cartesian(xlim = manual_limits$dist_urban, ylim = ylims, expand = FALSE)
      # scale_x_continuous(breaks = c(0, 10000, 20000, 30000), labels = c(0, 10, 20, 30))
  } else if (which == "total_rcs") {
    xlabel <- x_labels[which]
    p <- p +
      coord_cartesian(xlim = manual_limits$total_rcs, ylim = ylims, expand = FALSE)
      # scale_x_continuous(breaks = c(250, 500, 750), labels = c(250*4, 500*4, 750*4))
  } else {
    xlabel <- paste("% Coverage")
    p <- p + 
      scale_x_continuous(breaks = c(0, 0.25, 0.5, 0.75, 1), labels = c("0%", "25%", "50%", "75%", "100%"))
  }
  
  p <- p +
    theme_classic(base_size = 10) +
    labs(x = xlabel, y = expression(paste("Flight response ", log[10], "(VIR)"))) +
    theme(axis.line.y.left = element_line(color = colors[which], size = 1),
          axis.line.y.right = element_line(color = "grey"),
          axis.ticks.y.left = element_line(color = colors[which], size = 1),
          axis.ticks.y.right = element_line(color = "grey"),
          axis.title.y = element_text(),
          axis.title.y.right = element_blank())
  
  if (!which == "forests") {
    p <- p +
      theme(axis.title.y.left = element_blank())
  }
  
  if (!which == "waterbodies") {
    p <- p +
      theme(axis.title.x = element_blank())
  }

  p
}

predictors <- c("dist_urban", "total_rcs")
plots <- mapply(function(predictors) {plot_mboost_pdp(modelci, data_cleaned, which = predictors, confidence.intervals, colors, ylims, bvi_biol, 
                                                      offset = FALSE)}, 
                predictors = predictors,
                SIMPLIFY = FALSE)

saveRDS(plots, file = "data/plots/paper/fig_rcs_disturban.RDS")

plots[[1]] + plots[[2]] +
plot_layout(widths = c(1, 1), nrow = 1) & 
theme(axis.text.x = element_text(size = 10),
      axis.title.x = element_text(size = 12),
      axis.text.y = element_text(size = 10)) -> plots_rcs_disturban

ggsave(filename = "data/plots/paper/fig_rcs_disturban.pdf", plot = plots_rcs_disturban, width = 14, height = 5, dpi = 300, units = "cm")

plots_rcs_disturban