Analisis Spasio-Temporal Kasus Stunting Menggunakan Generalized Linear Mixed Model Negative Binomial di Provinsi Nusa Tenggara Timur Tahun 2021-2025

Nama	Raihanah Rafidah
NPM	140720250009
Program Studi	Magister Statistikan Terapan
Mata Kuliah	Epidemiologi

ABSTRAK

Stunting merupakan salah satu permasalahan kesehatan masyarakat yang memiliki pola penyebaran dipengaruhi oleh aspek spasial dan temporal. Penelitian ini bertujuan menganalisis pola spasio-temporal kasus stunting di Provinsi Nusa Tenggara Timur (NTT) periode 2021-2025 serta membangun model yang mampu mengakomodasi overdispersi dan heterogenitas wilayah. Data yang digunakan adalah data sekunder jumlah kasus stunting dan populasi balita pada 21 kabupaten/kota di Provinsi NTT selama tahun 2021-2025. Analisis dilakukan melalui statistik deskriptif, Incidence Rate (IR) dan Standardized Incidence Ratio (SIR), pembentukan matriks ketetanggaan menggunakan K-Nearest Neighbor (KNN), pengujian autokorelasi spasio-temporal dengan Global Moran’s I dan Local Moran’s I (LISA), uji overdispersi, serta pemodelan Negative Binomial dan Generalized Linear Mixed Model (GLMM) Negative Binomial. Hasil penelitian menunjukkan adanya autokorelasi spasial yang signifikan, yaitu Moran’s I sebesar 0,189 dengan beberapa klaster High-High dan Low-Low. Data juga mengalami overdispersi (φ) sebesar 532, sehingga model Negative Binomial lebih sesuai dibandingkan Poisson. Selanjutnya, model GLMM Negative Binomial berhasil mengidentifikasi pengaruh spasial yang signifikan melalui variabel lag kasus (p = 0,0319) setelah memasukkan efek acak kabupaten/kota dan tahun. Hasil penelitian menunjukkan bahwa pendekatan GLMM Negative Binomial mampu menggambarkan karakteristik spasio-temporal kasus stunting secara lebih komprehensif dan dapat menjadi dasar dalam penyusunan kebijakan penanggulangan stunting yang lebih tepat sasaran.

Kata kunci: stunting, spasio-temporal, Global Moran’s I, LISA, GLMM Negative Binomial, Nusa Tenggara Timur.

PENDAHULUAN

Stunting merupakan salah satu permasalahan kesehatan masyarakat yang menjadi perhatian global karena berdampak terhadap kualitas sumber daya manusia, produktivitas ekonomi, serta pembangunan berkelanjutan. Menurut WHO, stunting adalah kondisi gagal tumbuh pada anak akibat kekurangan gizi kronis yang ditandai dengan nilai height-for-age z-score kurang dari -2 standar deviasi berdasarkan standar pertumbuhan anak WHO. Dampak stunting tidak hanya terjadi pada masa kanak-kanak, tetapi juga berpengaruh terhadap perkembangan kognitif, prestasi belajar, kapasitas kerja, hingga meningkatkan risiko penyakit tidak menular pada usia dewasa. Oleh karena itu, penurunan prevalensi stunting menjadi salah satu indikator penting dalam pencapaian Sustainable Development Goals (SDGs), khususnya tujuan kedua (Zero Hunger) dan tujuan ketiga (Good Health and Well-being).

Indonesia masih menghadapi tantangan yang cukup besar dalam upaya penanggulangan stunting. Berdasarkan hasil Survei Status Gizi Indonesia (SSGI) dan Survei Kesehatan Indonesia (SKI), prevalensi stunting nasional mengalami penurunan dari 24,4% pada tahun 2021 menjadi 21,6% pada tahun 2022 dan kembali menurun menjadi 21,5% pada tahun 2023. Meskipun menunjukkan tren yang membaik, angka tersebut masih berada di atas ambang batas yang ditetapkan WHO, yaitu kurang dari 20%. Pemerintah Indonesia sendiri menargetkan prevalensi stunting turun hingga 14% sebagai bagian dari strategi percepatan penurunan stunting nasional. Namun, pencapaian target tersebut masih menghadapi berbagai tantangan, terutama karena adanya ketimpangan kondisi antarwilayah.

Salah satu provinsi yang masih memiliki beban stunting cukup tinggi adalah Provinsi Nusa Tenggara Timur (NTT). Karakteristik geografis berupa wilayah kepulauan, keterbatasan akses pelayanan kesehatan, tingkat kemiskinan yang relatif tinggi, serta perbedaan kondisi sosial ekonomi antar kabupaten/kota menyebabkan distribusi kasus stunting di provinsi ini tidak merata. Berdasarkan data yang dianalisis dalam penelitian ini, jumlah kasus stunting di Provinsi NTT menunjukkan pola yang dinamis selama periode 2021–2025. Total kasus tercatat sebanyak 78.124 kasus pada tahun 2021, menurun menjadi 74.756 kasus pada tahun 2022, kemudian mencapai titik terendah sebesar 61.395 kasus pada tahun 2023. Namun demikian, jumlah kasus kembali meningkat menjadi 62.218 kasus pada tahun 2024 dan 62.713 kasus pada tahun 2025. Selain mengalami perubahan secara temporal, distribusi kasus juga menunjukkan variasi yang cukup besar antar kabupaten/kota, sehingga mengindikasikan adanya karakteristik spasial yang perlu dikaji lebih lanjut.

Dalam analisis epidemiologi, data jumlah kasus penyakit umumnya berupa data cacah (count data) yang sering kali tidak memenuhi asumsi distribusi Poisson akibat terjadinya overdispersi, yaitu kondisi ketika varians lebih besar daripada nilai rata-rata. Selain itu, distribusi kasus penyakit sering memperlihatkan autokorelasi spasial, yaitu kecenderungan wilayah yang berdekatan memiliki tingkat kejadian penyakit yang serupa karena dipengaruhi oleh karakteristik lingkungan, sosial ekonomi, maupun akses pelayanan kesehatan yang relatif sama. Apabila karakteristik tersebut diabaikan, model statistik yang dihasilkan berpotensi memberikan estimasi parameter yang bias serta kesimpulan inferensial yang kurang akurat. Oleh karena itu, diperlukan pendekatan analisis yang mampu mengakomodasi data cacah, overdispersi, serta ketergantungan spasial dan temporal secara simultan.

Beberapa penelitian terdahulu telah menerapkan berbagai metode statistik untuk menganalisis kasus stunting maupun penyakit berbasis wilayah. Regresi Poisson banyak digunakan untuk menganalisis data kejadian penyakit, namun memiliki keterbatasan ketika data mengalami overdispersi. Untuk mengatasi permasalahan tersebut, Hilbe (2011) dan Cameron dan Trivedi (2013) merekomendasikan penggunaan regresi Negative Binomial karena mampu mengakomodasi varians yang lebih besar daripada rata-rata. Penelitian lain memanfaatkan Global Moran’s I dan Local Moran’s I (LISA) untuk mengidentifikasi pola autokorelasi spasial dan klaster wilayah berisiko tinggi, sehingga dapat diketahui adanya pengelompokan kejadian penyakit secara geografis. Sementara itu, perkembangan metode statistik modern menunjukkan bahwa Generalized Linear Mixed Model (GLMM) mampu mengakomodasi heterogenitas data melalui penambahan efek acak (random effects) sehingga memberikan estimasi parameter yang lebih stabil pada data longitudinal maupun data spasial.

Meskipun demikian, sebagian besar penelitian sebelumnya masih menganalisis aspek spasial atau temporal secara terpisah. Penelitian yang mengombinasikan analisis Incidence Rate (IR), Standardized Incidence Ratio (SIR), pembentukan matriks ketetanggaan K-Nearest Neighbor (KNN), pengujian Global Moran’s I dan Local Moran’s I, serta pemodelan GLMM Negative Binomial dalam satu kerangka analisis terpadu masih relatif terbatas, khususnya pada kasus stunting di Provinsi Nusa Tenggara Timur. Selain itu, sebagian besar penelitian terdahulu belum mengakomodasi karakteristik data yang mengalami overdispersi sekaligus mempertimbangkan heterogenitas spasial dan temporal melalui efek acak. Kondisi tersebut menunjukkan adanya research gap, yaitu belum optimalnya pendekatan analisis yang mampu menggambarkan pola penyebaran stunting secara komprehensif baik dari aspek ruang maupun waktu.

Berdasarkan kesenjangan tersebut, penelitian ini menawarkan pendekatan analisis spasio-temporal yang mengintegrasikan analisis deskriptif, perhitungan IR dan SIR, pembentukan matriks ketetanggaan spasio-temporal menggunakan KNN, pengujian autokorelasi spasial melalui Global Moran’s I dan Local Moran’s I (LISA), serta pemodelan GLMM Negative Binomial. Pendekatan ini diharapkan mampu mengidentifikasi wilayah dengan risiko stunting tinggi, mengungkap pola pengelompokan spasial, serta menghasilkan model statistik yang lebih sesuai dengan karakteristik data stunting yang bersifat spasio-temporal dan mengalami overdispersi.

Penelitian ini bertujuan untuk menganalisis pola spasio-temporal kasus stunting di Provinsi Nusa Tenggara Timur periode 2021-2025, mengidentifikasi adanya autokorelasi spasial dan klaster wilayah berisiko tinggi, serta membangun model GLMM Negative Binomial yang mampu menjelaskan variasi jumlah kasus stunting dengan mempertimbangkan efek spasial dan temporal. Secara teoretis, penelitian ini diharapkan dapat memperkaya pengembangan metode analisis spasio-temporal pada data epidemiologi berbentuk data cacah. Secara praktis, hasil penelitian diharapkan dapat menjadi dasar bagi pemerintah daerah dan pemangku kebijakan dalam menyusun strategi intervensi penanggulangan stunting yang lebih tepat sasaran, efektif, dan berbasis karakteristik wilayah.

DATA DAN METODE

Data dan Variabel

Penelitian ini menggunakan data sekunder yang bersifat agregat pada tingkat kabupaten/kota di Provinsi Nusa Tenggara Timur (NTT) selama periode tahun 2021–2025. Data yang digunakan terdiri atas jumlah kasus stunting dan jumlah populasi berisiko pada masing-masing kabupaten/kota untuk setiap tahun pengamatan. Selain itu, penelitian ini memanfaatkan data spasial berupa shapefile batas administrasi kabupaten/kota Provinsi NTT yang digunakan untuk analisis spasial dan visualisasi disease mapping.

Unit observasi dalam penelitian ini merupakan kombinasi antara dimensi wilayah dan waktu (spatio-temporal panel data), yaitu 21 kabupaten/kota yang diamati selama lima tahun berturut-turut. Dengan demikian, jumlah observasi yang dianalisis sebanyak 105 unit pengamatan. Kabupaten Malaka tidak disertakan dalam analisis karena keterbatasan ketersediaan data selama periode penelitian sehingga observasi difokuskan pada wilayah yang memiliki data lengkap.

Variabel utama yang digunakan dalam penelitian ini adalah jumlah kasus stunting sebagai variabel respons dan jumlah populasi berisiko sebagai dasar perhitungan ukuran epidemiologi. Berdasarkan kedua variabel tersebut dihitung beberapa indikator epidemiologi, yaitu Incidence Rate (IR), Expected Cases, dan Standardized Incidence Ratio (SIR). Selain itu, dibentuk variabel spatial lag yang merepresentasikan rata-rata jumlah kasus pada wilayah-wilayah tetangga berdasarkan matriks ketetanggaan spasial.

Analisis Deskriptif

Analisis deskriptif merupakan tahap awal dalam penelitian epidemiologi yang bertujuan menggambarkan distribusi penyakit berdasarkan dimensi waktu (time), tempat (place), dan populasi (person). Menurut Leon Gordis, epidemiologi deskriptif digunakan untuk mengidentifikasi pola kejadian penyakit serta menghasilkan hipotesis mengenai faktor-faktor yang mungkin berhubungan dengan kejadian tersebut. Dalam epidemiologi spasial, analisis deskriptif berfungsi untuk memberikan gambaran awal mengenai variasi kasus antarwilayah dan perubahan kasus dari waktu ke waktu sebelum dilakukan analisis spasial yang lebih kompleks. Analisis deskriptif pada penelitian ini digunakan untuk mengevaluasi karakteristik dasar kasus stunting di Provinsi NTT selama periode 2021–2025. Hasil analisis ini memberikan informasi awal mengenai pola distribusi kasus, tren temporal, dan variasi geografis yang menjadi dasar dalam analisis epidemiologi spasial berikutnya.

ukuran yang digunakan meliputi:

\[ \begin{aligned} \bar{Y} &= \frac{\sum_{i=1}^{n} Y_i}{n} \\ Y_{min} &= \min(Y_i) \\ Y_{max} &= \max(Y_i) \end{aligned} \]

dengan total kasus tahunan;

\[Y_t = \sum_{i=1}^{n} Y_{it}\]

Incidence Rate (IR)

Incidence Rate (IR) merupakan ukuran epidemiologi yang menggambarkan tingkat kejadian suatu penyakit pada populasi berisiko dalam periode tertentu. Menurut Kenneth Rothman, incidence rate mengukur laju terjadinya kasus baru dalam suatu populasi sehingga dapat digunakan untuk membandingkan tingkat risiko antarwilayah yang memiliki ukuran populasi berbeda. Kabupaten/kota di NTT memiliki jumlah penduduk yang berbeda sehingga penggunaan jumlah kasus absolut dapat menghasilkan interpretasi yang bias. Oleh karena itu incidence rate digunakan untuk memperoleh ukuran risiko yang telah mempertimbangkan ukuran populasi sehingga memungkinkan perbandingan risiko stunting antarwilayah secara lebih objektif.

secara matematis adalah sebagai berikut.

\[IR_i = \frac{Y_i}{N_i} \times 1000\]

Keterangan:

\(Y_i\) = jumlah kasus stunting wilayah ke-i

\(N_i\) = jumlah populasi berisiko wilayah ke-i

Standardized Incidence Ratio (SIR)

Standardized Incidence Ratio (SIR) merupakan ukuran risiko relatif yang membandingkan jumlah kasus yang diamati dengan jumlah kasus yang diharapkan berdasarkan tingkat kejadian rata-rata populasi referensi. Menurut Andrew Lawson, SIR merupakan indikator yang paling umum digunakan dalam disease mapping karena mampu menunjukkan wilayah yang memiliki risiko lebih tinggi atau lebih rendah dibandingkan wilayah referensi. SIR digunakan karena dapat menghilangkan pengaruh perbedaan ukuran populasi antarwilayah. Dalam penelitian ini SIR digunakan untuk mengidentifikasi kabupaten/kota yang memiliki risiko stunting relatif lebih tinggi dibandingkan rata-rata Provinsi NTT.

Tingkat Kejadian Keseluruhan:

\[r = \frac{\sum_{i} Y_i}{\sum_{i} N_i}\]

Kasus Harapan:

\[E_i = N_i \times r\]

SIR:

\[SIR_i = \frac{Y_i}{E_i}\]

Matriks Bobot Spasial

Menurut Luc Anselin, matriks bobot spasial merupakan representasi matematis dari hubungan geografis antar lokasi yang menjadi dasar dalam pengukuran autokorelasi spasial. Provinsi NTT memiliki karakteristik wilayah kepulauan sehingga pendekatan berbasis batas administrasi seperti Queen atau Rook Contiguity kurang optimal. Oleh karena itu metode KNN dipilih karena hubungan spasial ditentukan berdasarkan jarak antarwilayah sehingga lebih sesuai untuk menggambarkan interaksi geografis antar kabupaten/kota di NTT.

Bobot spasial KNN: \[w_{ij} = \begin{cases} 1, & j \in N_k(i) \\ 0, & \text{lainnya} \end{cases}\]

dengan \(k = 3\).

Karena penelitian ini melibatkan dimensi ruang dan waktu, maka dibentuk matriks bobot spasio-temporal menggunakan produk Kronecker antara matriks identitas waktu dan matriks bobot spasial.

Matriks spasio-temporal: \[W_{ST} = I_T \otimes W\] dengan:

\(I_T\) = matriks identitas waktu

\(W\) = matriks bobot spasial

Global Moran’s I

Autokorelasi spasial menunjukkan adanya ketergantungan nilai suatu variabel pada lokasi yang berdekatan. Konsep ini didasarkan pada Hukum Pertama Geografi Tobler yang menyatakan bahwa lokasi yang berdekatan cenderung memiliki karakteristik yang lebih mirip dibandingkan lokasi yang berjauhan. Global Moran’s I merupakan statistik yang paling umum digunakan untuk mengukur autokorelasi spasial secara keseluruhan (global) dalam suatu wilayah studi (Anselin, 1995).

Global Moran’s I digunakan untuk menguji apakah distribusi kasus stunting menunjukkan pola pengelompokan spasial yang signifikan selama periode pengamatan. Hasil pengujian ini menjadi dasar untuk menentukan apakah analisis spasial lanjutan diperlukan.

secara matematis adalah sebagai berikut.

\[I = \frac{n}{S_0} \frac{\sum_{i} \sum_{j} w_{ij}(y_i - \bar{y})(y_j - \bar{y})}{\sum_{i} (y_i - \bar{y})^2}\] dengan: \[S_0 = \sum_{i} \sum_{j} w_{ij}\]

Local Moran’s I (LISA)

Global Moran’s I hanya memberikan informasi mengenai keberadaan autokorelasi secara keseluruhan tanpa menunjukkan lokasi spesifik yang membentuk pola tersebut. Oleh karena itu digunakan Local Indicators of Spatial Association (LISA) yang diperkenalkan oleh Anselin (1995). Metode ini memungkinkan identifikasi hotspot, coldspot, dan pencilan spasial (spatial outlier) kasus stunting pada setiap wilayah kabupaten/kota di NTT.

secara matematis adalah sebagai berikut.

\[I_i = z_i \sum_{j} w_{ij} z_j\] dengan: \[z_i = \frac{y_i - \bar{y}}{s}\]

Uji Overdispersi

Regresi Poisson merupakan model standar untuk data cacahan yang mengasumsikan bahwa nilai rata-rata sama dengan varians. Menurut Joseph Hilbe, data epidemiologi sering menunjukkan varians yang lebih besar daripada rata-ratanya sehingga terjadi kondisi yang disebut overdispersi. Pengujian overdispersi dilakukan untuk memastikan kesesuaian penggunaan model Poisson. Apabila ditemukan overdispersi maka model Negative Binomial lebih tepat digunakan karena mampu mengakomodasi variabilitas data yang lebih besar.

Asumsi Poisson: \[E(Y) = Var(Y)\]

Statistik overdispersi: \[\phi = \frac{\text{Deviance}}{\text{df Residual}}\]

Kriteria: \[\phi > 1\]

menunjukkan adanya overdispersi.

Negative Binomial

Regresi Negative Binomial merupakan pengembangan dari regresi Poisson yang dirancang untuk menangani data cacahan yang mengalami overdispersi. Menurut Hilbe (2011), model ini menambahkan parameter dispersi sehingga varians dapat melebihi nilai rata-rata. Model Negative Binomial digunakan karena data stunting berupa data cacahan dan hasil uji menunjukkan adanya overdispersi. Selain itu model ini mampu mengevaluasi pengaruh faktor spasial melalui variabel spatial lag serta mempertimbangkan perbedaan risiko antarwilayah dan antarperiode.

secara matematis adalah sebagai berikut.

\[Y_{it} \sim NB(\mu_{it}, \theta)\]

dengan: \[Var(Y_{it}) = \mu_{it} + \frac{\mu_{it}^2}{\theta}\]

sehingga, \[\log(\mu_{it}) = \beta_0 + \beta_1 LagKasus_{it} + \alpha_i + \gamma_t + \log(Pop_{it})\]

Generalized Linear Mixed Model (GLMM) Negative Binomial

GLMM merupakan perluasan dari Generalized Linear Model yang menggabungkan efek tetap (fixed effects) dan efek acak (random effects). Menurut Alain Zuur, GLMM sangat sesuai untuk data longitudinal dan hierarkis karena mampu menangkap heterogenitas yang tidak dapat dijelaskan oleh kovariat yang diamati. Data penelitian memiliki struktur panel spasio-temporal berupa pengamatan berulang pada kabupaten/kota selama lima tahun. Oleh karena itu digunakan GLMM Negative Binomial untuk mengakomodasi heterogenitas antarwilayah dan antarwaktu yang tidak dapat dijelaskan oleh variabel yang diamati.

secara matematis adalah sebagai berikut.

\[Y_{it} \sim NB(\mu_{it}, \theta)\] \[\log(\mu_{it}) = \beta_0 + \beta_1 LagKasus_{it} + u_i + v_t + \log(Pop_{it})\]

dengan: \[u_i \sim N(0, \sigma_u^2)\] \[v_t \sim N(0, \sigma_v^2)\]

Akaike Information Criterion (AIC)

Ukuran yang digunakan untuk membandingkan beberapa model statistik berdasarkan keseimbangan antara tingkat kecocokan model dan kompleksitas model. Menurut Hirotugu Akaike, model terbaik adalah model yang mampu menjelaskan data dengan baik menggunakan jumlah parameter yang seminimal mungkin. AIC digunakan untuk membandingkan model Negative Binomial dan GLMM Negative Binomial. Model dengan nilai AIC terkecil dipilih sebagai model terbaik karena dianggap memiliki keseimbangan paling optimal antara ketepatan prediksi dan kompleksitas model.

secara matematis adalah sebagai berikut.

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

dengan

\(L\) = likelihood model

\(k\) = jumlah parameter model

HASIL DAN PEMBAHASAN

Statistik Deskriptif Data: Analisis statistik deskriptif dilakukan untuk memberikan gambaran umum mengenai karakteristik data jumlah kasus stunting pada 21 kabupaten/kota di Provinsi Nusa Tenggara Timur selama periode tahun 2021–2025. Statistik deskriptif meliputi ukuran pemusatan, ukuran penyebaran, serta bentuk distribusi data. Hasil analisis disajikan pada tabel berikut.

Berdasarkan Tabel 1, jumlah observasi yang dianalisis sebanyak 105, yang merupakan kombinasi data dari 21 kabupaten/kota selama lima tahun pengamatan (2021-2025). Rata-rata jumlah kasus stunting sebesar 3.230,53 kasus dengan simpangan baku sebesar 2.452,58 kasus. Nilai simpangan baku yang relatif besar, yaitu sekitar 75,9% dari nilai rata-rata yang menunjukkan bahwa jumlah kasus stunting memiliki variasi yang cukup tinggi antarwilayah maupun antarperiode pengamatan. Hal ini mengindikasikan bahwa distribusi kasus stunting di Provinsi NTT bersifat heterogen dan tidak tersebar secara merata.

Nilai minimum jumlah kasus stunting tercatat sebesar 549 kasus, sedangkan nilai maksimum mencapai 13.123 kasus, sehingga diperoleh rentang data (range) sebesar 12.574 kasus. Perbedaan yang cukup besar antara nilai minimum dan maksimum menunjukkan adanya kesenjangan jumlah kasus yang signifikan antar kabupaten/kota. Dengan kata lain, terdapat wilayah yang memiliki jumlah kasus stunting jauh lebih tinggi dibandingkan wilayah lainnya. Kondisi ini mengindikasikan adanya variasi spasial yang kuat sehingga diperlukan analisis lebih lanjut untuk mengetahui apakah perbedaan tersebut membentuk pola pengelompokan secara geografis.

Selain ukuran pemusatan dan penyebaran, bentuk distribusi data juga dapat dilihat melalui nilai skewness dan kurtosis. Nilai skewness sebesar 1,93 yang menunjukkan bahwa distribusi jumlah kasus stunting menceng ke kanan (positively skewed). Artinya, sebagian besar kabupaten/kota memiliki jumlah kasus yang relatif rendah hingga sedang, sedangkan hanya beberapa wilayah yang memiliki jumlah kasus sangat tinggi sehingga menarik distribusi ke arah kanan. Pola ini umum dijumpai pada data epidemiologi, yaitu di mana hanya sedikit wilayah yang menjadi pusat konsentrasi kasus.

Selanjutnya, nilai kurtosis sebesar 4,18 mengindikasikan bahwa distribusi data bersifat leptokurtik, yaitu memiliki puncak distribusi yang lebih tinggi dibandingkan distribusi normal serta ekor distribusi yang lebih panjang. Kondisi tersebut menunjukkan adanya beberapa observasi dengan nilai yang jauh lebih besar daripada sebagian besar data lainnya (extreme values). Keberadaan nilai ekstrem ini memperkuat dugaan bahwa distribusi kasus stunting antarwilayah tidak homogen dan berpotensi dipengaruhi oleh karakteristik wilayah tertentu.

Secara keseluruhan, hasil statistik deskriptif menunjukkan bahwa data kasus stunting di Provinsi Nusa Tenggara Timur memiliki tingkat variasi yang tinggi, distribusi yang tidak simetris, serta adanya nilai-nilai ekstrem. Karakteristik tersebut mengindikasikan bahwa data memiliki heterogenitas yang cukup besar sehingga pendekatan analisis yang hanya mengandalkan statistik deskriptif belum mampu menjelaskan pola hubungan antarwilayah. Oleh karena itu, penelitian ini dilanjutkan dengan analisis epidemiologi spasial melalui perhitungan Incidence Rate (IR), Standardized Incidence Ratio (SIR), analisis Global Moran’s I, Local Moran’s I (LISA), serta pemodelan Negatife Binomial dan Generalized Linear Mixed Model (GLMM) untuk mengevaluasi pengaruh aspek spasial dan temporal terhadap distribusi kasus stunting di Provinsi Nusa Tenggara Timur.

Berdasarkan Gambar, terlihat bahwa distribusi kasus stunting di Provinsi Nusa Tenggara Timur (NTT) selama periode 2021–2025 menunjukkan pola spasial yang tidak merata antar kabupaten/kota. Perbedaan intensitas warna pada peta menunjukkan adanya variasi jumlah kasus stunting yang cukup besar di setiap wilayah. Kabupaten Timor Tengah Selatan dan Sumba Barat Daya secara konsisten memiliki jumlah kasus yang relatif tinggi dibandingkan kabupaten lainnya, yang ditunjukkan dengan warna yang lebih terang pada sebagian besar periode pengamatan. Sebaliknya, beberapa kabupaten di wilayah Flores bagian tengah dan utara cenderung memiliki jumlah kasus yang lebih rendah.

Jika diamati dari tahun ke tahun, pola persebaran kasus tidak mengalami perubahan yang signifikan. Wilayah yang memiliki jumlah kasus tinggi pada awal periode pengamatan tetap menjadi wilayah dengan jumlah kasus tinggi pada tahun-tahun berikutnya. Kondisi tersebut menunjukkan adanya persistensi spasial, yaitu kecenderungan suatu wilayah mempertahankan tingkat kejadian stunting yang relatif sama sepanjang waktu. Hal ini mengindikasikan bahwa faktor-faktor lokal, seperti kondisi sosial ekonomi, kepadatan penduduk, akses terhadap pelayanan kesehatan, sanitasi, dan karakteristik geografis, masih berperan dalam membentuk pola penyebaran kasus stunting di Provinsi NTT.

Secara umum, peta sebaran ini memberikan indikasi awal bahwa kasus stunting memiliki karakteristik spasial sehingga diperlukan analisis autokorelasi spasial pada tahap berikutnya untuk menguji apakah pengelompokan wilayah tersebut signifikan secara statistik.

Berdasarkan Gambar, perkembangan total kasus stunting di Provinsi Nusa Tenggara Timur selama periode 2021-2025 menunjukkan pola yang fluktuatif dengan kecenderungan menurun pada awal periode pengamatan. Total kasus stunting tercatat sebesar 78.124 kasus pada tahun 2021, kemudian menurun menjadi 74.756 kasus pada tahun 2022. Penurunan yang lebih besar terjadi pada tahun 2023, yaitu menjadi 61.395 kasus, atau sekitar 17,9% dibandingkan tahun sebelumnya. Kondisi ini mengindikasikan adanya perbaikan dalam upaya penanganan stunting melalui berbagai program intervensi yang telah dilaksanakan.

Namun demikian, setelah mencapai titik terendah pada tahun 2023, jumlah kasus kembali mengalami peningkatan meskipun relatif kecil. Total kasus meningkat menjadi 62.218 kasus pada tahun 2024 dan kembali naik menjadi 62.713 kasus pada tahun 2025. Peningkatan tersebut menunjukkan bahwa meskipun program penurunan stunting telah memberikan dampak positif, keberhasilan yang dicapai belum sepenuhnya stabil sehingga masih diperlukan upaya pengendalian yang berkelanjutan.

Secara keseluruhan, tren temporal menunjukkan bahwa kasus stunting di Provinsi NTT mengalami penurunan yang cukup signifikan selama periode 2021-2023, kemudian mengalami kenaikan ringan pada periode 2024–2025. Pola tersebut mengindikasikan bahwa dinamika kasus stunting dipengaruhi oleh faktor waktu, sehingga analisis yang mengintegrasikan dimensi waktu dan spasial menjadi penting untuk memahami perubahan risiko stunting secara lebih komprehensif.

Standardized Incidence Ratio (SIR): SIR digunakan untuk menggambarkan tingkat risiko relatif kejadian stunting pada setiap kabupaten/kota setelah memperhitungkan ukuran populasi berisiko. Berbeda dengan jumlah kasus absolut, nilai SIR menunjukkan apakah suatu wilayah memiliki jumlah kasus yang lebih tinggi atau lebih rendah dibandingkan jumlah kasus yang diharapkan berdasarkan rata-rata risiko di seluruh Provinsi Nusa Tenggara Timur. Nilai SIR > 1 menunjukkan bahwa jumlah kasus yang diamati lebih tinggi daripada jumlah kasus yang diharapkan, sedangkan SIR < 1 menunjukkan jumlah kasus yang lebih rendah daripada yang diharapkan.

Hasil pemetaan SIR menunjukkan bahwa distribusi risiko stunting di Provinsi Nusa Tenggara Timur bersifat heterogen. Walaupun secara umum jumlah kasus stunting mengalami penurunan selama periode penelitian, tidak semua wilayah mengalami penurunan risiko yang sama. Beberapa kabupaten tetap memiliki nilai SIR di atas satu secara konsisten, yang mengindikasikan bahwa jumlah kasus yang terjadi masih lebih tinggi dibandingkan jumlah kasus yang diharapkan berdasarkan rata-rata provinsi.

Temuan ini menunjukkan pentingnya penggunaan SIR dibandingkan hanya mengandalkan jumlah kasus absolut. Kabupaten dengan jumlah kasus yang besar belum tentu memiliki risiko relatif yang tinggi apabila wilayah tersebut juga memiliki populasi berisiko yang besar. Sebaliknya, kabupaten dengan jumlah kasus yang tidak terlalu besar dapat memiliki nilai SIR tinggi apabila jumlah kasus yang terjadi jauh melebihi jumlah kasus yang diharapkan berdasarkan ukuran populasinya. Oleh karena itu, SIR mampu memberikan gambaran risiko yang lebih objektif karena telah melakukan standardisasi terhadap ukuran populasi (Lawson, 2018).

Selain itu, pola konsentrasi nilai SIR yang relatif tetap pada beberapa wilayah selama lima tahun pengamatan mengindikasikan adanya faktor-faktor lokal yang memengaruhi kejadian stunting. Faktor tersebut dapat berupa kondisi sosial ekonomi masyarakat, tingkat kemiskinan, akses terhadap pelayanan kesehatan, kualitas sanitasi, kondisi geografis, maupun karakteristik lingkungan yang berbeda antarwilayah. Menurut Elliott et al. (2000), pola risiko yang konsisten pada lokasi tertentu merupakan salah satu indikasi adanya variasi spasial yang tidak bersifat acak, sehingga memerlukan analisis autokorelasi spasial untuk mengetahui apakah wilayah dengan risiko tinggi cenderung berdekatan secara geografis.

Oleh karena itu, hasil pemetaan SIR pada penelitian ini menjadi dasar untuk melanjutkan analisis menggunakan Global Moran’s I dan Local Moran’s I (LISA). Kedua metode tersebut digunakan untuk menguji apakah pola risiko relatif stunting yang terlihat pada peta benar-benar membentuk klaster spasial yang signifikan secara statistik, sehingga dapat memberikan informasi yang lebih kuat dalam penentuan wilayah prioritas intervensi kesehatan di Provinsi Nusa Tenggara Timur.

Analisis Autokorelasi Spatio-Temporal: Analisis autokorelasi spasio-temporal dilakukan untuk mengetahui apakah distribusi kasus stunting di Provinsi Nusa Tenggara Timur (NTT) memiliki pola pengelompokan (cluster) yang signifikan secara spasial selama periode 2021–2025. Pengujian dilakukan menggunakan Global Moran’s I untuk mengidentifikasi autokorelasi spasial secara keseluruhan, kemudian dilanjutkan dengan Local Moran’s I (LISA) untuk mengidentifikasi lokasi-lokasi yang menjadi pusat klaster risiko tinggi maupun rendah.

Berdasarkan Tabel 3 diperoleh nilai Global Moran’s I sebesar 0,189 dengan nilai harapan (Expected Index) sebesar -0,010 dan p-value sebesar 0,003. Nilai Moran’s I yang bernilai positif menunjukkan bahwa distribusi kasus stunting di Provinsi Nusa Tenggara Timur memiliki kecenderungan membentuk pola pengelompokan spasial (positive spatial autocorrelation). Selain itu, nilai p-value yang lebih kecil dari 0,05 menunjukkan bahwa autokorelasi spasial tersebut signifikan secara statistik, sehingga hipotesis nol yang menyatakan tidak terdapat autokorelasi spasial ditolak.

Nilai Moran’s I sebesar 0,189 menunjukkan tingkat autokorelasi yang positif tetapi relatif lemah hingga sedang. Dengan demikian, kabupaten/kota yang memiliki jumlah kasus stunting tinggi cenderung bertetangga dengan wilayah yang juga memiliki jumlah kasus relatif tinggi, sedangkan wilayah dengan jumlah kasus rendah cenderung berada di sekitar wilayah dengan jumlah kasus rendah. Meskipun demikian, besarnya koefisien yang masih jauh dari nilai satu menunjukkan bahwa pengelompokan tersebut tidak terjadi secara sempurna pada seluruh wilayah penelitian.

Hasil ini mengindikasikan bahwa distribusi kasus stunting di Provinsi NTT tidak terjadi secara acak (random), melainkan dipengaruhi oleh keterkaitan antarwilayah yang berdekatan. Kondisi tersebut menguatkan dugaan bahwa faktor geografis, kondisi sosial ekonomi, akses pelayanan kesehatan, maupun karakteristik lingkungan yang saling berdekatan berkontribusi terhadap pola penyebaran kasus stunting.

Garis regresi pada Moran Scatter Plot memiliki kemiringan positif yang menunjukkan adanya hubungan positif antara jumlah kasus stunting pada suatu kabupaten/kota dengan rata-rata jumlah kasus pada wilayah tetangganya. Semakin tinggi jumlah kasus pada suatu wilayah, semakin tinggi pula kecenderungan wilayah di sekitarnya memiliki jumlah kasus yang relatif tinggi.

Sebagian besar titik pengamatan berada di sekitar pusat diagram, yang menunjukkan bahwa mayoritas kabupaten/kota memiliki jumlah kasus yang mendekati rata-rata provinsi. Namun demikian, terdapat beberapa observasi yang berada jauh dari pusat diagram, seperti observasi 20, 41, 59, 62, 83, dan 104, yang menunjukkan adanya wilayah dengan karakteristik spasial yang berbeda dibandingkan wilayah lainnya.

Keberadaan titik-titik tersebut menunjukkan bahwa tidak seluruh wilayah mengikuti pola spasial yang sama. Beberapa wilayah memiliki jumlah kasus yang jauh lebih tinggi atau lebih rendah dibandingkan tetangganya sehingga memerlukan analisis lokal untuk mengetahui tipe pengelompokan yang terbentuk.

Secara keseluruhan, Moran Scatter Plot memperkuat hasil uji statistik Global Moran’s I bahwa distribusi kasus stunting di Provinsi NTT menunjukkan autokorelasi spasial positif yang signifikan, meskipun tingkat pengelompokannya masih tergolong sedang.

Berdasarkan Gambar di atas terlihat bahwa sebagian besar kabupaten/kota di Provinsi Nusa Tenggara Timur termasuk dalam kategori Not Significant yang ditunjukkan dengan warna abu-abu. Hal ini menunjukkan bahwa sebagian besar wilayah tidak memiliki autokorelasi lokal yang signifikan selama periode penelitian.

Meskipun demikian, terdapat beberapa kabupaten yang secara konsisten membentuk klaster High–High, yaitu wilayah dengan jumlah kasus stunting tinggi yang dikelilingi oleh wilayah-wilayah dengan jumlah kasus tinggi pula. Klaster ini tampak terutama berada di bagian timur Pulau Flores pada tahun 2021, 2022, 2024, dan 2025. Keberadaan klaster tersebut mengindikasikan bahwa wilayah-wilayah tersebut merupakan hotspot stunting yang memiliki risiko relatif lebih tinggi dibandingkan wilayah lainnya.

Selain klaster High–High, pada tahun 2021 dan 2022 juga terlihat satu wilayah yang termasuk dalam kategori Low–High. Klaster ini menunjukkan adanya kabupaten dengan jumlah kasus relatif rendah, tetapi dikelilingi oleh wilayah-wilayah yang memiliki jumlah kasus tinggi. Kondisi tersebut mengindikasikan adanya perbedaan karakteristik lokal yang menyebabkan wilayah tersebut memiliki tingkat kejadian stunting yang lebih rendah dibandingkan lingkungan sekitarnya.

Menariknya, pada tahun 2023 hampir seluruh wilayah termasuk kategori Not Significant yang menunjukkan bahwa pola pengelompokan lokal pada tahun tersebut relatif melemah. Kondisi ini sejalan dengan hasil analisis deskriptif sebelumnya yang menunjukkan bahwa tahun 2023 merupakan periode dengan jumlah kasus stunting terendah selama penelitian. Penurunan jumlah kasus tersebut tampaknya diikuti dengan berkurangnya konsentrasi spasial sehingga klaster signifikan hampir tidak terbentuk.

Pada tahun 2024 dan 2025, klaster High–High kembali muncul pada wilayah yang hampir sama dengan tahun-tahun sebelumnya. Hal ini menunjukkan bahwa meskipun terjadi penurunan jumlah kasus secara keseluruhan, beberapa wilayah tetap mempertahankan tingkat risiko yang relatif tinggi dari waktu ke waktu. Pola tersebut mengindikasikan adanya persistensi spasio-temporal, yaitu kecenderungan suatu wilayah untuk tetap menjadi daerah berisiko tinggi selama beberapa periode pengamatan berturut-turut.

Uji Overdispersi Model Poisson: Model Poisson digunakan untuk mengetahui apakah data jumlah kasus stunting memenuhi asumsi equidispersion, yaitu kondisi ketika nilai rata-rata sama dengan nilai varians. Apabila asumsi tersebut tidak terpenuhi, maka model Poisson menjadi kurang sesuai sehingga diperlukan model alternatif yang mampu mengakomodasi varians yang lebih besar daripada rata-rata. Model Poisson pada penelitian ini menggunakan jumlah kasus stunting sebagai variabel respon, faktor tahun sebagai variabel penjelas, serta logaritma jumlah populasi sebagai offset, sehingga model mengestimasi laju kejadian (incidence rate) pada masing-masing kabupaten/kota.

Berdasarkan Tabel 5 diperoleh nilai Null Deviance sebesar 57.522, sedangkan Residual Deviance sebesar 53.186. Penurunan nilai deviance tersebut menunjukkan bahwa penambahan variabel tahun mampu meningkatkan kemampuan model dalam menjelaskan variasi jumlah kasus stunting dibandingkan model tanpa prediktor. Namun demikian, perhatian utama pada penelitian ini adalah nilai Dispersion Parameter (φ) yang diperoleh sebesar 532. Nilai tersebut jauh lebih besar dibandingkan satu (φ >> 1), sehingga menunjukkan adanya overdispersi yang sangat kuat pada data kasus stunting.

Secara teoritis, model Poisson mengasumsikan bahwa nilai rata-rata sama dengan nilai varians. Ketika varians jauh lebih besar daripada rata-rata, maka asumsi tersebut tidak lagi terpenuhi. Akibatnya, simpangan baku parameter menjadi terestimasi lebih kecil daripada kondisi sebenarnya (underestimated standard error), sehingga uji signifikansi parameter pada Tabel 4 berpotensi menghasilkan kesimpulan yang terlalu optimistis (inflated Type I error).

Oleh karena itu, model Poisson dinilai kurang sesuai untuk menganalisis data kasus stunting pada penelitian ini. Sebagai alternatif, digunakan negatif binomial yang menambahkan parameter dispersi sehingga mampu mengakomodasi varians yang lebih besar daripada rata-rata Cameron & Trivedi (2013). Selanjutnya, untuk mengakomodasi ketergantungan spasial dan temporal yang masih tersisa antar kabupaten dan antar tahun, analisis dilanjutkan menggunakan GLMM Negatif Binomial dengan efek acak spasial dan temporal.

Pemodelan

Pada penelitian ini dibandingkan dua model, yaitu:

Model 1 : Regresi Binomial Negatif dengan efek tetap (Fixed Effect).
Model 2 : Generalized Linear Mixed Model (GLMM) Negative Binomial dengan efek acak kabupaten/kota dan tahun.

Model 1

Model binomial negatif yang dibangun pada penelitian ini menggunakan fungsi link log dengan offset log(populasi), sehingga model yang diperoleh adalah sebagai berikut.

\[ \begin{aligned}\log(\mu_i) = {} & \log(populasi_i) - 1.948 - 0.072(Belu) - 0.361(Ende) + 0.326(Flores\ Timur) \\& + 0.108(Kota\ Kupang) - 0.103(Kupang) - 0.092(Lembata) + 0.039(Manggarai) \\& - 0.138(Manggarai\ Barat) - 0.355(Manggarai\ Timur) - 0.360(Nagekeo) - 0.294(Ngada) \\& + 0.286(Rote\ Ndao) + 0.284(Sabu\ Raijua) + 0.133(Sikka) + 0.292(Sumba\ Barat) \\& + 0.866(Sumba\ Barat\ Daya) - 0.362(Sumba\ Tengah) + 0.027(Sumba\ Timur) \\& + 0.606(Timor\ Tengah\ Selatan) + 0.360(Timor\ Tengah\ Utara) - 0.185(Tahun2022) \\& - 0.341(Tahun2023) - 0.153(Tahun2024) - 0.028(Tahun2025) + 0.000(LagKasus)\end{aligned} \]

Berdasarkan Tabel 6, sebagian besar perbedaan jumlah kasus stunting dipengaruhi oleh karakteristik kabupaten/kota. Kabupaten Sumba Barat Daya memiliki koefisien terbesar (β = 0,866; p < 0,001), diikuti Timor Tengah Selatan (β = 0,606; p < 0,001), Flores Timur (β = 0,326; p = 0,0077), Sumba Barat (β = 0,292; p = 0,016), Rote Ndao (β = 0,286; p = 0,014), dan Sabu Raijua (β = 0,284; p = 0,014). Nilai koefisien positif menunjukkan bahwa wilayah-wilayah tersebut memiliki laju kejadian stunting yang lebih tinggi dibandingkan wilayah referensi setelah memperhitungkan ukuran populasi. Sebaliknya, Kabupaten Ende, Manggarai Timur, Nagekeo, Ngada, dan Sumba Tengah memiliki koefisien negatif yang signifikan (p < 0,05). Hal tersebut menunjukkan bahwa laju kejadian stunting pada kabupaten tersebut relatif lebih rendah dibandingkan wilayah referensi.

Variabel waktu juga menunjukkan adanya perubahan kejadian stunting. Tahun 2022, 2023, dan 2024 memiliki koefisien negatif yang signifikan, ini mengindikasikan adanya penurunan jumlah kasus dibandingkan tahun 2021. Penurunan terbesar terjadi pada tahun 2023 (β = −0,341), sejalan dengan hasil analisis deskriptif yang menunjukkan penurunan total kasus stunting pada tahun tersebut. Sementara itu, tahun 2025 tidak berbeda secara signifikan dengan tahun 2021 (p = 0,663), sehingga mengindikasikan bahwa jumlah kasus mulai kembali meningkat.

Menariknya, variabel lag_kasus memiliki koefisien positif tetapi tidak signifikan (p = 0,201). Hal ini menunjukkan bahwa hubungan spasial antarwilayah belum dapat dijelaskan secara memadai apabila seluruh kabupaten diasumsikan memiliki karakteristik yang sama (fixed effect).

Model 2

Untuk mengakomodasi heterogenitas antarwilayah dan antarperiode, dilakukan pengembangan model menggunakan Generalized Linear Mixed Model (GLMM) Negative Binomial dengan efek acak kabupaten/kota dan tahun. Model yang diperoleh adalah sebagai berikut.

\[ \log(\mu_{ij}) = \log(\textit{populasi}_{ij}) - 2.108 + 0.000(\textit{LagKasus}_{ij}) + u_i + v_j \]

dengan

\(\begin{aligned} u_i & \text{ merupakan efek acak kabupaten/kota,} \\ v_j & \text{ merupakan efek acak tahun.} \end{aligned}\)

Model ini memungkinkan setiap kabupaten memiliki tingkat risiko dasar (baseline risk) yang berbeda serta mengakomodasi variasi antarperiode pengamatan. Berdasarkan Tabel 6 diperoleh bahwa intersep model bernilai −2,108 dengan nilai p < 0,001 yang menunjukkan bahwa model secara keseluruhan signifikan. Berbeda dengan model sebelumnya, variabel lag_kasus kini memiliki pengaruh yang signifikan terhadap jumlah kasus stunting (p = 0,0319). Walaupun koefisiennya sangat kecil akibat skala data yang besar, arah koefisien yang positif menunjukkan bahwa peningkatan jumlah kasus pada wilayah-wilayah yang bertetangga akan meningkatkan risiko kejadian stunting pada suatu kabupaten.

Hasil Perbandingan

Berdasarkan Tabel 7, model Negative Binomial Global memiliki 27 parameter, nilai AIC sebesar 1.641,02, BIC sebesar 1.712,67, dan Log-Likelihood sebesar -793,51. Sementara itu, model GLMM Negative Binomial (Mixed Effect) hanya menggunakan 5 parameter, dengan nilai AIC sebesar 1.687,10, BIC sebesar 1.700,37, dan Log-Likelihood sebesar -838,55.

Apabila ditinjau berdasarkan nilai AIC, model Negative Binomial memiliki nilai yang lebih kecil dibandingkan model GLMM Negative Binomial, yaitu 1.641,02 < 1.687,10. Menurut Akaike (1974), model dengan nilai AIC yang lebih rendah memiliki keseimbangan yang lebih baik antara tingkat kecocokan model (goodness of fit) dan kompleksitas model. Hasil ini menunjukkan bahwa berdasarkan kriteria AIC, model Negative Binomial memberikan kecocokan yang lebih baik terhadap data.

Apabila ditinjau dari nilai BIC, model GLMM Negative Binomial memiliki nilai yang lebih rendah dibandingkan model Negative Binomial , yaitu 1.700,37 < 1.712,67. Berbeda dengan AIC, BIC memberikan penalti yang lebih besar terhadap jumlah parameter sehingga lebih mengutamakan model yang lebih sederhana (parsimonious). Hal ini menunjukkan bahwa GLMM Negative Binomial mampu menghasilkan model yang relatif lebih efisien karena hanya menggunakan lima parameter tetapi masih mampu menjelaskan variasi data secara memadai.

Perbedaan tersebut juga terlihat dari jumlah parameter yang digunakan. Model Negative Binomial memasukkan seluruh efek tetap kabupaten/kota dan tahun sehingga jumlah parameter yang diestimasi menjadi lebih banyak. Sebaliknya, model GLMM hanya mengestimasi efek tetap untuk variabel lag_kasus, sedangkan variasi antar kabupaten/kota dan antar tahun dimodelkan sebagai efek acak (random effects). Pendekatan ini menghasilkan model yang lebih ringkas tanpa harus mengestimasi koefisien untuk setiap wilayah secara terpisah.

Dilihat dari nilai Log-Likelihood, model Negative Binomial menghasilkan nilai yang lebih besar (lebih mendekati nol), yaitu -793,51, dibandingkan GLMM sebesar -838,55. Nilai tersebut menunjukkan bahwa secara matematis model Negative Binomial memiliki kecocokan yang sedikit lebih baik terhadap data yang diamati. Akan tetapi, ukuran ini belum mempertimbangkan kompleksitas model sehingga perlu dievaluasi bersama AIC dan BIC.

Secara keseluruhan, hasil evaluasi menunjukkan bahwa kedua model memiliki keunggulan yang berbeda. Negative Binomial Global memberikan kecocokan model (fit) yang lebih tinggi berdasarkan AIC dan Log-Likelihood, sedangkan GLMM Negative Binomial menghasilkan model yang lebih sederhana berdasarkan BIC serta lebih mampu mengakomodasi heterogenitas spasial dan temporal melalui efek acak. Selain itu, pada model GLMM variabel lag_kasus terbukti signifikan (p = 0,0319), sedangkan pada model Negative Binomial Global variabel tersebut tidak signifikan (p = 0,201). Temuan ini menunjukkan bahwa setelah variasi antarwilayah dan antarwaktu dimodelkan sebagai efek acak, hubungan spasial antarwilayah menjadi lebih jelas.

Dengan demikian, pemilihan model terbaik bergantung pada tujuan analisis. Jika tujuan utama adalah memperoleh model dengan tingkat kecocokan statistik tertinggi, maka Negative Binomial lebih unggul berdasarkan nilai AIC dan Log-Likelihood. Namun, apabila tujuan penelitian adalah memperoleh model yang mampu merepresentasikan struktur data spatio-temporal secara lebih realistis dengan mempertimbangkan heterogenitas antar kabupaten/kota dan antar tahun, maka GLMM Negative Binomial lebih tepat digunakan, meskipun nilai AIC sedikit lebih besar. Dalam konteks penelitian epidemiologi spasial, kemampuan GLMM dalam mengakomodasi efek acak dan mengidentifikasi pengaruh spasial yang signifikan menjadikan model ini lebih representatif untuk menjelaskan pola penyebaran kasus stunting di Provinsi Nusa Tenggara Timur periode 2021-2025.

Incidence Rate Rasio

Hasil menunjukkan bahwa Kabupaten Sumba Barat Daya memiliki risiko tertinggi dengan IRR sekitar 2,38, berarti laju kejadian stunting sekitar 2,38 kali lebih tinggi dibandingkan wilayah referensi. Kabupaten Timor Tengah Selatan (IRR ≈ 1,83), Timor Tengah Utara (IRR ≈ 1,43), Flores Timur (IRR ≈ 1,39), serta Sumba Barat, Rote Ndao, dan Sabu Raijua (IRR ≈ 1,34) juga memiliki risiko yang lebih tinggi.

Sebaliknya, Ende, Nagekeo, Manggarai Timur, Ngada, dan Sumba Tengah memiliki IRR < 1 yang menunjukkan laju kejadian stunting lebih rendah dibandingkan wilayah referensi. Variabel waktu menunjukkan bahwa tahun 2022-2024 mengalami penurunan laju kejadian dibandingkan tahun 2021 dengan penurunan terbesar pada tahun 2023 (IRR ≈ 0,71), sedangkan tahun 2025 tidak berbeda signifikan dari tahun dasar.

Sementara itu, variabel lag_kasus memiliki nilai IRR yang mendekati 1, sehingga pengaruh wilayah tetangga belum signifikan pada model Negative Binomial. Secara keseluruhan, hasil ini menunjukkan bahwa variasi risiko stunting lebih dipengaruhi oleh karakteristik masing-masing kabupaten dengan Sumba Barat Daya dan Timor Tengah Selatan sebagai wilayah prioritas dalam upaya penanggulangan stunting.

Efek Acak Spasial

Gambar di atas menunjukkan distribusi efek acak spasial yang diestimasi menggunakan model GLMM Negative Binomial. Efek acak merepresentasikan variasi risiko kasus stunting antar kabupaten/kota yang tidak dapat dijelaskan oleh variabel tetap (fixed effects) dalam model. Nilai efek acak yang positif menunjukkan adanya kecenderungan jumlah kasus stunting yang lebih tinggi dari yang diperkirakan model, sedangkan nilai negatif menunjukkan kecenderungan jumlah kasus yang lebih rendah.

Secara umum, sebagian besar kabupaten/kota di Provinsi NTT memiliki efek acak yang relatif kecil, namun terdapat beberapa wilayah dengan nilai efek acak positif yang cukup tinggi, seperti Belu, Timor Tengah Selatan, dan Manggarai Timur yang mengindikasikan masih adanya faktor lokal yang meningkatkan risiko stunting di luar variabel yang dianalisis. Sebaliknya, beberapa wilayah seperti Sumba Barat Daya, Flores Timur, dan Timor Tengah Utara menunjukkan efek acak negatif, sehingga risiko kasus stunting di wilayah tersebut relatif lebih rendah setelah memperhitungkan pengaruh variabel dalam model.

Temuan ini menunjukkan bahwa masih terdapat heterogenitas spasial antar kabupaten/kota di Provinsi NTT. Oleh karena itu, selain mempertimbangkan faktor-faktor yang telah dimasukkan ke dalam model, diperlukan kebijakan intervensi yang disesuaikan dengan karakteristik lokal masing-masing wilayah untuk menurunkan kejadian stunting secara lebih efektif.

DAFTAR PUSTAKA

Badan Pusat Statistik Provinsi Nusa Tenggara Timur, 2021-2025. Jumlah Balita Stunting menurut Kabupaten/Kota. Kupang: BPS Provinsi NTT. Tersedia di: https://ntt.bps.go.id/id/statistics-table/2/MTQ4OSMy/jumlah-balita-stunting-menurut-kabupaten-kota.html
Anselin, L. (1995). Local Indicators of Spatial Association—LISA. Geographical Analysis, 27(2), 93–115. https://doi.org/10.1111/j.1538-4632.1995.tb00338.x
Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). Springer. https://doi.org/10.1007/978-1-4614-7618-4
Burnham, K. P., & Anderson, D. R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach (2nd ed.). Springer. https://doi.org/10.1007/b97636
Pinheiro, J. C., & Bates, D. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer. https://doi.org/10.1007/b98882
Rothman, K. J., Greenland, S., & Lash, T. L. (2021). Modern Epidemiology (4th ed.). Wolters Kluwer. ISBN: 9781975169274.