josiahparry / sfdep Goto Github PK

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
}

guerry <- sfdep::guerry %>% 
  st_set_crs(27572) %>% 
  select(code_dept, crime_pers, crime_prop)

guerry_nb <- guerry %>% 
  mutate(nb = st_contiguity(geometry),
         wt = st_weights(nb))

guerry_lisa <- guerry_nb %>% 
  mutate(local_moran_crime_pers = local_moran(crime_pers, nb, wt, nsim = 199), # error occurs here
         local_moran_crime_prop = local_moran(crime_prop, nb, wt, nsim = 199)) # tried with another column as well, just in case

Error in `stopifnot()`:
! Problem while computing `local_moran_crime_pers = local_moran(crime_pers, nb, wt, nsim = 199)`.
Caused by error in `names(object) <- nm`:
! 'names' attribute [9] must be the same length as the vector [8]

library(sf)
library(sfdep)
library(simplevis)

# replicate the guerry dataset 10 times
x <- purrr::map_dfr(1:10, ~guerry)  %>%  
  dplyr::select(code_dept, crime_pers) %>% 
  # create an indicator for time period
  mutate(time_period = sort(rep(1:10, 85)), 
         # add some noise 
         crime_pers = crime_pers * runif(850, max = 2))

x

# check what this looks like for one time period
x %>% 
  filter(time_period == 1) %>% 
  st_set_crs(27572) %>% 
  simplevis::gg_sf_col(col_var = crime_pers)

# test using emerging_hotspot_analysis() which first requires creating a spacetime object...
spt <- as_spacetime(x, "code_dept", "time_period")

ehsa_guerry <-
  emerging_hotspot_analysis(spt,
                            "crime_pers",
                            include_gi = TRUE, # note the include_gi argument doesn't seem to work.
                            threshold = 0.05)  %>% 
  mutate(classification = factor(
    classification,
    levels = c(
      "consecutive hotspot",
      "sporadic hotspot",
      "new hotspot",
      "new coldspot",
      "sporadic coldspot",
      "consecutive coldspot",
      "persistent coldspot",
      "no pattern detected"
      ))) 

rcartocolor::carto_pal(n = 9, name = "Prism") %>% 
  scales::show_col()

p <- c(rev(c("#CC503E", "#E17C05", "#EDAD08", "#73AF48", "#0F8554", "#38A6A5", "#1D6996")),"grey")

guerry %>% 
  left_join(ehsa_guerry, by = c("code_dept" = "location")) %>% 
  st_set_crs(27572) %>% 
  simplevis::gg_sf_col(col_var = classification,
                       pal = p,
                       title = "Emerging hotspots - using sfdep::emerging_hotspot_analysis()")

# compare what the output from emerging_hotspot_analysis() looks like relative to doing it step-by-step----
# identify neighbours and assign weights----
x_ <- x %>%
  mutate(nb = include_self(st_contiguity(geometry)),
         wt = st_weights(nb)) %>% 
  split(list(.$time_period))

# calculate local_gstar----
guerry_gi_star <- x_ %>% 
  map(. %>% mutate(nb = include_self(st_contiguity(geometry)),
                   wt = st_weights(nb)) %>% transmute(gi_star = local_gstar(crime_pers, nb, wt))) %>% 
  map(. %>% tidyr::unnest(gi_star)) %>% 
  map2(names(.), ~ mutate(.x, time_period = .y)) %>% 
  map(. %>% mutate(code_dept = x_[[1]]$code_dept)) %>% 
  map(. %>% relocate(code_dept, everything())) %>% 
  bind_rows()

# map it for one time period----
guerry_gi_star %>% 
  st_set_crs(27572) %>% 
  filter(time_period == "1") %>% 
  simplevis::gg_sf_col(col_var = gi_star)
  
# and now calculate the trends----
guerry_gi_star_l <- guerry_gi_star %>% 
  st_drop_geometry() %>% 
  split(list(.$code_dept))

mk <- guerry_gi_star_l %>%
  furrr::future_map( ~ .x %>%
                       pull(gi_star) %>%
                       mmkh() %>% as.data.frame() %>% rownames_to_column(var = "output") %>% pivot_wider(names_from = "output",
                                                                                                         values_from = ".")
  )

mk_ <- mk %>%
  bind_rows(.id = "code_dept") %>%
  dplyr::select(-c(`Original Z`, `old P.value`, old.variance)) %>%
  janitor::clean_names() %>% 
  mutate(code_dept = as.factor(code_dept))

mk_

guerry_gi_star_mk <- guerry_gi_star %>%
  left_join(mk_, by = c("code_dept"))

# categories adapted from https://pro.arcgis.com/en/pro-app/latest/tool-reference/space-time-pattern-mining/learnmoreemerging.htm----
guerry_gi_star_mk_cat <- guerry_gi_star_mk %>%
  mutate(
    hotspot = case_when(
      gi_star <= -2.58 ~ "Cold spot 99%",
      gi_star > -2.58 & gi_star <= -1.645 ~ "Cold spot 90%",
      gi_star > -1.645 & gi_star <= -0.9542 ~ "Cold spot 66%",
      gi_star >= 0.9542 & gi_star < 1.645 ~ "Hot spot 66%",
      gi_star >= 1.645 & gi_star < 2.58 ~ "Hot spot 90%",
      gi_star >= 2.58 ~ "Hot spot 99%",
      TRUE ~ "Not Significant"
    )
  ) %>% # adapted from https://crd230.github.io/lab7.html#Local_spatial_autocorrelation and adapted to align this with IPCC thresholds
  mutate(hotspot = factor(
    hotspot,
    levels = c(
      "Cold spot 99%",
      "Cold spot 90%",
      "Cold spot 66%",
      "Not Significant",
      "Hot spot 66%",
      "Hot spot 90%",
      "Hot spot 99%"
    )
  )) %>%
  group_by(code_dept) %>%
  arrange(code_dept, time_period) %>%
  mutate(
    hot = case_when(hotspot %in% c(# "Hot spot 66%", # remove to aligned with sfdep::emerging_hotspot_analysis() which seems to use a default threshold of 0.01.
      "Hot spot 90%",
      "Hot spot 99%") ~ 1,
      TRUE ~ 0),
    cold = case_when(hotspot %in% c(# "Cold spot 66%", # remove to aligned with sfdep::emerging_hotspot_analysis() which seems to use a default threshold of 0.01.
      "Cold spot 90%",
      "Cold spot 99%") ~ 1,
      TRUE ~ 0),
    prev_hot = lag(hot),
    prev_cold = lag(cold),
    n_hot = sum(hot),
    p_hot = n_hot / length(hot),
    n_cold = sum(cold),
    p_cold = n_cold / length(cold),
    next_h = case_when(hot != lag(hot) ~ 1),
    n_changes_hot = sum(next_h, na.rm = TRUE),
    next_c = case_when(cold != lag(cold) ~ 1),
    n_changes_cold = sum(next_c, na.rm = TRUE),
    cat = case_when(
      last(hot) == 1 &
        n_hot == 1 ~ "new hotspot",
      # A location that is a statistically significant hot spot for the final time step and has never been a statistically significant hot spot before.
      last(hot) == 1 &
        n_hot == 2 &
      p_hot < 0.9 ~ "consecutive hotspot",
      # A location with a single uninterrupted run of at least two statistically significant hot spot bins in the final time-step intervals. The location has never been a statistically significant hot spot prior to the final hot spot run and less than 90 percent of all bins are statistically significant hot spots.
      last(hot) == 1 &
        n_changes_hot >= 3 &
        p_hot < 0.9 &
        n_cold == 0 ~ "sporadic hotspot",
      # A statistically significant hot spot for the final time-step interval with a history of also being an on-again and off-again hot spot (i.e. there are not more than two consecutive hotspots). Less than 90 percent of the time-step intervals have been statistically significant hot spots and none of the time-step intervals have been statistically significant cold spots.
      p_hot >= 0.9 &
        last(hot) == 1 &
        tau > 0 &
        new_p_value <= 0.333 ~ "intensifying hotspot",
      # A location that has been a statistically significant hot spot for 90 percent of the time-step intervals, including the final time step. In addition, the intensity of clustering of high counts in each time step is increasing overall and that increase is statistically significant.
      p_hot >= 0.9 &
        new_p_value > 0.333 ~ "persistent hotspot",
      # A location that has been a statistically significant hot spot for 90 percent of the time-step intervals with no discernible trend in the intensity of clustering over time.
      p_hot >= 0.9 &
        new_p_value <= 0.333 &
        tau < 0 &
        last(hot == 1) ~ "diminishing hotspot",
      # A location that has been a statistically significant hot spot for 90 percent of the time-step intervals, including the final time step. In addition, the intensity of clustering in each time step is decreasing overall and that decrease is statistically significant.
      last(hot == 0) & p_hot >= 0.9 ~ "historic hotspot",
      # The most recent time period is not hot, but at least 90 percent of the time-step intervals have been statistically significant hot spots.
      last(cold) == 1 &
        n_cold == 1 ~ "new coldspot",
      # A location that is a statistically significant cold spot for the final time step and has never been a statistically significant cold spot before.
      last(hot) == 1 &
        n_hot == 2 &
      p_cold < 0.9 ~ "consecutive coldspot",
      # A location with a single uninterrupted run of at least two statistically significant cold spot bins in the final time-step intervals. The location has never been a statistically significant cold spot prior to the final cold spot run and less than 90 percent of all bins are statistically significant cold spots.
      last(cold) == 1 &
        n_changes_cold >= 3 &
        p_cold < 0.9 &
        n_hot == 0 ~ "sporadic coldspot",
      # A statistically significant cold spot for the final time-step interval with a history of also being an on-again and off-again cold spot (i.e. there are not more than two consecutive coldspots). Less than 90 percent of the time-step intervals have been statistically significant cold spots and none of the time-step intervals have been statistically significant hot spots.
      p_cold >= 0.9 &
        last(cold) == 1 &
        tau < 0 &
        new_p_value <= 0.333 ~ "intensifying coldspot",
      # A location that has been a statistically significant cold spot for 90 percent of the time-step intervals, including the final time step. In addition, the intensity of clustering of low counts in each time step is increasing overall and that increase is statistically significant.
      p_cold >= 0.9 &
        new_p_value > 0.333 ~ "persistent coldspot",
      # A location that has been a statistically significant cold spot for 90 percent of the time-step intervals with no discernible trend in the intensity of clustering over time.
      p_cold >= 0.9 &
        new_p_value <= 0.333 &
        tau < 0 &
        last(cold == 1) ~ "diminishing coldspot",
      # A location that has been a statistically significant cold spot for 90 percent of the time-step intervals, including the final time step. In addition, the intensity of clustering in each time step is decreasing overall and that decrease is statistically significant.
      last(cold == 0) & p_cold >= 0.9 ~ "historic coldspot"
    ),
    # The most recent time period is not cold, but at least 90 percent of the time-step intervals have been statistically significant cold spots.
  ) %>%
  mutate(cat = replace_na(cat, "no pattern detected")) %>%
  ungroup() %>%
  mutate(cat = factor(
    cat,
    levels = c(
      "consecutive hotspot",
      "sporadic hotspot",
      "new hotspot",
      "new coldspot",
      "sporadic coldspot",
      "consecutive coldspot",
      "persistent coldspot",
      "no pattern detected"
    )
  )) 

# map it----
guerry_gi_star_mk_cat %>% 
  group_by(code_dept) %>% 
  slice(1) %>% 
  st_set_crs(27572) %>% 
  simplevis::gg_sf_col(col_var = cat,
                             pal = p,
                             title = "Emerging hotspots - manually using sfdep::local_gstar, mann_kendall, etc.")

# mmm... so there appear to be lots of differences between the two methods... Though it is unclear how emerging_hotspot_analysis() categorises the data so a more robust comparision is to look at the tau and p-values (noting that the include_gi = TRUE argument doesn't seem to be working so we can't compare gi_star values)

mk <- guerry_gi_star_l %>%
  furrr::future_map( ~ .x %>%
                       values_from = ".")
                       mutate(mk = mann_kendall(gi_star, alternative = "greater")))

mk_ <- mk %>%
  bind_rows(.id = "code_dept") %>%
  janitor::clean_names() %>% 
  mutate(code_dept = as.factor(code_dept))

mk_ %>% 
  group_by(code_dept) %>% 
  slice(1) %>% 
  unnest() %>% 
  dplyr::select(code_dept, gi_star, p_value, tau) %>% 
  left_join(ehsa_guerry, by = c("code_dept" = "location")) %>% # note I'm pretty sure that emerging_hotspot_analyais() calls MannKendall() which returns a two-sided p-value.
  ungroup() %>% 
  summarise(sbs_n_p_less_05 = sum(p_value.x < 0.05),
            ehsa_n_p_less_05 = sum((p_value.y/2) < 0.05), # I'm dividing by 2 to turn a two-sided p-value into a one-sided p-value
            dif_tau = sum(tau.x != tau.y),
            n_sbs_positive_tau = sum(tau.x > 0),
            n_sbs_negative_tau = sum(tau.x < 0),
            n_ehsa_positive_tau = sum(tau.y > 0),
            n_ehsa_negative_tau = sum(tau.y < 0),
            n_diff_pos_neg = sum(tau.x > 0 & tau.y < 0 | tau.x < 0 & tau.y > 0 ),
            n_dif_ehsa_sig = sum((p_value.y/2) < 0.05 & (tau.x > 0 & tau.y < 0 | tau.x < 0 & tau.y > 0 )))

# note it is unclear whether the ehsa_guerry p_value relates to the gi_star or mannKendall result... I suspect is relates to the gi_star...

# So there are 16 code_dept(s) that ehsa finds have a significant p-value (< 0.05) but that have different tau values. This seems to suggest that over half the estimated slopes might be wrong??? That said it is apparent from below that there are only five polygons/code_dept that have a statistically significant (p-value <= 0.1) trend and of these five none had more than 0.5 of the timeseries that were hotspots (they needed to have at least 0.9 of the timesteps as hotspots). Hence none of them were intensifying or diminishing.

guerry_gi_star_mk_cat %>% 
  filter(new_p_value <= 0.05) %>% 
  # distinct(code_dept)
  View()

# ok so there are only five polygons/code_dept that have a statistically significant (p-value <= 0.1) trend and of these five none had more than 0.5 of the timeseries that were hotspots. Hence none of them were intensifying or diminishing.