Welcome to Spatial Disease Mapping!

In epidemiology, disease mapping helps us visualize geographic variations in disease risk. Instead of just looking at raw case counts, we use statistical models to smooth out random fluctuations and get a clear, reliable picture of where risk is truly elevated.
We will be using R-INLA (Integrated Nested Laplace Approximations), a powerful tool that allows us to fit complex Bayesian spatial models much faster than traditional simulation methods.

setwd("C:/Users/lenovo/Desktop/Rpubs/Rinla")
# 1. Install and Load Required Packages

library(INLA) #Our core engine for Bayesian modeling
library(spData)
library(sf) #Simple Features: The modern standard for working with geographic boundaries (shapes)
library(spdep) #A package used to figure out which geographic areas touch each other
library(ggplot2)
library(report)

Import Spatial Data from a Local Shapefile

To make a map, R needs to know the geometric boundaries of the areas we are studying (such as counties or sub-counties). This information is stored in a file called a Shapefile (.shp).

Note: For a shapefile to load properly, make sure its sibling files (.shx, .dbf, and .prj) are kept in the exact same folder!

We will also create a region_id. R-INLA doesn’t understand text names (like “Nairobi” or “Mombasa”) for its spatial calculations; it requires a unique, sequential tracking number (1, 2, 3…) for each region.

# 2. Import Spatial Data from Local Shapefile
# Replace "kenya_counties.shp" or "NAirobi1.shp" with your actual file path.
# Ensure that associated files (.shx, .dbf, .prj) are in the same directory.
kenya_map <- st_read("ken_admin1.shp")
## Reading layer `ken_admin1' from data source 
##   `C:\Users\lenovo\Desktop\Rpubs\Rinla\ken_admin1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 47 features and 21 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 33.91029 ymin: -4.67688 xmax: 41.90602 ymax: 5.414124
## Geodetic CRS:  WGS 84
# Verify the spatial geometry by plotting a basic outline
plot(st_geometry(kenya_map), col = "grey", main = "Imported Study Area Geometry")

# Create a unique integer ID required by R-INLA for spatial index mapping
kenya_map$region_id <- 1:nrow(kenya_map)

Handle Disease Data (Observed vs. Expected Cases)

In real-world disease mapping, we look at two main numbers for every area:
- Observed Cases: The actual number of sickness counts recorded.

- Expected Cases: The baseline number of cases we expect to see purely based on the local population size.
For example, a crowded city naturally yields more cases than a quiet village, even if the underlying risk is the same.
For learning purposes, we will simulate this data here using a Poisson distribution (the standard statistical distribution for counting rare events like diseases).

set.seed(42)
n_areas <- nrow(kenya_map)
kenya_map$expected <- runif(n_areas, 10, 50) # Simulated baseline population expectations

# Generate simulated true relative risk and sample observed cases using a Poisson distribution

relative_risk <- exp(rnorm(n_areas, mean = 0, sd = 0.2)) 
kenya_map$observed <- rpois(n_areas, lambda = kenya_map$expected * relative_risk)

Create a Spatial Neighborhood Graph

Diseases don’t care about administrative borders; an outbreak in one county easily spills over into the next. Therefore, spatial models assume that neighboring areas (contiguous) behave similarly.
Here, we use poly2nb to scan our map and build a “neighborhood directory” listing exactly which areas share a border. We then convert this directory into a .graph file format that R-INLA can easily read.

# Step A: Identify which regions are contiguous neighbors
nb <- poly2nb(kenya_map)

# Step B: Save this neighborhood structure into an INLA-readable graph file
nb2INLA("kenya.graph", nb)

Specify and Fit the BYM Disease Mapping Model

Now we build our statistical model using the famous Besag-York-Mollié (BYM) formula. The BYM model is brilliant because it splits your map’s hidden patterns into two helpful components:
- Spatially Structured Variation: Patterns that look smooth across borders (clumped risks shared with neighbors).
- Unstructured Variation: Pure local noise, quirks, or random luck unique to that specific area.
We use an offset(log(expected)) to ensure the model automatically adjusts for different population sizes across the map regions.

# Define our formula: observed cases depend on a baseline constant (1) and a spatial BYM effect (f)
formula <- observed ~ 1 + f(region_id, model = "bym", graph = "kenya.graph")

# Run the Bayesian engine
model_result <- inla(
  formula,
  family = "poisson",
  data = as.data.frame(kenya_map),
  offset = log(expected),
  control.predictor = list(compute = TRUE),# Tells INLA to calculate predictions for us
  control.compute = list(dic = TRUE, waic = TRUE) # Computes diagnostic values to assess model fit
)

Extract Fitted Relative Risks and Generate Risk Map

The model has finished running! Our goal now is to extract the Relative Risk (RR).
- An RR of 1.0 means the area has an completely average risk.
- An RR of 1.5 means the area has a 50% higher risk of disease than expected.
- An RR of 0.8 means the risk is 20% lower.
We extract these average predictions (mean) and use ggplot2 to draw a clear map. We use a vibrant color scale (YlOrRd - Yellow to Orange to Red) so that high-risk hot-spots immediately pop out to policy-makers.

# Extract the mean posterior values of the fitted relative risks
kenya_map$fitted_risk <- model_result$summary.fitted.values[, "mean"]

# Produce an automatic textual interpretation of the calculated risk values
report(kenya_map$fitted_risk)
## x: n = 47, Mean = 32.43, SD = 12.15, Median = 32.90, MAD = 9.00, range: [9.76,
## 59.45], Skewness = 8.45e-03, Kurtosis = -0.18, 0 missing
# Plot the final spatial risk map
ggplot(kenya_map) +
  geom_sf(aes(fill = fitted_risk)) +
  # Using scale_fill_distiller with direction = 1 ensures light yellow is low risk and dark red is high risk
  scale_fill_distiller(palette = "YlOrRd", direction = 1) +
  theme_minimal() +
  labs(
    title = "Estimated Disease Relative Risk Mapping",
    subtitle = "Fitted via R-INLA spatial BYM model using local shapefile data",
    fill = "Relative Risk"
  )

# Save the plot as a PNG file in your working directory
ggsave(
  filename = "kenya_disease_risk_map.png", 
  plot = last_plot(), # Automatically grabs the last ggplot you generated
  width = 8,          # Width in inches
  height = 5,         # Height in inches
  dpi = 300           # Sharp, high-resolution quality
)