r-spatial / spdep

Spatial Dependence: Weighting Schemes and Statistics
https://r-spatial.github.io/spdep/
116 stars 26 forks source link

Set `sym = TRUE` on dnearneigh when lower = 0 #159

Closed JosiahParry closed 2 weeks ago

JosiahParry commented 2 weeks ago

When using a fixed distance band to identify neighbors we can be sure that the neighbors matrix is symmetrix when the lower bound is set to zero. We should set the attribute sym = TRUE in these cases.

rsbivand commented 2 weeks ago

Thanks!

Now: https://github.com/r-spatial/spdep/blob/0e1001961f1e2d7911af949bc8441345cdcd31dc/R/dnearneigh.R#L125-L126 by default sets attr(z, "sym") <- TRUE unless symtest is TRUE, in which case symmetry is checked. My expectation would be that d_ij == d_ji, but that the bounds argument might be set such that this did not hold.

Do you have a use case where the d_ij == d_ji equality breaks down where d1 != 0? Or in other cases? I didn't check after s2 arrived, but I'd expect symmetry there too.

JosiahParry commented 2 weeks ago

Dang both you and @edzer are fast! You're correct. This is a misinterpretation of a message from spacemakeR.

library(spdep)

nb <- dnearneigh(spData::cycle_hire, 0, 0.75)

# shows the symmetric attribute
attr(nb, "sym")
#> [1] TRUE

# we get TRUE
is.symmetric.nb(nb)
#> [1] TRUE

# create a rowstandardized spatial weights matrix 
listw <- nb2listw(nb, style = "W")

# calculate MEMs
# we get asymmetric weights matrix
res <- spacemakeR::scores.listw(listw)
#> listw not symmetric, (w+t(w)) used in the place of w

This confused me because since the neighbor matrix is symmetric, I assumed the SWM would be symmetric as well. But that is only true in the case of the indicator / binary weights matrix. The corresponding weights matrix will almost never be symmetric.

Time to take a break and eat breakfast. 😅

rsbivand commented 2 weeks ago

See also: https://github.com/r-spatial/spatialreg/blob/bfb1ebe1b5b15f2eb20d2fb24dd9bc93f69dbe84/R/jacobian_setup.R#L270-L337 and friends using Ord's (1975, p. 125, appendix C) result that the eigenvalues of DC are same as those of D^{1/2}CD^{1/2} if D is diagonal with the row sums of C on the principal diagonal: dii=1(sum{j=1,n} C_{ij}), and C is symmetric (doesn't need to be binary, but needs to be symmetric).

Did Stéphane release the package?

JosiahParry commented 2 weeks ago

Thank you! I think this math is a bit beyond me, unfortunately!

I believe it was written by Stéphane Dray based on the authorship at https://rdrr.io/rforge/spacemakeR/

As completely off topic but fun to share:

I do struggle to follow along with the requirements of the package. Notably the

  1. the need for symmetric matrix
  2. doubly centering

The approach taken in spacemakeR requires a dense matrix due to the symmetry and the double centering. Though I'm not convinced that it notably improved the eigenvectors. Here is an example of creating the first 15 eigenvectors using a partial eigenvector from irlba. It's ridiculously fast and scalable. Using this approach and including the resultant vectors into an OLS model does noticeably reduce the spatial autocorrelation

Dray '06 writes

PCoA is usually computed on a distance matrix but Gower (1966) showed that this analysis can also be computed from a similarity matrix

In general, I think we can think of a SWM as a form of similarity matrix which is followed up as

A more interesting interpretation is to consider S to be a spatial weighting matrix (Bavaud, 1998, Tiefelsdorf et al., 1999), which indicates the strength of the potential interaction among the spatial units. In the PCoA of D (Fig. 1a), all sites are considered to be neighbours except those corresponding to the largest distance, for which the similarity is null

I think leaning into the sparsity of spatial weights matrix emphasized

"Our interpretation of PCNM base functions as a particular case of MEM generalizes the original approach because “the use of a generalised weighting matrix [. . .] allows the investigator to choose a set of weights which he deems appropriate from prior consider- ations. This allows great flexibility” (Cliff and Ord, 1973, p. 12)." (p 485, Dray et al 2006)

# install.packages("tripack")
# install.packages('ade4')
# install.packages("spacemakeR", repos="http://R-Forge.R-project.org")
library(sfdep)
library(spdep)
library(spData)
library(spacemakeR)
library(ggplot2)
library(patchwork)

# using my approach without symmetry and doing partial eigen 
partial_mem <- function(listw, n = 15) {
  requireNamespace("spatialreg")
  m <- as(listw, "CsparseMatrix")
  irlba::partial_eigen(m, n, symmetric = is.symmetric.nb(listw$neighbours))
}

visualize_partial_mems <- function(x, listw, n = 15) {
  mems <- partial_mem(listw, n)
  plots <- vector("list", n)

  plots <- lapply(1:n, \(i) {
    z <- mems$vectors[,i]
    mi <- moran.test(z, listw)
    ggplot(x, aes(color = z)) +
      geom_sf() +
        scale_color_viridis_c() +
        labs(title = paste0("Moran's I: ", round(mi$estimate[[1]], 3))) +
        theme_void() + 
        theme(legend.position = "none")
  })

  list(
    plots = patchwork::wrap_plots(plots),
    evs = mems
  )
}

# symmetric 
nb <- st_knn(cycle_hire$geometry, 10)
wt <- st_kernel_weights(nb, cycle_hire$geometry, "gaussian", TRUE)
listw <- nb2listw(nb, wt, "B")
res <- visualize_partial_mems(cycle_hire, listw, 10)

res$plots

rsbivand commented 2 weeks ago

Just this now: did you also look at https://cran.r-project.org/package=adespatial, which he wrote after spacemakeR if I remember correctly?