R/user_add_neighbors.R

In RSTr: Gibbs Samplers for Discrete Bayesian Spatiotemporal Models

#' Add neighbors to adjacency information
#'
#' Modifies \code{adjacency} to indicate that \code{neighs} they should be treated as neighbors.
#'
#' \code{add_neighbors()} is useful when adjacency information generated by \code{spdep::poly2nb()} indicates lone regions without links/neighbors, particularly in island counties such as the Hawaiian islands, Nantucket in Massachusetts, or San Juan in Washington. Note that \code{add_neighbors()} marks all listed counties as adjacent, so if you have a set of chaining counties where the first may not be connected to the last, several instances of \code{add_neighbors()} will be needed.
#'
#' @param adjacency Adjacency information generated by \code{spdep::poly2nb()}.
#' @param neighs A vector of regions to mark as adjacent. Accepts a vector of indices or names assigned to \code{adjacency}.
#' @returns A modified \code{adjacency} \code{list}.
#' @examples
#' if (requireNamespace("sf", quietly = TRUE) &&
#'     requireNamespace("spdep", quietly = TRUE)) {
#'
#'   mamap <- sf::st_as_sf(mamap[order(mamap$GEOID), ])
#'   ma_adj <- spdep::poly2nb(mamap)
#'   new_neighs <- c(1, 4, 10) # attach regions 1, 4, and 10
#'   ma_adj <- add_neighbors(ma_adj, new_neighs)
#'
#'   # Add neighbors by FIPS code instead of index
#'   ma_adj <- suppressWarnings(spdep::poly2nb(mamap))
#'   names(ma_adj) <- mamap$GEOID
#'   ma_adj <- add_neighbors(ma_adj, neighs = c("25001", "25007", "25019"))
#'
#'   ma_adj <- suppressWarnings(spdep::poly2nb(mamap))
#'   ma_adj <- add_neighbors(ma_adj, c(1, 4)) # only attach 1 and 4
#'   ma_adj <- add_neighbors(ma_adj, c(4, 10)) # only attach 4 and 10
#' }
#' @export
add_neighbors <- function(adjacency, neighs) {
  if (is.character(neighs)) {
    if (!all(neighs %in% names(adjacency))) {
      stop(
        "Not all `neighs` match names in `adjacency`. Check `neighs %in% names(adjacency)` for FALSE values."
      )
    }
    neighs <- match(neighs, names(adjacency))
  }
  if (any(neighs < 1 | neighs > length(adjacency))) {
    stop(
      "Not all indices in `neigh` are within `adjacency`. check `neighs %in% seq_along(adjacency)` for FALSE values."
    )
  }
  adjacency <- as_nb(adjacency)
  adjacency_card <- spdep::card(adjacency)
  for (neighbor in neighs) {
    if (adjacency_card[neighbor] == 0) {
      adjacency[[neighbor]] <- integer(length = 0)
    }
    adjacency[[neighbor]] <- union(
      adjacency[[neighbor]],
      setdiff(neighs, neighbor)
    )
  }
  as_nb(adjacency)
}

as_nb <- function(adjacency) {
  adjacency <- lapply(adjacency, as.integer)
  class(adjacency) <- "nb"
  adjacency
}