Analisis Epidemiologi Tuberkulosis di Indonesia

Insidens, Pemetaan Risiko, Hotspot Spasial, dan Determinan antarprovinsi, 2022-2024

library(tidyverse)
library(readxl)
library(sf)
library(spdep)
library(ggrepel)
library(knitr)
library(kableExtra)
library(MASS)
library(car)
library(scales)
library(corrplot)
library(CARBayes)    # model spasial CAR Poisson

Ringkasan Proyek

Studi ekologi tingkat provinsi terhadap pola spasial dan determinan kejadian tuberkulosis di Indonesia periode 2022-2024. Analisis mencakup pemetaan tematik kejadian (incidence rate & SIR), evaluasi autokorelasi spasial (Moran’s I global, LISA, Getis-Ord Gi*), dan identifikasi determinan melalui regresi data cacah (Negative Binomial dan model CAR Poisson spasial) terhadap faktor sosioekonomi, demografi, serta sistem kesehatan.

1 Abstrak

Latar belakang. Tuberkulosis (TB) masih menjadi tantangan utama kesehatan masyarakat, dan Indonesia secara konsisten tergolong negara dengan beban TB tertinggi kedua di dunia. Karena itu, faktor kejadian TB menunjukkan variasi antarwilayah yang dipengaruhi kondisi demografi, sosial ekonomi, lingkungan, dan akses pelayanan kesehatan, pemahaman atas pola spasial dan determinannya penting untuk mendukung pengendalian yang tepat sasaran.

Metode. Penelitian ini menganalisis data agregat tingkat provinsi (34 provinsi yang memiliki data lengkap dan konsisten dengan batas wilayah) selama 2022-2024, bersumber dari publikasi BPS dan Profil Kesehatan Indonesia. Kejadian TB dideskripsikan melalui incidence rate (IR) dan standardized incidence ratio (SIR), divisualisasikan dalam peta tematik. Autokorelasi spasial global dievaluasi dengan indeks Moran’s I, klaster lokal dengan Local Indicators of Spatial Association (LISA), dan intensitas pengelompokan dengan statistik Getis-Ord Gi*. Determinan jumlah kasus dimodelkan menggunakan regresi data cacah dengan offset jumlah penduduk terhadap lima prediktor: indeks pembangunan manusia (IPM), persentase penduduk miskin, akses sanitasi layak, kepadatan penduduk, dan rasio Puskesmas per kecamatan. Pemilihan model ditentukan secara empiris: uji overdispersi mengarahkan penggunaan regresi Negative Binomial, dan uji autokorelasi spasial pada residual menentukan kebutuhan model spasial. Bila residual tidak berautokorelasi, Negative Binomial dinilai memadai; bila residual berautokorelasi spasial, digunakan model Conditional Autoregressive (CAR) Poisson dengan efek acak spasial.

Hasil. Total kasus TB ternotifikasi meningkat dari 724.309 (2022) menjadi 842.383 (2024), dengan rata-rata incidence rate (IR) antarprovinsi naik dari 245,3 menjadi 286,3 per 100.000 penduduk. Berdasarkan SIR 2022, sdebanyak sembilan provinsi memiliki risiko di atas rata-rata nasional, tertinggi di DKI Jakarta (SIR=1,83), Gorontalo (1,51), dan Jawa Barat (1,44). Indeks Moran’s I positif dan signifikan sepanjang periode (I=0,233–0,273; p<0,05), menunjukkan autokorelasi spasial yang konsisten — provinsi ber-IR tinggi cenderung mengelompok secara geografis. Pada model Negative Binomial 2022, rasio Puskesmas per kecamatan berasosiasi signifikan dengan penurunan insidens (IRR=0,644; 95% CI: 0,456–0,918) dan kepadatan penduduk dengan kenaikan insidens (p=0,002), sementara IPM, kemiskinan, dan sanitasi tidak signifikan. Uji autokorelasi residual mengarahkan penggunaan model CAR Poisson untuk 2024 (residual berautokorelasi, p=0,044), yang hasilnya konsisten dengan model Negative Binomial.

Kesimpulan. Insidens TB di Indonesia menunjukkan tren meningkat dan pengelompokan spasial yang menetap antarprovinsi sepanjang 2022–2024, dengan beban tertinggi terkonsentrasi di provinsi padat penduduk seperti DKI Jakarta dan sekitarnya. Rasio Puskesmas per kecamatan yang lebih tinggi berasosiasi dengan insidens lebih rendah, mengindikasikan peran penting ketersediaan fasilitas kesehatan primer. Temuan ini mendukung strategi pengendalian TB berbasis wilayah yang menargetkan provinsi hotspot dan penguatan layanan kesehatan primer.

Kata kunci: tuberkulosis; epidemiologi spasial; incidence rate; Indonesia; autokorelasi spasial; Negative Binomial; model CAR.

2 Pendahuluan

Tuberkulosis (TB) merupakan salah satu penyakit menular yang masih menjadi tantangan utama kesehatan masyarakat di dunia. Penyakit yang disebabkan oleh Mycobacterium tuberculosis ini ditularkan melalui udara dan menjadi salah satu penyebab utama kematian akibat penyakit infeksi. Indonesia termasuk negara dengan beban TB tertinggi kedua di dunia dan secara konsisten berada dalam kelompok negara prioritas pengendalian TB menurut Organisasi Kesehatan Dunia (WHO). Tingginya jumlah kasus dan kematian akibat TB menunjukkan bahwa upaya pengendalian penyakit ini masih memerlukan perhatian serius, terutama dalam mengidentifikasi wilayah berisiko tinggi serta faktor-faktor yang berkontribusi terhadap tingginya kejadian penyakit.

Dalam epidemiologi, distribusi penyakit dipelajari berdasarkan dimensi orang (person), tempat (place), dan waktu (time). Dimensi tempat berperan penting karena kejadian TB sering menunjukkan variasi antarwilayah yang dipengaruhi oleh kondisi demografi, sosial ekonomi, lingkungan, maupun akses terhadap pelayanan kesehatan. Sebagai negara kepulauan dengan karakteristik antarprovinsi yang sangat beragam, Indonesia berpotensi mengalami perbedaan beban TB antarwilayah, sehingga pemahaman atas pola distribusi geografisnya menjadi penting untuk perencanaan program pengendalian yang efektif dan tepat sasaran.

Pendekatan epidemiologi spasial memungkinkan identifikasi pola persebaran penyakit dan hubungan antarwilayah yang berdekatan secara geografis. Salah satu metode yang umum digunakan adalah analisis autokorelasi spasial global menggunakan indeks Moran’s I; nilai yang positif dan signifikan menunjukkan wilayah berdekatan cenderung memiliki tingkat kejadian serupa sehingga membentuk pola pengelompokan (spatial clustering). Selain itu, karena jumlah kasus TB merupakan data cacah (count data) yang sering mengalami overdispersi (varians melebihi rata-rata), model Negative Binomial lebih sesuai dibandingkan Poisson standar untuk mengestimasi hubungan antara faktor sosial, demografi, dan kesehatan dengan jumlah kasus secara lebih akurat.

Beberapa penelitian sebelumnya menerapkan pemodelan spasiotemporal Bayesian kompleks pada data TB tingkat kabupaten/kota dengan jumlah unit spasial besar. Pada tingkat provinsi dengan unit wilayah lebih terbatas, regresi Negative Binomial memberikan alternatif yang lebih sederhana namun tetap memadai untuk mengidentifikasi determinan kejadian TB. Berdasarkan latar belakang tersebut, penelitian ini bertujuan untuk: (1) menggambarkan distribusi spasial kasus TB antarprovinsi di Indonesia; (2) mengevaluasi keberadaan autokorelasi spasial menggunakan indeks Moran’s I; dan (3) mengidentifikasi determinan yang berhubungan dengan jumlah kasus TB menggunakan regresi data cacah.

Pertanyaan penelitian.

1. Bagaimana distribusi spasial kasus tuberkulosis antarprovinsi di Indonesia?
2. Apakah terdapat autokorelasi spasial pada distribusi kasus tuberkulosis antarprovinsi?
3. Faktor-faktor apa yang berhubungan dengan variasi jumlah kasus tuberkulosis pada tingkat provinsi di Indonesia?

3 Data dan Metode

3.1 Sumber dan Unit Data

sumber: Data dikumpulkan dari publikasi BPS dan Profil Kesehatan Indonesia yang diakses pada tanggal 12 Juni 2026. Unit: provinsi (34), 3 tahun.

Variabel	Definisi	Satuan
Kasus	Jumlah kasus TB ternotifikasi (L+P)	kasus
Penduduk	Jumlah penduduk	jiwa
IR	Insidens rate = Kasus/Penduduk x 100.000	per 100.000
IPM	Indeks Pembangunan Manusia	indeks
Miskin	Persentase penduduk miskin	%
Sanitasi	Persentase sanitasi layak	%
Kepadatan	Kepadatan penduduk	per km2
Rasio_Pkm	Rasio Puskesmas per kecamatan	rasio
JKN	Total peserta aktif JKN	jiwa

3.2 Pembacaan Data (dari GitHub)

user   <- "Winaliaagwil"
repo   <- "analisis_tb_indonesia"
branch <- "main"
base   <- sprintf("https://raw.githubusercontent.com/%s/%s/%s/", user, repo, branch)

url_xlsx <- paste0(base, "Data%20TB%202022_2024.xlsx")
url_zip  <- paste0(base, "indonesia_34_provinsi.zip")

download.file(url_xlsx, "DataTB.xlsx", mode = "wb", quiet = TRUE)
download.file(url_zip,  "shp.zip",     mode = "wb", quiet = TRUE)
unzip("shp.zip", exdir = "shp_folder")

rename_map <- c(
  "Wilayah / Provinsi"             = "Provinsi",
  "IPM"                            = "IPM",
  "Persentase Penduduk Miskin"     = "Miskin",
  "Persentase Sanitasi Layak"      = "Sanitasi",
  "Kepadatan penduduk (per km2)"   = "Kepadatan",
  "Rasio Puskesmas per Kecamatan"  = "Rasio_Pkm",
  "Total Peserta Aktif JKN (Jiwa)" = "JKN",
  "IR TB per 100000 Penduduk"      = "IR",
  "Jumlah Kasus TB"                = "Kasus",
  "Jumlah Penduduk"                = "Penduduk"
)

baca_sheet <- function(path, sheet) {
  df <- read_excel(path, sheet = sheet)
  names(df) <- dplyr::recode(names(df), !!!rename_map)
  prov <- str_squish(toupper(df$Provinsi))
  prov[prov == "DI YOGYAKARTA"]        <- "DAERAH ISTIMEWA YOGYAKARTA"
  prov[prov == "KEP. BANGKA BELITUNG"] <- "KEPULAUAN BANGKA BELITUNG"
  df$Provinsi <- prov
  df$Tahun <- as.integer(sheet)
  df %>%
    dplyr::select(Provinsi, Tahun, Kasus, Penduduk, IR,
                  IPM, Miskin, Sanitasi, Kepadatan, Rasio_Pkm, JKN)
}

tb_all <- bind_rows(lapply(c("2022","2023","2024"),
                           function(s) baca_sheet("DataTB.xlsx", s)))
tb22 <- filter(tb_all, Tahun == 2022)
tb23 <- filter(tb_all, Tahun == 2023)
tb24 <- filter(tb_all, Tahun == 2024)

indo <- st_read(list.files("shp_folder", pattern="\\.shp$", full.names=TRUE)[1], quiet=TRUE)
indo$Provinsi <- str_squish(toupper(indo$Provinsi))

sf::sf_use_s2(FALSE)
indo <- sf::st_make_valid(indo)

tidak_cocok <- setdiff(indo$Provinsi, tb22$Provinsi)
if (length(tidak_cocok)>0) message("Provinsi peta tanpa data: ", paste(tidak_cocok, collapse=", "))

# daftar data per tahun untuk dipakai ulang
data_tahun <- list("2022"=tb22, "2023"=tb23, "2024"=tb24)

Tahun 2022

724,309

Total kasus

245.3

Rata-rata IR

34

Provinsi

Dki Jakarta

IR tertinggi

Tahun 2023

810,882

Total kasus

282.7

Rata-rata IR

34

Provinsi

Papua

IR tertinggi

Tahun 2024

842,383

Total kasus

286.3

Rata-rata IR

34

Provinsi

Papua

IR tertinggi

3.3 Metode Statistik

Incidence Rate (IR) dan Standardized Incidence Ratio (SIR): \[ IR=\frac{O_i}{P_i}\times 100{,}000, \qquad E_i=\bar r\times P_i, \qquad SIR_i=\frac{O_i}{E_i} \] dengan \(\bar r\) rate nasional. SIR > 1 = risiko di atas rata-rata nasional.

Pola & hotspot spasial. Moran’s I global (autokorelasi), LISA (klaster lokal), dan Getis-Ord Gi* (hotspot/coldspot). Pembobot k-nearest neighbours (k=4) dipakai karena Indonesia kepulauan.

Determinan. Karena outcome berupa cacahan kasus dengan populasi sebagai offset, dipakai regresi count. Overdispersi diuji (Poisson vs Negative Binomial), autokorelasi spasial residual diuji, dan model spasial (CAR Poisson) disertakan sebagai analisis sensitivitas.

# ============================================================
# Fungsi-fungsi pembuat peta agar bisa dipakai ulang per tahun
# ============================================================

# bobot spasial (sama untuk semua tahun karena geometri tetap)
coords  <- suppressWarnings(st_coordinates(st_centroid(indo)))
k       <- 4
nb_k    <- knn2nb(knearneigh(coords, k=k))
lw      <- nb2listw(nb_k, style="W")
lw_self <- nb2listw(include.self(nb_k), style="W")

# --- Peta IR ---
peta_ir <- function(df, th){
  peta <- indo %>% left_join(df, by="Provinsi")
  top5 <- peta %>% filter(!is.na(IR)) %>% slice_max(IR, n=5)
  co <- as.data.frame(st_coordinates(suppressWarnings(st_centroid(top5))))
  co$label <- paste0(str_to_title(top5$Provinsi), "\n", round(top5$IR,1))
  ggplot(peta) +
    geom_sf(aes(fill=IR), color="white", linewidth=0.15) +
    scale_fill_viridis_c(option="plasma", direction=-1, name="IR TB\n(per 100.000)", na.value="grey85") +
    geom_point(data=co, aes(X,Y), size=1.3) +
    ggrepel::geom_label_repel(data=co, aes(X,Y,label=label), size=3, fontface="bold",
      alpha=0.9, box.padding=0.8, min.segment.length=0, segment.size=0.5,
      arrow=arrow(length=unit(0.015,"npc")), seed=42) +
    labs(title=paste0("Angka Insidensi TB per Provinsi, ", th),
         caption=paste0("Sumber: Profil Kesehatan Indonesia ", th)) +
    theme_minimal(base_size=12) +
    theme(axis.text=element_blank(), axis.title=element_blank(),
          axis.ticks=element_blank(), panel.grid=element_blank())
}

# --- Hitung SIR ---
hitung_sir <- function(df){
  rn <- sum(df$Kasus)/sum(df$Penduduk)
  df %>% mutate(Expected=rn*Penduduk, SIR=Kasus/Expected)
}

# --- Peta SIR ---
peta_sir_fun <- function(df, th){
  sir <- hitung_sir(df)
  peta <- indo %>% left_join(sir, by="Provinsi")
  ggplot(peta) +
    geom_sf(aes(fill=SIR), color="white", linewidth=0.15) +
    scale_fill_gradient2(low="#2166ac", mid="white", high="#b2182b", midpoint=1,
                         name="SIR", na.value="grey85") +
    labs(title=paste0("Standardized Incidence Ratio (SIR) TB, ", th),
         caption="SIR > 1 (merah) = risiko di atas rata-rata nasional") +
    theme_minimal(base_size=12) +
    theme(axis.text=element_blank(), axis.title=element_blank(),
          axis.ticks=element_blank(), panel.grid=element_blank())
}

# --- Tabel SIR ---
tabel_sir <- function(df, th){
  hitung_sir(df) %>% arrange(desc(SIR)) %>%
    transmute(Provinsi=str_to_title(Provinsi), Observed=Kasus,
              Expected=round(Expected,0), SIR=round(SIR,2)) %>%
    head(10) %>%
    kable(caption=paste0("Sepuluh provinsi SIR tertinggi (", th, ").")) %>%
    kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))
}

# --- Peta LISA ---
peta_lisa <- function(df, th){
  peta <- indo %>% left_join(df, by="Provinsi")
  x <- peta$IR; xs <- as.numeric(scale(x)); lagxs <- lag.listw(lw, xs)
  locm <- localmoran(x, lw, zero.policy=TRUE); sig <- locm[,5] < 0.05
  peta$klaster <- "Tidak signifikan"
  peta$klaster[sig & xs>0 & lagxs>0] <- "High-High"
  peta$klaster[sig & xs<0 & lagxs<0] <- "Low-Low"
  peta$klaster[sig & xs>0 & lagxs<0] <- "High-Low"
  peta$klaster[sig & xs<0 & lagxs>0] <- "Low-High"
  peta$klaster <- factor(peta$klaster,
    levels=c("High-High","Low-Low","High-Low","Low-High","Tidak signifikan"))
  ggplot(peta) +
    geom_sf(aes(fill=klaster), color="white", linewidth=0.15) +
    scale_fill_manual(values=c("High-High"="#b2182b","Low-Low"="#2166ac",
      "High-Low"="#ef8a62","Low-High"="#67a9cf","Tidak signifikan"="grey88"),
      name="Klaster LISA", drop=FALSE) +
    labs(title=paste0("Klaster Spasial LISA Insidens TB, ", th),
         caption="High-High = hotspot; Low-Low = coldspot") +
    theme_minimal(base_size=12) +
    theme(axis.text=element_blank(), axis.title=element_blank(),
          axis.ticks=element_blank(), panel.grid=element_blank())
}

# --- Peta Hotspot Getis-Ord Gi* ---
peta_getis <- function(df, th){
  peta <- indo %>% left_join(df, by="Provinsi")
  gi <- as.numeric(localG(peta$IR, lw_self, zero.policy=TRUE))
  peta$hot <- cut(gi, breaks=c(-Inf,-2.58,-1.96,-1.65,1.65,1.96,2.58,Inf),
    labels=c("Coldspot 99%","Coldspot 95%","Coldspot 90%","Tidak signifikan",
             "Hotspot 90%","Hotspot 95%","Hotspot 99%"))
  warna_g <- c("Coldspot 99%"="#2166ac","Coldspot 95%"="#67a9cf","Coldspot 90%"="#d1e5f0",
    "Tidak signifikan"="grey88","Hotspot 90%"="#fddbc7","Hotspot 95%"="#ef8a62","Hotspot 99%"="#b2182b")
  ggplot(peta) +
    geom_sf(aes(fill=hot), color="white", linewidth=0.15) +
    scale_fill_manual(values=warna_g, name="Getis-Ord Gi*", drop=FALSE) +
    labs(title=paste0("Hotspot Insidens TB (Getis-Ord Gi*), ", th),
         caption="Merah = hotspot (kumpulan IR tinggi); biru = coldspot") +
    theme_minimal(base_size=12) +
    theme(axis.text=element_blank(), axis.title=element_blank(),
          axis.ticks=element_blank(), panel.grid=element_blank())
}

4 Hasil Analisis

4.1 Analisis Deskriptif

4.1.1 Tren Nasional

tb_all %>% group_by(Tahun) %>% summarise(Total=sum(Kasus)) %>%
  ggplot(aes(Tahun, Total)) +
  geom_col(fill="#1a6985", width=0.6) +
  geom_text(aes(label=comma(Total)), vjust=-0.5, size=3.5) +
  scale_x_continuous(breaks=2022:2024) +
  scale_y_continuous(labels=comma, expand=expansion(mult=c(0,0.12))) +
  labs(x=NULL, y="Total kasus TB") + theme_minimal(base_size=12)

Total kasus TB nasional, 2022-2024.

Interpretasi. Total kasus TB nasional bergerak dari 724,309 (2022) menjadi 842,383 (2024), menunjukkan tren yang meningkat selama periode pengamatan. Perubahan ini dapat mencerminkan dinamika penemuan kasus (notifikasi) maupun beban penyakit yang sesungguhnya.

4.1.2 Provinsi IR Tertinggi

plot_top10 <- function(df, th){
  df %>% slice_max(IR, n=10) %>%
    ggplot(aes(reorder(str_to_title(Provinsi), IR), IR)) +
    geom_col(fill="#e07b39") +
    geom_text(aes(label=round(IR,0)), hjust=-0.2, size=3) +
    coord_flip() + scale_y_continuous(expand=expansion(mult=c(0,0.12))) +
    labs(x=NULL, y=paste0("IR TB per 100.000 (", th, ")")) + theme_minimal(base_size=12)
}

4.1.2.1 2022

plot_top10(tb22, 2022)

Sepuluh provinsi IR tertinggi (2022).

4.1.2.2 2023

plot_top10(tb23, 2023)

Sepuluh provinsi IR tertinggi (2023).

4.1.2.3 2024

plot_top10(tb24, 2024)

Sepuluh provinsi IR tertinggi (2024).

4.1.3 Ringkasan Statistik

tb_all %>% group_by(Tahun) %>%
  summarise(Total_kasus=sum(Kasus), IR_min=round(min(IR),1),
            IR_median=round(median(IR),1), IR_mean=round(mean(IR),1),
            IR_max=round(max(IR),1)) %>%
  kable(caption="Tabel 1. Ringkasan kasus & IR per tahun.") %>%
  kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))

Tabel 1. Ringkasan kasus & IR per tahun.

Tahun	Total_kasus	IR_min	IR_median	IR_mean	IR_max
2022	724309	105.7	231.8	245.3	480.2
2023	810882	120.4	264.2	282.7	570.2
2024	842383	121.9	265.4	286.3	615.7

Interpretasi. Rata-rata IR antarprovinsi berkisar 245.3-286.3 per 100.000 penduduk sepanjang 2022-2024. Rentang nilai yang lebar (minimum sekitar 105.7 hingga maksimum 615.7) mengindikasikan disparitas insidens yang besar antarprovinsi, yang menjadi dasar perlunya analisis spasial untuk mengetahui apakah provinsi ber-IR tinggi cenderung mengelompok secara geografis.

4.1.4 Ringkasan Variabel per Tahun

Statistik deskriptif (min, median, rata-rata, maks) untuk seluruh variabel pada tiap tahun.

ringkas_var <- function(df, th){
  df %>%
    dplyr::select(Kasus, IR, IPM, Miskin, Sanitasi, Kepadatan, Rasio_Pkm) %>%
    pivot_longer(everything(), names_to="Variabel", values_to="Nilai") %>%
    group_by(Variabel) %>%
    summarise(Min=round(min(Nilai),1), Median=round(median(Nilai),1),
              Mean=round(mean(Nilai),1), Maks=round(max(Nilai),1), .groups="drop") %>%
    kable(caption=paste0("Ringkasan variabel, ", th, ".")) %>%
    kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))
}

4.1.4.1 2022

ringkas_var(tb22, 2022)

Ringkasan variabel, 2022.

Variabel	Min	Median	Mean	Maks
IPM	61.4	72.2	72.0	81.7
IR	105.7	231.8	245.3	480.2
Kasus	1738.0	8793.0	21303.2	184406.0
Kepadatan	10.0	100.5	772.9	17031.0
Miskin	4.5	8.5	10.3	26.8
Rasio_Pkm	0.7	1.4	1.6	7.2
Sanitasi	40.3	81.7	81.0	96.2

4.1.4.2 2023

ringkas_var(tb23, 2023)

Ringkasan variabel, 2023.

Variabel	Min	Median	Mean	Maks
IPM	63.0	73.9	73.8	83.6
IR	120.4	264.2	282.7	570.2
Kasus	1982.0	9980.0	23849.5	212613.0
Kepadatan	9.0	103.5	783.0	17153.0
Miskin	4.3	8.4	10.1	26.0
Rasio_Pkm	0.9	1.4	1.4	2.1
Sanitasi	43.0	83.4	82.6	96.4

4.1.4.3 2024

ringkas_var(tb24, 2024)

Ringkasan variabel, 2024.

Variabel	Min	Median	Mean	Maks
IPM	67.7	74.5	74.7	84.2
IR	121.9	265.4	286.3	615.7
Kasus	2051.0	9960.0	24776.0	225841.0
Kepadatan	9.6	105.1	774.3	16685.9
Miskin	3.8	7.5	9.2	21.1
Rasio_Pkm	0.9	1.4	1.4	2.1
Sanitasi	72.8	84.6	85.0	96.8

4.2 Peta Insidens Rate (IR)

4.2.1 2022

peta_ir(tb22, 2022)

4.2.2 2023

peta_ir(tb23, 2023)

4.2.3 2024

peta_ir(tb24, 2024)

Interpretasi. Peta menunjukkan provinsi dengan warna lebih pekat memiliki insidens TB lebih tinggi. Tiga provinsi dengan IR tertinggi relatif konsisten antar tahun (mis. pada 2022: Dki Jakarta, Gorontalo, Jawa Barat), mengindikasikan beban TB yang terkonsentrasi pada wilayah-wilayah tertentu. Klik tab tiap tahun untuk membandingkan perubahan sebaran.

4.3 Standardized Incidence Ratio (SIR)

4.3.1 2022

tabel_sir(tb22, 2022)

Sepuluh provinsi SIR tertinggi (2022).

Provinsi	Observed	Expected	SIR
Dki Jakarta	54025	29582	1.83
Gorontalo	4786	3169	1.51
Jawa Barat	184406	127891	1.44
Papua	15822	11457	1.38
Banten	42429	31935	1.33
Sulawesi Utara	8784	7005	1.25
Papua Barat	3393	3053	1.11
Kepulauan Riau	5825	5525	1.05
Sumatera Utara	41057	40245	1.02
Sumatera Barat	14844	14788	1.00

peta_sir_fun(tb22, 2022)

4.3.2 2023

tabel_sir(tb23, 2023)

Sepuluh provinsi SIR tertinggi (2023).

Provinsi	Observed	Expected	SIR
Papua	6187	3179	1.95
Dki Jakarta	60362	33220	1.82
Papua Barat	2668	1658	1.61
Jawa Barat	212613	146205	1.45
Banten	52887	36537	1.45
Gorontalo	5000	3624	1.38
Sulawesi Utara	10192	7794	1.31
Sumatera Utara	49425	45332	1.09
Sulawesi Tenggara	8580	8069	1.06
Kepulauan Riau	6633	6384	1.04

peta_sir_fun(tb23, 2023)

4.3.3 2024

tabel_sir(tb24, 2024)

Sepuluh provinsi SIR tertinggi (2024).

Provinsi	Observed	Expected	SIR
Papua	6787	3305	2.05
Dki Jakarta	63403	33096	1.92
Papua Barat	2892	1728	1.67
Jawa Barat	225841	153865	1.47
Banten	56130	38623	1.45
Gorontalo	4964	3751	1.32
Sulawesi Utara	10082	7932	1.27
Sumatera Utara	50409	46897	1.07
Kepulauan Riau	7100	6812	1.04
Sulawesi Selatan	29154	28569	1.02

peta_sir_fun(tb24, 2024)

Interpretasi. SIR membandingkan jumlah kasus teramati dengan yang diharapkan bila provinsi memiliki risiko setara rata-rata nasional. Pada 2022, sebanyak 10 dari 34 provinsi memiliki SIR > 1 (risiko di atas rata-rata nasional, ditandai warna merah pada peta), sedangkan sisanya di bawah rata-rata. Provinsi dengan SIR tinggi menjadi prioritas perhatian karena bebannya melampaui ekspektasi populasi.

4.4 Analisis Pola Spasial: Moran’s I Global

moran_tahun <- function(df, th){
  peta <- indo %>% left_join(df, by="Provinsi")
  mc <- moran.mc(peta$IR, lw, nsim=999, zero.policy=TRUE)
  data.frame(Tahun=th, Moran_I=round(unname(mc$statistic),3), p_value=round(mc$p.value,3))
}
moran_tab <- rbind(moran_tahun(tb22,2022), moran_tahun(tb23,2023), moran_tahun(tb24,2024))
moran_tab %>%
  kable(caption="Tabel 2. Moran's I global IR per tahun (999 permutasi).") %>%
  kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))

Tabel 2. Moran’s I global IR per tahun (999 permutasi).

Tahun	Moran_I	p_value
2022	0.233	0.015
2023	0.273	0.008
2024	0.262	0.003

Interpretasi. Nilai Moran’s I berkisar 0.233-0.273 dan signifikan pada ketiga tahun (p < 0,05). Nilai positif menandakan autokorelasi spasial positif: provinsi dengan IR tinggi cenderung berdekatan dengan provinsi ber-IR tinggi lainnya, membentuk pola pengelompokan (spatial clustering). Konsistensi nilai antar tahun menunjukkan pola spasial yang stabil sepanjang periode, sehingga keberadaan klaster ini bukan kejadian sesaat melainkan karakteristik yang menetap.

Cara baca. Moran’s I mendekati +1 = pengelompokan spasial positif; mendekati 0 = acak. Jika p-value < 0,05, autokorelasi spasial signifikan.

4.5 Peta Klaster LISA (Local Moran)

4.5.1 2022

peta_lisa(tb22, 2022)

4.5.2 2023

peta_lisa(tb23, 2023)

4.5.3 2024

peta_lisa(tb24, 2024)

Interpretasi. Pada 2022, analisis LISA mengidentifikasi 1 provinsi klaster High-High (Banten) — yaitu provinsi ber-IR tinggi yang dikelilingi provinsi ber-IR tinggi (hotspot), dan 1 provinsi klaster Low-Low (coldspot). Klaster High-High menandai wilayah prioritas intervensi karena bebannya tinggi dan menyebar secara terkonsentrasi. Bandingkan antar tahun melalui tab.

Cara baca LISA. High-High = hotspot; Low-Low = coldspot; High-Low / Low-High = outlier spasial.

4.6 Peta Hotspot Getis-Ord Gi*

4.6.1 2022

peta_getis(tb22, 2022)

4.6.2 2023

peta_getis(tb23, 2023)

4.6.3 2024

peta_getis(tb24, 2024)

Interpretasi. Statistik Getis-Ord Gi* pada 2022 mengidentifikasi 4 provinsi sebagai hotspot (skor-Z > 1,65; kumpulan IR tinggi) dan 4 provinsi sebagai coldspot (IR rendah terkelompok). Provinsi hotspot mencakup: Maluku Utara, Dki Jakarta, Jawa Barat, Banten. Hotspot ini memperkuat temuan LISA dan menandai wilayah dengan intensitas pengelompokan kasus tertinggi yang layak diprioritaskan.

Cara baca hotspot. Gi* tinggi positif = provinsi tersebut dan tetangganya sama-sama ber-IR tinggi (hotspot). Gi* negatif = kumpulan IR rendah (coldspot). Tingkat 90/95/99% menunjukkan keyakinan statistik.

4.7 Analisis Determinan: Pemodelan Count

Outcome: jumlah kasus TB dengan offset(log(Penduduk)). Prediktor (5 covariate): IPM, % penduduk miskin, % sanitasi layak, kepadatan penduduk, dan rasio Puskesmas per kecamatan. Pemilihan model (NB vs model CAR spasial) ditentukan oleh uji autokorelasi spasial pada residual tiap tahun.

4.7.1 Eksplorasi Korelasi antar Prediktor

Matriks korelasi antar enam prediktor (tanpa IR, karena IR bukan prediktor melainkan turunan dari outcome) untuk memeriksa potensi multikolinearitas.

plot_corr <- function(df){
  vars <- df %>% dplyr::select(IPM, Miskin, Sanitasi, Kepadatan, Rasio_Pkm)
  corrplot(cor(vars), method="color", type="upper", addCoef.col="black",
           number.cex=0.7, tl.col="black", tl.srt=45)
}

4.7.1.1 2022

plot_corr(tb22)

Matriks korelasi antar prediktor (2022).

4.7.1.2 2023

plot_corr(tb23)

Matriks korelasi antar prediktor (2023).

4.7.1.3 2024

plot_corr(tb24)

Matriks korelasi antar prediktor (2024).

4.7.2 Uji Overdispersi & Autokorelasi Residual (per Tahun)

f_cov <- Kasus ~ IPM + Miskin + Sanitasi + Kepadatan + Rasio_Pkm

diagnostik <- function(df, th){
  mp <- glm(f_cov, family=poisson("log"), offset=log(Penduduk), data=df)
  disp <- sum(residuals(mp, type="pearson")^2)/mp$df.residual
  mnb <- glm.nb(update(f_cov, .~. + offset(log(Penduduk))), data=df)
  rdf <- df %>% mutate(.r=residuals(mnb, type="pearson")) %>% dplyr::select(Provinsi, .r)
  peta <- indo %>% left_join(rdf, by="Provinsi")
  mc <- moran.mc(peta$.r, lw, nsim=999, zero.policy=TRUE)
  data.frame(Tahun=th, Dispersi=round(disp,1),
             Moran_residual=round(unname(mc$statistic),3),
             p_value=round(mc$p.value,3),
             Model=ifelse(mc$p.value<0.05, "Model CAR", "NB"))
}
tabel_diag <- rbind(diagnostik(tb22,2022), diagnostik(tb23,2023), diagnostik(tb24,2024))
tabel_diag %>%
  kable(caption="Tabel 3. Uji overdispersi & autokorelasi residual NB per tahun.") %>%
  kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))

Tabel 3. Uji overdispersi & autokorelasi residual NB per tahun.

Tahun	Dispersi	Moran_residual	p_value	Model
2022	1790.5	0.071	0.165	NB
2023	1920.8	0.145	0.062	NB
2024	2343.4	0.156	0.052	NB

Interpretasi. Rasio dispersi jauh di atas 1 (1790.5-2343.4) pada semua tahun menegaskan adanya overdispersi, sehingga Poisson standar tidak memadai dan Negative Binomial menjadi titik awal yang tepat. Selanjutnya, uji Moran’s I pada residual menentukan kebutuhan model spasial: Tahun 2022, 2023, 2024 (residual acak, p > 0,05) cukup dimodelkan dengan Negative Binomial. Logika berbasis bukti ini memastikan pemilihan model didasarkan pada karakteristik data, bukan asumsi.

Logika pemilihan model. Rasio dispersi >> 1 pada ketiga tahun menegaskan penggunaan Negative Binomial (bukan Poisson). Uji Moran’s I residual kemudian menentukan: bila p > 0,05 (tidak ada autokorelasi tersisa) NB memadai; bila p < 0,05 digunakan model CAR spasial. Hasil: 2022 -> NB (residual acak), 2023 & 2024 -> Model CAR (residual berautokorelasi).

4.7.3 Model 2022: Negative Binomial

4.7.3.1 Ringkasan & IRR

m_nb22 <- glm.nb(update(f_cov, .~. + offset(log(Penduduk))), data=tb22)
summary(m_nb22)

## 
## Call:
## glm.nb(formula = update(f_cov, . ~ . + offset(log(Penduduk))), 
##     data = tb22, init.theta = 18.50380797, link = log)
## 
## Coefficients:
##                Estimate  Std. Error z value Pr(>|z|)   
## (Intercept) -3.92546584  1.31222101  -2.991  0.00278 **
## IPM         -0.01346706  0.01616374  -0.833  0.40475   
## Miskin      -0.00370327  0.01099262  -0.337  0.73620   
## Sanitasi    -0.00685009  0.00416233  -1.646  0.09982 . 
## Kepadatan    0.00019950  0.00006372   3.131  0.00174 **
## Rasio_Pkm   -0.43954797  0.18689015  -2.352  0.01868 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(18.5038) family taken to be 1)
## 
##     Null deviance: 57.060  on 33  degrees of freedom
## Residual deviance: 34.297  on 28  degrees of freedom
## AIC: 642.43
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  18.50 
##           Std. Err.:  4.46 
## 
##  2 x log-likelihood:  -628.43

irr22 <- exp(cbind(IRR=coef(m_nb22), confint(m_nb22)))
round(irr22,3) %>%
  kable(caption="Tabel 4. IRR & 95% CI - Negative Binomial 2022.") %>%
  kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))

Tabel 4. IRR & 95% CI - Negative Binomial 2022.

	IRR	2.5 %	97.5 %
(Intercept)	0.020	0.001	0.291
IPM	0.987	0.954	1.021
Miskin	0.996	0.975	1.018
Sanitasi	0.993	0.985	1.001
Kepadatan	1.000	1.000	1.000
Rasio_Pkm	0.644	0.456	0.918

Interpretasi. IRR menyatakan perubahan rate TB per kenaikan satu satuan prediktor (dengan prediktor lain konstan). Faktor yang berasosiasi signifikan (95% CI tidak memuat 1): kepadatan penduduk (IRR=1.000, meningkatkan rate TB); rasio Puskesmas per kecamatan (IRR=0.644, menurunkan rate TB). Prediktor dengan IRR di atas 1 berasosiasi dengan insidens lebih tinggi, sedangkan di bawah 1 dengan insidens lebih rendah.

4.7.3.2 Multikolinearitas (VIF)

vif22 <- car::vif(m_nb22)
data.frame(VIF=round(vif22,2)) %>%
  kable(caption="Tabel 5. VIF prediktor (2022).") %>%
  kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))

Tabel 5. VIF prediktor (2022).

	VIF
IPM	2.42
Miskin	2.06
Sanitasi	1.01
Kepadatan	20.83
Rasio_Pkm	22.21

Interpretasi. Nilai VIF melebihi 10, menandakan multikolinearitas serius (tertinggi pada Rasio_Pkm = 22.21). Dengan demikian, estimasi koefisien model dapat dianggap ditafsirkan dengan hati-hati.

Interpretasi IRR. IRR=exp(koef). IRR>1 = kenaikan covariate berasosiasi dengan kenaikan rate TB; IRR<1 = penurunan. Signifikan bila 95% CI tidak memuat 1.

4.7.4 Model CAR (Poisson Spasial): 2023 & 2024

Karena residual NB pada 2023 dan 2024 menunjukkan autokorelasi spasial, digunakan model CAR (Conditional Autoregressive) dengan likelihood Poisson dan efek acak spasial (formulasi Leroux). Pada model ini, overdispersi tidak ditangani lewat parameter dispersi (seperti NB), melainkan diserap oleh efek acak spasial dan tak terstruktur — pendekatan yang konsisten dengan Farkhan et al. yang menemukan Poisson berstruktur spasial memadai. Estimasi Bayesian via CARBayes::S.CARleroux (MCMC); koefisien di-eksponen menjadi rate ratio.

# matriks ketetanggaan biner (dipakai kedua tahun)
W <- nb2mat(nb_k, style="B", zero.policy=TRUE)
W <- pmax(W, t(W))

fit_car <- function(df){
  CARBayes::S.CARleroux(
    formula = update(f_cov, .~. + offset(log(Penduduk))),
    family  = "poisson", W = W, data = df,
    burnin = 20000, n.sample = 120000, thin = 10, verbose = FALSE)
}
ringkas_car <- function(m, th){
  tab <- exp(m$summary.results[, c("Median","2.5%","97.5%")])
  round(tab,3) %>% as.data.frame() %>%
    kable(caption=paste0("Rate Ratio (Model CAR) & 95% CrI - ", th, ".")) %>%
    kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))
}

4.7.4.1 2023

m_car23 <- fit_car(tb23)
ringkas_car <- function(m, th){
  tab <- exp(m$summary.results[, c("Mean","2.5%","97.5%")])
  colnames(tab) <- c("RateRatio","CI_2.5%","CI_97.5%")
  round(tab, 3) %>% as.data.frame() %>%
    kable(caption=paste0("Rate Ratio (Model CAR) & 95% CrI - ", th, ".")) %>%
    kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))
}

4.7.4.2 2024

m_car24 <- fit_car(tb24)
ringkas_car(m_car24, 2024)

Rate Ratio (Model CAR) & 95% CrI - 2024.

	RateRatio	CI_2.5%	CI_97.5%
(Intercept)	0.384	0.068	4.398
IPM	0.949	0.924	0.983
Miskin	0.999	0.976	1.015
Sanitasi	0.990	0.965	1.006
Kepadatan	1.000	1.000	1.000
Rasio_Pkm	0.834	0.351	1.825
tau2	1.200	1.089	1.504
rho	1.137	1.003	1.658

Interpretasi. Setelah memperhitungkan efek acak spasial melalui model CAR, faktor yang berasosiasi (95% CrI tidak memuat 1) adalah: 2023 — IPM, penduduk miskin, sanitasi layak, kepadatan penduduk, NA, NA; 2024 — IPM, kepadatan penduduk, NA, NA. Konsistensi arah rate ratio dengan model NB 2022 memperkuat keandalan temuan determinan. Interval kredibel yang lebar wajar mengingat hanya 34 area, sehingga interpretasi difokuskan pada arah asosiasi, bukan presisi titik estimasi.

Catatan. Model CAR Bayesian untuk 34 area bersifat eksploratif; interval kredibel cenderung lebar. Fokuskan interpretasi pada arah & konsistensi rate ratio dengan model NB 2022, serta faktor yang kredibel intervalnya tidak memuat 1.

5 Pembahasan

Sesuai pedoman bagian F. Bahas: temuan utama (IR, SIR, hotspot, determinan signifikan); makna epidemiologis; implikasi kebijakan; keterbatasan (data agregat/ekologis, 34 area membatasi prediktor & model spasial Bayesian, 3 tahun, kemungkinan under-reporting). Bandingkan dengan paper acuan: di sana Poisson + BYM lebih baik (514 kabupaten dengan random effect spasial), sedangkan pada skala 34 provinsi tanpa random effect, NB lebih tepat untuk mengakomodasi overdispersi. Analisis determinan difokuskan pada 2022 karena pola spasial stabil sepanjang periode (Moran’s I 0,23-0,27; p < 0,05).

Catatan data JKN. Variabel cakupan Jaminan Kesehatan Nasional tidak dimasukkan ke dalam model karena ditemukan inkonsistensi antara data peserta JKN dan jumlah penduduk: proporsi cakupan melebihi 100% pada beberapa provinsi (mencapai sekitar 387% di Papua pada 2023-2024). Hal ini kemungkinan disebabkan oleh ketidaksesuaian batas wilayah administratif dan periode pendataan, khususnya pasca-pemekaran provinsi Papua tahun 2022. Memasukkan variabel dengan inkonsistensi sebesar ini berisiko membiaskan estimasi, sehingga variabel ini dikeluarkan demi menjaga validitas model.

(Tuliskan pembahasan di sini.)

6 Kesimpulan

(Ringkasan temuan & rekomendasi proporsional.)

7 Daftar Pustaka

1. (Sumber kasus TB & treatment - SITB/Profil Kesehatan Kemenkes.)
2. (Sumber penduduk - BPS.)
3. (Sumber covariate - BPS: IPM, kemiskinan, sanitasi, kepadatan; JKN.)
4. (Shapefile batas provinsi.)
5. Farkhan A, dkk. Spatiotemporal epidemiology and associated risk factors of TB incidence and mortality in Indonesia 2017-2022. Population Health Metrics. 2026.
6. (Referensi metode: Moran’s I, LISA, Getis-Ord, Negative Binomial, CAR.)

8 Lampiran: Kontribusi Anggota

Anggota	Kontribusi
Sinta Septi Pangastuti	(…)
Winalia Agwil	(…)

---

Informasi sesi R

sessionInfo()

## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Dutch_Belgium.utf8  LC_CTYPE=Dutch_Belgium.utf8   
## [3] LC_MONETARY=Dutch_Belgium.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Dutch_Belgium.utf8    
## 
## time zone: Europe/Brussels
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] CARBayes_6.1.1   Rcpp_1.1.1-1.1   corrplot_0.95    scales_1.4.0    
##  [5] car_3.1-5        carData_3.0-6    MASS_7.3-65      kableExtra_1.4.0
##  [9] knitr_1.51       ggrepel_0.9.8    spdep_1.4-2      spData_2.3.5    
## [13] sf_1.1-1         readxl_1.5.0     lubridate_1.9.5  forcats_1.0.1   
## [17] stringr_1.6.0    dplyr_1.2.1      purrr_1.2.2      readr_2.2.0     
## [21] tidyr_1.3.2      tibble_3.3.1     ggplot2_4.0.3    tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.3.0          deldir_2.0-4       s2_1.1.11          rlang_1.2.0       
##  [5] magrittr_2.0.5     otel_0.2.0         e1071_1.7-17       compiler_4.5.2    
##  [9] png_0.1-9          systemfonts_1.3.2  vctrs_0.7.3        quantreg_6.1      
## [13] pkgconfig_2.0.3    wk_0.9.5           shape_1.4.6.1      fastmap_1.2.0     
## [17] mcmc_0.9-8         labeling_0.4.3     leafem_0.2.5       rmarkdown_2.31    
## [21] tzdb_0.5.0         MatrixModels_0.5-4 xfun_0.58          satellite_1.0.6   
## [25] glmnet_5.0         cachem_1.1.0       jsonlite_2.0.0     CARBayesdata_3.0  
## [29] mapview_2.11.4     terra_1.9-27       parallel_4.5.2     R6_2.6.1          
## [33] bslib_0.11.0       stringi_1.8.7      RColorBrewer_1.1-3 GGally_2.4.0      
## [37] boot_1.3-32        jquerylib_0.1.4    cellranger_1.1.0   iterators_1.0.14  
## [41] base64enc_0.1-6    igraph_2.3.2       splines_4.5.2      Matrix_1.7-4      
## [45] timechange_0.4.0   tidyselect_1.2.1   rstudioapi_0.19.0  abind_1.4-8       
## [49] yaml_2.3.12        codetools_0.2-20   lattice_0.22-7     withr_3.0.2       
## [53] S7_0.2.2           coda_0.19-4.1      evaluate_1.0.5     survival_3.8-3    
## [57] ggstats_0.13.0     units_1.0-1        proxy_0.4-29       xml2_1.5.2        
## [61] pillar_1.11.1      KernSmooth_2.23-26 foreach_1.5.2      stats4_4.5.2      
## [65] generics_0.1.4     sp_2.2-1           truncnorm_1.0-9    hms_1.1.4         
## [69] class_7.3-23       glue_1.8.1         tools_4.5.2        SparseM_1.84-2    
## [73] dotCall64_1.2      grid_4.5.2         MCMCpack_1.7-1     crosstalk_1.2.2   
## [77] raster_3.6-32      Formula_1.2-5      cli_3.6.6          textshaping_1.0.5 
## [81] spam_2.11-4        viridisLite_0.4.3  svglite_2.2.2      gtable_0.3.6      
## [85] sass_0.4.10        digest_0.6.39      classInt_0.4-11    htmlwidgets_1.6.4 
## [89] farver_2.1.2       htmltools_0.5.9    lifecycle_1.0.5    leaflet_2.2.3