LINK DASHBOARD https://chrisrafael.shinyapps.io/Dashboard_Epidemiologi/

Abstrak

Penyakit Tuberkulosis (TBC), merupakan sebuah penyakit yang disebabkan oleh bakteri Mycobacterium tuberculosis dan menyerangparu-paru seorang individu. TBC merupakan salah satu tantangan serius yang dihadapi sistem kesehatan di seluruh dunia, terutama di Indonesia dimana kasus TBC sedang menjadi masalah besar di wilayah Asia Tenggara. Pada tingkat negara, didapatkan bahwa provinsi Jawa Barat merupakan provinsi dengan jumlah kasus TBC tertinggi di seluruh Indonesia. Karena itu, dalam rangka membantu pemerintah pusat serta lokal dalam mengidentifikasi serta mencegah penyebaran TBC akan dilakukan analisis epidemiologis dalam bentuk Insidensi Rate, mengidentifikasi autokorelasi spasial pada insidensi rate, serta melihat tren insidensi rate pada Jawa Barat dan juga dilakukan case-control study untuk melihat faktor-faktor pada Individu yang dapat memudahkan terjangkitnya TBC. Ddiapatkan bahwa Insidensi Rate TBC di Jawa Barat terkonsentrasi pada titik-titik perkotaan seperti Kota Cirebon, Kota Bogor serta didapatkan trend Insidensi Rate yang meningkat dari 2021-2024. Didapatkan bahwa tidak terdapat autokorelasi cukup signifikan pada Insidensi Rate. didapatkan bahwa faktor individu seperti status merokok serta riwayat diabetes mempunyai peran dalan menentukan apakah sebuah individu akan terjangkit atau tidak.
Kata Kunci: TBC, Jawa Barat, Epidemiologi

1 Pendahuluan

1.1 Latar Belakang

TBC, merupakan penyakit penyebar yang disebabkan oleh bakteri Mycobacterium tuberculosis. Melalui udara, bakteria bisa ditransmisi dari host ke host dimana akan memasuki tubuh dan menjadi dorman. Dimana sebuah bakteria memiliki peluang 5-10% akan reaktivisasi setelah dorman didalam tubuh. Mycobacterium tuberculosis akan menyerang paru-paru, namun bisa menyerang organ lainnya. Dampaknya pada seorang individu adalah kepucatan pada kulit, batuk berdarah, demam, infeksi, penurunan berat badan,serta juga kematian (Fogel, 2015).

Tuberkulosis (TBC) masih menjadi salah satu tantangan kesehatan masyarakat paling serius di dunia. Walaupun berbagai kemajuan telah dicapai dalam hal penanganan TBC, World Health Organization (WHO) tetap mencatat TBC sebagai salah satu penyebab utama kematian akibat penyakit menular, dengan jutaan kasus baru muncul setiap tahunnya (WHO, 2024). Pada 2021, tercatat bahwa 45% dari kasus TBC mayoritas dari area Asia Tenggara, diikuti oeh Afrika sebesar 23%. Dimana pada 2022, sekitar 2 mmiliyar atau sekitar 25% dari populasi dunia sedang menderita dari TBC. (Alsayed dan Gunosewoyo, 2023). Penyebaran penyakit ini pun tidak merata, di mana sebagian besar kasus terkonsentrasi di negara-negara berkembang, termasuk Indonesia.

Di tingkat nasional, Indonesia memiliki tantangan TBC yang sangat besar. Selama beberapa tahun terakhir, Indonesia secara konsisten berada pada posisi kedua atau ketiga dengan jumlah kasus TBC tertinggi di dunia (Kementerian Kesehatan RI, 2024). Kondisi ini menjadikan pengendalian TBC sebagai prioritas utama dalam kebijakan kesehatan nasional, sebagaimana tercantum dalam berbagai strategi dan rencana aksi untuk mencapai target eliminasi TBC pada tahun 2030.

Dalam kerangka tersebut, Provinsi Jawa Barat menjadi wilayah dengan perhatian khusus. Sebagai provinsi berpenduduk terbanyak di Indonesia—lebih dari 49 juta jiwa (BPS Provinsi Jawa Barat, 2023)—Jawa Barat menanggung konsekuensi demografis yang besar terhadap tingginya kasus TBC. Data kesehatan menunjukkan bahwa provinsi ini secara konsisten menjadi penyumbang utama jumlah kasus TBC nasional (Dinas Kesehatan Provinsi Jawa Barat, 2023). Kepadatan penduduk yang tinggi, pesatnya urbanisasi, serta berbagai tantangan sosial dan ekonomi di wilayah perkotaan maupun pedesaan menjadikan pengendalian TBC di Jawa Barat sebagai tugas yang kompleks dan mendesak.

1.2 Rumusan Masalah

Berdasarkan latar tersebut, beberapa masalah utama yang diidentifikasi adalah berapa angka prevalensi kasus tuberkulosis (TBC) di Provinsi Jawa Barat pada tahun 2021-2024? Bagaimana distribusi prevalensi tuberkulosis (TBC) berdasarkan Kabupaten/Kota di Provinsi Jawa Barat? Apakah faktor risiko biologis dalam bentuk jenis kelamin serta faktor seperti status merokok dan riwayat diabetes memiliki pengaruh signifikan?

1.3 Tujuan Penelitian

Penelitian ini bertujuan untuk mendeskripsikan jumlah kasus TBC total di Jawa Barat serta prevalensinya pada populasi Jawa Barat sebagai bentuk dukungan data pada rspon pemerintah serta institusi kesehatan pada kasus TBC. Serta dilakukan riset case-control untuk menunjukkan adanya pengaruh faktor risiko yang membuat beberapa individu lebih rawan terhadap terjangkit dalam membantu antisipasi dan pencegahan TBC.

2 Tinjauan Pustaka

2.1 Epidemiologic Triad (Agent-Host-Environment)

Untuk memahami dinamika penularan penyakit menular seperti tuberkulosis (TBC), kerangka kerja konseptual yang paling fundamental adalah model Segitiga Epidemiologi atau Epidemiologic Triad (Gordis, 2018). Model ini bersifat dinamis yang kejadian dan distribusi penyakit dalam populasi bukanlah disebabkan oleh satu faktor tunggal, melainkan hasil dari interaksi kompleks dan berkelanjutan antara tiga komponen utama: Agent, Host, dan Environment. Kejadian penyakit (TBC) terjadi ketika keseimbangan antara ketiga faktor ini terganggu. Misalnya, kerentanan pejamu meningkat (misal, akibat gizi buruk) atau lingkungan berubah menjadi lebih mendukung transmisi (misal, urbanisasi padat). Dalam konteks TBC secara keseluruhan di Jawa Barat, ketiga elemen ini berinteraksi secara konstan.

2.1.1 Agent (Agen)

Agen merujuk pada mikroorganisme penyebab penyakit. Dalam kasus TBC, agennya adalah bakteri Mycobacterium tuberculosis. Dalam konteks TBC keseluruhan, karakteristik agen yang relevan meliputi:

  • Infektivitas: Kemampuan agen untuk masuk dan berkembang biak dalam pejamu.

  • Patogenisitas: Kemampuan agen untuk menimbulkan penyakit klinis (TBC aktif) setelah berhasil menginfeksi. Penting diingat bahwa tidak semua orang yang terinfeksi M. tuberculosis akan jatuh sakit (menjadi TBC Laten).

  • Virulensi: Tingkat keparahan atau keganasan penyakit yang disebabkan oleh agen.

  • Transmisibilitas: Kemudahan agen untuk ditularkan. M. tuberculosis utamanya ditularkan melalui airborne (via droplet nuclei), yang memungkinkannya bertahan lama di udara dalam ruangan tertutup dan terhirup oleh pejamu baru (Kementerian Kesehatan RI, 2023).

2.1.2 Host (Pejamu)

Host merujuk pada manusia yang rentan (suspektibel) terhadap infeksi. Faktor-faktor intrinsik pada pejamu sangat menentukan apakah paparan terhadap agen akan menyebabkan infeksi (TBC laten) dan apakah infeksi tersebut akan berkembang menjadi penyakit aktif (TBC aktif). Faktor pejamu yang krusial untuk TBC meliputi:

  • Faktor Imunologi: Ini adalah faktor terpenting. Sistem kekebalan tubuh yang lemah (imunokompromais) akibat kondisi seperti koinfeksi HIV, diabetes melitus, atau malnutrisi (gizi buruk), secara drastis meningkatkan risiko transisi dari TBC laten ke TBC aktif (WHO, 2024).

  • Faktor Perilaku: Perilaku yang meningkatkan paparan atau menurunkan pertahanan tubuh. Contoh utamanya adalah kebiasaan merokok, yang merusak mekanisme pertahanan paru-paru dan meningkatkan risiko TBC aktif.

  • Faktor Demografis: Usia (kelompok usia sangat muda dan lanjut usia lebih rentan) dan riwayat kontak erat sebelumnya dengan penderita TBC.

  • Faktor Genetik: Adanya predisposisi genetik tertentu yang mungkin memengaruhi respons imun individu terhadap bakteri.

2.1.3 Environment (Lingkungan)

Environment merujuk pada kondisi ekstrinsik di luar pejamu yang mempengaruhi paparan dan kerentanan. Lingkungan bertindak sebagai mediator yang memfasilitasi atau menghambat transmisi agen.

  • Lingkungan Fisik: Ini sangat sentral dalam penularan TBC. Meliputi:
    1. Kepadatan hunian, di mana tinggal di rumah yang sempit dan penuh sesak (overcrowding) meningkatkan intensitas dan durasi paparan;
      1. Kualitas ventilasi udara yang buruk, yang menyebabkan droplet nuclei terakumulasi di dalam ruangan; dan (3) Kurangnya paparan sinar matahari langsung (sinar UV dapat membunuh M. tuberculosis).
  • Lingkungan Sosial-Ekonomi: Faktor ini seringkali menjadi akar masalah (root cause) dari kerentanan. Kemiskinan berkorelasi langsung dengan malnutrisi (menurunkan status imun host) dan kondisi hunian yang buruk (faktor environment). Tingkat pendidikan dan pengetahuan yang rendah dapat menghambat perilaku pencarian layanan kesehatan, sementara stigma sosial dapat menyebabkan penderita menunda pengobatan (Dinas Kesehatan Provinsi Jawa Barat, 2023).

Interaksi ketiga komponen inilah yang menentukan tingkat endemisitas TBC di suatu wilayah.

2.2 Desain Studi Epidemiologi

Desain studi adalah cetak biru atau arsitektur penelitian yang memandu peneliti dalam mengumpulkan dan menganalisis data untuk menjawab pertanyaan penelitian (Budiarto, 2012). Secara garis besar, desain studi dibagi menjadi deskriptif (mendeskripsikan distribusi penyakit) dan analitik (menguji hipotesis hubungan sebab-akibat).Desain analitik observasional yang umum digunakan meliputi:

2.2.1 Studi Cross-sectional

Dalam studi cross-sectional, pengukuran variabel paparan (faktor risiko) dan variabel outcome (penyakit/TBC) dilakukan secara bersamaan (simultaneously) pada satu titik waktu (Rothman, 2012). Desain ini dapat diibaratkan sebagai “potret” atau snapshot dari populasi pada saat itu.

Kelebihan:

  • Efisien: Relatif cepat untuk dilaksanakan dan biaya lebih murah dibandingkan desain analitik lainnya.

  • Mengukur Prevalensi: Merupakan desain terbaik untuk mengukur frekuensi penyakit (prevalensi) di populasi, sehingga cocok untuk menggambarkan beban penyakit (burden of disease).

  • Hipoetetis Awal: Baik untuk mengeksplorasi asosiasi dan menghasilkan hipotesis awal untuk penelitian kohort atau eksperimen lebih lanjut.

Kekurangan:

  • Masalah Temporal : Kekurangan utamanya adalah ketidakmampuan menentukan mana yang lebih dulu terjadi (temporalitas). Sulit dipastikan apakah paparan menyebabkan penyakit atau penyakit yang menyebabkan paparan (misalnya, apakah status gizi buruk menyebabkan TBC, atau TBC aktif menyebabkan gizi buruk?).

  • Tidak Mengukur Insidensi: Tidak dapat digunakan untuk mengukur laju kasus baru (insidensi), sehingga tidak dapat menghitung Risk Ratio (RR).

  • Bias Prevalensi: Dapat terjadi bias jika paparan (faktor risiko) berhubungan dengan durasi penyakit.

2.2.2 Studi Kohort (Cohort)

Desain ini mengikuti sekelompok individu (kohort) yang bebas dari penyakit pada awal studi, yang dibagi berdasarkan status paparan (terpapar dan tidak terpapar). Kohort tersebut kemudian diikuti dari waktu ke waktu (follow-up) untuk melihat siapa yang berkembang menjadi sakit.

Kelebihan:

  • Temporalitas Jelas: Dapat menentukan urutan waktu (paparan terjadi sebelum penyakit), sehingga menjadi desain observasional terkuat untuk menentukan penyebab (causality).

  • Mengukur Insidensi & RR: Dapat menghitung angka insidensi (kasus baru) secara langsung dan menghitung Risk Ratio (RR).

  • Paparan Langka: Baik untuk meneliti paparan yang langka.

Kekurangan:

  • Mahal dan Lama: Membutuhkan biaya sangat besar dan waktu yang sangat lama, terutama untuk penyakit dengan masa laten panjang seperti TBC.

  • Loss to Follow-up: Risiko subjek penelitian hilang atau keluar dari studi selama periode follow-up, yang dapat menyebabkan bias.

  • Penyakit Langka: Tidak efisien untuk meneliti penyakit yang langka.

2.2.3 Studi Kasus-Kontrol (Case-Control)

Desain ini dimulai dengan mengidentifikasi kelompok kasus (individu yang sakit TBC) dan kelompok kontrol (individu yang tidak sakit TBC). Peneliti kemudian melihat ke belakang (retrospektif) untuk membandingkan riwayat paparan faktor risiko antara kedua kelompok tersebut.

Kelebihan:

  • Efisien untuk Penyakit Langka: Sangat cocok untuk meneliti penyakit yang jarang terjadi atau memiliki masa laten panjang (seperti TBC).

  • Cepat dan Murah: Relatif cepat dan murah dibandingkan kohort, karena tidak memerlukan follow-up.

  • Multiple Paparan: Dapat meneliti berbagai kemungkinan faktor risiko (paparan) secara bersamaan.

Kekurangan:

  • Bias Recall (Daya Ingat): Sangat rentan terhadap bias daya ingat, karena subjek (terutama kasus) mungkin mengingat paparan masa lalu secara berbeda dibandingkan kelompok kontrol.

  • Bias Seleksi: Sulit untuk memilih kelompok kontrol yang benar-benar sebanding (representatif) dengan populasi asal kasus.

  • Tidak Mengukur Prevalensi/Insidensi: Tidak dapat menghitung prevalensi atau insidensi. Hanya menghasilkan Odds Ratio (OR).

2.3 Ukuran Asosiasi dan Frekuensi

Dalam epidemiologi, ukuran digunakan untuk mengkuantifikasi frekuensi penyakit dan kekuatan hubungan antara paparan dan penyakit. Sering kali digunakan untuk melihat adanya signifikansi faktor diantara variabel-variabel saat melakukan sebuah riset seperti case-control study. Sering kali digunakan analisis menggunakan metode Analisis Data Kategorik.

3 Metodologi Penelitian

3.1 Sumber Data

Data yang digunakan dalam penelitian ini merupakan data sekunder yang diperoleh dari sumber resmi yaitu Website Dinas Kesehatan Provinsi Jawa Barat dan Badan Pusat Statistik (BPS).

3.2 Variabel

Variabel yang digunakan adalah jumlah kasus Tuberkulosis (TBC) di setiap kabupaten/kota di Provinsi Jawa Barat dan jumlah penduduk di setiap Kabupaten/Kota di Provinsi Jawa Barat pada tahun 2021-2024.

3.3 Metode Analisis

Analisis yang digunakan terdiri atas perhitungan prevalensi dan autokorelasi spasial. Analisis deskriptif dilakukan untuk memvisualisasikan persebaran kasus TBC di Jawa Barat. Autokorelasi spasial dihitung menggunakan Global Moran’s I dan Geary’s C untuk mendeteksi pola sebaran secara keseluruhan dan Local Moran’s I (LISA) untuk mengidentifikasi klaster signifikan antarwilayah.

3.3.1 Prevalensi

Prevalensi adalah ukuran proporsi individu dalam suatu populasi yang memiliki penyakit (kasus existing, baik baru maupun lama) pada satu titik waktu tertentu (point prevalence) atau selama periode waktu tertentu (period prevalence).

\[ \text{Prevalensi} = \frac{\text{Jumlah kasus (baru + lama) pada waktu } t}{\text{Total populasi pada waktu } t} \]

Ukuran ini menggambarkan beban penyakit (burden of disease) di masyarakat.

3.3.2 Insidensi

Insidensi mengukur jumlah kasus baru suatu penyakit yang terjadi dalam populasi berisiko selama periode waktu tertentu. Insidensi mencerminkan risiko seseorang untuk terkena penyakit. Ini biasanya diukur dalam studi kohort. Rumus untuk Cumulative Incidence (CI) atau Risk adalah:

\[ \text{Insidensi Kumulatif (Risk)} = \frac{\text{Jumlah kasus baru selama periode waktu tertentu}} {\text{Jumlah populasi berisiko}} \]

3.3.3 Attack rate (AR)

Attack Rate adalah ukuran frekuensi yang digunakan dalam epidemiologi untuk menggambarkan proporsi individu dalam suatu populasi yang menjadi sakit (terinfeksi) selama periode waktu tertentu.

\[ AR = \frac{\text{Jumlah kasus baru selama periode wabah}} {\text{Jumlah populasi yang berisiko}} \times 100\% \]

3.3.4 Odds Ratio (OR)

OR adalah ukuran asosiasi yang digunakan dalam studi kasus-kontrol dan cross-sectional. OR membandingkan odds (peluang) terjadinya penyakit pada kelompok terpapar dengan odds terjadinya penyakit pada kelompok tidak terpapar.

  • OR = 1 → tidak ada asosiasi antara paparan dan penyakit.

  • OR > 1 → paparan berhubungan dengan peningkatan odds penyakit (faktor risiko).

  • OR < 1 → paparan berhubungan dengan penurunan odds penyakit (faktor protektif).

3.3.5 Risk Ratio (RR) atau Relative Risk

RR adalah ukuran asosiasi yang digunakan dalam studi kohort. RR membandingkan risiko (insidensi) penyakit pada kelompok terpapar dengan risiko (insidensi) penyakit pada kelompok tidak terpapar.

\[ RR = \frac{\text{Insidensi pada kelompok terpapar}} {\text{Insidensi pada kelompok tidak terpapar}} = \frac{a / (a + b)}{c / (c + d)} \]

3.3.6 Attributable Risk (AR)

AR (atau Selisih Risiko) adalah ukuran yang menunjukkan kelebihan risiko penyakit pada kelompok terpapar yang dapat diatribusikan (disebabkan) oleh paparan tersebut. Ukuran ini juga memerlukan data insidensi dari studi kohort.

\[ AR = (\text{Insidensi terpapar}) - (\text{Insidensi tidak terpapar}) \]

3.4 Alur Kerja Penelitian

Penelitian dimulai dengan pengumpulan data-data sekunder, yaitu populasi total kabupaten serta kota di Jawa Barat dari 2021-2024 kemudian total kasus penyakit TBC yang akan diikuti oleh analisis deskriptif berdasarkan epidemiologi. Analisis akan dilakukan pada program lunak Rstudio. Pertama akan dilakukan plottingan peta menggunakan jumlah total kasus TBC kemudian melakukan analisis deskriptif. Kemudian menghitung Prevalensi Rate berdasarkan populasi pada semua kabupaten/kota pada interval 2021-2024 dimana kemudian akan dilakukan plottingan pada peta kemudian melihat trend data pada semua variabel dari 2021-2024 dan kemudian analisis spasial autokorelasi pada prevalensi rate. Kemudian akan dilakukan eksperimen berdasarkan faktor risiko host.

4 Hasil dan Pembahasan

4.1 Tabel Hasil

Berdasarkan perhitungan yang dilakukan didapatkan hasil tabel berikut yang mencangkup data populasi dalam interval 4 tahun serta jumlah total kasus TBC pada masing-masing kabupaten. Dilakukan perhitungan Prevalensi Rate pada masing-masing kabupaten serta kota untuk melihat tingkat keparahan atau beban TBC yang diberikan pada populasi di masing-masing kabupaten/kota dimana didapatkan tabel sebagai berikut:

# Package
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')`
## Loading required package: 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(sp)
## Warning: package 'sp' was built under R version 4.3.3
library(sf)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(spData)
library(geodata)
## Loading required package: terra
## terra 1.7.46
library(skimr)
library(psych)
## Warning: package 'psych' was built under R version 4.3.3
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:terra':
## 
##     describe, distance, rescale
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.3.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(grid)
## 
## Attaching package: 'grid'
## The following object is masked from 'package:terra':
## 
##     depth
library(viridis) 
## Warning: package 'viridis' was built under R version 4.3.3
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.3.3
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.3.3
## Loading required package: Matrix
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:psych':
## 
##     logit
## The following object is masked from 'package:dplyr':
## 
##     recode
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:terra':
## 
##     time<-
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dplyr)

# Set up Data
setwd("C:\\Users\\HP\\Documents\\KULIAH\\SEMESTER 5\\Epidem")
indo_adm2 <- st_read("IDN_AdminBoundaries_candidate.gdb",
                     layer = "idn_admbnda_adm2_bps_20200401")
## Reading layer `idn_admbnda_adm2_bps_20200401' from data source 
##   `C:\Users\HP\Documents\KULIAH\SEMESTER 5\Epidem\IDN_AdminBoundaries_candidate.gdb' 
##   using driver `OpenFileGDB'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received a polygon with more than 100 parts.  The
## processing may be really slow.  You can skip the processing by setting
## METHOD=SKIP.
## Simple feature collection with 522 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
## Geodetic CRS:  WGS 84
jabar <- indo_adm2[indo_adm2$admin1Name_en == "Jawa Barat", ]
jabar <- jabar[jabar$admin2Name_en != "Waduk Cirata", ]
data <- (read.csv("ceremony.csv", sep=",", dec="."))
jabar$id <- c(1:27)
jabar_sf <- st_as_sf(jabar)
data$id <- c(1:27)
jabar_merged <- jabar_sf %>%
  left_join(data, by = "id")

# Prevalence Rate
jabar_merged <- jabar_merged %>%
  mutate(
    PR_2021 = (Kasus.TBC.2021 / Populasi.2021) * 100000,
    PR_2022 = (Kasus.TBC.2022 / Populasi.2022) * 100000,
    PR_2023 = (Kasus.TBC.2023 / Populasi.2023) * 100000,
    PR_2024 = (Kasus.TBC.2024 / Populasi.2024) * 100000
  )

# Tabel Hasil
df_tbc <- jabar_merged %>%
  st_drop_geometry() %>%  
  select(
    Kabupaten.Kota,
    Populasi.2021, Populasi.2022, Populasi.2023, Populasi.2024,
    Kasus.TBC.2021, Kasus.TBC.2022, Kasus.TBC.2023, Kasus.TBC.2024, PR_2021, PR_2022, PR_2023, PR_2024
  )
print(df_tbc)
##      Kabupaten.Kota Populasi.2021 Populasi.2022 Populasi.2023 Populasi.2024
## 1           BANDUNG       3652400       3687250       3721110       3753120
## 2     BANDUNG BARAT       1808420       1834230       1859640       1884190
## 3            BEKASI       3148740       3193840       3237420       3273870
## 4             BOGOR       5484150       5556310       5627020       5682300
## 5            CIAMIS       1234830       1243320       1251540       1259230
## 6           CIANJUR       2500640       2529810       2558140       2584990
## 7           CIREBON       2301330       2331360       2360440       2387960
## 8             GARUT       2613530       2648950       2683670       2716950
## 9         INDRAMAYU       1851730       1873400       1894330       1914040
## 10         KARAWANG       2465570       2496190       2526000       2554380
## 11     KOTA BANDUNG       2461410       2484150       2506600       2528160
## 12      KOTA BANJAR        202720        205140        207510        209790
## 13      KOTA BEKASI       2568020       2598070       2627210       2644060
## 14       KOTA BOGOR       1050920       1060940       1070720       1078350
## 15      KOTA CIMAHI        574450        582650        590780        598700
## 16     KOTA CIREBON        335810        338940        341980        344850
## 17       KOTA DEPOK       2081130       2113620       2145400       2163640
## 18    KOTA SUKABUMI        350150        355420        360640        365740
## 19 KOTA TASIKMALAYA        723100        732480        741760        750730
## 20         KUNINGAN       1175950       1189010       1201760       1213930
## 21       MAJALENGKA       1315010       1328010       1340620       1352540
## 22      PANGANDARAN        425590        428600        431460        434100
## 23       PURWAKARTA       1008930       1023180       1037070       1050340
## 24           SUBANG       1620700       1635560       1649820       1663160
## 25         SUKABUMI       2747450       2775310       2802400       2828020
## 26         SUMEDANG       1159260       1168840       1178240       1187130
## 27      TASIKMALAYA       1876890       1892220       1907050       1920920
##    Kasus.TBC.2021 Kasus.TBC.2022 Kasus.TBC.2023 Kasus.TBC.2024   PR_2021
## 1            5708          10650          12697          14012 156.28080
## 2            1736           3043           4178           4569  95.99540
## 3            4798           8379          13515          15458 152.37841
## 4           11946          21516          27690          29110 217.82774
## 5            1606           2946           3042           3477 130.05839
## 6            4673           7059           9327          12197 186.87216
## 7            3383           6981           8238           9587 147.00195
## 8            4786           7920           8615           9516 183.12397
## 9            1701           2729           5148           6092  91.86004
## 10           4528           8020          12868          13474 183.64922
## 11           8919          14541          18327          18627 362.35329
## 12            269            742           1175           1519 132.69534
## 13           6051          10747          14279          13640 235.62901
## 14           1440           7740          10354          11418 137.02280
## 15           1782           4321           4602           4513 310.20977
## 16           1910           2760           4127           4122 568.77401
## 17           4042           6549           7849           8422 194.22141
## 18           1223           2383           3431           3477 349.27888
## 19           1476           2838           4749           4693 204.12115
## 20           1651           2441           3705           4045 140.39713
## 21           1723           3162           4318           4517 131.02562
## 22            410            734            954            968  96.33685
## 23           2410           4589           5607           4687 238.86692
## 24           2912           4769           5666           6321 179.67545
## 25           4805           7491          10950          12411 174.88944
## 26           1372           2501           3020           3930 118.35136
## 27           1995           3110           3528           4881 106.29286
##     PR_2022   PR_2023   PR_2024
## 1  288.8331  341.2154  373.3427
## 2  165.9007  224.6671  242.4915
## 3  262.3488  417.4621  472.1629
## 4  387.2354  492.0900  512.2926
## 5  236.9462  243.0605  276.1211
## 6  279.0328  364.6008  471.8393
## 7  299.4390  349.0027  401.4724
## 8  298.9864  321.0156  350.2457
## 9  145.6710  271.7584  318.2797
## 10 321.2896  509.4220  527.4861
## 11 585.3511  731.1498  736.7809
## 12 361.7042  566.2378  724.0574
## 13 413.6532  543.5043  515.8733
## 14 729.5417  967.0129 1058.8399
## 15 741.6116  778.9702  753.7999
## 16 814.3034 1206.7957 1195.3023
## 17 309.8476  365.8525  389.2514
## 18 670.4744  951.3642  950.6753
## 19 387.4509  640.2340  625.1249
## 20 205.2968  308.2978  333.2153
## 21 238.1006  322.0898  333.9642
## 22 171.2552  221.1097  222.9901
## 23 448.5037  540.6578  446.2365
## 24 291.5821  343.4314  380.0596
## 25 269.9158  390.7365  438.8583
## 26 213.9728  256.3145  331.0505
## 27 164.3572  184.9978  254.0970

4.2 Peta Deskriptif

Pertama dilakukan pemlottingan peta berdasarkan dari jumlah total kasus TBC yang terdapat pada masing-masing kabupaten/kota yang ada di Jawa Barat. Dimana akan dilakukan analisis data pada 4 tahun terakhir dari 2021 sampai 2024. Berikut didapatkan peta Total Kasus TBC di Jawa Barat pada 2021.

library(scales)
## Warning: package 'scales' was built under R version 4.3.3
## 
## Attaching package: 'scales'
## The following object is masked from 'package:viridis':
## 
##     viridis_pal
## The following objects are masked from 'package:psych':
## 
##     alpha, rescale
## The following object is masked from 'package:terra':
## 
##     rescale
# Fungsi bikin label interval (benar)
interval_labels <- function(breaks) {
  labs <- c()
  for (i in seq_along(breaks)[-length(breaks)]) {
    labs[i] <- paste0(
      format(breaks[i], big.mark = ","),
      " – ",
      format(breaks[i + 1], big.mark = ",")
    )
  }
  labs
}

# Tentukan breaks (misalnya per 5.000)
max_val <- max(jabar_merged$Kasus.TBC.2024, na.rm = TRUE)
breaks_interval <- seq(0, ceiling(max_val / 5000) * 5000, by = 5000)

# Ambil label dengan satu lebih sedikit dari breaks
label_interval <- interval_labels(breaks_interval)

# === Fungsi plotting untuk tiap tahun ===
plot_tbc_interval <- function(data, var, year) {
  ggplot(data) +
    geom_sf(aes(fill = .data[[var]]),
            color = "white", linewidth = 0.3) +
    scale_fill_viridis_c(
      option = "plasma",
      direction = -1,
      name = "Jumlah Kasus TBC",
      breaks = breaks_interval[-1], # buang elemen pertama (biar jumlah sama)
      labels = label_interval
    ) +
    labs(
      title = paste("Jumlah Kasus TBC di Jawa Barat", year),
      caption = "Sumber: Open Data Jabar"
    ) +
    theme_minimal(base_size = 14) +
    theme(
      panel.background = element_rect(fill = "white", color = NA),
      panel.grid.major = element_line(color = "gray90"),
      panel.grid.minor = element_line(color = "gray95"),
      plot.title = element_text(face = "bold", hjust = 0.5),
      legend.title = element_text(face = "bold"),
      legend.text = element_text(size = 9)
    )
}

# === Jalankan untuk semua tahun ===
plot_tbc_interval(jabar_merged, "Kasus.TBC.2021", 2021)

plot_tbc_interval(jabar_merged, "Kasus.TBC.2022", 2022)

plot_tbc_interval(jabar_merged, "Kasus.TBC.2023", 2023)

plot_tbc_interval(jabar_merged, "Kasus.TBC.2024", 2024)

Bisa dilihat tingginya konsentrasi kasus TBC pada Kabupaten Bogor serta pada Kabupaten Bandung dan juga di Kota Bandung. Serta dilihat bahwa ada pola pada persebaran pada total kasus TBC di Jawa Barat pada tahun 2021, dimana Jawa Barat pada bagian barat sering kali memiliki jumlah kasus yang lebih tinggi dari bagian timur. Namun, hal tersebut bisa dijelaskan dengan adanya homogenitas populasi pada masing-masing kabupaten/kota, sehingga dibutuhkan perhitungan lebih menjelaskan. Tahun 2022, juga kurang lebih sama dengan tahun 2021 dengan perbedaan yang minim seperti di kabupaten Tasikmalaya. Untuk data jumlah kasus 2023 terlihat juga perubahan minim seperti peningkatan kasus pada Kabupaten Bekasi, Kabupaten Karawang, dan Kota Bekasi. Untuk 2024, terlihat jelas peningkatan kasus TBC pada mayoritas Kabupaten/Kota di Jawa Barat dengan Kabupaten Ciamis sebagai pengecualian dimana kasusnya berkurang. Namun berdasarkan dari angka yang diberikan bisa dilihat sebuah peningkatan signifikan dari tahun ke tahun di semua kabupaten serta kota. Untuk melihat lebih dalam akan dilakukan analisis deskriptif pada masing-masing tahun mengenai jumlah total kasus TBC:

library(psych)
library(dplyr)

describe(
  jabar_merged %>%
    st_drop_geometry() %>%         # buang kolom geometry dulu
    select(Kasus.TBC.2021:Kasus.TBC.2024)
)
##                vars  n    mean      sd median trimmed     mad min   max range
## Kasus.TBC.2021    1 27 3305.74 2650.67   1995 2943.96 1359.54 269 11946 11677
## Kasus.TBC.2022    2 27 5950.41 4575.08   4589 5353.39 3270.62 734 21516 20782
## Kasus.TBC.2023    3 27 7850.33 5933.35   5607 7122.30 3835.49 954 27690 26736
## Kasus.TBC.2024    4 27 8506.78 6240.98   6092 7802.57 3877.00 968 29110 28142
##                skew kurtosis      se
## Kasus.TBC.2021 1.51     2.20  510.12
## Kasus.TBC.2022 1.57     2.73  880.47
## Kasus.TBC.2023 1.48     2.33 1141.87
## Kasus.TBC.2024 1.37     2.01 1201.08

Dari deskriptif bisa dilihat karakteristik pada data jumlah total kasus TBC di Jawa Barat, namun hal yang perlu diperhatikan merupakan peningkatan rata-rata serta median dan juga nilai minimum dan maximum dari tahun ke tahun menunjukkan sebuah tren naik, untuk hal tersebut akan didalami pada trend line. Didapatkan juga visualisasi histogram jumlah kasus TBC di Jawa Barat:

library(ggplot2)
library(scales)

plot_histogram_tbc <- function(data, var, tahun, warna, bins = 10) {
  x <- data[[var]]
  
  # Hitung breaks manual
  breaks <- pretty(range(x, na.rm = TRUE), n = bins)
  
  ggplot(data, aes(x = x)) +
    geom_histogram(
      breaks = breaks,
      fill = warna,
      color = "white"
    ) +
    scale_x_continuous(
      breaks = breaks,     # ini bikin label batas bawah & atas
      labels = comma       # biar format ribuan
    ) +
    labs(
      title = paste("Distribusi Jumlah Kasus TBC Tahun", tahun),
      x = "Batas Interval Jumlah Kasus",
      y = "Frekuensi"
    ) +
    theme_minimal(base_size = 13) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
}

# Plot satu per satu
plot_histogram_tbc(jabar_merged, "Kasus.TBC.2021", 2021, "#5E4FA2")

plot_histogram_tbc(jabar_merged, "Kasus.TBC.2022", 2022, "#3288BD")

plot_histogram_tbc(jabar_merged, "Kasus.TBC.2023", 2023, "#66C2A5")

plot_histogram_tbc(jabar_merged, "Kasus.TBC.2024", 2024, "#F46D43")

Dilihat bahwa distribusi kebanyakan berbentuk right-skewed, mayoritas kabupaten memiliki jumlah kasus yang relatif rendah dari rata-rata.

Salah satu tantangan yang dihadapi saat melakukan analisis adalah heterogenitas populasi, dimana kasus yang banyak merupakan indikasi bahwa karena ada lebih banyak populasi sehingga akan lebih banyak kasus. Karena itu, digunakannya prevalence rate untuk menangani hal tersebut dalam bentuk persentase. Didapatkan plot map seperti berikut:

interval_labels <- function(breaks) {
  labs <- c()
  for (i in seq_along(breaks)[-length(breaks)]) {
    labs[i] <- paste0(
      format(breaks[i], big.mark = ","),
      " – ",
      format(breaks[i + 1], big.mark = ",")
    )
  }
  labs
}


max_val <- max(
  jabar_merged$PR_2021,
  jabar_merged$PR_2022,
  jabar_merged$PR_2023,
  jabar_merged$PR_2024,
  na.rm = TRUE
)

breaks_interval <- seq(0, ceiling(max_val / 100) * 100, by = 100)
label_interval <- interval_labels(breaks_interval)

plot_pr_interval <- function(data, var, year) {
  ggplot(data) +
    geom_sf(aes(fill = .data[[var]]),
            color = "white", linewidth = 0.3) +
    scale_fill_viridis_c(
      option = "plasma",
      direction = -1,
      name = "Prevalence Rate (PR) 100000 Orang",
      breaks = breaks_interval[-1],  
      labels = label_interval
    ) +
    labs(
      title = paste("Prevalence Rate (PR) TBC di Jawa Barat", year),
      caption = "Sumber: Open Data Jabar"
    ) +
    theme_minimal(base_size = 14) +
    theme(
      panel.background = element_rect(fill = "white", color = NA),
      panel.grid.major = element_line(color = "gray90"),
      panel.grid.minor = element_line(color = "gray95"),
      plot.title = element_text(face = "bold", hjust = 0.5),
      legend.title = element_text(face = "bold"),
      legend.text = element_text(size = 9)
    )
}


plot_pr_interval(jabar_merged, "PR_2021", 2021)

plot_pr_interval(jabar_merged, "PR_2022", 2022)

plot_pr_interval(jabar_merged, "PR_2023", 2023)

plot_pr_interval(jabar_merged, "PR_2024", 2024)

Bisa dilihat sekarang dengan jelas, sebuah perbedaan signifikan. Dilihat pada data tahun 2021 bahwa pola dimana bagian timur memiliki prevalensi yang lebih rendah daripada di barat. Juga bisa dilihat bahwa ialah pada di pusat populasi, yaitu di kota, seperti di Kota Bandung, Kota Bogor, Kota Cirebon, Kota Banjar, Kota Tasikmalaya, dan Kota Sukabumi memiliki prevalensi yang tinggi dibanding yang berkabupaten. Jika dilihat dari tahun ke tahun juga bisa diamati berkembangnya prevalensi TBC di Jawa Barat dengan meningkatnya prevalensi di seluruh kabupaten Jawa Barat dan juga secara signifikan di Kota Bogor, namun pada 2024 bisa diamati bahwa persebaran tersebut memang tetap mengikuti distribusi barat yang tinggi serta timur yang rendah.

library(psych)
library(dplyr)

describe(
  jabar_merged %>%
    st_drop_geometry() %>%         
    select(PR_2021:PR_2024)
)
##         vars  n   mean     sd median trimmed    mad    min     max   range skew
## PR_2021    1 27 193.53 103.40 174.89  178.53  63.66  91.86  568.77  476.91 1.94
## PR_2022    2 27 359.36 188.73 298.99  340.72 130.84 145.67  814.30  668.63 1.09
## PR_2023    3 27 476.04 257.45 365.85  446.66 187.16 185.00 1206.80 1021.80 1.20
## PR_2024    4 27 505.03 251.10 438.86  474.62 156.63 222.99 1195.30  972.31 1.23
##         kurtosis    se
## PR_2021     4.09 19.90
## PR_2022     0.00 36.32
## PR_2023     0.65 49.55
## PR_2024     0.66 48.32

Sama seperti pada deskriptif jumlah kasus bisa dilihat dengan jelas kenaikan rata-rata serta median dan juga nilai nimum dan maksimum dari tahun ke tahun menunjukkan adanya peningkatan prevalensi data dari 2021-2024.

4.3 Trendline Data

Untuk melihat pola data terutama berhubungan dengan waktu diperlukan visualisasi data dari jangka waktu tersebut sebagai bentuk asesmen. Karena itu disajikan data jumlah kasus total di Jawa Barat serta rata-rata prevalensi rate di Bandung sebagai berikut:

library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.3
## 
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack
## The following object is masked from 'package:terra':
## 
##     extract
trend_kasus <- jabar_merged %>%
  st_drop_geometry() %>%   
  summarise(
    Total_2021 = sum(Kasus.TBC.2021, na.rm = TRUE),
    Total_2022 = sum(Kasus.TBC.2022, na.rm = TRUE),
    Total_2023 = sum(Kasus.TBC.2023, na.rm = TRUE),
    Total_2024 = sum(Kasus.TBC.2024, na.rm = TRUE)
  ) %>%
  pivot_longer(cols = everything(), names_to = "Tahun", values_to = "Total_Kasus") %>%
  mutate(Tahun = as.numeric(gsub("Total_", "", Tahun)))

# Plot
ggplot(trend_kasus, aes(x = Tahun, y = Total_Kasus)) +
  geom_col(fill = "#3288BD") +
  geom_line(aes(group = 1), color = "darkred", size = 1.2) +
  geom_point(color = "darkred", size = 3) +
  labs(
    title = "Tren Total Kasus TBC Jawa Barat (2021–2024)",
    x = "Tahun",
    y = "Jumlah Kasus"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Trendline Prevalensi
trend_prev <- jabar_merged %>%
  st_drop_geometry() %>%
  summarise(
    MeanPR_2021 = mean(PR_2021, na.rm = TRUE),
    MeanPR_2022 = mean(PR_2022, na.rm = TRUE),
    MeanPR_2023 = mean(PR_2023, na.rm = TRUE),
    MeanPR_2024 = mean(PR_2024, na.rm = TRUE)
  ) %>%
  pivot_longer(cols = everything(), names_to = "Tahun", values_to = "Mean_Prevalensi") %>%
  mutate(Tahun = as.numeric(gsub("MeanPR_", "", Tahun)))

ggplot(trend_prev, aes(x = Tahun, y = Mean_Prevalensi)) +
  geom_col(fill = "#F46D43") +
  geom_line(aes(group = 1), color = "darkblue", size = 1.2) +
  geom_point(color = "darkblue", size = 3) +
  labs(
    title = "Tren Rata-rata Prevalensi TBC Jawa Barat (2021–2024)",
    x = "Tahun",
    y = "Prevalensi (per 100.000 penduduk)"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

Dilihat dengan jelas pada trendline yang diberikan total kasus TBC di Jawa Barat serta rata-rata prevalensi TBC di Jawa Barat bahwa Jawa Barat sedang menghadapi sebuah krisis kesehatan dalam bentuk peningkatan kasus TBC serta banyaknya beban populasi yang menderita TBC. Sebuah krisis TBC sedang terjadi di Jawa Barat, dan dibutuhkan respon yang cepat serta efektif.

4.4 Autokorelasi Spasial

Setelah dilakukan analisis deskriptif menggunakan prevalensi rate, bisa dilakukan analisis lebih dalam, terutama pada efek spasial. Dilakukan analisis untuk melihat besarnya efek dependensi spasial pada penyebaran TBC, dan bagaimana hal tersebut bisa diambil sebagai konsiderasi dalam merespon ataupun menghambat penyebarannya TBC di masa yang akan datang. Akan digunakan uji Global Moran’s I dan Geary’s C untuk Uji Global. Uji Local Moran dan Uji Ger-Ord akan digunakan untuk Uji Lokal.

Pertama akan dibuatkan weight menggunakan row-based proportion dengan jarak queen, serta pembuatan Peta Ketentanggaan per Kabupaten/Kota.

# Analisis Autokorelasi
nb <- poly2nb(as_Spatial(jabar_merged), queen = TRUE)
lwW <- nb2listw(nb, style = "W", zero.policy = TRUE)
lwB <- nb2listw(nb, style = "B", zero.policy = TRUE)

jabar_sp <- as_Spatial(jabar_merged)
coords <- coordinates(jabar_sp)
plot(jabar_sp, col = "gray90", border = "gray60", main = "Peta Ketetanggaan Kabupaten/Kota Jawa Barat")
points(coords[,1], coords[,2], pch = 19, cex = 0.7, col = "blue")
text(coords[,1], coords[,2], labels = jabar_sp$admin2Name_en, cex = 0.6, pos = 3)
plot(nb, coords, add = TRUE, col = "red", lwd = 1.5)

Kemudian bisa dilakukan pengujian autokorelasi spasial, pertama adalah pembuatan function untuk mempermudah analisis

analyze_spatial <- function(data_sf, var_col, label) {
  cat("\n\n==============================\n", label, "\n==============================\n")
  
  x_raw <- data_sf[[var_col]]
  
  # === Uji autokorelasi spasial global (asimptotik) ===
  moran_res <- moran.test(x_raw, listw = lwW, zero.policy = TRUE)
  geary_res <- geary.test(x_raw, listw = lwW, zero.policy = TRUE)
  
  print(moran_res)
  print(geary_res)
  
  # === LISA (Local Moran’s I) ===
  x <- scale(x_raw)[, 1]
  lagx <- lag.listw(lwW, x, zero.policy = TRUE)
  lisa <- localmoran(x, lwW, alternative = "two.sided", zero.policy = TRUE)
  lisa_df <- as.data.frame(lisa)
  names(lisa_df) <- c("Ii","Ei","Vi","Zi","Pi.two.sided")
  
  alpha <- 0.05
  quad <- dplyr::case_when(
    x >= 0 & lagx >= 0 ~ "High-High",
    x < 0 & lagx < 0 ~ "Low-Low",
    x >= 0 & lagx < 0 ~ "High-Low (Outlier)",
    x < 0 & lagx >= 0 ~ "Low-High (Outlier)"
  )
  
  LISA_sf <- bind_cols(data_sf, lisa_df) |>
    mutate(quad = ifelse(Pi.two.sided <= alpha, quad, "Not significant"))
  
  # === Local Geary’s C ===
  localC_vals <- localC(x_raw, lwW, zero.policy = TRUE)
  Geary_sf <- mutate(data_sf, LocalC = as.numeric(localC_vals))
  
  # === Getis–Ord Gi* ===
  Wb <- listw2mat(lwB)
  Wb_star <- Wb; diag(Wb_star) <- 1
  G_star_raw <- as.numeric(Wb_star %*% x_raw) / sum(x_raw)
  Gz <- localG(x_raw, listw = lwW, zero.policy = TRUE)
  
  G_sf <- mutate(data_sf,
                 z_Gistar = as.numeric(Gz),
                 hotcold = case_when(
                   z_Gistar >= 1.96 ~ "Hot spot (p<0.05)",
                   z_Gistar <= -1.96 ~ "Cold spot (p<0.05)",
                   TRUE ~ "Not significant"
                 ))
  
  # === Visualisasi ===
  p1 <- ggplot(LISA_sf) +
    geom_sf(aes(fill = quad), color = "white", size = 0.2) +
    scale_fill_manual(values = c(
      "High-High" = "#d73027",
      "Low-Low" = "#4575b4",
      "High-Low (Outlier)" = "#fdae61",
      "Low-High (Outlier)" = "#74add1",
      "Not significant" = "grey85"
    )) +
    labs(title = paste("Local Moran's I (LISA) —", label),
         fill = "Kategori") +
    theme_minimal()
  
  p2 <- ggplot(Geary_sf) +
    geom_sf(aes(fill = LocalC), color = "white", size = 0.2) +
    scale_fill_viridis_c(option = "magma", direction = -1,
                         name = "Local Geary’s C") +
    labs(title = paste("Local Geary’s C —", label)) +
    theme_minimal()
  
  p3 <- ggplot(G_sf) +
    geom_sf(aes(fill = hotcold), color = "white", size = 0.2) +
    scale_fill_manual(values = c(
      "Hot spot (p<0.05)" = "#b2182b",
      "Cold spot (p<0.05)" = "#2166ac",
      "Not significant" = "grey85"
    )) +
    labs(title = paste("Getis–Ord Gi* —", label), fill = NULL) +
    theme_minimal()
  
  print(p1); print(p2); print(p3)
}

Dilakukan analisis spasial autokorelasi untuk prevalensi tahun 2021 dengan hasil sebagai berikut:

analyze_spatial(jabar_merged, "PR_2021", "Prevalensi Rate 2021")
## 
## 
## ==============================
##  Prevalensi Rate 2021 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW    
## 
## Moran I statistic standard deviate = 0.45122, p-value = 0.3259
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.01746724       -0.03846154        0.01536363 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW   
## 
## Geary C statistic standard deviate = 1.281, p-value = 0.1001
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.71732764        1.00000000        0.04869638

Didapatkan bahwa prevalensi rate 2021 tidak memiliki autokorelasi global yang signifikan, dengan Moran’s I sebesar 0.02 dan Geary sebesar 0.7. Menunjukkan bahwa pola TBC pada saat tersebut bisa dikatakan acak. Namun secara lokal terdapat cold spot dimana area dengan prevalensi TBC yang rendah disekelilingi yang rendah yakni di Kabupaten Majalengka.

Dilakukan analisis spasial autokorelasi untuk prevalensi tahun 2022 dengan hasil sebagai berikut:

analyze_spatial(jabar_merged, "PR_2022", "Prevalensi Rate 2022")
## 
## 
## ==============================
##  Prevalensi Rate 2022 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW    
## 
## Moran I statistic standard deviate = 0.84992, p-value = 0.1977
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.07889033       -0.03846154        0.01906451 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW   
## 
## Geary C statistic standard deviate = 1.803, p-value = 0.03569
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.69138634        1.00000000        0.02929659

Seperti sebelumnya, tidak terdapat autokorelasi signifikan pada prevalensi TBC di tahun 2022. Namun uji Geary mendapatkan hasil yang signifikan dengan nilai sebesar 0.69, walaupun Moran’s I tidak dengan nilai sebesar 0.07. Secara lokal didapatkan juga cold spot di Kabupaten Majalengka.

Dilakukan analisis spasial autokorelasi untuk prevalensi tahun 2023 dengan hasil sebagai berikut:

analyze_spatial(jabar_merged, "PR_2023", "Prevalensi Rate 2023")
## 
## 
## ==============================
##  Prevalensi Rate 2023 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW    
## 
## Moran I statistic standard deviate = 0.58868, p-value = 0.278
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.04155478       -0.03846154        0.01847590 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW   
## 
## Geary C statistic standard deviate = 1.5209, p-value = 0.06415
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.72632061        1.00000000        0.03238205

Seperti sebelumnya tidak terdapat autokorelasi global signifikan berdasarkan pengujian Global Moran’s serta Geary C, dengan nilai sebesar 0.04 untuk Moran dan 0.7 untuk Geary. Namun Geary juga mendapatkan hasil yang hampir signifikan. Dan untuk lokal bisa didapatkan Kabupaten Majalengka sebagai coldspot seperti sebelumnya namun juga sekarang terdapat Kabupaten Sumedang sebagai coldspot juga.

Dilakukan analisis spasial autokorelasi untuk prevalensi tahun 2024 dengan hasil sebagai berikut:

analyze_spatial(jabar_merged, "PR_2024", "Prevalensi Rate 2024")
## 
## 
## ==============================
##  Prevalensi Rate 2024 
## ==============================
## 
##  Moran I test under randomisation
## 
## data:  x_raw  
## weights: lwW    
## 
## Moran I statistic standard deviate = 0.60773, p-value = 0.2717
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.04413265       -0.03846154        0.01847043 
## 
## 
##  Geary C test under randomisation
## 
## data:  x_raw 
## weights: lwW   
## 
## Geary C statistic standard deviate = 1.6571, p-value = 0.04875
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.70168080        1.00000000        0.03241075

Seperti sebelumnya, tidak terdapat autokorelasi signifikan pada prevalensi TBC di tahun 2022. Namun uji Geary mendapatkan hasil yang signifikan dengan nilai sebesar 0.70, walaupun Moran’s I tidak dengan nilai sebesar 0.04. Secara lokal didapatkan juga cold spot di Kabupaten Majalengka.

4.5 Interpretasi Epidemiologis

Hasil deskriptif menunjukkan peningkatan jumlah kasus dan prevalensi TBC di Jawa Barat dari 2021 sampai 2024, dengan konsentrasi beban di daerah timur di Jawa Barat yang berdekatan dengan area industri seperti di Bogor, Bekasi, dan Karawang. Serta juga didapatkan prevalensi rendah terutama di daerah barat di Jawa Barat, dimana Kabupaten seperti Majalengka dan kepada batas tertentu, Kabupaten Sumedang memiliki prevalensi yang rendah. Serta diidentifikasikannya area-area perkotaan seperti Kota Bogor, Kota Bandung, dan pusat urba lainnya yang memiliki prevalensi yang sangatlah tinggi dibanding kabupateb lainnya. Analisis spasial global tidak menunjukkan autokorelasi signifikan secara global menunjukan bahwa TBC tidak memiliki pengaruh spasial yang terlalu signifkman, tetapi analisis lokal mengidentifikasi hotspot dan coldspot yang relevan untuk intervensi program. Dibutuhkan riset lebih dalam mengenai deteksinya TBC, pelaporan, serta juga analisis untuk mendapatkan informasi terbagus dalam melaksanakan kebijakan

5 Simulasi Desain Eksperimen

5.1 Studi Kasus-Kontrol: Faktor Risiko Tuberkulosis (TBC)

Studi kasus-kontrol adalah studi observasional retrospektif yang memulai dengan memilih kelompok kasus (penderita) dan kontrol (bukan penderita), lalu menelusuri paparan di masa lalu untuk menguji hipotesis hubungan paparan–penyakit. Pendekatan ini cocok untuk penyakit dengan insidens relatif rendah dan hipotesis awal tentang faktor risiko. Hasil analisis utama adalah Odds Ratio (OR) sebagai ukuran asosiasi antara paparan dan penyakit diikuti oleh pemodelan regresi logistik.

Berdasarkan riset Kevelaitiene, Davidaviciene, dan Danila (2022) mendapatkan faktor risiko pada sebuah host untuk dalam faktor terkenanya penyakit TBC yaitu penggunaan rokok dan adanya riwayat menderita diabetes. Serta itu faktor biologis seperti kelamin juga memberikan faktor terhadap terjangkitnya TBC atau tidak, dimana lelaki memiliki rata-rata lebih tinggi daripada wanita, namun alasan yang memberikan hal tersebut masih dipertanyakan. Karena itu, tujuan riset adalah melakukan konfirmasi terhadap hipotesis yang diberikan oleh (Kevelaitiene et al. 2022) serta melihat lebih dalam apa yang memengaruhi seseorang terjangkit penyakit TBC.

5.2 Populasi Target dan Sampling

Populasi yang menjadi unit eksperimen merupakan penduduk dewasa (>18 tahun) pada wilayah Jawa Barat. Sampling akan dilakukan dengan mengambil daftar pasien yang sedang menderita TBC yang ada di rumah sakit/klinik di Jawa Barat, diprioritaskan pasien-pasien yang relatif baru menderita penyakit TBC untuk mengurangu adanya Berkson Bias. Untuk memastikan benarinya sampling serta mengurangi bias juga diberlakukan skrining ketat bagi yang menderita untuk memastikan bahwa individu benari-benar terdiagnosa menderita TBC oleh dokter. Untuk kontrol, individu yang tidak menderita TBC ataupun punya riwayat TBC, akan dipilih berdasarkan dari individu dari komunitas sama contohnya adalah tetangga dari individu yang menderita ataupun teman/rekannya. Untuk mengurangi bias maka akan dilakukan frequency matching untuk ememastikan persamaan karakteristik pada kontrol serta yang menderita TBC dalam aspek etnis dan umur untuk mengurangi bias. Serta diberlakukan skrining untuk memastikan tidak adanya riwayat menderita TBC.

5.3 Variabel Utama

Status TBC (Dependen) Pasien dewasa yang menderita dari penyakit TBC melalui diagnosis pada fasilitas kesehatan (sudah terkonfirmasi). Kontrol berupa individu dewasa yang tidak menderita serta mempunyai riwayat TBC.

Variabel Independen (Faktor Risiko Host)

  1. Merokok (Ya/Tidak) – riwayat merokok aktif sebelum diagnosis TBC.
  2. Diabetes (Ya/Tidak) – adanya riwayat diabetes melitus.
  3. Kelamin (Pria/Wanita)

5.4 Desain Studi dan Langkah Penelitian

  1. Penentuan Kasus dan Kontrol: Kasus: Semua penderita TBC terkonfirmasi dalam periode tertentu di Jawa Barat.

Kontrol: Individu dari populasi sama tanpa tanpa TBC serta riwayat TBC dipilih berdasarkan frequency matching.

  1. Pengukuran Paparan: Melalui wawancara dilihat status perokok dan riwayat diabetes. Ditanyakan pada kontrol serta yang menderita TBC.

  2. Langkah Analisis: Setelah data siap, analisis yang pertama dilakukan berupa asosiasi awal (uji Odds Ratio) antara TBC dan setiap variabel, lalu analisis dengan regresi logistik untuk mengontrol variabel confounder.

5.5 Potential Bias

  1. Selection Bias Adanya bias pada pemilihan unit observasi terutama pada yang menderita TBC karena kemungkinan yang masuk pada rumah sakit tidak mewakili semua penderita TBC terutama yang sedang rawat inap sehingga melakukan distorsi pada OR. Serta juga pada pemolihan kontrol yang bisa tidak mewakili sumber populasi sehingga membuat ambigu pada analisis data.

  2. Prevalence-Incidence Bias Memasuki pasien yang sudah lama (prevalent) dimana bisa mudah hilang dari sampling karena sudah meninggal, menyebabkan underestimasi OR.

  3. Recall Bias Melaporkan paparan faktor risiko yang lebih intens daripada kontrol sehingga menimbulkan bias pada asosiasi dimana odds akan terlihat lebih besar karena kurangnya intensitas pada kontrol.

  4. Measurement Error Adanya permasalahn mendata terutama mengenai faktor risiko dimana ada kasus underreporting karena tekanan sosial sehingga banyak yang tidak mengaku dirinya sebagai perokok ataupun riwayat diabetes sehingga terjadinya underreporting.

  5. Confounding Variabel Adanya variabel ketiga yang memengaruhi peluangnya seseorang terkena penyakit TBC serta pada variabel independen seperti variabel sosio-ekonomik, HIV, dan juga pekerjaan yang tidak dimasukkan model.

5.6 Analisis Simulasi

Sudah didapatkan data sebagai berikut:

library(dplyr)
library(broom)
## Warning: package 'broom' was built under R version 4.3.3
library(pscl)
## Warning: package 'pscl' was built under R version 4.3.3
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
library(ggplot2)

set.seed(42)


n_case <- 250
n_ctrl <- 250
n <- n_case + n_ctrl

Gender <- sample(c("Male", "Female"), n, replace = TRUE, prob = c(0.55, 0.45))
Smoking <- rbinom(n, 1, ifelse(Gender == "Male", 0.55, 0.35)) 
Diabetes <- rbinom(n, 1, 0.15)  

b0 <- -2.5
b1 <- 1.2     
b2 <- 0.9     
b3 <- 0.6     

logit_p <- b0 + b1*Smoking + b2*Diabetes + b3*(Gender=="Male")
p <- plogis(logit_p)


TBC <- rbinom(n, 1, p)

TBC_case_idx <- which(TBC == 1)
TBC_ctrl_idx <- which(TBC == 0)
set.seed(42)
idx <- c(sample(TBC_case_idx, n_case, replace = TRUE),
          sample(TBC_ctrl_idx, n_ctrl, replace = TRUE))

data_tbc <- tibble(
  TBC = TBC[idx],
  Smoking = Smoking[idx],
  Diabetes = Diabetes[idx],
  Gender = factor(Gender[idx], levels = c("Female","Male"))
)

table(data_tbc$TBC)
## 
##   0   1 
## 250 250
table(data_tbc$Gender, data_tbc$TBC)
##         
##            0   1
##   Female 122  80
##   Male   128 170

Menghitung OR

calc_or <- function(var){
  tab <- table(data_tbc[[var]], data_tbc$TBC)
  or <- (tab[2,2] / tab[1,2]) / (tab[2,1] / tab[1,1])
  or
}

cat("Crude OR Smoking:", round(calc_or("Smoking"), 2), "\n")
## Crude OR Smoking: 4.11
cat("Crude OR Diabetes:", round(calc_or("Diabetes"), 2), "\n")
## Crude OR Diabetes: 2.32
cat("Crude OR Gender (Male vs Female):", round(calc_or("Gender"), 2), "\n")
## Crude OR Gender (Male vs Female): 2.03

Model Logistik

model <- glm(TBC ~ Smoking + Diabetes + Gender, data = data_tbc, family = binomial)
summary(model)
## 
## Call:
## glm(formula = TBC ~ Smoking + Diabetes + Gender, family = binomial, 
##     data = data_tbc)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.1097     0.1823  -6.089 1.14e-09 ***
## Smoking       1.3243     0.1968   6.730 1.69e-11 ***
## Diabetes      0.7912     0.2752   2.876  0.00403 ** 
## GenderMale    0.4317     0.1993   2.166  0.03033 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 693.15  on 499  degrees of freedom
## Residual deviance: 621.06  on 496  degrees of freedom
## AIC: 629.06
## 
## Number of Fisher Scoring iterations: 4
ll_full <- as.numeric(logLik(model))
ll_null <- as.numeric(logLik(glm(TBC ~ 1, data = data_tbc, family = binomial)))
mcfadden <- 1 - (ll_full / ll_null)
cat("McFadden’s pseudo-R² =", round(mcfadden, 3), "\n")
## McFadden’s pseudo-R² = 0.104
result <- broom::tidy(model) %>%
  mutate(
    OR = exp(estimate),
    OR_low = exp(estimate - 1.96 * std.error),
    OR_high = exp(estimate + 1.96 * std.error),
    p.value = signif(p.value, 3)
  ) %>%
  select(term, OR, OR_low, OR_high, p.value)

print(result)
## # A tibble: 4 × 5
##   term           OR OR_low OR_high  p.value
##   <chr>       <dbl>  <dbl>   <dbl>    <dbl>
## 1 (Intercept) 0.330  0.231   0.471 1.14e- 9
## 2 Smoking     3.76   2.56    5.53  1.69e-11
## 3 Diabetes    2.21   1.29    3.78  4.03e- 3
## 4 GenderMale  1.54   1.04    2.28  3.03e- 2
plot_df <- result %>%
  filter(term != "(Intercept)") %>%
  mutate(term = case_when(
    term == "Smoking" ~ "Smoking (Yes vs No)",
    term == "Diabetes" ~ "Diabetes (Yes vs No)",
    term == "GenderMale" ~ "Male vs Female",
    TRUE ~ term
  ))

ggplot(plot_df, aes(x = OR, y = term)) +
  geom_point(size = 3, color = "#2E86AB") +
  geom_errorbarh(aes(xmin = OR_low, xmax = OR_high), height = 0.2, color = "#2E86AB") +
  geom_vline(xintercept = 1, linetype = "dashed", color = "gray60") +
  scale_x_log10() +
  labs(
    title = "Adjusted Odds Ratios (Studi Kasus–Kontrol TBC)",
    x = "Odds Ratio (log scale)",
    y = ""
  ) +
  theme_minimal(base_size = 13)

5.7 Kesimpulan Eksperimen

Orang yang merokok memiliki odds 4,1.1 kali lebih besar untuk mendapatkan TBC dibandingkan yang tidak. Serta juga orang yang mempunyai riwayat diabetes memiliki odds 2.32 kali lebih besar untuk terkena penyakit TBC dibanding yang tidak. Dan lelaki memiliki odds 2 kali lebih besar dari wanita untuk terkena TBC.

Berdasarkan dari hasil regresi logistik, individu yang merokok tetap memiliki odds 3.76 kali lebih besar untuk menderita TBC serta juga diabetes dengan odds 2.21 kali lebih besar daripada yang tidak ada riwayat diabetes untuk terkena TBC. Laki-laki memiliki odds 1.54 kali lebih tinggi dibanding perempuan, setelah memperhitungkan smoking dan diabetes.

6 Kesimpulan dan Saran

Berdasarkan riset yang dilakukan dibutuhkan sebuah riset besar-besaran terhadap deteksi penyakit TBC yang ada di Jawa Barat dikarenakan meningkatnya kasus serta prevalensi TBC pada populasi Jawa Barat. Karena itu dibutuhkan koordinasi besar antar pemerintah lokal serta institusi kesehatan untuk mencrgah penyebaran TBC dengan konsiderasi faktor risiko pada host berupa kelamin, perokok, dan riwayat diabetes untuk memudahkan pencegahan serta antisipasi penyebaran TBC.

7 Daftar Pustaka

Badan Pusat Statistik Provinsi Jawa Barat. (2023). Provinsi Jawa Barat dalam angka 2023. BPS Provinsi Jawa Barat.

Budiarto, E. (2012). Metodologi penelitian kedokteran: Sebuah pengantar. EGC.

Celentano, D. D., & Szklo, M. (2018). Gordis epidemiology (6th ed.). Elsevier.

Dinas Kesehatan Provinsi Jawa Barat. (2023). Profil kesehatan Provinsi Jawa Barat tahun 2023. Dinas Kesehatan Provinsi Jawa Barat.

Nicole Fogel, Tuberculosis: A disease without boundaries, Tuberculosis, Volume 95, Issue 5, 2015, Pages 527-531, ISSN 1472-9792, https://doi.org/10.1016/j.tube.2015.05.017.

Kementerian Kesehatan RI. (2023). Pedoman nasional pelayanan kedokteran: Tata laksana tuberkulosis. Kementerian Kesehatan RI.

Kementerian Kesehatan RI. (2024). Profil kesehatan Indonesia 2023. Kementerian Kesehatan RI.

Kėvelaitienė, K., Davidavičienė, V. E., & Danila, E. (2025). Tuberculosis treatment failure: what are the risk factors? A comprehensive literature review. Multidisciplinary Respiratory Medicine, 20. https://doi.org/10.5826/mrm.2025.1030

Rothman, K. J. (2012). Epidemiology: An introduction (2nd ed.). Oxford University Press.

World Health Organization. (2024). Global tuberculosis report 2024. World Health Organization.

8 Lampiran

LINK DASHBOARD https://chrisrafael.shinyapps.io/Dashboard_Epidemiologi/

TABEL DATA

setwd("C:\\Users\\HP\\Documents\\KULIAH\\SEMESTER 5\\Epidem")
data <- (read.csv("ceremony.csv", sep=",", dec=".")); data
##      Kabupaten.Kota Populasi.2021 Populasi.2022 Populasi.2023 Populasi.2024
## 1           BANDUNG       3652400       3687250       3721110       3753120
## 2     BANDUNG BARAT       1808420       1834230       1859640       1884190
## 3            BEKASI       3148740       3193840       3237420       3273870
## 4             BOGOR       5484150       5556310       5627020       5682300
## 5            CIAMIS       1234830       1243320       1251540       1259230
## 6           CIANJUR       2500640       2529810       2558140       2584990
## 7           CIREBON       2301330       2331360       2360440       2387960
## 8             GARUT       2613530       2648950       2683670       2716950
## 9         INDRAMAYU       1851730       1873400       1894330       1914040
## 10         KARAWANG       2465570       2496190       2526000       2554380
## 11     KOTA BANDUNG       2461410       2484150       2506600       2528160
## 12      KOTA BANJAR        202720        205140        207510        209790
## 13      KOTA BEKASI       2568020       2598070       2627210       2644060
## 14       KOTA BOGOR       1050920       1060940       1070720       1078350
## 15      KOTA CIMAHI        574450        582650        590780        598700
## 16     KOTA CIREBON        335810        338940        341980        344850
## 17       KOTA DEPOK       2081130       2113620       2145400       2163640
## 18    KOTA SUKABUMI        350150        355420        360640        365740
## 19 KOTA TASIKMALAYA        723100        732480        741760        750730
## 20         KUNINGAN       1175950       1189010       1201760       1213930
## 21       MAJALENGKA       1315010       1328010       1340620       1352540
## 22      PANGANDARAN        425590        428600        431460        434100
## 23       PURWAKARTA       1008930       1023180       1037070       1050340
## 24           SUBANG       1620700       1635560       1649820       1663160
## 25         SUKABUMI       2747450       2775310       2802400       2828020
## 26         SUMEDANG       1159260       1168840       1178240       1187130
## 27      TASIKMALAYA       1876890       1892220       1907050       1920920
##    Kasus.TBC.2021 Kasus.TBC.2022 Kasus.TBC.2023 Kasus.TBC.2024
## 1            5708          10650          12697          14012
## 2            1736           3043           4178           4569
## 3            4798           8379          13515          15458
## 4           11946          21516          27690          29110
## 5            1606           2946           3042           3477
## 6            4673           7059           9327          12197
## 7            3383           6981           8238           9587
## 8            4786           7920           8615           9516
## 9            1701           2729           5148           6092
## 10           4528           8020          12868          13474
## 11           8919          14541          18327          18627
## 12            269            742           1175           1519
## 13           6051          10747          14279          13640
## 14           1440           7740          10354          11418
## 15           1782           4321           4602           4513
## 16           1910           2760           4127           4122
## 17           4042           6549           7849           8422
## 18           1223           2383           3431           3477
## 19           1476           2838           4749           4693
## 20           1651           2441           3705           4045
## 21           1723           3162           4318           4517
## 22            410            734            954            968
## 23           2410           4589           5607           4687
## 24           2912           4769           5666           6321
## 25           4805           7491          10950          12411
## 26           1372           2501           3020           3930
## 27           1995           3110           3528           4881

FILE PETA https://drive.google.com/drive/folders/15nmgQnOc58jNLIaXt10To6aZyqdRigT2?usp=sharing