cannot classify Bivariate Moran Hotspots

[JosiahParry]

I am trying to use the hotspot() function to classify Bivariate Local Moran output but this results in an error

library(sf)
library(spdep)
columbus <- st_read(system.file("shapes/columbus.shp", package="spData"))

nb <- poly2nb(columbus)
listw <- nb2listw(nb)
set.seed(1)
(res <- localmoran_bv(columbus$CRIME, columbus$INC, listw, nsim = 499))

hotspot(res, Prname = "Pr(folded) Sim", cutoff = 0.05, p.adjust = "none")

Error in `UseMethod()`:
! no applicable method for 'droplevels' applied to an object of class "logical"

[rsbivand]

@JosiahParry the object returned by localmoran_bv does not have an "quadr" attribute; in localmoran it is constructed as: https://github.com/r-spatial/spdep/blob/1fb494364507cc9539901271c202cad69e61f324/R/localmoran.R#L45-L56. Would you like to suggest how it might be added? The pysal variant is like GeoDa (split on zero for z and lz), median was a speculation. For the present I'm adding a better error message.

[rsbivand]

Now says:

> hotspot(res, Prname = "Pr(folded) Sim", cutoff = 0.05, p.adjust = "none")
Error in hotspot.localmoran(res, Prname = "Pr(folded) Sim", cutoff = 0.05,  : 
  object has no quadr attribute

[JosiahParry]

Thanks Roger! My apologies for the delay in getting back to you I was stuck in airports for far too long!

I was having a think about this and how we could calculate the quadrant of the cluster. Given that the variables are not guaranteed to be on the same scale, I think one has to center and scale them so that they can be effectively compared. Here's how I might go about it

library(spdep)
#> Loading required package: spData
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://nowosad.github.io/drat/', type='source')`
#> Loading required package: sf
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(sfdep)
library(ggplot2)

# constants
zero.policy=NULL
NAOK=FALSE

# SWM
listw <- nb2listw(
  poly2nb(guerry)
)

# vars 
x <- guerry$literacy
y <- guerry$instruction

# note the replacement of z with scale(y)
ly <- lag.listw(listw, scale(y), zero.policy=zero.policy, NAOK=NAOK) 
lbs <- c("Low", "High") 

quadr_ps <- interaction(
  cut(scale(x), c(-Inf, 0, Inf), labels=lbs),
  cut(ly, c(-Inf, 0, Inf), labels=lbs),
  sep="-"
) 

# Moran plot
data.frame(
  x = scale(x),
  ly,
  clust = quadr_ps
) |> 
  ggplot(aes(x, ly, color = clust)) + 
  geom_point(size = 3) 

# calculate bv moran
bv_res <- localmoran_bv(x, y, listw, 499, TRUE) |> 
  as.data.frame()

# using 0.05 for example purposes
sig <- bv_res[["Pr(z != E(Ibvi)) Sim"]] < 0.05

# removing "insignificant" clusers 
quadr_ps[!sig] <- NA

# visualize
ggplot() +
  geom_sf(aes(fill = quadr_ps), guerry, lwd = 0.15, color = "black") +
  scale_fill_viridis_d()

[rsbivand]

@JosiahParry Could you please review https://github.com/r-spatial/spdep/commit/0152f720569d8a9b5e0c278a9d4c1e5651e4b504 and PR #156

[rsbivand]

@JosiahParry Should I modify moran.plot to take an optional y argument for the bivariate case?

[rsbivand]

https://github.com/r-spatial/spdep/pull/156/commits/895c5204a5d88bd3cb36333b17013e72ecb056dd adds bi-variate Moran plot.