Analisis Epidemiologi Kejadian 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)    

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, sebanyak 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 di dunia dan secara konsisten berada dalam kelompok negara prioritas pengendalian TB menurut World Health Organization (WHO, 2023). Tingginya jumlah kasus dan kematian akibat TB menunjukkan bahwa upaya pengendalian penyakit ini masih memerlukan perhatian yang serius, terutama dalam mengidentifikasi wilayah dengan risiko tinggi dan faktor-faktor yang berkontribusi terhadap tingginya kejadian penyakit. Dalam epidemiologi, distribusi penyakit dipelajari berdasarkan dimensi orang (person), tempat (place), dan waktu (time). Dimensi tempat memiliki peran penting karena kejadian TB seringkali menunjukkan variasi antarwilayah yang dipengaruhi oleh kondisi demografi, sosial ekonomi, lingkungan, maupun akses terhadap pelayanan kesehatan. Sebagai negara kepulauan yang terdiri atas 38 provinsi dengan karakteristik yang sangat beragam, Indonesia memiliki potensi terjadinya perbedaan beban TB antarwilayah. Oleh karena itu, pemahaman mengenai pola distribusi geografis TB menjadi penting untuk mendukung perencanaan program pengendalian yang lebih efektif dan tepat sasaran. Pendekatan epidemiologi spasial menggunakan indeks Moran’s I dilakukan untuk mengidentifikasi pola pengelompokan (spatial clustering) kasus Tuberkulosis (TB) akibat faktor risiko regional. Selain aspek spasial, analisis determinan juga diperlukan menggunakan pendekatan data cacah (count data). Karena data epidemiologi TB sering mengalami overdispersi (varians melebihi rata-rata), model regresi negative binomial dipilih karena lebih akurat dibanding regresi Poisson standar. Pada skala provinsi dengan jumlah unit wilayah yang terbatas, model ini menjadi alternatif yang lebih sederhana namun memadai dibandingkan pendekatan Bayesian hierarkis yang kompleks. Penelitian ini bertujuan untuk menganalisis distribusi spasial, mengevaluasi autokorelasi, dan mengidentifikasi determinan kasus TB tingkat provinsi di Indonesia guna mendukung kebijakan pengendalian berbasis wilayah yang lebih efektif.

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 Ringkasan Data TB 2022-2024

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). Menggunakan pembobot k-nearest neighbours (k=4) 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 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).

Interpretasi. DKI Jakarta tercatat sebagai IR TB tertinggi pada 2022 (480 per 100.000), diikuti Gorontalo (397) dan Jawa Barat (379). Sepuluh provinsi teratas didominasi wilayah di Jawa dan Indonesia timur, menunjukkan beban insidens yang terkonsentrasi pada provinsi tertentu.

4.1.2.2 2023

plot_top10(tb23, 2023)

Sepuluh provinsi IR tertinggi (2023).

Interpretasi (2023). Papua tertinggi (570 per 100.000), menunjukkan pergeseran dari DKI Jakarta tahun sebelumnya.

4.1.2.3 2024

plot_top10(tb24, 2024)

Sepuluh provinsi IR tertinggi (2024).

Interpretasi (2024). Papua tetap tertinggi (616 per 100.000), mengindikasikan beban insidens yang menetap di wilayah timur.

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

Interpretasi. Pada 2022, IR TB antarprovinsi berkisar 105,7-480,2 per 100.000 (rata-rata 245,3), menunjukkan disparitas yang lebar. Jumlah kasus sangat bervariasi (1.738-184.406), sejalan dengan rentang kepadatan penduduk yang ekstrem (10-17.031 per km²), sehingga populasi perlu diperhitungkan sebagai offset dalam pemodelan.

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

Interpretasi. Pada 2023, rata-rata IR meningkat menjadi 282,7 per 100.000 dengan nilai maksimum mencapai 570,2, menandakan kenaikan beban insidens dibanding tahun sebelumnya. Sebaran covariate sosioekonomi (IPM, kemiskinan, sanitasi) relatif stabil antartahun.

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

Interpretasi. Pada 2024, rata-rata IR kembali naik menjadi 286,3 per 100.000 (maksimum 615,7), melanjutkan tren peningkatan. Disparitas antarprovinsi yang konsisten lebar sepanjang periode menegaskan perlunya analisis spasial dan penyesuaian populasi dalam model.

4.2 Peta Insidens Rate (IR)

4.2.1 2022

peta_ir(tb22, 2022)

Interpretasi. Pada 2022, IR TB tertinggi terkonsentrasi di provinsi padat di Pulau Jawa — DKI Jakarta (480), Jawa Barat (379), dan Banten (349) — serta beberapa wilayah Indonesia timur seperti Gorontalo (397) dan Papua (363). Kedekatan warna gelap antarprovinsi di Jawa bagian barat mengisyaratkan potensi pengelompokan spasial yang diuji melalui Moran’s I.

4.2.2 2023

peta_ir(tb23, 2023)

Interpretasi. Pada 2023, beban insidens tertinggi bergeser ke Papua (570), menandai meningkatnya IR di Indonesia timur, sementara provinsi padat di Jawa tetap berada pada kategori tinggi. Pola sebaran yang serupa dengan tahun sebelumnya mengindikasikan konsentrasi geografis yang menetap.

4.2.3 2024

peta_ir(tb24, 2024)

Interpretasi. Pada 2024, Papua kembali tercatat sebagai IR tertinggi (616), mempertegas beban TB yang menetap di wilayah timur, diiringi insidens yang konsisten tinggi di provinsi-provinsi Jawa. Kestabilan pola spasial sepanjang 2022-2024 memperkuat temuan autokorelasi spasial yang signifikan.

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-2024, 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.") %>%
  kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))

Tabel 2. Moran’s I global IR per tahun.

Tahun	Moran_I	p_value
2022	0.233	0.013
2023	0.273	0.011
2024	0.262	0.006

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)

Interpretasi. Pada 2022, klaster High-High (hotspot) berada di Jawa bagian barat (merah) dan klaster Low-Low (coldspot) di Jawa bagian timur (biru). Pada tahap ini, konsentrasi beban TB masih berpusat di Pulau Jawa.

4.5.2 2023

peta_lisa(tb23, 2023)

Interpretasi. Pada 2023, klaster High-High bergeser ke Indonesia timur (Papua dan Papua Barat, merah), sementara Pulau Jawa berubah menjadi klaster Low-Low (coldspot). Terdapat pula outlier Low-High di sekitar Maluku — wilayah ber-IR rendah yang dikelilingi tetangga ber-IR tinggi. Pergeseran ini menandai perpindahan pusat beban TB dari barat ke timur.

4.5.3 2024

peta_lisa(tb24, 2024)

Interpretasi. Pada 2024, klaster High-High tetap di Papua dan Papua Barat, dan Jawa konsisten sebagai klaster Low-Low. Pola ini mempertegas perpindahan dan menetapnya hotspot TB di kawasan timur Indonesia pada dua tahun terakhir periode pengamatan.

4.6 Peta Hotspot Getis-Ord Gi*

4.6.1 2022

peta_getis(tb22, 2022)

Interpretasi. Pada 2022, analisis Getis-Ord Gi* mengidentifikasi hotspot (kumpulan IR tinggi) di Pulau Jawa, mempertegas temuan klaster LISA. Hotspot menandai wilayah dengan intensitas pengelompokan kasus tertinggi yang layak diprioritaskan dalam intervensi pengendalian TB.

4.6.2 2023

peta_getis(tb23, 2023)

Interpretasi. Pada 2023, hotspot Getis-Ord Gi* justru terkonsentrasi di Indonesia bagian timur — terutama Papua dan Papua Barat (Hotspot 95%) — menandai pergeseran dari pola 2022 yang berpusat di Jawa. Sementara itu, sebagian Pulau Jawa bagian tengah-timur teridentifikasi sebagai coldspot. Pergeseran ini sejalan dengan naiknya IR Papua menjadi tertinggi pada 2023 (570 per 100.000).

4.6.3 2024

peta_getis(tb24, 2024)

Interpretasi. Pada 2024, hotspot Getis-Ord Gi* di Papua menjadi tingkat tertinggi (Hotspot 99%), mempertegas konsentrasi beban TB di Indonesia bagian timur yang mulai tampak pada 2023. Pulau Jawa bagian tengah-timur konsisten menjadi coldspot. Penguatan hotspot Papua sejalan dengan IR tertinggi yang kembali tercatat di Papua (616 per 100.000).

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).

Interpretasi. Sebagian besar prediktor menunjukkan korelasi rendah hingga sedang. Namun terdapat korelasi yang sangat tinggi antara Kepadatan penduduk dan Rasio Puskesmas (r = 0,97), mengindikasikan potensi multikolinearitas serius antar kedua variabel. Korelasi sedang juga terlihat antara IPM dengan kemiskinan (r = -0,67) dan IPM dengan rasio Puskesmas (r = 0,55). Temuan ini perlu dikonfirmasi melalui nilai VIF pada tahap pemodelan; bila VIF tinggi, salah satu dari variabel yang berkorelasi kuat sebaiknya dipertimbangkan untuk dikeluarkan.

4.7.1.2 2023

plot_corr(tb23)

Matriks korelasi antar prediktor (2023).

Interpretasi. Pada 2023, korelasi antar prediktor seluruhnya tergolong rendah hingga sedang (|r| < 0,7), sehingga tidak ada indikasi multikolinearitas yang mengkhawatirkan. Korelasi tertinggi terdapat antara IPM dan persentase kemiskinan (r = -0,69), yang secara substansi wajar karena daerah dengan IPM tinggi cenderung memiliki tingkat kemiskinan lebih rendah. Berbeda dengan 2022, korelasi Kepadatan-Rasio Puskesmas pada 2023 hanya -0,24, menguatkan dugaan adanya anomali data pada tahun 2022.

4.7.1.3 2024

plot_corr(tb24)

Matriks korelasi antar prediktor (2024).

Interpretasi. Pada 2024, seluruh korelasi antar prediktor tergolong rendah hingga sedang (|r| < 0,6), menandakan tidak ada masalah multikolinearitas. Korelasi tertinggi tetap antara IPM dan kemiskinan (r = -0,58), yang wajar secara substansi. Konsistensi pola korelasi pada 2023 dan 2024 (Kepadatan-Rasio Puskesmas keduanya sekitar -0,24) memperkuat dugaan bahwa korelasi ekstrem pada 2022 merupakan anomali akibat pencilan DKI Jakarta, bukan hubungan yang sesungguhnya.

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.176	NB
2023	1920.8	0.145	0.060	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 model yang tepat. Selanjutnya, uji Moran’s I pada residual menentukan kebutuhan model spasial: Tahun 2022, 2023, 2024 (residual acak, p > 0,10) 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,10 (tidak ada autokorelasi pada residual) NB memadai; bila p < 0,10 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,4) %>%
  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.0197	0.0013	0.2912
IPM	0.9866	0.9543	1.0210
Miskin	0.9963	0.9748	1.0184
Sanitasi	0.9932	0.9845	1.0015
Kepadatan	1.0002	1.0001	1.0003
Rasio_Pkm	0.6443	0.4560	0.9182

Interpretasi. Berdasarkan pemodelan Negative Binomial tahun 2022, faktor signifikan yang memengaruhi variable kejadian TB adalah Kepadatan Penduduk (sebagai faktor pemicu/peningkat laju kejadian) dengan IRR= 1,002; CI 1,001–1,003) dan Rasio Puskesmas (sebagai faktor pengendali/penurun laju kejadian) dengan IRR=0,644; 95% CI 0,456–0,9182. Hasil ini menyarankan bahwa intervensi kebijakan pemerintah sebaiknya difokuskan pada pemerataan sarana kesehatan (Puskesmas), terutama di wilayah-wilayah dengan tingkat kepadatan penduduk yang tinggi.

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.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.

4.7.4 Model CAR (Poisson Spasial): 2023 & 2024

Karena residual NB pada 2023 dan 2024 menunjukkan autokorelasi spasial, digunakan model CAR (Conditional Autoregressive). Pada model CAR, overdispersi tidak ditangani melalui parameter dispersi (seperti pada Negative Binomial), melainkan diserap oleh efek acak spasial dan tak terstruktur. Formulasi Leroux (Leroux et al., 2000) yang digunakan di sini secara eksplisit memisahkan parameter overdispersi dari kekuatan dependensi spasial.

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

fit_car <- function(df){
   #set.seed(2606)
  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("Mean","2.5%","97.5%")])
  round(tab,3) %>% as.data.frame() %>%
    kable(caption=paste0("Rate Ratio (Model CAR) & 95% CI - ", th, ".")) %>%
    kable_styling(full_width=FALSE, bootstrap_options=c("striped","hover"))
}

4.7.4.1 2023

m_car23 <- fit_car(tb23)
ringkas_car(m_car23, 2023)

Rate Ratio (Model CAR) & 95% CI - 2023.

	Mean	2.5%	97.5%
(Intercept)	0.391	0.035	3.092
IPM	0.963	0.932	0.985
Miskin	0.984	0.963	1.004
Sanitasi	0.986	0.975	0.997
Kepadatan	1.000	1.000	1.000
Rasio_Pkm	0.516	0.203	0.928
tau2	1.210	1.081	1.522
rho	1.268	1.010	2.039

Interpretasi (2023). nilai parameter\(\rho\) dan \(\tau^2\)yang signifikan menegaskan bahwa aspek spasial (konteks wilayah dan efek ketetanggaan) memegang peranan krusial, sehingga strategi penanganan masalah ini wajib dirancang secara regional terintegrasi, bukan hanya fokus per daerah secara terisolasi

4.7.4.2 2024

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

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

	Mean	2.5%	97.5%
(Intercept)	0.271	0.014	3.837
IPM	0.958	0.937	0.982
Miskin	0.981	0.957	1.022
Sanitasi	0.991	0.976	1.006
Kepadatan	1.000	1.000	1.000
Rasio_Pkm	0.680	0.343	1.268
tau2	1.217	1.099	1.555
rho	1.104	1.002	1.520

Interpretasi (2024). Hasil pemodelan CAR tahun 2024 menunjukkan konsistensi dengan pola tahun sebelumnya (2023), di mana IPM menjadi satu-satunya faktor struktural yang secara signifikan mampu menekan laju kejadian. Selain itu, bertahannya nilai parameter\(\rho\) dan \(\tau^2\)yang signifikan menegaskan bahwa aspek spasial (konteks wilayah dan efek ketetanggaan) memegang peranan krusial, sehingga strategi penanganan masalah ini wajib dirancang secara regional terintegrasi, bukan hanya fokus per daerah secara terisolasi

5 Pembahasan

5.1 Temuan Utama dan Maknanya

Penelitian ini menunjukkan bahwa angka kejadian tuberkulosis di Indonesia meningkat sepanjang 2022-2024, dengan total kasus ternotifikasi naik dari 724.309 menjadi 842.383 dan rata-rata incidence rate (IR) antarprovinsi dari 245,3 menjadi 286,3 per 100.000 penduduk. Analisis SIR memjelaskan bahwa beban tidak merata: pada 2022 sebanyak 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 global yang positif dan signifikan pada ketiga tahun (I=0,233-0,273; p<0,05) menegaskan adanya autokorelasi spasial: provinsi dengan IR tinggi cenderung berdekatan secara geografis, membentuk pengelompokan yang tidak acak. Keberadaan autokorelasi ini konsisten sepanjang periode pengamatan.

5.2 Pergeseran Spasial Kejadian TB dari Barat ke Timur

Temuan paling menonjol dari penelitian ini adalah pergeseran lokasi pusat kejadian TB. Meskipun autokorelasi spasial global konsisten signifikan, analisis klaster lokal (LISA) dan hotspot (Getis-Ord Gi*) menunjukkan bahwa lokasi klaster berpindah seiring waktu. Pada 2022, hotspot (klaster High-High) terkonsentrasi di Pulau Jawa bagian barat, dengan Jawa bagian timur sebagai coldspot. Namun pada 2023, hotspot bergeser ke Indonesia bagian timur, terutama Papua dan Papua Barat (Hotspot 95%), sementara Jawa berubah menjadi coldspot. Pada 2024, hotspot Papua menguat hingga tingkat tertinggi (Hotspot 99%) dan meluas ke wilayah sekitarnya. Pergeseran ini sejalan dengan perpindahan provinsi ber-IR tertinggi dari DKI Jakarta (2022) ke Papua (2023 dan 2024).

5.3 Determinan Insidens TB

Pemilihan model dilakukan dengan melakukan pengujian overdispersi untuk mengonfirmasi penggunaan regresi Negative Binomial (rasio dispersi jauh di atas 1 pada ketiga tahun), dan uji autokorelasi spasial residual menentukan kebutuhan model spasial. Pada 2022 residual tidak berautokorelasi sehingga model NB memadai, sedangkan pada 2023 dan 2024 residual menunjukkan autokorelasi spasial sehingga digunakan model CAR (Poisson spasial, formulasi Leroux).

Pada model Negative Binomial (NB) 2022, rasio Puskesmas berasosiasi signifikan dengan penurunan laju kasus TB. Hal ini mengindikasikan bahwa ketersedian fasilitas kesehatan primer yang memadai berbanding lurus dengan rendahnya insidens TB. Sementara itu, pada model CAR 2023, variable IPM, sanitasi dan rasio puskesmas merupakan variabel yang signifikan menekan insidens TB. Model CAR 2024, hanya variabel IPM yang signifikan menekan insidens TB. Hasil ini menunjukkan bahwa faktor struktural (IPM) dan ketersediaan fasilitas kesehatan (rasio Puskesmas) merupakan determinan penting dalam mengendalikan laju kejadian TB di Indonesia.

Catatan data. Variabel cakupan Jaminan Kesehatan Nasional (JKN) 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). Selain itu, nilai rasio Puskesmas DKI Jakarta pada 2022 teridentifikasi sebagai pencilan yang memengaruhi struktur korelasi pada tahun tersebut, sehingga estimasi koefisien kepadatan dan rasio Puskesmas pada model 2022 ditafsirkan dengan hati-hati.

6 Daftar Pustaka

1. Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x
2. Blangiardo, M., & Cameletti, M. (2015). Spatial and spatio-temporal Bayesian models with R-INLA. John Wiley & Sons. https://doi.org/10.1002/9781118950203
3. Cameron, A. C., & Trivedi, P. K. (2013). Regression analysis of count data (2nd ed.). Cambridge University Press. https://doi.org/10.1017/CBO9781139013567
4. Farkhan, A., Pakasi, T. T., Sulistyo, S., et al. (2026). Spatiotemporal epidemiology and associated risk factors of tuberculosis incidence and mortality in Indonesia 2017–2022: A nationwide space-time hierarchical analysis. Population Health Metrics, 24, Artikel 12. https://doi.org/10.1186/s12963-026-00458-5
5. Hilbe, J. M. (2011). Negative binomial regression (2nd ed.). Cambridge University Press. https://doi.org/10.1017/CBO9780511973420
6. Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2004). Applied linear statistical models (5th ed.). McGraw-Hill Irwin. ISBN: 978-0073108742
7. Mirahmadizadeh, A., Heiran, A., & Delam, H. (2022). Correlation between human development index and its components with COVID-19 indices: A global level ecologic study. BMC Public Health, 22(1), Artikel 1549. https://doi.org/10.1186/s12889-022-13959-1
8. Munayco, C. V., Mújica, O., Del Granado, M., Barceló, A., & Jiang, W. (2015). Social determinants and inequalities in tuberculosis incidence in Latin America and the Caribbean. Revista Panamericana de Salud Pública, 38(3), 177–185. https://doi.org/10.1590/S1020-49892015000300005
9. Ord, J. K., & Getis, A. (1995). Local spatial autocorrelation statistics: Distributional issues and an application. Geographical Analysis, 27(4), 286–306. https://doi.org/10.1111/j.1538-4632.1995.tb00912.x
10. Spence, D. P., Hotchkiss, J., Williams, C. S. D., & Davies, P. D. O. (1993). Tuberculosis and poverty. BMJ, 307(6907), 759–761. https://doi.org/10.1136/bmj.307.6907.759
11. World Health Organization. (2023). Global tuberculosis report 2023. World Health Organization. https://iris.who.int/handle/10665/373828

7 Lampiran: Kontribusi Anggota

Anggota	Kontribusi
Sinta Septi Pangastuti	Pengumpulan Data dan Pembuatan Makalah
Winalia Agwil	Pembuatan file Rmd dan Analisis

---

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