1. Persiapan Data

1. 1 Instalasi & load package
required <- c("tidyverse","caret","e1071","randomForest","rpart","rpart.plot",
"mgcv","Metrics","gridExtra","readr","scales")
new <- required[!(required %in% installed.packages()[,"Package"])]
if(length(new)) install.packages(new)
invisible(lapply(required, library, character.only = TRUE))
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.2
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: lattice
## 
## 
## Attaching package: 'caret'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     lift
## 
## 
## 
## Attaching package: 'e1071'
## 
## 
## The following object is masked from 'package:ggplot2':
## 
##     element
## 
## 
## randomForest 4.7-1.2
## 
## Type rfNews() to see new features/changes/bug fixes.
## 
## 
## Attaching package: 'randomForest'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
## 
## 
## Loading required package: nlme
## 
## 
## Attaching package: 'nlme'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## 
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
## 
## 
## Attaching package: 'Metrics'
## 
## 
## The following objects are masked from 'package:caret':
## 
##     precision, recall
## 
## 
## 
## Attaching package: 'gridExtra'
## 
## 
## The following object is masked from 'package:randomForest':
## 
##     combine
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## 
## 
## Attaching package: 'scales'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
1.2 Set seed & path file
set.seed(123)
train_path <- "D:/Kuliah/Statistika Lingkungan/kualitasair.xlsx - Training (1).csv"
test_path  <- "D:/Kuliah/Statistika Lingkungan/kualitasair.xlsx - Testing (1).csv"
1.3 Baca file CSV & tampilkan 5 baris pertama
train_df <- read_csv(train_path)
## Rows: 300 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Lokasi, Status
## dbl (5): pH, DO, BOD, TSS, Suhu
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
test_df  <- read_csv(test_path)
## Rows: 75 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Lokasi
## dbl (5): pH, DO, BOD, TSS, Suhu
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
cat("Training data:", nrow(train_df), "baris\n")
## Training data: 300 baris
cat("Testing data:", nrow(test_df), "baris\n")
## Testing data: 75 baris
glimpse(train_df)
## Rows: 300
## Columns: 7
## $ Lokasi <chr> "S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9", "S10", "S…
## $ pH     <dbl> 7.6855, 6.7177, 7.1816, 7.3164, 7.2021, 6.9469, 7.7558, 6.9527,…
## $ DO     <dbl> NA, 5.7236, 4.8906, 6.1339, 7.7853, 8.4222, 4.9232, 6.4859, 7.3…
## $ BOD    <dbl> 1.7136, 1.4402, 2.7274, 3.1398, 1.1778, 3.2324, NA, 4.0358, 3.1…
## $ TSS    <dbl> 43.1415, 44.2963, NA, 41.0104, 48.0967, 48.5610, 49.0343, 51.81…
## $ Suhu   <dbl> 26.7972, 27.7284, 26.0255, 29.6639, 26.4099, 28.6809, 29.7409, …
## $ Status <chr> "Tercemar ringan", "Tercemar ringan", "Tercemar ringan", "Terce…
print("===== 5 BARIS PERTAMA DATA TRAINING =====")
## [1] "===== 5 BARIS PERTAMA DATA TRAINING ====="
print(head(train_df, 5))
## # A tibble: 5 × 7
##   Lokasi    pH    DO   BOD   TSS  Suhu Status         
##   <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <chr>          
## 1 S1      7.69 NA     1.71  43.1  26.8 Tercemar ringan
## 2 S2      6.72  5.72  1.44  44.3  27.7 Tercemar ringan
## 3 S3      7.18  4.89  2.73  NA    26.0 Tercemar ringan
## 4 S4      7.32  6.13  3.14  41.0  29.7 Tercemar ringan
## 5 S5      7.20  7.79  1.18  48.1  26.4 baik

2. Data Cleaning & Eksplorasi (Soal 1)

Perubahan sudah diterapkan: bagian Data Cleaning & Eksplorasi (Soal 1) kini mencakup cek missing, imputasi sederhana, deteksi outlier & winsorize, serta standarisasi Status. Nama kolom juga dinormalisasi dan variabel numerik diidentifikasi.

2.2 Normalisasi nama kolom & identifikasi variabel numerik
# Normalisasi nama kolom (hapus spasi berlebih dan ganti dengan underscore)
names(train_df) <- trimws(names(train_df))
names(train_df) <- gsub("\\s+", "_", names(train_df))

names(test_df) <- trimws(names(test_df))
names(test_df) <- gsub("\\s+", "_", names(test_df))

# Deteksi variabel numerik

num_vars <- intersect(c("pH","DO","BOD","TSS","Suhu"), names(train_df))
cat("Variabel numerik yang terdeteksi:", paste(num_vars, collapse = ", "), "\n")
## Variabel numerik yang terdeteksi: pH, DO, BOD, TSS, Suhu
2.3 Cek Missing Value (Training & Testing)
# Cek missing value di training
cat("===== Missing value per kolom (TRAINING) =====\n")
## ===== Missing value per kolom (TRAINING) =====
missing_train <- sapply(train_df, function(x) sum(is.na(x)))
print(missing_train)
## Lokasi     pH     DO    BOD    TSS   Suhu Status 
##      0      0     23     22     24      0      0
# Cek missing value di testing
cat("\n===== Missing value per kolom (TESTING) =====\n")
## 
## ===== Missing value per kolom (TESTING) =====
missing_test <- sapply(test_df, function(x) sum(is.na(x)))
print(missing_test)
## Lokasi     pH     DO    BOD    TSS   Suhu 
##      0      0      8      9      7      0
2.4 Tampilkan baris yang memiliki missing (preview)
# Training
cat("\n===== Baris dengan Missing (TRAINING) =====\n")
## 
## ===== Baris dengan Missing (TRAINING) =====
rows_with_na_train <- train_df[rowSums(is.na(train_df)) > 0, ]
cat("Jumlah baris dengan minimal 1 NA di TRAINING:", nrow(rows_with_na_train), "\n")
## Jumlah baris dengan minimal 1 NA di TRAINING: 60
if(nrow(rows_with_na_train) > 0) print(head(rows_with_na_train, 10)) else print("Tidak ada baris dengan NA di training.")
## # A tibble: 10 × 7
##    Lokasi    pH    DO   BOD   TSS  Suhu Status         
##    <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <chr>          
##  1 S1      7.69 NA     1.71  43.1  26.8 Tercemar ringan
##  2 S3      7.18  4.89  2.73  NA    26.0 Tercemar ringan
##  3 S7      7.76  4.92 NA     49.0  29.7 Tercemar ringan
##  4 S10     6.97  5.80 NA     NA    23.8 BAIK           
##  5 S16     7.32 NA     2.66  35.9  28.6 Baik           
##  6 S28     6.12 NA     3.28  56.2  27.3 Tercemar ringan
##  7 S43     7.38 NA     3.74  42.0  23.4 Tercemar ringan
##  8 S45     6.32 NA     3.04  56.8  27.9 Tercemar ringan
##  9 S47     6.59  4.55  1.17  NA    27.2 Tercemar ringan
## 10 S48     7.72  6.10 NA     51.8  30.0 Baik
# Testing
cat("\n===== Baris dengan Missing (TESTING) =====\n")
## 
## ===== Baris dengan Missing (TESTING) =====
rows_with_na_test <- test_df[rowSums(is.na(test_df)) > 0, ]
cat("Jumlah baris dengan minimal 1 NA di TESTING:", nrow(rows_with_na_test), "\n")
## Jumlah baris dengan minimal 1 NA di TESTING: 21
if(nrow(rows_with_na_test) > 0) print(head(rows_with_na_test, 10)) else print("Tidak ada baris dengan NA di testing.")
## # A tibble: 10 × 6
##    Lokasi    pH    DO   BOD   TSS  Suhu
##    <chr>  <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 S303    7.02  6.59  3.00  NA    26.6
##  2 S304    7.37 NA     3.50  39.0  26.7
##  3 S312    6.62  5.42 NA     51.1  29.0
##  4 S318    7.08 NA    NA     59.4  25.4
##  5 S319    6.38  5.19 NA     40.9  26.3
##  6 S322    7.09  5.20 NA     50.1  26.7
##  7 S323    7.03  3.81 NA     51.2  26.6
##  8 S324    7.00  5.71  3.13  NA    26.2
##  9 S326    6.59 NA     3.73  51.2  27.5
## 10 S329    6.58  5.67 NA     47.1  28.0
2.5 Tangani Missing Value (Imputasi Median untuk Numerik)

Median dipilih sebagai metode imputasi karena nilai median tidak terpengaruh oleh outlier dan distribusi data yang skewed. Dengan demikian, imputasi menggunakan median menghasilkan estimasi yang lebih stabil dan representatif untuk variabel numerik dibanding mean.

# Imputasi median pada TRAINING
cat("===== Imputasi Missing Value (TRAINING) =====\n")
## ===== Imputasi Missing Value (TRAINING) =====
for(v in num_vars){
  if (v %in% names(train_df)) {
    if(any(is.na(train_df[[v]]))){
      med <- median(train_df[[v]], na.rm = TRUE)
      train_df[[v]][is.na(train_df[[v]])] <- med
      cat("Imputasi median untuk", v, "=", med, "\n")
    }
  }
}
## Imputasi median untuk DO = 5.9909 
## Imputasi median untuk BOD = 3.0661 
## Imputasi median untuk TSS = 49.52205
if("Status" %in% names(train_df)){
  cat("Nilai NA pada Status (TRAINING):", sum(is.na(train_df$Status)), "\n")
}
## Nilai NA pada Status (TRAINING): 0
# Imputasi median pada TESTING
cat("\n===== Imputasi Missing Value (TESTING) =====\n")
## 
## ===== Imputasi Missing Value (TESTING) =====
for(v in num_vars){
  if (v %in% names(test_df)) {
    if(any(is.na(test_df[[v]]))){
      med <- median(test_df[[v]], na.rm = TRUE)
      test_df[[v]][is.na(test_df[[v]])] <- med
      cat("Imputasi median untuk", v, "=", med, "\n")
    }
  }
}
## Imputasi median untuk DO = 5.644 
## Imputasi median untuk BOD = 3.06165 
## Imputasi median untuk TSS = 49.56945
if("Status" %in% names(test_df)){
  cat("Nilai NA pada Status (TESTING):", sum(is.na(test_df$Status)), "\n")
}
2.6 Deteksi Outlier (IQR Rule)
cat("\n===== Deteksi Outlier (TRAINING) =====\n")
## 
## ===== Deteksi Outlier (TRAINING) =====
outlier_train <- sapply(num_vars, function(v){
  x <- train_df[[v]]
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  sum((x < (Q1 - 1.5*IQR)) | (x > (Q3 + 1.5*IQR)), na.rm = TRUE)
})
print(outlier_train)
##   pH   DO  BOD  TSS Suhu 
##    4    4    5    5    2
cat("\n===== Deteksi Outlier (TESTING) =====\n")
## 
## ===== Deteksi Outlier (TESTING) =====
outlier_test <- sapply(num_vars, function(v){
  x <- test_df[[v]]
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  sum((x < (Q1 - 1.5*IQR)) | (x > (Q3 + 1.5*IQR)), na.rm = TRUE)
})
print(outlier_test)
##   pH   DO  BOD  TSS Suhu 
##    2    3    2    2    1
2.6 Winsorizing

Winsorizing = batasi nilai ekstrem (outlier) ke batas tertentu supaya tidak mengacaukan analisis, tapi data tetap dipertahankan.

# Winsorize pada TRAINING
for(v in num_vars){
  x <- train_df[[v]]
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  lower <- Q1 - 1.5*IQR
  upper <- Q3 + 1.5*IQR
  train_df[[v]] <- pmin(pmax(x, lower), upper)
}

# Winsorize pada TESTING
for(v in num_vars){
  x <- test_df[[v]]
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  lower <- Q1 - 1.5*IQR
  upper <- Q3 + 1.5*IQR
  test_df[[v]] <- pmin(pmax(x, lower), upper)
}
cat("Winsorizing selesai untuk variabel numerik.\n")
## Winsorizing selesai untuk variabel numerik.
2.7 Standarisasi Penulisan Kategori Status
standardize_status <- function(df){
  if("Status" %in% names(df)){
    df$Status <- as.character(df$Status)
    df$Status <- trimws(tolower(df$Status))
    df$Status[df$Status %in% c("1","baik","baik=1","baik = 1")] <- "Baik"
    df$Status[df$Status %in% c("2","tercemar ringan","tercemar_ringan","tercemar ringan=2")] <- "Tercemar_Ringan"
    df$Status[df$Status %in% c("3","tercemar berat","tercemar_berat","tercemar berat=3")] <- "Tercemar_Berat"
    df$Status <- factor(df$Status, levels = c("Baik","Tercemar_Ringan","Tercemar_Berat"))
  }
  return(df)
}

train_df <- standardize_status(train_df)
test_df  <- standardize_status(test_df)

cat("Distribusi Status TRAINING:\n")
## Distribusi Status TRAINING:
print(table(train_df$Status, useNA = "always"))
## 
##            Baik Tercemar_Ringan  Tercemar_Berat            <NA> 
##              72             221               7               0
cat("\nDistribusi Status TESTING:\n")
## 
## Distribusi Status TESTING:
print(table(test_df$Status, useNA = "always"))
## Warning: Unknown or uninitialised column: `Status`.
## 
## <NA> 
##    0
2.8 Ringkasan statistik (setelah pembersihan)
cat("\n===== Ringkasan Statistik Numerik (TRAINING) =====\n")
## 
## ===== Ringkasan Statistik Numerik (TRAINING) =====
print(summary(train_df[num_vars]))
##        pH              DO             BOD              TSS       
##  Min.   :5.697   Min.   :3.615   Min.   :0.8513   Min.   :27.28  
##  1st Qu.:6.670   1st Qu.:5.413   1st Qu.:2.4599   1st Qu.:44.28  
##  Median :6.988   Median :5.991   Median :3.0661   Median :49.52  
##  Mean   :6.990   Mean   :5.977   Mean   :3.0041   Mean   :49.68  
##  3rd Qu.:7.318   3rd Qu.:6.611   3rd Qu.:3.5323   3rd Qu.:55.62  
##  Max.   :8.290   Max.   :8.409   Max.   :5.1409   Max.   :72.62  
##       Suhu      
##  Min.   :22.77  
##  1st Qu.:26.62  
##  Median :28.01  
##  Mean   :28.12  
##  3rd Qu.:29.46  
##  Max.   :33.73
cat("\n===== Ringkasan Statistik Numerik (TESTING) =====\n")
## 
## ===== Ringkasan Statistik Numerik (TESTING) =====
print(summary(test_df[num_vars]))
##        pH              DO             BOD             TSS       
##  Min.   :5.514   Min.   :4.005   Min.   :1.368   Min.   :24.95  
##  1st Qu.:6.615   1st Qu.:5.354   1st Qu.:2.574   1st Qu.:43.35  
##  Median :6.965   Median :5.644   Median :3.062   Median :49.57  
##  Mean   :6.972   Mean   :5.801   Mean   :3.034   Mean   :49.16  
##  3rd Qu.:7.349   3rd Qu.:6.253   3rd Qu.:3.378   3rd Qu.:55.62  
##  Max.   :8.449   Max.   :7.602   Max.   :4.584   Max.   :74.02  
##       Suhu      
##  Min.   :23.39  
##  1st Qu.:26.79  
##  Median :28.28  
##  Mean   :28.37  
##  3rd Qu.:29.88  
##  Max.   :34.52

Data numerik sudah bersih, outlier sudah dikurangi, missing value telah diimputasi menggunakan median, dan variabel Status telah distandarisasi menjadi tiga kategori: Baik, Tercemar_ringan, Tercemar_berat. Distribusi numerik pada training dan testing cukup serupa, menandakan dataset konsisten untuk model selanjutnya.

3. Klasifikasi Status Kualitas Air (Soal 2)

3.1 Persiapan Data & Split Internal
features <- intersect(c("pH","DO","BOD","TSS","Suhu"), names(train_df))
cat("Fitur numerik yang digunakan:", paste(features, collapse=", "), "\n")
## Fitur numerik yang digunakan: pH, DO, BOD, TSS, Suhu
# Pastikan target sebagai faktor
train_df$Status <- factor(train_df$Status, levels = c("Baik","Tercemar_Ringan","Tercemar_Berat"))

Model nanti akan belajar hubungan antara lima variabel fisik/kimia air dengan kategori status kualitas air. Output menunjukkan dataset sudah siap untuk tahap modeling.

3.2 Preprocessing Data (Center & Scale + Split Internal)

Kode ini memecah data menjadi train (80%) dan validation (20%), memastikan semua fitur numerik, lalu menormalkan (center & scale) berdasarkan train saja supaya model belajar dari data training tanpa bocor ke validation atau testing.

# Split train -> train/validation (80/20)

train_index <- caret::createDataPartition(train_df$Status, p = 0.8, list = FALSE)
train_set <- train_df[train_index, ]
val_set   <- train_df[-train_index, ]
cat("Baris train:", nrow(train_set), " | Baris validation:", nrow(val_set), "\n")
## Baris train: 241  | Baris validation: 59
# Pastikan semua kolom numerik

for(v in features){
train_set[[v]] <- as.numeric(train_set[[v]])
val_set[[v]]   <- as.numeric(val_set[[v]])
test_df[[v]]   <- as.numeric(test_df[[v]])
}

# Buat objek preprocessing hanya dari training

preProc <- caret::preProcess(train_set[, features], method = c("center","scale"))
train_scaled <- predict(preProc, train_set[, features]) %>% as.data.frame()
val_scaled   <- predict(preProc, val_set[, features]) %>% as.data.frame()
test_scaled  <- predict(preProc, test_df[, features]) %>% as.data.frame()

# Tambahkan kembali target

train_scaled$Status <- train_set$Status
val_scaled$Status   <- val_set$Status

sebagian besar data digunakan untuk melatih model (train), sisanya untuk mengevaluasi kinerja model (validation).

3.3 Model 1 — SVM (Radial Kernel)

Model SVM (Support Vector Machine) digunakan untuk mengklasifikasikan status kualitas air berdasarkan fitur numerik. Kernel radial dipilih agar model bisa menangkap hubungan non-linear antar variabel.

set.seed(123)
svm_model <- e1071::svm(Status ~ ., data = train_scaled, kernel = "radial", probability = TRUE)
svm_pred  <- predict(svm_model, val_scaled)
conf_svm  <- caret::confusionMatrix(svm_pred, val_scaled$Status)
cat("\n=== CONFUSION MATRIX: SVM ===\n")
## 
## === CONFUSION MATRIX: SVM ===
print(conf_svm)
## Confusion Matrix and Statistics
## 
##                  Reference
## Prediction        Baik Tercemar_Ringan Tercemar_Berat
##   Baik               8               1              0
##   Tercemar_Ringan    6              43              1
##   Tercemar_Berat     0               0              0
## 
## Overall Statistics
##                                           
##                Accuracy : 0.8644          
##                  95% CI : (0.7502, 0.9396)
##     No Information Rate : 0.7458          
##     P-Value [Acc > NIR] : 0.02096         
##                                           
##                   Kappa : 0.5913          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: Baik Class: Tercemar_Ringan Class: Tercemar_Berat
## Sensitivity               0.5714                 0.9773               0.00000
## Specificity               0.9778                 0.5333               1.00000
## Pos Pred Value            0.8889                 0.8600                   NaN
## Neg Pred Value            0.8800                 0.8889               0.98305
## Prevalence                0.2373                 0.7458               0.01695
## Detection Rate            0.1356                 0.7288               0.00000
## Detection Prevalence      0.1525                 0.8475               0.00000
## Balanced Accuracy         0.7746                 0.7553               0.50000

Model SVM digunakan untuk mengklasifikasikan status kualitas air dengan kernel radial agar dapat menangkap hubungan non-linear antar variabel. Akurasi keseluruhan sekitar 86%, menunjukkan model cukup baik, terutama untuk kelas Tercemar_Ringan. Namun, sensitivitas untuk kelas Baik rendah (57%) dan Tercemar_Berat tidak terdeteksi sama sekali, menandakan model kesulitan memprediksi kelas minoritas. Kappa 0.59 menunjukkan kesepakatan moderat antara prediksi dan aktual.

3.4 Model 2 — Decision Tree

Decision Tree adalah model klasifikasi berbentuk pohon yang memprediksi kelas dengan membagi data berdasarkan aturan sederhana dari variabel input. Model ini mudah diinterpretasi dan menangani data numerik maupun kategorik.

set.seed(123)
tree_model <- rpart::rpart(Status ~ ., data = train_scaled, method = "class")
tree_pred  <- predict(tree_model, val_scaled, type="class")
conf_tree  <- caret::confusionMatrix(tree_pred, val_scaled$Status)
cat("\n=== CONFUSION MATRIX: Decision Tree ===\n")
## 
## === CONFUSION MATRIX: Decision Tree ===
print(conf_tree)
## Confusion Matrix and Statistics
## 
##                  Reference
## Prediction        Baik Tercemar_Ringan Tercemar_Berat
##   Baik              13               0              0
##   Tercemar_Ringan    1              44              1
##   Tercemar_Berat     0               0              0
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9661          
##                  95% CI : (0.8829, 0.9959)
##     No Information Rate : 0.7458          
##     P-Value [Acc > NIR] : 6.696e-06       
##                                           
##                   Kappa : 0.9075          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: Baik Class: Tercemar_Ringan Class: Tercemar_Berat
## Sensitivity               0.9286                 1.0000               0.00000
## Specificity               1.0000                 0.8667               1.00000
## Pos Pred Value            1.0000                 0.9565                   NaN
## Neg Pred Value            0.9783                 1.0000               0.98305
## Prevalence                0.2373                 0.7458               0.01695
## Detection Rate            0.2203                 0.7458               0.00000
## Detection Prevalence      0.2203                 0.7797               0.00000
## Balanced Accuracy         0.9643                 0.9333               0.50000

Model Decision Tree pada dataset ini mencapai akurasi sekitar 96,6% dengan Kappa 0,91, menunjukkan kesepakatan yang sangat baik antara prediksi dan nilai aktual. Sensitivitas dan spesifikasi tinggi untuk kelas Baik dan Tercemar_Ringan menandakan model mampu mengenali dua kelas dominan dengan baik. Namun, kelas Tercemar_Berat masih tidak terdeteksi karena jumlah observasi minoritas sangat sedikit, sehingga model kesulitan mempelajari pola untuk kelas ini.

3.5 Model 3 — Random Forest

Random Forest adalah ensemble dari banyak decision tree yang digabungkan untuk meningkatkan akurasi dan mengurangi overfitting. Setiap pohon dilatih pada subset data dan fitur acak, kemudian prediksi akhir ditentukan lewat majority voting.

set.seed(123)
rf_model <- randomForest::randomForest(Status ~ ., data = train_scaled, ntree=500, importance=TRUE)
rf_pred  <- predict(rf_model, val_scaled)
conf_rf  <- caret::confusionMatrix(rf_pred, val_scaled$Status)
cat("\n=== CONFUSION MATRIX: Random Forest ===\n")
## 
## === CONFUSION MATRIX: Random Forest ===
print(conf_rf)
## Confusion Matrix and Statistics
## 
##                  Reference
## Prediction        Baik Tercemar_Ringan Tercemar_Berat
##   Baik              13               0              0
##   Tercemar_Ringan    1              44              1
##   Tercemar_Berat     0               0              0
## 
## Overall Statistics
##                                           
##                Accuracy : 0.9661          
##                  95% CI : (0.8829, 0.9959)
##     No Information Rate : 0.7458          
##     P-Value [Acc > NIR] : 6.696e-06       
##                                           
##                   Kappa : 0.9075          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: Baik Class: Tercemar_Ringan Class: Tercemar_Berat
## Sensitivity               0.9286                 1.0000               0.00000
## Specificity               1.0000                 0.8667               1.00000
## Pos Pred Value            1.0000                 0.9565                   NaN
## Neg Pred Value            0.9783                 1.0000               0.98305
## Prevalence                0.2373                 0.7458               0.01695
## Detection Rate            0.2203                 0.7458               0.00000
## Detection Prevalence      0.2203                 0.7797               0.00000
## Balanced Accuracy         0.9643                 0.9333               0.50000

Model Random Forest pada dataset ini mencapai akurasi sekitar 96,6% dengan Kappa 0,91, menunjukkan kesepakatan yang sangat baik antara prediksi dan nilai aktual. Model mampu mengklasifikasikan kelas Baik dan Tercemar_Ringan dengan sensitivitas dan spesifikasi tinggi, menandakan performa baik pada kelas dominan. Namun, kelas Tercemar_Berat tidak terdeteksi karena jumlah observasi sangat sedikit, sehingga model kesulitan mempelajari pola untuk kelas minoritas.

3.6 Evaluasi & Perbandingan Akurasi
# Perbandingan akurasi

accs <- c(
SVM  = conf_svm$overall["Accuracy"],
Tree = conf_tree$overall["Accuracy"],
RF   = conf_rf$overall["Accuracy"]
)
cat("\n=== PERBANDINGAN AKURASI ===\n")
## 
## === PERBANDINGAN AKURASI ===
print(accs)
##  SVM.Accuracy Tree.Accuracy   RF.Accuracy 
##     0.8644068     0.9661017     0.9661017
best_model_name <- names(accs)[which.max(accs)]
cat("\nModel terbaik berdasarkan akurasi validasi:", best_model_name, "\n")
## 
## Model terbaik berdasarkan akurasi validasi: Tree.Accuracy
# Interpretasi akurasi

if(best_model_name == "SVM"){ acc <- round(conf_svm$overall["Accuracy"],3)
} else if(best_model_name == "Tree"){ acc <- round(conf_tree$overall["Accuracy"],3)
} else { acc <- round(conf_rf$overall["Accuracy"],3) }

cat("\nInterpretasi Akurasi Model:\n")
## 
## Interpretasi Akurasi Model:
cat(paste0(
"Model klasifikasi terbaik adalah ", best_model_name,
" dengan akurasi sekitar ", acc*100, "%.\n",
"Artinya, model mampu mengklasifikasikan status kualitas air dengan benar pada sekitar ",
round(acc*100,1), "% data validasi.\n",
"Akurasi di atas 80% sudah menunjukkan model cukup baik.\n"
))
## Model klasifikasi terbaik adalah Tree.Accuracy dengan akurasi sekitar 96.6%.
## Artinya, model mampu mengklasifikasikan status kualitas air dengan benar pada sekitar 96.6% data validasi.
## Akurasi di atas 80% sudah menunjukkan model cukup baik.
3.7 Prediksi Status untuk Data Testing (75 baris)

Pada tahap ini, model klasifikasi terbaik yang sudah dilatih (SVM, Decision Tree, Random Forest) digunakan untuk memprediksi Status pada data testing yang belum memiliki label. Model membaca fitur numerik (pH, DO, BOD, TSS, Suhu) dan menghasilkan prediksi untuk masing-masing baris. Hasil prediksi kemudian ditambahkan ke dataframe testing dan disimpan ke file CSV di folder:

D:/Kuliah/Statistika Lingkungan/prediksi_status_testing.csv

# Prediksi status menggunakan ketiga model
pred_svm_test  <- predict(svm_model, test_scaled)
pred_tree_test <- predict(tree_model, test_scaled, type="class")
pred_rf_test   <- predict(rf_model, test_scaled)

# Tambahkan prediksi ke dataframe testing
test_df <- test_df %>%
  mutate(
    Pred_SVM  = pred_svm_test,
    Pred_Tree = pred_tree_test,
    Pred_RF   = pred_rf_test
  )

# Simpan hasil prediksi ke folder spesifik
write.csv(test_df, 
          "D:/Kuliah/Statistika Lingkungan/prediksi_status_testing.csv", 
          row.names = FALSE)

cat("\n✅ Prediksi Status untuk data testing telah disimpan di 'D:/Kuliah/Statistika Lingkungan/prediksi_status_testing.csv'.\n")
## 
## ✅ Prediksi Status untuk data testing telah disimpan di 'D:/Kuliah/Statistika Lingkungan/prediksi_status_testing.csv'.
3.8 Variable Importance (Random Forest)
# 3.8 Variable Importance (Random Forest)
# Hitung importance
rf_importance <- randomForest::importance(rf_model)
print(rf_importance)
##            Baik Tercemar_Ringan Tercemar_Berat MeanDecreaseAccuracy
## pH    0.7833028      -0.7618183    -1.87529165           -0.4983741
## DO   60.4601245      57.0806754     9.22156783           69.3062844
## BOD  59.4636042      53.2867129    11.07591019           64.9630909
## TSS   0.6254071      -1.4311747    -3.52774642           -1.1503440
## Suhu -1.6014858      -2.3524093    -0.02718638           -2.7117089
##      MeanDecreaseGini
## pH           5.528646
## DO          41.579246
## BOD         37.280048
## TSS          6.301945
## Suhu         6.110864
# Visualisasi
library(ggplot2)
importance_df <- data.frame(
  Feature = rownames(rf_importance),
  Importance = rf_importance[,1]  # %IncMSE atau MeanDecreaseGini
)

ggplot(importance_df, aes(x = reorder(Feature, Importance), y = Importance)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Variable Importance (Random Forest)", 
       x = "Fitur", 
       y = "Importance") +
  theme_minimal()

getwd()
## [1] "D:/Kuliah/Statistika Lingkungan"

4. Prediksi Variabel DO (Soal 3)

4.1 Persiapan Data
# Fitur numerik untuk prediksi DO
features_do <- intersect(c("pH","BOD","TSS","Suhu"), names(train_df))
cat("Fitur numerik untuk prediksi DO:", paste(features_do, collapse=", "), "\n")
## Fitur numerik untuk prediksi DO: pH, BOD, TSS, Suhu
# Split train & validation (80/20)
set.seed(123)
trainIndex <- caret::createDataPartition(train_df$DO, p = 0.8, list = FALSE)
train_set <- train_df[trainIndex, ]
val_set   <- train_df[-trainIndex, ]

# Pastikan semua fitur numeric
for(v in features_do){
  train_set[[v]] <- as.numeric(train_set[[v]])
  val_set[[v]]   <- as.numeric(val_set[[v]])
}

Data dibagi untuk melatih model (train) dan mengevaluasi performa (validation). Semua kolom numerik agar model regresi bisa dijalankan.

4.2 Model 1 — Regresi Linear

Regresi linear digunakan untuk memodelkan hubungan linier antara variabel dependen (DO) dengan variabel independen (pH, BOD, TSS, Suhu). Model ini mencoba memperkirakan DO sebagai kombinasi linear dari keempat variabel prediktor.

# Fit Linear Model
lm_model <- lm(DO ~ pH + BOD + TSS + Suhu, data = train_set)

# Print ringkasan model
summary(lm_model)
## 
## Call:
## lm(formula = DO ~ pH + BOD + TSS + Suhu, data = train_set)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.42805 -0.52690  0.01684  0.59689  2.42337 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.6418662  1.2762782   5.204 4.23e-07 ***
## pH          -0.0417877  0.1255561  -0.333    0.740    
## BOD          0.1529793  0.0793183   1.929    0.055 .  
## TSS          0.0009295  0.0066100   0.141    0.888    
## Suhu        -0.0315591  0.0295869  -1.067    0.287    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9597 on 236 degrees of freedom
## Multiple R-squared:   0.02,  Adjusted R-squared:  0.003392 
## F-statistic: 1.204 on 4 and 236 DF,  p-value: 0.3097
# Prediksi pada data validasi
lm_pred <- predict(lm_model, newdata = val_set)
print(head(lm_pred))
##        1        2        3        4        5        6 
## 5.777254 5.747557 5.893798 6.115169 5.719706 5.969021
# Evaluasi performa
lm_mse <- mean((lm_pred - val_set$DO)^2)
lm_rmse <- sqrt(lm_mse)
lm_r2 <- summary(lm_model)$r.squared
cat("\n--- Evaluasi Linear Model ---\n")
## 
## --- Evaluasi Linear Model ---
cat("R² =", round(lm_r2,4), " | MSE =", round(lm_mse,4), " | RMSE =", round(lm_rmse,4), "\n")
## R² = 0.02  | MSE = 0.8406  | RMSE = 0.9168

Berdasarkan hasil, R² hanya 0,02, artinya model hanya mampu menjelaskan sekitar 2% variasi DO. MSE dan RMSE relatif tinggi (0,84 dan 0,92), menunjukkan prediksi kurang akurat. Tidak ada variabel yang signifikan pada tingkat 5%, meskipun BOD mendekati signifikan (p ≈ 0,055). Secara keseluruhan, regresi linear ini kurang efektif untuk memprediksi DO pada dataset ini.

4.3 Model Regresi Spline (GAM)

Regresi spline atau Generalized Additive Model (GAM) digunakan untuk menangkap hubungan non-linear antara DO dengan prediktor (pH, BOD, TSS, Suhu) menggunakan fungsi spline. Model ini lebih fleksibel dibanding regresi linear karena bisa menyesuaikan pola non-linear pada data.

library(mgcv)

# Fit GAM model dengan spline
gam_model <- gam(DO ~ s(pH) + s(BOD) + s(TSS) + s(Suhu), data = train_set, method = "REML")

# Print ringkasan GAM
summary(gam_model)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DO ~ s(pH) + s(BOD) + s(TSS) + s(Suhu)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.97520    0.06105   97.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df     F p-value
## s(pH)   1.000  1.000 0.085   0.770
## s(BOD)  3.637  4.539 1.958   0.101
## s(TSS)  1.000  1.000 0.001   0.981
## s(Suhu) 1.000  1.000 1.564   0.212
## 
## R-sq.(adj) =  0.0281   Deviance explained =  5.5%
## -REML = 338.58  Scale est. = 0.89822   n = 241
# Prediksi pada validasi
gam_pred <- predict(gam_model, newdata = val_set)
print(head(gam_pred))
##        1        2        3        4        5        6 
## 5.777867 5.613550 5.880721 6.128478 5.804982 5.939335
# Evaluasi performa GAM
gam_mse <- mean((gam_pred - val_set$DO)^2)
gam_rmse <- sqrt(gam_mse)
sst <- sum((val_set$DO - mean(train_set$DO))^2)
sse_gam <- sum((val_set$DO - gam_pred)^2)
gam_r2 <- 1 - (sse_gam / sst)
cat("\n--- Evaluasi GAM (Spline) ---\n")
## 
## --- Evaluasi GAM (Spline) ---
cat("R² =", round(gam_r2,4), " | MSE =", round(gam_mse,4), " | RMSE =", round(gam_rmse,4), "\n")
## R² = -0.0622  | MSE = 0.8201  | RMSE = 0.9056

Hasil menunjukkan R² sangat rendah (-0,0622), artinya model GAM tidak mampu menjelaskan variasi DO dengan baik. Nilai MSE dan RMSE (0,82 dan 0,91) masih tinggi, sehingga prediksi kurang akurat. Semua term smooth tidak signifikan secara statistik (p > 0,05), menunjukkan tidak ada prediktor yang memberikan pengaruh nyata pada DO dalam bentuk non-linear pada dataset ini. Secara keseluruhan, GAM tidak lebih baik daripada regresi linear untuk dataset ini.

4.4 Perbandingan Model
eval_df <- data.frame(
  Model = c("Linear Regression", "GAM (Spline)"),
  R2 = c(lm_r2, gam_r2),
  MSE = c(lm_mse, gam_mse),
  RMSE = c(lm_rmse, gam_rmse)
)
print(eval_df)
##               Model          R2       MSE      RMSE
## 1 Linear Regression  0.02000236 0.8406028 0.9168440
## 2      GAM (Spline) -0.06219665 0.8200622 0.9055728

Dari tabel, regresi linear memiliki R² = 0,02, MSE = 0,84, RMSE = 0,92, sedangkan GAM (Spline) memiliki R² negatif (-0,062), MSE = 0,82, RMSE = 0,91.

Interpretasinya, kedua model sama-sama kurang mampu menjelaskan variasi DO, tapi GAM sedikit lebih rendah performanya secara R². Secara praktis, kedua model menunjukkan prediksi DO kurang akurat, dan regresi linear tetap lebih sederhana untuk digunakan.

4.5 Visualisasi Prediksi vs Aktual
library(ggplot2)

# Linear
p1 <- ggplot(val_set, aes(x = DO, y = lm_pred)) +
  geom_point(color = "purple") +
  geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
  labs(title="Linear Regression: Prediksi vs Aktual DO", x="DO Aktual", y="DO Prediksi") +
  theme_minimal()
print(p1)

# GAM
p2 <- ggplot(val_set, aes(x = DO, y = gam_pred)) +
  geom_point(color = "orange") +
  geom_abline(slope=1, intercept=0, color="red", linetype="dashed") +
  labs(title="GAM (Spline): Prediksi vs Aktual DO", x="DO Aktual", y="DO Prediksi") +
  theme_minimal()
print(p2)

4.6 Variabel Paling Berpengaruh
# Linear: koefisien standar
lm_std <- lm(scale(DO) ~ scale(pH) + scale(BOD) + scale(TSS) + scale(Suhu), data = train_set)
coef_std <- summary(lm_std)$coefficients
cat("\n--- Koefisien Standar Linear ---\n")
## 
## --- Koefisien Standar Linear ---
print(coef_std)
##                  Estimate Std. Error       t value   Pr(>|t|)
## (Intercept)  1.080157e-17 0.06430631  1.679706e-16 1.00000000
## scale(pH)   -2.144972e-02 0.06444829 -3.328206e-01 0.73956518
## scale(BOD)   1.245784e-01 0.06459269  1.928677e+00 0.05497044
## scale(TSS)   9.077519e-03 0.06455671  1.406131e-01 0.88829555
## scale(Suhu) -6.883079e-02 0.06452927 -1.066660e+00 0.28721536
# GAM: signifikansi smooth terms
cat("\n--- Signifikansi Variabel pada GAM ---\n")
## 
## --- Signifikansi Variabel pada GAM ---
summary(gam_model)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## DO ~ s(pH) + s(BOD) + s(TSS) + s(Suhu)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.97520    0.06105   97.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df     F p-value
## s(pH)   1.000  1.000 0.085   0.770
## s(BOD)  3.637  4.539 1.958   0.101
## s(TSS)  1.000  1.000 0.001   0.981
## s(Suhu) 1.000  1.000 1.564   0.212
## 
## R-sq.(adj) =  0.0281   Deviance explained =  5.5%
## -REML = 338.58  Scale est. = 0.89822   n = 241

Berdasarkan hasil analisis regresi linear dan GAM, variabel yang paling berpengaruh terhadap kadar DO adalah BOD, meskipun secara statistik belum signifikan pada α = 0,05. Variabel lain seperti pH, TSS, dan Suhu menunjukkan pengaruh yang relatif kecil. Hal ini menunjukkan bahwa peningkatan BOD cenderung berkaitan dengan perubahan DO, namun model saat ini hanya menjelaskan sebagian kecil variasi DO, sehingga kemungkinan terdapat faktor lain yang juga memengaruhi.