R/plot_moran.R

In BlasBenito/spatialRF: Easy Spatial Modeling with Random Forest

#' @title Plots a Moran's I test of model residuals
#' @description Plots the results of spatial autocorrelation tests for a variety of functions within the package. The x axis represents the Moran's I estimate, the y axis contains the values of the distance thresholds, the dot sizes represent the p-values of the Moran's I estimate, and the red dashed line represents the theoretical null value of the Moran's I estimate.
#' @usage
#' plot_moran(
#'   model,
#'   point.color = viridis::viridis(
#'     100,
#'     option = "F",
#'     direction = -1
#'    ),
#'   line.color = "gray30",
#'   option = 1,
#'   ncol = 1,
#'   verbose = TRUE
#' )
#' @param model A model fitted with [rf()], [rf_repeat()], or [rf_spatial()], or a data frame generated by [moran()]. Default: `NULL`
#' @param point.color Colors of the plotted points. Can be a single color name (e.g. "red4"), a character vector with hexadecimal codes (e.g. "#440154FF" "#21908CFF" "#FDE725FF"), or function generating a palette (e.g. `viridis::viridis(100)`). Default: `viridis::viridis(100, option = "F")`
#' @param line.color Character string, color of the line produced by `ggplot2::geom_smooth()`. Default: `"gray30"`
#' @param option Integer, type of plot. If `1` (default) a line plot with Moran's I and p-values across distance thresholds is returned. If `2`, scatterplots of residuals versus lagged residuals per distance threshold and their corresponding slopes are returned. In models fitted with [rf_repeat()], the residuals and lags of the residuals are computed from the median residuals across repetitions. Option `2` is disabled if `x` is a data frame generated by [moran()].
#' @param ncol Number of columns of the plot. Only relevant when `option = 2`. Argument `ncol` of \link[patchwork]{wrap_plots}.
#' @param verbose Logical, if `TRUE`, the resulting plot is printed, Default: `TRUE`
#' @return A ggplot.
#' @seealso [moran()], [moran_multithreshold()]
#' @examples
#' if(interactive()){
#'
#'  #loading example data
#'  data(plant_richness_df)
#'  data(distance.matrix)
#'
#'  #fitting a random forest model
#'  rf.model <- rf(
#'    data = plant_richness_df,
#'    dependent.variable.name = "richness_species_vascular",
#'    predictor.variable.names = colnames(plant_richness_df)[5:21],
#'    distance.matrix = distance_matrix,
#'    distance.thresholds = c(0, 1000, 2000),
#'    n.cores = 1,
#'    verbose = FALSE
#'  )
#'
#'  #Incremental/multiscale Moran's I
#'  plot_moran(rf.model)
#'
#'  #Moran's scatterplot
#'  plot_moran(rf.model, option = 2)
#'
#' }
#' @rdname plot_moran
#' @export
#' @importFrom ggplot2 ggplot aes geom_hline geom_point geom_line xlab ylab ggtitle theme labs scale_colour_manual
#' @export
plot_moran <- function(
  model = NULL,
  point.color = viridis::viridis(
    100,
    option = "F",
    direction = -1
  ),
  line.color = "gray30",
  option = 1,
  ncol = 1,
  verbose = TRUE
  ){

  #checking option
  if(!(option %in% c(1, 2)) |
     inherits(model, "data.frame")){
    option <- 1
  }

  #option 1
  if(option == 1){

    #declaring variables
    distance.threshold <- NULL
    moran.i <- NULL
    moran.i.null <- NULL
    p.value.binary <- NULL
    repetition <- NULL

    if(!inherits(model, "data.frame")){
      #importance from rf_repeat
      if(inherits(model, "rf_repeat")){
        x <- model$residuals$autocorrelation$per.repetition
      } else {
        x <- model$residuals$autocorrelation$per.distance
      }
    } else {
      x <- model
    }

    #adding binary p.value
    x$p.value.binary <- "< 0.05"
    x[x$p.value >= 0.05, "p.value.binary"] <- ">= 0.05"
    x$p.value.binary <- factor(x$p.value.binary, levels = c("< 0.05", ">= 0.05"))

    #plotting rf
    if(!("repetition" %in% colnames(x))){

      p1 <- ggplot2::ggplot(data = x) +
        ggplot2::aes(
          x = distance.threshold,
          y = moran.i,
          size = p.value.binary,
          fill = moran.i
        ) +
        ggplot2::geom_line(size = 1, color = line.color) +
        ggplot2::geom_point(pch = 21) +
        ggplot2::scale_fill_gradientn(colors = point.color) +
        ggplot2::geom_hline(
          yintercept = x$moran.i.null[1],
          col = line.color,
          linetype = "dashed"
        ) +
        ggplot2::scale_size_manual(
          breaks = c("< 0.05", ">= 0.05"),
          values = c(2.5, 5),
          drop = FALSE
        ) +
        ggplot2::scale_x_continuous(breaks = x$distance.threshold) +
        ggplot2::xlab("Distance thresholds") +
        ggplot2::ylab("Moran's I of model residuals") +
        ggplot2::ggtitle("Multiscale Moran's I") +
        ggplot2::theme_bw() +
        ggplot2::theme(legend.position = "bottom") +
        ggplot2::labs(size = "Moran's I p-value") +
        ggplot2::theme(plot.title = element_text(hjust = 0.5)) +
        ggplot2::guides(fill = "none")

    } else {

      p1 <- ggplot2::ggplot(data = x) +
        ggplot2::aes(
          x = distance.threshold,
          y = moran.i,
          group = repetition,
          size = p.value.binary,
          fill = moran.i
        ) +
        ggplot2::geom_line(
          size = 1,
          color = line.color,
          alpha = 0.5
          ) +
        ggplot2::geom_point(
          pch = 21,
          alpha = 0.7
        ) +
        ggplot2::scale_fill_gradientn(colors = point.color) +
        ggplot2::geom_hline(
          yintercept = x$moran.i.null[1],
          col = line.color,
          linetype = "dashed"
        ) +
        ggplot2::scale_size_manual(
          breaks = c("< 0.05", ">= 0.05"),
          values = c(2.5, 5),
          drop = FALSE
        ) +
        ggplot2::scale_x_continuous(breaks = x$distance.threshold) +
        ggplot2::xlab("Distance thresholds") +
        ggplot2::ylab("Moran's I of model residuals") +
        ggplot2::ggtitle("Multiscale Moran's I") +
        ggplot2::theme_bw() +
        ggplot2::theme(legend.position = "bottom") +
        ggplot2::labs(size = "Moran's I p-value") +
        ggplot2::theme(plot.title = element_text(hjust = 0.5)) +
        ggplot2::guides(fill = "none")

    }

    if(verbose == TRUE){
      suppressWarnings(print(p1))
    }

    return(p1)

  }#end of option == 1

  #option 2
  if(option == 2){

    residual <- NULL
    residual.lag <- NULL

    #getting residuals (median if model is rf_repeat())
    m.residuals <- model$residuals$values

    #getting distance matrix
    if(!is.null(model$ranger.arguments$distance.matrix)){
      distance.matrix <- model$ranger.arguments$distance.matrix
    } else {
      stop("distance.matrix not found in model$ranger.arguments")
    }

    #getting distance thresholds
    if(!is.null(model$ranger.arguments$distance.thresholds)){
      distance.thresholds <- model$ranger.arguments$distance.thresholds
    } else {
      stop("distance.matrix not found in x$ranger.arguments")
    }

    #getting Moran's I results
    m.moran <- model$residuals$autocorrelation$per.distance

    #iterating through distance thresholds
    loop.out <- foreach::foreach(distance.threshold = distance.thresholds) %do% {

      #matrix of weights
      distance.weights <- weights_from_distance_matrix(
        distance.matrix = distance.matrix,
        distance.threshold = distance.threshold
      )

      #computing weighted mean of the lag
      residuals.lag <- apply(
        X = distance.weights,
        MARGIN = 1,
        FUN = function(x){x %*% m.residuals}
      )

      #subset moran i
      moran.i <- round(m.moran[m.moran$distance.threshold == distance.threshold, "moran.i"], 3)
      p.value <- round(m.moran[m.moran$distance.threshold == distance.threshold, "p.value"], 4)

      #label for facet
      facet.label <- paste0(
        "Distance = ",
        distance.threshold,
        ", I = ",
        moran.i,
        ", p = ",
        p.value
      )

      #out data frame
      out.df <- data.frame(
        distance.threshold = rep(facet.label, length(m.residuals)),
        residual = m.residuals,
        residual.lag = residuals.lag
      )

    }

    #to df
    plot.df <- do.call(
      "rbind",
      loop.out
    )

    #plot
    p2 <- ggplot2::ggplot(data = plot.df) +
      ggplot2::theme_bw() +
      ggplot2::facet_wrap(
        "distance.threshold",
        ncol = ncol) +
      ggplot2::aes(
        x = residual,
        y = residual.lag,
        color = residual
      ) +
      ggplot2::geom_point() +
      ggplot2::scale_color_gradientn(colors = rev(point.color)) +
      ggplot2::geom_hline(
        yintercept = 0,
        color = "black",
        linetype = "dotted"
      ) +
      ggplot2::geom_smooth(
        method = "lm",
        color = line.color,
        fill = "gray50",
        se = TRUE,
        alpha = 0.2,
        formula = y ~ x,
        size = 0.6
      ) +
      ggplot2::coord_fixed(expand = FALSE) +
      ggplot2::xlab("Residuals") +
      ggplot2::ylab("Lagged residuals") +
      ggplot2::theme(legend.position = "none") +
      ggplot2::ggtitle("Moran plots") +
      ggplot2::theme(plot.title = element_text(hjust = 0.5))


    if(verbose == TRUE){
      suppressWarnings(print(p2))
    }

    return(p2)

  }

}