Laporan Epidemiologi UTS Firgie

Link Dashboard: https://firgiealdiansyahf.shinyapps.io/Dashboard_Epidemiologi/

Pendahuluan

Latar Belakang

Tuberculosis (TBC) adalah penyakit infeksi yang disebabkan oleh bakteri Mycobacterium tuberculosis yang umumnya menyerang paru-paru, namun dapat pula menyerang organ lain. Penyakit ini tetap menjadi masalah kesehatan masyarakat global dan nasional karena tingginya kasus yang terjadi. Penularan TBC dapat terjadi melalui udara saat penderita batuk, bersin, dan berbicara [1] oleh karena itu resiko penyebaran nya sangat besar terutama di lingkungan dengan kepadatan penduduk yang tinggi. Menurut Global Tuberculosis Report 2024 oleh World Health Organization (WHO), terdapat sekitar 10,6 juta kasus baru TBC di seluruh dunia dan 1,3 juta kematian akibat TBC [1]. Kondisi tersebut menunjukkan bahwa TBC masih menjadi tantangan serius bagi sistem kesehatan global, termasuk di Indonesia.

Indonesia termasuk salah satu negara dengan beban TBC tertinggi di dunia, menempati urutan kedua setelah India. WHO memperkirakan terdapat sekitar 1 juta kasus baru TBC setiap tahun di Indonesia, atau sekitar 394 kasus per 100.000 penduduk [1]. Laporan Kementerian Kesehatan RI tahun 2023 menunjukkan bahwa TBC menjadi salah satu dari sepuluh penyebab kematian tertinggi di Indonesia, bersama dengan stroke dan penyakit jantung [2]. Beban TBC di Indonesia tidak hanya berdampak pada aspek medis, tetapi juga pada aspek sosial dan ekonomi masyarakat, karena penyakit ini cenderung berkaitan dengan kemiskinan, kepadatan hunian, dan rendahnya tingkat pendidikan[3][4].

Provinsi Jawa Barat merupakan salah satu wilayah dengan beban TBC tertinggi di Indonesia [5]. Sebagai provinsi berpenduduk terbesar dengan lebih dari 49 juta jiwa, variasi kondisi sosial ekonomi antar kabupaten/kota memengaruhi tingkat penularan penyakit. Daerah dengan kepadatan tinggi dan kondisi lingkungan yang kurang baik cenderung memiliki jumlah kasus lebih besar dibanding wilayah lain[6][7]. Hal ini sejalan dengan hasil penelitian yang menunjukkan bahwa determinan sosial seperti kepadatan penduduk, status sosial ekonomi, dan kualitas hunian memiliki hubungan yang signifikan terhadap angka kejadian TBC di Indonesia [8]. Oleh karena itu, analisis epidemiologi terhadap distribusi spasial TBC di Jawa Barat menjadi penting untuk mengetahui pola sebaran penyakit dan faktor wilayah yang berpengaruh. Hasil penelitian ini diharapkan dapat menjadi dasar bagi upaya intervensi dan kebijakan pengendalian TBC dalam rangka mendukung target Indonesia Bebas TBC 2030. [5]

Rumusan Masalah

Berdasarkan latar belakang yang telah diuraikan, permasalahan utama dalam penelitian ini adalah tingginya angka kejadian Tuberkulosis (TBC) di Provinsi Jawa Barat yang masih menjadi tantangan dalam upaya pengendalian penyakit menular. Perbedaan jumlah kasus antar kabupaten/kota menunjukkan bahwa terdapat variasi kondisi wilayah yang dapat memengaruhi tingkat penularan TBC. Faktor-faktor seperti kepadatan penduduk, kondisi lingkungan tempat tinggal, serta aspek sosial ekonomi masyarakat diduga berperan dalam meningkatkan risiko terjadinya penyakit ini. Oleh karena itu, penelitian ini difokuskan untuk mengidentifikasi tingkat kejadian TBC di berbagai wilayah di Jawa Barat dan menganalisis hubungan antara faktor-faktor wilayah tersebut dengan tingginya kasus TBC guna mendukung perumusan kebijakan pengendalian yang lebih tepat dan efektif.

Tujuan

Penelitian ini bertujuan untuk menggambarkan kondisi epidemiologi penyakit Tuberkulosis (TBC) di Provinsi Jawa Barat berdasarkan data kasus yang tersedia. Secara khusus, penelitian ini bertujuan untuk: 1. Mengetahui tingkat kejadian dan distribusi kasus TBC di setiap kabupaten/kota di Provinsi Jawa Barat.

2. Menganalisis adanya perbedaan laju kasus TBC antar wilayah di Provinsi Jawa Barat.

3. Mengidentifikasi faktor-faktor wilayah seperti kepadatan penduduk dan kondisi lingkungan dengan tingginya angka kejadian TBC.

Tinjauan Pustaka

Agent-Host-Environment

Dalam epidemiologi penyakit menular, seperti TBC dipengaruhi oleh tiga komponen utama, yaitu agent-host-environment: 1. Agent penyebab penyakit TBC adalah bakteri Mycobacterium tuberculosis. Penularan terjadi melalui udara saat penderita batuk, bersin, atau berbicara. Droplet yang mengandung bakteri dapat terhirup oleh orang lain, kemudian menginfeksi jaringan paru-paru.

2. Host adalah individu yang terinfeksi atau berisiko terinfeksi. Faktor yang memengaruhi kerentanan seseorang terhadap TBC antara lain status daya tahan tubuh yang lemah serta adanya penyakit penyerta seperti HIV/AIDS atau diabetes melitus. Selain itu, kelompok usia produktif cenderung lebih rentan karena sering berinteraksi di lingkungan padat dan tempat kerja tertutup.

3. Environment (lingkungan) berperan penting dalam penularan TBC. Lingkungan dengan kepadatan penduduk tinggi dan tempat tinggal tidak layak huni meningkatkan risiko penyebaran bakteri TBC.

Ukuran Asosiasi

Dalam studi epidemiologi, ukuran asosiasi digunakan untuk mengetahui hubungan antara faktor risiko dengan kejadian penyakit. Ukuran asosiasi menggambarkan seberapa besar pengaruh suatu variabel bebas terhadap variabel terikat.

Beberapa ukuran asosiasi yang umum digunakan antara lain relative risk (RR), odds ratio (OR), dan attributable risk (AR). Nilai OR atau RR yang lebih besar dari 1 menunjukkan bahwa faktor tersebut meningkatkan risiko penyakit, sedangkan nilai di bawah 1 menurunkan risiko penyakit.

Dalam konteks penelitian TBC di Jawa Barat, ukuran asosiasi dapat digunakan untuk melihat hubungan antara faktor-faktor wilayah seperti kepadatan penduduk, kemiskinan, atau kondisi lingkungan dengan jumlah kasus TBC yang terjadi. Dengan demikian, ukuran asosiasi membantu peneliti menilai sejauh mana faktor tersebut berperan dalam penyebaran penyakit TBC.

Desain Studi

Desain studi yang digunakan dalam penelitian ini adalah cross-sectional. Studi cross-sectional merupakan rancangan penelitian observasional yang dilakukan dengan cara mengamati variabel paparan (faktor risiko) dan variabel hasil (kejadian penyakit) pada waktu yang bersamaan. Tujuannya adalah untuk mengetahui hubungan antara faktor risiko dan penyakit tanpa menelusuri urutan waktu kejadian.

Kelebihan desain ini adalah pelaksanaannya relatif cepat dan dapat dilakukan menggunakan data sekunder. Selain itu, desain ini cocok digunakan untuk menggambarkan distribusi penyakit dan faktor risikonya dalam suatu populasi pada periode tertentu. Kelemahannya adalah tidak dapat menentukan hubungan sebab-akibat secara langsung karena waktu paparan dan kejadian penyakit diukur bersamaan.

Dalam penelitian ini, desain cross-sectional digunakan untuk menggambarkan kondisi epidemiologi Tuberkulosis (TBC) di Jawa Barat berdasarkan data kasus per kabupaten/kota. Penelitian ini juga bertujuan untuk menilai hubungan antara faktor-faktor seperti kepadatan penduduk dan kondisi lingkungan dengan tingkat kejadian TBC pada periode tahun pengamatan tertentu.

Metodologi

Sumber Data

Data yang digunakan merupakan data sekunder yang bersumber dari Open Data Jawa Barat. Data diperoleh dari publikasi resmi pemerintah Daerah Jawa Barat yang dapat diakses secara umum. Data yang digunakan mencakup data jumlah TBC, kepadatan penduduk, dan persentase rumah tangga yang menempati tempat huni yang layak pada kota/kabupaten yang berjumlah 27 wilayah di Jawa Barat.

Variabel

Variabel independen (Y) dan variabel dependen (X) yang digunakan adalah sebagai berikut:

Variabel	Keterangan
Y	Jumlah Kasus TBC
X1	Kepadatan Penduduk
X2	Persentase Rumah Tangga yang Menempati Rumah Layak Huni

Metode Analisis

Penelitian ini menggunakan desain studi Cross-Sectional dengan pendekatan Regresi Binomial Negatif. Penelitian ini bertujuan utnuk melihat pengaruh dari varibel prediktor yang memengaruhi variabel respon. Tahapan analisis yang dilakukan sebagai berikut :

1. Eksplorasi Data : Dilakukan eksplorasi data dengan melakukan analisis deskriptif, membuat peta sebaran

2. Uji Autokorelasi Spasial : Mengukur adanya ketergantungan spasial menggunakan statistik Moran’s I dan menggunakan local Moran’s I (LISA) untuk mendeteksi klaster wilayah

3. Menghitung Prevalensi untuk mengetahui proporsi individu yang memiliki penyakit pada waktu tertentu

4. Menggunakan Regresi Binomial sebagai model regresi terbaik

5. Menghitung Prevalensi Ratio sebagai ukuran asosiasi yang digunakan

Alur Kerja

1. Pengumpulan Data : Data diperoleh dari Open Data Jabar dengan variabel yang digunakan yaitu jumlah kasus TBC, kepadatan penduduk, dan persentase rumah tangga yang menghuni rumah yang layak

2. Eksplorasi Data : Dilakukan analisis deskriptif dan menvisualisasikan sebaran TBC

3. Uji Autokorelasi Spasial : Menguji adanya pola ketergantungan spasial secara global dan lokal

4. Menghitung Prevalensi : Untuk melihat proporsi individu yang memiliki penyakit pada waktu tertentu

5.Memilih Analisis Terbaik : Dengan melihat dispersi dan AIC yang dihasilkan Menghitung Prevalence Ratio : Untuk melihat hubungan antara faktor paparan dan penyakit

Hasil dan Pembahasan

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.3

## Warning: package 'ggplot2' was built under R version 4.3.3

## Warning: package 'tidyr' was built under R version 4.3.3

## Warning: package 'readr' was built under R version 4.3.3

## Warning: package 'purrr' was built under R version 4.3.3

## Warning: package 'dplyr' was built under R version 4.3.3

## Warning: package 'stringr' was built under R version 4.3.3

## Warning: package 'forcats' was built under R version 4.3.3

## Warning: package 'lubridate' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(MASS)

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select

library(broom)

## Warning: package 'broom' was built under R version 4.3.3

library(dplyr)
library(sf)

## Warning: package 'sf' was built under R version 4.3.3

## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(spdep)

## Warning: package 'spdep' was built under R version 4.3.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.3.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(tmap)
library(geodata)

## Loading required package: terra

## Warning: package 'terra' was built under R version 4.3.3

## terra 1.8.29
## 
## Attaching package: 'terra'
## 
## The following object is masked from 'package:MASS':
## 
##     area
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

setwd("C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Epidemiologi")
data <- read.csv("DATA UTS EPIDEMIOLOGI.csv", sep = ";", dec =".")
data

str(data)

## 'data.frame':    27 obs. of  5 variables:
##  $ KABKOT  : chr  "BANDUNG" "BANDUNG BARAT" "BEKASI" "BOGOR" ...
##  $ Populasi: int  3753120 1884190 3273870 5682300 1259230 2584990 2387960 2716950 1914040 2554380 ...
##  $ Y       : int  14012 4569 15458 29110 3477 12197 9587 9516 6092 13474 ...
##  $ X1      : int  2156 1468 2617 1899 789 712 2228 876 922 1335 ...
##  $ X2      : num  54.2 38.9 63.8 46.7 67.2 ...

# Ubah nama kolom agar mudah dipakai
data <- data %>%
  rename(
    kabkot = KABKOT,
    population = Populasi,
    cases = Y,
    density = X1,
    housing = X2
  )
# 2. Hitung Prevalensi TBC per 100.000 penduduk
data <- data %>%
  mutate(
    prevalence_per100k = (cases / population) * 100000,
    housing_prop = housing / 100   # ubah persen jadi proporsi
  )

# Lihat hasil
head(data)

# --- 1. Ambil shapefile Indonesia ---
indo_map <- readRDS("C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Epidemiologi/gadm41_IDN_2_pk.rds") %>%
  st_as_sf()

## Samain nama kab/kota di shapefile
jabar <- indo_map %>%
  filter(NAME_1 == "Jawa Barat") %>%
  filter(!NAME_2 %in% c("Waduk Cirata", "Waduk Jatiluhur")) %>%
  mutate(KABKOT = str_squish(toupper(NAME_2)))



# Samain juga di data CSV (buat jaga-jaga)
data <- data %>%
  mutate(
    KABKOT = str_squish(toupper(kabkot))
  )


# --- 4. Join shapefile dan data ---
data_spatial <- left_join(jabar, data, by = "KABKOT")

# --- 5. Cek kesesuaian nama ---
setdiff(unique(jabar$KABKOT), unique(data$KABKOT))

## character(0)

Tabel Hasil

data_deskriptif <- data %>%
  dplyr::select(cases, population, density, housing) %>%
  dplyr::summarise(
    across(everything(), 
      list(
        mean = ~mean(., na.rm = TRUE),
        median = ~median(., na.rm = TRUE),
        sd = ~sd(., na.rm = TRUE),
        min = ~min(., na.rm = TRUE),
        max = ~max(., na.rm = TRUE)
      ),
      .names = "{.col}_{.fn}"
    )
  )

print(data_deskriptif)

##   cases_mean cases_median cases_sd cases_min cases_max population_mean
## 1   8506.778         6092 6240.981       968     29110         1864637
##   population_median population_sd population_min population_max density_mean
## 1           1884190       1228053         209790        5682300     3910.926
##   density_median density_sd density_min density_max housing_mean housing_median
## 1           1468   4668.116         385       15176     58.72407          61.91
##   housing_sd housing_min housing_max
## 1   17.09696       32.83       85.78

Moran’s

# --- 5. Hitung Moran’s I Global ---
# Konversi sf → Spatial
map_sp <- as_Spatial(data_spatial)

# Buat matriks tetangga (queen contiguity)
nb <- poly2nb(map_sp)

# Buat bobot spasial (row-standardized)
lw <- nb2listw(nb, style = "W")

# Uji Moran's I global
moran_global <- moran.test(data_spatial$cases, lw)
moran_global

## 
##  Moran I test under randomisation
## 
## data:  data_spatial$cases  
## weights: lw    
## 
## Moran I statistic standard deviate = 2.9347, p-value = 0.00167
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.35148462       -0.04000000        0.01779563

Hasil uji Moran’s I menunjukkan nilai Moran I = 0.3515 dengan p-value = 0.00167, yang berarti terdapat autokorelasi spasial positif yang signifikan pada kasus TBC di Jawa Barat. Artinya, daerah dengan jumlah kasus TBC tinggi cenderung berdekatan dengan daerah lain yang juga memiliki kasus tinggi, dan sebaliknya. Dengan demikian, penyebaran kasus TBC di Jawa Barat tidak bersifat acak, melainkan menunjukkan pola spasial yang mengelompok (clustered).

Peta Kasus Sebaran

tmap_mode("plot")

## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`

tm_shape(data_spatial) +
  tm_polygons(
    col = "cases",           # pakai kolom 'cases' bukan 'Y'
    palette = "Reds",        # warna gradasi merah
    title = "Jumlah Kasus TBC"
  ) +
  tm_borders(lwd = 0.7, col = "gray40") +
  tm_layout(
    title = "Sebaran Kasus TBC di Kabupaten/Kota Jawa Barat (2024)",
    legend.outside = TRUE,
    frame = FALSE
  )

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Peta Klaster LISA

# --- 8. Hitung Local Moran’s I (LISA) ---
local_moran <- localmoran(st_drop_geometry(data_spatial)$cases, lw)
head(local_moran)

##            Ii        E.Ii     Var.Ii        Z.Ii Pr(z != E(Ii))
## 1 -0.03867622 -0.02966099 0.08017623 -0.03183863    0.974600737
## 2 -0.14554706 -0.01949157 0.06556356 -0.49230075    0.622506751
## 3  1.05549883 -0.05775959 1.41500884  0.93587190    0.349339144
## 4  1.80599497 -0.04838894 0.36582139  3.06595023    0.002169794
## 5  1.47678414 -0.44997849 0.56976067  2.55259539    0.010692363
## 6  0.70568331 -0.03086110 0.12960434  2.04592245    0.040764001

# --- 9. Klasifikasi LISA (HH, LL, HL, LH) ---

# 1️⃣ Ambil nilai kasus dan rata-ratanya
cases <- st_drop_geometry(data_spatial)$cases
mean_val <- mean(cases, na.rm = TRUE)

# 2️⃣ Hitung nilai lag spasial (rata-rata tetangga)
lag_cases <- lag.listw(lw, cases)

# 3️⃣ Buat klasifikasi berdasarkan kuadran LISA
lisa_cluster <- ifelse(cases > mean_val & lag_cases > mean(lag_cases), "High-High",
                ifelse(cases < mean_val & lag_cases < mean(lag_cases), "Low-Low",
                ifelse(cases > mean_val & lag_cases < mean(lag_cases), "High-Low",
                ifelse(cases < mean_val & lag_cases > mean(lag_cases), "Low-High", "Undefined"))))

# 4️⃣ Tambahkan ke data spasial
data_spatial$lisa_cluster <- lisa_cluster

# 5️⃣ Cek distribusi klaster
table(data_spatial$lisa_cluster)

## 
## High-High  High-Low  Low-High   Low-Low 
##         7         4         4        11

# 6️⃣ Visualisasikan peta LISA
tm_shape(data_spatial) +
  tm_polygons(
    col = "lisa_cluster",
    palette = c("red", "blue", "pink", "lightblue", "gray"),
    title = "Klaster LISA (Local Moran’s I)"
  ) +
  tm_borders() +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE
  )

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'

Analisis Regresi Poisson

# 3. Model Cross-Sectional → Prevalence Rate Ratio (PRR) via Poisson Regression
model_pois <- glm(
  cases ~ density + housing_prop + offset(log(population)),
  family = poisson(link = "log"),
  data = data
)

summary(model_pois)

## 
## Call:
## glm(formula = cases ~ density + housing_prop + offset(log(population)), 
##     family = poisson(link = "log"), data = data)
## 
## Coefficients:
##                Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)  -5.482e+00  8.087e-03 -677.955   <2e-16 ***
## density       4.172e-05  4.125e-07  101.141   <2e-16 ***
## housing_prop -1.339e-01  1.339e-02   -9.999   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 27033  on 26  degrees of freedom
## Residual deviance: 17191  on 24  degrees of freedom
## AIC: 17484
## 
## Number of Fisher Scoring iterations: 4

Dapat dilihat terdapat overdispersi karena nilai residual deviance sebesar 17191 dan nilai df sebesar 24, oleh karena itu digunakan binomial negatif

Binomial Negatif

# 4. Konversi koefisien ke PRR (Prevalence Rate Ratio)
PRR <- tidy(model_pois, exponentiate = TRUE, conf.int = TRUE)
PRR

# 5. Cek Overdispersion
disp <- model_pois$deviance / model_pois$df.residual
disp

## [1] 716.3092

# Jika overdispersi → gunakan Negative Binomial
if (disp > 1.5) {
  model_nb <- glm.nb(
    cases ~ density + housing_prop + offset(log(population)),
    data = data
  )
  
  summary(model_nb)
}

## 
## Call:
## glm.nb(formula = cases ~ density + housing_prop + offset(log(population)), 
##     data = data, init.theta = 8.967353605, link = log)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.686e+00  2.551e-01 -22.289   <2e-16 ***
## density       6.687e-05  1.432e-05   4.671    3e-06 ***
## housing_prop  1.558e-01  3.911e-01   0.398     0.69    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(8.9674) family taken to be 1)
## 
##     Null deviance: 48.868  on 26  degrees of freedom
## Residual deviance: 27.527  on 24  degrees of freedom
## AIC: 501.31
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  8.97 
##           Std. Err.:  2.40 
## 
##  2 x log-likelihood:  -493.31

Berdasarkan hasil uji dispersi, diperoleh nilai sebesar 8.967 (>1), yang menunjukkan bahwa data mengalami overdispersi sehingga asumsi model Poisson tidak terpenuhi. Oleh karena itu, digunakan model Regresi Binomial Negatif sebagai alternatif yang lebih sesuai. Model ini memperhitungkan parameter dispersi tambahan untuk mengatasi variansi yang lebih besar dari rataan. Hasil analisis menunjukkan bahwa kepadatan penduduk berpengaruh positif dan signifikan terhadap prevalensi TBC (p = 3e-06), sedangkan rumah layak huni tidak menunjukkan pengaruh yang signifikan (p = 0.69).

Interpretasi

Hasil analisis menunjukkan bahwa wilayah dengan kasus tinggi cenderung berdekatan dengan wilayah lain yang juga memiliki kasus tinggi, terutama di kawasan perkotaan seperti Kota Bandung, Bekasi, dan Depok. Faktor kepadatan penduduk berhubungan positif dengan jumlah kasus TBC, sedangkan persentase rumah layak huni berhubungan negatif. Hal ini menunjukkan bahwa lingkungan padat dan kondisi hunian yang buruk meningkatkan risiko penularan, sementara hunian yang layak berperan protektif. Secara epidemiologis, hasil ini menegaskan bahwa faktor sosial dan lingkungan memiliki peran penting terhadap penyebaran TBC di Jawa Barat. Oleh karena itu, pengendalian TBC perlu difokuskan pada wilayah dengan kepadatan tinggi dan kualitas perumahan yang rendah guna mendukung target Indonesia Bebas TBC 2030.

Kesimpulan dan Saran

Kesimpulan

Penelitian ini menunjukkan bahwa penyebaran kasus Tuberkulosis (TBC) di Provinsi Jawa Barat memiliki pola mengelompok, dengan wilayah berpenduduk padat seperti Kota Bandung, Bekasi, dan Depok memiliki jumlah kasus tertinggi. Hasil analisis juga menunjukkan bahwa kepadatan penduduk berhubungan positif dengan jumlah kasus TBC, sedangkan persentase rumah layak huni berhubungan negatif, yang berarti kondisi hunian yang baik dapat menurunkan risiko penularan.

Secara epidemiologis, temuan ini menegaskan bahwa faktor sosial dan lingkungan berperan penting terhadap penyebaran TBC di Jawa Barat. Upaya pengendalian tidak hanya perlu difokuskan pada aspek medis, tetapi juga pada perbaikan kondisi lingkungan dan sosial masyarakat.

Saran

1. Pemerintah daerah perlu meningkatkan upaya perbaikan kualitas hunian dan pengawasan kesehatan lingkungan di wilayah padat penduduk.

2. Dinas Kesehatan Provinsi Jawa Barat disarankan memperkuat surveilans dan edukasi masyarakat terkait pencegahan penularan TBC.

3. Penelitian selanjutnya dapat menambahkan variabel sosial ekonomi lain seperti tingkat pendidikan atau status pekerjaan untuk memperkaya analisis determinan TBC di tingkat wilayah.

Daftar Pustaka

[1] World Health Organization. (2024). Global Tuberculosis Report 2024: The Second National TB Inventory Study in Indonesia. https://www.who.int/teams/global-programme-on-tuberculosis-and-lung-health/tb-reports/global-tuberculosis-report-2024

[2] Kementerian Kesehatan Republik Indonesia. (2024). Laporan Program Penanggulangan Tuberkulosis Tahun 2023. https://www.tbindonesia.or.id/wp-content/uploads/2024/12/Laporan-Program-Penanggulangan-TBC-2023_Final.pdf

[3] Kementerian Kesehatan. (2025, 11 April). Gerakan Indonesia Akhiri TBC. https://kemkes.go.id/id/indonesias-movement-to-end-tb

[4] Kementerian Kesehatan Republik Indonesia. (2024). Laporan Hasil Studi Inventori Tuberkulosis Indonesia 2023-2024. https://cdn.who.int/media/docs/default-source/searo/indonesia/non-who-publications/laporan-hasil-studi-inventori-tb-indonesia-2023-2024.pdf?sfvrsn=e041d479_5&download=true

[5] Rahayuningrum, I. O., & Sulistyani, S. (2024). DETERMINAN SOSIAL KESEHATAN PENYAKIT TUBERKULOSIS DI INDONESIA. INFOKES: Jurnal Ilmiah Rekam Medis dan Informatika Kesehatan, 14(1), 1–6. https://ojs.udb.ac.id/infokes/article/view/3438

[6] Zeanova, H., Taniwan, P., Putri, R. A., & Faidah, D. Y. (2024). ANALISIS FAKTOR PENYEBAB PENYAKIT TBC DI JAWA BARAT MENGGUNAKAN REGRESI BINOMIAL NEGATIF. ResearchGate. https://www.researchgate.net/publication/387640440_ANALISIS_FAKTOR_PENYEBAB_PENYAKIT_TBC_DI_JAWA_BARAT_MENGGUNAKAN_REGRESI_BINOMIAL_NEGATIF

[7] Kustanto, A. (2025). Socioeconomic Determinants and Tuberculosis Burden: A Panel Data Analysis of Indonesian Provinces [MPRA Paper No. 126466]. Munich Personal RePEc Archive. https://mpra.ub.uni-muenchen.de/126466/

[8] Badan Pusat Statistik. Kepadatan penduduk berdasarkan kabupaten/kota di Jawa Barat [Dataset]. Open Data Jabar. https://opendata.jabarprov.go.id/id/dataset/kepadatan-penduduk-berdasarkan-kabupatenkota-di-jawa-barat

[9] Badan Pusat Statistik. Persentase rumah tangga yang menempati rumah layak huni berdasarkan kabupaten/kota di Jawa Barat [Dataset]. Open Data Jabar. https://opendata.jabarprov.go.id/id/dataset/persentase-rumah-tangga-yang-menempati-rumah-layak-huni-berdasarkan-kabupatenkota-di-jawa-barat

[10] Dinas Kesehatan. Jumlah kasus penyakit tuberkulosis berdasarkan kabupaten/kota di Jawa Barat [Dataset]. Open Data Jabar. https://opendata.jabarprov.go.id/id/dataset/jumlah-kasus-penyakit-tuberkulosis-berdasarkan-kabupatenkota-di-jawa-barat

Lampiran

library(tidyverse)
library(MASS)
library(broom)
library(dplyr)
library(sf)
library(spdep)
library(tmap)
library(geodata)

setwd("C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Epidemiologi")
data <- read.csv("DATA UTS EPIDEMIOLOGI.csv", sep = ";", dec =".")
data

str(data)

## 'data.frame':    27 obs. of  5 variables:
##  $ KABKOT  : chr  "BANDUNG" "BANDUNG BARAT" "BEKASI" "BOGOR" ...
##  $ Populasi: int  3753120 1884190 3273870 5682300 1259230 2584990 2387960 2716950 1914040 2554380 ...
##  $ Y       : int  14012 4569 15458 29110 3477 12197 9587 9516 6092 13474 ...
##  $ X1      : int  2156 1468 2617 1899 789 712 2228 876 922 1335 ...
##  $ X2      : num  54.2 38.9 63.8 46.7 67.2 ...

# Ubah nama kolom agar mudah dipakai
data <- data %>%
  rename(
    kabkot = KABKOT,
    population = Populasi,
    cases = Y,
    density = X1,
    housing = X2
  )

# 2. Hitung Prevalensi TBC per 100.000 penduduk
data <- data %>%
  mutate(
    prevalence_per100k = (cases / population) * 100000,
    housing_prop = housing / 100   # ubah persen jadi proporsi
  )

# Lihat hasil
head(data)

# --- 1. Ambil shapefile Indonesia ---
indo_map <- readRDS("C:/Users/Fizky Firdzansyah/Downloads/Semester 5/Epidemiologi/gadm41_IDN_2_pk.rds") %>%
  st_as_sf()

## Samain nama kab/kota di shapefile
jabar <- indo_map %>%
  filter(NAME_1 == "Jawa Barat") %>%
  filter(!NAME_2 %in% c("Waduk Cirata", "Waduk Jatiluhur")) %>%
  mutate(KABKOT = str_squish(toupper(NAME_2)))



# Samain juga di data CSV (buat jaga-jaga)
data <- data %>%
  mutate(
    KABKOT = str_squish(toupper(kabkot))
  )


# --- 4. Join shapefile dan data ---
data_spatial <- left_join(jabar, data, by = "KABKOT")

# --- 5. Cek kesesuaian nama ---
setdiff(unique(jabar$KABKOT), unique(data$KABKOT))

## character(0)

tmap_mode("plot")

## ℹ tmap modes "plot" - "view"

tm_shape(data_spatial) +
  tm_polygons(
    col = "cases",           # pakai kolom 'cases' bukan 'Y'
    palette = "Reds",        # warna gradasi merah
    title = "Jumlah Kasus TBC"
  ) +
  tm_borders(lwd = 0.7, col = "gray40") +
  tm_layout(
    title = "Sebaran Kasus TBC di Kabupaten/Kota Jawa Barat (2024)",
    legend.outside = TRUE,
    frame = FALSE
  )

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
## "brewer.reds"Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

# --- 5. Hitung Moran’s I Global ---
# Konversi sf → Spatial
map_sp <- as_Spatial(data_spatial)

# Buat matriks tetangga (queen contiguity)
nb <- poly2nb(map_sp)

# Buat bobot spasial (row-standardized)
lw <- nb2listw(nb, style = "W")

# Uji Moran's I global
moran_global <- moran.test(data_spatial$cases, lw)
moran_global

## 
##  Moran I test under randomisation
## 
## data:  data_spatial$cases  
## weights: lw    
## 
## Moran I statistic standard deviate = 2.9347, p-value = 0.00167
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.35148462       -0.04000000        0.01779563

# --- 8. Hitung Local Moran’s I (LISA) ---
local_moran <- localmoran(st_drop_geometry(data_spatial)$cases, lw)
head(local_moran)

##            Ii        E.Ii     Var.Ii        Z.Ii Pr(z != E(Ii))
## 1 -0.03867622 -0.02966099 0.08017623 -0.03183863    0.974600737
## 2 -0.14554706 -0.01949157 0.06556356 -0.49230075    0.622506751
## 3  1.05549883 -0.05775959 1.41500884  0.93587190    0.349339144
## 4  1.80599497 -0.04838894 0.36582139  3.06595023    0.002169794
## 5  1.47678414 -0.44997849 0.56976067  2.55259539    0.010692363
## 6  0.70568331 -0.03086110 0.12960434  2.04592245    0.040764001

# --- 9. Klasifikasi LISA (HH, LL, HL, LH) ---

# 1️⃣ Ambil nilai kasus dan rata-ratanya
cases <- st_drop_geometry(data_spatial)$cases
mean_val <- mean(cases, na.rm = TRUE)

# 2️⃣ Hitung nilai lag spasial (rata-rata tetangga)
lag_cases <- lag.listw(lw, cases)

# 3️⃣ Buat klasifikasi berdasarkan kuadran LISA
lisa_cluster <- ifelse(cases > mean_val & lag_cases > mean(lag_cases), "High-High",
                ifelse(cases < mean_val & lag_cases < mean(lag_cases), "Low-Low",
                ifelse(cases > mean_val & lag_cases < mean(lag_cases), "High-Low",
                ifelse(cases < mean_val & lag_cases > mean(lag_cases), "Low-High", "Undefined"))))

# 4️⃣ Tambahkan ke data spasial
data_spatial$lisa_cluster <- lisa_cluster

# 5️⃣ Cek distribusi klaster
table(data_spatial$lisa_cluster)

## 
## High-High  High-Low  Low-High   Low-Low 
##         7         4         4        11

# 6️⃣ Visualisasikan peta LISA
tm_shape(data_spatial) +
  tm_polygons(
    col = "lisa_cluster",
    palette = c("red", "blue", "pink", "lightblue", "gray"),
    title = "Klaster LISA (Local Moran’s I)"
  ) +
  tm_borders() +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE
  )

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values') to fill.scale
## = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'

# --- 1. Pastikan library corrplot sudah di-load ---
library(corrplot)

## Warning: package 'corrplot' was built under R version 4.3.3

## corrplot 0.95 loaded

# --- 2. Pilih variabel numerik ---
data_num <- data %>% 
  dplyr::select(cases, density, housing)

# --- 3. Hitung matriks korelasi ---
cor_matrix <- cor(data_num, method = "pearson", use = "complete.obs")

# --- 4. Plot hasil korelasi ---
corrplot(
  cor_matrix, 
  method = "color",       # bisa juga "circle", "number", dsb.
  type = "upper",         # hanya tampilkan segitiga atas
  order = "hclust",       # urutkan berdasarkan clustering
  addCoef.col = "black",  # tampilkan nilai korelasi
  tl.col = "black",       # warna label
  tl.srt = 45,            # rotasi label
  col = colorRampPalette(c("blue", "white", "red"))(200)  # gradasi warna
)

# 3. Model Cross-Sectional → Prevalence Rate Ratio (PRR) via Poisson Regression
model_pois <- glm(
  cases ~ density + housing_prop + offset(log(population)),
  family = poisson(link = "log"),
  data = data
)

summary(model_pois)

## 
## Call:
## glm(formula = cases ~ density + housing_prop + offset(log(population)), 
##     family = poisson(link = "log"), data = data)
## 
## Coefficients:
##                Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)  -5.482e+00  8.087e-03 -677.955   <2e-16 ***
## density       4.172e-05  4.125e-07  101.141   <2e-16 ***
## housing_prop -1.339e-01  1.339e-02   -9.999   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 27033  on 26  degrees of freedom
## Residual deviance: 17191  on 24  degrees of freedom
## AIC: 17484
## 
## Number of Fisher Scoring iterations: 4

# 4. Konversi koefisien ke PRR (Prevalence Rate Ratio)
PRR <- tidy(model_pois, exponentiate = TRUE, conf.int = TRUE)
PRR

# 5. Cek Overdispersion
disp <- model_pois$deviance / model_pois$df.residual
disp

## [1] 716.3092

# Jika overdispersi → gunakan Negative Binomial
if (disp > 1.5) {
  model_nb <- glm.nb(
    cases ~ density + housing_prop + offset(log(population)),
    data = data
  )
  
  summary(model_nb)
}

## 
## Call:
## glm.nb(formula = cases ~ density + housing_prop + offset(log(population)), 
##     data = data, init.theta = 8.967353605, link = log)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -5.686e+00  2.551e-01 -22.289   <2e-16 ***
## density       6.687e-05  1.432e-05   4.671    3e-06 ***
## housing_prop  1.558e-01  3.911e-01   0.398     0.69    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(8.9674) family taken to be 1)
## 
##     Null deviance: 48.868  on 26  degrees of freedom
## Residual deviance: 27.527  on 24  degrees of freedom
## AIC: 501.31
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  8.97 
##           Std. Err.:  2.40 
## 
##  2 x log-likelihood:  -493.31

PRR <- tidy(model_nb, exponentiate = TRUE, conf.int = TRUE)
PRR

# 6. Tampilkan tabel prevalensi tertinggi
dplyr::select(data, kabkot, prevalence_per100k) %>%
  arrange(desc(prevalence_per100k)) %>%
  head(10)