## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  eval = FALSE
)

## ----rj-data------------------------------------------------------------------
# library(treeSS)
# data(rj_mortality)
# data(rj_tree)
# 
# c(rows           = nrow(rj_mortality),
#   municipalities = length(unique(rj_mortality$region_id)),
#   tree_nodes     = nrow(rj_tree),
#   total_deaths   = sum(rj_mortality$cases),
#   total_births   = sum(unique(
#     rj_mortality[, c("region_id", "live_births")]
#   )$live_births))
# #> rows           municipalities  tree_nodes  total_deaths  total_births
# #> 1440           89              622         2981          219124

## ----rj-scan------------------------------------------------------------------
# result_rj <- treespatial_scan(
#   rj_mortality,
#   cases       = cases,
#   population  = live_births,
#   region_id   = region_id,
#   x           = x,
#   y           = y,
#   node_id     = node_id,
#   tree        = rj_tree,
#   max_pop_pct = 0.50,
#   nsim        = 999, seed = 2016,
#   n_cores     = 4L
# )
# print(result_rj)

## ----rj-spatial-only----------------------------------------------------------
# rj_agg <- unique(rj_mortality[, c("region_id", "x", "y", "live_births")])
# rj_agg$cases <- tapply(rj_mortality$cases,
#                         rj_mortality$region_id,
#                         sum)[as.character(rj_agg$region_id)]
# 
# result_spatial <- circular_scan(
#   rj_agg,
#   cases       = cases,
#   population  = live_births,
#   region_id   = region_id,
#   x           = x,
#   y           = y,
#   max_pop_pct = 0.50,
#   nsim        = 999, seed = 2016
# )
# result_spatial$most_likely_cluster$region_ids

## ----rj-filter----------------------------------------------------------------
# filter_clusters(result_rj)[1:5,
#                             c("node_id", "n_regions",
#                               "cases", "expected", "llr")]

## ----rj-getregions------------------------------------------------------------
# # Most likely cluster only (single map)
# cr1 <- get_cluster_regions(result_rj, n_clusters = 1, overlap = FALSE)
# 
# # Top-2 distinct clusters (faceted map)
# cr2 <- get_cluster_regions(result_rj, n_clusters = 2, overlap = TRUE)

## ----rj-plot------------------------------------------------------------------
# library(ggplot2)
# library(sf)
# 
# data(rj_map)   # sf polygons: ibge_code, NAME, geometry
# 
# cr2 <- merge(get_cluster_regions(result_rj, n_clusters = 2,
#                                   overlap = TRUE),
#              unique(rj_mortality[, c("region_id", "ibge_code")]),
#              by = "region_id")
# mun_facet <- merge(rj_map, cr2, by = "ibge_code")
# 
# ggplot(mun_facet) +
#   geom_sf(aes(fill = factor(cluster)), color = "gray40",
#           alpha = 0.75) +
#   scale_fill_manual(values = c("1" = "#C44E52", "2" = "#4C72B0"),
#                     na.value = "gray95", na.translate = FALSE) +
#   facet_wrap(~ panel, nrow = 1) +
#   theme_minimal() +
#   theme(legend.position = "none")

## ----rj-seq-------------------------------------------------------------------
# seq_rj <- sequential_scan(
#   rj_mortality,
#   cases       = cases,
#   population  = live_births,
#   region_id   = region_id,
#   x           = x,
#   y           = y,
#   node_id     = node_id,
#   tree        = rj_tree,
#   max_iter    = 5L, alpha = 0.05,
#   nsim        = 999, seed = 2016,
#   max_pop_pct = 0.50, n_cores = 4L
# )
# print(seq_rj)

## ----rj-seq-table-------------------------------------------------------------
# seq_rj$clusters[, c("iteration", "node_id", "n_regions",
#                     "llr", "pvalue", "significant")]