The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Epidemiology with pkmapr

library(pkmapr)
library(spdep)
library(dplyr)

The problem: Non-contiguous units

Gilgit-Baltistan and Azad Kashmir share no land boundaries with other provinces. Under standard contiguity rules, they would have zero neighbours — breaking most spatial statistics.

pk_neighbors() handles this automatically.

Build spatial weights

districts <- get_districts()
w <- pk_neighbors(districts, disputed = "exclude")

The returned object contains: - w$nb: neighbours list - w$listw: row-standardised weights (ready for spdep)

Global Moran’s I (synthetic data)

set.seed(2023)

# Create spatially autocorrelated synthetic data
n <- nrow(districts)
rho <- 0.6
W_mat <- listw2mat(w$listw)
epsilon <- rnorm(n, mean = 50, sd = 10)
y <- as.numeric(solve(diag(n) - rho * W_mat) %*% epsilon)

districts$synthetic_rate <- y

# Test for spatial autocorrelation
moran_result <- moran.test(districts$synthetic_rate, w$listw)
print(moran_result)

A significant p-value indicates spatial clustering.

Local Indicators of Spatial Association (LISA)

Identify local clusters and outliers:

# Calculate local Moran's I
lisa <- localmoran(districts$synthetic_rate, w$listw)

# Add results to districts
districts$lisa_i <- lisa[, 1]  # local I
districts$lisa_p <- lisa[, 5]  # p-value

# Classify clusters (simplified)
districts$cluster <- case_when(
  districts$lisa_p > 0.05 ~ "Not significant",
  districts$synthetic_rate > mean(districts$synthetic_rate) &
    districts$lisa_i > 0 ~ "High-High",
  districts$synthetic_rate < mean(districts$synthetic_rate) &
    districts$lisa_i > 0 ~ "Low-Low",
  districts$lisa_i < 0 ~ "Outlier",
  TRUE ~ "Other"
)

# Map the clusters
ggplot2::ggplot(districts) +
  ggplot2::geom_sf(ggplot2::aes(fill = cluster)) +
  ggplot2::theme_void() +
  ggplot2::labs(title = "Spatial Clusters")

Hotspot detection (Getis-Ord Gi*)

# Calculate Gi* for the synthetic data
gi_star <- localG(districts$synthetic_rate, w$listw)

districts$gi_star <- as.numeric(gi_star)
districts$hotspot <- ifelse(districts$gi_star > 1.96, "Hotspot",
                            ifelse(districts$gi_star < -1.96, "Coldspot", "Not significant"))

pk_map(districts, fill = "hotspot", title = "Hotspots (p < 0.05)")

Disputed boundary handling

The Line of Control creates analytical ambiguity. pk_neighbors() makes your decision explicit:

# Exclude (default) — GB/AJK get nearest neighbour fallback
w_exclude <- pk_neighbors(districts, disputed = "exclude")

# Include — treat disputed boundaries as shared
w_include <- pk_neighbors(districts, disputed = "include")

# Flag — document which units are affected
w_flag <- pk_neighbors(districts, disputed = "flag")
print(w_flag$boundary_note)

Complete workflow example

# 1. Get data
districts <- get_districts()

# 2. Build weights with explicit dispute handling
w <- pk_neighbors(districts, disputed = "exclude")

# 3. Attach your real case data (replace with your own)
# districts <- pk_join(districts, your_case_data, by = "district_code")

# 4. Calculate rates (example)
# districts <- districts |>
#   mutate(rate = cases / population * 100000)

# 5. Test for spatial autocorrelation
# moran.test(districts$rate, w$listw)

# 6. Identify clusters
# lisa <- localmoran(districts$rate, w$listw)
# districts$lisa_cluster <- attr(lisa, "quadr")$mean

# 7. Map results
# pk_map(districts, fill = "lisa_cluster")

References

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.