Hubungan Asosiasi Multidimensi antara Jenis Sarana Air Minum dan Klasifikasi Wilayah terhadap Kualitas Air Minum Rumah Tangga di Indonesia
Algafuri Decembra Nawa & Zhafran Fakhry Sani
09 June 2026
Abstrak
1 Pendahuluan
Air minum yang aman merupakan kebutuhan dasar manusia sekaligus pilar utama dalam menjaga kesehatan masyarakat secara berkelanjutan. Secara global, penyediaan akses terhadap air minum layak terus mengalami peningkatan, namun kesenjangan antara ketersediaan sarana yang “layak” secara fisik dan keamanan air yang “aman” secara mikrobiologis tetap menjadi persoalan krusial yang belum terselesaikan.
Pada tahun 2022, WHO dan UNICEF melaporkan bahwa sedikitnya 1,7 miliar penduduk dunia masih mengonsumsi air dari sumber yang terkontaminasi feses. Ketimpangan geografis membuat akses air minum perkotaan dan perdesaan menjadi terlihat jelas, terutama antara wilayah perkotaan yang mencatat akses air minum aman sebesar 81% dan wilayah perdesaan yang tertinggal di angka 62%. Sebagai negara kepulauan terbesar dengan populasi lebih dari 270 juta jiwa, Indonesia menghadapi tantangan serupa dalam pemenuhan komitmen global target Sustainable Development Goals (SDGs) 6.1.1, yaitu menyediakan akses merata terhadap air minum yang aman dan terjangkau bagi seluruh lapisan masyarakat pada tahun 2030.
Konteks Indonesia — SKAM-RT 2020: Hasil Studi Kualitas Air Minum Rumah Tangga (SKAM-RT) tahun 2020 yang diselenggarakan oleh Balitbangkes Kementerian Kesehatan RI menunjukkan realitas yang masih jauh dari target tersebut. Meskipun capaian akses SAM layak secara nasional telah menyentuh angka 93%, hanya 18,1% rumah tangga yang memiliki air minum yang benar-benar aman secara bakteriologis dari kontaminasi bakteri E. coli. Selain itu, akses air minum aman di perkotaan (15%) bernilai hampir dua kali lipat dibandingkan wilayah perdesaan (8%). Kontaminasi bakteriologis ini bukan sekadar masalah teknis infrastruktur, melainkan ancaman kesehatan publik yang serius mengingat diare akibat air tidak higienis tetap menjadi salah satu penyebab utama kematian balita di Indonesia.
Fenomena kontaminasi air ini melibatkan interaksi multidimensional dari berbagai karakteristik kategorik, seperti jenis sarana air minum yang digunakan dan klasifikasi wilayah tempat tinggal. Sejumlah penelitian yang ada umumnya menggunakan pendekatan analisis bivariat (seperti uji Chi-Square independensi biasa) atau model regresi logistik. Pendekatan tersebut menempatkan Kualitas Air secara kaku sebagai variabel respons yang dipengaruhi oleh variabel prediktor secara searah (sebab-akibat). Padahal, dalam kajian kemasyarakatan dan kesehatan masyarakat, hubungan antarvariabel kategorik tersebut sering kali bersifat timbal balik dan membentuk suatu struktur asosiasi yang kompleks di dalam tabel kontingensi multidimensi. Selain itu, uji bivariat parsial mengabaikan keberadaan variabel ketiga, sehingga berisiko menimbulkan bias akibat lewatnya efek interaksi bersyarat.
Penelitian ini menerapkan pendekatan Model Log-Linear Tiga Arah dengan memanfaatkan basis data berskala nasional dari SKAM-RT 2020. Berbeda dengan model regresif atau bivariat biasa yang membedakan variabel dependen dan independen secara kaku, Model Log-Linear memperlakukan ketiga variabel Klasifikasi Wilayah (X), Jenis Sarana Air Minum (Y), dan Status Kualitas Air Bakteriologis (Z) secara setara tanpa pembedaan respons-prediktor. Fokus utama analisis ini bukan untuk memprediksi pengaruh acak satu arah, melainkan untuk membedah struktur hubungan asosiasi, menguji pola independensi (baik independensi total, joint, maupun kondisional), serta mengidentifikasi keberadaan interaksi dua arah maupun interaksi tiga arah di antara ketiga variabel tersebut secara simultan menggunakan kriteria Goodness-of-Fit berbasis Deviance (G²) dan Pearson Chi-Square (X²). Melalui pendekatan pemodelan log-linear hierarkis ini, penelitian diharapkan dapat membuktikan secara empiris apakah kekuatan hubungan antara jenis sarana air dengan kualitas bakteriologis air bersifat konstan atau berubah (bergantung) pada latar geografis wilayahnya melalui pembuktian signifikansi interaksi [XYZ]. Hasil pemetaan struktur asosiasi homogen [XY][XZ][YZ] ini diharapkan dapat menjadi landasan ilmiah berbasis bukti (evidence-based policy) bagi para pengambil kebijakan untuk menentukan prioritas intervensi infrastruktur sanitasi yang tepat sasaran demi mempercepat pencapaian target air minum aman SDGs 2030 di Indonesia.
2 Metode Penelitian
2.1 Sumber Data dan Karakteristik Variabel
Data yang digunakan merupakan data sekunder yang bersumber dari Studi Kualitas Air Minum Rumah Tangga (SKAM-RT) 2020 yang dilaksanakan oleh Badan Penelitian dan Pengembangan Kesehatan (Balitbangkes) Kementerian Kesehatan RI. Desain penelitian yang digunakan adalah potong lintang (cross-sectional) dengan cakupan nasional yang tersebar di 34 provinsi di Indonesia. Total sampel riil domestik yang dianalisis dalam penelitian ini sebanyak N = 18.020 rumah tangga.
Analisis data dilakukan menggunakan tiga variabel kategorik utama yang dibentuk menjadi sebuah struktur tabel kontingensi tiga arah (X × Y × Z). Definisi operasional beserta skala pengukuran dari masing-masing variabel disajikan pada Tabel 1 berikut.
Tabel 1.1. Karakteristik Variabel
2.1.1 Data Observasi Tabel Kontingensi 2 × 2 × 3
# ── LOAD PACKAGES ─────────────────────────────────────────────────────────────
required_packages <- c("dplyr", "tidyr", "ggplot2", "nnet", "broom", "knitr",
"scales")
missing_pkg <- required_packages[!sapply(required_packages, requireNamespace,
quietly = TRUE)]
if (length(missing_pkg) > 0) install.packages(missing_pkg)
suppressPackageStartupMessages({
library(dplyr); library(tidyr); library(ggplot2)
library(nnet); library(broom); library(knitr); library(scales)
})
# ── INPUT DATA RIIL SKAM-RT 2020 ──────────────────────────────────────────────
data_riil <- data.frame(
Daerah = c(
rep("Perkotaan", 6), rep("Perdesaan", 6)
),
Jenis_SAM = c(
rep("Air Terlindungi", 3), rep("Air Tidak Terlindungi", 3),
rep("Air Terlindungi", 3), rep("Air Tidak Terlindungi", 3)
),
Kualitas_Air = rep(c("Aman", "Layak", "Tidak Layak"), 4),
Jumlah = c(
1102, 5925, 173, # Perkotaan + Terlindungi
1, 24, 175, # Perkotaan + Tidak Terlindungi
898, 8526, 1396, # Perdesaan + Terlindungi
2, 116, 1702 # Perdesaan + Tidak Terlindungi
)
)
data_riil <- data_riil %>%
mutate(
Kualitas_Air = factor(Kualitas_Air,
levels = c("Aman", "Layak", "Tidak Layak")),
Jenis_SAM = factor(Jenis_SAM,
levels = c("Air Terlindungi", "Air Tidak Terlindungi")),
Daerah = factor(Daerah, levels = c("Perkotaan", "Perdesaan"))
)Tabel 1.2. Data Variabel
2.2 Metode Analisis Data
Metode analisis data yang digunakan untuk melihat interaksi multidimensi antarvariabel adalah Model Log-Linear Tiga Arah (Three-Way Log-Linear Model). Model log-linear memperlakukan seluruh variabel (X, Y, Z) secara simetris sebagai variabel respons tanpa adanya partisi kaku berupa variabel dependen-independen. Pemodelan ini bertujuan utama untuk mengestimasi parameter intensitas hubungan asosiasi serta pola interaksi yang terjadi di dalam sel tabel kontingensi multidimensi. Langkah-langkah pemodelan yang dilakukan sebagai berikut:
2.2.1 Analisis Deskriptif Kontingensi Terstrata dan Pembentukan Margin
Menganalisis frekuensi observasi nyata (nijk) pada tabel kontingensi berukuran I × J × K (2 × 2 × 3), di mana i=1,2 (indeks klasifikasi wilayah X), j=1,2 (indeks jenis sarana air minum Y), dan k=1,2,3 (indeks status kualitas air mikrobiologis Z). Berdasarkan teori struktur data hitung (count data), langkah awal dilakukan dengan membentuk margin satu arah melalui penjumlahan sel:
Margin satu arah:
ni++ = Σj Σk nijk ; n+j+ = Σi
Σk nijk ; n++k = Σi Σj nijk
Serta pembentukan margin dua arah (nij+, ni+k, n+jk) yang merepresentasikan sub-total frekuensi berpasangan antarvariabel.
2.2.2 Pemodelan Mutual Independence dan Perhitungan Ekspektasi Manual
Membangun model dasar awal berupa Model Independensi Total (Mutual Independence) yang mengasumsikan tidak terdapat hubungan asosiasi sama sekali di antara ketiga variabel ([X][Y][Z]). Konsep probabilitas gabungan independen dengan nilai frekuensi harapan teoritis (μijk) di bawah kondisi H₀. Maka, independensi total dihitung secara manual dengan rumusan:
μijk = (ni++ × n+j+ × n++k) / N²
Di mana N melambangkan ukuran sampel keseluruhan objek penelitian (N = 18.020).
2.2.3 Estimasi Parameter Menggunakan Metode GLM Poisson
Estimasi parameter menggunakan metode Generalized Linear Models (GLM) di mana frekuensi sel nijk diasumsikan berdistribusi Poisson secara independen dengan rata-rata μijk. Fungsi tautan (link function) logaritma natural digunakan untuk mentransformasikan model multiplikatif menjadi model aditif linear sebagai berikut:
Model Mutual Independence [X][Y][Z]:
ln(μijk) =
λ + λi^X + λj^Y + λk^Z
Model Asosiasi Homogen [XY][XZ][YZ]:
ln(μijk) = λ + λi^X + λj^Y + λk^Z + λij^XY + λik^XZ + λjk^YZ
Model Jenuh (Saturated Model) [XYZ]:
ln(μijk) = λ + λi^X + λj^Y + λk^Z + λij^XY + λik^XZ + λjk^YZ +
λijk^XYZ
Di mana λ menyatakan efek rata-rata umum; λiX, λjY, λkZ menyatakan efek utama tunggal; λijXY, λikXZ, λjkYZ menyatakan parameter interaksi asosiasi dua arah; dan λijkXYZ menyatakan parameter interaksi tiga arah tingkat tinggi. Estimasi parameter menggunakan pembatasan kategori acuan (corner-point constraints) dengan menetapkan efek kategori terakhir bernilai 0.
2.2.4 Pemilihan Model Terbaik dengan Backward Elimination
Prosedur seleksi hierarkis dilakukan untuk memperoleh model parsimoni terbaik dengan mengevaluasi nilai Deviance (G²) global yang bersumber dari rasio fungsi likelihood (Likelihood Ratio Test / LRT). Rumusan G² global dituliskan sebagai berikut:
Deviance global (G²):
G² = 2 Σi Σj Σk [ nijk
ln(nijk / μijk) ]
Pengujian dilakukan bertahap menggunakan selisih deviasi antarmodel. Untuk uji interaksi tiga arah (Homogeneous vs Saturated), hipotesis yang diuji adalah:
H₀: λijkXYZ = 0 (Model Asosiasi Homogen memadai) vs H₁: λijkXYZ ≠ 0 (Model Jenuh diperlukan). Statistik uji ditentukan melalui:
ΔG² = G²(Model Homogen) − G²(Model Saturated) ; Δdf = df(Model Homogen) − df(Model Saturated)
Jika ΔG² ≤ χ²(α, Δdf), maka H₀ gagal ditolak, yang berarti interaksi tiga arah tidak signifikan secara statistik sehingga model disederhanakan menjadi model asosiasi homogen. Metode yang sama diterapkan untuk menguji signifikansi parsial interaksi dua arah dengan membandingkan model homogen (full model) terhadap model bersyarat yang dikurangi (reduced model):
ΔG² = G²(Model Reduced) − G²(Model Homogen)
2.2.5 Estimasi Model Final Menggunakan Algoritma Iterative Proportional Fitting (IPF)
Memvalidasi kestabilan nilai parameter harapan menggunakan algoritma non-parametrik Iterative Proportional Fitting (IPF) yang secara khusus dirancang untuk model log-linear. Algoritma ini melakukan penyesuaian nilai harapan m̂ijk secara berulang (siklis) langsung ke margin observasi konstan (nij+, ni+k, n+jk) tanpa estimasi parameter λ terlebih dahulu. Proses iterasi (m) terus berjalan hingga mencapai batas konvergensi terkecil (ε = 10⁻¹⁰) yang diukur dari kesalahan maksimum sel:
Max Error(m) = max( |m̂ij+(m) − nij+|, |m̂i+k(m) − ni+k|, |m̂+jk(m) − n+jk| ) < ε
2.2.6 Evaluasi Diagnostik Residual Model
Ketepatan kecocokan model pada tingkat individu sel dievaluasi menggunakan komponen sisaan Pearson Terstandardisasi (Standardized Pearson Residual). Berdasarkan formulasi sisaan untuk data counts, persamaannya dirumuskan sebagai:
eijk = (nijk − μ̂ijk) / √( μ̂ijk (1 − hijk) )
Di mana hijk melambangkan nilai diagonal matriks leverage (topi). Sel dengan nilai |eijk| > 2 mengindikasikan adanya penyimpangan kecocokan lokal (outlier sel). Selain itu, kontribusi komponen deviance per sel (G²ijk) terhadap total deviance dihitung melalui formula:
G²ijk = 2 [ nijk ln(nijk / μ̂ijk) − (nijk − μ̂ijk) ]
2.2.7 Hubungan Melalui Odds Ratio (OR) Kondisional
Dengan konsep asosiasi parsial, interpretasi kekuatan hubungan log-linear dianalisis menggunakan Odds Ratio (OR) kondisional (bersyarat) pada tingkat strata variabel ketiga yang dibuat konstan. Rumus OR kondisional untuk hubungan variabel X dan Y pada strata Z=k dirumuskan sebagai:
OR kondisional asosiasi XY pada strata Z = k:
θXY(k) = (μ11k × μ22k) / (μ12k × μ21k)
Dalam struktur Model Asosiasi Homogen, sifat mutlak yang berlaku adalah nilai OR kondisional bernilai seragam/konstan di seluruh tingkatan strata pengondisinya (θXY(1) = θXY(2) = … = θXY(K)). Nilai asosiasi konstan ini memiliki hubungan matematis langsung dengan parameter interaksi dua arah model aditifnya, yaitu θXY = exp(λijXY).
3 Hasil dan Pembahasan
3.1 Pemilihan Model Log-Linear Terbaik
Langkah awal dalam analisis data kategori tiga dimensi adalah menentukan struktur hubungan terbaik antar variabel: Daerah (X), Jenis Sarana Air Minum (Y), dan Kualitas Air (Z). Evaluasi dilakukan terhadap sembilan model log-linear potensial, mulai dari model independensi lengkap [X][Y][Z] hingga model jenuh (saturated model) [XYZ]. Nilai Akaike Information Criterion (AIC) dan Bayesian Information Criterion (BIC) dari masing-masing model dirangkum pada tabel berikut.
3.1.1 Perbandingan Sembilan Model Kandidat
# ── KONSTRUKSI TABEL KONTINGENSI ──────────────────────────────────────────────
ct <- xtabs(Jumlah ~ Daerah + Jenis_SAM + Kualitas_Air, data = data_riil)
# ── FIT SEMBILAN MODEL LOG-LINEAR MENGGUNAKAN GLM POISSON ────────────────────
# Bentuk data flat untuk GLM
df_flat <- as.data.frame(ct)
df_flat <- df_flat[df_flat$Freq > 0 | TRUE, ] # pertahankan sel nol
names(df_flat)[names(df_flat) == "Freq"] <- "n"
# Definisi formula masing-masing model
formulas <- list(
"[X][Y][Z]" = n ~ Daerah + Jenis_SAM + Kualitas_Air,
"[X][YZ]" = n ~ Daerah + Jenis_SAM * Kualitas_Air,
"[Y][XZ]" = n ~ Jenis_SAM + Daerah * Kualitas_Air,
"[Z][XY]" = n ~ Kualitas_Air + Daerah * Jenis_SAM,
"[XY][XZ]" = n ~ Daerah * Jenis_SAM + Daerah * Kualitas_Air,
"[XY][YZ]" = n ~ Daerah * Jenis_SAM + Jenis_SAM * Kualitas_Air,
"[XZ][YZ]" = n ~ Daerah * Kualitas_Air + Jenis_SAM * Kualitas_Air,
"[XY][XZ][YZ]" = n ~ Daerah * Jenis_SAM + Daerah * Kualitas_Air +
Jenis_SAM * Kualitas_Air,
"[XYZ] (Saturated)" = n ~ Daerah * Jenis_SAM * Kualitas_Air
)
fits <- lapply(formulas, function(f) glm(f, data = df_flat, family = poisson()))
# ── RINGKASAN AIC & BIC ────────────────────────────────────────────────────────
model_summary <- data.frame(
Model = names(formulas),
Spesifikasi = names(formulas),
Deviance = sapply(fits, deviance),
df = sapply(fits, function(m) m$df.residual),
AIC = sapply(fits, AIC),
BIC = sapply(fits, BIC)
)
kable(model_summary[, c("Model","Deviance","df","AIC","BIC")],
digits = 3,
caption = "Tabel 1. Perbandingan Nilai AIC dan BIC antar Model Log-Linear",
col.names = c("Model Log-Linear", "Deviance (G²)", "df", "AIC", "BIC"))| Model Log-Linear | Deviance (G²) | df | AIC | BIC | |
|---|---|---|---|---|---|
| [X][Y][Z] | [X][Y][Z] | 8421.207 | 7 | 8518.684 | 8521.108 |
| [X][YZ] | [X][YZ] | 1696.012 | 5 | 1797.489 | 1800.883 |
| [Y][XZ] | [Y][XZ] | 6765.102 | 5 | 6866.579 | 6869.974 |
| [Z][XY] | [Z][XY] | 7579.426 | 6 | 7678.903 | 7681.813 |
| [XY][XZ] | [XY][XZ] | 5923.322 | 4 | 6026.799 | 6030.678 |
| [XY][YZ] | [XY][YZ] | 854.231 | 4 | 957.708 | 961.588 |
| [XZ][YZ] | [XZ][YZ] | 39.908 | 3 | 145.385 | 149.749 |
| [XY][XZ][YZ] | [XY][XZ][YZ] | 18.920 | 2 | 126.397 | 131.246 |
| [XYZ] (Saturated) | [XYZ] (Saturated) | 0.000 | 0 | 111.477 | 117.296 |
Temuan: Berdasarkan tabel di atas, model jenuh [XYZ] menghasilkan nilai AIC minimum (111,477) dan BIC minimum (117,296). Untuk menguji apakah interaksi tiga faktor [XYZ] signifikan secara statistik, dilakukan uji selisih deviance (likelihood ratio test) antara model asosiasi homogen [XY][XZ][YZ] dan model jenuh [XYZ].
3.2 Uji Signifikansi Interaksi Tiga Faktor [XYZ]
# ── UJI SELISIH DEVIANCE: [XY][XZ][YZ] vs [XYZ] ──────────────────────────────
m_homogen <- fits[["[XY][XZ][YZ]"]]
m_saturated <- fits[["[XYZ] (Saturated)"]]
delta_G2 <- deviance(m_homogen) - deviance(m_saturated)
delta_df <- m_homogen$df.residual - m_saturated$df.residual
p_value <- pchisq(delta_G2, df = delta_df, lower.tail = FALSE)
tbl_uji <- data.frame(
Perbandingan = "[XY][XZ][YZ] vs [XYZ]",
Delta_G2 = round(delta_G2, 4),
Delta_df = delta_df,
P_Value = format(p_value, scientific = TRUE, digits = 4)
)
kable(tbl_uji,
caption = "Tabel 2. Pengujian Interaksi Tiga Faktor ([XY][XZ][YZ] vs [XYZ])",
col.names = c("Perbandingan Model", "ΔG²", "Δdf", "P-Value"))| Perbandingan Model | ΔG² | Δdf | P-Value |
|---|---|---|---|
| [XY][XZ][YZ] vs [XYZ] | 18.9199 | 2 | 7.791e-05 |
#>
#> ΔG² = 18.9199
#> Δdf = 2
#> P-Value = 7.791036e-05Interpretasi: Hasil pengujian menunjukkan bahwa ΔG² = 18,9199 dengan Δdf = 2 memberikan p-value sebesar 0,0001 (< 0,05). Dengan demikian, hipotesis nol ditolak, yang berarti interaksi tiga faktor [XYZ] signifikan secara statistik. Hal ini menandakan bahwa pola hubungan dua faktor (misalnya antara karakteristik daerah dan kualitas air) tidak konstan (heterogen) di seluruh tingkatan variabel ketiga. Meskipun model jenuh merupakan model terbaik secara matematis, model asosiasi homogen [XY][XZ][YZ] tetap dianalisis secara mendalam sebagai baseline untuk mengidentifikasi sel-sel spesifik yang memicu munculnya interaksi tingkat tinggi tersebut.
3.3 Analisis Konvergensi dan Komparasi Algoritma
Pada model log-linear non-saturasi seperti model asosiasi homogen [XY][XZ][YZ], nilai harapan frekuensi sel (μijk) tidak memiliki bentuk eksplisit (closed-form) sehingga memerlukan algoritma iteratif. Penelitian ini membandingkan efisiensi algoritma Iterative Proportional Fitting (IPF) dengan pendekatan Generalized Linear Model (GLM) berbasis fungsi link log distribusi Poisson.
# ── IMPLEMENTASI MANUAL ALGORITMA IPF ─────────────────────────────────────────
# Inisialisasi dari nilai seragam (semua sel = 1) agar riwayat konvergensi
# bermakna secara visual — bukan dari observasi yang sudah "dekat"
ipf_loglinear <- function(n_obs, max_iter = 100, tol = 1e-8) {
# Inisialisasi: semua sel bernilai 1 (uniform start)
m <- array(1, dim = dim(n_obs), dimnames = dimnames(n_obs))
marg_xy <- apply(n_obs, c(1,2), sum)
marg_xz <- apply(n_obs, c(1,3), sum)
marg_yz <- apply(n_obs, c(2,3), sum)
hist_err <- numeric(max_iter)
for (iter in seq_len(max_iter)) {
# Langkah 1 — sesuaikan ke margin XY
fitted_xy <- apply(m, c(1,2), sum)
for (k in seq_len(dim(m)[3]))
m[,,k] <- m[,,k] * (marg_xy / pmax(fitted_xy, 1e-12))
# Langkah 2 — sesuaikan ke margin XZ
fitted_xz <- apply(m, c(1,3), sum)
for (j in seq_len(dim(m)[2]))
m[,j,] <- m[,j,] * (marg_xz / pmax(fitted_xz, 1e-12))
# Langkah 3 — sesuaikan ke margin YZ
fitted_yz <- apply(m, c(2,3), sum)
for (i in seq_len(dim(m)[1]))
m[i,,] <- m[i,,] * (marg_yz / pmax(fitted_yz, 1e-12))
# Hitung maximum absolute error ke ketiga margin target
err_xy <- max(abs(apply(m, c(1,2), sum) - marg_xy))
err_xz <- max(abs(apply(m, c(1,3), sum) - marg_xz))
err_yz <- max(abs(apply(m, c(2,3), sum) - marg_yz))
hist_err[iter] <- max(err_xy, err_xz, err_yz)
if (hist_err[iter] < tol) {
hist_err <- hist_err[seq_len(iter)]
break
}
}
list(fitted = m, history = hist_err)
}
# Jalankan IPF pada tabel kontingensi 2×2×3
n_array <- ct
ipf_result <- ipf_loglinear(n_array)
# ── RIWAYAT KONVERGENSI ────────────────────────────────────────────────────────
hist_df <- data.frame(
Iterasi = seq_along(ipf_result$history),
MaxError = ipf_result$history
)
kable(hist_df,
digits = 4,
caption = "Tabel 3. Riwayat Iterasi dan Maximum Error Margin Algoritma IPF",
col.names = c("Iterasi", "Maximum Error Margin"))| Iterasi | Maximum Error Margin |
|---|---|
| 1 | 485.4461 |
| 2 | 11.1570 |
| 3 | 5.1512 |
| 4 | 1.9671 |
| 5 | 0.7414 |
| 6 | 0.2786 |
| 7 | 0.1046 |
| 8 | 0.0393 |
| 9 | 0.0147 |
| 10 | 0.0055 |
| 11 | 0.0021 |
| 12 | 0.0008 |
| 13 | 0.0003 |
| 14 | 0.0001 |
| 15 | 0.0000 |
| 16 | 0.0000 |
| 17 | 0.0000 |
| 18 | 0.0000 |
| 19 | 0.0000 |
| 20 | 0.0000 |
| 21 | 0.0000 |
| 22 | 0.0000 |
| 23 | 0.0000 |
| 24 | 0.0000 |
tol_line <- min(hist_df$MaxError) * 1.5 # garis toleransi di dekat nilai konvergen
ggplot(hist_df, aes(x = Iterasi, y = MaxError)) +
geom_line(color = "#17a2b8", linewidth = 1.3) +
geom_point(color = "#0a7c8c", size = 3.5, shape = 21,
fill = "#a8e6ef", stroke = 1.2) +
geom_hline(yintercept = tol_line, linetype = "dashed",
color = "#b84f5a", linewidth = 0.8) +
annotate("text",
x = max(hist_df$Iterasi) * 0.55,
y = tol_line * 2.5,
label = paste0("Konvergen pada iterasi ke-", nrow(hist_df)),
color = "#b84f5a", size = 3.5, fontface = "italic") +
scale_x_continuous(breaks = hist_df$Iterasi) +
scale_y_log10(labels = scales::scientific_format()) +
labs(
title = "Kurva Konvergensi Algoritma IPF",
subtitle = paste0("Model Asosiasi Homogen [XY][XZ][YZ] — ",
nrow(hist_df), " iterasi hingga konvergen"),
x = "Iterasi ke-",
y = "Maximum Error Margin (Skala Log)"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", color = "#0a7c8c", size = 13),
plot.subtitle = element_text(color = "#64748b", size = 10),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "#e0f6fa"),
axis.text = element_text(color = "#1a2e35")
)
Gambar 2. Kurva Konvergensi Algoritma IPF
# ── KOMPARASI FITTED VALUES: IPF vs GLM ───────────────────────────────────────
# Reorder IPF fitted array ke urutan baris df_flat (Daerah × Jenis_SAM × Kualitas_Air)
ipf_arr <- ipf_result$fitted # array [Daerah, Jenis_SAM, Kualitas_Air]
fitted_ipf_vec <- mapply(function(d, s, k)
ipf_arr[as.character(d), as.character(s), as.character(k)],
df_flat$Daerah, df_flat$Jenis_SAM, df_flat$Kualitas_Air)
fitted_glm_vec <- fitted(m_homogen)
komparasi_df <- data.frame(
Daerah = df_flat$Daerah,
Jenis_SAM = df_flat$Jenis_SAM,
Kualitas = df_flat$Kualitas_Air,
Freq = df_flat$n,
Fitted_IPF = round(fitted_ipf_vec, 4),
Fitted_GLM = round(fitted_glm_vec, 4),
Selisih = round(abs(fitted_ipf_vec - fitted_glm_vec), 8)
)
kable(komparasi_df,
digits = 4,
caption = "Tabel 4. Komparasi Nilai Harapan antara Algoritma IPF dan GLM",
col.names = c("Daerah","Jenis SAM","Kualitas Air","Freq","Fitted IPF",
"Fitted GLM","Selisih |IPF-GLM|"))| Daerah | Jenis SAM | Kualitas Air | Freq | Fitted IPF | Fitted GLM | Selisih |IPF-GLM| |
|---|---|---|---|---|---|---|
| Perkotaan | Air Terlindungi | Aman | 1102 | 1101.6742 | 1101.6742 | 0 |
| Perdesaan | Air Terlindungi | Aman | 898 | 898.3258 | 898.3258 | 0 |
| Perkotaan | Air Tidak Terlindungi | Aman | 1 | 1.3258 | 1.3258 | 0 |
| Perdesaan | Air Tidak Terlindungi | Aman | 2 | 1.6742 | 1.6742 | 0 |
| Perkotaan | Air Terlindungi | Layak | 5925 | 5905.8014 | 5905.8014 | 0 |
| Perdesaan | Air Terlindungi | Layak | 8526 | 8545.1986 | 8545.1986 | 0 |
| Perkotaan | Air Tidak Terlindungi | Layak | 24 | 43.1986 | 43.1986 | 0 |
| Perdesaan | Air Tidak Terlindungi | Layak | 116 | 96.8014 | 96.8014 | 0 |
| Perkotaan | Air Terlindungi | Tidak Layak | 173 | 192.5244 | 192.5244 | 0 |
| Perdesaan | Air Terlindungi | Tidak Layak | 1396 | 1376.4756 | 1376.4756 | 0 |
| Perkotaan | Air Tidak Terlindungi | Tidak Layak | 175 | 155.4756 | 155.4756 | 0 |
| Perdesaan | Air Tidak Terlindungi | Tidak Layak | 1702 | 1721.5244 | 1721.5244 | 0 |
Temuan Penting — Ekuivalensi Matematis: Hasil komparasi menunjukkan bahwa selisih mutlak antara Fitted_IPF dan Fitted_GLM adalah tepat 0 untuk keseluruhan 12 sel pengamatan. Hal ini secara empiris membuktikan kesetaraan matematis (mathematical equivalence) kedua metode dalam mengestimasi parameter model log-linear non-saturasi, sehingga memberikan jaminan reliabilitas tinggi pada nilai parameter yang dihasilkan.
3.4 Uji Kebaikan Suai Global Model Asosiasi Homogen
Uji kebaikan suai (goodness-of-fit) digunakan untuk mengevaluasi apakah model asosiasi homogen [XY][XZ][YZ] mampu merepresentasikan struktur data observasi secara keseluruhan.
# ── STATISTIK GOODNESS-OF-FIT MODEL [XY][XZ][YZ] ─────────────────────────────
G2_hom <- deviance(m_homogen)
df_hom <- m_homogen$df.residual
p_G2 <- pchisq(G2_hom, df = df_hom, lower.tail = FALSE)
# Pearson Chi-Square manual
pearson_x2 <- sum((df_flat$n - fitted(m_homogen))^2 / fitted(m_homogen))
p_X2 <- pchisq(pearson_x2, df = df_hom, lower.tail = FALSE)
gof_tbl <- data.frame(
Statistik = c("Deviance (G²)", "Pearson (X²)"),
Nilai = round(c(G2_hom, pearson_x2), 4),
df = c(df_hom, df_hom),
P_Value = format(c(p_G2, p_X2), scientific = TRUE, digits = 4)
)
kable(gof_tbl,
caption = "Tabel 5. Statistik Kebaikan Suai Model Asosiasi [XY][XZ][YZ]",
col.names = c("Statistik Uji","Nilai","Derajat Bebas","P-Value"))| Statistik Uji | Nilai | Derajat Bebas | P-Value |
|---|---|---|---|
| Deviance (G²) | 18.9199 | 2 | 7.791e-05 |
| Pearson (X²) | 17.5195 | 2 | 1.569e-04 |
Interpretasi: Hasil pada tabel di atas menunjukkan bahwa baik statistik Deviance (G² = 18,9199) maupun Pearson (X² = 17,5195) menghasilkan p-value yang sangat signifikan (< 0,05). Akibatnya, hipotesis nol yang menyatakan model asosiasi homogen cocok dengan data terpaksa ditolak. Penolakan global ini mengindikasikan adanya fluktuasi lokal yang ekstrem atau deviasi dari asumsi kehomogenan odds rasio pada sel-sel tertentu.
3.5 Analisis Diagnostik Sisaan dan Kontribusi Sel
Untuk menelusuri sumber ketidakcocokan model [XY][XZ][YZ], dilakukan analisis diagnostik sisaan (residuals) dan kontribusi deviance per sel.
# ── PEARSON & DEVIANCE RESIDUALS ──────────────────────────────────────────────
obs <- df_flat$n
fitted_ <- fitted(m_homogen)
pearson_res <- (obs - fitted_) / sqrt(fitted_)
dev_res <- sign(obs - fitted_) * sqrt(2 * abs(
ifelse(obs > 0, obs * log(obs / fitted_), 0) - (obs - fitted_)
))
# Leverage diagonal (hii) via hatvalues (base stats)
h_ii <- hatvalues(m_homogen)
std_res <- pearson_res / sqrt(pmax(1 - h_ii, 1e-8))
resid_df <- df_flat %>%
mutate(
Fitted = round(fitted_, 4),
Pearson_Res = round(pearson_res, 4),
Std_Res = round(std_res, 4),
Dev_Res = round(dev_res, 4)
) %>%
dplyr::select(Daerah, Jenis_SAM, Kualitas_Air, n, Fitted,
Pearson_Res, Std_Res, Dev_Res)
kable(resid_df,
caption = "Tabel 6. Analisis Sisaan Lengkap Model Asosiasi Homogen",
col.names = c("Daerah","Jenis SAM","Kualitas Air","Freq","Fitted",
"Pearson Res","Std Res","Dev Res"))| Daerah | Jenis SAM | Kualitas Air | Freq | Fitted | Pearson Res | Std Res | Dev Res |
|---|---|---|---|---|---|---|---|
| Perkotaan | Air Terlindungi | Aman | 1102 | 1101.6742 | 0.0098 | 0.3803 | 0.0098 |
| Perdesaan | Air Terlindungi | Aman | 898 | 898.3258 | -0.0109 | -0.3803 | -0.0109 |
| Perkotaan | Air Tidak Terlindungi | Aman | 1 | 1.3258 | -0.2829 | -0.3803 | -0.2959 |
| Perdesaan | Air Tidak Terlindungi | Aman | 2 | 1.6742 | 0.2518 | 0.3803 | 0.2442 |
| Perkotaan | Air Terlindungi | Layak | 5925 | 5905.8014 | 0.2498 | 4.1434 | 0.2497 |
| Perdesaan | Air Terlindungi | Layak | 8526 | 8545.1986 | -0.2077 | -4.1434 | -0.2078 |
| Perkotaan | Air Tidak Terlindungi | Layak | 24 | 43.1986 | -2.9210 | -4.1434 | -3.1914 |
| Perdesaan | Air Tidak Terlindungi | Layak | 116 | 96.8014 | 1.9513 | 4.1434 | 1.8916 |
| Perkotaan | Air Terlindungi | Tidak Layak | 173 | 192.5244 | -1.4071 | -4.1819 | -1.4320 |
| Perdesaan | Air Terlindungi | Tidak Layak | 1396 | 1376.4756 | 0.5263 | 4.1819 | 0.5250 |
| Perkotaan | Air Tidak Terlindungi | Tidak Layak | 175 | 155.4756 | 1.5658 | 4.1819 | 1.5347 |
| Perdesaan | Air Tidak Terlindungi | Tidak Layak | 1702 | 1721.5244 | -0.4706 | -4.1819 | -0.4715 |
# ── KONTRIBUSI DEVIANCE PER SEL ───────────────────────────────────────────────
dev_contrib <- 2 * (ifelse(obs > 0, obs * log(obs / fitted_), 0) - (obs - fitted_))
pct_contrib <- dev_contrib / sum(dev_contrib) * 100
kontrib_df <- df_flat %>%
mutate(
Fitted = round(fitted_, 4),
Dev_Res = round(dev_res, 4),
Kontribusi = round(dev_contrib, 5),
Pct_G2 = round(pct_contrib, 2)
) %>%
arrange(desc(Kontribusi)) %>%
dplyr::select(Daerah, Jenis_SAM, Kualitas_Air, n, Fitted, Dev_Res, Kontribusi, Pct_G2)
kable(kontrib_df,
caption = "Tabel 7. Kontribusi dan Proporsi Sisa Deviance per Sel terhadap G² Global",
col.names = c("Daerah","Jenis SAM","Kualitas Air","Freq","Fitted",
"Dev Res","Kontribusi G²","Proporsi G² (%)"))| Daerah | Jenis SAM | Kualitas Air | Freq | Fitted | Dev Res | Kontribusi G² | Proporsi G² (%) |
|---|---|---|---|---|---|---|---|
| Perkotaan | Air Tidak Terlindungi | Layak | 24 | 43.1986 | -3.1914 | 10.18502 | 53.83 |
| Perdesaan | Air Tidak Terlindungi | Layak | 116 | 96.8014 | 1.8916 | 3.57828 | 18.91 |
| Perkotaan | Air Tidak Terlindungi | Tidak Layak | 175 | 155.4756 | 1.5347 | 2.35521 | 12.45 |
| Perkotaan | Air Terlindungi | Tidak Layak | 173 | 192.5244 | -1.4320 | 2.05057 | 10.84 |
| Perdesaan | Air Terlindungi | Tidak Layak | 1396 | 1376.4756 | 0.5250 | 0.27564 | 1.46 |
| Perdesaan | Air Tidak Terlindungi | Tidak Layak | 1702 | 1721.5244 | -0.4715 | 0.22227 | 1.17 |
| Perkotaan | Air Tidak Terlindungi | Aman | 1 | 1.3258 | -0.2959 | 0.08755 | 0.46 |
| Perkotaan | Air Terlindungi | Layak | 5925 | 5905.8014 | 0.2497 | 0.06234 | 0.33 |
| Perdesaan | Air Tidak Terlindungi | Aman | 2 | 1.6742 | 0.2442 | 0.05963 | 0.32 |
| Perdesaan | Air Terlindungi | Layak | 8526 | 8545.1986 | -0.2078 | 0.04317 | 0.23 |
| Perdesaan | Air Terlindungi | Aman | 898 | 898.3258 | -0.0109 | 0.00012 | 0.00 |
| Perkotaan | Air Terlindungi | Aman | 1102 | 1101.6742 | 0.0098 | 0.00010 | 0.00 |
Temuan Diagnostik: Terdapat empat sel utama yang
memiliki nilai sisaan standar melampaui batas kritis ±1,96 sekaligus
mendominasi sisa deviance model:
1. Sel S5 (Perkotaan -
Air Tidak Terlindungi - Layak): Memiliki sisaan standar -4,1434
dan menyumbang 53,83% dari total deviance. Frekuensi aktual (24) jauh
lebih kecil daripada estimasi model (43,1986), mengindikasikan model
secara signifikan overestimate terhadap kelompok ini.
2. Sel S11 (Perdesaan - Air Tidak Terlindungi - Layak):
Memiliki sisaan standar 4,1434 dan menyumbang 18,91% deviance. Model
bersifat underestimate karena frekuensi aktual (116) melampaui
nilai harapan (96,8014).
3. Sel S6 dan S3:
Masing-masing menyumbang 12,45% dan 10,84% terhadap deviance global
dengan sisaan standar ekstrem masing-masing sebesar 4,1819 dan
-4,1819.
Kombinasi dari keempat sel pencilan lokal (local
outliers) ini menyumbang 96,03% dari total
interaksi sisa yang tidak terjelaskan. Hal ini membuktikan bahwa
ketidakcocokan model homogeneous association berakar dari perilaku
spesifik kelompok rumah tangga pengguna “Air Tidak Terlindungi” lintas
daerah, serta disparitas kualitas “Tidak Layak” di wilayah
perkotaan.
``` r
resid_df_plot <- resid_df %>%
mutate(
Sel = paste(Daerah, Jenis_SAM, Kualitas_Air, sep = "\n"),
Outlier = abs(Std_Res) > 2
)
ggplot(resid_df_plot, aes(x = reorder(Sel, abs(Std_Res)), y = Std_Res,
fill = Outlier)) +
geom_col(width = 0.65, color = "white") +
geom_hline(yintercept = c(-2, 2), linetype = "dashed",
color = "#b84f5a", linewidth = 0.8) +
scale_fill_manual(values = c("TRUE" = "#0a7c8c", "FALSE" = "#a8e6ef"),
labels = c("Normal", "Outlier Sel (|e| > 2)")) +
coord_flip() +
labs(title = "Residual Standar per Sel Kontingensi",
subtitle = "Garis merah putus-putus = batas kritis ±2",
x = "Kombinasi Sel (X × Y × Z)", y = "Residual Standar",
fill = "Status Sel") +
theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", color = "#0a7c8c"),
plot.subtitle = element_text(color = "#64748b"),
legend.position = "bottom")
Gambar 3. Peta Panas Residual Standar per Sel (|e| > 2 = outlier)
3.6 Analisis Asosiasi Parsial
Untuk mengevaluasi kekuatan kepentingan relatif dari interaksi dua faktor, dilakukan uji asosiasi parsial dengan menghapus satu per satu komponen interaksi dari model [XY][XZ][YZ]. Hasil dari uji ini disajikan pada tabel berikut.
# ── UJI ASOSIASI PARSIAL: hapus masing-masing efek dua arah ──────────────────
m_no_YZ <- glm(n ~ Daerah * Jenis_SAM + Daerah * Kualitas_Air,
data = df_flat, family = poisson())
m_no_XZ <- glm(n ~ Daerah * Jenis_SAM + Jenis_SAM * Kualitas_Air,
data = df_flat, family = poisson())
m_no_XY <- glm(n ~ Daerah * Kualitas_Air + Jenis_SAM * Kualitas_Air,
data = df_flat, family = poisson())
delta_YZ <- deviance(m_no_YZ) - deviance(m_homogen)
delta_XZ <- deviance(m_no_XZ) - deviance(m_homogen)
delta_XY <- deviance(m_no_XY) - deviance(m_homogen)
parsial_df <- data.frame(
Interaksi_Dihapus = c("YZ (Jenis SAM × Kualitas Air)",
"XZ (Daerah × Kualitas Air)",
"XY (Daerah × Jenis SAM)"),
Model_Terreduksi = c("[XY][XZ]", "[XY][YZ]", "[XZ][YZ]"),
Delta_G2 = round(c(delta_YZ, delta_XZ, delta_XY), 4),
Signifikansi = c("Sangat Signifikan", "Sangat Signifikan", "Signifikan")
)
kable(parsial_df,
caption = "Tabel 8. Hasil Uji Dampak Penghapusan Komponen Interaksi Dua Faktor",
col.names = c("Interaksi Dihapus", "Model Terreduksi", "Peningkatan Deviance (ΔG²)", "Signifikansi Efek"))| Interaksi Dihapus | Model Terreduksi | Peningkatan Deviance (ΔG²) | Signifikansi Efek |
|---|---|---|---|
| YZ (Jenis SAM × Kualitas Air) | [XY][XZ] | 5904.4017 | Sangat Signifikan |
| XZ (Daerah × Kualitas Air) | [XY][YZ] | 835.3116 | Sangat Signifikan |
| XY (Daerah × Jenis SAM) | [XZ][YZ] | 20.9877 | Signifikan |
Tabel 8 mengindikasikan bahwa seluruh komponen interaksi dua faktor bermakna secara statistik. Penjelasan mendalam mengenai dampak masing-masing interaksi dijabarkan sebagai berikut:
3.6.1 Analisis Penghapusan Interaksi YZ
Penghapusan efek YZ menghasilkan lonjakan deviance terbesar, yaitu ΔG² = 5.904,40. Ini menunjukkan bahwa hubungan antara jenis sarana air minum dan kualitas air yang dihasilkan adalah efek paling dominan di dalam sistem variabel ini. Penghapusan interaksi YZ sangat memengaruhi sel S12 dan S11 di pedesaan (menyumbang gabungan > 63%). Hal ini menegaskan bahwa kualitas air (Layak vs Tidak Layak) sangat bergantung pada jenis sarana air yang digunakan, khususnya bagi masyarakat di wilayah pedesaan.
3.6.2 Analisis Penghapusan Interaksi XZ
Penghapusan efek XZ menghasilkan peningkatan deviance ΔG² = 835,31. Komponen spasial (geografis) memiliki andil besar terhadap higienitas air. Sel S3 (Perkotaan - Terlindungi - Tidak Layak) menyerap 55,10% dampak pemotongan efek XZ. Ini menunjukkan adanya fenomena spesifik di perkotaan, di mana meskipun sarana air minum dikategorikan terlindungi, proporsi air yang berstatus tidak layak berbeda secara drastis dibandingkan wilayah pedesaan.
3.6.3 Analisis Penghapusan Interaksi XY
Meskipun memberikan nilai terkecil (ΔG² = 20,99), interaksi antara wilayah tempat tinggal dan jenis kepemilikan sarana air tetap signifikan. Dampak penghapusan interaksi XY terfokus pada pengguna air tidak terlindungi yang berkualitas layak di perkotaan (S5) dan perdesaan (S11), merefleksikan disparitas struktural aksesibilitas infrastruktur air bersih antar wilayah geografis.
3.7 Interpretasi Parameter dan Odds Ratio Kondisional
Analisis tahap akhir dilakukan dengan menginterpretasikan parameter hubungan asosiasi melalui besaran nilai Conditional Odds Ratio (OR) untuk memahami magnitudo dampak riil di lapangan.
3.7.1 Odds Ratio Kondisional XY
Di bawah model asosiasi homogen [XY][XZ][YZ], parameter interaksi dua faktor dipaksa memiliki nilai odds rasio konstan lintas strata variabel ketiga. Model mengestimasi nilai OR_Fitted_Final sebesar 1,5487 untuk seluruh strata kualitas air.
ORXY = exp(λijXY) = 1,5487
Nilai ini mengindikasikan bahwa secara umum, rumah tangga yang berada di wilayah Perdesaan memiliki kecenderungan 1,5487 kali lebih tinggi untuk menggunakan sarana air tidak terlindungi dibandingkan dengan rumah tangga di perkotaan.
3.7.2 Analisis Odds Ratio Kondisional XZ
Analisis odds rasio bersyarat pada hubungan letak wilayah (Daerah) dengan Kualitas Air (difokuskan pada kategori Aman vs Layak) dipisahkan berdasarkan jenis sarana air minum yang digunakan. Pada kelompok Air Terlindungi, nilai OR = 1,7659 (95% CI: 1,6072; 1,9402) — signifikan karena tidak mencakup angka 1. Ini menunjukkan bahwa pada fasilitas air yang sudah terlindungi, rumah tangga di perdesaan justru memiliki peluang 1,7659 kali lebih besar untuk menikmati air berkategori “Aman” relatif terhadap kualitas “Layak”. Pada kelompok Air Tidak Terlindungi, nilai OR = 2,4167, namun interval kepercayaannya sangat lebar (95% CI: 0,2106; 27,7378) dan mencakup angka 1, sehingga hubungan ini tidak signifikan secara statistik akibat keterbatasan ukuran sampel (sparsity).
3.7.3 Analisis Odds Ratio Kondisional YZ
Bagian terakhir mengukur kontribusi protektif dari fasilitas Jenis Sarana Air Minum terhadap higienitas air (Aman vs Layak) yang dikelompokkan berdasarkan wilayah geografis. Di wilayah Perdesaan, nilai OR = 6,1088 (95% CI: 1,5072; 24,7606) — bernilai signifikan. Ini membuktikan bahwa di pedesaan, pemanfaatan SAM Terlindungi mampu meningkatkan odds perolehan air bersih dengan kualitas “Aman” sebesar 6,11 kali lipat dibandingkan jika menggunakan sarana tidak terlindungi. Di wilayah Perkotaan, nilai OR estimasi bernilai 4,4638, namun interval kepercayaan (95% CI: 0,6032; 33,0311) mencakup nilai 1 (tidak signifikan). Hal ini mengindikasikan bahwa di perkotaan, intervensi pengubahan jenis sarana air minum saja belum tentu menjamin lompatan kualitas air dari layak menjadi aman tanpa adanya perbaikan pada sistem sanitasi lingkungan makro secara terintegrasi.
# ── PERHITUNGAN OR KONDISIONAL DARI FITTED VALUES ─────────────────────────────
# Susun array fitted values: dimensi [Daerah, Jenis_SAM, Kualitas_Air]
mu <- array(fitted(m_homogen), dim = dim(ct),
dimnames = dimnames(ct))
# ── OR XY: Daerah × Jenis SAM, dikontrol pada strata Z = "Aman"
# θXY(k) = (μ[Perkotaan,Terlindungi,k] × μ[Perdesaan,TidakTerlindungi,k]) /
# (μ[Perkotaan,TidakTerlindungi,k] × μ[Perdesaan,Terlindungi,k])
# Dalam model asosiasi homogen, nilai ini KONSTAN di semua k → ambil strata k=1 ("Aman")
or_XY <- (mu["Perkotaan","Air Terlindungi","Aman"] *
mu["Perdesaan","Air Tidak Terlindungi","Aman"]) /
(mu["Perkotaan","Air Tidak Terlindungi","Aman"] *
mu["Perdesaan","Air Terlindungi","Aman"])
# ── OR XZ: Daerah × Kualitas Air (Aman vs Tidak Layak), dikontrol pada strata Y = "Air Terlindungi"
# θXZ(j) = (μ[Perkotaan,j,Aman] × μ[Perdesaan,j,TidakLayak]) /
# (μ[Perkotaan,j,TidakLayak] × μ[Perdesaan,j,Aman])
# Konstan di semua j → ambil strata j=1 ("Air Terlindungi")
or_XZ <- (mu["Perkotaan","Air Terlindungi","Aman"] *
mu["Perdesaan","Air Terlindungi","Tidak Layak"]) /
(mu["Perkotaan","Air Terlindungi","Tidak Layak"] *
mu["Perdesaan","Air Terlindungi","Aman"])
# ── OR YZ: Jenis SAM × Kualitas Air (Aman vs Tidak Layak), dikontrol pada strata X = "Perkotaan"
# θYZ(i) = (μ[i,Terlindungi,Aman] × μ[i,TidakTerlindungi,TidakLayak]) /
# (μ[i,Terlindungi,TidakLayak] × μ[i,TidakTerlindungi,Aman])
# Konstan di semua i → ambil strata i=1 ("Perkotaan")
or_YZ <- (mu["Perkotaan","Air Terlindungi","Aman"] *
mu["Perkotaan","Air Tidak Terlindungi","Tidak Layak"]) /
(mu["Perkotaan","Air Terlindungi","Tidak Layak"] *
mu["Perkotaan","Air Tidak Terlindungi","Aman"])
or_tbl <- data.frame(
Asosiasi = c("XY: Daerah ↔ Jenis SAM",
"XZ: Daerah ↔ Kualitas Air (Aman vs Tidak Layak)",
"YZ: Jenis SAM ↔ Kualitas Air (Terlindungi:Aman vs TT:Tidak Layak)"),
OR = round(c(or_XY, or_XZ, or_YZ), 4),
Interpretasi = c(
paste0("Perdesaan ", round(or_XY, 2), "\u00d7 lebih mungkin menggunakan SAM Tidak Terlindungi"),
paste0("Perdesaan ", round(or_XZ, 2), "\u00d7 lebih mungkin memiliki air Tidak Layak vs Aman"),
paste0("SAM Tidak Terlindungi ", round(or_YZ, 2), "\u00d7 lebih mungkin air Tidak Layak vs Aman")
)
)
kable(or_tbl,
caption = "Tabel 10. Odds Ratio Kondisional Model Asosiasi Homogen [XY][XZ][YZ]",
col.names = c("Pasangan Asosiasi","OR Kondisional","Interpretasi"))| Pasangan Asosiasi | OR Kondisional | Interpretasi |
|---|---|---|
| XY: Daerah ↔︎ Jenis SAM | 1.5487 | Perdesaan 1.55× lebih mungkin menggunakan SAM Tidak Terlindungi |
| XZ: Daerah ↔︎ Kualitas Air (Aman vs Tidak Layak) | 8.7680 | Perdesaan 8.77× lebih mungkin memiliki air Tidak Layak vs Aman |
| YZ: Jenis SAM ↔︎ Kualitas Air (Terlindungi:Aman vs TT:Tidak Layak) | 671.0619 | SAM Tidak Terlindungi 671.06× lebih mungkin air Tidak Layak vs Aman |
671.06× OR YZ: SAM Tidak Terlindungi → Tidak Layak
8.77× OR XZ: Perdesaan → Risiko Air Tidak Layak
1.55× OR XY: Perdesaan → Pengguna SAM Tidak Terlindungi
3.8 Visualisasi Probabilitas Prediksi Model Multinomial
# ── MODEL REGRESI LOGISTIK MULTINOMIAL (PENDUKUNG) ────────────────────────────
model_multi <- multinom(Kualitas_Air ~ Jenis_SAM * Daerah,
weights = Jumlah, data = data_riil,
trace = FALSE)
# Grid prediksi
grid_pred <- expand.grid(
Jenis_SAM = levels(data_riil$Jenis_SAM),
Daerah = levels(data_riil$Daerah)
)
prob_pred <- predict(model_multi, newdata = grid_pred, type = "probs")
df_vis <- cbind(grid_pred, prob_pred) %>%
pivot_longer(cols = c("Aman", "Layak", "Tidak Layak"),
names_to = "Kualitas_Air",
values_to = "Probabilitas") %>%
mutate(
Kualitas_Air = factor(Kualitas_Air,
levels = c("Aman", "Layak", "Tidak Layak")),
Jenis_SAM_lbl = ifelse(Jenis_SAM == "Air Terlindungi",
"Terlindungi", "Tidak Terlindungi")
)
ggplot(df_vis, aes(x = Jenis_SAM_lbl, y = Probabilitas, fill = Kualitas_Air)) +
geom_bar(stat = "identity", position = "dodge",
color = "white", alpha = 0.92, width = 0.68) +
facet_wrap(~Daerah, labeller = labeller(Daerah = c(
"Perkotaan" = "Perkotaan", "Perdesaan" = "Perdesaan"
))) +
scale_fill_manual(
values = c("Aman" = "#0a7c8c",
"Layak" = "#d18b2f",
"Tidak Layak"= "#b84f5a")
) +
scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
geom_text(aes(label = percent(Probabilitas, accuracy = 0.1)),
position = position_dodge(width = 0.68),
vjust = -0.4, size = 3, color = "#1a2e35") +
labs(
title = "Probabilitas Prediksi Status Kualitas Air Rumah Tangga",
subtitle = "Model Multinomial Berbobot — Data Sampel Riil SKAM-RT 2020",
x = "Jenis Sarana Air Minum (SAM)", y = "Probabilitas Prediksi",
fill = "Kualitas Air (Z)"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 13, color = "#0a7c8c"),
plot.subtitle = element_text(color = "#64748b", size = 10),
strip.text = element_text(face = "bold", size = 11, color = "#1a2e35"),
strip.background = element_rect(fill = "#e0f6fa", color = "#17a2b8"),
panel.grid.major.x = element_blank(),
legend.position = "bottom"
)
Gambar 4. Probabilitas Prediksi Status Kualitas Air per Kombinasi Wilayah × Jenis SAM
3.9 Pembahasan: Homogenitas Hubungan Infrastruktur dan Higienitas Air (YZ)
Asosiasi dua arah antara Jenis Sarana Air Minum dan Kualitas Air Bakteriologis (YZ) menunjukkan nilai OR kondisional yang sangat ekstrem, mencapai angka 671,9314. Angka ini menggambarkan kerentanan biologis luar biasa pada rumah tangga yang mengandalkan infrastruktur tidak terlindungi (seperti sumur gali terbuka atau air sungai).
Penemuan bahwa interaksi tingkat tinggi tiga arah [XYZ] tidak signifikan secara statistik memiliki implikasi ilmiah yang sangat mendalam. Hal ini membuktikan bahwa besarnya nilai risiko kontaminasi biologis pada sarana yang tidak terlindungi bersifat homogen, seragam, dan konstan, baik di kawasan perkotaan metropolitan maupun di pelosok perdesaan. Karakteristik fisik sarana air itu sendiri yang menjadi benteng utama perlindungan bakteriologis — bukan letak geografis rumah tangga.
3.10 Pembahasan: Ketimpangan Spasial dan Faktor Kewilayahan (XZ dan XY)
Analisis asosiasi [XZ] dan [XY] mengungkap ketimpangan struktural yang nyata antara perdesaan dan perkotaan di Indonesia:
Asosiasi XY (Spasial–Infrastruktur): Rumah tangga perdesaan memiliki kecenderungan frekuensi 1,8163 kali lebih tinggi untuk menggunakan jenis SAM yang tidak terlindungi dibandingkan kawasan perkotaan. Ini mencerminkan keterbatasan pembangunan infrastruktur air bersih di wilayah perdesaan.
Asosiasi XZ (Spasial–Higienitas): Wilayah perdesaan secara independen berasosiasi kuat dengan keterbatasan higienitas air minum, dengan kecenderungan odds 8,7271 kali lebih tinggi berada pada status “Tidak Layak” dibandingkan “Aman”. Ini mengindikasikan adanya faktor eksternal perdesaan seperti buruknya manajemen sanitasi total berbasis masyarakat (STBM) dan limpasan limbah domestik pertanian.
3.11 Pembahasan: Implikasi Kebijakan terhadap Target SDGs 2030
Temuan penelitian ini memberikan basis ilmiah untuk intervensi berbasis bukti (evidence-based policy). Karena kekuatan asosiasi YZ bersifat konstan di seluruh wilayah, program transformasi SAM tidak terlindungi menjadi sarana komunal terlindungi memiliki urgensi strategis yang sama besarnya di kota maupun di desa. Penuntasan masalah akses air tidak layak di perdesaan akan secara langsung memotong pola asosiasi buruk dalam tabel kontingensi domestik.
4 Kesimpulan
Berdasarkan hasil analisis data kategori tiga dimensi menggunakan pemodelan log-linear terhadap variabel Daerah (X), Jenis Sarana Air Minum (Y), dan Kualitas Air (Z), diperoleh beberapa kesimpulan utama sebagai berikut:
Struktur Model Terbaik dan Interaksi Tingkat Tinggi: Model jenuh (saturated model) [XYZ] merupakan model terbaik berdasarkan kriteria informasi minimum (AIC = 111,477; BIC = 117,296). Uji selisih deviance formal mengonfirmasi bahwa interaksi tiga faktor [XYZ] signifikan secara statistik (ΔG² = 18,9199; p = 0,0001). Hal ini membuktikan secara empiris bahwa pola hubungan dua faktor bersifat heterogen, di mana odds rasio antar-variabel berubah secara dinamis bergantung pada tingkatan (strata) variabel ketiga.
Validitas Komputasi Algoritma: Pendekatan Iterative Proportional Fitting (IPF) dan Generalized Linear Model (GLM) menunjukkan kesetaraan matematis (mathematical equivalence) yang mutlak dalam mengestimasi nilai harapan frekuensi sel model log-linear non-saturasi, ditunjukkan oleh nilai selisih estimasi yang tepat sebesar nol. Algoritma IPF terbukti efisien dengan mencapai konvergensi pada iterasi ke-10 dengan batas toleransi kesalahan (maximum error margin) yang sangat ketat sebesar 0,00553.
Pemicu Ketidakcocokan Model Homogen (Local Outliers): Penolakan global terhadap model asosiasi homogen [XY][XZ][YZ] disebabkan oleh anomali lokal pada sel-sel spesifik. Sel S5 (Perkotaan - Air Tidak Terlindungi - Layak) menjadi komponen paling disruptif dengan menyumbang 53,83% terhadap total sisa deviance global, diikuti oleh sel S11 (Perdesaan - Air Tidak Terlindungi - Layak) sebesar 18,91%. Kedua sel pencilan ini mencerminkan adanya fluktuasi ekstrem pada kelompok pengguna air tidak terlindungi yang tidak mampu diakomodasi oleh asumsi odds rasio konstan.
Hierarki Kekuatan Asosiasi Parsial: Uji asosiasi parsial menunjukkan bahwa interaksi antara Jenis Sarana Air Minum dan Kualitas Air (YZ) memiliki kontribusi paling dominan dalam sistem variabel dengan peningkatan deviance terbesar saat dihapus (ΔG² = 5.904,40). Interaksi Daerah dan Kualitas Air (XZ) menempati posisi dominan kedua (ΔG² = 835,31), sedangkan interaksi Daerah dan Jenis Sarana Air Minum (XY) memiliki pengaruh terkecil namun tetap signifikan (ΔG² = 20,99).
Magnitudo Dampak melalui Rasio Odds Kondisional:
- Disparitas Spasial Akses Air: Rumah tangga di wilayah perdesaan secara umum memiliki kecenderungan 1,55 kali lebih tinggi untuk menggunakan sarana air minum tidak terlindungi dibandingkan dengan rumah tangga di perkotaan.
- Efektivitas Sarana Terlindungi di Perdesaan: Di wilayah perdesaan, pemanfaatan Jenis Sarana Air Minum Terlindungi secara signifikan bertindak sebagai faktor protektif, meningkatkan odds perolehan air bersih berkategori “Aman” sebesar 6,11 kali lipat (95% CI: 1,5072; 24,7606) dibandingkan dengan sarana tidak terlindungi.
- Tantangan Sanitasi Perkotaan: Di wilayah perkotaan, pengaruh jenis sarana air minum terhadap peningkatan kualitas air dari layak menjadi aman tidak signifikan secara statistik (95% CI mencakup nilai 1). Hal ini mengindikasikan bahwa intervensi infrastruktur sarana fisik saja di perkotaan belum menjamin keamanan kualitas air tanpa adanya mitigasi terhadap pencemaran lingkungan makro.
Rekomendasi Kebijakan: Kementerian Kesehatan RI dan Kementerian PUPR perlu mengintegrasikan perluasan jaringan air perpipaan dan revitalisasi sumur terlindungi dengan memprioritaskan infrastruktur di wilayah perdesaan guna memutus rantai asosiasi buruk demi tercapainya target SDGs 6.1.1 tahun 2030.
Referensi
Agresti, A. An Introduction to Categorical Data Analysis, 2nd ed.; John Wiley & Sons, Inc.: Hoboken, New Jersey, 2007; pp. 200–245.
Dobson, A.J. An Introduction to Statistical Modelling; Chapman and Hall: London, 1983; pp. 95–110.
Irianto, J.; dkk. Studi Kualitas Air Minum Rumah Tangga (SKAM-RT) di Indonesia 2020; Badan Penelitian dan Pengembangan Kesehatan (Balitbangkes), Kementerian Kesehatan RI: Jakarta, 2020; pp. 45–80.
Fleiss, J.L.; Levin, B.; Paik, M.C. Statistical Methods for Rates and Proportions, 3rd ed.; John Wiley & Sons, Inc.: Hoboken, New Jersey, 2003; pp. 112–135.
Hosmer, D.W.; Lemeshow, S.; Sturdivant, R.X. Applied Logistic Regression, 3rd ed.; John Wiley & Sons, Inc.: Hoboken, New Jersey, 2013; pp. 210–240.
World Health Organization. Drinking-water Fact Sheet; World Health Organization: Geneva, 2023. Available online: https://www.who.int (accessed on Day Month Year).
Salamah; Ayuningrum. Analisis Faktor Sanitasi dan Sumber Air Minum yang Mempengaruhi Insiden Diare pada Balita di Jawa Timur. Jurnal Sains dan Seni ITS 2015, 4, 233–240.
Fikri; dkk. Kualitas Air Minum Rumah Tangga di Indonesia Berdasarkan Parameter Fisik, Kimia, dan Mikrobiologi. Jurnal Penelitian Inovatif 2025, 5, 112–125.
World Health Organization; United Nations Children’s Fund. Progress on Household Drinking Water, Sanitation and Hygiene 2000-2022: Special Focus on Gender; World Health Organization and United Nations Children’s Fund: New York, 2023; pp. 12–35.
Santos, T.M.; et al. E. coli contamination of drinking water sources in rural and urban settings: an analysis of 38 nationally representative household surveys (2014–2021). Journal of Water and Health 2023, 21, 485–502.