KustaJawaTengah2025

LAPORAN PROYEK UJIAN AKHIR SEMESTER
ANALISIS EPIDEMIOLOGI KASUS PENYAKIT
KUSTA DI PROVINSI JAWA TENGAH
EPIDEMIOLOGI 2025

Disusun Oleh:
Farah Nadira
NPM: 140610230015

Dosen Mata Kuliah:
Dr. I Gede Nyoman Mindra Jaya, M.Si.
NIP: 19800603 200312 1 002

PROGRAM STUDI S-1 STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
JATINANGOR
2025

ABSTRAK

Indonesia masih menempati peringkat ketiga kasus kusta tertinggi di dunia, dengan Provinsi Jawa Tengah sebagai salah satu penyumbang kasus terbesar yang menunjukkan tren stagnasi. Penelitian ini bertujuan untuk menganalisis determinan sosial-ekonomi dan lingkungan terhadap jumlah kasus baru kusta di 35 kabupaten/kota di Jawa Tengah tahun 2024. Melalui pendekatan segitiga epidemiologi, penelitian ini menguji pengaruh kepadatan penduduk, persentase penduduk miskin, rata-rata lama sekolah (RLS), rumah layak huni, akses sanitasi, akses air minum, dan jumlah fasilitas kesehatan. Analisis dilakukan menggunakan model Regresi Binomial Negatif karena ditemukan overdispersi pada data. Berdasarkan kriteria kebaikan model (AIC terkecil), dipilih model StepAIC sebagai model terbaik. Hasil menunjukkan bahwa variabel Penduduk Miskin (\(p<0.01\)), Akses Sanitasi Layak (\(p<0.001\)), Kepadatan Penduduk (\(p<0.05\)), dan Rata-Rata Lama Sekolah ($ p<0.001 $) berpengaruh signifikan terhadap jumlah kasus kusta di Jawa Tengah. Rendahnya tingkat pendidikan (RLS) dan kondisi ekonomi menjadi hambatan utama dalam literasi kesehatan dan deteksi dini, sementara kepadatan penduduk serta sanitasi yang buruk mempercepat transmisi bakteri. Penelitian ini merekomendasikan intervensi terpadu yang berfokus pada peningkatan literasi kesehatan masyarakat dan perbaikan infrastruktur lingkungan di wilayah dengan kepadatan tinggi.

Kata Kunci: Kusta, Jawa Tengah, Morbus hansen, Mycobacterium leprae, Regresi Binomial Negatif.

PENDAHULUAN

Latar Belakang

Kusta atau Morbus hansen merupakan penyakit infeksi yang disebabkan oleh bakteri Mycobacterium leprae, bakteri ini merupakan mikroorganisme tahan asam yang memiliki ketertarikan khusus terhadap jaringan saraf manusia. Bakteri ini menginvasi sel saraf dan memicu peradangan yang menyebabkan mati rasa (anestesi), kelemahan otot, hingga kelumpuhan. Karena hilangnya sensasi rasa sakit, penderita sering kali tidak menyadari luka atau trauma fisik yang dialaminya, yang pada akhirnya mengakibatkan infeksi sekunder, pengeroposan tulang, dan mutilasi fisik yang dikenal sebagai cacat permanen.

Kusta pada dasarnya bukan merupakan penyakit yang mematikan ataupun berbahaya karena Mycobacterium leprae sendiri memiliki sifat pertumbuhan yang sangat lambat dan masa inkubasi yang panjang yaitu sekitar 2 hingga 5 tahun. Selain itu, transmisi atau penyebaran penyakit kusta sendiri dapat segera diputus hanya dengan satu dosis pertama Multi-Drug Therapy (MDT), MDT ini sendiri disediakan oleh pemerintah secara gratis di rumah sakit atau puskesmas di Indonesia. Namun, kemudahan medis ini tidak sebanding dengan realitas di lapangan. Di tahun 2023, Indonesia sebagai negara dengan kepadatan penduduk diperkirakan mencapai 158 jiwa per km^2 masih menempati peringkat ketiga kusta terbanyak di dunia dengan 17.251 kasus kusta terdaftar.

Jawa Tengah sebagai provinsi terpadat ketiga di Indonesia pun diketahui menjadi provinsi dengan jumlah kasus Kusta tertinggi ketiga di Indonesia, berdasarkan data yang diambil dari Open Data Jawa Tengah yakni sebanyak 1.424 kasus tercatat di tahun 2024, meningkat 37.3% dari tahun 2022 yang hanya berjumlah 932 kasus. Di tahun 2024 terjadi penurunan jumlah kasus baru kusta dari tahun 2023 sebanyak 0.39% (5 orang). Meskipun terdapat penurunan, angka ini terhitung sangat kecil dan menunjukkan bahwa tren kasus baru kusta di Jawa Tengah masih cenderung sejajar dan belum menunjukkan tanda-tanda eliminasi yang signifikan. Kesenjangan antara kemudahan yang ada dan kondisi faktual ini memunculkan tanda tanya besar dan menunjukkan bahwa tingginya kasus kusta di Jawa Tengah tidak semata-mata dipengaruhi oleh aspek medis, melainkan juga oleh aspek lain seperti sosial ekonomi, lingkungan, dan pendidikan.

Secara spesifik, aspek lingkungan dan ekonomi berkaitan dengan bagaimana kepadatan penduduk dan persentase penduduk miskin menciptakan ruang interaksi fisik yang intim dan berkelanjutan, yang merupakan syarat utama transmisi bakteri Mycobacterium leprae. Kondisi ini semakin berisiko apabila infrastruktur dasar seperti akses sanitasi layak, air minum layak, dan kondisi rumah layak huni tidak terpenuhi karena lingkungan yang lembap dan kurang higienis dapat memperpanjang daya tahan bakteri di luar tubuh manusia, sehingga memperbesar peluang terjadinya penularan dalam rumah tangga.

Di sisi lain, efektivitas pemutusan rantai penularan sangat bergantung pada kapasitas layanan melalui jumlah fasilitas kesehatan serta tingkat literasi kesehatan masyarakat yang direpresentasikan oleh Rata-rata Lama Sekolah (RLS). Pendidikan menjadi faktor krusial karena kusta tidak menunjukkan gejala sakit di awal infeksi, sehingga tanpa pengetahuan yang memadai, penderita cenderung mengabaikan gejala tersebut hingga mencapai tahap kecacatan permanen. Kondisi tersebut menegaskan urgensi penelitian kasus kusta di Jawa Tengah, dengan menggunakan model regresi binomial negatif bertujuan untuk mengidentifikasi faktor-faktor yang signifikan terhadap kasus baru kusta serta wilayah kabupaten/kota mana yang memerlukan perhatian khusus.

Rumusan Masalah

Berdasarkan uraian latar belakang tersebut, didapat rumusan masalah dalam penelitian ini antara lain:

1. Bagaimana gambaran kasus Kusta di Provinsi Jawa Tengah?
2. Bagaimana proses terjadinya penyakit Kusta ditinjau berdasarkan model agent-host-environment di Provinsi Jawa Tengah?
3. Wilayah mana yang perlu mendapatkan perhatian khusus dan menjadi prioritas utama dalam upaya pencegahan dan pengendalian Kusta di Provinsi Jawa Tengah?
4. Apa faktor yang sekiranya mempengaruhi kenaikan kasus baru Kusta di Provinsi Jawa Tengah?

Tujuan Penelitian

Dari rumusan masalah di atas, tujuan penelitian yang akhirnya dijadikan dasar dalam penelitian ini antara lain:

1. Mendeskripsikan gambaran kasus Kusta di Provinsi Jawa Tengah.
2. Menjelaskan proses terjadinya Kusta berdasarkan pendekatan agent-host-environment.
3. Mengetahui wilayah yang perlu mendapatkan perhatian khusus dan menjadi prioritas utama dalam upaya pencegahan dan pengendalian Kusta di Provinsi Jawa Tengah.
4. Mengetahui faktor yang berpengaruh signifikan terhadap kenaikan kasus baru Kusta di Provinsi Jawa Tengah.

TINJAUAN PUSTAKA

Konsep Dasar Epidemiologi

Epidemiologi adalah cabang ilmu kesehatan masyarakat yang mempelajari kejadian penyakit dalam suatu populasi, termasuk seberapa sering penyakit muncul (frekuensi), bagaimana pola penyebarannya berdasarkan waktu, tempat, dan karakteristik orang (distribusi), serta faktor-faktor yang memengaruhinya (determinan). Secara etimologi, istilah ini berasal dari bahasa Yunani epi (pada), demos (penduduk), dan logos (ilmu). Dalam praktiknya, epidemiologi digunakan untuk mengidentifikasi penyebab penyakit, memahami mekanisme penularan, serta merancang upaya pencegahan dan pengendalian agar masalah kesehatan dapat diminimalkan secara efektif di tingkat populasi.

Model Agent-Host-Environment

Model segitiga epidemiologi menjelaskan bahwa munculnya penyakit terjadi karena interaksi antara tiga komponen utama: agent, host, dan environment. Agent adalah penyebab penyakit, bisa berupa virus, bakteri, parasit, atau zat kimia yang menimbulkan gangguan kesehatan. Kemampuan agent untuk menular, bertahan hidup, dan memicu reaksi pada tubuh manusia menentukan seberapa besar potensi penularannya. Host adalah individu yang bisa terinfeksi, dengan tingkat kerentanan yang dipengaruhi oleh faktor seperti usia, status gizi, daya tahan tubuh, genetik, dan perilaku. Sementara itu, environment atau lingkungan mencakup faktor fisik (iklim, suhu, kelembapan), biologis (keberadaan vektor), dan sosial (kepadatan penduduk, kebersihan, perilaku masyarakat) yang dapat mendukung atau menghambat penularan penyakit.

Ketiga komponen ini saling berinteraksi dan menentukan apakah suatu penyakit bisa muncul atau tidak. Penyakit akan terjadi jika agent yang patogen bertemu dengan host yang rentan dalam lingkungan yang kondusif. Sebaliknya, jika salah satu faktor tidak mendukung, misalnya host memiliki kekebalan kuat atau lingkungan tidak sesuai untuk agent penyakit bisa dicegah.

Ukuran Epidemiologi

Di dalam studi epidemiologi terdapat dua jenis ukuran, yakni ukuran frekuensi yang digunakan untuk mengetahui frekuensi atau seberapa sering suatu penyakit terjadi dan ukuran asosiasi yang digunakan untuk mengetahui seberapa kuat hubungan faktor risiko dengan kejadian penyakit tersebut.

Baik ukuran frekuensi maupun ukuran asosiasi sangat penting untuk membantu pihak terkait dalam memahami pola penyakit dan menentukan prioritas dari penyelesaian faktor-faktor yang paling erat hubungannya.

Ukuran Frekuensi

Ukuran frekuensi menggambarkan seberapa sering suatu penyakit atau kejadian kesehatan terjadi pada populasi tertentu dalam periode waktu tertentu.

1. Prevalensi menunjukkan berapa banyak orang yang sakit (baik kasus baru maupun lama) pada satu waktu tertentu. Biasanya digunakan untuk mengukur beban penyakit di masyarakat.
   …(1)

2. Cumulative Incidence (Insidensi Kumulatif) menunjukkan proporsi individu yang menjadi kasus baru selama periode waktu tertentu. Ukuran ini digunakan untuk mengetahui risiko seseorang tertular atau terkena penyakit dalam periode tertentu.
   x 100%

3. Incidence Rate (Angka Insidensi) menunjukkan laju terjadinya kasus baru suatu penyakit dalam periode tertentu di antara populasi yang berisiko. Ukuran ini menggambarkan kecepatan munculnya kasus baru dalam suatu waktu.
   …(2)

4. Case Fatality Rate (CFR) menunjukkan persentase penderita yang meninggal akibat penyakit tersebut. CFR sering digunakan untuk mengukur tingkat keparahan penyakit, misalnya kalau CFR-nya tinggi berarti penyakit itu cenderung berbahaya
   

Ukuran Asosiasi

Ukuran asosiasi digunakan untuk menilai hubungan antara paparan (exposure) dan kejadian penyakit (outcome). Tujuannya adalah untuk mengetahui apakah suatu faktor benar-benar berperan sebagai penyebab atau faktor risiko penyakit.

1. Relative Risk (RR) digunakan pada penelitian kohort untuk membandingkan risiko penyakit antara kelompok terpapar dan tidak terpapar.
   

2. Odds Ratio (OR) digunakan pada penelitian kasus-kontrol untuk membandingkan peluang paparan antara kelompok kasus dan kontrol.
   

3. Attributable Risk (AR) atau risk difference menunjukkan selisih risiko antara kelompok terpapar dan tidak terpapar, menggambarkan seberapa besar proporsi kejadian penyakit yang bisa dicegah jika paparan dihilangkan.
   

4. Population Attributable Risk (PAR) menunjukkan seberapa besar proporsi penyakit di populasi yang disebabkan oleh paparan tertentu. Jadi, ukuran ini bantu nunjukin dampak faktor risiko terhadap keseluruhan populasi.
   

Desain Studi

Desain studi epidemiologi digunakan untuk meneliti hubungan antara paparan (exposure) dan kejadian penyakit (outcome) dalam suatu populasi. Secara umum, desain ini dibagi menjadi dua kelompok besar, yaitu studi observasional, di mana peneliti tidak melakukan intervensi terhadap subjek, dan studi eksperimental, di mana peneliti secara langsung memberikan perlakuan atau intervensi untuk melihat pengaruhnya terhadap outcome.

1. Cross-sectional study atau studi potong lintang
   Desain ini mengukur paparan dan penyakit pada waktu yang sama, sehingga memberikan gambaran tentang kondisi populasi pada saat tertentu. Studi ini biasanya digunakan untuk memperkirakan prevalensi penyakit atau faktor risiko, serta untuk mengeksplorasi hubungan awal antara variabel. Kelebihan studi ini adalah prosesnya cepat dan biayanya relatif murah, cocok untuk survei kesehatan masyarakat. Namun, karena tidak dapat menentukan urutan waktu antara paparan dan penyakit, studi ini tidak bisa memastikan hubungan sebab-akibat dan kurang tepat digunakan untuk penyakit langka.

2. Case-control study atau studi kasus-kontrol
   Desain yang membandingkan individu yang memiliki penyakit (kasus) dengan individu yang tidak memiliki penyakit (kontrol), kemudian menelusuri apakah mereka pernah terpapar faktor risiko di masa lalu. Studi ini bersifat retrospektif karena dimulai dari outcome menuju exposure. Ukuran asosiasi yang digunakan adalah Odds Ratio (OR). Keunggulan desain ini adalah efisien untuk penyakit yang jarang terjadi dan tidak membutuhkan waktu lama, namun memiliki kelemahan berupa potensi bias recall dan bias seleksi karena bergantung pada ingatan atau catatan masa lalu.

3. Cohort study atau studi kohort
   Desain ini menjelaskan di mana sekelompok individu yang terpapar dan tidak terpapar diikuti dalam jangka waktu tertentu untuk melihat apakah mereka mengalami penyakit atau tidak. Studi ini bisa dilakukan secara prospektif (ke depan) maupun retrospektif (menggunakan data historis). Ukuran asosiasi yang umum digunakan adalah Relative Risk (RR). Keunggulan studi kohort adalah kemampuannya menentukan hubungan sebab-akibat serta mengukur beberapa outcome sekaligus. Namun, kekurangannya adalah membutuhkan waktu lama, biaya besar, dan risiko kehilangan partisipan selama periode pengamatan.

4. Randomized Controlled Trial (RCT)
   RCT termasuk dalam kategori studi eksperimental. Pada desain ini, peserta dibagi secara acak menjadi dua kelompok, yaitu kelompok intervensi dan kelompok kontrol. Tujuannya adalah menilai efektivitas suatu intervensi atau terapi dengan meminimalkan bias melalui randomisasi dan kontrol. RCT dianggap sebagai desain dengan tingkat bukti paling kuat dalam menentukan hubungan sebab-akibat, meskipun dalam pelaksanaannya membutuhkan biaya besar, waktu lama, serta pengawasan etika yang ketat.

Pemodelan dengan Binomial Negatif

Model Regresi Binomial Negatif

Regresi Binomial Negatif adalah model statistik untuk data cacah (count data) yang digunakan ketika terjadi overdispersi, yaitu kondisi di mana varians data (\(\sigma^2\)) lebih besar dari rata-ratanya (\(\mu\)). Model ini memperbaiki kelemahan Regresi Poisson dengan menambahkan parameter dispersi (\(\alpha\)) agar estimasi lebih akurat. Hubungan antara variabel dependen dan independen dalam model ini dinyatakan sebagai berikut:

\[ ln(\mu)=\beta_0+\beta_1X_1+...+\beta_kX_k \]

Di mana variansnya didefinisikan sebagai:

\[ Var(Y)=\mu+\alpha\mu^2 \]

Pengujian Asumsi

Uji Multikolinearitas

Uji multikolinearitas digunakan untuk menguji apakah dalam model regresi terdapat korelasi yang tinggi atau sempurna antar variabel prediktor. Pengujian tersebut dapat diperiksa dengan melihat VIF (variance Inflation Factor) dan koefisien korelasi dari masing-masing variabel prediktor. Statistik uji yang digunakan adalah:

\[ VIF_k(u_i,v_i) = \frac{1}{(R^2_k(u_i,v_i))} \]

Keputusan yang diambil untuk uji multikolinearitas yaitu terjadi penolakan H0 ketika VIF > 5 atau terdapat multikolinearitas pada model Regresi Poisson.

Uji Overdispersi

Uji overdispersi dilakukan untuk memastikan apakah varians data melampaui nilai rata-ratanya, yang menjadi syarat utama penggunaan regresi Binomial Negatif. Secara statistik, kondisi ini dideteksi melalui nilai dispersi (\(\Phi = \frac{X^2}{df}\)); jika dispersion rate atau \(\Phi > 1\) atau p-value < 0,05, maka \(H_0\) ditolak dan data dinyatakan mengalami overdispersi. Jika data dinyatakan overdispersi, model regresi poisson tidak bisa digunakan dan digunakan model regresi binomial negatif sebagai alternatifnya.

Pemilihan Model Terbaik

Untuk mengevaluasi performa atau akurasi dari suatu model dan membandingkan antar model yang ada, diperlukan metrik untuk memilih model terbaik dari beberapa model yang ada. Metrik yang digunakan antara lain:

• AIC (Akaike Information Criterion) dan BIC (Bayesian Information Criterion), digunakan untuk mengukur efisiensi model. Model terbaik adalah yang memiliki nilai AIC dan BIC paling kecil.

  \[ AIC = -2ln(L)+2k \]

  \[ BIC = -2ln(L)+k.ln(n)\]

• McFadden’s Adjusted R-Square, digunakan untuk mengukur seberapa besar variabel independen mampu menjelaskan variasi variabel dependen, berbeda dengan \(R^2\) pada regresi linear, McFadden’s Adjusted R-Square memberikan penalti terhadap penambahan variabel yang tidak signifikan agar model tetap efisien. Ini dihitung dengan rumus:

  \[R_{adj}^2= 1-\frac{ln(L_{full})-k}{ln(L_{null})}\]

  Nilai \(R_{adj}^2\) yang terbesar yang paling mendekati 1 menunjukkan bahwa model memiliki kemampuan prediksi yang semakin baik dalam menjelaskan fenomena yang diteliti.

METODOLOGI PENELITIAN

Sumber Data

Data yang digunakan merupakan data sekunder yang berkaitan dengan Penyakit Kusta di Provinsi Jawa Tengah dengan 35 kabupaten/kota tahun 2024 yang bersumber dari sumber resmi seperti Badan Pusat Statistik (BPS) dan Open Data Provinsi Jawa Tengah.

Variabel Penelitian

Penelitian ini menggunakan beberapa variabel yang berkaitan dengan kejadian Penyakit Kusta di Provinsi Jawa Tengah. Variabel-variabel tersebut dibagi menjadi dua kelompok, yaitu variabel dependen (terikat) dan variabel independen (bebas). Variabel dependen dalam penelitian ini adalah jumlah kasus baru Kusta di setiap kabupaten/kota di Jawa Tengah yang menggambarkan tingkat penambahan kejadian penyakit di setiap tahunnya dan menjadi fokus utama analisis.

Sementara itu, variabel independen meliputi beberapa faktor lingkungan, sosial ekonomi, dan pendidikan yang diduga berhubungan dengan penyebaran Kusta di Jawa Tengah. Variabel-variabel tersebut antara lain:

1. Kepadatan penduduk, dipilih karena secara epidemiologis kusta membutuhkan kontak erat dalam durasi yang lama untuk dapat berpindah dari inang ke inang lainnya. Di wilayah dengan kepadatan tinggi, jarak fisik antarmanusia menjadi sangat terbatas, sehingga frekuensi terpapar percikan cairan pernapasan (droplet) penderita meningkat secara signifikan. Variabel ini digunakan untuk melihat sejauh mana kepadatan ruang hidup di Jawa Tengah mempermudah sirkulasi bakteri M. leprae di tengah masyarakat.

2. Persentase penduduk miskin, diambil karena kemiskinan sering kali berbanding lurus dengan keterbatasan luas hunian dan buruknya kualitas ruang mobilitas. Penduduk miskin cenderung tinggal dalam lingkungan yang sesak dan berbagi ruang tidur yang sama dengan banyak anggota keluarga, kondisi ini memenuhi syarat utama penularan kusta yaitu kontak erat yang lama. Selain itu, variabel ini penting untuk membedah bagaimana keterbatasan ekonomi melemahkan daya tahan tubuh (imun) melalui asupan gizi yang tidak memadai.

3. Persentase rumah layak huni, didasarkan pada karakteristik biologis bakteri kusta yang bersifat thermolabile atau sensitif terhadap suhu panas dan sinar matahari (UV). Rumah yang tidak layak huni biasanya memiliki ventilasi yang buruk dan pencahayaan alami yang minim, sehingga menciptakan suasana interior yang gelap dan lembap. Kondisi mikroklimat seperti ini diduga menjadi faktor kunci yang memungkinkan bakteri bertahan hidup lebih lama di luar tubuh manusia dan menginfeksi penghuni rumah lainnya.

4. Akses sanitasi dan air minum layak, dipilih karena kebersihan lingkungan dan personal higiene memiliki peran dalam menjaga integritas perlindungan kulit. Sanitasi dan air bersih yang buruk memaksa masyarakat hidup dalam kondisi yang kurang higienis, yang dapat mempermudah masuknya bakteri melalui luka kecil pada kulit atau kontak dengan permukaan yang tercemar. Variabel ini digunakan untuk melihat apakah kegagalan infrastruktur dasar di suatu wilayah berkontribusi pada langgengnya penyakit kusta.

5. Jumlah fasilitas kesehatan digunakan untuk merepresentasikan kemudahan akses masyarakat terhadap layanan diagnosis dan pengobatan Multi-Drug Therapy (MDT). Secara teoritis, keberadaan faskes yang cukup seharusnya dapat memutus rantai penularan dengan cepat karena penderita segera ditangani. Variabel ini digunakan untuk menguji apakah ketersediaan sarana kesehatan di Jawa Tengah sudah cukup efektif dalam menjaring kasus kusta sebelum penderita menularkan infeksinya lebih luas ke komunitas.

6. Variabel rata-rata lama sekolah dipilih sebagai indikator tingkat literasi kesehatan dan kapasitas individu dalam mengenali risiko penyakit. Kusta memiliki tantangan unik karena tidak menimbulkan rasa sakit di awal infeksi, sehingga pengetahuan yang didapat dari pendidikan menjadi sangat krusial agar masyarakat tidak mengabaikan gejala bercak mati rasa. Variabel ini penting untuk membuktikan apakah tingkat pendidikan berpengaruh terhadap kecepatan penemuan kasus dan pencegahan kecacatan tingkat dua yang sering terjadi akibat keterlambatan diagnosa.

Variabel-variabel tersebut dianalisis untuk mengetahui apakah variabel independen dari aspek kondisi lingkungan, sosial ekonomi, dan pendidikan tersebut berpengaruh terhadap jumlah kasus baru Kusta di Jawa Tengah. Dengan mengetahui variabel yang berpengaruh, diharapkan hasil penelitian ini dapat memberikan gambaran faktor risiko utama yang berkontribusi terhadap tingginya angka jumlah kasus baru kusta di wilayah tersebut.

Kemudian digunakan juga variabel bantu untuk mencari ukuran frekuensi dari penyakit Kusta di Provinsi Jawa Tengah maupun untuk memperkaya interpretasi epidemiologi dari penyakit, adapun variabel tersebut antara lain:

1. Jumlah penduduk yang digunakan untuk menghitung ukuran frekuensi prevalensi sesuai dengan rumus …(1).
2. Jumlah kasus kusta tercatat yang digunakan untuk menghitung ukuran frekuensi prevalensi sesuai dengan rumus …(1)
3. Jumlah kasus kusta basah yang digunakan untuk menghitung ukuran frekuensi prevalensi kusta basah sesuai dengan rumus …(1)
4. Jumlah kasus kusta kering yang digunakan untuk menghitung ukuran frekuensi prevalensi kusta kering sesuai dengan rumus …(1)
5. Jumlah kasus baru kusta cacat tingkat 0 untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.
6. Jumlah kasus baru kusta cacat tingkat 2 untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.
7. Jumlah pasien selesai beobat di tahun 2024 untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.
8. Persentase pasien kusta basah menyelesaikan pengobatan tepat waktu untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.
9. Jumlah kasus kusta kering anak untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.
10. Jumlah kasus kusta kering dewasa untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.
11. Jumlah kasus kusta basah anak untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.
12. Jumlah kasus kusta basah dewasa untuk memperkaya gambaran mengenai kondisi penyakit kusta di Jawa Tengah.

Metode Analisis

Penelitian ini menggunakan pendekatan analisis deskriptif epidemiologi untuk memahami pola penyebaran dan faktor yang berhubungan dengan kejadian penyakit kusta di Provinsi Jawa Tengah. Analisis deskriptif dipilih karena sesuai untuk menggambarkan karakteristik penyakit berdasarkan waktu, tempat, dan populasi tanpa memanipulasi variabel. Analisis dilakukan melalui ukuran frekuensi untuk melihat seberapa besar penyebaran kasus kusta pada tiap wilayah dengan menghitung prevalensi dan insidensi, ini membantu menggambarkan tingkat kejadian dan keparahan penyakit di populasi.

Selain itu, digunakan statistika deskriptif untuk memperjelas distribusi dan karakteristik data, serta visualisasi peta dan grafik untuk memperlihatkan pola spasial antarwilayah. Dilakukan juga pemodelan dengan menggunakan distribusi binomial negatif untuk mengetahui bagaimana dan seberapa besar pengaruh variabel independen pada variabel dependen yang diteliti, dalam perjalanannya dilakukan juga

Alur Kerja

Penelitian ini dimulai dengan tahap mengkaji proses dan penyebab terjadinya penyakit kusta agar dapat mengetahui gambarannya berdasarkan model agent, host, dan environment. Kemudian dicari variabel-variabel yang dianggap berhubungan berdasarkan gambaran yang didapat, pengumpulan data-data sekunder dari variabel yang digunakan bersumber dari Badan Pusat Statistik (BPS) dan Open Data Jawa Tengah. Data yang dikumpulkan meliputi Jumlah Kasus Baru Kusta, Kepadatan Penduduk, Persentase Penduduk Miskin, Persentase Rumah Layak Huni, Persentase Rumah dengan Akses Sanitasi Layak, Persentase Rumah dengan Air Minum Layak, Jumlah Fasilitas Kesehatan, dan Rata-Rata Lama Sekolah per kabupaten/kota di Jawa Tengah. Selain itu dilakukan juga pengumpulan variabel-variabel bantu untuk memperkaya interpretasi gambaran kasus.

Selanjutnya dilakukan tahap pengolahan data, seperti pengecekan kesesuaian format, penyatuan nama wilayah, serta pembuatan variabel turunan agar data siap dianalisis. Setelah data bersih, dilakukan analisis statistika deskriptif untuk memberikan gambaran umum mengenai sebaran data, seperti nilai rata-rata, minimum, maksimum, dan persebaran kasus Kusta di tiap kabupaten/kota.

Kemudian dilanjutkan dengan perhitungan ukuran frekuensi untuk melihat pola penyebaran Kusta di wilayah Jawa Tengah. Ukuran frekuensi yang digunakan meliputi prevalensi, yang menunjukkan proporsi penduduk yang terjangkit Kusta dan angka insidensi, yang menggambarkan laju munculnya kasus baru.

Kemudian dilakukan visualisasi dan interpretasi hasil. Data yang telah dianalisis disajikan dalam bentuk grafik dan peta distribusi guna menggambarkan pola spasial antarwilayah. Dari hasil tersebut, dilakukan interpretasi untuk memahami kecenderungan, wilayah dengan kasus tinggi, serta faktor risiko yang mungkin berperan dalam penyebaran kasus Kusta di Jawa Tengah.

Terakhir, dilakukan pemodelan dengan regresi binomial negatif untuk mengetahui variabel mana yang berpengaruh signifikan terhadap kenaikan kasus baru Kusta di Provinsi Jawa Tengah, nantinya variabel yang berpengaruh tersebut dapat diberikan perhatian khusus untuk menurunkan kasus kusta di Jawa Tengah.

HASIL DAN PEMBAHASAN

Deskripsi Kasus

Jumlah kasus baru Kusta di Provinsi Jawa Tengah pada tahun 2024 terhitung cukup tinggi dengan Kabupaten Tegal sebagai kabutapen/kota dengan jumlah kasus baru tertinggi di angka sekitar 265 kasus, kemudian setelahnya ada Kabupaten Brebes di angka 195 kasus, dua kabupaten ini memiliki jumlah kasus baru yang sangat tinggi dan cukup jauh jika dibandingkan dengan kabutapen/kota lainnya di Provinsi Jawa Tengah.

Dapat dilihat dari boxplot di atas, mayoritas kabupaten/kota berada di kisaran 35 sampai 12 kasus jumlah DBD yang berarti utamanya sebaran data jumlah kasus baru Kusta di Provinsi Jawa Tengah berada pada rentang tersebut. Namun, dapat dilihat terdapat empat outlier dengan jumlah kasus yang sangat tinggi yakni Kab. Tegal (265 kasus), Kab. Brebes (195 kasus), Kab. Pekalongan (81 kasus), dan Kab. Pemalang (73 kasus).

Faktor Agent, Host, Environment

Kusta merupakan penyakit infeksi yang disebabkan oleh bakteri Mycobactrium leprae dan ditularkan melalui percikan ecairan hidung (droplets) atau kontak erat dalam waktu yang lama dengan penderita yang belum diobati. Penularan penyakit ini dapat dijelaskan menggunakan model agent host environment, di mana ketiga faktor tersebut saling berinteraksi dan berperan dalam terjadinya penyakit.

1. Agent (Agen penyebab), agen penyebab penyakit kusta adalah bakteri Mycobactrium leprae yang kemudian terbagi menjadi dua jenis kusta, yakni kusta basah (Multibacillary/MB) dan kusta kering (Paucibacillary/PB). Bakteri ini bersifat intraseluler obligat yang menyerang saraf tepi, kulit, dan saluran pernapasan atas, berinkubasi selama 2-5 tahun, kemudian ditularkan secara langsung melalui percikan cairan hidung (droplets) atau kontak erat dalam waktu yang lama dengan penderita yang belum diobati.

2. Host (Pejamu), host atau pejamu utama penyakit ini adalah manusia. Risiko tertinggi terdapat pada masyarakat yang tinggal di wilayah padat penduduk. Beberapa faktor yang mempengaruhi kerentanan tubuh host seperti status gizi, jenis kelamin, usia produktif, dan tingkat imun yang dipengaruhi oleh faktor genetik juga dapat mempengaruhi terjangkitnya host (manusia) dari terinfeksi penyakit kusta.

3. Environment (Lingkungan), faktor lingkungan berperan penting dalam penyebaran kusta, terutama dengan karakteristik alur infeksi dan penyebaran penyakit kusta yang telah disebutkan sebelumnya. Lingkungan dengan penampungan air tercemar, tempat yang lembab, dan kepadatan penduduk yang tinggi menjadi tempat yang ideal bagi perkembangbiakan bakteri. Pengetahuan yang rendah akan penyakit kusta dan stigma masyarakat juga memengaruhi proses deteksi dan pemulihan penyakit, penderita yang tidak segera diobati juga akan menyebarkan bakteri ke penderita yang lain. mempercepat perkembangbiakan nyamuk.

Ketiga faktor tersebut saling berhubungan dalam siklus penularan penyakit kusta. Upaya pengendalian penyakit perlu dilakukan lebih giat dan merata, serta peningkatan pemahaman tentang penyakit kusta agar penderita langsung mendapatkan pengobatan dan stigma tentang penyakit kusta tidak lagi ada.

Diagram Hubungan Faktor Risiko dan Penyakit

Berikut merupakan diagram dari model segitiga epidemiologi (Agent, Host, Environment) yang menjelaskan bagaimana faktor lingkungan, host (manusia), dan agent (virus dengue) saling berinteraksi hingga menyebabkan kasus Kusta.

Dari diagram di atas, dapat dilihat bahwa ketiga faktor berperan dan saling berhubungan terhadap terjadinya kasus penyakit kusta. Untuk Environment sendiri terbagi ke dalam beberapa variabel yang memungkinkan, antara lain:

1. Kepadatan Penduduk Tinggi & Penduduk Miskin: Kusta menular melalui droplet (percikan cairan hidung/mulut) saat kontak erat dan lama. Lingkungan yang padat dan kumuh mempersempit ruang gerak sehingga interaksi antarmanusia menjadi sangat intens, yang mempercepat transmisi bakteri.
2. Rumah Tidak Layak Huni dan Rumah dengan Akses Sanitasi Buruk: Bakteri kusta menyukai tempat yang lembap dan kurang sinar matahari. Rumah dengan sanitasi yang buruk akan membuat spot lembab di beberapa titik (akibat kebocoran dan sebagainya). Rumah yang lembab dengan ventilasi dan pencahayaan buruk memperpanjang masa hidup bakteri di lingkungan.
3. Air Minum Tercemar: Air merupakan reservoir alami bagi bakteri kusta, jika suatu sumber air tercemar bakteri kusta maka potensi masyarakat terjangkit akan semakin besar.
4. Rata-rata Lama Sekolah Rendah: Kurangnya pengetahuan akan penyakit kusta membuat masyarakat tidak mengenali gejala awal kusta (seperti bercak putih dan mati rasa), hal ini akan membuat pemerintah setempat kesulitan dalam mendeteksi kasus Kusta hingga akhirnya sampai pada kecacatan (kusta cacat tingkat 2). Akibatnya, penderita tidak segera mencari pengobatan dan penyebaran bakteri masih terus berjalan (karena penderita belum meminum obat).
5. Stigma Masyarakat: Stigma negatif membuat penderita merasa malu atau takut dikucilkan, sehingga mereka cenderung menyembunyikan penyakitnya. Hal ini menyebabkan bakteri terus menyebar di lingkungan tanpa terdeteksi dan terjadi infeksi yang lebih parah pada penderita.

Kemudian untuk faktor yang memengaruhi kerentanan host akan menentukan apakah seseorang yang terpapar bakteri akan benar-benar sakit atau tetap sehat, dapat dijelaskan sebagai berikut:

1. Tingkat Imun dan Gizi Buruk: Sekitar 95% manusia sebenarnya memiliki kekebalan alami terhadap kusta. Namun, gizi buruk dan sistem imun yang lemah membuat tubuh gagal melawan bakteri, sehingga infeksi berkembang menjadi penyakit.
2. Usia Produktif: Memiliki mobilitas tinggi, sehingga risiko terpapar dengan berbagai sumber penularan di luar rumah lebih besar.
3. Jenis Kelamin Laki-laki: Secara statistik, kusta seringkali lebih banyak ditemukan pada laki-laki, yang diduga berkaitan dengan faktor gaya hidup, mobilitas, dan paparan lingkungan kerja.

Keseluruhan variabel-variabel dalam diagram tersebut berujung pada peningkatan jumlah kasus kusta, sehingga perlu penanganan lebih lanjut dari variabel-variabel terkait agar jumlah kasus terinfeksi penyakit kusta dapat ditekan seminimal mungkin.

Ukuran Epidemiologi

Ukuran Frekuensi

Akibat keterbatasan data yang ada, ukuran frekuensi yang dapat dihitung pada analisis kali ini hanya insidensi dan prevalensi.

Insidensi Kusta (Tipe Basah dan Kering)

Hasil perhitungan menunjukkan bahwa Kabupaten Tegal memiliki insidensi kusta tertinggi di Jawa Tengah, yaitu sebesar 15,82 per 100.000 penduduk. Angka ini menunjukkan bahwa dari setiap 100.000 penduduk di Kabupaten Tegal, terdapat sekitar 16 orang yang baru terdeteksi sakit kusta selama tahun 2024. Selanjutnya, Kota Pekalongan (13,08) dan Kabupaten Brebes (9,44) juga menunjukkan insidensi tinggi, yang mencerminkan beban penularan aktif yang berat di wilayah tersebut. Sebaliknya, Kabupaten Temanggung mencatatkan angka insidensi terendah, yaitu hanya 0,12 per 100.000 penduduk, yang memperlihatkan perbedaan beban risiko penularan yang sangat ekstrem antar wilayah.

Prevalensi

Prevalensi (Tipe Basah dan Kering)

Hasil perhitungan menunjukkan bahwa Kabupaten Brebes memiliki prevalensi kusta tertinggi di Jawa Tengah, yaitu mencapai 15,49 per 100.000 penduduk. Angka ini memberikan gambaran bahwa terdapat sekitar 15 hingga 16 penderita kusta yang sedang dalam masa pengobatan di antara 100.000 penduduk di wilayah tersebut. Wilayah lain dengan prevalensi tinggi adalah Kabupaten Tegal (12,78) dan Kota Pekalongan (12,77), yang menunjukkan beban kerja fasilitas kesehatan yang lebih besar dalam menangani pasien dibandingkan wilayah seperti Kabupaten Temanggung dengan prevalensi hanya 0,25.

Prevalensi (Tipe Basah)

Hasil perhitungan menunjukkan bahwa prevalensi kusta tipe basah atau Multibasiler (MB) mendominasi secara signifikan, dengan angka tertinggi di Kabupaten Brebes sebesar 15,05 per 100.000 penduduk yang berarti terdapat sekitar 15 penderita kusta tipe basah yang sedang dalam masa pengobatan di antara 100.000 penduduk di wilayah tersebut. Angka ini mengindikasikan bahwa mayoritas penderita kusta di wilayah tersebut memiliki muatan bakteri yang tinggi, yang selanjutnya diikuti oleh Kabupaten Tegal (12,6) dan Kota Pekalongan (12,46). Tingginya prevalensi kusta basah ini menjadi perhatian serius karena penderita tipe MB merupakan sumber penularan utama yang dapat terus menyebarkan bakteri di lingkungan sekitarnya.

Prevalensi (Tipe Kering)

Hasil perhitungan menunjukkan bahwa prevalensi kusta tipe kering atau Pausibasiler (PB) di Jawa Tengah berada pada tingkat yang sangat rendah, dengan angka tertinggi ditemukan di Kota Tegal sebesar 0,7 per 100.000 penduduk yang berarti terdapat sekitar 1 penderita kusta tipe kering yang sedang dalam masa pengobatan di antara 100.000 penduduk di wilayah tersebut. Wilayah berikutnya dengan prevalensi PB tertinggi adalah Kabupaten Pekalongan sebesar 0,69 dan Kabupaten Brebes sebesar 0,44. Rendahnya angka ini dibandingkan dengan tipe MB mengindikasikan bahwa sebagian besar penderita kemungkinan baru terdeteksi setelah penyakitnya berlanjut ke tahap yang lebih menular, atau adanya keterlambatan dalam deteksi dini di masyarakat.

Ukuran Asosiasi

Akibat keterbatasan data, tidak ada ukuran asosiasi yang dapat dihitung dari data yang ada.

Visualisasi

Peta Jumlah Kasus Kusta (Distribusi Kasus)

Peta tersebut menunjukkan sebaran jumlah kasus tercatat di wilayah Jawa Tengah, di mana terdapat perbedaan intensitas kasus yang mencolok antara area pesisir utara dan wilayah lainnya. Berdasarkan gradasi warna dan ukuran simbol, wilayah Tegal dan area di sisi baratnya menjadi titik dengan beban kasus tertinggi, yang ditandai dengan lingkaran besar berwarna merah hingga oranye tua yang mengindikasikan angka di atas 150 hingga mendekati 300 kasus. Sementara itu, sebagian besar wilayah Jawa Tengah lainnya, terutama di bagian tengah dan selatan seperti sekitar Surakarta, Salatiga, dan Magelang, menunjukkan angka yang jauh lebih rendah (di bawah 50 kasus) dengan simbol titik kuning kecil yang tersebar cukup merata.

Peta tersebut menunjukkan bahwa beban kasus baru cacat tingkat 0 di Jawa Tengah terkonsentrasi di wilayah Tegal dengan angka mencapai sekitar 200 kasus (simbol merah), serta tren yang cukup tinggi di sepanjang pesisir utara dibandingkan wilayah tengah dan selatan yang mayoritas berada di bawah 50 kasus. Tingginya angka pada kategori cacat tingkat 0 ini secara klinis mengindikasikan bahwa upaya deteksi penyakit di wilayah-wilayah tersebut berjalan cukup cepat dan responsif, karena kasus telah berhasil ditemukan sejak dini sebelum penderita mengalami kecatatan fisik yang tampak. Kondisi ini menunjukkan efektivitas penemuan kasus baru di lapangan, meskipun terdapat kesenjangan beban kasus yang mencolok antar daerah di dalam provinsi tersebut.

Peta sebaran tersebut menunjukkan bahwa konsentrasi tertinggi kasus baru dengan cacat tingkat 2 di Jawa Tengah berada di wilayah Tegal, yang ditandai dengan lingkaran merah pekat paling besar dengan jumlah mencapai sekitar 20 kasus. Wilayah lain di sepanjang pesisir utara, seperti area Pekalongan dan wilayah di sebelah timur Semarang, juga menunjukkan beban kasus yang cukup signifikan dengan rentang antara 10 hingga 15 kasus (simbol oranye). Munculnya angka pada kategori cacat tingkat 2 ini mengindikasikan bahwa masih banyak penderita yang baru ditemukan atau melaporkan diri setelah penyakitnya berkembang ke tahap lanjut, sehingga sudah menimbulkan kecacatan fisik yang terlihat secara permanen. Sementara itu, wilayah bagian tengah dan selatan Jawa Tengah secara umum memiliki angka yang jauh lebih rendah, didominasi oleh titik kuning kecil yang menunjukkan jumlah kasus baru cacat tingkat 2 berada di kisaran 0 sampai 5 kasus.

Peta Jumlah Penduduk Provinsi Jawa Tengah

Peta Jumlah Penduduk tersebut menggambarkan kepadatan populasi di berbagai kabupaten/kota di wilayah Jawa Tengah. Berdasarkan gradasi warna dari kuning hingga merah pekat, wilayah dengan jumlah penduduk terbanyak (mencapai angka 1,8 juta hingga 2 juta jiwa) ditandai dengan lingkaran merah besar, seperti yang terlihat di wilayah Cilacap dan sekitarnya. Sementara itu, wilayah dengan simbol lingkaran oranye hingga kuning menunjukkan populasi yang lebih moderat hingga rendah, berkisar antara 200 ribu hingga 1,2 juta jiwa. Data ini berfungsi sebagai pembanding penting; misalnya, jika suatu daerah memiliki jumlah penduduk yang sangat besar namun jumlah kasus kesehatannya rendah, maka beban penyakit di sana relatif lebih terkendali dibandingkan daerah berpenduduk sedikit namun memiliki jumlah kasus yang tinggi.

Peta Insidensi Kusta

Peta tersebut menunjukkan bahwa tingkat kejadian kasus baru paling tinggi terkonsentrasi di wilayah Tegal dan Pekalongan, yang ditandai dengan simbol lingkaran merah besar dengan nilai 13 hingga 16. Angka ini memiliki arti terdapat sekitar 13 sampai 16 kasus baru per 100.000 penduduk di wilayah tersebut. Secara epidemiologi, angka ini menunjukkan tingkat risiko penularan yang tinggi dan beban penyakit yang berat di kedua daerah tersebut dibandingkan dengan wilayah Jawa Tengah lainnya yang mayoritas hanya memiliki angka insidensi di bawah 2. Munculnya angka insidensi yang tinggi ini mempertegas bahwa meskipun suatu daerah mungkin memiliki jumlah penduduk yang berbeda, proporsi temuan kasus baru di Tegal dan Pekalongan tetap yang paling signifikan secara statistik.

Peta Prevalensi Kusta

Peta tersebut menunjukkan bahwa angka prevalensi tertinggi terkonsentrasi di wilayah pesisir utara, khususnya di area Brebes dan Tegal yang ditandai dengan simbol lingkaran merah hingga oranye tua dengan nilai antara 13 hingga 15. Angka 13 hingga 15 ini mengindikasikan bahwa terdapat sekitar 13 hingga 15 kasus per 100.000 penduduk yang sedang menderita atau dalam masa pengobatan di wilayah tersebut pada satu titik waktu tertentu. Sebaliknya, sebagian besar wilayah Jawa Tengah lainnya, terutama di bagian tengah dan selatan, didominasi oleh titik kuning kecil dengan angka prevalensi rendah di bawah 2, yang menunjukkan bahwa jumlah penderita aktif di wilayah tersebut jauh lebih sedikit secara proporsional terhadap jumlah penduduknya.

Peta prevalensi menunjukkan beban penyakit yang didominasi secara mutlak oleh tipe kusta basah, terutama di wilayah koridor utara seperti Brebes, Tegal, dan Pekalongan. Wilayah Brebes dan Tegal dkk tercatat memiliki kasus MB tertinggi sekaligus kasus PB tertinggi di Jawa Tengah, namun angka prevalensinya sangat timpang; sebagai contoh, di Kabupaten Brebes angka MB mencapai 15, sedangkan PB-nya hanya 0,44. Ketimpangan yang drastis meski sama-sama memiliki kasus kusta MB dan PB tertinggi ini menunjukkan bahwa kusta basah (MB) memang sangat banyak dan jauh lebih mendominasi dibandingkan kusta kering (PB) di wilayah tersebut. Hal ini kemungkinan besar disebabkan karena tipe MB memiliki daya tular yang jauh lebih tinggi dan lebih mudah menularkan bakteri kepada orang lain di lingkungan sekitarnya dibandingkan tipe PB, sehingga menciptakan konsentrasi kasus yang jauh lebih besar di daerah-daerah tersebut meski sama-sama menempati urutan tertinggi.

Desain Studi Epidemiologi

Desain Studi

Dilakukan penelitian dengan menggunakan desain cross-sectional (potong lintang) di mana unit analisis adalah kabupaten/kota di Provinsi Jawa Tengah pada tahun 2024. Penelitian ini dilakukan untuk mengetahui hubungan antara variabel Kepadatan Penduduk, Persentase Penduduk Miskin, Persentase Rumah Layak Huni, Persentase Rumah dengan Akses Sanitasi Layak, Persentase Rumah dengan Air Minum Layak, Jumlah Fasilitas Kesehatan, dan Rata-Rata Lama Sekolah terhadap jumlah kasus baru kusta di Provinsi Jawa Tengah.

Variabel Utama

Variabel utama (respons) yang digunakan adalah variabel Jumlah Kasus Baru Kusta per kabupaten/kota di Provinsi Jawa Tengah dengan total 35 kabupaten/kota. Sedangkan variabel lainnya yang kemudian menjadi variabel prediktor adalah variabel Kepadatan Penduduk, Persentase Penduduk Miskin, Persentase Rumah Layak Huni, Persentase Rumah dengan Akses Sanitasi Layak, Persentase Rumah dengan Air Minum Layak, Jumlah Fasilitas Kesehatan, dan Rata-Rata Lama Sekolah terhadap jumlah kasus baru kusta di Provinsi Jawa Tengah. Bentuk datanya antara lain sebagai berikut

## # A tibble: 6 × 9
##   Wilayah      `Jumlah Kasus Baru` Persentase Penduduk …¹ Persentase Rumah Lay…²
##   <chr>                      <dbl>                  <dbl>                  <dbl>
## 1 Kab. Cilacap                  30                   10.7                   60.7
## 2 Kab. Banyum…                  12                   12.0                   67  
## 3 Kab. Purbal…                  20                   14.2                   63.3
## 4 Kab. Banjar…                   6                   14.7                   42.8
## 5 Kab. Kebumen                  24                   15.7                   74.1
## 6 Kab. Purwor…                   8                   10.9                   72.2
## # ℹ abbreviated names: ¹​`Persentase Penduduk Miskin`,
## #   ²​`Persentase Rumah Layak Huni`
## # ℹ 5 more variables: `Persentase Akses Sanitasi Layak` <dbl>,
## #   `Persentase Akses Air Minum Layak` <dbl>,
## #   `Jumlah Fasilitas Kesehatan` <dbl>, `Kepadatan Penduduk (per km^2)` <dbl>,
## #   `Rata-Rata Lama Sekolah (IPM) (Tahun)` <dbl>

Semua data yang digunakan bersumber dari Badan Pusat Statistik (BPS) dan Open Data Jawa Tengah dalam periode waktu yang sama (tahun 2024), agar hasil analisis tetap konsisten dan dapat dibandingkan antarwilayah.

Populasi dan Sampling

Populasi penelitian mencakup 35 kabupaten/kota di Provinsi Jawa Tengah, di mana setiap wilayah menjadi satu unit analisis. Karena semua wilayah yang ada dimasukkan dalam penelitian, metode yang digunakan adalah total sampling, sehingga gambaran yang dihasilkan mewakili keseluruhan provinsi tanpa harus melakukan generalisasi dari sampel kecil. Semua data diperoleh dari sumber resmi pemerintah pada periode yang sama agar hasilnya tetap konsisten.

Potensi Bias

1. Ecological Fallacy (Ekologi Keliru) Potensi bias pertama dan yang paling utama dalam penelitian ini adalah ecological fallacy, yaitu kesalahan logika dalam mengambil kesimpulan tentang individu berdasarkan data tingkat populasi. Karena data yang dianalisis adalah agregat tingkat kabupaten/kota di Jawa Tengah, hubungan signifikan antara Rata-Rata Lama Sekolah (RLS) dan jumlah kasus kusta tidak serta-merta berarti bahwa setiap individu yang memiliki pendidikan rendah pasti akan menderita kusta. Hubungan ini hanya berlaku pada tingkat wilayah, di mana wilayah dengan pendidikan rata-rata rendah cenderung memiliki beban kasus yang lebih tinggi, namun faktor risiko di tingkat personal bisa saja berbeda.
2. Underreporting Bias (Bias Laporan Rendah) Penelitian ini sangat bergantung pada data sekunder yang berasal dari laporan fasilitas kesehatan, sehingga terdapat potensi underreporting bias atau fenomena gunung es. Rendahnya angka insidensi di beberapa wilayah seperti Kabupaten Temanggung (0,12) belum tentu mencerminkan kondisi bebas kusta yang sesungguhnya, melainkan bisa disebabkan oleh lemahnya sistem surveilans, adanya stigma sosial yang membuat penderita enggan memeriksakan diri, atau kurangnya tenaga kesehatan yang terlatih dalam deteksi dini. Akibatnya, model Regresi Binomial Negatif mungkin mengestimasi koefisien berdasarkan angka yang lebih rendah dari kenyataan di lapangan (kasus yang tidak terdeteksi).
3. Omitted Variable Bias (Bias Variabel yang Terlewat) Masih terdapat potensi omitted variable bias di mana ada faktor-faktor krusial lain yang tidak masuk dalam model karena keterbatasan data. Faktor-faktor seperti status gizi spesifik (mikronutrien), riwayat vaksinasi BCG, faktor genetik populasi lokal, atau mobilitas penduduk antar provinsi tidak dihitung dalam penelitian ini.
4. Bias Temporal (Cross-sectional) Penelitian ini menggunakan desain cross-sectional pada tahun 2024, yang menimbulkan potensi bias dalam menentukan hubungan sebab-akibat yang pasti. Sebagai contoh, hubungan antara Penduduk Miskin dan kusta bersifat dua arah (kemiskinan meningkatkan risiko kusta karena gizi buruk, namun penyakit kusta juga dapat menyebabkan seseorang jatuh miskin karena kecacatan dan hilangnya pekerjaan). Karena data diambil pada satu titik waktu yang sama, penelitian ini tidak bisa membuktikan secara absolut mana yang menjadi penyebab awal (kausalitas).

Pemodelan Penyakit

Pemodelan dilakukan menggunakan variabel-variabel yang telah disebutkan pada desain studi sebelumnya. Karena variabel Y yang digunakan (jumlah kasus baru kusta) merupakan data diskrit atau hasil hitung (data count) yang tidak mungkin bernilai negatif, sehingga tidak bisa menggunakan regresi linear biasa. Model awal yang paling standar untuk data count adalah Regresi Poisson yang bekerja dengan asumsi dasar bahwa nilai rata-rata (mean) harus sama dengan nilai variansnya (equidispersion).
Namun, dalam data kesehatan masyarakat seperti kusta di Jawa Tengah, asumsi ini jarang terpenuhi karena penyebaran kasus antar wilayah cenderung sangat timpang, di mana beberapa wilayah memiliki kasus yang sangat melonjak sementara wilayah lain sangat rendah. Jika setelah diuji overdispersi didapat nilai varians tidak sama dengan rata-ratanya, maka tidak dapat menggunakan model Regresi Poisson dan dialternatifkan menggunakan model Regresi Binomial Negatif. Regresi Binomial Negatif dipilih karena memiliki parameter tambahan yang disebut parameter dispersi (\(\alpha\)) untuk mengakomodasi keragaman data yang ekstrem tersebut. Model ini jauh lebih fleksibel dalam menangani data yang memiliki nilai nol yang banyak atau lonjakan kasus yang sangat tinggi di beberapa titik.

Uji Asumsi

Terlebih dahulu membuat model awal standar, model Regresi Poisson untuk diuji asumsi sebagai berikut:

model_pois <- glm(
  Jumlah_Kasus_Baru ~
    Persentase_Penduduk_Miskin +
    Persentase_Rumah_Layak_Huni +
    Persentase_Akses_Sanitasi_Layak +
    Persentase_Akses_Air_Minum_Layak +
    Jumlah_Fasilitas_Kesehatan +
    Kepadatan_Penduduk +
    Rata_Rata_Lama_Sekolah,
  family = poisson(link = "log"),
  data = data
)

Uji Multikolinearitas

Setelah menentukan model dasar, dilakukan pengujian multikolinearitas untuk memastikan bahwa tidak terdapat hubungan linear yang sangat kuat antar variabel independen yang dapat mengganggu stabilitas hasil estimasi model.

vif(model_pois)

##       Persentase_Penduduk_Miskin      Persentase_Rumah_Layak_Huni 
##                         2.049445                         3.575208 
##  Persentase_Akses_Sanitasi_Layak Persentase_Akses_Air_Minum_Layak 
##                         2.698691                         2.470035 
##       Jumlah_Fasilitas_Kesehatan               Kepadatan_Penduduk 
##                         1.604797                         1.943696 
##           Rata_Rata_Lama_Sekolah 
##                         2.384187

Berdasarkan uji multikolinearitas terhadap model, seluruh variabel independen menunjukkan nilai VIF yang berada jauh di bawah ambang batas 10. Nilai VIF terendah ditemukan pada variabel Jumlah Fasilitas Kesehatan sebesar 1,60, diikuti oleh Kepadatan Penduduk (1,94) dan Persentase Penduduk Miskin (2,04). Variabel lainnya seperti Rata-Rata Lama Sekolah (2,38), Persentase Akses Air Minum Layak (2,47), serta Persentase Akses Sanitasi Layak (2,69) juga menunjukkan angka yang rendah. Nilai VIF tertinggi hanya sebesar 3,57 yang ditemukan pada variabel Persentase Rumah Layak Huni, yang masih bisa dikategorikan aman.

Hasil tersebut memberikan kesimpulan logis bahwa tidak terdapat masalah multikolinearitas dalam penelitian ini. Karena tidak ada hubungan korelasi yang ekstrem antar variabel bebas, maka setiap variabel dapat diinterpretasikan secara independen. Hal ini memperkuat validitas hasil pada tahap selanjutnya, di mana pengaruh dari faktor ekonomi, pendidikan, dan lingkungan terhadap kasus kusta dapat dibedakan secara jelas tanpa adanya bias statistik antar variabel.

Uji Overdispersi

Kemudian dilakukan uji Overdispersi untuk memeriksa pemenuhan asumsi dasar pada model Regresi Poisson, hasil ini akan menentukan metode apa yang sebaiknya digunakan antara metode awal standar Regresi Poisson dan Regresi Binomial Negatif.

residual_deviance <- model_pois$deviance
df_residual <- model_pois$df.residual
dispersion_ratio <- residual_deviance / df_residual
print(dispersion_ratio)

## [1] 16.49831

Berdasarkan hasil pengujian terhadap model Poisson, ditemukan indikasi overdispersi yang sangat kuat. Hal ini terlihat dari nilai dispersion ratio sebesar 16,49 yang dihitung melalui perbandingan residual deviance dan degree of freedom.

dispersiontest(model_pois)

## 
##  Overdispersion test
## 
## data:  model_pois
## z = 3.1995, p-value = 0.0006884
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   14.93393

Hasil ini diperkuat dengan uji formal Cameron-Trivedi melalui fungsi dispersiontest() yang menghasilkan nilai estimasi dispersi sebesar 14,93 dengan n nilai \(p-value = 0,00068\) (\(p < 0,05\)). Karena kedua indikator secara konsisten menunjukkan bahwa nilai dispersi jauh lebih besar dari 1, maka asumsi equidispersion pada model Poisson dinyatakan tidak terpenuhi, sehingga penggunaan model Regresi Binomial Negatif menjadi mutlak diperlukan.

Model Regresi Binomial Negatif

Model Penuh

model_nb_full <- glm.nb(
  Jumlah_Kasus_Baru ~
    Persentase_Penduduk_Miskin +
    Persentase_Rumah_Layak_Huni +
    Persentase_Akses_Sanitasi_Layak +
    Persentase_Akses_Air_Minum_Layak +
    Jumlah_Fasilitas_Kesehatan +
    Kepadatan_Penduduk +
    Rata_Rata_Lama_Sekolah,
  data = data
)
summary(model_nb_full)

## 
## Call:
## glm.nb(formula = Jumlah_Kasus_Baru ~ Persentase_Penduduk_Miskin + 
##     Persentase_Rumah_Layak_Huni + Persentase_Akses_Sanitasi_Layak + 
##     Persentase_Akses_Air_Minum_Layak + Jumlah_Fasilitas_Kesehatan + 
##     Kepadatan_Penduduk + Rata_Rata_Lama_Sekolah, data = data, 
##     init.theta = 2.248684369, link = log)
## 
## Coefficients:
##                                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       1.191e+01  4.649e+00   2.561  0.01043 *  
## Persentase_Penduduk_Miskin       -1.767e-01  5.721e-02  -3.089  0.00201 ** 
## Persentase_Rumah_Layak_Huni       2.939e-02  1.855e-02   1.584  0.11314    
## Persentase_Akses_Sanitasi_Layak   2.878e-02  1.959e-02   1.469  0.14182    
## Persentase_Akses_Air_Minum_Layak -3.645e-02  4.678e-02  -0.779  0.43596    
## Jumlah_Fasilitas_Kesehatan        2.339e-04  2.594e-04   0.902  0.36727    
## Kepadatan_Penduduk                2.132e-04  7.954e-05   2.680  0.00736 ** 
## Rata_Rata_Lama_Sekolah           -1.051e+00  1.778e-01  -5.912 3.38e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2487) family taken to be 1)
## 
##     Null deviance: 91.904  on 34  degrees of freedom
## Residual deviance: 37.106  on 27  degrees of freedom
## AIC: 306.42
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.249 
##           Std. Err.:  0.582 
## Warning while fitting theta: alternation limit reached 
## 
##  2 x log-likelihood:  -288.418

Dari hasil pemodelan regresi binomial negatif dengan seluruh variabel yang digunakan, didapatkan bahwa variabel Kepadatan Penduduk, Rata-Rata Lama Sekolah, dan Persentase Penduduk Miskin memiliki pengaruh yang signifikan secara statistik terhadap jumlah kasus baru kusta (p-value < 0,05). Sementara itu, faktor lingkungan fisik seperti akses sanitasi, air minum, dan rumah layak huni tidak menunjukkan pengaruh signifikan dalam model ini.

Model dengan Variabel Signifikan

Kemudian dicoba melakukan pemodelan hanya dengan variabel signifikan saja dari model full, yakni Kepadatan Penduduk, Rata-Rata Lama Sekolah, dan Persentase Penduduk Miskin. Adapun hasilnya adalah sebagai berikut:

model_nb_sig <- glm.nb(
  Jumlah_Kasus_Baru ~
    Persentase_Penduduk_Miskin +
    Kepadatan_Penduduk +
    Rata_Rata_Lama_Sekolah,
  data = data
)
summary(model_nb_sig)

## 
## Call:
## glm.nb(formula = Jumlah_Kasus_Baru ~ Persentase_Penduduk_Miskin + 
##     Kepadatan_Penduduk + Rata_Rata_Lama_Sekolah, data = data, 
##     init.theta = 1.530575644, link = log)
## 
## Coefficients:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 1.086e+01  1.771e+00   6.132 8.69e-10 ***
## Persentase_Penduduk_Miskin -1.190e-01  6.087e-02  -1.955   0.0506 .  
## Kepadatan_Penduduk          1.365e-04  8.259e-05   1.653   0.0983 .  
## Rata_Rata_Lama_Sekolah     -7.869e-01  1.653e-01  -4.759 1.95e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.5306) family taken to be 1)
## 
##     Null deviance: 64.920  on 34  degrees of freedom
## Residual deviance: 38.066  on 31  degrees of freedom
## AIC: 312.08
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.531 
##           Std. Err.:  0.369 
## 
##  2 x log-likelihood:  -302.077

Setelah dilakukan penyederhanaan model dengan hanya menyertakan variabel yang sebelumnya signifikan, hasil pemodelan menunjukkan perubahan pada tingkat signifikansi dan kekuatan pengaruh variabel. Dapat dilihat model full memiliki nilai AIC yang lebih baik dibandingkan model dengan variabel signifikan (306.42 vs 312.08), ini menunjukkan bahwa model pertama (model lengkap) sebenarnya memberikan performa yang sedikit lebih baik dalam menjelaskan data. Hal ini menandakan bahwa variabel lingkungan (sanitasi, air, rumah) tetap memberikan kontribusi informasi bagi model meskipun secara individu mereka tidak signifikan.

Model StepAIC

Karena itu, dilakukan percobaan untuk mencari kombinasi variabel terbaik dengan menggunakan metode Stepwise AIC, adapun hasilnya adalah sebagai berikut:

model_nb_step <- stepAIC(model_nb_full, direction = "both")

## Start:  AIC=304.42
## Jumlah_Kasus_Baru ~ Persentase_Penduduk_Miskin + Persentase_Rumah_Layak_Huni + 
##     Persentase_Akses_Sanitasi_Layak + Persentase_Akses_Air_Minum_Layak + 
##     Jumlah_Fasilitas_Kesehatan + Kepadatan_Penduduk + Rata_Rata_Lama_Sekolah

##                                    Df    AIC
## - Persentase_Akses_Air_Minum_Layak  1 303.08
## - Jumlah_Fasilitas_Kesehatan        1 303.19
## <none>                                304.42
## - Persentase_Akses_Sanitasi_Layak   1 304.50
## - Persentase_Rumah_Layak_Huni       1 304.99
## - Kepadatan_Penduduk                1 308.94
## - Persentase_Penduduk_Miskin        1 312.87
## - Rata_Rata_Lama_Sekolah            1 329.80

## 
## Step:  AIC=303.08
## Jumlah_Kasus_Baru ~ Persentase_Penduduk_Miskin + Persentase_Rumah_Layak_Huni + 
##     Persentase_Akses_Sanitasi_Layak + Jumlah_Fasilitas_Kesehatan + 
##     Kepadatan_Penduduk + Rata_Rata_Lama_Sekolah

##                                    Df    AIC
## - Jumlah_Fasilitas_Kesehatan        1 301.61
## - Persentase_Rumah_Layak_Huni       1 303.01
## <none>                                303.08
## - Persentase_Akses_Sanitasi_Layak   1 303.16
## + Persentase_Akses_Air_Minum_Layak  1 304.42
## - Kepadatan_Penduduk                1 306.94
## - Persentase_Penduduk_Miskin        1 310.88
## - Rata_Rata_Lama_Sekolah            1 328.36
## 
## Step:  AIC=301.61
## Jumlah_Kasus_Baru ~ Persentase_Penduduk_Miskin + Persentase_Rumah_Layak_Huni + 
##     Persentase_Akses_Sanitasi_Layak + Kepadatan_Penduduk + Rata_Rata_Lama_Sekolah

##                                    Df    AIC
## - Persentase_Rumah_Layak_Huni       1 301.07
## <none>                                301.61
## - Persentase_Akses_Sanitasi_Layak   1 303.06
## + Jumlah_Fasilitas_Kesehatan        1 303.08
## + Persentase_Akses_Air_Minum_Layak  1 303.19
## - Kepadatan_Penduduk                1 305.06
## - Persentase_Penduduk_Miskin        1 310.27
## - Rata_Rata_Lama_Sekolah            1 331.28
## 
## Step:  AIC=301.07
## Jumlah_Kasus_Baru ~ Persentase_Penduduk_Miskin + Persentase_Akses_Sanitasi_Layak + 
##     Kepadatan_Penduduk + Rata_Rata_Lama_Sekolah
## 
##                                    Df    AIC
## <none>                                301.07
## + Persentase_Rumah_Layak_Huni       1 301.61
## + Jumlah_Fasilitas_Kesehatan        1 303.01
## + Persentase_Akses_Air_Minum_Layak  1 303.05
## - Kepadatan_Penduduk                1 304.22
## - Persentase_Penduduk_Miskin        1 309.61
## - Persentase_Akses_Sanitasi_Layak   1 310.08
## - Rata_Rata_Lama_Sekolah            1 330.15

summary(model_nb_step)

## 
## Call:
## glm.nb(formula = Jumlah_Kasus_Baru ~ Persentase_Penduduk_Miskin + 
##     Persentase_Akses_Sanitasi_Layak + Kepadatan_Penduduk + Rata_Rata_Lama_Sekolah, 
##     data = data, init.theta = 2.070701281, link = log)
## 
## Coefficients:
##                                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      9.011e+00  1.634e+00   5.515 3.48e-08 ***
## Persentase_Penduduk_Miskin      -1.686e-01  5.406e-02  -3.120 0.001810 ** 
## Persentase_Akses_Sanitasi_Layak  5.251e-02  1.372e-02   3.828 0.000129 ***
## Kepadatan_Penduduk               1.844e-04  7.791e-05   2.367 0.017933 *  
## Rata_Rata_Lama_Sekolah          -1.072e+00  1.714e-01  -6.256 3.96e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.0707) family taken to be 1)
## 
##     Null deviance: 85.379  on 34  degrees of freedom
## Residual deviance: 37.107  on 30  degrees of freedom
## AIC: 303.07
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.071 
##           Std. Err.:  0.525 
## 
##  2 x log-likelihood:  -291.066

Dari hasil percobaan beberapa kombinasi variabel, didapatkan bahwa kombinasi dengan variabel Persentase Penduduk Miskin, Persentase Akses Sanitasi Layak, Kepadatan Penduduk, dan Rata-Rata Lama Sekolah. Untuk variabel akses sanitasi layak yang sebelumnya tidak signifikan saat di model full, diikutsertakan pada model ini karena menghasilkan nilai AIC terkecil. Nilai theta yang lebih besar dari 1 juga menunjukkan bahwa model ini berhasil menangani masalah overdispersion (keragaman data yang sangat tinggi) yang tidak bisa ditangani oleh regresi Poisson biasa.

Perbandingan Model dan Evaluasi Performa

Kemudian dilakukan perbandingan dari ketiga model yang telah dibuat untuk mengetahui model mana yang merupakan model terbaik. Evaluasi performa dilakukan dengan menggunakan metrik AIC, BIC, dan Adjusted McFadden \(R^2\) dengan hasil sebagai berikut:

## fitting null model for pseudo-r2
## fitting null model for pseudo-r2
## fitting null model for pseudo-r2

##                       AIC      BIC Adj_R2_McFadden
## Full_Model       306.4179 320.4161      0.10625673
## Significant_Only 312.0768 319.8535      0.06393098
## StepAIC_Model    303.0657 312.3978      0.09805171

Berdasarkan nilai AIC dan BIC terkecil, model terbaik adalah Model StepAIC dengan AIC = 303,07 dan BIC = 312,39. Model ini memiliki nilai AIC dan BIC terkecil yang menunjukkan keseimbangan terbaik antara kesesuaian model (goodness of fit) dan kompleksitas model. Meskipun Full Model memiliki McFadden R² yang sedikit lebih tinggi dari model terpilih (0,106 vs 0,098), selisih nilai tersebut relatif kecil. Hal ini menunjukkan bahwa penambahan variabel pada Full Model hanya memberikan peningkatan kemampuan penjelasan yang terbatas, sementara nilai AIC dan BIC yang lebih rendah pada model Stepwise mengindikasikan bahwa model ini lebih efisien dalam menjelaskan data.

Model Terbaik

Model terbaik terpilih adalah model StepAIC dengan variabel Persentase Penduduk Miskin, Persentase Akses Sanitasi Layak, Kepadatan Penduduk, dan Rata-Rata Lama Sekolah, karena model ini menggunakan link function berupa log, maka bentuk persamaannya adalah:

\[ ln(/mu) =9,0107-0,1686(PendudukMiskin)+0,0525(AksesSanitasiLayak)+0,0002(KepadatanPenduduk)-1,0720(RataRataLamaSekolah) \]

• Variabel Persentase Penduduk Miskin memiliki koefisien negatif sebesar -0,1686 dan berpengaruh signifikan terhadap model. Artinya, setiap peningkatan satu persen penduduk miskin, jumlah kasus kusta yang terlaporkan diprediksi mengalami penurunan sebesar 15,5%. Secara teori, kemiskinan berkaitan dengan lingkungan kumuh yang mempercepat transmisi, namun hasil negatif ini menunjukkan adanya pengaruh faktor lain kemungkinan seperti stigma masyarakat yang menyebabkan penderita tidak terdeteksi oleh sistem kesehatan (under-reporting).

• Variabel Persentase Akses Sanitasi Layak memiliki koefisien positif sebesar 0,0525 yang signifikan. Artinya, setiap peningkatan satu persen akses sanitasi layak, jumlah kasus terdeteksi diprediksi meningkat sebesar 5,4%. Meskipun hasilnya berbanding terbalik dengan logika di awal, terdapat kemungkinan bahwa wilayah dengan sanitasi layak umumnya memiliki bangunan yang rapat dan berpotensi menghalangi masuknya sinar matahari dan menciptakan area mikro yang lembap. Mengingat bakteri kusta sangat menyukai tempat lembap dan air merupakan reservoir alami baginya, lingkungan seperti ini justru dapat memperpanjang masa hidup bakteri di luar tubuh manusia. Selain itu, wilayah dengan sanitasi baik biasanya lebih aktif dalam pelaporan kesehatan, sehingga kasus yang sebelumnya tersembunyi menjadi lebih mudah ditemukan.

• Variabel Kepadatan Penduduk memiliki koefisien positif sebesar 0,00018 yang signifikan. Artinya, setiap kenaikan satu jiwa per km², jumlah kasus kusta diprediksi meningkat sebesar 0,018%. Hal ini sangat berkaitan dengan cara penularan kusta yang melalui droplet (percikan cairan hidung/mulut) saat kontak erat dan lama. Lingkungan yang memiliki kepadatan penduduk tinggi mempersempit ruang gerak sehingga interaksi antarmanusia menjadi sangat intens. Kondisi ini mempermudah transmisi bakteri dari penderita ke orang sehat, sehingga mempercepat penyebaran penyakit di wilayah tersebut.

• Variabel Rata-Rata Lama Sekolah (RLS) memiliki koefisien negatif dan memiliki pengaruh paling signifikan (-1,0720). Artinya, setiap peningkatan satu unit rata-rata lama sekolah, maka jumlah kasus baru kusta diprediksi akan mengalami penurunan sebesar 65,8%. Hal ini menunjukkan bahwa tingkat pendidikan yang lebih tinggi meningkatkan kesadaran masyarakat terhadap penyakit kusta baik dalam hal pencegahan maupun deteksi dini (mengenali bercak dan mati rasa), dan kemauan untuk segera berobat sebelum terjadi penularan lebih lanjut.

KESIMPULAN DAN SARAN

Kesimpulan

Gambaran kasus kusta di Provinsi Jawa Tengah menunjukkan tren yang masih fluktuatif dengan persebaran yang tidak merata, di mana beban kasus tertinggi masih mendominasi wilayah pesisir utara (Pantura) khususnya daerah Tegal, Brebes, Pekalongan dan sebagian wilayah timur. Secara epidemiologi, proses terjadinya penyakit ini merupakan hasil interaksi kompleks antara bakteri Mycobacterium leprae (agent), tingkat imunitas dan mobilitas yang memengaruhi kerentanan manusia (host), dan kondisi lingkungan fisik rumah yang lembap atau padat penghuni (environment), angka rata-rata lama sekolah sebagai indikator tingkat literasi ataupun pengetahuan masyarakat juga memengaruhi angka munculnya kasus baru. Dari sisi spasial, wilayah Kabupaten Brebes, Tegal, dan sekitarnya teridentifikasi sebagai prioritas utama yang membutuhkan perhatian khusus karena statusnya sebagai wilayah dengan transmisi tinggi. Kenaikan kasus baru di wilayah-wilayah tersebut secara signifikan dipengaruhi oleh faktor variabel rata-rata lama sekolah, kepadatan penduduk, persentase penduduk miskin, dan rumah dengan akses sanitasi layak yang memengaruhi angka kasus baru kusta di Provinsi Jawa Tengah.

Saran

Disarankan kepada Pemerintah Provinsi Jawa Tengah untuk memprioritaskan alokasi sumber daya dan program pemeriksaan kusta secara aktif pada wilayah-wilayah di zona prioritas tinggi seperti Tegal, Brebes, dan Pekalongan. Intervensi harus difokuskan pada perbaikan variabel signifikan hasil pemodelan, yaitu dengan meningkatkan kualitas lingkungan perumahan melalui program bedah rumah sehat dan memperluas cakupan deteksi dini melalui pelacakan kontak erat secara intensif. Selain itu, edukasi mengenai tanda-tanda awal penyakit kusta, penyebab terjadinya penyakit kusta, dan penghapusan stigma di masyarakat perlu diperkuat untuk memastikan setiap penderita mendapatkan pengobatan tepat waktu guna mencegah keterlambatan deteksi yang menyebabkan penyebaran penyakit lebih luas dan kecacatan permanen.

REFERENSI

Mungroo, M., Khan, N., & Siddiqui, R. (2020). Mycobacterium leprae: Pathogenesis, diagnosis, and treatment options.. Microbial pathogenesis, 104475 . https://doi.org/10.1016/j.micpath.2020.104475.

Dharmawan, Y., Fuady, A., Korfage, I., & Richardus, J. (2021). Individual and community factors determining delayed leprosy case detection: A systematic review. PLoS Neglected Tropical Diseases, 15. https://doi.org/10.1371/journal.pntd.0009651.

Ulfa, Y., Soleh, A., & Sartono, B. (2021). Handling of Overdispersion in the Poisson Regression Model with Negative Binomial for the Number of New Cases of Leprosy in Java. , 5, 1-13. https://doi.org/10.29244/ijsa.v5i1p1-13.

Saifudin, T., Rifada, M., Makhbubah, K., & Ramadhanty, D. (2025). LEPROSY CASE MODELING IN EAST JAVA USING SPATIAL REGRESSION WITH QUEEN CONTIGUITY WEIGHTING. BAREKENG: Jurnal Ilmu Matematika dan Terapan. https://doi.org/10.30598/barekengvol19iss3pp2141-2154.

Collin, S., Lima, A., Heringer, S., Sanders, V., Pessotti, H., & Deps, P. (2023). Systematic Review of Hansen Disease Attributed to Mycobacterium lepromatosis. Emerging Infectious Diseases, 29, 1376 - 1385. https://doi.org/10.3201/eid2907.230024.

Dharmawan, Y., Korfage, I., Abqari, U., Widjanarko, B., & Richardus, J. (2023). Measuring leprosy case detection delay and associated factors in Indonesia: a community-based study. BMC Infectious Diseases, 23. https://doi.org/10.1186/s12879-023-08552-x.

Pescarini, J., Strina, A., Nery, J., Skalinski, L., De Andrade, K., Penna, M., Brickley, E., Rodrigues, L., Barreto, M., & Penna, G. (2018). Socioeconomic risk markers of leprosy in high-burden countries: A systematic review and meta-analysis. PLoS Neglected Tropical Diseases, 12. https://doi.org/10.1371/journal.pntd.0006622.

Rahmawati, D., & Bimanto, H. (2021). Analysis of Factors Affecting Leprosy Cases in East Java Province with Spatial Autoregressive Model (SAR). Jurnal Berkala Kesehatan, 7, 59-64. https://doi.org/10.20527/jbk.v7i1.10285.

Radra, H., Yuziani, Y., & Topik, M. (2025). Evaluasi Penerapan Terapi MDT pada Kasus Kusta di Kota Lhokseumawe. JURNAL RISET RUMPUN ILMU KESEHATAN. https://doi.org/10.55606/jurrikes.v4i1.4542.

Badan Pusat Statistik Kota Surakarta. (2024). Jumlah penduduk menurut kabupaten/kota di Provinsi Jawa Tengah, 2021-2023. https://surakartakota.bps.go.id/id/statistics-table/2/NDAyIzI=/jumlah-pendudukmenurut-kabupaten-kota-di-provinsi-jawa-tengah.html

Badan Pusat Statistik Provinsi Jawa Tengah. (2024). Indikator strategis IPM: Rata-rata lama sekolah menurut kabupaten/kota di Provinsi Jawa Tengah, 2021-2023. https://jateng.bps.go.id/id/statistics-table/2/MjQxNCMy/-indikator-strategis-ipm-rata-rata-lama-sekolah-menurut-kabupaten-kota-di-provinsi-jawa-tengah.html

Badan Pusat Statistik Provinsi Jawa Tengah. (2024). Jumlah rumah sakit umum, rumah sakit khusus, puskesmas, klinik pratama, dan posyandu menurut kabupaten/kota di Provinsi Jawa Tengah, 2023. https://jateng.bps.go.id/id/statistics-table/3/YmlzemNGUkNVblZLVVhOblREWnZXbkEzWld0eVVUMDkjMw==/number-of-general-hospital–specialized-hospital–public-health-center–primary-clinic–and-integrated-health-post-by-regency-municipality-in-jawa-tengah-province–2018.html

Badan Pusat Statistik Provinsi Jawa Tengah. (2024). Kemiskinan (P0) menurut kabupaten/kota di Provinsi Jawa Tengah (Persen), 2022-2024. https://jateng.bps.go.id/id/statistics-table/2/MzQjMg==/kemiskinan-menurut-kabupaten-kota-di-provinsi-jawa-tengah.html

Badan Pusat Statistik Provinsi Jawa Tengah. (2024). Persentase rumah tangga yang memiliki akses terhadap air minum layak menurut kabupaten/kota di Provinsi Jawa Tengah, 2021-2023. https://jateng.bps.go.id/id/statistics-table/2/MjQ4MCMy/persentase-rumah-tangga-yang-memiliki-akses-terhadap-air-minum-layak-menurut-kabupaten-kota-di-provinsi-jawa-tengah.html

Badan Pusat Statistik Provinsi Jawa Tengah. (2024). Persentase rumah tangga yang memiliki akses terhadap sanitasi layak menurut kabupaten/kota di Provinsi Jawa Tengah, 2021-2023. https://jateng.bps.go.id/id/statistics-table/2/MjQ3OSMy/persentase-rumah-tangga-yang-memiliki-akses-terhadap-sanitasi-layak-menurut-kabupaten-kota-di-provinsi-jawa-tengah.html

Badan Pusat Statistik Provinsi Jawa Tengah. (2024). Persentase rumah tangga yang menempati rumah layak huni di Provinsi Jawa Tengah, 2021-2023. https://jateng.bps.go.id/id/statistics-table/2/MjQyOCMy/persentase-rumahtangga-yang-menempati-rumah-layak-huni-di-provinsi-jawa-tengah.html

Badan Pusat Statistik Provinsi Jawa Tengah. (2025). Jumlah penduduk, laju pertumbuhan penduduk, distribusi persentase penduduk, kepadatan penduduk, rasio jenis kelamin penduduk menurut kabupaten/kota di Provinsi Jawa Tengah, 2024. https://jateng.bps.go.id/id/statistics-table/3/V1ZSbFRUY3lTbFpEYTNsVWNGcDZjek53YkhsNFFUMDkjMw==/jumlah-penduduk–laju-pertumbuhan-penduduk–distribusi-persentase-penduduk–kepadatan-penduduk–rasio-jenis-kelamin-penduduk-menurut-kabupaten-kota-di-provinsi-jawa-tengah–2025.html?year=2024

Satu Data Provinsi Jawa Tengah. (2024). Jumlah kasus terdaftar dan angka prevalensi penyakit kusta tahun 2024. Dinas Kesehatan Provinsi Jawa Tengah. https://data.jatengprov.go.id/dataset/jumlah-kasus-terdaftar-dan-angka-prevelensi-penyakit-kusta-tahun-2024

LAMPIRAN

Link Video Presentasi
https://bit.ly/PresentasiKustaJawaTengahEpidemiologi2025

Link Dashboard
https://farahnadira.shinyapps.io/KustaJawaTengahEpidemiologi2025/

Syntax Dashboard

# 01_data_loading.R - Data Loading and Helper Functions
suppressPackageStartupMessages({
  library(shiny)
  library(shinydashboard)
  library(DT)
  library(dplyr)
  library(ggplot2)
  library(readr)
  library(readxl)
  library(stringr)
  library(scales)
  library(bslib)
  library(tidyr)
  library(leaflet)
  library(plotly)
  library(MASS)
  library(car)
})

# SETTING Data
DATA_FILE <- "Dataset Kusta Jawa Tengah 2024.xlsx"
DATA_TREN_FILE <- "Data Tren Kasus Baru Kusta 2021-2024.xlsx"
safe_num <- function(x) {
  if (is.numeric(x)) return(as.numeric(x))
  x2 <- as.character(x)
  x2 <- stringr::str_replace_all(x2, "\\s+", "")
  x2 <- stringr::str_replace_all(x2, "\\.", "")
  x2 <- stringr::str_replace_all(x2, ",", ".")
  x2 <- stringr::str_replace_all(x2, "[^0-9eE+\\-\\.]", "")
  suppressWarnings(as.numeric(x2))
}

# Load data
initial_data <- tryCatch({
  if (!file.exists(DATA_FILE)) {
    stop("File data utama tidak ditemukan")
  }
  
  df <- readxl::read_excel(DATA_FILE)
  
  names(df) <- stringr::str_trim(names(df))
  df <- df[, !grepl("^Unnamed", names(df))]
  
  numeric_cols <- c(
    "Jumlah Kasus Baru", "Jumlah Kasus Tercatat",
    "Jumlah Kasus Baru Cacat Tingkat 0", "Jumlah Baru Kasus Cacat Tingkat 2",
    "Jumlah Pasien selesai berobat di tahun 2024",
    "Persentase Pasien Kusta Basah menyelesaikan pengobatan tepat waktu",
    "Jumlah Penduduk",
    "Jumlah Kasus Kusta Kering Anak", "Jumlah Kasus Kusta Kering Dewasa",
    "Jumlah Kasus Kusta Kering", "Jumlah Kasus Kusta Basah Anak",
    "Jumlah Kasus Kusta Basah Dewasa", "Jumlah Kasus Kusta Basah",
    "Jumlah Kasus Kusta Anak", "Jumlah Kasus Kusta Dewasa",
    "Kepadatan Penduduk (per km^2)", "Persentase Penduduk Miskin",
    "Persentase Rumah Layak Huni", "Persentase Akses Sanitasi Layak",
    "Persentase Akses Air Minum Layak", "Jumlah Fasilitas Kesehatan",
    "Rata-Rata Lama Sekolah (IPM) (Tahun)"
  )
  
  for (col in numeric_cols) {
    if (col %in% names(df)) {
      df[[col]] <- safe_num(df[[col]])
    }
  }
  
  df <- df %>%
    dplyr::mutate(
      Insidensi = round((`Jumlah Kasus Baru` / `Jumlah Penduduk`) * 100000, 2),
      Prevalensi = round((`Jumlah Kasus Tercatat` / `Jumlah Penduduk`) * 100000, 2),
      Prevalensi_Kusta_Basah = round((`Jumlah Kasus Kusta Basah` / `Jumlah Penduduk`) * 100000, 2),
      Prevalensi_Kusta_Kering = round((`Jumlah Kasus Kusta Kering` / `Jumlah Penduduk`) * 100000, 2),
      Persen_Basah = round((`Jumlah Kasus Kusta Basah` / `Jumlah Kasus Tercatat`) * 100, 2)
    )
  
  df
}, error = function(e) {
  warning(e$message)
  NULL
})

data_tren <- tryCatch({
  if (!file.exists(DATA_TREN_FILE)) {
    stop("File data tren tidak ditemukan")
  }
  
  df <- readxl::read_excel(DATA_TREN_FILE)
  
  names(df)[1] <- "Wilayah"
  df$Wilayah <- stringr::str_trim(df$Wilayah)
  
  year_cols <- names(df)[names(df) != "Wilayah"]
  
  for (col in year_cols) {
    df[[col]] <- safe_num(df[[col]])
  }
  
  message("Tren data loaded successfully: ", nrow(df), " rows")
  df
}, error = function(e) {
  warning(e$message)
  NULL
})

# Koordinat Jawa Tengah
coords_jateng <- data.frame(
  Wilayah = c("Kab. Cilacap", "Kab. Banyumas", "Kab. Purbalingga", "Kab. Banjarnegara", 
              "Kab. Kebumen", "Kab. Purworejo", "Kab. Wonosobo", "Kab. Magelang",
              "Kab. Boyolali", "Kab. Klaten", "Kab. Sukoharjo", "Kab. Wonogiri",
              "Kab. Karanganyar", "Kab. Sragen", "Kab. Grobogan", "Kab. Blora",
              "Kab. Rembang", "Kab. Pati", "Kab. Kudus", "Kab. Jepara",
              "Kab. Demak", "Kab. Semarang", "Kab. Temanggung", "Kab. Kendal",
              "Kab. Batang", "Kab. Pekalongan", "Kab. Pemalang", "Kab. Tegal",
              "Kab. Brebes", "Kota Magelang", "Kota Surakarta", "Kota Salatiga",
              "Kota Semarang", "Kota Pekalongan", "Kota Tegal"),
  lat = c(-7.7261, -7.5200, -7.3881, -7.3858, -7.6671, -7.6700, -7.3661, -7.4700,
          -7.5403, -7.7117, -7.6800, -7.8144, -7.5971, -7.4261, -7.0300, -6.9697,
          -6.7089, -6.7500, -6.8047, -6.5928, -6.8906, -7.1950, -7.3200, -6.9167,
          -6.9056, -6.8889, -6.8983, -6.8694, -6.8700, -7.4797, -7.5755, -7.3319,
          -6.9933, -6.8886, -6.8694),
  lng = c(109.0128, 109.2400, 109.3581, 109.6872, 109.6544, 110.0200, 109.9028, 110.2172,
          110.5961, 110.6061, 110.8200, 110.9261, 110.9775, 111.0200, 110.8000, 111.4189,
          111.3425, 111.0300, 110.8406, 110.6686, 110.6400, 110.4200, 110.1700, 110.2028,
          110.0350, 109.6739, 109.3758, 109.1403, 108.8400, 110.2175, 110.8244, 110.4981,
          110.4203, 109.6758, 109.1403)
)

# ----------------------------
# 02_ui.R - User Interface Definition

theme_bs <- bs_theme(
  bootswatch = "minty",
  bg = "#FAFAFA",
  fg = "#0f3d3e",
  primary = "#2C7873"
)

ui <- dashboardPage(
  skin = "black",
  dashboardHeader(title = span("Dashboard Analisis Kasus Kusta Provinsi Jawa Tengah", style = "font-weight:700;")),
  
  dashboardSidebar(
    width = 300,
    sidebarMenu(
      id = "tabs",
      menuItem("Data", tabName = "data", icon = icon("table")),
      menuItem("Statistika Deskriptif", tabName = "statdes", icon = icon("chart-bar")),
      menuItem("Ukuran Frekuensi", tabName = "freq", icon = icon("percent")),
      menuItem("Interpretasi Epidemiologi", tabName = "interp", icon = icon("stethoscope")),
      menuItem("Peta Visualisasi", tabName = "peta", icon = icon("globe")),
      menuItem("Tren Waktu", tabName = "tren", icon = icon("line-chart")),
      menuItem("Pemodelan", tabName = "pemodelan", icon = icon("calculator"))
    )
  ),
  
  dashboardBody(
    theme = theme_bs,
    tags$head(tags$style(HTML("
      body, .content-wrapper, .right-side { background-color: #FAFAFA; }
      .box { border-radius: 10px; border: 1px solid #ECECEC; box-shadow: 0 2px 6px rgba(0,0,0,0.04); }
      .box .box-title { color: #115e59; font-weight: 700; }
      .value-box { border-radius: 8px; }
      .small-box { border-radius: 8px; }
      .interpretation-box {
        background: #e8f5e9;
        border-left: 4px solid #4caf50;
        padding: 15px;
        margin: 10px 0;
        border-radius: 5px;
      }
    "))),
    
    tabItems(
      # TAB DATA
      tabItem(tabName = "data",
              fluidRow(
                valueBoxOutput("vb_total_wilayah", width = 2),
                valueBoxOutput("vb_total_kasus_baru", width = 2),
                valueBoxOutput("vb_total_kasus_tercatat", width = 2),
                valueBoxOutput("vb_total_kasus_basah", width = 2),
                valueBoxOutput("vb_total_kasus_kering", width = 2)
              ),
              fluidRow(
                box(width = 12,
                    tags$div(class = "interpretation-box",
                             HTML("<b>Interpretasi:</b> 96.9% kasus tercatat merupakan kasus kusta basah, yang menunjukkan bahwa pencegahan penyebaran penyakit masih belum optimal. Kusta basah (MB/Multibacillary) lebih menular dibanding kusta kering, sehingga memerlukan perhatian khusus dalam upaya pengendalian.")
                    )
                )
              ),
              fluidRow(
                box(title = "Data - Indikator Kusta", status = "info", solidHeader = TRUE, width = 12,
                    DTOutput("table_kusta")
                )
              ),
              fluidRow(
                box(title = "Data - Indikator Sosial Ekonomi dan Kesehatan", status = "warning", solidHeader = TRUE, width = 12,
                    DTOutput("table_sosek")
                )
              )
      ),
      
      # TAB STATISTIKA DESKRIPTIF
      tabItem(tabName = "statdes",
              fluidRow(
                box(title = "Pilih Variabel", status = "primary", solidHeader = TRUE, width = 8,
                    uiOutput("select_var_ui")
                ),
                box(title = "Filter Tahun", status = "info", solidHeader = TRUE, width = 4,
                    uiOutput("year_filter_statdes")
                )
              ),
              fluidRow(
                box(title = "Histogram", status = "success", solidHeader = TRUE, width = 6,
                    plotOutput("histogram_plot", height = "400px")
                ),
                box(title = "Boxplot", status = "info", solidHeader = TRUE, width = 6,
                    plotlyOutput("boxplot_plot", height = "400px")
                )
              ),
              fluidRow(
                box(title = "Info Data Terpilih", status = "warning", solidHeader = TRUE, width = 6,
                    htmlOutput("clicked_info")
                ),
                box(title = "Distribusi Kategori Kasus", status = "primary", solidHeader = TRUE, width = 6,
                    selectInput("pie_var", "Pilih Kategori:",
                                choices = c("Kusta Kering: Anak vs Dewasa" = "kering",
                                            "Kusta Basah: Anak vs Dewasa" = "basah",
                                            "Kusta Kering vs Basah" = "tipe"),
                                selected = "tipe"),
                    plotlyOutput("pie_chart", height = "300px")
                )
              ),
              fluidRow(
                box(title = "Statistika Deskriptif", status = "primary", solidHeader = TRUE, width = 12,
                    DTOutput("statdesc_tbl")
                )
              )
      ),
      
      # TAB UKURAN FREKUENSI
      tabItem(tabName = "freq",
              fluidRow(
                box(title = "Ringkasan Ukuran Frekuensi (agregat)", status = "primary", solidHeader = TRUE, width = 12,
                    htmlOutput("freq_summary")
                )
              ),
              fluidRow(
                box(title = "Insidensi (per 100,000 penduduk)", 
                    status = "success", solidHeader = TRUE, width = 6,
                    plotlyOutput("incidence_plot", height = "600px")
                ),
                box(title = "Prevalensi (per 100,000 penduduk)", 
                    status = "info", solidHeader = TRUE, width = 6,
                    plotlyOutput("prevalensi_plot", height = "600px")
                )
              ),
              fluidRow(
                box(title = "Prevalensi Kusta Basah (per 100,000 penduduk)", 
                    status = "warning", solidHeader = TRUE, width = 6,
                    plotlyOutput("prevalensi_basah_plot", height = "600px")
                ),
                box(title = "Prevalensi Kusta Kering (per 100,000 penduduk)", 
                    status = "primary", solidHeader = TRUE, width = 6,
                    plotlyOutput("prevalensi_kering_plot", height = "600px")
                )
              ),
              fluidRow(
                box(title = "Tabel Ukuran Frekuensi per Wilayah", status = "info", solidHeader = TRUE, width = 12,
                    DTOutput("freq_table")
                )
              )
      ),
      
      # TAB INTERPRETASI EPIDEMIOLOGI
      tabItem(tabName = "interp",
              fluidRow(
                box(title = "Ringkasan Ukuran Epidemiologi", status = "primary", solidHeader = TRUE, width = 12,
                    htmlOutput("epi_summary")
                )
              ),
              fluidRow(
                box(title = "Pilih Ukuran untuk Visualisasi", status = "primary", solidHeader = TRUE, width = 12,
                    uiOutput("epi_var_plot_ui")
                )
              ),
              fluidRow(
                box(title = "Grafik Perbandingan Top 15 Wilayah", status = "success", solidHeader = TRUE, width = 12,
                    plotOutput("epi_comparison_plot", height = "600px")
                )
              ),
              fluidRow(
                box(title = "Interpretasi Epidemiologi", status = "info", solidHeader = TRUE, width = 12,
                    htmlOutput("epi_interpretation")
                )
              )
      ),
      
      # TAB PETA VISUALISASI
      tabItem(tabName = "peta",
              fluidRow(
                box(title = "Kontrol Peta", status = "primary", solidHeader = TRUE, width = 12,
                    fluidRow(
                      column(4, uiOutput("map_var_ui")),
                      column(4, uiOutput("map_year_ui")),
                      column(4, 
                             selectInput("color_low", "Warna Rendah:", 
                                         choices = c("Kuning" = "yellow", "Putih" = "white", "Biru Muda" = "lightblue"),
                                         selected = "yellow"),
                             selectInput("color_high", "Warna Tinggi:",
                                         choices = c("Merah" = "red", "Merah Tua" = "darkred", "Oranye" = "orange"),
                                         selected = "red"))
                    )
                )
              ),
              fluidRow(
                box(title = "Peta Visualisasi Jawa Tengah", status = "info", solidHeader = TRUE, width = 9,
                    leafletOutput("map_leaflet", height = "600px"),
                    tags$div(style = "margin-top: 15px; padding: 10px; background: #f0f0f0; border-radius: 5px;",
                             HTML("<b>Keterangan:</b> Warna yang lebih tua dan ukuran circle yang lebih besar menunjukkan nilai yang lebih tinggi.")
                    )
                ),
                box(title = "Informasi Peta", status = "warning", solidHeader = TRUE, width = 3,
                    htmlOutput("map_info")
                )
              ),
              fluidRow(
                box(title = "Interpretasi Pola Spasial", status = "success", solidHeader = TRUE, width = 12,
                    htmlOutput("spatial_interpretation")
                )
              )
      ),
      
      # TAB TREN WAKTU
      tabItem(tabName = "tren",
              fluidRow(
                box(title = "Grafik Tren Kasus Baru Kusta Jawa Tengah 2021-2024", 
                    status = "primary", solidHeader = TRUE, width = 12,
                    plotlyOutput("tren_plot", height = "500px")
                )
              ),
              fluidRow(
                box(title = "Pilih Wilayah untuk Perbandingan", 
                    status = "info", solidHeader = TRUE, width = 12,
                    uiOutput("tren_wilayah_ui")
                )
              ),
              fluidRow(
                box(title = "Tren Kasus Baru Kusta per Wilayah", 
                    status = "success", solidHeader = TRUE, width = 12,
                    plotlyOutput("tren_wilayah_plot", height = "500px")
                )
              ),
              fluidRow(
                box(title = "Tabel Data Tren", status = "warning", solidHeader = TRUE, width = 12,
                    DTOutput("tren_table")
                )
              ),
              fluidRow(
                box(title = "Interpretasi Tren", status = "primary", solidHeader = TRUE, width = 12,
                    htmlOutput("tren_interpretation")
                )
              )
      ),
      
      # TAB MODELLING
      tabItem(tabName = "pemodelan",
              fluidRow(
                valueBoxOutput("vb_jumlah_observasi", width = 4),
                valueBoxOutput("vb_variabel_x", width = 4),
                valueBoxOutput("vb_variabel_y", width = 4)
              ),
              fluidRow(
                box(title = "Informasi Dataset", status = "primary", solidHeader = TRUE, width = 12,
                    htmlOutput("model_info")
                )
              ),
              fluidRow(
                box(title = "Uji Overdispersi", status = "warning", solidHeader = TRUE, width = 6,
                    htmlOutput("uji_asumsi")
                ),
                box(title = "Uji Multikolinearitas", status = "info", solidHeader = TRUE, width = 6,
                    htmlOutput("uji_multikol")
                )
              ),
              fluidRow(
                box(title = "Model Full", status = "primary", solidHeader = TRUE, width = 4,
                    htmlOutput("model_full_summary"),
                    tags$div(style = "margin-top: 15px; padding: 10px; background: #fff3cd; border-radius: 5px;",
                             htmlOutput("model_full_note")
                    )
                ),
                box(title = "Model Signifikan", status = "success", solidHeader = TRUE, width = 4,
                    htmlOutput("model_sign_summary"),
                    tags$div(style = "margin-top: 15px; padding: 10px; background: #d4edda; border-radius: 5px;",
                             HTML("<b>Catatan:</b> Merupakan model dengan hanya menggunakan variabel yang signifikan dari model full.")
                    )
                ),
                box(title = "Model StepAIC", status = "info", solidHeader = TRUE, width = 4,
                    htmlOutput("model_step_summary"),
                    tags$div(style = "margin-top: 15px; padding: 10px; background: #d1ecf1; border-radius: 5px;",
                             HTML("<b>Catatan:</b> Merupakan model dengan kombinasi variabel dengan AIC terkecil dari beberapa kombinasi variabel yang dicoba dengan menggunakan metode StepAIC. Variabel Rumah dengan Akses Sanitasi Layak (%) (yang dianggap tidak signifikan pada model full) turut diikutsertakan karena menghasilkan nilai AIC terkecil dari seluruh kombinasi StepAIC yang dicoba.")
                    )
                )
              ),
              fluidRow(
                box(title = "Perbandingan Model", status = "danger", solidHeader = TRUE, width = 12,
                    DTOutput("model_comparison"),
                    tags$div(style = "margin-top: 15px; padding: 10px; background: #f8d7da; border-radius: 5px;",
                             htmlOutput("model_comparison_interp")
                    )
                )
              ),
              fluidRow(
                box(title = "Bentuk Model Terpilih", status = "success", solidHeader = TRUE, width = 12,
                    htmlOutput("model_final_equation"),
                    htmlOutput("model_interpretation")
                )
              )
      )
    )
  )
)

# ----------------------------
# 03a_server_part1.R - Server Logic (Part 1)

server_part1 <- function(input, output, session, std_df, tren_df) {
  
  # ============ TAB DATA ============
  output$vb_total_wilayah <- renderValueBox({
    df <- std_df()
    val <- if(is.null(df)) 0 else nrow(df)
    valueBox(val, "Jumlah Kabupaten/Kota", icon = icon("map-marker"), color = "aqua")
  })
  
  output$vb_total_kasus_baru <- renderValueBox({
    df <- std_df()
    val <- if(is.null(df)) 0 else comma(sum(df$`Jumlah Kasus Baru`, na.rm = TRUE))
    valueBox(val, "Jumlah Kasus Baru", icon = icon("user-plus"), color = "green")
  })
  
  output$vb_total_kasus_tercatat <- renderValueBox({
    df <- std_df()
    val <- if(is.null(df)) 0 else comma(sum(df$`Jumlah Kasus Tercatat`, na.rm = TRUE))
    valueBox(val, "Jumlah Kasus Tercatat", icon = icon("clipboard-list"), color = "orange")
  })
  
  output$vb_total_kasus_basah <- renderValueBox({
    df <- std_df()
    val <- if(is.null(df)) 0 else comma(sum(df$`Jumlah Kasus Kusta Basah`, na.rm = TRUE))
    valueBox(val, "Jumlah Kasus Kusta Basah", icon = icon("droplet"), color = "red")
  })
  
  output$vb_total_kasus_kering <- renderValueBox({
    df <- std_df()
    val <- if(is.null(df)) 0 else comma(sum(df$`Jumlah Kasus Kusta Kering`, na.rm = TRUE))
    valueBox(val, "Jumlah Kasus Kusta Kering", icon = icon("sun"), color = "yellow")
  })
  
  output$table_kusta <- renderDT({
    df <- std_df()
    if (is.null(df)) return(NULL)
    
    df_display <- df[, c("Wilayah", "Jumlah Kasus Baru", "Jumlah Kasus Baru Cacat Tingkat 0",
                         "Jumlah Baru Kasus Cacat Tingkat 2", "Jumlah Pasien selesai berobat di tahun 2024",
                         "Jumlah Penduduk", "Persentase Pasien Kusta Basah menyelesaikan pengobatan tepat waktu",
                         "Jumlah Kasus Kusta Kering Anak", "Jumlah Kasus Kusta Kering Dewasa",
                         "Jumlah Kasus Kusta Kering", "Jumlah Kasus Kusta Basah Anak",
                         "Jumlah Kasus Kusta Basah Dewasa", "Jumlah Kasus Kusta Basah",
                         "Jumlah Kasus Kusta Anak", "Jumlah Kasus Kusta Dewasa", "Jumlah Kasus Tercatat")]
    
    datatable(df_display, options = list(pageLength = 10, scrollX = TRUE), rownames = FALSE) %>%
      formatRound(columns = "Jumlah Penduduk", digits = 0)
  })
  
  output$table_sosek <- renderDT({
    df <- std_df()
    if (is.null(df)) return(NULL)
    
    df_display <- df[, c("Wilayah", "Kepadatan Penduduk (per km^2)", "Persentase Penduduk Miskin",
                         "Persentase Rumah Layak Huni", "Persentase Akses Sanitasi Layak",
                         "Persentase Akses Air Minum Layak", "Jumlah Fasilitas Kesehatan",
                         "Rata-Rata Lama Sekolah (IPM) (Tahun)")]
    
    datatable(df_display, options = list(pageLength = 10, scrollX = TRUE), rownames = FALSE) %>%
      formatRound(columns = c("Kepadatan Penduduk (per km^2)", "Persentase Penduduk Miskin",
                              "Persentase Rumah Layak Huni", "Persentase Akses Sanitasi Layak",
                              "Persentase Akses Air Minum Layak", "Rata-Rata Lama Sekolah (IPM) (Tahun)"), 
                  digits = 2)
  })
  
  # ============ TAB STATISTIKA DESKRIPTIF ============
  output$select_var_ui <- renderUI({
    df <- std_df()
    if (is.null(df)) return(NULL)
    
    num_vars <- names(df)[sapply(df, is.numeric)]
    num_vars <- num_vars[!num_vars %in% c("Wilayah")]
    
    selectInput("select_var", "Pilih Variabel:", choices = num_vars, selected = "Jumlah Kasus Baru")
  })
  
  output$year_filter_statdes <- renderUI({
    df <- std_df()
    df_tren <- tren_df()
    req(input$select_var)
    
    if(input$select_var == "Jumlah Kasus Baru" && !is.null(df_tren)) {
      year_cols <- names(df_tren)[names(df_tren) != "Wilayah"]
      selectInput("year_statdes", "Pilih Tahun:", 
                  choices = c("2024 (Data Utama)", year_cols), 
                  selected = "2024 (Data Utama)")
    } else {
      tags$p("Tahun: 2024", style = "margin-top: 25px; font-weight: bold;")
    }
  })
  
  output$histogram_plot <- renderPlot({
    df <- std_df()
    df_tren <- tren_df()
    req(df, input$select_var)
    
    var_name <- input$select_var
    
    if(var_name == "Jumlah Kasus Baru" && !is.null(df_tren) && !is.null(input$year_statdes)) {
      if(input$year_statdes != "2024 (Data Utama)") {
        data_plot <- df_tren[[input$year_statdes]]
        data_plot <- data_plot[!is.na(data_plot)]
        var_label <- paste("Jumlah Kasus Baru", input$year_statdes)
      } else {
        data_plot <- df[[var_name]][!is.na(df[[var_name]])]
        var_label <- var_name
      }
    } else {
      data_plot <- df[[var_name]][!is.na(df[[var_name]])]
      var_label <- var_name
    }
    
    ggplot(data.frame(x = data_plot), aes(x = x)) +
      geom_histogram(bins = 15, fill = "#2C7873", color = "white", alpha = 0.8) +
      theme_minimal(base_size = 14) +
      labs(title = paste("Histogram:", var_label), x = var_label, y = "Frekuensi") +
      theme(plot.title = element_text(face = "bold", size = 16), panel.grid.minor = element_blank())
  }, res = 96)
  
  clicked_point <- reactiveVal(NULL)
  
  output$boxplot_plot <- renderPlotly({
    df <- std_df()
    df_tren <- tren_df()
    req(df, input$select_var)
    
    var_name <- input$select_var
    
    if(var_name == "Jumlah Kasus Baru" && !is.null(df_tren) && !is.null(input$year_statdes)) {
      if(input$year_statdes != "2024 (Data Utama)") {
        data_plot <- data.frame(
          Wilayah = df_tren$Wilayah,
          value = df_tren[[input$year_statdes]]
        ) %>% filter(!is.na(value))
        var_label <- paste("Jumlah Kasus Baru", input$year_statdes)
      } else {
        data_plot <- df %>% filter(!is.na(!!sym(var_name)))
        data_plot <- data.frame(Wilayah = data_plot$Wilayah, value = data_plot[[var_name]])
        var_label <- var_name
      }
    } else {
      data_plot <- df %>% filter(!is.na(!!sym(var_name)))
      data_plot <- data.frame(Wilayah = data_plot$Wilayah, value = data_plot[[var_name]])
      var_label <- var_name
    }
    
    plot_ly(source = "boxplot_source") %>%
      add_trace(
        data = data_plot, y = ~value, x = 0, type = "box", name = "",
        marker = list(color = "#2C7873"), customdata = ~Wilayah,
        hovertemplate = paste('<b>%{customdata}</b><br>Nilai: %{y}<br><extra></extra>')
      ) %>%
      add_markers(
        data = data_plot, y = ~value, x = 0,
        marker = list(size = 8, color = "#2C7873", opacity = 0),
        customdata = ~Wilayah,
        hovertemplate = paste('<b>%{customdata}</b><br>Nilai: %{y}<br><extra></extra>'),
        showlegend = FALSE
      ) %>%
      layout(
        title = paste("Boxplot:", var_label),
        yaxis = list(title = var_label),
        xaxis = list(title = "", showticklabels = FALSE),
        showlegend = FALSE
      )
  })
  
  observeEvent(event_data("plotly_click", source = "boxplot_source"), {
    df <- std_df()
    df_tren <- tren_df()
    req(df, input$select_var)
    
    click <- event_data("plotly_click", source = "boxplot_source")
    if(!is.null(click)) {
      clicked_y <- click$y
      var_name <- input$select_var
      
      if(var_name == "Jumlah Kasus Baru" && !is.null(df_tren) && !is.null(input$year_statdes)) {
        if(input$year_statdes != "2024 (Data Utama)") {
          matching_row <- df_tren %>%
            filter(abs(!!sym(input$year_statdes) - clicked_y) < 0.5) %>%
            slice(1)
          matching_row$nilai_terpilih <- matching_row[[input$year_statdes]]
        } else {
          matching_row <- df %>%
            filter(!is.na(!!sym(var_name))) %>%
            filter(abs(!!sym(var_name) - clicked_y) < 0.5) %>%
            slice(1)
          matching_row$nilai_terpilih <- matching_row[[var_name]]
        }
      } else {
        matching_row <- df %>%
          filter(!is.na(!!sym(var_name))) %>%
          filter(abs(!!sym(var_name) - clicked_y) < 0.5) %>%
          slice(1)
        matching_row$nilai_terpilih <- matching_row[[var_name]]
      }
      
      if(nrow(matching_row) > 0) {
        clicked_point(matching_row)
      }
    }
  })
  
  output$clicked_info <- renderUI({
    df <- std_df()
    req(df, input$select_var)
    
    var_name <- input$select_var
    data_valid <- df %>% filter(!is.na(!!sym(var_name)))
    clicked <- clicked_point()
    
    if(!is.null(clicked)) {
      HTML(paste0(
        "<p><b>Variabel:</b> ", var_name, "</p>",
        "<p><b>Wilayah Terpilih:</b> ", clicked$Wilayah, "</p>",
        "<p><b>Nilai:</b> ", round(clicked$nilai_terpilih, 2), "</p>",
        "<hr>",
        "<p><b>Total Data Valid:</b> ", nrow(data_valid), "</p>"
      ))
    } else {
      HTML(paste0(
        "<p><b>Variabel:</b> ", var_name, "</p>",
        "<p><b>Jumlah Data Valid:</b> ", nrow(data_valid), "</p>",
        "<p><i>Klik pada titik boxplot untuk melihat detail wilayah</i></p>"
      ))
    }
  })
  
  output$pie_chart <- renderPlotly({
    df <- std_df()
    req(df, input$pie_var)
    
    total_df <- df %>%
      summarise(
        Kering_Anak = sum(`Jumlah Kasus Kusta Kering Anak`, na.rm = TRUE),
        Kering_Dewasa = sum(`Jumlah Kasus Kusta Kering Dewasa`, na.rm = TRUE),
        Basah_Anak = sum(`Jumlah Kasus Kusta Basah Anak`, na.rm = TRUE),
        Basah_Dewasa = sum(`Jumlah Kasus Kusta Basah Dewasa`, na.rm = TRUE),
        Total_Kering = sum(`Jumlah Kasus Kusta Kering`, na.rm = TRUE),
        Total_Basah = sum(`Jumlah Kasus Kusta Basah`, na.rm = TRUE)
      )
    
    if(input$pie_var == "kering") {
      pie_data <- data.frame(
        Kategori = c("Anak", "Dewasa"),
        Jumlah = c(total_df$Kering_Anak, total_df$Kering_Dewasa)
      )
      title_text <- "Distribusi Kusta Kering"
    } else if(input$pie_var == "basah") {
      pie_data <- data.frame(
        Kategori = c("Anak", "Dewasa"),
        Jumlah = c(total_df$Basah_Anak, total_df$Basah_Dewasa)
      )
      title_text <- "Distribusi Kusta Basah"
    } else {
      pie_data <- data.frame(
        Kategori = c("Kering", "Basah"),
        Jumlah = c(total_df$Total_Kering, total_df$Total_Basah)
      )
      title_text <- "Distribusi Tipe Kusta"
    }
    
    pie_data$Persen <- round(pie_data$Jumlah / sum(pie_data$Jumlah) * 100, 1)
    
    plot_ly(pie_data, labels = ~Kategori, values = ~Jumlah, type = 'pie',
            textposition = 'inside', textinfo = 'label+percent', hoverinfo = 'text',
            text = ~paste(Kategori, '<br>', Jumlah, 'kasus<br>', Persen, '%'),
            marker = list(colors = c('#66c2a5', '#fc8d62'), line = list(color = '#FFFFFF', width = 2))) %>%
      layout(title = title_text, showlegend = TRUE)
  })
  
  output$statdesc_tbl <- renderDT({
    df <- std_df()
    if (is.null(df)) return(NULL)
    
    num_vars <- names(df)[sapply(df, is.numeric)]
    num_vars <- num_vars[!num_vars %in% c("Wilayah")]
    
    stat_df <- data.frame(
      Variabel = character(), N = integer(), Mean = numeric(),
      Median = numeric(), SD = numeric(), Min = numeric(), Max = numeric(),
      stringsAsFactors = FALSE
    )
    
    for(var in num_vars) {
      vals <- df[[var]][!is.na(df[[var]])]
      if(length(vals) > 0) {
        stat_df <- rbind(stat_df, data.frame(
          Variabel = var, N = length(vals), Mean = mean(vals, na.rm = TRUE),
          Median = median(vals, na.rm = TRUE), SD = sd(vals, na.rm = TRUE),
          Min = min(vals, na.rm = TRUE), Max = max(vals, na.rm = TRUE)
        ))
      }
    }
    
    datatable(stat_df, options = list(pageLength = 15, scrollX = TRUE), rownames = FALSE) %>%
      formatRound(columns = c("Mean", "Median", "SD", "Min", "Max"), digits = 2)
  })
  
  # ============ TAB UKURAN FREKUENSI ============
  output$freq_summary <- renderUI({
    df <- std_df()
    if (is.null(df)) return(HTML("<i>Data tidak tersedia</i>"))
    
    total_kasus_baru <- sum(df$`Jumlah Kasus Baru`, na.rm = TRUE)
    total_kasus_tercatat <- sum(df$`Jumlah Kasus Tercatat`, na.rm = TRUE)
    total_penduduk <- sum(df$`Jumlah Penduduk`, na.rm = TRUE)
    
    incidence_agregat <- round((total_kasus_baru / total_penduduk) * 100000, 2)
    prevalensi_agregat <- round((total_kasus_tercatat / total_penduduk) * 100000, 2)
    
    HTML(paste0(
      "<div style='font-size:16px; line-height:2;'>",
      "<table style='width:100%; border-collapse: collapse;'>",
      "<tr style='background:#f0f0f0;'><th style='padding:10px; text-align:left; border:1px solid #ddd;'>Ukuran Frekuensi</th>",
      "<th style='padding:10px; text-align:left; border:1px solid #ddd;'>Nilai (per 100,000 penduduk)</th></tr>",
      "<tr><td style='padding:8px; border:1px solid #ddd;'><b>Insidensi Agregat</b></td>",
      "<td style='padding:8px; border:1px solid #ddd;'>", incidence_agregat, "</td></tr>",
      "<tr style='background:#f9f9f9;'><td style='padding:8px; border:1px solid #ddd;'><b>Prevalensi Agregat</b></td>",
      "<td style='padding:8px; border:1px solid #ddd;'>", prevalensi_agregat, "</td></tr>",
      "</table></div>"
    ))
  })
  
  output$incidence_plot <- renderPlotly({
    df <- std_df()
    req(df)
    plot_data <- df %>% filter(!is.na(Insidensi)) %>% arrange(desc(Insidensi))
    plot_ly(plot_data, x = ~Insidensi, y = ~reorder(Wilayah, Insidensi),
            type = "bar", orientation = "h", marker = list(color = "#2C7873"),
            text = ~paste("Insidensi:", round(Insidensi, 2)), hoverinfo = "text+y") %>%
      layout(xaxis = list(title = "Insidensi (per 100,000 penduduk)"),
             yaxis = list(title = ""), margin = list(l = 150))
  })
  
  output$prevalensi_plot <- renderPlotly({
    df <- std_df()
    req(df)
    plot_data <- df %>% filter(!is.na(Prevalensi)) %>% arrange(desc(Prevalensi))
    plot_ly(plot_data, x = ~Prevalensi, y = ~reorder(Wilayah, Prevalensi),
            type = "bar", orientation = "h", marker = list(color = "#5bc0de"),
            text = ~paste("Prevalensi:", round(Prevalensi, 2)), hoverinfo = "text+y") %>%
      layout(xaxis = list(title = "Prevalensi (per 100,000 penduduk)"),
             yaxis = list(title = ""), margin = list(l = 150))
  })
  
  output$prevalensi_basah_plot <- renderPlotly({
    df <- std_df()
    req(df)
    plot_data <- df %>% filter(!is.na(Prevalensi_Kusta_Basah)) %>% arrange(desc(Prevalensi_Kusta_Basah))
    plot_ly(plot_data, x = ~Prevalensi_Kusta_Basah, y = ~reorder(Wilayah, Prevalensi_Kusta_Basah),
            type = "bar", orientation = "h", marker = list(color = "#d9534f"),
            text = ~paste("Prevalensi Basah:", round(Prevalensi_Kusta_Basah, 2)), hoverinfo = "text+y") %>%
      layout(xaxis = list(title = "Prevalensi Kusta Basah (per 100,000 penduduk)"),
             yaxis = list(title = ""), margin = list(l = 150))
  })
  
  output$prevalensi_kering_plot <- renderPlotly({
    df <- std_df()
    req(df)
    plot_data <- df %>% filter(!is.na(Prevalensi_Kusta_Kering)) %>% arrange(desc(Prevalensi_Kusta_Kering))
    plot_ly(plot_data, x = ~Prevalensi_Kusta_Kering, y = ~reorder(Wilayah, Prevalensi_Kusta_Kering),
            type = "bar", orientation = "h", marker = list(color = "#f0ad4e"),
            text = ~paste("Prevalensi Kering:", round(Prevalensi_Kusta_Kering, 2)), hoverinfo = "text+y") %>%
      layout(xaxis = list(title = "Prevalensi Kusta Kering (per 100,000 penduduk)"),
             yaxis = list(title = ""), margin = list(l = 150))
  })
  
  output$freq_table <- renderDT({
    df <- std_df()
    if (is.null(df)) return(NULL)
    df_display <- df[, c("Wilayah", "Jumlah Kasus Baru", "Jumlah Kasus Tercatat", "Jumlah Penduduk",
                         "Insidensi", "Prevalensi", "Prevalensi_Kusta_Basah", "Prevalensi_Kusta_Kering")]
    datatable(df_display, options = list(pageLength = 15, scrollX = TRUE), rownames = FALSE) %>%
      formatRound(columns = c("Insidensi", "Prevalensi", "Prevalensi_Kusta_Basah", "Prevalensi_Kusta_Kering"), digits = 2)
  })
}

# 03b_server_part2.R - Server Logic (Part 2)

server_part2 <- function(input, output, session, std_df, tren_df) {
  
  # ============ TAB INTERPRETASI EPIDEMIOLOGI ============
  output$epi_summary <- renderUI({
    df <- std_df()
    if (is.null(df)) return(HTML("<i>Data tidak tersedia</i>"))
    
    top_incidence <- df %>% filter(!is.na(Insidensi)) %>% arrange(desc(Insidensi)) %>% slice(1)
    top_prevalensi <- df %>% filter(!is.na(Prevalensi)) %>% arrange(desc(Prevalensi)) %>% slice(1)
    top_pre_basah <- df %>% filter(!is.na(Prevalensi_Kusta_Basah)) %>% arrange(desc(Prevalensi_Kusta_Basah)) %>% slice(1)
    top_pre_kering <- df %>% filter(!is.na(Prevalensi_Kusta_Kering)) %>% arrange(desc(Prevalensi_Kusta_Kering)) %>% slice(1)
    
    HTML(paste0(
      "<div style='font-size:15px; line-height:2;'>",
      "<h4>Ukuran Frekuensi Kusta - Jawa Tengah 2024</h4>",
      "<table style='width:100%; border-collapse: collapse;'>",
      "<tr style='background:#f0f0f0;'><th style='padding:10px; text-align:left; border:1px solid #ddd;'>Ukuran</th>",
      "<th style='padding:10px; text-align:left; border:1px solid #ddd;'>Wilayah Tertinggi</th></tr>",
      "<tr><td style='padding:8px; border:1px solid #ddd;'><b>Insidensi</b></td>",
      "<td style='padding:8px; border:1px solid #ddd;'>", top_incidence$Wilayah, " (", round(top_incidence$Insidensi, 2), " per 100,000)</td></tr>",
      "<tr style='background:#f9f9f9;'><td style='padding:8px; border:1px solid #ddd;'><b>Prevalensi</b></td>",
      "<td style='padding:8px; border:1px solid #ddd;'>", top_prevalensi$Wilayah, " (", round(top_prevalensi$Prevalensi, 2), " per 100,000)</td></tr>",
      "<tr><td style='padding:8px; border:1px solid #ddd;'><b>Prevalensi Kasus Basah</b></td>",
      "<td style='padding:8px; border:1px solid #ddd;'>", top_pre_basah$Wilayah, " (", round(top_pre_basah$Prevalensi_Kusta_Basah, 2), " per 100,000)</td></tr>",
      "<tr style='background:#f9f9f9;'><td style='padding:8px; border:1px solid #ddd;'><b>Prevalensi Kasus Kering</b></td>",
      "<td style='padding:8px; border:1px solid #ddd;'>", top_pre_kering$Wilayah, " (", round(top_pre_kering$Prevalensi_Kusta_Kering, 2), " per 100,000)</td></tr>",
      "</table></div>"
    ))
  })
  
  output$epi_var_plot_ui <- renderUI({
    selectInput("epi_var_plot", "Pilih Variabel:", 
                choices = c("Insidensi" = "Insidensi", "Prevalensi" = "Prevalensi",
                            "Prevalensi Basah" = "Prevalensi_Kusta_Basah",
                            "Prevalensi Kering" = "Prevalensi_Kusta_Kering"),
                selected = "Insidensi")
  })
  
  output$epi_comparison_plot <- renderPlot({
    df <- std_df()
    req(df, input$epi_var_plot)
    
    var_name <- input$epi_var_plot
    plot_data <- df %>% filter(!is.na(!!sym(var_name))) %>%
      arrange(desc(!!sym(var_name))) %>% slice_head(n = 15)
    
    var_labels <- c("Insidensi" = "Insidensi (per 100,000 penduduk)",
                    "Prevalensi" = "Prevalensi (per 100,000 penduduk)",
                    "Prevalensi_Kusta_Basah" = "Prevalensi Kasus Basah (per 100,000 penduduk)",
                    "Prevalensi_Kusta_Kering" = "Prevalensi Kasus Kering (per 100,000 penduduk)")
    
    ggplot(plot_data, aes(x = reorder(Wilayah, !!sym(var_name)), y = !!sym(var_name))) +
      geom_col(fill = "#2C7873", alpha = 0.8) +
      geom_text(aes(label = round(!!sym(var_name), 2)), hjust = -0.2, size = 3.5) +
      coord_flip() + theme_minimal(base_size = 13) +
      labs(title = paste("Top 15 Wilayah berdasarkan", var_labels[var_name]),
           x = "Wilayah", y = var_labels[var_name]) +
      theme(plot.title = element_text(face = "bold", size = 15), panel.grid.major.y = element_blank())
  }, res = 96)
  
  output$epi_interpretation <- renderUI({
    df <- std_df()
    if (is.null(df)) return(HTML("<i>Data tidak tersedia</i>"))
    
    total_kasus_baru <- sum(df$`Jumlah Kasus Baru`, na.rm = TRUE)
    total_penduduk <- sum(df$`Jumlah Penduduk`, na.rm = TRUE)
    total_kasus_tercatat <- sum(df$`Jumlah Kasus Tercatat`, na.rm = TRUE)
    total_basah <- sum(df$`Jumlah Kasus Kusta Basah`, na.rm = TRUE)
    
    incidence_agregat <- round((total_kasus_baru / total_penduduk) * 100000, 2)
    prevalensi_agregat <- round((total_kasus_tercatat / total_penduduk) * 100000, 2)
    pct_basah <- round((total_basah / total_kasus_tercatat) * 100, 1)
    
    HTML(paste0(
      "<div style='font-size:14px; line-height:1.9;'>",
      "<h4>Interpretasi Epidemiologi</h4>",
      "<p><b>Jumlah Kasus Baru Kusta:</b> ", comma(total_kasus_baru), " kasus</p>",
      "<p><b>Insidensi (Agregat):</b> ", incidence_agregat, " per 100,000 penduduk</p>",
      "<p><b>Prevalensi (Agregat):</b> ", prevalensi_agregat, " per 100,000 penduduk</p>",
      "<hr>",
      "<p style='text-align:justify;'><i><b>Interpretasi:</b> Angka insidensi sebesar ", incidence_agregat, 
      " per 100,000 penduduk menunjukkan bahwa dari setiap 100,000 penduduk di Jawa Tengah, terdapat sekitar ", 
      round(incidence_agregat, 0), " kasus baru kusta yang terdiagnosis pada tahun 2024. ",
      "Angka prevalensi sebesar ", prevalensi_agregat, " per 100,000 penduduk menunjukkan beban penyakit kusta bahwa dari setiap 100,000 penduduk, terdapat sekitar ",
      round(prevalensi_agregat, 0), " kasus kusta yang tercatat pada tahun 2024. Dari total kasus tercatat, ", pct_basah, 
      "% merupakan kasus kusta basah (multibacillary), yang menandakan bahwa transmisi penyakit masih aktif ",
      "dan memerlukan upaya pencegahan dan deteksi dini yang lebih intensif untuk mengurangi penyebaran penyakit.</i></p>",
      "</div>"
    ))
  })
  
  # ============ TAB PETA VISUALISASI ============
  output$map_var_ui <- renderUI({
    df <- std_df()
    req(df)
    
    num_vars <- names(df)[sapply(df, is.numeric)]
    priority_vars <- c("Insidensi", "Prevalensi", "Jumlah Kasus Baru", "Jumlah Kasus Tercatat",
                       "Persentase Akses Sanitasi Layak", "Persentase Rumah Layak Huni")
    available_priority <- intersect(priority_vars, num_vars)
    other_vars <- setdiff(num_vars, priority_vars)
    ordered_vars <- c(available_priority, other_vars)
    
    if (length(ordered_vars) == 0) {
      return(tags$div("Tidak ada variabel numerik untuk divisualisasikan."))
    }
    selectInput("map_var", "Variabel:", choices = ordered_vars, selected = ordered_vars[1])
  })
  
  output$map_year_ui <- renderUI({
    df_tren <- tren_df()
    req(input$map_var)
    
    if(input$map_var == "Jumlah Kasus Baru" && !is.null(df_tren)) {
      year_cols <- names(df_tren)[names(df_tren) != "Wilayah"]
      selectInput("map_year", "Pilih Tahun:", 
                  choices = c("2024", year_cols), selected = "2024")
    } else {
      NULL
    }
  })
  
  output$map_leaflet <- renderLeaflet({
    req(input$map_var)
    df <- std_df()
    df_tren <- tren_df()
    req(df)
    
    if(input$map_var == "Jumlah Kasus Baru" && !is.null(df_tren) && !is.null(input$map_year)) {
      if(input$map_year != "2024") {
        df_for_map <- data.frame(
          Wilayah = df_tren$Wilayah,
          value = df_tren[[input$map_year]]
        )
      } else {
        df_for_map <- data.frame(
          Wilayah = df$Wilayah,
          value = df[[input$map_var]]
        )
      }
    } else {
      df_for_map <- data.frame(
        Wilayah = df$Wilayah,
        value = df[[input$map_var]]
      )
    }
    
    df_map <- merge(coords_jateng, df_for_map, by = "Wilayah", all.x = TRUE)
    
    vals <- df_map$value
    valid_vals <- vals[!is.na(vals)]
    
    if (length(valid_vals) == 0) {
      return(leaflet() %>% addTiles() %>% setView(lng = 110.42, lat = -7.15, zoom = 8))
    }
    
    output$map_info <- renderUI({
      year_info <- if(input$map_var == "Jumlah Kasus Baru" && !is.null(input$map_year)) {
        paste0("<li><strong>Tahun:</strong> ", input$map_year, "</li>")
      } else ""
      
      HTML(paste0(
        "<h5>Statistik Peta:</h5>",
        "<ul style='font-size:13px;'>", year_info,
        "<li><strong>Min:</strong> ", comma(round(min(valid_vals, na.rm=TRUE), 2)), "</li>",
        "<li><strong>Maks:</strong> ", comma(round(max(valid_vals, na.rm=TRUE), 2)), "</li>",
        "<li><strong>Rata-rata:</strong> ", comma(round(mean(valid_vals, na.rm=TRUE), 2)), "</li>",
        "</ul>"
      ))
    })
    
    pal <- colorNumeric(
      palette = colorRampPalette(c(input$color_low, input$color_high))(10),
      domain = df_map$value,
      na.color = "#cccccc"
    )
    
    max_val <- max(valid_vals, na.rm = TRUE)
    df_map$radius <- ifelse(is.na(df_map$value), 5, 5 + (df_map$value / max_val) * 15)
    
    leaflet(df_map) %>%
      addProviderTiles("CartoDB.Positron") %>%
      setView(lng = 110.42, lat = -7.15, zoom = 8) %>%
      addCircleMarkers(
        lng = ~lng, lat = ~lat, radius = ~radius,
        color = ~pal(value), fillColor = ~pal(value),
        fillOpacity = 0.7, stroke = TRUE, weight = 2,
        popup = ~paste0("<b>", Wilayah, "</b><br>",
                        input$map_var, ": ", 
                        ifelse(is.na(value), "Tidak ada data", round(value, 2)))
      ) %>%
      addLegend(
        position = "bottomright", pal = pal, values = ~value,
        title = input$map_var, opacity = 0.7, na.label = "Tidak ada data"
      )
  })
  
  output$spatial_interpretation <- renderUI({
    req(input$map_var)
    df <- std_df()
    req(df)
    
    var_name <- input$map_var
    
    if(!(var_name %in% names(df))) return(HTML("<p>Variabel tidak tersedia untuk interpretasi.</p>"))
    
    vals <- df[[var_name]]
    valid_vals <- vals[!is.na(vals)]
    
    if (length(valid_vals) == 0) {
      return(HTML("<p>Tidak dapat memberikan interpretasi karena tidak ada data valid.</p>"))
    }
    
    q <- quantile(valid_vals, probs = c(0.25, 0.75), na.rm = TRUE)
    high_threshold <- q[2]
    low_threshold <- q[1]
    
    high_regions <- df[!is.na(vals) & vals >= high_threshold, ]
    low_regions <- df[!is.na(vals) & vals <= low_threshold, ]
    
    high_names <- if(nrow(high_regions) > 0) {
      paste(head(high_regions$Wilayah, 5), collapse = ", ")
    } else "Tidak ada"
    
    low_names <- if(nrow(low_regions) > 0) {
      paste(head(low_regions$Wilayah, 5), collapse = ", ")
    } else "Tidak ada"
    
    HTML(paste0(
      "<h4>Interpretasi Pola Spasial: ", tools::toTitleCase(gsub("_", " ", var_name)), "</h4>",
      "<p><strong>Wilayah TINGGI</strong> (Q3 = ", comma(round(high_threshold, 2)), "): ", high_names, "</p>",
      "<p><strong>Wilayah RENDAH</strong> (Q1 = ", comma(round(low_threshold, 2)), "): ", low_names, "</p>",
      "<p style='text-align:justify;'><i>Ukuran dan warna lingkaran pada peta menunjukkan besaran nilai. ",
      "Wilayah dengan circle lebih besar dan warna lebih tua memiliki nilai lebih tinggi.</i></p>"
    ))
  })
  
  # ============ TAB TREN WAKTU ============
  output$tren_plot <- renderPlotly({
    df_tren <- tren_df()
    req(df_tren)
    
    year_cols <- names(df_tren)[names(df_tren) != "Wilayah"]
    df_long <- df_tren[, c("Wilayah", year_cols)] %>%
      pivot_longer(cols = all_of(year_cols), names_to = "Tahun", values_to = "Kasus")
    
    df_summary <- df_long %>% group_by(Tahun) %>% summarise(Total_Kasus = sum(Kasus, na.rm = TRUE))
    
    plot_ly(df_summary, x = ~Tahun, y = ~Total_Kasus, type = "scatter", mode = "lines+markers",
            line = list(color = "#2C7873", width = 3), marker = list(size = 10, color = "#2C7873"),
            text = ~paste("Tahun:", Tahun, "<br>Total:", comma(Total_Kasus)), hoverinfo = "text",
            name = "Total Provinsi") %>%
      layout(title = "Tren Kasus Baru Kusta Provinsi Jawa Tengah 2021-2024",
             xaxis = list(title = "Tahun"), yaxis = list(title = "Jumlah Kasus Baru"), showlegend = TRUE)
  })
  
  output$tren_wilayah_ui <- renderUI({
    df_tren <- tren_df()
    req(df_tren)
    
    selectizeInput("tren_wilayah", "Pilih Wilayah (bisa lebih dari 1):", 
                   choices = df_tren$Wilayah, selected = df_tren$Wilayah[1],
                   multiple = TRUE, options = list(maxItems = 10))
  })
  
  output$tren_wilayah_plot <- renderPlotly({
    df_tren <- tren_df()
    req(df_tren, input$tren_wilayah)
    
    if(length(input$tren_wilayah) == 0) return(NULL)
    
    year_cols <- names(df_tren)[names(df_tren) != "Wilayah"]
    
    df_selected <- df_tren %>% filter(Wilayah %in% input$tren_wilayah)
    df_selected <- df_selected[, c("Wilayah", year_cols)] %>%
      pivot_longer(cols = all_of(year_cols), names_to = "Tahun", values_to = "Kasus")
    
    df_avg <- df_tren[, year_cols] %>%
      summarise(across(everything(), ~mean(.x, na.rm = TRUE))) %>%
      pivot_longer(everything(), names_to = "Tahun", values_to = "Kasus") %>%
      mutate(Wilayah = "Rata-rata Provinsi")
    
    df_plot <- bind_rows(df_selected, df_avg)
    
    p <- plot_ly()
    
    for(wilayah in input$tren_wilayah) {
      df_w <- df_plot %>% filter(Wilayah == wilayah)
      p <- p %>% add_trace(data = df_w, x = ~Tahun, y = ~Kasus, type = "scatter", mode = "lines+markers",
                           name = wilayah, line = list(width = 2), marker = list(size = 8),
                           text = ~paste(Wilayah, "<br>Tahun:", Tahun, "<br>Kasus:", Kasus), hoverinfo = "text")
    }
    
    df_avg_plot <- df_plot %>% filter(Wilayah == "Rata-rata Provinsi")
    p <- p %>% add_trace(data = df_avg_plot, x = ~Tahun, y = ~Kasus, type = "scatter", mode = "lines+markers",
                         name = "Rata-rata Provinsi", line = list(width = 3, dash = "dash", color = "black"),
                         marker = list(size = 10, color = "black", symbol = "diamond"),
                         text = ~paste("Rata-rata Provinsi<br>Tahun:", Tahun, "<br>Rata-rata:", round(Kasus, 1)),
                         hoverinfo = "text")
    
    p %>% layout(title = "Perbandingan Tren Kasus Baru Kusta per Wilayah",
                 xaxis = list(title = "Tahun"), yaxis = list(title = "Jumlah Kasus Baru"),
                 showlegend = TRUE, legend = list(x = 1.05, y = 1))
  })
  
  output$tren_table <- renderDT({
    df_tren <- tren_df()
    req(df_tren)
    
    year_cols <- names(df_tren)[names(df_tren) != "Wilayah"]
    df_display <- df_tren[, c("Wilayah", year_cols)]
    datatable(df_display, options = list(pageLength = 15, scrollX = TRUE), rownames = FALSE)
  })
  
  output$tren_interpretation <- renderUI({
    df_tren <- tren_df()
    req(df_tren)
    
    year_cols <- names(df_tren)[names(df_tren) != "Wilayah"]
    first_year <- year_cols[1]
    last_year <- year_cols[length(year_cols)]
    
    interpretations <- c()
    
    for(i in 1:nrow(df_tren)) {
      wilayah <- df_tren$Wilayah[i]
      val_first <- df_tren[[first_year]][i]
      val_last <- df_tren[[last_year]][i]
      
      if(!is.na(val_first) && !is.na(val_last) && val_first > 0) {
        pct_change <- round(((val_last - val_first) / val_first) * 100, 1)
        
        if(abs(pct_change) < 10) {
          trend <- "Tren Stabil"
          color <- "blue"
        } else if(pct_change > 0) {
          trend <- "Tren Meningkat"
          color <- "red"
        } else {
          trend <- "Tren Menurun"
          color <- "green"
        }
        
        interpretations <- c(interpretations, paste0(
          "<li style='color:", color, ";'><b>", wilayah, "</b>: ", trend, 
          " (", ifelse(pct_change > 0, "+", ""), pct_change, "%)</li>"
        ))
      }
    }
    
    HTML(paste0(
      "<div style='font-size:14px; line-height:1.8;'>",
      "<h4>Interpretasi Tren ", first_year, "-", last_year, "</h4>",
      "<ul>", paste(head(interpretations, 10), collapse = ""), "</ul>",
      "</div>"
    ))
  })
}

# 03c_server_part3.R - Server Logic (Modelling Section)

server_part3 <- function(input, output, session, model_data) {
  
  output$vb_jumlah_observasi <- renderValueBox({
    df <- model_data()
    val <- if(is.null(df)) 0 else nrow(df)
    valueBox(val, "Jumlah Observasi", icon = icon("database"), color = "aqua")
  })
  
  output$vb_variabel_y <- renderValueBox({
    valueBox(1, "Jumlah Variabel Y", icon = icon("bullseye"), color = "yellow")
  })
  
  output$vb_variabel_x <- renderValueBox({
    df <- model_data()
    val <- if(is.null(df)) 0 else (ncol(df) - 1)
    valueBox(val, "Jumlah Variabel X", icon = icon("list"), color = "green")
  })
  
  output$model_info <- renderUI({
    df <- model_data()
    req(df)
    x_vars <- names(df)[names(df) != "Y"]
    HTML(paste0(
      "<div style='font-size:15px; line-height:2;'>",
      "<p><b>Dataset:</b> 5px; line-height:2;'>",
      "<p><b>Dataset:</b> Data yang digunakan merupakan data sekunder yang berkaitan dengan Kasus Baru Penyakit Kusta di Provinsi Jawa Tengah dengan 35 kabupaten/kota tahun 2024 yang bersumber dari sumber resmi seperti Badan Pusat Statistik (BPS) dan Open Data Provinsi Jawa Tengah.</p>",
      "<p><b>Variabel Y:</b> Jumlah Kasus Baru</p>",
      "<p><b>Variabel X:</b> ", paste(x_vars, collapse = ", "), "</p>",
      "</div>"
    ))
  })
  
  output$uji_asumsi <- renderUI({
    df <- model_data()
    req(df)
    
    model_poisson <- glm(Y ~ ., data = df, family = poisson(link = "log"))
    residual_deviance <- model_poisson$deviance
    df_residual <- model_poisson$df.residual
    dispersion_ratio <- residual_deviance / df_residual
    
    interpretation <- if(dispersion_ratio > 1) {
      paste0("<span style='color:red;'><b>Terdapat overdispersi</b></span> (ratio = ", 
             round(dispersion_ratio, 2), " > 1), asumsi model Poisson tidak terpenuhi. Menggunakan model Binomial Negatif lebih sesuai.")
    } else {
      paste0("<span style='color:green;'><b>Tidak terdapat overdispersi</b></span> (ratio = ", 
             round(dispersion_ratio, 2), " ≤ 1). Model Poisson dapat digunakan.")
    }
    
    HTML(paste0(
      "<div style='font-size:14px; line-height:1.8;'>",
      "<h5>Uji Overdispersi</h5>",
      "<p><b>Residual Deviance:</b> ", round(residual_deviance, 2), "</p>",
      "<p><b>DF:</b> ", df_residual, "</p>",
      "<p><b>Dispersion Ratio:</b> ", round(dispersion_ratio, 2), "</p>",
      "<p>", interpretation, "</p>",
      "</div>"
    ))
  })
  
  output$uji_multikol <- renderUI({
    df <- model_data()
    req(df)
    
    model_nb <- glm.nb(Y ~ ., data = df)
    vif_values <- vif(model_pois)
    
    vif_df <- data.frame(
      Variabel = names(vif_values),
      VIF = round(vif_values, 2),
      Status = ifelse(vif_values > 10, "Multikolinearitas Tinggi", 
                      ifelse(vif_values > 5, "Multikolinearitas Sedang", "Tidak ada"))
    )
    
    vif_table <- paste0(
      "<table style='width:100%; border-collapse: collapse; margin-top:10px;'>",
      "<tr style='background:#f0f0f0;'><th style='padding:8px; border:1px solid #ddd;'>Variabel</th>",
      "<th style='padding:8px; border:1px solid #ddd;'>VIF</th>",
      "<th style='padding:8px; border:1px solid #ddd;'>Status</th></tr>",
      paste(apply(vif_df, 1, function(row) {
        color <- if(as.numeric(row[2]) > 10) "red" else if(as.numeric(row[2]) > 5) "orange" else "green"
        paste0("<tr><td style='padding:8px; border:1px solid #ddd;'>", row[1], "</td>",
               "<td style='padding:8px; border:1px solid #ddd;'>", row[2], "</td>",
               "<td style='padding:8px; border:1px solid #ddd; color:", color, ";'><b>", row[3], "</b></td></tr>")
      }), collapse = ""),
      "</table>"
    )
    
    HTML(paste0("<div style='font-size:14px; line-height:1.8;'>",
                "<h5>Uji Multikolinearitas (VIF)</h5>",
                vif_table, "</div>"))
  })
  
  output$model_full_summary <- renderUI({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    coef_summary <- summary(model_full)$coefficients
    
    coef_df <- data.frame(
      Variabel = rownames(coef_summary),
      Estimate = round(coef_summary[, 1], 4),
      p_value = round(coef_summary[, 4], 4),
      Sig = ifelse(coef_summary[, 4] < 0.001, "***",
                   ifelse(coef_summary[, 4] < 0.01, "**",
                          ifelse(coef_summary[, 4] < 0.05, "*", "ns")))
    )
    
    coef_table <- paste0(
      "<table style='width:100%; border-collapse: collapse; font-size:11px;'>",
      "<tr style='background:#f0f0f0;'><th style='padding:5px; border:1px solid #ddd;'>Variabel</th>",
      "<th style='padding:5px; border:1px solid #ddd;'>Estimate</th>",
      "<th style='padding:5px; border:1px solid #ddd;'>P-value</th>",
      "<th style='padding:5px; border:1px solid #ddd;'>Sig</th></tr>",
      paste(apply(coef_df, 1, function(row) {
        paste0("<tr><td style='padding:5px; border:1px solid #ddd; font-size:10px;'>", row[1], "</td>",
               "<td style='padding:5px; border:1px solid #ddd;'>", row[2], "</td>",
               "<td style='padding:5px; border:1px solid #ddd;'>", row[3], "</td>",
               "<td style='padding:5px; border:1px solid #ddd;'><b>", row[4], "</b></td></tr>")
      }), collapse = ""),
      "</table>"
    )
    
    HTML(paste0("<div style='font-size:12px;'>", 
                "<p><b>AIC:</b> ", round(AIC(model_full), 2), "</p>",
                coef_table, "</div>"))
  })
  
  output$model_full_note <- renderUI({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    coef_summary <- summary(model_full)$coefficients
    non_sig <- rownames(coef_summary)[coef_summary[, 4] >= 0.05 & rownames(coef_summary) != "(Intercept)"]
    
    if(length(non_sig) > 0) {
      HTML(paste0("<b>Variabel tidak signifikan:</b> ", paste(non_sig, collapse = ", ")))
    } else {
      HTML("<b>Semua variabel signifikan</b>")
    }
  })
  
  output$model_sign_summary <- renderUI({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    coef_summary <- summary(model_full)$coefficients
    sig_vars <- rownames(coef_summary)[coef_summary[, 4] < 0.05 & rownames(coef_summary) != "(Intercept)"]
    
    if(length(sig_vars) > 0) {
      formula_sig <- as.formula(paste("Y ~", paste(sig_vars, collapse = " + ")))
      model_sig <- glm.nb(formula_sig, data = df)
      coef_summary_sig <- summary(model_sig)$coefficients
      
      coef_df <- data.frame(
        Variabel = rownames(coef_summary_sig),
        Estimate = round(coef_summary_sig[, 1], 4),
        p_value = round(coef_summary_sig[, 4], 4),
        Sig = ifelse(coef_summary_sig[, 4] < 0.001, "***",
                     ifelse(coef_summary_sig[, 4] < 0.01, "**",
                            ifelse(coef_summary_sig[, 4] < 0.05, "*", "ns")))
      )
      
      coef_table <- paste0(
        "<table style='width:100%; border-collapse: collapse; font-size:11px;'>",
        "<tr style='background:#f0f0f0;'><th style='padding:5px; border:1px solid #ddd;'>Variabel</th>",
        "<th style='padding:5px; border:1px solid #ddd;'>Estimate</th>",
        "<th style='padding:5px; border:1px solid #ddd;'>P-value</th>",
        "<th style='padding:5px; border:1px solid #ddd;'>Sig</th></tr>",
        paste(apply(coef_df, 1, function(row) {
          paste0("<tr><td style='padding:5px; border:1px solid #ddd; font-size:10px;'>", row[1], "</td>",
                 "<td style='padding:5px; border:1px solid #ddd;'>", row[2], "</td>",
                 "<td style='padding:5px; border:1px solid #ddd;'>", row[3], "</td>",
                 "<td style='padding:5px; border:1px solid #ddd;'><b>", row[4], "</b></td></tr>")
        }), collapse = ""),
        "</table>"
      )
      
      HTML(paste0("<div style='font-size:12px;'>", 
                  "<p><b>AIC:</b> ", round(AIC(model_sig), 2), "</p>",
                  coef_table, "</div>"))
    } else {
      HTML("<p>Tidak ada variabel signifikan.</p>")
    }
  })
  
  output$model_step_summary <- renderUI({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    model_step <- stepAIC(model_full, direction = "both", trace = FALSE)
    coef_summary_step <- summary(model_step)$coefficients
    
    coef_df <- data.frame(
      Variabel = rownames(coef_summary_step),
      Estimate = round(coef_summary_step[, 1], 4),
      p_value = round(coef_summary_step[, 4], 4),
      Sig = ifelse(coef_summary_step[, 4] < 0.001, "***",
                   ifelse(coef_summary_step[, 4] < 0.01, "**",
                          ifelse(coef_summary_step[, 4] < 0.05, "*", "ns")))
    )
    
    coef_table <- paste0(
      "<table style='width:100%; border-collapse: collapse; font-size:11px;'>",
      "<tr style='background:#f0f0f0;'><th style='padding:5px; border:1px solid #ddd;'>Variabel</th>",
      "<th style='padding:5px; border:1px solid #ddd;'>Estimate</th>",
      "<th style='padding:5px; border:1px solid #ddd;'>P-value</th>",
      "<th style='padding:5px; border:1px solid #ddd;'>Sig</th></tr>",
      paste(apply(coef_df, 1, function(row) {
        paste0("<tr><td style='padding:5px; border:1px solid #ddd; font-size:10px;'>", row[1], "</td>",
               "<td style='padding:5px; border:1px solid #ddd;'>", row[2], "</td>",
               "<td style='padding:5px; border:1px solid #ddd;'>", row[3], "</td>",
               "<td style='padding:5px; border:1px solid #ddd;'><b>", row[4], "</b></td></tr>")
      }), collapse = ""),
      "</table>"
    )
    
    HTML(paste0("<div style='font-size:12px;'>", 
                "<p><b>AIC:</b> ", round(AIC(model_step), 2), "</p>",
                coef_table, "</div>"))
  })
  
  output$model_comparison <- renderDT({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    coef_summary <- summary(model_full)$coefficients
    sig_vars <- rownames(coef_summary)[coef_summary[, 4] < 0.05 & rownames(coef_summary) != "(Intercept)"]
    
    if(length(sig_vars) > 0) {
      formula_sig <- as.formula(paste("Y ~", paste(sig_vars, collapse = " + ")))
      model_sig <- glm.nb(formula_sig, data = df)
    } else {
      model_sig <- model_full
    }
    
    model_step <- stepAIC(model_full, direction = "both", trace = FALSE)
    
    mcfadden_r2 <- function(model) {
      null_model <- glm.nb(Y ~ 1, data = df)
      1 - (logLik(model)[1] / logLik(null_model)[1])
    }
    
    comparison_df <- data.frame(
      Model = c("Full", "Signifikan", "StepAIC"),
      AIC = c(AIC(model_full), AIC(model_sig), AIC(model_step)),
      BIC = c(BIC(model_full), BIC(model_sig), BIC(model_step)),
      McFadden_R2 = c(mcfadden_r2(model_full), mcfadden_r2(model_sig), mcfadden_r2(model_step))
    )
    
    datatable(comparison_df, options = list(dom = 't', pageLength = 3), rownames = FALSE) %>%
      formatRound(columns = c("AIC", "BIC", "McFadden_R2"), digits = 3)
  })
  
  output$model_comparison_interp <- renderUI({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    model_step <- stepAIC(model_full, direction = "both", trace = FALSE)
    
    HTML(paste0(
      "<b>Interpretasi:</b> Berdasarkan nilai AIC dan BIC terkecil, model terpilih adalah <b>Model StepAIC</b> dengan AIC = ", 
      round(AIC(model_step), 2), " dan BIC = ", round(BIC(model_step), 2), 
      ". Model ini memiliki nilai AIC dan BIC terkecil yang menunjukkan keseimbangan terbaik antara kesesuaian model (goodness of fit) dan kompleksitas model. Meskipun Full Model memiliki McFadden R² yang sedikit lebih tinggi dari model terpilih(0,106 vs 0,098), selisih nilai tersebut relatif kecil. Hal ini menunjukkan bahwa penambahan variabel pada Full Model hanya memberikan peningkatan kemampuan penjelasan yang terbatas, sementara nilai AIC dan BIC yang lebih rendah pada model Stepwise mengindikasikan bahwa model ini lebih efisien dalam menjelaskan data."
    ))
  })
  
  output$model_final_equation <- renderUI({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    model_step <- stepAIC(model_full, direction = "both", trace = FALSE)
    coefs <- coef(model_step)
    var_names <- names(coefs)[-1]
    
    equation_parts <- paste0(round(coefs[-1], 4), " × ", var_names)
    equation <- paste0("log(Y) = ", round(coefs[1], 4), " + ", paste(equation_parts, collapse = " + "))
    
    HTML(paste0(
      "<div style='font-size:14px; line-height:1.8; background:#f0f9ff; padding:15px; border-radius:5px;'>",
      "<h5>Bentuk Model Terpilih (StepAIC)</h5>",
      "<p style='font-family:monospace; background:white; padding:10px; border-radius:3px;'>", 
      equation, "</p>",
      "</div>"
    ))
  })
  
  output$model_interpretation <- renderUI({
    df <- model_data()
    req(df)
    
    model_full <- glm.nb(Y ~ ., data = df)
    model_step <- stepAIC(model_full, direction = "both", trace = FALSE)
    coefs <- coef(model_step)
    sig_coefs <- summary(model_step)$coefficients
    
    interpretations <- c()
    for(i in 2:length(coefs)) {
      var_name <- names(coefs)[i]
      estimate <- coefs[i]
      p_value <- sig_coefs[i, 4]
      
      effect <- ifelse(estimate > 0, "meningkat", "menurun")
      pct_change <- round((exp(estimate) - 1) * 100, 2)
      
      interpretations <- c(interpretations, paste0(
        "<li><b>", var_name, "</b>: Setiap kenaikan 1 unit ", var_name,
        " maka jumlah kasus baru kusta akan", effect, 
        " sebesar ", abs(pct_change), "% (p-value = ", round(p_value, 4), ")</li>"
      ))
    }
    
    HTML(paste0(
      "<div style='font-size:14px; line-height:1.8; margin-top:15px;'>",
      "<h5>Interpretasi Koefisien</h5>",
      "<ul>", paste(interpretations, collapse = ""), "</ul>",
      "</div>"
    ))
  })
}

server <- function(input, output, session) {
  # DEBUG: Print data ke console
  cat("\n=== CHECKING DATA ===\n")
  cat("initial_data exists:", exists("initial_data"), "\n")
  cat("initial_data is null:", is.null(initial_data), "\n")
  if (!is.null(initial_data)) {
    cat("Number of rows:", nrow(initial_data), "\n")
    cat("Total Kasus Baru:", sum(initial_data$`Jumlah Kasus Baru`, na.rm = TRUE), "\n")
  }
  
  # Reactive data
  std_df <- reactive({ initial_data })
  tren_df <- reactive({ data_tren })
  
  # Model data
  model_data <- reactive({
    df <- std_df()
    req(df)
    
    df_model <- data.frame(
      Y = df$`Jumlah Kasus Baru`,
      Penduduk_Miskin = df$`Persentase Penduduk Miskin`,
      Rumah_Layak_Huni = df$`Persentase Rumah Layak Huni`,
      Akses_Sanitasi_Layak = df$`Persentase Akses Sanitasi Layak`,
      Akses_Air_Minum_Layak = df$`Persentase Akses Air Minum Layak`,
      Jumlah_Fasilitas_Kesehatan = df$`Jumlah Fasilitas Kesehatan`,
      Kepadatan_Penduduk = df$`Kepadatan Penduduk (per km^2)`,
      Rata_Rata_Lama_Sekolah = df$`Rata-Rata Lama Sekolah (IPM) (Tahun)`
    )
    
    df_model[complete.cases(df_model), ]
  })
  # Panggil semua fungsi server part
  server_part1(input, output, session, std_df, tren_df)
  server_part2(input, output, session, std_df, tren_df)
  server_part3(input, output, session, model_data)
}

# Run App
shinyApp(ui = ui, server = server)

ACKNOWLEDGEMENT

AI digunakan untuk membantu membuat syntax dashboard agar tampilannya terlihat lebih menarik, membantu memunculkan ide, sebagai partner diskusi jika ditemukan hal-hal atau pernyataan di jurnal yang kemudian memunculkan banyak pertanyaan baru di otak saya, serta membantu membuat kalimat interpretasi yang efektif dan mudah dimengerti saat membacanya. Adapun untuk interpretasi aslinya merupakan interpretasi yang saya pahami sendiri dan kemudian dijabarkan kalimatnya agar lebih dapat dipahami.