Pemodelan Regresi Logistik Biner untuk Identifikasi Faktor Risiko dan Prediksi Kontaminasi Lahan Industri di China
Salsabila Najwa Effianti & Fathia Salwadifina
Program Studi Statistika, Fakultas MIPA, Universitas Padjadjaran
2026-06-10
required_packages <- c("dplyr", "ggplot2", "broom", "knitr", "scales",
"readxl", "ResourceSelection", "DescTools", "car",
"tidyr", "kableExtra")
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0) install.packages(missing_packages)
invisible(lapply(required_packages, library, character.only = TRUE))1 Pendahuluan
Latar Belakang: Pencemaran lahan industri merupakan salah satu permasalahan lingkungan yang serius, terutama di negara berkembang yang mengalami urbanisasi pesat seperti China. Jumlah lahan terkontaminasi di China meningkat hingga 325% antara tahun 1990 dan 2018, dengan 43.676 lokasi terkontaminasi teridentifikasi dari 83.498 lahan perusahaan.
Tujuan penelitian:
- Mengidentifikasi faktor risiko yang secara signifikan memengaruhi kontaminasi lahan industri.
- Membangun model prediktif berbasis regresi logistik biner yang parsimonious dengan validasi asumsi komprehensif.
- Menginterpretasikan besar dan arah pengaruh masing-masing faktor melalui nilai Odds Ratio (OR).
2 Landasan Metodologi
2.1 Model Regresi Logistik Biner
\[g(x) = \ln\!\left(\frac{\pi(x)}{1-\pi(x)}\right) = \alpha + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p \tag{1}\]
\[\pi(x) = \frac{e^{g(x)}}{1 + e^{g(x)}} \tag{2}\]
2.2 Odds Ratio
\[OR = \exp(\hat{\beta}_p) \tag{3}\]
OR > 1 → meningkatkan peluang kontaminasi | OR < 1 → efek protektif
2.3 Likelihood Ratio Test (LRT)
H0: Variabel yang dihapus tidak berkontribusi signifikan terhadap model
H1: Variabel yang dihapus berkontribusi signifikan terhadap model
\[G = -2\ln\!\left(\frac{L_{\text{without variable}}}{L_{\text{with variable}}}\right) \sim \chi^2_{df} \tag{4}\]
2.4 Uji Multikolinearitas — VIF
H0: Tidak terdapat masalah multikolinearitas antar variabel independen
H1: Terdapat masalah multikolinearitas antar variabel independen
\[VIF_j = \frac{1}{1 - R_j^2} \tag{5}\]
Ambang batas: GVIF¹/(²ᵈᶠ) > √10 ≈ 3,16 (kategorik) | VIF > 10 (numerik)
2.5 Uji Linearitas Logit — Box-Tidwell
H0: γ₁ = 0 (asumsi linearitas logit terpenuhi)
H1: γ₁ ≠ 0 (asumsi linearitas logit dilanggar)
\[\ln\!\left(\frac{\pi(x)}{1-\pi(x)}\right) = \alpha + \beta_1 X_1 + \gamma_1 \bigl[X_1 \cdot \ln(X_1)\bigr] \tag{6}\]
2.6 Uji Kecocokan Model — Hosmer-Lemeshow
\[C = \sum_{g=1}^{G} \frac{(O_g - n_g \hat{\pi}_g)^2}{n_g \hat{\pi}_g (1 - \hat{\pi}_g)} \sim \chi^2_{G-2} \tag{7}\]
Model fit apabila p-value > 0,05 | G = 10 kelompok | df = G − 2 = 8
2.7 Indeks Youden — Threshold Optimal
\[J = \text{Sensitivitas} + \text{Spesifisitas} - 1 \tag{8}\]
2.8 Metrik Evaluasi Klasifikasi
\[\text{Akurasi} = \frac{TP+TN}{TP+TN+FP+FN} \tag{9}\]
\[\text{Sensitivitas} = \frac{TP}{TP+FN} \tag{10}\]
\[\text{Spesifisitas} = \frac{TN}{TN+FP} \tag{11}\]
\[\text{Presisi} = \frac{TP}{TP+FP} \tag{12}\]
\[F1\text{-Score} = 2 \times \frac{\text{Presisi} \times \text{Sensitivitas}}{\text{Presisi}+\text{Sensitivitas}} \tag{13}\]
\[\text{Balanced Accuracy} = \frac{\text{Sensitivitas}+\text{Spesifisitas}}{2} \tag{14}\]
2.9 Nagelkerke Pseudo-\(R^2\)
Nagelkerke Pseudo-\(R^2\) digunakan untuk mengukur kemampuan model regresi logistik dalam menjelaskan variasi pada variabel respons.
\[R^2_{\text{Nagelkerke}} = \frac{1 - \left(\dfrac{L_0}{L_M}\right)^{2/n}}{1 - (L_0)^{2/n}} \tag{15}\]
Dengan:
\(L_0\) : Nilai likelihood dari model tanpa variabel prediktor.
\(L_M\) : Nilai likelihood dari model yang sedang diuji.
\(n\) : Total ukuran sampel data penelitian.
3 Data dan Prapemrosesan
3.1 Sumber Data
Dataset yang digunakan adalah Urban Contaminated Sites along
China’s Urbanization yang dipublikasikan oleh Li (2023) melalui
Mendeley Data (DOI: 10.17632/r4y2vcpfmx.1). Dataset ini
mencakup 2.005 lokasi industri dengan 15
variabel (1 respon + 14 prediktor).
3.2 Baca Data
# Ganti path sesuai lokasi file kamu
raw_contamination <- readxl::read_excel("C:/Users/Salsabila Najwa/Downloads/Input dataset.xlsx")
names(raw_contamination) <- trimws(names(raw_contamination))
ringkasan_data <- data.frame(
Keterangan = c("Jumlah observasi (lokasi industri)", "Jumlah variabel"),
Nilai = c(nrow(raw_contamination), ncol(raw_contamination))
)
knitr::kable(ringkasan_data,
caption = "Ukuran Dataset Kontaminasi Lahan Perkotaan China") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE)| Keterangan | Nilai |
|---|---|
| Jumlah observasi (lokasi industri) | 2005 |
| Jumlah variabel | 16 |
3.3 Kamus Variabel
kamus_variabel <- data.frame(
`Kolom Asli` = c(
"Presence of contamination", "Duration", "Starting time",
"Industry class", "Violations", "Precipitation (mm)",
"Temperature", "Impervious surfaces", "Scale",
"Pollutant mobility", "Pollutant volatility",
"Persistent organic pollutants",
"Soil texture (Soil clay content in percentage)",
"Soil erosion", "Wind speed (m/s)"
),
`Variabel Analisis` = c(
"status_kontaminasi", "durasi_operasi", "tahun_mulai",
"jenis_industri", "jumlah_pelanggaran", "curah_hujan",
"suhu", "permukaan_kedap", "skala_industri",
"mobilitas_polutan", "volatilitas_polutan",
"polutan_organik_persisten",
"tekstur_tanah", "erosi_tanah", "kecepatan_angin"
),
`Keterangan` = c(
"Status kontaminasi lahan (RESPON Y)", "Lama operasi industri (tahun)",
"Tahun mulai beroperasi", "Jenis/kelas industri",
"Jumlah pelanggaran lingkungan", "Total curah hujan tahunan (mm)",
"Rata-rata suhu tahunan (°C)", "Ada/tidaknya permukaan kedap air",
"Skala ukuran industri", "Tingkat mobilitas polutan dalam tanah",
"Tingkat volatilitas polutan", "Ada/tidaknya polutan organik persisten",
"Kandungan liat tanah (%)", "Tingkat erosi tanah",
"Rata-rata kecepatan angin (m/s)"
),
`Tipe` = c(
"Respon biner", "Numerik", "Numerik",
"Kategorik (banyak level)", "Numerik", "Numerik",
"Numerik", "Kategorik (biner)", "Kategorik (ordinal)",
"Numerik", "Numerik", "Kategorik (biner)",
"Numerik", "Kategorik (ordinal)", "Numerik"
),
check.names = FALSE
)
knitr::kable(kamus_variabel,
caption = "Kamus Variabel Dataset Kontaminasi Lahan") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE)| Kolom Asli | Variabel Analisis | Keterangan | Tipe |
|---|---|---|---|
| Presence of contamination | status_kontaminasi | Status kontaminasi lahan (RESPON Y) | Respon biner |
| Duration | durasi_operasi | Lama operasi industri (tahun) | Numerik |
| Starting time | tahun_mulai | Tahun mulai beroperasi | Numerik |
| Industry class | jenis_industri | Jenis/kelas industri | Kategorik (banyak level) |
| Violations | jumlah_pelanggaran | Jumlah pelanggaran lingkungan | Numerik |
| Precipitation (mm) | curah_hujan | Total curah hujan tahunan (mm) | Numerik |
| Temperature | suhu | Rata-rata suhu tahunan (°C) | Numerik |
| Impervious surfaces | permukaan_kedap | Ada/tidaknya permukaan kedap air | Kategorik (biner) |
| Scale | skala_industri | Skala ukuran industri | Kategorik (ordinal) |
| Pollutant mobility | mobilitas_polutan | Tingkat mobilitas polutan dalam tanah | Numerik |
| Pollutant volatility | volatilitas_polutan | Tingkat volatilitas polutan | Numerik |
| Persistent organic pollutants | polutan_organik_persisten | Ada/tidaknya polutan organik persisten | Kategorik (biner) |
| Soil texture (Soil clay content in percentage) | tekstur_tanah | Kandungan liat tanah (%) | Numerik |
| Soil erosion | erosi_tanah | Tingkat erosi tanah | Kategorik (ordinal) |
| Wind speed (m/s) | kecepatan_angin | Rata-rata kecepatan angin (m/s) | Numerik |
3.4 Preprocessing & Rekoding
Catatan preprocessing
jenis_industri: Penggabungan kategori dilakukan menggunakan ambang batas frekuensi minimum sebesar 31 observasi. Sebanyak 18 jenis industri dengan frekuensi ≥ 31 dipertahankan sebagai kategori tersendiri, sedangkan 42 jenis industri dengan frekuensi < 31 digabungkan ke dalam kategori “Industri Lainnya” sehingga menghasilkan 427 observasi pada kategori tersebut.
collapse_rare <- function(x, min_n = 31, other_label = "Industri Lainnya") {
x_chr <- as.character(x)
keep <- names(which(table(x_chr) >= min_n))
factor(ifelse(x_chr %in% keep, x_chr, other_label))
}
erosi_levels <- c("Slight", "Mild", "Moderate", "Strong", "Very strong", "Severe")
contamination <- raw_contamination %>%
transmute(
status_kontaminasi = factor(`Presence of contamination`, levels = c("No", "Yes")),
kontaminasi_biner = as.integer(`Presence of contamination` == "Yes"),
durasi_operasi = as.numeric(Duration),
tahun_mulai = as.numeric(`Starting time`),
jumlah_pelanggaran = as.numeric(Violations),
curah_hujan = as.numeric(`Precipitation (mm)`),
suhu = as.numeric(`Temperature (℃）`),
mobilitas_polutan = as.numeric(`Pollutant mobility`),
volatilitas_polutan = as.numeric(`Pollutant volatility`),
tekstur_tanah = as.numeric(`Soil texture (Soil clay content in percentage)`),
kecepatan_angin = as.numeric(`Wind speed (m/s)`),
log_mobilitas_polutan = log(as.numeric(`Pollutant mobility`) + 1),
permukaan_kedap = factor(`Impervious surfaces`,
levels = c("No","Yes"), labels = c("Tidak","Ya")),
polutan_organik_persisten = factor(`Persistent organic pollutants`,
levels = c("No","Yes"), labels = c("Tidak","Ya")),
skala_industri = factor(Scale,
levels = c("Micro","Small","Medium","Large","Extra-large"),
labels = c("Mikro","Kecil","Menengah","Besar","Sangat Besar")),
erosi_tanah = factor(`Soil erosion`, levels = erosi_levels,
labels = c("Ringan","Agak Ringan","Sedang",
"Kuat","Sangat Kuat","Parah")),
jenis_industri = {
x_chr <- as.character(`Industry class`)
freq <- table(x_chr)
keep <- names(freq[freq >= 31])
factor(ifelse(x_chr %in% keep, x_chr, "Industri Lainnya"))
}
) %>%
droplevels() %>%
mutate(jenis_industri = relevel(jenis_industri, ref = "Industri Lainnya"))3.4.1 Verifikasi
Preprocessing jenis_industri
freq_industri <- sort(table(contamination$jenis_industri), decreasing = TRUE)
freq_df <- as.data.frame(freq_industri)
names(freq_df) <- c("Jenis Industri", "Frekuensi")
knitr::kable(freq_df,
caption = "Frekuensi Jenis Industri Setelah Penggabungan") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE)| Jenis Industri | Frekuensi |
|---|---|
| Industri Lainnya | 427 |
| Metal surface treatment and heat treatment processing | 209 |
| Organic chemical raw material manufacturing | 134 |
| Cotton printing and dyeing processing | 133 |
| Chemical reagent and additive manufacturing | 129 |
| Pesticide manufacturing | 126 |
| Leather tanning and processing | 121 |
| Coking | 104 |
| Crude oil processing and petroleum product manufacturing | 100 |
| Iron ore mining and dressing | 84 |
| Hazardous waste treatment | 67 |
| Lead-zinc smelting | 52 |
| Petroleum extraction | 52 |
| Iron smelting | 50 |
| Steel smelting | 49 |
| Chemical pharmaceutical raw material manufacturing | 48 |
| Lead-zinc ore mining and dressing | 47 |
| Gold ore mining and dressing | 40 |
| Copper smelting | 33 |
## Baseline jenis_industri: Industri Lainnya4 Eksplorasi Data
4.1 Distribusi Variabel Respon
distribusi_y <- contamination %>%
count(status_kontaminasi, name = "Jumlah") %>%
mutate(Proporsi = scales::percent(Jumlah / sum(Jumlah), accuracy = 0.1)) %>%
rename("Status Kontaminasi" = status_kontaminasi)
knitr::kable(distribusi_y,
caption = "Distribusi Kelas Respon (Y)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"),
full_width = FALSE)| Status Kontaminasi | Jumlah | Proporsi |
|---|---|---|
| No | 801 | 40.0% |
| Yes | 1204 | 60.0% |
ggplot(contamination, aes(x = status_kontaminasi, fill = status_kontaminasi)) +
geom_bar(width = 0.62, color = "white", linewidth = 0.8) +
geom_text(stat = "count", aes(label = after_stat(count)),
vjust = -0.4, fontface = "bold") +
scale_fill_manual(values = c("No" = "#2d6a2d", "Yes" = "#eab308")) +
labs(title = "Distribusi Status Kontaminasi Lahan",
x = NULL, y = "Jumlah Lokasi") +
theme_minimal(base_size = 12) +
theme(legend.position = "none")
Distribusi variabel respon: Dataset terdiri atas 1.204 lokasi terkontaminasi (60,0%) dan 801 lokasi tidak terkontaminasi (40,0%). Distribusi kelas menunjukkan ketidakseimbangan yang relatif moderat dengan selisih proporsi sebesar 20 poin persentase. Oleh karena itu, evaluasi model perlu mempertimbangkan kondisi tersebut sehingga balanced accuracy digunakan sebagai metrik yang lebih representatif dibandingkan akurasi biasa.
ringkasan_numerik <- contamination %>%
group_by(status_kontaminasi) %>%
summarise(
Jumlah = n(),
`Rata-rata Durasi (th)` = round(mean(durasi_operasi, na.rm = TRUE), 2),
`Tahun Mulai (rata-rata)` = round(mean(tahun_mulai, na.rm = TRUE), 2),
`Rata-rata Pelanggaran` = round(mean(jumlah_pelanggaran, na.rm = TRUE), 2),
`Median Curah Hujan (mm)` = round(median(curah_hujan, na.rm = TRUE), 2),
`Rata-rata Suhu (°C)` = round(mean(suhu, na.rm = TRUE), 2),
`Rata-rata Mobilitas Pol.` = round(mean(mobilitas_polutan, na.rm = TRUE), 2),
`Rata-rata Volatilitas Pol.` = round(mean(volatilitas_polutan, na.rm = TRUE), 2),
`Rata-rata Tekstur Tanah` = round(mean(tekstur_tanah, na.rm = TRUE), 2),
`Rata-rata Kec. Angin (m/s)` = round(mean(kecepatan_angin, na.rm = TRUE), 2),
.groups = "drop"
) %>%
rename("Status Kontaminasi" = status_kontaminasi)
ringkasan_numerik %>%
tidyr::pivot_longer(-`Status Kontaminasi`, names_to = "Statistik", values_to = "Nilai") %>%
tidyr::pivot_wider(names_from = `Status Kontaminasi`, values_from = Nilai) %>%
knitr::kable(caption = "Ringkasan Numerik per Kelas Respon") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Statistik | No | Yes |
|---|---|---|
| Jumlah | 801.00 | 1204.00 |
| Rata-rata Durasi (th) | 13.36 | 24.92 |
| Tahun Mulai (rata-rata) | 2004.07 | 1992.11 |
| Rata-rata Pelanggaran | 1.00 | 2.34 |
| Median Curah Hujan (mm) | 8584.68 | 12222.13 |
| Rata-rata Suhu (°C) | 13.55 | 14.26 |
| Rata-rata Mobilitas Pol. | 0.59 | 0.60 |
| Rata-rata Volatilitas Pol. | 0.57 | 0.64 |
| Rata-rata Tekstur Tanah | 26.11 | 26.14 |
| Rata-rata Kec. Angin (m/s) | 21.49 | 22.13 |
4.2 Visualisasi Eksplorasi
# Scatter: durasi vs curah hujan
ggplot(contamination, aes(x = durasi_operasi, y = curah_hujan,
color = status_kontaminasi)) +
geom_point(alpha = 0.45, size = 1.8) +
scale_color_manual(values = c("No" = "#2d6a2d", "Yes" = "#eab308")) +
scale_y_continuous(labels = scales::comma) +
labs(title = "Durasi Operasi vs Curah Hujan",
subtitle = "Warna = status kontaminasi lahan",
x = "Durasi operasi (tahun)", y = "Curah hujan tahunan (mm)",
color = "Status Kontaminasi") +
theme_minimal(base_size = 12)
Scatter plot — Durasi Operasi vs Curah Hujan: Lokasi terkontaminasi (kuning) cenderung terkonsentrasi pada durasi operasi yang lebih panjang (>20 tahun) dan pada curah hujan lebih tinggi, sementara lokasi bersih (hijau) lebih banyak di bawah 20 tahun operasi. Pola ini memberikan indikasi awal bahwa durasi operasi dan curah hujan berpotensi menjadi prediktor signifikan.
4.3 Contoh 8 Data Pertama
contamination %>%
select(status_kontaminasi, jenis_industri, durasi_operasi,
tahun_mulai, skala_industri, permukaan_kedap, erosi_tanah,
curah_hujan, suhu) %>%
head(8) %>%
knitr::kable(caption = "Contoh 8 Data Pertama") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| status_kontaminasi | jenis_industri | durasi_operasi | tahun_mulai | skala_industri | permukaan_kedap | erosi_tanah | curah_hujan | suhu |
|---|---|---|---|---|---|---|---|---|
| Yes | Petroleum extraction | 21 | 1992 | Mikro | Tidak | Ringan | 13775.359 | 19.7 |
| No | Petroleum extraction | 16 | 2000 | Mikro | Tidak | Ringan | 18130.318 | 21.5 |
| Yes | Petroleum extraction | 44 | 1978 | Mikro | Tidak | Ringan | 7784.797 | 14.3 |
| No | Petroleum extraction | 23 | 1979 | Mikro | Ya | Ringan | 4900.000 | 13.5 |
| No | Petroleum extraction | 3 | 2011 | Besar | Tidak | Kuat | 4841.117 | 9.8 |
| Yes | Petroleum extraction | 36 | 1986 | Mikro | Ya | Ringan | 5840.898 | 13.4 |
| Yes | Petroleum extraction | 42 | 1980 | Menengah | Ya | Ringan | 5885.018 | 13.4 |
| Yes | Petroleum extraction | 46 | 1976 | Besar | Ya | Ringan | 4303.062 | 7.2 |
5 Split Data Training & Testing
set.seed(42)
stratified_split <- function(y, prop = 0.8) {
idx_by_class <- split(seq_along(y), y)
unlist(lapply(idx_by_class,
function(idx) sample(idx, floor(length(idx) * prop))),
use.names = FALSE)
}
train_id <- stratified_split(contamination$kontaminasi_biner, 0.8)
train_data <- contamination[train_id, ]
test_data <- contamination[-train_id, ]
bind_rows(
train_data %>% count(status_kontaminasi) %>% mutate(Data = "Training"),
test_data %>% count(status_kontaminasi) %>% mutate(Data = "Testing")
) %>%
group_by(Data) %>%
mutate(Proporsi = scales::percent(n / sum(n), accuracy = 0.1)) %>%
ungroup() %>%
rename("Status Kontaminasi" = status_kontaminasi, Jumlah = n) %>%
knitr::kable(caption = "Distribusi Kelas Training dan Testing") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Status Kontaminasi | Jumlah | Data | Proporsi |
|---|---|---|---|
| No | 640 | Training | 39.9% |
| Yes | 963 | Training | 60.1% |
| No | 161 | Testing | 40.0% |
| Yes | 241 | Testing | 60.0% |
6 Pembangunan Model Secara Iteratif
Pengembangan model dilakukan secara iteratif melalui tiga tahap untuk memperbaiki pelanggaran asumsi dan menghasilkan model yang parsimonious.
6.1 Model v1 — Full Model Awal (14 Variabel)
6.1.1 Estimasi Model v1
contamination_fit1 <- glm(
kontaminasi_biner ~
durasi_operasi + tahun_mulai + jumlah_pelanggaran +
curah_hujan + suhu + kecepatan_angin +
tekstur_tanah + mobilitas_polutan + volatilitas_polutan +
permukaan_kedap + polutan_organik_persisten +
skala_industri + erosi_tanah + jenis_industri,
data = train_data,
family = binomial(link = "logit")
)
data.frame(
Keterangan = c("Jumlah observasi training", "Jumlah variabel prediktor",
"Null deviance", "Residual deviance", "Derajat bebas residual", "AIC"),
Nilai = c(nobs(contamination_fit1), length(contamination_fit1$coefficients) - 1,
round(contamination_fit1$null.deviance, 3),
round(contamination_fit1$deviance, 3),
contamination_fit1$df.residual, round(AIC(contamination_fit1), 3))
) %>%
knitr::kable(caption = "Ringkasan Kecocokan Model v1 (Full — 14 variabel)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 1603.000 |
| Jumlah variabel prediktor | 38.000 |
| Null deviance | 2156.699 |
| Residual deviance | 1411.060 |
| Derajat bebas residual | 1564.000 |
| AIC | 1489.060 |
6.1.2 Uji Multikolinearitas (VIF) — Model v1
vif_v1 <- car::vif(contamination_fit1)
vif_mat_v1 <- as.data.frame(vif_v1)
var_kategorik <- c("permukaan_kedap", "polutan_organik_persisten",
"skala_industri", "erosi_tanah", "jenis_industri")
tabel_vif_num_v1 <- vif_mat_v1[!rownames(vif_mat_v1) %in% var_kategorik, , drop = FALSE] %>%
transmute(VIF = round(GVIF, 6),
Keterangan = ifelse(GVIF > 10, "Bermasalah", "Aman"))
tabel_vif_kat_v1 <- vif_mat_v1[rownames(vif_mat_v1) %in% var_kategorik, , drop = FALSE] %>%
transmute(GVIF = round(GVIF, 6),
Df = Df,
`GVIF^(1/(2Df))` = round(`GVIF^(1/(2*Df))`, 6),
Keterangan = ifelse(`GVIF^(1/(2*Df))` > sqrt(10), "Bermasalah", "Aman"))
knitr::kable(tabel_vif_num_v1,
caption = "VIF Variabel Numerik — Model v1 (cutoff > 10)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| VIF | Keterangan | |
|---|---|---|
| durasi_operasi | 1.992835 | Aman |
| tahun_mulai | 2.036131 | Aman |
| jumlah_pelanggaran | 1.106675 | Aman |
| curah_hujan | 3.027219 | Aman |
| suhu | 3.387764 | Aman |
| kecepatan_angin | 1.464081 | Aman |
| tekstur_tanah | 1.574873 | Aman |
| mobilitas_polutan | 2.716300 | Aman |
| volatilitas_polutan | 85.554369 | Bermasalah |
knitr::kable(tabel_vif_kat_v1,
caption = "GVIF Variabel Kategorik — Model v1 (cutoff GVIF^(1/(2Df)) > √10 ≈ 3,16)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| GVIF | Df | GVIF^(1/(2Df)) | Keterangan | |
|---|---|---|---|---|
| permukaan_kedap | 1.196102 | 1 | 1.093664 | Aman |
| polutan_organik_persisten | 15.095518 | 1 | 3.885295 | Bermasalah |
| skala_industri | 1.663734 | 4 | 1.065701 | Aman |
| erosi_tanah | 1.352301 | 5 | 1.030641 | Aman |
| jenis_industri | 4636.502745 | 18 | 1.264266 | Aman |
Temuan v1 — VIF: Dua variabel bermasalah:
volatilitas_polutan: VIF = 85,55 > 10 → multikolinearitas beratpolutan_organik_persisten: GVIF^(1/(2Df)) = 3,89 > √10 ≈ 3,16 → melebihi batas GVIF kategorikKeputusan: Hapus kedua variabel tersebut dari model → model v2.
6.1.3 Uji Linearitas Logit (Box-Tidwell) — Model v1
train_bt1 <- train_data %>%
mutate(
dur_log = durasi_operasi * log(durasi_operasi + 1),
hujan_log = curah_hujan * log(curah_hujan + 1),
suhu_log = suhu * log(abs(suhu) + 1),
angin_log = kecepatan_angin * log(kecepatan_angin + 1),
tanah_log = tekstur_tanah * log(tekstur_tanah + 1),
mob_log = mobilitas_polutan * log(mobilitas_polutan + 1),
vol_log = volatilitas_polutan * log(volatilitas_polutan + 1)
)
bt_fit1 <- glm(
kontaminasi_biner ~
durasi_operasi + dur_log + curah_hujan + hujan_log +
suhu + suhu_log + kecepatan_angin + angin_log +
tekstur_tanah + tanah_log + mobilitas_polutan + mob_log +
volatilitas_polutan + vol_log,
data = train_bt1, family = binomial
)
broom::tidy(bt_fit1) %>%
filter(grepl("_log$", term)) %>%
transmute(
`Variabel (interaksi log)` = term,
Koefisien = round(estimate, 4),
`p-value` = signif(p.value, 3),
Kesimpulan = ifelse(p.value < 0.05, "Asumsi DILANGGAR", "Asumsi terpenuhi")
) %>%
knitr::kable(caption = "Uji Box-Tidwell v1: Linearitas Logit Variabel Kontinu") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Variabel (interaksi log) | Koefisien | p-value | Kesimpulan |
|---|---|---|---|
| dur_log | 0.0144 | 0.63200 | Asumsi terpenuhi |
| hujan_log | 0.0000 | 0.60900 | Asumsi terpenuhi |
| suhu_log | 0.0506 | 0.19700 | Asumsi terpenuhi |
| angin_log | 0.0228 | 0.21000 | Asumsi terpenuhi |
| tanah_log | 0.0101 | 0.76200 | Asumsi terpenuhi |
| mob_log | -199.1905 | 0.00433 | Asumsi DILANGGAR |
| vol_log | -0.6576 | 0.06170 | Asumsi terpenuhi |
Temuan v1 — Box-Tidwell:
mobilitas_polutanmelanggar asumsi linearitas logit (p = 0,004 < 0,05). Keputusan: Transformasilog(mobilitas_polutan + 1)→ dipakai sebagailog_mobilitas_polutanpada v2 & v3.
6.1.4 Uji Signifikansi Per Variabel (LRT) — Model v1
## Analysis of Deviance Table (Type II tests)
##
## Response: kontaminasi_biner
## LR Chisq Df Pr(>Chisq)
## durasi_operasi 88.443 1 < 2.2e-16 ***
## tahun_mulai 8.100 1 0.004428 **
## jumlah_pelanggaran 149.124 1 < 2.2e-16 ***
## curah_hujan 23.646 1 1.158e-06 ***
## suhu 0.168 1 0.682219
## kecepatan_angin 3.576 1 0.058637 .
## tekstur_tanah 8.383 1 0.003788 **
## mobilitas_polutan 5.270 1 0.021696 *
## volatilitas_polutan 0.016 1 0.899957
## permukaan_kedap 0.325 1 0.568672
## polutan_organik_persisten 8.165 1 0.004270 **
## skala_industri 1.459 4 0.833867
## erosi_tanah 3.628 5 0.604182
## jenis_industri 97.302 18 6.894e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.1.5 Prediksi & Evaluasi — Model v1
p_train_v1 <- predict(contamination_fit1, newdata = train_data, type = "response")
p_test_v1 <- predict(contamination_fit1, newdata = test_data, type = "response")
roc_test_v1 <- roc_points(test_data$kontaminasi_biner, p_test_v1) %>%
mutate(data = "Testing v1")
auc_test_v1 <- auc_value(roc_test_v1)
met_v1_def <- classification_metrics(test_data$kontaminasi_biner, p_test_v1, 0.5)
confusion_matrix_tabel(test_data$kontaminasi_biner, p_test_v1, 0.5) %>%
knitr::kable(caption = "Confusion Matrix Testing (Threshold = 0,50) — v1") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Prediksi: Terkontaminasi | Prediksi: Bersih | Sum | |
|---|---|---|---|
| Aktual: Terkontaminasi | 192 | 49 | 241 |
| Aktual: Bersih | 52 | 109 | 161 |
| Sum | 244 | 158 | 402 |
format_metrics_indonesia(met_v1_def) %>%
mutate(across(where(is.numeric), round, 3)) %>%
knitr::kable(caption = "Metrik Evaluasi Testing (Threshold = 0,50) — v1") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Threshold | Akurasi | Error Rate | Sensitivity | Specificity | Presisi | F1-Score | Balanced Accuracy |
|---|---|---|---|---|---|---|---|
| 0.5 | 0.749 | 0.251 | 0.797 | 0.677 | 0.787 | 0.792 | 0.737 |
# Kurva ROC v1
roc_train_v1 <- roc_points(train_data$kontaminasi_biner, p_train_v1) %>%
mutate(data = "Training v1")
optimal_v1 <- roc_train_v1 %>%
filter(is.finite(threshold)) %>% arrange(desc(youden), desc(sensitivity)) %>% slice(1)
threshold_opt_v1 <- optimal_v1$threshold[1]
ggplot(bind_rows(roc_train_v1, roc_test_v1), aes(x = fpr, y = sensitivity, color = data)) +
geom_path(linewidth = 1.1) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#6c757d") +
geom_point(data = optimal_v1, aes(x = fpr, y = sensitivity),
inherit.aes = FALSE, color = "#854d0e", fill = "#eab308",
shape = 21, size = 4, stroke = 1.2) +
coord_equal() +
scale_color_manual(values = c("Training v1" = "#2d6a2d", "Testing v1" = "#eab308")) +
labs(title = "Kurva ROC — Model v1 (Full — 14 variabel)",
subtitle = paste0("AUC Training = ", round(auc_value(roc_train_v1), 3),
" | AUC Testing = ", round(auc_test_v1, 3),
" | Threshold Optimal = ", round(threshold_opt_v1, 3)),
x = "False Positive Rate (1 - Specificity)",
y = "Sensitivity (True Positive Rate)", color = "Data") +
theme_minimal(base_size = 12) + theme(legend.position = "bottom")
data.frame(Data = c("Training v1","Testing v1"),
AUC = round(c(auc_value(roc_train_v1), auc_test_v1), 3)) %>%
knitr::kable(caption = "Nilai AUC — Model v1") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Data | AUC |
|---|---|
| Training v1 | 0.863 |
| Testing v1 | 0.841 |
AUC v1: AUC sebesar 0,863 (training) dan 0,841 (testing) menunjukkan kemampuan diskriminasi yang baik. Namun, model v1 masih belum memenuhi seluruh asumsi yang diperlukan karena terdapat multikolinearitas dan pelanggaran linearitas logit.
6.2 Model v2 — Full Model Terkoreksi (12 Variabel)
Perubahan dari v1:
- [DIHAPUS]
volatilitas_polutan→ VIF = 85,55 > 10 (multikolinearitas berat)- [DIHAPUS]
polutan_organik_persisten→ GVIF^(1/(2Df)) = 3,89 > √10 ≈ 3,16 (melebihi batas GVIF kategorik)- [DITRANSFORMASI]
mobilitas_polutan→log(x+1)=log_mobilitas_polutan
6.2.1 Estimasi Model v2
contamination_fit2 <- glm(
kontaminasi_biner ~
durasi_operasi + tahun_mulai + jumlah_pelanggaran +
curah_hujan + suhu + kecepatan_angin +
tekstur_tanah +
log_mobilitas_polutan + # [v2] sudah log-transform
# volatilitas_polutan # [v2] DIHAPUS: VIF = 85,55 > 10
# polutan_organik_persisten # [v2] DIHAPUS: GVIF^(1/(2Df)) = 3,89 > sqrt(10)
permukaan_kedap +
skala_industri + erosi_tanah + jenis_industri,
data = train_data,
family = binomial(link = "logit")
)
data.frame(
Keterangan = c("Jumlah observasi training", "Jumlah variabel prediktor",
"Null deviance", "Residual deviance", "Derajat bebas residual", "AIC"),
Nilai = c(nobs(contamination_fit2), length(contamination_fit2$coefficients) - 1,
round(contamination_fit2$null.deviance, 3),
round(contamination_fit2$deviance, 3),
contamination_fit2$df.residual, round(AIC(contamination_fit2), 3))
) %>%
knitr::kable(caption = "Ringkasan Kecocokan Model v2 (12 variabel)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 1603.000 |
| Jumlah variabel prediktor | 36.000 |
| Null deviance | 2156.699 |
| Residual deviance | 1419.257 |
| Derajat bebas residual | 1566.000 |
| AIC | 1493.257 |
6.2.2 Uji Multikolinearitas (VIF) — Model v2
vif_v2 <- car::vif(contamination_fit2)
vif_mat_v2 <- as.data.frame(vif_v2)
var_kategorik <- c("permukaan_kedap", "skala_industri", "erosi_tanah", "jenis_industri")
tabel_vif_num_v2 <- vif_mat_v2[!rownames(vif_mat_v2) %in% var_kategorik, , drop = FALSE] %>%
transmute(VIF = round(GVIF, 6),
Keterangan = ifelse(GVIF > 10, "Bermasalah", "Aman"))
tabel_vif_kat_v2 <- vif_mat_v2[rownames(vif_mat_v2) %in% var_kategorik, , drop = FALSE] %>%
transmute(GVIF = round(GVIF, 6),
Df = Df,
`GVIF^(1/(2Df))` = round(`GVIF^(1/(2*Df))`, 6),
Keterangan = ifelse(`GVIF^(1/(2*Df))` > sqrt(10), "Bermasalah", "Aman"))
knitr::kable(tabel_vif_num_v2,
caption = "VIF Variabel Numerik — Model v2 (cutoff > 10)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| VIF | Keterangan | |
|---|---|---|
| durasi_operasi | 1.994930 | Aman |
| tahun_mulai | 2.040736 | Aman |
| jumlah_pelanggaran | 1.102698 | Aman |
| curah_hujan | 3.015168 | Aman |
| suhu | 3.383602 | Aman |
| kecepatan_angin | 1.460018 | Aman |
| tekstur_tanah | 1.574520 | Aman |
| log_mobilitas_polutan | 2.424257 | Aman |
knitr::kable(tabel_vif_kat_v2,
caption = "GVIF Variabel Kategorik — Model v2 (cutoff GVIF^(1/(2Df)) > √10 ≈ 3,16)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| GVIF | Df | GVIF^(1/(2Df)) | Keterangan | |
|---|---|---|---|---|
| permukaan_kedap | 1.193966 | 1 | 1.092688 | Aman |
| skala_industri | 1.642345 | 4 | 1.063979 | Aman |
| erosi_tanah | 1.336795 | 5 | 1.029453 | Aman |
| jenis_industri | 6.536135 | 18 | 1.053532 | Aman |
Setelah penghapusan
volatilitas_polutandanpolutan_organik_persisten, seluruh variabel memiliki VIF < 10 dan GVIF^(1/(2Df)) < √10 ≈ 3,16 sehingga masalah multikolinearitas telah teratasi sepenuhnya.
6.2.3 Uji Linearitas Logit (Box-Tidwell) — Model v2
train_bt2 <- train_data %>%
mutate(
dur_log = durasi_operasi * log(durasi_operasi + 1),
hujan_log = curah_hujan * log(curah_hujan + 1),
suhu_log = suhu * log(abs(suhu) + 1),
angin_log = kecepatan_angin * log(kecepatan_angin + 1),
tanah_log = tekstur_tanah * log(tekstur_tanah + 1)
)
bt_fit2 <- glm(
kontaminasi_biner ~
durasi_operasi + dur_log + curah_hujan + hujan_log +
suhu + suhu_log + kecepatan_angin + angin_log +
tekstur_tanah + tanah_log,
data = train_bt2, family = binomial
)
broom::tidy(bt_fit2) %>%
filter(grepl("_log$", term)) %>%
transmute(
`Variabel (interaksi log)` = term,
Koefisien = round(estimate, 4),
`p-value` = signif(p.value, 3),
Kesimpulan = ifelse(p.value < 0.05, "Asumsi DILANGGAR", "Asumsi terpenuhi")
) %>%
knitr::kable(caption = "Uji Box-Tidwell v2: Linearitas Logit (Variabel Tersisa)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Variabel (interaksi log) | Koefisien | p-value | Kesimpulan |
|---|---|---|---|
| dur_log | 0.0168 | 0.580 | Asumsi terpenuhi |
| hujan_log | 0.0000 | 0.328 | Asumsi terpenuhi |
| suhu_log | 0.0402 | 0.298 | Asumsi terpenuhi |
| angin_log | 0.0254 | 0.154 | Asumsi terpenuhi |
| tanah_log | 0.0043 | 0.896 | Asumsi terpenuhi |
Semua p-value > 0,05 → asumsi linearitas logit terpenuhi untuk semua variabel numerik di v2.
6.2.4 Uji Signifikansi Per Variabel (LRT) — Model v2
## Analysis of Deviance Table (Type II tests)
##
## Response: kontaminasi_biner
## LR Chisq Df Pr(>Chisq)
## durasi_operasi 86.410 1 < 2.2e-16 ***
## tahun_mulai 8.286 1 0.003996 **
## jumlah_pelanggaran 148.947 1 < 2.2e-16 ***
## curah_hujan 25.144 1 5.320e-07 ***
## suhu 0.232 1 0.630183
## kecepatan_angin 3.546 1 0.059688 .
## tekstur_tanah 9.163 1 0.002470 **
## log_mobilitas_polutan 8.364 1 0.003828 **
## permukaan_kedap 0.267 1 0.605459
## skala_industri 1.584 4 0.811610
## erosi_tanah 3.891 5 0.565142
## jenis_industri 92.770 18 4.579e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.2.5 Koefisien & Odds Ratio — Model v2
coef_v2 <- broom::tidy(contamination_fit2) %>%
filter(term != "(Intercept)") %>%
mutate(
odds_ratio = exp(estimate),
ci_low = exp(estimate - 1.96 * std.error),
ci_high = exp(estimate + 1.96 * std.error),
Signifikansi = case_when(
p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*", p.value < 0.1 ~ ".",
TRUE ~ "tidak signifikan"
)
) %>%
arrange(p.value) %>%
transmute(`Variabel/Level` = term,
`Odds Ratio` = round(odds_ratio, 3),
`CI 95%` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
`p-value` = signif(p.value, 3),
Signifikansi = Signifikansi)
knitr::kable(coef_v2, caption = "Koefisien & Odds Ratio Lengkap — Model v2 (diurutkan p-value)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Variabel/Level | Odds Ratio | CI 95% | p-value | Signifikansi |
|---|---|---|---|---|
| jumlah_pelanggaran | 1.549 | 1.421 - 1.688 | 0.00e+00 | *** |
| durasi_operasi | 1.140 | 1.107 - 1.173 | 0.00e+00 | *** |
| curah_hujan | 1.000 | 1 - 1 | 7.00e-07 | *** |
| jenis_industriHazardous waste treatment | 0.077 | 0.027 - 0.223 | 2.10e-06 | *** |
| jenis_industriCrude oil processing and petroleum product manufacturing | 3.517 | 1.758 - 7.036 | 3.77e-04 | *** |
| jenis_industriOrganic chemical raw material manufacturing | 3.045 | 1.634 - 5.676 | 4.54e-04 | *** |
| jenis_industriCotton printing and dyeing processing | 0.365 | 0.207 - 0.644 | 4.97e-04 | *** |
| tekstur_tanah | 0.968 | 0.948 - 0.989 | 2.64e-03 | ** |
| log_mobilitas_polutan | 0.420 | 0.233 - 0.759 | 4.01e-03 | ** |
| tahun_mulai | 0.964 | 0.941 - 0.989 | 4.26e-03 | ** |
| jenis_industriLead-zinc smelting | 2.634 | 1.198 - 5.791 | 1.59e-02 |
|
| jenis_industriPesticide manufacturing | 1.977 | 1.073 - 3.644 | 2.88e-02 |
|
| jenis_industriCoking | 1.923 | 0.977 - 3.785 | 5.84e-02 | . |
| kecepatan_angin | 1.016 | 0.999 - 1.032 | 6.19e-02 | . |
| jenis_industriLeather tanning and processing | 1.803 | 0.961 - 3.382 | 6.62e-02 | . |
| jenis_industriPetroleum extraction | 0.393 | 0.14 - 1.102 | 7.59e-02 | . |
| jenis_industriSteel smelting | 2.307 | 0.848 - 6.271 | 1.01e-01 | tidak signifikan |
| erosi_tanahSangat Kuat | 3.948 | 0.676 - 23.06 | 1.27e-01 | tidak signifikan |
| jenis_industriChemical reagent and additive manufacturing | 0.700 | 0.375 - 1.308 | 2.63e-01 | tidak signifikan |
| erosi_tanahAgak Ringan | 1.322 | 0.808 - 2.164 | 2.67e-01 | tidak signifikan |
| jenis_industriIron ore mining and dressing | 0.674 | 0.317 - 1.435 | 3.07e-01 | tidak signifikan |
| skala_industriKecil | 1.219 | 0.829 - 1.792 | 3.15e-01 | tidak signifikan |
| jenis_industriIron smelting | 1.524 | 0.63 - 3.686 | 3.50e-01 | tidak signifikan |
| jenis_industriLead-zinc ore mining and dressing | 1.500 | 0.582 - 3.868 | 4.01e-01 | tidak signifikan |
| permukaan_kedapYa | 0.928 | 0.698 - 1.234 | 6.06e-01 | tidak signifikan |
| suhu | 0.988 | 0.94 - 1.039 | 6.30e-01 | tidak signifikan |
| jenis_industriMetal surface treatment and heat treatment processing | 1.117 | 0.668 - 1.866 | 6.73e-01 | tidak signifikan |
| jenis_industriGold ore mining and dressing | 1.217 | 0.442 - 3.349 | 7.03e-01 | tidak signifikan |
| jenis_industriCopper smelting | 1.077 | 0.425 - 2.73 | 8.76e-01 | tidak signifikan |
| erosi_tanahKuat | 1.088 | 0.374 - 3.162 | 8.77e-01 | tidak signifikan |
| jenis_industriChemical pharmaceutical raw material manufacturing | 0.945 | 0.373 - 2.396 | 9.05e-01 | tidak signifikan |
| skala_industriMenengah | 0.979 | 0.641 - 1.493 | 9.20e-01 | tidak signifikan |
| erosi_tanahSedang | 0.969 | 0.517 - 1.816 | 9.22e-01 | tidak signifikan |
| skala_industriSangat Besar | 0.982 | 0.612 - 1.575 | 9.40e-01 | tidak signifikan |
| skala_industriBesar | 0.991 | 0.657 - 1.495 | 9.67e-01 | tidak signifikan |
| erosi_tanahParah | 1279.112 | 0 - 3.42645693810851e+279 | 9.82e-01 | tidak signifikan |
Keempat variabel tersebut, yaitu suhu, permukaan kedap, skala industri, dan erosi tanah, dieliminasi pada model v3 karena tidak memenuhi kriteria yang ditetapkan dalam proses seleksi variabel. Selain itu, variabel erosi tanah pada kategori “Parah” menunjukkan indikasi quasi-complete separation, yang berpotensi menyebabkan ketidakstabilan dalam estimasi parameter model.
6.2.6 Prediksi & Evaluasi — Model v2
p_train_v2 <- predict(contamination_fit2, newdata = train_data, type = "response")
p_test_v2 <- predict(contamination_fit2, newdata = test_data, type = "response")
confusion_matrix_tabel(test_data$kontaminasi_biner, p_test_v2, 0.5) %>%
knitr::kable(caption = "Confusion Matrix Testing (Threshold = 0,50) — v2") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Prediksi: Terkontaminasi | Prediksi: Bersih | Sum | |
|---|---|---|---|
| Aktual: Terkontaminasi | 191 | 50 | 241 |
| Aktual: Bersih | 53 | 108 | 161 |
| Sum | 244 | 158 | 402 |
met_v2_def <- classification_metrics(test_data$kontaminasi_biner, p_test_v2, 0.5)
format_metrics_indonesia(met_v2_def) %>%
mutate(across(where(is.numeric), round, 3)) %>%
knitr::kable(caption = "Metrik Evaluasi Testing (Threshold = 0,50) — v2") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Threshold | Akurasi | Error Rate | Sensitivity | Specificity | Presisi | F1-Score | Balanced Accuracy |
|---|---|---|---|---|---|---|---|
| 0.5 | 0.744 | 0.256 | 0.793 | 0.671 | 0.783 | 0.788 | 0.732 |
roc_train_v2 <- roc_points(train_data$kontaminasi_biner, p_train_v2) %>% mutate(data = "Training v2")
roc_test_v2 <- roc_points(test_data$kontaminasi_biner, p_test_v2) %>% mutate(data = "Testing v2")
auc_train_v2 <- auc_value(roc_train_v2); auc_test_v2 <- auc_value(roc_test_v2)
data.frame(Data = c("Training v2","Testing v2"),
AUC = round(c(auc_train_v2, auc_test_v2), 3)) %>%
knitr::kable(caption = "Nilai AUC — Model v2") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Data | AUC |
|---|---|
| Training v2 | 0.860 |
| Testing v2 | 0.837 |
ggplot(bind_rows(roc_train_v2, roc_test_v2), aes(x = fpr, y = sensitivity, color = data)) +
geom_path(linewidth = 1.1) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#6c757d") +
coord_equal() +
scale_color_manual(values = c("Training v2" = "#2d6a2d", "Testing v2" = "#eab308")) +
labs(title = "Kurva ROC — Model v2 (Full Diperbaiki — 12 variabel)",
subtitle = paste0("AUC Training = ", round(auc_train_v2, 3),
" | AUC Testing = ", round(auc_test_v2, 3)),
x = "False Positive Rate (1 - Specificity)",
y = "Sensitivity (True Positive Rate)", color = "Data") +
theme_minimal(base_size = 12) + theme(legend.position = "bottom")
AUC v2: Model v2 tetap mampu membedakan lokasi terkontaminasi dan tidak terkontaminasi dengan baik (AUC training = 0,860; testing = 0,837) setelah perbaikan masalah multikolinearitas dan linearitas logit. Namun, hasil uji Hosmer–Lemeshow (p = 0,074 < 0,05 tidak terpenuhi secara ketat) masih mengindikasikan adanya misfit, sehingga dilakukan eliminasi variabel yang tidak signifikan pada tahap berikutnya.
6.2.7 Uji Hosmer-Lemeshow & Nagelkerke R² — Model v2
hl_v2 <- ResourceSelection::hoslem.test(train_data$kontaminasi_biner, fitted(contamination_fit2), g = 10)
print(hl_v2)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train_data$kontaminasi_biner, fitted(contamination_fit2)
## X-squared = 14.31, df = 8, p-value = 0.07403nagelkerke_v2 <- round(DescTools::PseudoR2(contamination_fit2, which = "Nagelkerke"), 4)
cat("Nagelkerke R² v2:", nagelkerke_v2)## Nagelkerke R² v2: 0.4986Hosmer-Lemeshow v2: menghasilkan χ²(8) = 14,310 (p = 0,074), yang mengindikasikan bahwa model v2 masih belum sepenuhnya fit. Oleh karena itu, dilakukan eliminasi variabel yang tidak signifikan pada model v3.
Nagelkerke R² v2: sebesar 0,4986 menunjukkan bahwa model mampu menjelaskan sekitar 49,9% variasi status kontaminasi lahan.
6.3 Model v3 — Reduced Model / Model Final (8 Variabel)
Perubahan dari v2:
- [DIHAPUS]
suhu(p = 0,630)- [DIHAPUS]
permukaan_kedap(p = 0,606)- [DIHAPUS]
skala_industri(semua level p > 0,31)- [DIHAPUS]
erosi_tanah(semua level p > 0,13 + quasi-complete separation)- [DIPERTAHANKAN]
kecepatan_angin(p = 0,060 borderline, relevan substantif)
6.3.1 Estimasi Model v3
contamination_fit3 <- glm(
kontaminasi_biner ~
durasi_operasi + tahun_mulai + jumlah_pelanggaran +
curah_hujan + kecepatan_angin + tekstur_tanah +
log_mobilitas_polutan +
# polutan_organik_persisten # [DIHAPUS di v2] GVIF^(1/(2Df)) = 3,89 > sqrt(10)
jenis_industri,
data = train_data,
family = binomial(link = "logit")
)
data.frame(
Keterangan = c("Jumlah observasi training", "Jumlah variabel prediktor",
"Null deviance", "Residual deviance", "Derajat bebas residual", "AIC"),
Nilai = c(nobs(contamination_fit3), length(contamination_fit3$coefficients) - 1,
round(contamination_fit3$null.deviance, 3),
round(contamination_fit3$deviance, 3),
contamination_fit3$df.residual, round(AIC(contamination_fit3), 3))
) %>%
knitr::kable(caption = "Ringkasan Kecocokan Model v3 (Reduced / Final)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Keterangan | Nilai |
|---|---|
| Jumlah observasi training | 1603.000 |
| Jumlah variabel prediktor | 25.000 |
| Null deviance | 2156.699 |
| Residual deviance | 1425.287 |
| Derajat bebas residual | 1577.000 |
| AIC | 1477.287 |
6.3.2 Perbandingan AIC: v2 vs v3 & Uji LRT
aic_v2 <- AIC(contamination_fit2); aic_v3 <- AIC(contamination_fit3)
selisih_aic <- round(abs(aic_v2 - aic_v3), 3)
data.frame(
Model = c("v2 — Full Model (12 variabel)", "v3 — Reduced Model (8 variabel)"),
`Jumlah Variabel` = c(12, 8),
`Residual Deviance` = c(round(contamination_fit2$deviance, 3),
round(contamination_fit3$deviance, 3)),
`Df Residual` = c(contamination_fit2$df.residual, contamination_fit3$df.residual),
AIC = c(round(aic_v2, 3), round(aic_v3, 3)),
check.names = FALSE
) %>%
knitr::kable(caption = "Perbandingan Full Model (v2) vs Reduced Model (v3)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Model | Jumlah Variabel | Residual Deviance | Df Residual | AIC |
|---|---|---|---|---|
| v2 — Full Model (12 variabel) | 12 | 1419.257 | 1566 | 1493.257 |
| v3 — Reduced Model (8 variabel) | 8 | 1425.287 | 1577 | 1477.287 |
## Selisih AIC (v2 - v3): 15.97 -> ΔAIC > 10 -> bukti kuat memilih v3## --- LIKELIHOOD RATIO TEST (LRT): v3 vs v2 ---## Analysis of Deviance Table
##
## Model 1: kontaminasi_biner ~ durasi_operasi + tahun_mulai + jumlah_pelanggaran +
## curah_hujan + kecepatan_angin + tekstur_tanah + log_mobilitas_polutan +
## jenis_industri
## Model 2: kontaminasi_biner ~ durasi_operasi + tahun_mulai + jumlah_pelanggaran +
## curah_hujan + suhu + kecepatan_angin + tekstur_tanah + log_mobilitas_polutan +
## permukaan_kedap + skala_industri + erosi_tanah + jenis_industri
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 1577 1425.3
## 2 1566 1419.3 11 6.03 0.8714##
## p-value = 0.8714 > 0.05## -> Variabel yang dihapus TIDAK signifikan secara bersama-sama.## -> Reduced model (v3) VALID digunakan sebagai MODEL FINAL.LRT v3 vs v2: Hasil Likelihood Ratio Test (p = 0,871) menunjukkan bahwa penghapusan variabel suhu, permukaan kedap, skala industri, dan erosi tanah tidak menyebabkan penurunan kecocokan model yang signifikan. Oleh karena itu, model v3 dipilih sebagai model final berdasarkan prinsip parsimoni.
6.3.3 Uji Multikolinearitas (VIF) — Model v3
vif_v3 <- car::vif(contamination_fit3)
vif_mat_v3 <- as.data.frame(vif_v3)
var_kategorik_v3 <- c("skala_industri", "erosi_tanah", "jenis_industri")
tabel_vif_num_v3 <- vif_mat_v3[!rownames(vif_mat_v3) %in% var_kategorik_v3, , drop = FALSE] %>%
transmute(VIF = round(GVIF, 6),
Keterangan = ifelse(GVIF > 10, "Bermasalah", "Aman"))
tabel_vif_kat_v3 <- vif_mat_v3[rownames(vif_mat_v3) %in% var_kategorik_v3, , drop = FALSE] %>%
transmute(GVIF = round(GVIF, 6),
Df = Df,
`GVIF^(1/(2Df))` = round(`GVIF^(1/(2*Df))`, 6),
Keterangan = ifelse(`GVIF^(1/(2*Df))` > sqrt(10), "Bermasalah", "Aman"))
knitr::kable(tabel_vif_num_v3,
caption = "VIF Variabel Numerik — Model v3 (cutoff > 10)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| VIF | Keterangan | |
|---|---|---|
| durasi_operasi | 1.836481 | Aman |
| tahun_mulai | 1.772312 | Aman |
| jumlah_pelanggaran | 1.087497 | Aman |
| curah_hujan | 1.602981 | Aman |
| kecepatan_angin | 1.226603 | Aman |
| tekstur_tanah | 1.494252 | Aman |
| log_mobilitas_polutan | 2.391970 | Aman |
knitr::kable(tabel_vif_kat_v3,
caption = "GVIF Variabel Kategorik — Model v3 (cutoff GVIF^(1/(2Df)) > √10 ≈ 3,16)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| GVIF | Df | GVIF^(1/(2Df)) | Keterangan | |
|---|---|---|---|---|
| jenis_industri | 3.797707 | 18 | 1.037762 | Aman |
6.3.4 Uji Signifikansi Per Variabel (LRT) — Model v3
## Analysis of Deviance Table (Type II tests)
##
## Response: kontaminasi_biner
## LR Chisq Df Pr(>Chisq)
## durasi_operasi 91.853 1 < 2.2e-16 ***
## tahun_mulai 10.200 1 0.001405 **
## jumlah_pelanggaran 152.671 1 < 2.2e-16 ***
## curah_hujan 42.013 1 9.069e-11 ***
## kecepatan_angin 6.844 1 0.008894 **
## tekstur_tanah 10.597 1 0.001132 **
## log_mobilitas_polutan 9.254 1 0.002350 **
## jenis_industri 96.170 18 1.108e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.3.5 Koefisien, Odds Ratio & p-value — Model v3 (Final)
coef_v3 <- broom::tidy(contamination_fit3) %>%
filter(term != "(Intercept)") %>%
mutate(
odds_ratio = exp(estimate),
ci_low = exp(estimate - 1.96 * std.error),
ci_high = exp(estimate + 1.96 * std.error),
Signifikansi = case_when(
p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*", p.value < 0.1 ~ ".",
TRUE ~ "-"
)
) %>%
arrange(p.value) %>%
transmute(`Variabel/Level` = term,
`Odds Ratio` = round(odds_ratio, 3),
`CI 95% (Low - High)` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
`p-value` = signif(p.value, 3),
Signifikansi = Signifikansi)
knitr::kable(coef_v3, caption = "Koefisien & Odds Ratio — Model v3 (Reduced / Final)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Variabel/Level | Odds Ratio | CI 95% (Low - High) | p-value | Signifikansi |
|---|---|---|---|---|
| jumlah_pelanggaran | 1.549 | 1.422 - 1.687 | 0.00e+00 | *** |
| durasi_operasi | 1.137 | 1.106 - 1.169 | 0.00e+00 | *** |
| curah_hujan | 1.000 | 1 - 1 | 0.00e+00 | *** |
| jenis_industriHazardous waste treatment | 0.077 | 0.027 - 0.219 | 1.60e-06 | *** |
| jenis_industriCotton printing and dyeing processing | 0.353 | 0.204 - 0.61 | 1.96e-04 | *** |
| jenis_industriCrude oil processing and petroleum product manufacturing | 3.396 | 1.707 - 6.758 | 4.96e-04 | *** |
| jenis_industriOrganic chemical raw material manufacturing | 2.899 | 1.569 - 5.359 | 6.83e-04 | *** |
| tekstur_tanah | 0.967 | 0.947 - 0.987 | 1.23e-03 | ** |
| tahun_mulai | 0.963 | 0.941 - 0.986 | 1.56e-03 | ** |
| log_mobilitas_polutan | 0.405 | 0.226 - 0.728 | 2.49e-03 | ** |
| jenis_industriLead-zinc smelting | 2.783 | 1.282 - 6.044 | 9.67e-03 | ** |
| kecepatan_angin | 1.020 | 1.005 - 1.035 | 9.84e-03 | ** |
| jenis_industriPesticide manufacturing | 1.896 | 1.036 - 3.471 | 3.81e-02 |
|
| jenis_industriCoking | 1.969 | 1.006 - 3.854 | 4.81e-02 |
|
| jenis_industriPetroleum extraction | 0.394 | 0.146 - 1.059 | 6.47e-02 | . |
| jenis_industriLeather tanning and processing | 1.729 | 0.935 - 3.199 | 8.10e-02 | . |
| jenis_industriSteel smelting | 2.252 | 0.834 - 6.076 | 1.09e-01 |
|
| jenis_industriChemical reagent and additive manufacturing | 0.692 | 0.375 - 1.279 | 2.40e-01 |
|
| jenis_industriIron smelting | 1.565 | 0.652 - 3.754 | 3.16e-01 |
|
| jenis_industriIron ore mining and dressing | 0.698 | 0.334 - 1.457 | 3.38e-01 |
|
| jenis_industriLead-zinc ore mining and dressing | 1.544 | 0.602 - 3.957 | 3.66e-01 |
|
| jenis_industriGold ore mining and dressing | 1.257 | 0.467 - 3.384 | 6.51e-01 |
|
| jenis_industriCopper smelting | 1.110 | 0.442 - 2.791 | 8.24e-01 |
|
| jenis_industriChemical pharmaceutical raw material manufacturing | 0.906 | 0.364 - 2.252 | 8.32e-01 |
|
| jenis_industriMetal surface treatment and heat treatment processing | 1.038 | 0.641 - 1.68 | 8.80e-01 |
|
6.3.6 Prediksi & Confusion Matrix — Model v3
p_train_v3 <- predict(contamination_fit3, newdata = train_data, type = "response")
p_test_v3 <- predict(contamination_fit3, newdata = test_data, type = "response")
# Preview 8 data pertama
head(data.frame("Status Aktual" = as.character(test_data$status_kontaminasi),
"Peluang Terkontaminasi" = round(p_test_v3, 4),
check.names = FALSE), 8) %>%
knitr::kable(caption = "Contoh Peluang Prediksi v3 (8 Data Pertama)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Status Aktual | Peluang Terkontaminasi |
|---|---|
| No | 0.2109 |
| No | 0.2149 |
| Yes | 0.4692 |
| No | 0.2037 |
| No | 0.2727 |
| No | 0.1762 |
| Yes | 0.9728 |
| Yes | 0.9803 |
# Confusion Matrix threshold 0.50
confusion_matrix_tabel(test_data$kontaminasi_biner, p_test_v3, 0.5) %>%
knitr::kable(caption = "Confusion Matrix Testing (Threshold = 0,50) — v3") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Prediksi: Terkontaminasi | Prediksi: Bersih | Sum | |
|---|---|---|---|
| Aktual: Terkontaminasi | 190 | 51 | 241 |
| Aktual: Bersih | 48 | 113 | 161 |
| Sum | 238 | 164 | 402 |
met_v3_def <- classification_metrics(test_data$kontaminasi_biner, p_test_v3, 0.5)
format_metrics_indonesia(met_v3_def) %>%
mutate(across(where(is.numeric), \(x) round(x, 3))) %>%
knitr::kable(caption = "Metrik Evaluasi Testing (Threshold = 0,50) — v3") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Threshold | Akurasi | Error Rate | Sensitivity | Specificity | Presisi | F1-Score | Balanced Accuracy |
|---|---|---|---|---|---|---|---|
| 0.5 | 0.754 | 0.246 | 0.788 | 0.702 | 0.798 | 0.793 | 0.745 |
7 Evaluasi Kinerja Model v3
7.1 ROC & AUC — Model v3 (Final)
roc_train_v3 <- roc_points(train_data$kontaminasi_biner, p_train_v3) %>% mutate(data = "Training v3")
roc_test_v3 <- roc_points(test_data$kontaminasi_biner, p_test_v3) %>% mutate(data = "Testing v3")
auc_train_v3 <- auc_value(roc_train_v3); auc_test_v3 <- auc_value(roc_test_v3)
optimal_v3 <- roc_train_v3 %>%
filter(is.finite(threshold)) %>% arrange(desc(youden), desc(sensitivity)) %>% slice(1)
threshold_opt_v3 <- optimal_v3$threshold[1]
data.frame(Data = c("Training v3","Testing v3"),
AUC = round(c(auc_train_v3, auc_test_v3), 3)) %>%
knitr::kable(caption = "Nilai AUC — Model v3 (Final)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Data | AUC |
|---|---|
| Training v3 | 0.859 |
| Testing v3 | 0.847 |
## Threshold optimal v3 (Youden): 0.6589# Kurva ROC v3
ggplot(bind_rows(roc_train_v3, roc_test_v3), aes(x = fpr, y = sensitivity, color = data)) +
geom_path(linewidth = 1.2) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#6c757d") +
geom_point(data = optimal_v3, aes(x = fpr, y = sensitivity),
inherit.aes = FALSE, color = "#854d0e", fill = "#eab308",
shape = 21, size = 4.5, stroke = 1.3) +
coord_equal() +
scale_color_manual(values = c("Training v3" = "#2d6a2d", "Testing v3" = "#eab308")) +
labs(title = "Kurva ROC Model v3 (Reduced / Final)",
subtitle = paste0("AUC Training = ", round(auc_train_v3, 3),
" | AUC Testing = ", round(auc_test_v3, 3),
" | Threshold Optimal = ", round(threshold_opt_v3, 3)),
x = "False Positive Rate (1 - Specificity)",
y = "Sensitivity (True Positive Rate)", color = "Data") +
theme_minimal(base_size = 12) + theme(legend.position = "bottom")
AUC v3: Nilai AUC pada data training (0,859) dan testing (0,847) menunjukkan performa model yang konsisten, dengan selisih hanya 0,012 sehingga tidak terdapat indikasi overfitting yang berarti. Nilai AUC pengujian sebesar 0,847 termasuk kategori baik dan menunjukkan bahwa model mampu membedakan lokasi terkontaminasi dan tidak terkontaminasi dengan baik.
7.2 Threshold Default vs Threshold Optimal — Model v3
Threshold optimal ditentukan dengan memaksimalkan Indeks Youden (Persamaan 8):
\[J = \text{Sensitivitas} + \text{Spesifisitas} - 1\]
Metrik evaluasi yang digunakan (Persamaan 9–14):
| Metrik | Rumus |
|---|---|
| Akurasi | \((TP+TN)/(TP+TN+FP+FN)\) |
| Sensitivitas | \(TP/(TP+FN)\) |
| Spesifisitas | \(TN/(TN+FP)\) |
| Presisi | \(TP/(TP+FP)\) |
| F1-Score | \(2 \times (\text{Presisi} \times \text{Sensitivitas})/(\text{Presisi}+\text{Sensitivitas})\) |
| Balanced Accuracy | \((\text{Sensitivitas}+\text{Spesifisitas})/2\) |
met_v3_opt <- classification_metrics(test_data$kontaminasi_biner, p_test_v3, threshold_opt_v3)
confusion_matrix_tabel(test_data$kontaminasi_biner, p_test_v3, threshold_opt_v3) %>%
knitr::kable(caption = paste0("Confusion Matrix Testing (Threshold Optimal = ",
round(threshold_opt_v3, 3), ") — v3")) %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)| Prediksi: Terkontaminasi | Prediksi: Bersih | Sum | |
|---|---|---|---|
| Aktual: Terkontaminasi | 159 | 82 | 241 |
| Aktual: Bersih | 28 | 133 | 161 |
| Sum | 187 | 215 | 402 |
bind_rows(
met_v3_def %>% mutate(aturan = "Threshold 0.50"),
met_v3_opt %>% mutate(aturan = paste0("Threshold Optimal (", round(threshold_opt_v3, 3), ")"))
) %>%
format_metrics_indonesia() %>%
bind_cols(`Aturan Klasifikasi` = c("Threshold 0.50",
paste0("Threshold Optimal (", round(threshold_opt_v3, 3), ")")), .) %>%
select(`Aturan Klasifikasi`, everything()) %>%
mutate(across(where(is.numeric), round, 3)) %>%
knitr::kable(caption = "Perbandingan Metrik: Default vs Optimal — v3 (Final)") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Aturan Klasifikasi | Threshold | Akurasi | Error Rate | Sensitivity | Specificity | Presisi | F1-Score | Balanced Accuracy |
|---|---|---|---|---|---|---|---|---|
| Threshold 0.50 | 0.500 | 0.754 | 0.246 | 0.788 | 0.702 | 0.798 | 0.793 | 0.745 |
| Threshold Optimal (0.659) | 0.659 | 0.726 | 0.274 | 0.660 | 0.826 | 0.850 | 0.743 | 0.743 |
Perbandingan threshold: Threshold 0,50 menghasilkan sensitivitas lebih tinggi (79,7% vs 66,8%), sehingga lebih direkomendasikan untuk pengawasan lahan berisiko tinggi di mana meminimalkan lokasi kontaminasi yang terlewat menjadi prioritas utama. Dalam konteks identifikasi lahan terkontaminasi, meminimalkan lokasi terkontaminasi yang tidak terdeteksi menjadi prioritas utama, sehingga threshold 0,5 dinilai lebih sesuai untuk digunakan.
7.3 Distribusi Peluang Prediksi — Model v3
test_data %>%
mutate(peluang_kontaminasi = p_test_v3) %>%
ggplot(aes(x = peluang_kontaminasi, fill = status_kontaminasi)) +
geom_density(alpha = 0.55, color = "white", linewidth = 0.7) +
geom_vline(xintercept = threshold_opt_v3,
color = "#2d6a2d", linewidth = 1.2, linetype = "dashed") +
annotate("text", x = threshold_opt_v3 + 0.02, y = Inf,
label = paste0("Threshold = ", round(threshold_opt_v3, 3)),
vjust = 1.6, hjust = 0, color = "#2d6a2d", size = 3.5) +
scale_fill_manual(values = c("No" = "#2d6a2d", "Yes" = "#eab308")) +
labs(title = "Distribusi Peluang Prediksi Kontaminasi Lahan (Data Testing) — v3 Final",
subtitle = "Garis putus-putus = threshold optimal ROC",
x = "Peluang prediksi terkontaminasi", y = "Kepadatan",
fill = "Status Aktual") +
theme_minimal(base_size = 12)
Distribusi peluang prediksi: Distribusi peluang prediksi menunjukkan pemisahan yang cukup baik antara lokasi terkontaminasi dan tidak terkontaminasi. Lokasi terkontaminasi cenderung memiliki peluang prediksi yang lebih tinggi, sedangkan lokasi tidak terkontaminasi lebih banyak terkonsentrasi pada peluang yang rendah. Namun, masih terdapat area tumpang tindih di antara kedua distribusi yang mengindikasikan adanya observasi dengan karakteristik serupa, sehingga berpotensi menimbulkan kesalahan klasifikasi.
7.4 Kecocokan Model Final (Goodness-of-Fit)
hl_v3 <- ResourceSelection::hoslem.test(train_data$kontaminasi_biner, fitted(contamination_fit3), g = 10)
cat("--- UJI HOSMER-LEMESHOW — MODEL v3 ---\n")## --- UJI HOSMER-LEMESHOW — MODEL v3 ---##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train_data$kontaminasi_biner, fitted(contamination_fit3)
## X-squared = 12.556, df = 8, p-value = 0.128##
## --- OMNIBUS TEST (ANOVA LRT vs NULL) — MODEL v3 ---null_fit <- glm(kontaminasi_biner ~ 1, data = train_data, family = binomial)
omnibus_v3 <- anova(null_fit, contamination_fit3, test = "Chisq")
print(omnibus_v3)## Analysis of Deviance Table
##
## Model 1: kontaminasi_biner ~ 1
## Model 2: kontaminasi_biner ~ durasi_operasi + tahun_mulai + jumlah_pelanggaran +
## curah_hujan + kecepatan_angin + tekstur_tanah + log_mobilitas_polutan +
## jenis_industri
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 1602 2156.7
## 2 1577 1425.3 25 731.41 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nagelkerke_v3 <- round(DescTools::PseudoR2(contamination_fit3, which = "Nagelkerke"), 4)
cat("\n--- NAGELKERKE PSEUDO R² — MODEL v3 ---\n")##
## --- NAGELKERKE PSEUDO R² — MODEL v3 ---## Nagelkerke R² v3: 0.4954Hosmer-Lemeshow v3: nilai χ²(8) = 12,556 dengan p-value = 0,128, sehingga tidak terdapat perbedaan yang signifikan antara probabilitas prediksi model dan data observasi.
Omnibus LRT: χ² = 731,41, p < 0,001 menunjukkan bahwa model secara keseluruhan signifikan dan memberikan peningkatan yang nyata dibandingkan model tanpa prediktor.
Nagelkerke R² v3: sebesar 0,4954 mengindikasikan bahwa sekitar 49,5% variasi status kontaminasi lahan dapat dijelaskan oleh variabel-variabel dalam model.
Hasil Evaluasi Goodness-of-Fit:
- Hosmer-Lemeshow: p = 0,128 > 0,05 → model fit dengan data
- Omnibus LRT: χ² = 731,41, p < 0,001 → model secara keseluruhan signifikan
- Nagelkerke R²: 0,495 → model menjelaskan ~49,5% variasi kontaminasi
8 Perbandingan Akhir: v1 vs v2 vs v3
8.1 Tabel Perbandingan Lengkap
roc_test_v2_plot <- roc_points(test_data$kontaminasi_biner, p_test_v2) %>% mutate(data = "Testing v2")
met_v1_full <- classification_metrics(test_data$kontaminasi_biner, p_test_v1, 0.5)
met_v2_full <- classification_metrics(test_data$kontaminasi_biner, p_test_v2, 0.5)
met_v3_full <- classification_metrics(test_data$kontaminasi_biner, p_test_v3, 0.5)
data.frame(
Metrik = c("Jumlah Variabel","AIC","Nagelkerke R²","AUC Training","AUC Testing",
"Akurasi (threshold 0.5)","Sensitivity (threshold 0.5)",
"Specificity (threshold 0.5)","Balanced Accuracy","F1-Score",
"Hosmer-Lemeshow p-value"),
`v1 — Full (14 var)` = c(
length(contamination_fit1$coefficients) - 1, round(AIC(contamination_fit1), 3),
round(DescTools::PseudoR2(contamination_fit1, which = "Nagelkerke"), 4),
round(auc_value(roc_points(train_data$kontaminasi_biner, p_train_v1)), 3),
round(auc_test_v1, 3), round(met_v1_full$accuracy, 3),
round(met_v1_full$sensitivity, 3), round(met_v1_full$specificity, 3),
round(met_v1_full$balanced_accuracy, 3), round(met_v1_full$f1_score, 3),
round(ResourceSelection::hoslem.test(train_data$kontaminasi_biner,
fitted(contamination_fit1), g = 10)$p.value, 4)
),
`v2 — Full diperbaiki (12 var)` = c(
length(contamination_fit2$coefficients) - 1, round(AIC(contamination_fit2), 3),
round(DescTools::PseudoR2(contamination_fit2, which = "Nagelkerke"), 4),
round(auc_train_v2, 3), round(auc_test_v2, 3),
round(met_v2_full$accuracy, 3), round(met_v2_full$sensitivity, 3),
round(met_v2_full$specificity, 3), round(met_v2_full$balanced_accuracy, 3),
round(met_v2_full$f1_score, 3),
round(ResourceSelection::hoslem.test(train_data$kontaminasi_biner,
fitted(contamination_fit2), g = 10)$p.value, 4)
),
`v3 — Reduced / Final (8 var)` = c(
length(contamination_fit3$coefficients) - 1, round(AIC(contamination_fit3), 3),
round(DescTools::PseudoR2(contamination_fit3, which = "Nagelkerke"), 4),
round(auc_train_v3, 3), round(auc_test_v3, 3),
round(met_v3_full$accuracy, 3), round(met_v3_full$sensitivity, 3),
round(met_v3_full$specificity, 3), round(met_v3_full$balanced_accuracy, 3),
round(met_v3_full$f1_score, 3),
round(ResourceSelection::hoslem.test(train_data$kontaminasi_biner,
fitted(contamination_fit3), g = 10)$p.value, 4)
),
check.names = FALSE
) %>%
knitr::kable(caption = "Perbandingan Performa Lengkap: v1 vs v2 vs v3") %>%
kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE)| Metrik | v1 — Full (14 var) | v2 — Full diperbaiki (12 var) | v3 — Reduced / Final (8 var) |
|---|---|---|---|
| Jumlah Variabel | 38.0000 | 36.0000 | 25.0000 |
| AIC | 1489.0600 | 1493.2570 | 1477.2870 |
| Nagelkerke R² | 0.5029 | 0.4986 | 0.4954 |
| AUC Training | 0.8630 | 0.8600 | 0.8590 |
| AUC Testing | 0.8410 | 0.8370 | 0.8470 |
| Akurasi (threshold 0.5) | 0.7490 | 0.7440 | 0.7540 |
| Sensitivity (threshold 0.5) | 0.7970 | 0.7930 | 0.7880 |
| Specificity (threshold 0.5) | 0.6770 | 0.6710 | 0.7020 |
| Balanced Accuracy | 0.7370 | 0.7320 | 0.7450 |
| F1-Score | 0.7920 | 0.7880 | 0.7930 |
| Hosmer-Lemeshow p-value | 0.0185 | 0.0740 | 0.1280 |
8.2 Perbandingan Kurva ROC: v1 vs v2 vs v3
roc_compare_all <- bind_rows(
roc_test_v1 %>% mutate(data = "Testing v1 (Full — 14 var)"),
roc_test_v2_plot %>% mutate(data = "Testing v2 (Full diperbaiki — 12 var)"),
roc_test_v3 %>% mutate(data = "Testing v3 (Reduced — 8 var)")
)
ggplot(roc_compare_all, aes(x = fpr, y = sensitivity, color = data)) +
geom_path(linewidth = 1.2) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#6c757d") +
coord_equal() +
scale_color_manual(values = c(
"Testing v1 (Full — 14 var)" = "#86c287",
"Testing v2 (Full diperbaiki — 12 var)" = "#eab308",
"Testing v3 (Reduced — 8 var)" = "#2d6a2d"
)) +
labs(title = "Perbandingan Kurva ROC: v1 vs v2 vs v3 (Data Testing)",
subtitle = paste0("AUC v1 = ", round(auc_test_v1, 3),
" | AUC v2 = ", round(auc_test_v2, 3),
" | AUC v3 = ", round(auc_test_v3, 3)),
x = "False Positive Rate (1 - Specificity)",
y = "Sensitivity (True Positive Rate)", color = "Model") +
theme_minimal(base_size = 12) + theme(legend.position = "bottom")
Perbandingan kurva ROC: Ketiga model menunjukkan kemampuan klasifikasi yang baik dengan nilai AUC di atas 0,84. Model v3 memiliki nilai AUC tertinggi pada data testing (0,847), sedikit lebih tinggi dibandingkan model v2 (0,837). Hasil ini menunjukkan bahwa penyederhanaan model melalui penghapusan variabel yang kurang informatif tidak menurunkan performa model, bahkan memberikan peningkatan kecil pada kemampuan model dalam membedakan lokasi terkontaminasi dan tidak terkontaminasi.
9 Hasil & Interpretasi Model Final
9.1 Persamaan Model Final (v3)
\[ \widehat{\text{logit}}(\pi_i) = 72{,}13 + 0{,}129\,(\text{durasi operasi}) - 0{,}037\,(\text{tahun mulai}) \] \[ + 0{,}438\,(\text{jumlah pelanggaran}) + 0{,}0000926\,(\text{curah hujan}) + 0{,}020\,(\text{kecepatan angin}) \] \[ - 0{,}034\,(\text{tekstur tanah}) - 0{,}903\,(\text{log mobilitas polutan}) \] \[ + \beta_8\,(\text{jenis industri}) \]
9.2 Odds Ratio — Karakteristik Operasional Industri
Jumlah Pelanggaran Regulasi Lingkungan (OR = 1,549; CI 95%: 1,422–1,687; p < 0,001)
Prediktor paling kuat dalam model. Setiap tambahan satu pelanggaran regulasi lingkungan meningkatkan odds terjadinya kontaminasi sebesar 54,9%, dengan asumsi variabel lain tetap konstan. Hal ini menunjukkan bahwa efektivitas penegakan hukum lingkungan memiliki peran penting dalam menurunkan risiko kontaminasi.
Durasi Operasi (OR = 1,137; p < 0,001)
Setiap tambahan satu tahun masa operasi industri meningkatkan odds kontaminasi sebesar 13,7%. Hal ini berkaitan dengan akumulasi polutan dalam jangka panjang di area industri.
Tahun Mulai Beroperasi (OR = 0,963; p = 0,002)
Variabel ini menunjukkan hubungan negatif yang signifikan dengan status kontaminasi lahan. Setiap kenaikan satu tahun pada tahun mulai beroperasi dikaitkan dengan penurunan odds terjadinya kontaminasi sebesar 3,7%, dengan asumsi variabel lain konstan. Temuan ini mengindikasikan bahwa industri yang mulai beroperasi pada periode yang lebih baru cenderung memiliki risiko kontaminasi yang lebih rendah, yang kemungkinan berkaitan dengan peningkatan standar pengelolaan lingkungan serta pengetatan regulasi lingkungan pada periode yang lebih modern.
Jenis Industri (baseline: Industri Lainnya)
Sektor yang melibatkan proses kimia intensif dan peleburan logam memiliki risiko kontaminasi tertinggi:
Jenis Industri OR CI 95% p-value Organic chemical raw material mfg. 2,899 1,569–5,359 < 0,001 *** Crude oil processing & petroleum mfg. 3,396 1,707–6,758 < 0,001 *** Lead-zinc smelting 2,783 1,282–6,044 0,010 ** Pesticide manufacturing 1,896 1,036–3,471 0,038 * Coking 1,969 1,006–3,854 0,048 * Sebaliknya, sektor berikut menunjukkan risiko lebih rendah dibanding Industri Lainnya:
Jenis Industri OR Interpretasi Hazardous waste treatment 0,077 Peluang kontaminasi sekitar 92,3% lebih rendah dibandingkan industri lainnya. Cotton printing & dyeing 0,353 Peluang kontaminasi sekitar 64,7% lebih rendah dibandingkan industri lainnya.
9.3 Odds Ratio — Kondisi Lingkungan Biofisik
Curah Hujan (OR ≈ 1,000 per 1 mm; p < 0,001)
Nilai OR yang mendekati 1 untuk kenaikan 1 mm curah hujan menunjukkan bahwa perubahan kecil pada curah hujan tidak memberikan perubahan odds yang besar secara langsung. Hal ini dipengaruhi oleh penggunaan satuan yang sangat kecil, sehingga dampaknya lebih terlihat ketika perubahan curah hujan terjadi dalam skala yang lebih besar. Secara lingkungan, peningkatan curah hujan berpotensi meningkatkan risiko kontaminasi melalui proses infiltrasi, pencucian (leaching), dan limpasan permukaan yang dapat memperluas penyebaran polutan.
Kecepatan Angin (OR = 1,020; CI 95%: 1,005–1,035; p = 0,010)
Setiap kenaikan satu satuan kecepatan angin dikaitkan dengan peningkatan odds kontaminasi sebesar 2,0%, dengan asumsi variabel lain konstan.
Tekstur Tanah — Kandungan Liat (OR = 0,967; CI 95%: 0,947–0,987; p = 0,001)
Setiap peningkatan kandungan liat sebesar 1% menurunkan odds kontaminasi sebesar 3,3%. Hubungan protektif ini berkaitan dengan karakteristik fisik tanah berliat dalam membatasi pergerakan polutan.
Log Mobilitas Polutan (OR = 0,405; CI 95%: 0,226–0,728; p = 0,002)
Variabel ini menunjukkan hubungan negatif dengan status kontaminasi lahan. Setiap kenaikan satu unit pada log mobilitas polutan dikaitkan dengan penurunan odds terjadinya kontaminasi sebesar 59,5%, dengan asumsi variabel lain konstan. Temuan ini tidak serta-merta menunjukkan bahwa polutan dengan mobilitas tinggi lebih aman. Sebaliknya, polutan dengan mobilitas tinggi cenderung berpindah dari lapisan permukaan ke lapisan lebih dalam atau terbawa ke badan air, sehingga konsentrasinya pada titik sampling permukaan menjadi lebih rendah. Hasil ini lebih tepat dipahami sebagai indikasi perbedaan perilaku migrasi dan deteksi polutan.
10 Kesimpulan
1. Faktor risiko yang berpengaruh signifikan terhadap kontaminasi lahan industri Dari 14 variabel kandidat, terdapat 8 faktor yang terbukti berpengaruh signifikan terhadap probabilitas kontaminasi lahan industri di China. Dari sisi operasional, jumlah pelanggaran regulasi lingkungan merupakan prediktor terkuat (OR = 1,549), yang menunjukkan bahwa setiap tambahan satu pelanggaran dikaitkan dengan peningkatan odds kontaminasi sebesar 54,9%. Durasi operasi yang lebih lama juga meningkatkan risiko kontaminasi (OR = 1,137), sedangkan industri yang mulai beroperasi pada periode yang lebih akhir cenderung memiliki risiko yang lebih rendah (OR = 0,963), sejalan dengan perkembangan teknologi dan regulasi lingkungan yang semakin ketat. Dari sisi jenis industri, sektor crude oil processing & petroleum manufacturing (OR = 3,396), organic chemical raw material manufacturing (OR = 2,899), lead-zinc smelting (OR = 2,783), dan coking (OR = 1,969) menunjukkan risiko kontaminasi tertinggi dibanding industri lainnya. Dari sisi lingkungan biofisik, curah hujan dan kecepatan angin berasosiasi positif dengan kontaminasi, sedangkan kandungan liat tanah menunjukkan efek protektif.
2. Model prediktif regresi logistik biner yang parsimonious dengan performa baik Model final (v3) dengan 8 prediktor diperoleh melalui proses penyederhanaan model secara bertahap dengan mempertimbangkan multikolinearitas, linearitas logit, signifikansi variabel, dan goodness-of-fit. Hasil evaluasi menunjukkan bahwa model telah sesuai dengan data (Hosmer–Lemeshow p = 0,128) dan memiliki kemampuan yang baik dalam membedakan lokasi terkontaminasi dan tidak terkontaminasi (AUC testing = 0,847). Selain itu, nilai akurasi sebesar 75,4% dan Nagelkerke R² sebesar 0,495 menunjukkan bahwa model memiliki kemampuan prediksi yang baik serta mampu menjelaskan sekitar 49,5% variasi status kontaminasi lahan. Selisih AIC antara v2 (1493,257) dan v3 (1477,287) sebesar 15,97 juga mendukung pemilihan model yang lebih parsimoni.
3. Interpretasi pengaruh faktor melalui Odds Ratio Nilai Odds Ratio dari model final memberikan gambaran mengenai arah dan besaran pengaruh masing-masing faktor terhadap kontaminasi lahan. Hasil penelitian menunjukkan bahwa industri kimia dan metalurgi, terutama yang memiliki riwayat pelanggaran lingkungan dan telah beroperasi dalam jangka waktu panjang, merupakan kelompok yang perlu diprioritaskan dalam kegiatan pengawasan dan remediasi. Sementara itu, hubungan negatif yang ditemukan pada variabel mobilitas polutan (OR = 0,405) perlu diinterpretasikan secara hati-hati. Polutan dengan mobilitas tinggi cenderung bermigrasi ke lapisan tanah yang lebih dalam, sehingga konsentrasinya pada titik sampling permukaan menjadi lebih rendah. Dengan demikian, hasil ini lebih mencerminkan keterbatasan metode pengukuran daripada rendahnya risiko aktual.
Delapan faktor risiko signifikan:
jumlah_pelanggaran— prediktor terkuat (OR = 1,549; p < 0,001)durasi_operasi(OR = 1,137; p < 0,001)tahun_mulai(OR = 0,963; p = 0,002)jenis_industri— beberapa level signifikan, tertinggi: Crude oil & petroleum mfg. (OR = 3,396)curah_hujan(p < 0,001)kecepatan_angin(OR = 1,020; p = 0,010)tekstur_tanah— efek protektif (OR = 0,967; p = 0,001)log_mobilitas_polutan(OR = 0,405; p = 0,002)
Implikasi praktis: Implikasi praktis dari hasil penelitian ini menunjukkan bahwa prioritas pemantauan dan remediasi perlu difokuskan pada industri kimia dan metalurgi (organic chemical manufacturing, steel smelting, coking, dan pesticide manufacturing) yang memiliki riwayat pelanggaran lingkungan tinggi serta telah beroperasi dalam jangka waktu lama. Selain itu, hubungan negatif pada variabel mobilitas polutan mengindikasikan bahwa metode pemantauan yang hanya berbasis sampel tanah permukaan berpotensi tidak sepenuhnya menangkap risiko kontaminasi. Oleh karena itu, diperlukan pengembangan sistem pemantauan yang lebih komprehensif dengan mempertimbangkan karakteristik migrasi dan distribusi polutan di lingkungan.
Dataset: Urban Contaminated Sites along China’s Urbanization — Li, K. (2023). Mendeley Data. DOI: 10.17632/r4y2vcpfmx.1