Teaching Lab: Spatio-temporal Disease Analysis with R (Simulated Data)
PGS.TS. Trần Ngọc Đăng (template)
16-12-2025
Mục tiêu của lab
- Tạo dữ liệu giả lập (bản đồ + ca bệnh theo tuần)
- Đọc – kiểm tra – hợp nhất dữ liệu một cách “chắc tay” (có kiểm lỗi)
- Vẽ bản đồ theo thời gian (small multiples)
- Fit mô hình không gian–thời gian với INLA: Poisson + BYM2 (space) + RW1 (time)
- Trích xuất và vẽ Relative Risk (RR)
0. Chuẩn bị môi trường
0.1 Cài và nạp gói
pkgs <- c("sf","dplyr","tidyr","ggplot2","readr","lubridate","spdep")
to_install <- pkgs[!pkgs %in% rownames(installed.packages())]
if(length(to_install)) install.packages(to_install)
# INLA (repo riêng)
if(!requireNamespace("INLA", quietly = TRUE)) {
install.packages(
"INLA",
repos = c(getOption("repos"), INLA = "https://inla.r-inla-download.org/R/stable"),
dep = TRUE
)
}
library(sf)
library(dplyr)
library(tidyr)
library(ggplot2)
library(readr)
library(lubridate)
library(spdep)
library(INLA)0.2 Tạo thư mục dữ liệu
dir.create("data", showWarnings = FALSE)
list.files("data")## [1] "dengue_cases_weekly.csv" "hanoi.graph"
## [3] "hanoi_districts.dbf" "hanoi_districts.gpkg"
## [5] "hanoi_districts.prj" "hanoi_districts.shp"
## [7] "hanoi_districts.shx"1. Tạo dữ liệu giả lập
1.1 Tạo bản đồ “quận/huyện” giả lập
Lưu ý teaching point: Shapefile (.shp) bị giới hạn tên biến <= 10 ký tự (DBF) ⇒ dễ bị cắt tên cột.
Vì vậy trong lab này: - Vẫn ghi
.shpđể bạn thấy lỗi thật ngoài đời, - Đồng thời ghi thêm.gpkgđể có lựa chọn “sạch” cho bài thật.
bbox <- st_bbox(c(xmin = 0, ymin = 0, xmax = 12, ymax = 8))
grid <- st_make_grid(st_as_sfc(bbox), n = c(6,4), what = "polygons") # 24 polygons
hanoi_map <- st_sf(
district_id = 1:length(grid),
district_name = c(
"Ba Dinh","Hoan Kiem","Tay Ho","Long Bien","Cau Giay","Dong Da",
"Hai Ba Trung","Hoang Mai","Thanh Xuan","Ha Dong","Bac Tu Liem","Nam Tu Liem",
"Gia Lam","Dong Anh","Soc Son","Me Linh","Thanh Tri","Thuong Tin",
"Phu Xuyen","Quoc Oai","Thach That","Dan Phuong","Hoai Duc","Son Tay"
),
geometry = st_sfc(grid)
)
# Ghi SHP (có thể bị cắt tên cột)
st_write(hanoi_map, "data/hanoi_districts.shp", delete_dsn = TRUE, quiet = TRUE)
# Ghi GPKG (không bị cắt tên cột)
st_write(hanoi_map, "data/hanoi_districts.gpkg", delete_dsn = TRUE, quiet = TRUE)
hanoi_map1.2 Tạo ca bệnh theo tuần (CSV)
weeks <- 1:52
year <- 2023
district_cov <- hanoi_map |>
st_drop_geometry() |>
mutate(
pop_density = exp(rnorm(n(), log(12000), 0.5)),
log_pop_density = log(pop_density),
population = round(runif(n(), 120000, 420000))
)
# "hotspot" giả lập (nguy cơ cao hơn)
hot <- c("Hoan Kiem","Ba Dinh","Dong Da","Hai Ba Trung","Cau Giay","Thanh Xuan")
district_cov <- district_cov |>
mutate(spatial_risk = ifelse(district_name %in% hot, 0.5, 0))
# mùa vụ: đỉnh quanh tuần 30
season <- tibble(
week = weeks,
season_effect = 1.2 * exp(-0.5 * ((week - 30)/7)^2)
)
analysis_data <- expand_grid(
district_name = district_cov$district_name,
year = year,
week = weeks
) |>
left_join(district_cov, by = "district_name") |>
left_join(season, by = "week") |>
mutate(
base_rate = 8/100000, # baseline weekly incidence (simulated)
expected_cases = population * base_rate,
linpred = -0.2 + 0.25 * scale(log_pop_density)[,1] + spatial_risk + season_effect,
rr_true = exp(linpred),
mu = expected_cases * rr_true,
cases = rpois(n(), lambda = pmax(mu, 0.1))
) |>
select(district_name, year, week, cases, population, expected_cases, log_pop_density)
write_csv(analysis_data, "data/dengue_cases_weekly.csv")
glimpse(analysis_data)## Rows: 1,248
## Columns: 7
## $ district_name <chr> "Ba Dinh", "Ba Dinh", "Ba Dinh", "Ba Dinh", "Ba Dinh",…
## $ year <dbl> 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2023, …
## $ week <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…
## $ cases <int> 21, 25, 26, 30, 17, 29, 32, 24, 25, 39, 26, 17, 24, 26…
## $ population <dbl> 273615, 273615, 273615, 273615, 273615, 273615, 273615…
## $ expected_cases <dbl> 21.8892, 21.8892, 21.8892, 21.8892, 21.8892, 21.8892, …
## $ log_pop_density <dbl> 9.292431, 9.292431, 9.292431, 9.292431, 9.292431, 9.29…2. Đọc dữ liệu + kiểm tra chất lượng
2.1 Helper: chuẩn hoá tên cột khi đọc SHP bị cắt
Teaching point: Nếu đọc SHP mà không thấy
district_namevì bị cắt thànhdstrct_n(hoặc tương tự), ta tự động rename.
normalize_hanoi_map_cols <- function(x) {
nm <- names(x)
# Case 1: đúng chuẩn rồi
if(all(c("district_id","district_name") %in% nm)) return(x)
# Case 2: bị cắt như bạn gặp: dstrct_d, dstrct_n
if(all(c("dstrct_d","dstrct_n") %in% nm)) {
return(dplyr::rename(x, district_id = dstrct_d, district_name = dstrct_n))
}
# Case 3: thử đoán theo pattern (fallback)
id_guess <- nm[grepl("^dstrct", nm) & grepl("d$", nm)]
name_guess <- nm[grepl("^dstrct", nm) & grepl("n$", nm)]
if(length(id_guess) == 1 && length(name_guess) == 1) {
return(dplyr::rename(x, district_id = all_of(id_guess), district_name = all_of(name_guess)))
}
stop(
"Không tìm thấy cột district_id/district_name (hoặc phiên bản bị cắt). ",
"Hãy chạy names(hanoi_map) để xem tên cột thực tế."
)
}2.2 Đọc map + CSV
# Bạn có thể đổi sang .gpkg nếu muốn "sạch" (không bị cắt tên cột)
hanoi_map_raw <- st_read("data/hanoi_districts.shp", quiet = TRUE)
# hanoi_map_raw <- st_read("data/hanoi_districts.gpkg", quiet = TRUE)
dengue_data <- read_csv("data/dengue_cases_weekly.csv", show_col_types = FALSE)
cat("Map columns:\n"); print(names(hanoi_map_raw))## Map columns:## [1] "dstrct_d" "dstrct_n" "geometry"cat("\nCSV columns:\n"); print(names(dengue_data))##
## CSV columns:## [1] "district_name" "year" "week" "cases"
## [5] "population" "expected_cases" "log_pop_density"2.3 Chuẩn hoá tên cột & hợp nhất
hanoi_map <- normalize_hanoi_map_cols(hanoi_map_raw)
# Check khóa join
stopifnot("district_name" %in% names(hanoi_map))
stopifnot("district_name" %in% names(dengue_data))
analysis_data_sf <- hanoi_map |>
left_join(dengue_data, by = "district_name")
# Quick sanity checks
cat("Số quận/huyện:", nrow(hanoi_map), "\n")## Số quận/huyện: 24cat("Số dòng sau join (district-week):", nrow(analysis_data_sf), "\n")## Số dòng sau join (district-week): 1248cat("Tỷ lệ NA cases:", mean(is.na(analysis_data_sf$cases)), "\n")## Tỷ lệ NA cases: 0analysis_data_sf |>
select(district_name, year, week, cases, population) |>
head()3. Khám phá dữ liệu (EDA) theo kiểu giảng dạy
3.1 Phân bố ca bệnh
analysis_data_sf |>
st_drop_geometry() |>
ggplot(aes(x = cases)) +
geom_histogram(bins = 30) +
labs(title = "Distribution of weekly cases (simulated)", x = "Cases", y = "Count") +
theme_minimal()
3.2 Xu hướng theo tuần (toàn thành phố)
analysis_data_sf |>
st_drop_geometry() |>
group_by(week) |>
summarise(total_cases = sum(cases), .groups = "drop") |>
ggplot(aes(x = week, y = total_cases)) +
geom_line() +
labs(title = "Total cases by week (simulated)", x = "Week", y = "Total cases") +
theme_minimal()
4. Bản đồ theo thời gian (Small multiples)
Teaching tip: Chọn vài tuần đại diện để nhìn mùa vụ / cụm nóng.
library(sf)
library(dplyr)
library(ggplot2)
# 1) Đảm bảo analysis_data_sf là sf hợp lệ
stopifnot(inherits(analysis_data_sf, "sf"))
# 2) Nếu CRS đang NA hoặc "lạ", gán CRS phẳng (Web Mercator) để ggplot không đòi transform
if (is.na(st_crs(analysis_data_sf))) {
st_crs(analysis_data_sf) <- 3857
} else {
# đôi khi CRS có nhưng GDAL không hiểu -> gán lại 3857 cho chắc
st_crs(analysis_data_sf) <- 3857
}
weeks_show <- c(18, 22, 26, 30, 34, 38, 42, 46)
analysis_data_sf |>
filter(week %in% weeks_show) |>
ggplot() +
geom_sf(aes(fill = cases), color = "grey30", linewidth = 0.2) +
facet_wrap(~week, ncol = 4) +
labs(title = "Weekly cases by district (simulated)", fill = "Cases") +
theme_minimal()
5. Mô hình không gian–thời gian với INLA
5.1 Tạo adjacency graph (láng giềng)
# =========================
# 5.1 Tạo adjacency graph (láng giềng) cho INLA
# (robust: sf::st_touches -> sparse adjacency -> INLA graph)
# =========================
library(sf)
library(dplyr)
library(Matrix)
library(INLA)
map_unique <- hanoi_map |>
arrange(district_id) |>
mutate(district_id = as.integer(district_id)) # chắc chắn là int
# 1) Tạo danh sách láng giềng theo tiếp giáp biên (queen=FALSE là rook; TRUE ~ queen)
# Ở đây dùng touches: chung biên/điểm đều tính là láng giềng (queen-like)
nb_list <- sf::st_touches(map_unique)
# 2) Chuyển sang ma trận kề (sparse)
n <- nrow(map_unique)
i <- rep(seq_len(n), lengths(nb_list))
j <- unlist(nb_list)
A <- Matrix::sparseMatrix(
i = i, j = j,
x = 1,
dims = c(n, n)
)
# đối xứng + bỏ self-loop
A <- (A + t(A)) > 0
diag(A) <- 0
A <- Matrix::Matrix(A * 1, sparse = TRUE)
# 3) Ghi graph theo định dạng INLA
# inla.write.graph() nhận "adjacency matrix" (sparse) ổn định nhất
INLA::inla.write.graph(A, file = "data/hanoi.graph")
# 4) Đọc lại graph để kiểm tra
g <- INLA::inla.read.graph("data/hanoi.graph")
cat("Số node trong graph:", g$n, "\n")## Số node trong graph: 245.2 Chuẩn bị data cho INLA
analysis_inla <- analysis_data_sf |>
st_drop_geometry() |>
mutate(
district_idx = district_id,
week_idx = week
)
# Check
stopifnot(!any(is.na(analysis_inla$expected_cases)))
stopifnot(!any(is.na(analysis_inla$district_idx)))
stopifnot(!any(is.na(analysis_inla$week_idx)))
glimpse(analysis_inla)## Rows: 1,248
## Columns: 10
## $ district_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ district_name <chr> "Ba Dinh", "Ba Dinh", "Ba Dinh", "Ba Dinh", "Ba Dinh",…
## $ year <dbl> 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2023, 2023, …
## $ week <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…
## $ cases <dbl> 21, 25, 26, 30, 17, 29, 32, 24, 25, 39, 26, 17, 24, 26…
## $ population <dbl> 273615, 273615, 273615, 273615, 273615, 273615, 273615…
## $ expected_cases <dbl> 21.8892, 21.8892, 21.8892, 21.8892, 21.8892, 21.8892, …
## $ log_pop_density <dbl> 9.292431, 9.292431, 9.292431, 9.292431, 9.292431, 9.29…
## $ district_idx <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ week_idx <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,…5.3 Fit mô hình: Poisson + BYM2 + RW1
formula <- cases ~ 1 +
log_pop_density +
f(district_idx, model = "bym2", graph = g) +
f(week_idx, model = "rw1")
fit <- inla(
formula,
family = "poisson",
data = analysis_inla,
E = expected_cases,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)
fit$waic$waic## [1] 7789.551fit$dic$dic## [1] 7789.0315.4 Diễn giải nhanh output (teaching checklist)
# fixed effects
fit$summary.fixed# random effects (space/time)
names(fit$summary.random)## [1] "district_idx" "week_idx"6. Trích xuất & vẽ bản đồ RR
6.1 RR không gian (xấp xỉ từ random effect BYM2)
library(sf)
library(dplyr)
library(ggplot2)
# 0) Ép CRS hợp lệ cho map_unique (nếu NA hoặc CRS lỗi)
# Vì đây là lưới giả lập (tọa độ phẳng), ta gán CRS phẳng đơn giản.
if (is.na(st_crs(map_unique))) {
st_crs(map_unique) <- 3857
} else {
# nhiều máy CRS có nhưng GDAL không hiểu -> gán lại cho chắc
st_crs(map_unique) <- 3857
}
# 1) Lấy random effect theo district (BYM2)
# (giữ code của bạn, nhưng đảm bảo đúng kiểu)
re_sp <- fit$summary.random$district_idx |>
transmute(district_id = as.integer(ID), re_mean = mean)
# 2) Join + tính RR
rr_map <- map_unique |>
left_join(re_sp, by = "district_id") |>
mutate(RR = exp(re_mean))
# 3) Ép CRS lại lần nữa cho rr_map (để chắc chắn)
st_crs(rr_map) <- 3857
# 4) Vẽ
ggplot(rr_map) +
geom_sf(aes(fill = RR), color = "grey30", linewidth = 0.2) +
labs(title = "Spatial Relative Risk (approx.)", fill = "RR") +
theme_minimal()
6.2 Ví dụ “prediction” cho 1 tuần mục tiêu
target_week <- 40
pred_rows <- which(analysis_inla$week_idx == target_week)
pred_mean <- fit$summary.fitted.values$mean[pred_rows]
pred_map <- map_unique |>
left_join(
analysis_inla |> filter(week_idx == target_week) |> select(district_name, expected_cases),
by = "district_name"
) |>
mutate(
mu_pred = pred_mean,
RR_pred = mu_pred / expected_cases
)
ggplot(pred_map) +
geom_sf(aes(fill = RR_pred), color = "grey30", linewidth = 0.2) +
labs(title = paste0("Predicted RR (week ", target_week, ")"), fill = "RR") +
theme_minimal()
7. Bài tập (Exercises)
- Thay
rw1→rw2. So sánh WAIC và hình dạng xu hướng thời gian.
- Thêm covariate theo tuần (mưa/nhiệt độ giả lập) và xem hệ số.
- Thử mô hình Negative Binomial nếu bạn cố tình làm dữ liệu over-dispersion.
8. Troubleshooting nhanh
Lỗi join không thấy district_name
- Nếu dùng
.shp, tên cột có thể bị cắt (VD:dstrct_n).
- Trong lab này đã có
normalize_hanoi_map_cols()để tự sửa.
- Cách “sạch”: dùng
data/hanoi_districts.gpkg.