R/ecospat.ESM.responsePlot.R

In ecospat: Spatial Ecology Miscellaneous Methods

ecospat.ESM.responsePlot <- function(ESM.EnsembleModeling.output,
                                     ESM.modeling.output,
                                     fixed.var.metric = 'median') {
  
  models <- ESM.modeling.output$models
  weights <- ESM.EnsembleModeling.output$weights
  weights.EF <- ESM.EnsembleModeling.output$weights.EF
  resp <- ESM.EnsembleModeling.output$ESM.fit$resp.var
  data <- ESM.modeling.output$data@data.env.var
  
  min.data <- apply(data, 2, min)
  max.data <- apply(data, 2, max)
  
  data.fixed.all <- as.data.frame(t(apply(data, 2, fixed.var.metric)))
  data.fixed.all <- data.fixed.all[rep(1, each = 1000), ]
  
  var.names <- colnames(data)
  
  proj.fixed.list <- list()
  
  for (i in 1:ncol(data)) {
    
    data.fixed <- data.fixed.all
    data.fixed[, i] <- seq(min.data[i], max.data[i],
                           (max.data[i] - min.data[i]) / 999)
    
    # unique projection name to avoid directory reuse
    proj.fixed <- ecospat.ESM.Projection(
      ESM.modeling.output,
      new.env = data.fixed,
      name.env = paste0("data.fixed.", i)
    )
    
    proj.fixed.list[[i]] <- ecospat.ESM.EnsembleProjection(
      ESM.prediction.output = proj.fixed,
      ESM.EnsembleModeling.output = ESM.EnsembleModeling.output
    )
    
    proj.fixed.list[[i]] <- cbind(data.fixed[, i],
                                  proj.fixed.list[[i]])
  }
  
  names(proj.fixed.list) <- colnames(data)
  
  old.par <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(old.par))
  
  xs <- floor(sqrt(length(var.names)))
  ys <- ceiling(length(var.names) / xs)
  graphics::par(mfrow = c(xs, ys))
  
  ColModels <- c("#222E50","#007991","#439A86","#BCD8C1",
                 "#E9D985","#FF6978","#6D435A","#352D39",
                 "#6E8894","#FA7921","#FE9920")
  
  for (i in 1:ncol(data)) {
    
    if (length(models) == 1) {
      
      plot(proj.fixed.list[[i]][, 2] ~ proj.fixed.list[[i]][, 1],
           xlab = "", main = names(proj.fixed.list)[i],
           ylab = "predicted value",
           ylim = c(min(sapply(proj.fixed.list, function(x) min(x[, -1]))),
                    max(sapply(proj.fixed.list, function(x) max(x[, -1])))),
           type = "n", las = TRUE)
      
      points(proj.fixed.list[[i]][, 2] ~ proj.fixed.list[[i]][, 1],
             col = "red", lwd = 2, type = "l")
      
    } else {
      
      plot(proj.fixed.list[[i]]$EF ~ proj.fixed.list[[i]][, 1],
           xlab = "", main = names(proj.fixed.list)[i],
           ylab = "predicted value",
           ylim = c(min(sapply(proj.fixed.list, function(x) min(x[, -1]))),
                    max(sapply(proj.fixed.list, function(x) max(x[, -1])))),
           type = "n", las = TRUE)
      
      legend("topleft", legend = c("ensemble", models),
             fill = c("red", ColModels[1:length(models)]),
             box.lty = 0)
      
      points(proj.fixed.list[[i]]$EF ~ proj.fixed.list[[i]][, 1],
             col = "red", lwd = 2, type = "l")
      
      for (mod.i in models) {
        points(proj.fixed.list[[i]][, mod.i] ~ proj.fixed.list[[i]][, 1],
               col = ColModels[which(models == mod.i)],
               lwd = 2, type = "l")
      }
      
      rug(data[, i], col = "black")
    }
  }
  
  # ---- cleanup ----
  unlink(grep('proj_data.fixed.ESM.BIOMOD',
              list.files(paste0(getwd(), '/', 'ESM.BIOMOD.output_', 
                                ESM.EnsembleModeling.output$species),
                         recursive = TRUE, full.names = TRUE, include.dirs = TRUE),
              value = TRUE),
         recursive = TRUE)
  
  unlink(grep("ESM_Projections.data.fixed",
              list.files(paste0(getwd(), "/", 
                                "ESM.BIOMOD.output_", 
                                ESM.EnsembleModeling.output$species),
                         full.names = TRUE),
              value = TRUE))
  
  return(proj.fixed.list = proj.fixed.list)
}