inst/doc/introduction.R

In treeSS: Tree-Spatial Scan Statistic for Cluster Detection

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

## ----rj-data------------------------------------------------------------------
#  library(treeSS)
#  data(rj_mortality)
#  data(rj_tree)
#  
#  str(rj_mortality, max = 1)
#  #> 'data.frame': 1440 obs. of 9 variables: region_id, ibge_code,
#  #>   name, live_births, x, y, node_id, cases, ...

## ----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-filter----------------------------------------------------------------
#  filter_clusters(result_rj)[1:5, c("node_id", "n_regions",
#                                     "cases", "expected", "llr")]

## ----rj-getregions, eval = FALSE----------------------------------------------
#  # 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")

## ----chi-scan-----------------------------------------------------------------
#  data(chicago_crimes)
#  data(chicago_tree)
#  
#  result_chi <- treespatial_scan(
#    cases       = chicago_crimes$cases,
#    population  = chicago_crimes$pop_residential,
#    region_id   = chicago_crimes$region_id,
#    x           = chicago_crimes$x,
#    y           = chicago_crimes$y,
#    node_id     = chicago_crimes$node_id,
#    tree        = chicago_tree,
#    max_pop_pct = 0.25,
#    nsim        = 999, seed = 2023,
#    n_cores     = 4L
#  )
#  print(result_chi)

## ----chi-plot-----------------------------------------------------------------
#  data(chicago_map)
#  
#  cr1 <- merge(get_cluster_regions(result_chi, n_clusters = 1, overlap = FALSE),
#               unique(chicago_crimes[, c("region_id", "area_number", "name")]),
#               by = "region_id")
#  chi <- merge(chicago_map, cr1, by.x = "AREA_NUM", by.y = "area_number",
#                all.x = TRUE)
#  
#  ggplot(chi) +
#    geom_sf(aes(fill = factor(cluster)), color = "gray40", alpha = 0.75) +
#    scale_fill_manual(values = c("1" = "#C44E52"), na.value = "gray95",
#                      name = "Cluster") +
#    theme_minimal()

## ----rj-iter------------------------------------------------------------------
#  iter_rj <- iterative_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    = 5, alpha = 0.05,
#    nsim        = 999, seed = 2016,
#    max_pop_pct = 0.50, n_cores = 4L
#  )
#  print(iter_rj)

## ----rj-iter-table------------------------------------------------------------
#  iter_rj$clusters[, c("iteration", "node_id", "n_regions",
#                        "llr", "pvalue", "pvalue_adjusted",
#                        "significant")]

## ----rj-iter-map--------------------------------------------------------------
#  cr_it <- merge(get_cluster_regions(iter_rj, overlap = TRUE),
#                 unique(rj_mortality[, c("region_id", "ibge_code")]),
#                 by = "region_id")
#  mun_it <- merge(mun, cr_it, by.x = "code6", by.y = "ibge_code",
#                  all.x = TRUE)
#  
#  ggplot(mun_it) +
#    geom_sf(aes(fill = factor(cluster)), color = "gray40", alpha = 0.75) +
#    scale_fill_manual(values = c("1" = "#C44E52", "2" = "#4C72B0",
#                                  "3" = "#55A868", "4" = "#8172B2",
#                                  "5" = "#CCB974"),
#                      na.value = "gray95", na.translate = FALSE) +
#    facet_wrap(~ panel) +
#    theme_minimal() + theme(legend.position = "none")

## ----examples-dir-------------------------------------------------------------
#  list.files(system.file("examples", package = "treeSS"))
#  #> [1] "example_brazil_rj.R"  "example_chicago.R"
#  #> [2] "example_florida.R"    "example_london.R"