Add categorize LISA function for labeling lisa clusters

[JosiahParry]

This was previously implemented in sfweight and would be a good utility to make avaialble to people

function(x, x_lag, scale = TRUE) {

  cats <- character(length(x))

  if (scale) {
    x <- scale(x)
    x_lag <- scale(x_lag)

    cats[x > 0 & x_lag > 0] <- "HH"
    cats[x > 0 & x_lag < 0] <- "HL"
    cats[x < 0 & x_lag < 0] <- "LL"
    cats[x < 0 & x_lag > 0] <- "LH"

  }

  cats[x > mean(x) & x_lag > mean(x_lag)] <- "HH"
  cats[x > mean(x) & x_lag < mean(x_lag)] <- "HL"
  cats[x < mean(x) & x_lag < mean(x_lag)] <- "LL"
  cats[x < mean(x) & x_lag > mean(x_lag)] <- "LH"
  cats[cats == ""] <- NA
  cats
}

[denis-or]

Hi @JosiahParry . Thanks for the package.

Just a question: this code work for to bivariate_moran function too?

[JosiahParry]

Uhm. I have to think about it. I don't....know. You'd have to ensure that both are scaled. The bivariate moran is somewhat of a....not phenomenal statistic. The idea behind it is that you're comparing a variable x to the spatial lag of a variable y. So the inference about it is x is spatially dependent upon the neighborhood values of y.

But with that said. Yeah, if you replace x_lag with y_lag it should aboslutely work. I'd transform x and y into z scores first using scale() and then use this.

[JosiahParry]

@denis-or check out the documentation for Geoda at https://geodacenter.github.io/workbook/5b_global_adv/lab5b.html for a good example of the bivariate (and theory) behind it

[denis-or]

Oh right. Thanks. I saw this documentation and I liked it. I would like to use sfdep because it is possible to use pipe as well. lol. But thanks again.

[JosiahParry]

What do you mean? Just suggesting you read the documentation for the bivariate moran and bv moran plot on how to do it in R the above function template is agreat start

[denis-or]

What do you mean? Just suggesting you read the documentation for the bivariate moran and bv moran plot on how to do it in R the above function template is agreat start

Yes yes! You're right. It's a great start for sure. I even wrote something from the documentation. I would like to apologize, for this misunderstanding of mine. My intention was just to confirm if it is possible to use pipe for all code.

bv_moran <- function (x, y, nb, wt, nsim = 999, alpha = 0.05, scale = TRUE) 
{
  listw <- sfdep::recreate_listw(nb, wt)

  data_bv <- sfdep:::local_moran_bv_impl(x, y, listw, nsim = nsim)

  p_sim <- data_bv['p_sim']

  lag_y <- sfdep::st_lag(y, nb, wt)

  cats <- character(length(x))

  if (scale) {
    x <- scale(x)
    y <- scale(lag_y)

    cats[x > 0 & y > 0] <- "HH"
    cats[x > 0 & y < 0] <- "HL"
    cats[x < 0 & y < 0] <- "LL"
    cats[x < 0 & y > 0] <- "LH"
    cats[p_sim >= alpha] <- "NS"

  }

  cats[x > mean(x) & lag_y > mean(lag_y)] <- "HH"
  cats[x > mean(x) & lag_y < mean(lag_y)] <- "HL"
  cats[x < mean(x) & lag_y < mean(lag_y)] <- "LL"
  cats[x < mean(x) & lag_y > mean(lag_y)] <- "LH"
  cats[cats == ""] <- NA
  cats[p_sim >= alpha] <- "NS"

  data_bv$pysal_bv <- cats

  data_bv
}

g <- guerry %>% 
  mutate(
    nb = st_contiguity(geometry),
    wt = st_weights(nb),
    moran_bv = bv_moran(crime_pers, wealth, nb, wt)
  )

g %>% 
  tidyr::unnest(moran_bv) %>% 
  ggplot(aes(fill = pysal_bv)) +
  geom_sf() +
  geom_sf(lwd = 0.2, color = "black") +
  theme_void() +
  scale_fill_manual(
    name = "With sfdep",
    values = c("#1C4769", "#24975E",
               "#EACA97", "#B20016",
               "#FAFAFA"), drop = F)

queen_w <- rgeoda::queen_weights(guerry)

lisa <- rgeoda::local_bimoran(queen_w,
                      guerry[c("crime_pers","wealth")])

g2 <- mutate(guerry,
       bv_moran = factor(rgeoda::lisa_clusters(lisa),
                      levels = 0:6,
                      labels = rgeoda::lisa_labels(lisa)))
g2 %>%
ggplot(aes(fill = bv_moran)) +
  geom_sf() +
  geom_sf(lwd = 0.2, color = "black") +
  theme_void() +
  scale_fill_manual(
    name = "With rgeoda",
    values = c("#FAFAFA","#1C4769","#B20016",
                               "#EACA97","#24975E","#650A01",
                               "#231300"), drop = F)