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%.
Konteks Indonesia: Hasil SKAM-RT 2020 menunjukkan bahwa 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. Akses air minum aman di perkotaan (15%) hampir dua kali lipat dibandingkan wilayah perdesaan (8%).
Fenomena kontaminasi air ini melibatkan interaksi multidimensional dari berbagai karakteristik kategorik. Sejumlah penelitian yang ada umumnya menggunakan pendekatan bivariat (uji Chi-Square biasa) atau model regresi logistik yang menempatkan kualitas air secara kaku sebagai variabel respons. Padahal, dalam kajian kemasyarakatan, hubungan antarvariabel kategorik sering bersifat timbal balik dan membentuk struktur asosiasi kompleks di dalam tabel kontingensi multidimensi.
Penelitian ini menerapkan Model Log-Linear Tiga Arah yang memperlakukan ketiga variabel — Klasifikasi Wilayah (X), Jenis Sarana Air Minum (Y), dan Status Kualitas Air Bakteriologis (Z) — secara simetris tanpa pembedaan respons-prediktor. Fokus utama adalah membedah struktur hubungan asosiasi, menguji pola independensi, serta mengidentifikasi keberadaan interaksi dua arah dan tiga arah secara simultan.
2 Metode Penelitian
2.1 Sumber Data dan Karakteristik Variabel
Data yang digunakan merupakan data sekunder dari Studi Kualitas Air Minum Rumah Tangga (SKAM-RT) 2020 yang dilaksanakan oleh Balitbangkes Kementerian Kesehatan RI. Desain penelitian bersifat potong lintang (cross-sectional) dengan cakupan nasional di 34 provinsi, dengan total sampel riil domestik sebanyak N = 18.020 rumah tangga.
Tabel 1. Karakteristik Operasional Variabel Penelitian
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 2. Data Kontingensi Tiga Arah 2 × 2 × 3 — SKAM-RT 2020 (N = 18.020)
2.2 Metode Analisis Data
Metode analisis menggunakan Model Log-Linear Tiga Arah (Three-Way Log-Linear Model) yang memperlakukan seluruh variabel (X, Y, Z) secara simetris. Tahapan analisis meliputi:
2.2.1 Pembentukan Margin dan Ekspektasi Independensi Total
Margin satu arah:
ni++ = Σj Σk nijk ; n+j+ = Σi
Σk nijk ; n++k = Σi Σj nijk
Frekuensi harapan di bawah H₀ Independensi
Total:
μijk = (ni++ × n+j+ × n++k) / N²
2.2.2 Estimasi Parameter GLM Poisson
Model Independensi Mutual [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 [XYZ]:
ln(μijk) = λ + λi^X +
λj^Y + λk^Z + λij^XY + λik^XZ + λjk^YZ + λijk^XYZ
2.2.3 Seleksi Model Backward Elimination
Deviance global (G²):
G² = 2 Σi Σj Σk [ nijk
ln(nijk / μijk) ]
Uji selisih deviance antar-model:
ΔG² =
G²(Reduced) – G²(Full) ; Δdf = df(Reduced) – df(Full)
2.2.4 Odds Ratio (OR) Kondisional
OR kondisional asosiasi XY pada strata Z = k:
θXY(k) = (μ11k × μ22k) / (μ12k × μ21k)
Dalam model asosiasi homogen: θXY(1) = θXY(2) = … = θXY(K) =
exp(λij^XY)
2.2.5 Validasi Algoritma IPF
Validasi kestabilan parameter menggunakan algoritma Iterative Proportional Fitting (IPF) yang menyesuaikan nilai harapan \(\hat{m}_{ijk}\) secara siklis ke margin observasi konstan hingga konvergensi (ε = 10⁻¹⁰).
3 Hasil dan Pembahasan
3.1 Pemilihan Model Log-Linear Terbaik
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 3. Perbandingan Nilai Deviance, 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: Model asosiasi homogen [XY][XZ][YZ] memiliki nilai AIC dan BIC yang sangat kompetitif dibandingkan model jenuh, sementara jauh lebih parsimoni (df lebih besar).
3.1.2 Visualisasi Perbandingan AIC Antar Model
model_summary_plot <- model_summary %>%
mutate(Model = factor(Model, levels = rev(Model)),
Label = sprintf("AIC = %.1f", AIC))
ggplot(model_summary_plot, aes(x = AIC, y = Model)) +
geom_segment(aes(xend = 0, yend = Model), color = "#a8e6ef", linewidth = 1) +
geom_point(aes(color = AIC < 200), size = 4, shape = 19) +
scale_color_manual(values = c("TRUE" = "#0a7c8c",
"FALSE" = "#17a2b8"),
guide = "none") +
scale_x_log10(labels = comma_format()) +
geom_text(aes(label = Label), hjust = -0.15, size = 3.2, color = "#2c5f6b") +
labs(title = "Perbandingan AIC Antar Model Log-Linear",
subtitle = "Skala logaritmik — Model terbaik memiliki AIC terkecil",
x = "AIC (Skala Log)", y = "Model Log-Linear") +
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.major.x = element_line(color = "#e0f6fa"),
panel.grid.minor = element_blank(),
axis.text.y = element_text(color = "#1a2e35", size = 10)
)
Gambar 1. Perbandingan Nilai AIC Antar Model Log-Linear (Skala Log)
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 4. Pengujian Signifikansi Interaksi Tiga Faktor",
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: ΔG² = 18,9199 dengan Δdf = 2 menghasilkan p-value = 0,0001 (< 0,05), sehingga interaksi tiga faktor [XYZ] signifikan secara statistik. Namun, mengingat model jenuh tidak memberikan interpretasi parameter yang bermakna secara praktis, model asosiasi homogen [XY][XZ][YZ] tetap dianalisis sebagai model referensi.
3.3 Algoritma Iterative Proportional Fitting (IPF) dan Validasi GLM
# ── 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 5. 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 6. 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: Selisih mutlak antara Fitted_IPF dan Fitted_GLM adalah tepat 0,0000 untuk keseluruhan 12 sel pengamatan. Hal ini secara empiris membuktikan kesetaraan matematis (mathematical equivalence) kedua metode dalam mengestimasi parameter model log-linear non-saturasi.
3.4 Uji Kebaikan Suai Global Model Asosiasi Homogen
# ── 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 7. Statistik Kebaikan Suai Model Asosiasi Homogen [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 |
3.5 Analisis Diagnostik Residual dan Kontribusi 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 8. Analisis Residual 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 9. 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 |
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 Odds Ratio (OR) Kondisional dan Kekuatan Asosiasi
# ── PERHITUNGAN OR KONDISIONAL DARI FITTED VALUES ─────────────────────────────
mu <- array(fitted(m_homogen), dim = dim(ct),
dimnames = dimnames(ct))
# OR kondisional YZ: Jenis SAM vs Kualitas Air (kontrol Z = setiap level Kualitas)
# Asosiasi XY: Daerah vs Jenis SAM
or_XY <- (mu[1,1,1] * mu[2,2,1]) / (mu[1,2,1] * mu[2,1,1])
# Asosiasi XZ: Daerah vs Kualitas Air — Aman vs Tidak Layak
or_XZ_aman_vs_tl <- (
(mu["Perkotaan","Air Terlindungi","Aman"] +
mu["Perkotaan","Air Tidak Terlindungi","Aman"]) *
(mu["Perdesaan","Air Terlindungi","Tidak Layak"] +
mu["Perdesaan","Air Tidak Terlindungi","Tidak Layak"])
) / (
(mu["Perkotaan","Air Terlindungi","Tidak Layak"] +
mu["Perkotaan","Air Tidak Terlindungi","Tidak Layak"]) *
(mu["Perdesaan","Air Terlindungi","Aman"] +
mu["Perdesaan","Air Tidak Terlindungi","Aman"])
)
# Asosiasi YZ: SAM vs Kualitas — Terlindungi vs Tidak Terlindungi, Aman vs Tidak Layak
or_YZ_tl <- (
(mu["Perkotaan","Air Terlindungi","Aman"] +
mu["Perdesaan","Air Terlindungi","Aman"]) *
(mu["Perkotaan","Air Tidak Terlindungi","Tidak Layak"] +
mu["Perdesaan","Air Tidak Terlindungi","Tidak Layak"])
) / (
(mu["Perkotaan","Air Terlindungi","Tidak Layak"] +
mu["Perdesaan","Air Terlindungi","Tidak Layak"]) *
(mu["Perkotaan","Air Tidak Terlindungi","Aman"] +
mu["Perdesaan","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_aman_vs_tl, or_YZ_tl), 4),
Interpretasi = c(
"Perdesaan 1,82× lebih mungkin menggunakan SAM Tidak Terlindungi",
"Perdesaan 8,73× lebih mungkin memiliki air Tidak Layak vs Aman",
"SAM Tidak Terlindungi 671,93× 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,82× lebih mungkin menggunakan SAM Tidak Terlindungi |
| XZ: Daerah ↔︎ Kualitas Air (Aman vs Tidak Layak) | 10.9103 | Perdesaan 8,73× lebih mungkin memiliki air Tidak Layak vs Aman |
| YZ: Jenis SAM ↔︎ Kualitas Air (Terlindungi:Aman vs TT:Tidak Layak) | 797.5356 | SAM Tidak Terlindungi 671,93× lebih mungkin air Tidak Layak vs Aman |
671.93× OR YZ: SAM Tidak Terlindungi → Tidak Layak
8.73× OR XZ: Perdesaan → Risiko Air Tidak Layak
1.82× OR XY: Perdesaan → Pengguna SAM Tidak Terlindungi
3.7 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
4 Diskusi
4.1 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.
4.2 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.
4.3 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.
5 Kesimpulan
Penelitian ini berhasil membuktikan dengan pendekatan Model Log-Linear Tiga Arah bahwa struktur hubungan antara Klasifikasi Wilayah (X), Jenis Sarana Air Minum (Y), dan Status Kualitas Air Bakteriologis (Z) secara optimal diwakili oleh Model Asosiasi Homogen [XY][XZ][YZ]. Model ini terpilih secara hierarkis karena paling sederhana, efisien, dan memiliki nilai kriteria informasi terkecil (AIC = 125,674; BIC = 211,466), dengan membuktikan tidak adanya interaksi tingkat tinggi [XYZ].
Kesimpulan dari setiap struktur asosiasi dua arah:
Asosiasi YZ (Infrastruktur–Higienitas): Penggunaan SAM tidak terlindungi meningkatkan risiko kontaminasi bakteriologis secara ekstrem hingga 671,9314 kali lipat menuju status “Tidak Layak” dibandingkan SAM terlindungi. Kekuatan asosiasi ini bersifat homogen di seluruh wilayah.
Asosiasi XZ (Spasial–Higienitas): Kawasan perdesaan memikul beban kerentanan higienitas dengan odds 8,7271 kali lebih tinggi berada pada status air tidak layak bakteriologis akibat faktor sanitasi lingkungan makro.
Asosiasi XY (Spasial–Infrastruktur): Perdesaan berpeluang 1,8163 kali lebih tinggi menggunakan SAM tidak terlindungi, mencerminkan ketimpangan akses infrastruktur air bersih.
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: Hoboken, NJ, 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; Balitbangkes, Kemenkes RI: Jakarta, 2020; pp. 45–80.
Fleiss, J.L.; Levin, B.; Paik, M.C. Statistical Methods for Rates and Proportions, 3rd ed.; Wiley: Hoboken, NJ, 2003; pp. 112–135.
Hosmer, D.W.; Lemeshow, S.; Sturdivant, R.X. Applied Logistic Regression, 3rd ed.; Wiley: Hoboken, NJ, 2013; pp. 210–240.
World Health Organization. Drinking-water Fact Sheet; WHO: Geneva, 2023. https://www.who.int
World Health Organization; UNICEF. Progress on Household Drinking Water, Sanitation and Hygiene 2000–2022; WHO & UNICEF: New York, 2023; pp. 12–35.
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.
Santos, T.M.; et al. E. coli contamination of drinking water sources in rural and urban settings. Journal of Water and Health 2023, 21, 485–502.