## ----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(
#    cases       = rj_mortality$cases,
#    population  = rj_mortality$live_births,
#    region_id   = rj_mortality$region_id,
#    x           = rj_mortality$x,
#    y           = rj_mortality$y,
#    node_id     = rj_mortality$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(
#    cases       = rj_agg$cases,
#    population  = rj_agg$live_births,
#    region_id   = rj_agg$region_id,
#    x           = rj_agg$x,
#    y           = rj_agg$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(geobr)
#  library(sf)
#  
#  mun       <- read_municipality(code_muni = "RJ", year = 2016)
#  mun$code6 <- as.integer(substr(mun$code_muni, 1, 6))
#  
#  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(mun, cr2, by.x = "code6", by.y = "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(
#    cases       = rj_mortality$cases,
#    population  = rj_mortality$live_births,
#    region_id   = rj_mortality$region_id,
#    x           = rj_mortality$x,
#    y           = rj_mortality$y,
#    node_id     = rj_mortality$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")]