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:

spdep/R/localmoran.R

Lines 45 to 56 in 1fb4943

	lz <- lag.listw(listw, z, zero.policy=zero.policy, NAOK=NAOK)
	lbs <- c("Low", "High")
	quadr_ps <- interaction(cut(z, c(-Inf, 0, Inf), labels=lbs),
	cut(lz, c(-Inf, 0, Inf), labels=lbs), sep="-")
	lx <- lag.listw(listw, x, zero.policy=zero.policy, NAOK=NAOK)
	lxx <- mean(lx, na.rm=NAOK)
	quadr <- interaction(cut(x, c(-Inf, xx, Inf), labels=lbs),
	cut(lx, c(-Inf, lxx, Inf), labels=lbs), sep="-")
	xmed <- median(x, na.rm=NAOK)
	lxmed <- median(lx, na.rm=NAOK)
	quadr_med <- interaction(cut(x, c(-Inf, xmed, Inf), labels=lbs),
	cut(lx, c(-Inf, lxmed, Inf), labels=lbs), sep="-")

. 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 0152f72 and PR #156

[rsbivand]

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

[rsbivand]

895c520 adds bi-variate Moran plot.