UTS Statistika Lingkungan Semester 5

————————-

PERSIAPAN

————————-

rm(list = ls())
set.seed(123)          
options(stringsAsFactors = FALSE)

library(readr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(caret)

## Loading required package: lattice

library(e1071)          # svm
library(rpart)          # decision tree
library(rpart.plot)
library(randomForest)

## randomForest 4.7-1.2

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:ggplot2':
## 
##     margin

## The following object is masked from 'package:dplyr':
## 
##     combine

library(splines)        # bs()
library(Metrics)        # mse, rmse

## 
## Attaching package: 'Metrics'

## The following objects are masked from 'package:caret':
## 
##     precision, recall

library(reshape2)

————————-

BACA DATA

————————-

train <- read.csv("C:/Users/ASUS/Downloads/kualitasair.xlsx - Training.csv", na.strings = c("", "NA"))
test  <- read.csv("C:/Users/ASUS/Downloads/kualitasair.xlsx - Testing.csv",  na.strings = c("", "NA"))

# Cek data
cat("Dimensi training:", dim(train), "\n")

## Dimensi training: 300 7

cat("Dimensi testing :", dim(test), "\n")

## Dimensi testing : 75 6

str(train)

## 'data.frame':    300 obs. of  7 variables:
##  $ Lokasi: chr  "S1" "S2" "S3" "S4" ...
##  $ pH    : num  7.69 6.72 7.18 7.32 7.2 ...
##  $ DO    : num  NA 5.72 4.89 6.13 7.79 ...
##  $ BOD   : num  1.71 1.44 2.73 3.14 1.18 ...
##  $ TSS   : num  43.1 44.3 NA 41 48.1 ...
##  $ Suhu  : num  26.8 27.7 26 29.7 26.4 ...
##  $ Status: chr  "Tercemar ringan" "Tercemar ringan" "Tercemar ringan" "Tercemar ringan" ...

summary(train)

##     Lokasi                pH              DO             BOD        
##  Length:300         Min.   :5.503   Min.   :2.982   Min.   :0.3026  
##  Class :character   1st Qu.:6.670   1st Qu.:5.375   1st Qu.:2.3573  
##  Mode  :character   Median :6.988   Median :5.991   Median :3.0661  
##                     Mean   :6.989   Mean   :5.976   Mean   :3.0005  
##                     3rd Qu.:7.318   3rd Qu.:6.688   3rd Qu.:3.5781  
##                     Max.   :8.351   Max.   :9.229   Max.   :5.7962  
##                                     NA's   :23      NA's   :22      
##       TSS             Suhu          Status         
##  Min.   :24.65   Min.   :22.77   Length:300        
##  1st Qu.:43.73   1st Qu.:26.62   Class :character  
##  Median :49.52   Median :28.01   Mode  :character  
##  Mean   :49.70   Mean   :28.31                     
##  3rd Qu.:56.44   3rd Qu.:29.46                     
##  Max.   :76.34   Max.   :90.00                     
##  NA's   :24

————————-

SOAL 1: DATA CLEANING & EKSPLORASI

————————-

1.1 Identifikasi missing value, outlier, inkonsistensi kategori

cat("\n--- Missing values (training) ---\n")

## 
## --- Missing values (training) ---

print(colSums(is.na(train)))

## Lokasi     pH     DO    BOD    TSS   Suhu Status 
##      0      0     23     22     24      0      0

cat("\n--- Missing values (testing) ---\n")

## 
## --- Missing values (testing) ---

print(colSums(is.na(test)))

## Lokasi     pH     DO    BOD    TSS   Suhu 
##      0      0      8      9      7      0

# Jika ada missing: lihat baris contoh
if(any(colSums(is.na(train)) > 0)){
  print("Contoh baris dengan NA (training):")
  print(head(train[!complete.cases(train), ], 5))
}

## [1] "Contoh baris dengan NA (training):"
##    Lokasi     pH     DO    BOD     TSS    Suhu          Status
## 1      S1 7.6855     NA 1.7136 43.1415 26.7972 Tercemar ringan
## 3      S3 7.1816 4.8906 2.7274      NA 26.0255 Tercemar ringan
## 7      S7 7.7558 4.9232     NA 49.0343 29.7409 Tercemar ringan
## 10    S10 6.9686 5.8043     NA      NA 23.7831            BAIK
## 16    S16 7.3180     NA 2.6550 35.8807 28.5968            Baik

1.2 Standarisasi penulisan kategori ‘Status’

if("Status" %in% names(train)){
  # deteksi apakah numeric
  if(is.numeric(train$Status) || all(grepl("^[0-9]+$", as.character(train$Status)))){
    train$Status <- factor(train$Status, levels = c("1", "2", "3"),
                           labels = c("Baik", "Tercemar_ringan", "Tercemar_berat"))
  } else {
    # bersihkan spasi, huruf besar/kecil
    train$Status <- tolower(trimws(train$Status))
    train$Status <- gsub("\\s+", "_", train$Status)
    train$Status <- factor(train$Status)
  }
  cat("Levels Status (training):", levels(train$Status), "\n")
}

## Levels Status (training): baik tercemar_berat tercemar_ringan

memastikan kolom numerik benar tipe numeric

num_vars <- c("pH", "DO", "BOD", "TSS", "Suhu")
for(v in num_vars){
  if(v %in% names(train)){
    train[[v]] <- as.numeric(train[[v]])
  }
  if(v %in% names(test)){
    test[[v]]  <- as.numeric(test[[v]])
  }
}

1.3 Tangani missing value sederhana

na_counts <- colSums(is.na(train))
print(na_counts)

## Lokasi     pH     DO    BOD    TSS   Suhu Status 
##      0      0     23     22     24      0      0

Imputasi Median untuk Variabel Numerik

for(v in num_vars){
  if(na_counts[v] > 0){
    med <- median(train[[v]], na.rm = TRUE)
    train[[v]][is.na(train[[v]])] <- med
    cat("Imputed NA di", v, "dengan median =", med, "\n")
  }
}

## Imputed NA di DO dengan median = 5.9909 
## Imputed NA di BOD dengan median = 3.0661 
## Imputed NA di TSS dengan median = 49.52205

if("Status" %in% names(train)){
  before <- nrow(train)
  train <- train[!is.na(train$Status), ]
  cat("Menghapus", before - nrow(train), "baris karena Status NA\n")
}

## Menghapus 0 baris karena Status NA

1.4 Deteksi dan handling outlier (metode IQR)

detect_outliers <- function(x){
  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
  return(which(x < lower | x > upper))
}

outlier_index <- list()
for(v in num_vars){
  idx <- detect_outliers(train[[v]])
  outlier_index[[v]] <- idx
  cat("Outlier count for", v, ":", length(idx), "\n")
}

## Outlier count for pH : 4 
## Outlier count for DO : 4 
## Outlier count for BOD : 5 
## Outlier count for TSS : 5 
## Outlier count for Suhu : 2

Pilihan penanganan: winsorize (truncate) pada 1% dan 99% quantiles

winsorize <- function(x, lower.p = 0.01, upper.p = 0.99){
  lower <- quantile(x, lower.p, na.rm = TRUE)
  upper <- quantile(x, upper.p, na.rm = TRUE)
  x[x < lower] <- lower
  x[x > upper] <- upper
  return(x)
}

for(v in num_vars){
  train[[v]] <- winsorize(train[[v]], 0.01, 0.99)
  # juga untuk test agar konsisten
  if(v %in% names(test)){
    test[[v]] <- winsorize(test[[v]], 0.01, 0.99)
  }
}

1.5 Ringkasan statistik deskriptif setelah pembersihan

cat("\n--- Ringkasan Statistik Deskriptif (training setelah cleaning) ---\n")

## 
## --- Ringkasan Statistik Deskriptif (training setelah cleaning) ---

print(summary(train[num_vars]))

##        pH              DO             BOD              TSS       
##  Min.   :5.779   Min.   :3.830   Min.   :0.9563   Min.   :27.71  
##  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.989   Mean   :5.977   Mean   :3.0045   Mean   :49.68  
##  3rd Qu.:7.318   3rd Qu.:6.611   3rd Qu.:3.5323   3rd Qu.:55.62  
##  Max.   :8.106   Max.   :8.224   Max.   :5.0996   Max.   :72.23  
##       Suhu      
##  Min.   :23.52  
##  1st Qu.:26.62  
##  Median :28.01  
##  Mean   :28.12  
##  3rd Qu.:29.46  
##  Max.   :33.25

cat("\n--- Distribusi status ---\n")

## 
## --- Distribusi status ---

print(table(train$Status))

## 
##            baik  tercemar_berat tercemar_ringan 
##              72               7             221

Hasil statistik deskriptif menunjukkan bahwa nilai pH air berada pada kisaran normal (5,77–8,10) dengan rata-rata mendekati netral, menandakan kondisi air tidak terlalu asam atau basa. Nilai DO (3,83–8,22 mg/L) menunjukkan ketersediaan oksigen terlarut yang cukup baik, sementara BOD (0,96–5,09 mg/L) mengindikasikan tingkat pencemaran organik masih tergolong moderat. Variabel TSS (27,71–72,23 mg/L) dan Suhu (23,52–33,25°C) berada dalam kisaran wajar untuk perairan sungai tropis. Distribusi kategori Status Kualitas Air menunjukkan mayoritas sampel termasuk “Tercemar ringan” (221 data), sedangkan “Baik” hanya 72 dan “Tercemar berat” 7 data. Hal ini menggambarkan bahwa sebagian besar kondisi air masih tercemar ringan dan terdapat ketidakseimbangan jumlah antar kategori.

Visualisasi singkat (boxplot semua numerik)

train_melt <- melt(train[, c(num_vars)], variable.name = "Variable", value.name = "Value")

## No id variables; using all as measure variables

ggplot(train_melt, aes(x = Variable, y = Value)) +
  geom_boxplot() + ggtitle("Boxplot variabel numerik (setelah cleaning)") + theme_minimal()

Gambar boxplot di atas menampilkan distribusi lima variabel numerik, yaitu pH, DO, BOD, TSS, dan Suhu setelah dilakukan proses pembersihan data (data cleaning). Dari visualisasi tersebut terlihat bahwa sebagian besar variabel memiliki sebaran data yang relatif sempit, kecuali TSS yang menunjukkan variasi nilai cukup tinggi. Nilai median pH, DO, dan BOD tampak berada pada rentang rendah dengan jarak antar kuartil yang kecil, mengindikasikan bahwa data untuk ketiga variabel tersebut cukup homogen. Sementara itu, variabel TSS memiliki rentang yang jauh lebih besar dibandingkan variabel lainnya, menandakan adanya perbedaan yang signifikan dalam skala nilai dan kemungkinan masih terdapat variasi ekstrem dalam data. Variabel Suhu berada pada kisaran menengah dengan sebaran sedang, menunjukkan kondisi yang relatif stabil.

————————-

SOAL 2: KLASIFIKASI STATUS KUALITAS AIR

————————-

2.1 Variabel input: pH, DO, BOD, TSS, Suhu

predictors <- num_vars
response <- "Status"

2.2 Bagi data training internal untuk validasi model (misal 70:30)

set.seed(123)
trainIndex <- createDataPartition(train$Status, p = 0.7, list = FALSE)
train_part <- train[trainIndex, ]
val_part   <- train[-trainIndex, ]

cat("Train part:", nrow(train_part), "Validation part:", nrow(val_part), "\n")

## Train part: 211 Validation part: 89

table(train_part$Status)

## 
##            baik  tercemar_berat tercemar_ringan 
##              51               5             155

table(val_part$Status)

## 
##            baik  tercemar_berat tercemar_ringan 
##              21               2              66

Fungsi untuk menghitung confusion matrix & akurasi

eval_class <- function(true, pred){
  cm <- confusionMatrix(pred, true)
  print(cm)
  return(list(confusion = cm$table, accuracy = cm$overall["Accuracy"]))
}

2.3 Model Decision Tree (rpart)

model_tree <- rpart(Status ~ pH + DO + BOD + TSS + Suhu, data = train_part, method = "class")
pred_tree_val <- predict(model_tree, val_part, type = "class")
cat("\n--- Evaluasi Decision Tree (validasi) ---\n")

## 
## --- Evaluasi Decision Tree (validasi) ---

res_tree <- eval_class(val_part$Status, pred_tree_val)

## Confusion Matrix and Statistics
## 
##                  Reference
## Prediction        baik tercemar_berat tercemar_ringan
##   baik              18              0               1
##   tercemar_berat     0              0               0
##   tercemar_ringan    3              2              65
## 
## Overall Statistics
##                                          
##                Accuracy : 0.9326         
##                  95% CI : (0.859, 0.9749)
##     No Information Rate : 0.7416         
##     P-Value [Acc > NIR] : 3.598e-06      
##                                          
##                   Kappa : 0.816          
##                                          
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: baik Class: tercemar_berat Class: tercemar_ringan
## Sensitivity               0.8571               0.00000                 0.9848
## Specificity               0.9853               1.00000                 0.7826
## Pos Pred Value            0.9474                   NaN                 0.9286
## Neg Pred Value            0.9571               0.97753                 0.9474
## Prevalence                0.2360               0.02247                 0.7416
## Detection Rate            0.2022               0.00000                 0.7303
## Detection Prevalence      0.2135               0.00000                 0.7865
## Balanced Accuracy         0.9212               0.50000                 0.8837

Hasil confusion matrix menunjukkan bahwa model Decision Tree memiliki akurasi sebesar 93,26%, yang berarti model mampu memprediksi status kualitas air dengan tingkat ketepatan yang sangat baik. Nilai Kappa = 0,816 juga menandakan kesesuaian tinggi antara prediksi model dan data aktual. Model paling akurat dalam mengenali kelas “Tercemar ringan” dengan sensitivity 0,9848 dan specificity 0,7826, artinya model sangat baik mendeteksi kelas ini. Untuk kelas “Baik”, sensitivitasnya 0,8571 menunjukkan model cukup baik mengenali kondisi air yang baik. Namun, kelas “Tercemar berat” tidak terdeteksi (sensitivity = 0), kemungkinan karena jumlah datanya sangat sedikit sehingga model kesulitan mengenali pola dari kelas tersebut. Secara keseluruhan, model Decision Tree sudah bekerja efektif, terutama untuk dua kelas dominan (“Baik” dan “Tercemar ringan”).

rpart.plot(model_tree,
           type = 4,
           extra = 101,
           fallen.leaves = TRUE,
           box.palette = "RdYlGn",
           shadow.col = "gray",
           main = "Pohon Keputusan – Kualitas Air")

Gambar di atas menunjukkan hasil visualisasi pohon keputusan (decision tree) yang digunakan untuk mengklasifikasikan kualitas air berdasarkan beberapa parameter seperti DO (Dissolved Oxygen) dan BOD (Biochemical Oxygen Demand). Pohon keputusan ini memisahkan data berdasarkan nilai ambang tertentu untuk menghasilkan kategori kualitas air, yaitu baik dan tercemar ringan. Dari hasil model terlihat bahwa percabangan pertama terjadi pada variabel DO, di mana jika nilai DO lebih kecil dari 6 maka kualitas air cenderung masuk kategori tercemar ringan. Sebaliknya, jika nilai DO ≥ 6, maka percabangan berikutnya ditentukan oleh nilai BOD, di mana BOD < 3.1 diklasifikasikan sebagai baik, sedangkan BOD ≥ 3.1 masih tergolong tercemar ringan. Hasil akhir menunjukkan bahwa sebagian besar sampel air berada pada kategori tercemar ringan (ditandai dengan warna hijau), sementara sebagian kecil tergolong baik (warna kuning).

2.4 Model Random Forest

model_rf <- randomForest(Status ~ pH + DO + BOD + TSS + Suhu, data = train_part, ntree = 500)
pred_rf_val <- predict(model_rf, val_part)
cat("\n--- Evaluasi Random Forest (validasi) ---\n")

## 
## --- Evaluasi Random Forest (validasi) ---

res_rf <- eval_class(val_part$Status, pred_rf_val)

## Confusion Matrix and Statistics
## 
##                  Reference
## Prediction        baik tercemar_berat tercemar_ringan
##   baik              18              0               1
##   tercemar_berat     0              0               0
##   tercemar_ringan    3              2              65
## 
## Overall Statistics
##                                          
##                Accuracy : 0.9326         
##                  95% CI : (0.859, 0.9749)
##     No Information Rate : 0.7416         
##     P-Value [Acc > NIR] : 3.598e-06      
##                                          
##                   Kappa : 0.816          
##                                          
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: baik Class: tercemar_berat Class: tercemar_ringan
## Sensitivity               0.8571               0.00000                 0.9848
## Specificity               0.9853               1.00000                 0.7826
## Pos Pred Value            0.9474                   NaN                 0.9286
## Neg Pred Value            0.9571               0.97753                 0.9474
## Prevalence                0.2360               0.02247                 0.7416
## Detection Rate            0.2022               0.00000                 0.7303
## Detection Prevalence      0.2135               0.00000                 0.7865
## Balanced Accuracy         0.9212               0.50000                 0.8837

Hasil model Random Forest menunjukkan performa yang sama kuat dengan akurasi 93,26% dan nilai Kappa = 0,816, menandakan model memiliki tingkat kesesuaian yang tinggi antara hasil prediksi dan data aktual. Model ini sangat baik dalam mengenali kelas “Tercemar ringan” (sensitivity = 0,9848) serta cukup akurat pada kelas “Baik” (sensitivity = 0,8571). Namun, kelas “Tercemar berat” tidak berhasil terdeteksi (sensitivity = 0) karena datanya sangat sedikit, sehingga model kesulitan mempelajari pola dari kelas tersebut. Nilai specificity yang tinggi di semua kelas (> 0,78) menunjukkan kemampuan model yang baik dalam membedakan antar kategori. Secara keseluruhan, Random Forest memberikan hasil yang stabil dan akurat, terutama untuk dua kelas utama, dan memiliki kemampuan generalisasi yang baik terhadap data validasi.

2.5 Model SVM (e1071)

preproc <- preProcess(train_part[, predictors], method = c("center", "scale"))
train_scaled <- predict(preproc, train_part[, predictors])
val_scaled   <- predict(preproc, val_part[, predictors])

model_svm <- svm(x = train_scaled, y = train_part$Status, probability = FALSE)
pred_svm_val <- predict(model_svm, val_scaled)
cat("\n--- Evaluasi SVM (validasi) ---\n")

## 
## --- Evaluasi SVM (validasi) ---

res_svm <- eval_class(val_part$Status, pred_svm_val)

## Confusion Matrix and Statistics
## 
##                  Reference
## Prediction        baik tercemar_berat tercemar_ringan
##   baik              12              0               1
##   tercemar_berat     0              0               0
##   tercemar_ringan    9              2              65
## 
## Overall Statistics
##                                           
##                Accuracy : 0.8652          
##                  95% CI : (0.7763, 0.9283)
##     No Information Rate : 0.7416          
##     P-Value [Acc > NIR] : 0.003614        
##                                           
##                   Kappa : 0.5942          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: baik Class: tercemar_berat Class: tercemar_ringan
## Sensitivity               0.5714               0.00000                 0.9848
## Specificity               0.9853               1.00000                 0.5217
## Pos Pred Value            0.9231                   NaN                 0.8553
## Neg Pred Value            0.8816               0.97753                 0.9231
## Prevalence                0.2360               0.02247                 0.7416
## Detection Rate            0.1348               0.00000                 0.7303
## Detection Prevalence      0.1461               0.00000                 0.8539
## Balanced Accuracy         0.7784               0.50000                 0.7533

2.6 Bandingkan akurasi dan pilih model terbaik (berdasarkan validasi)

accuracies <- c(
  Decision_Tree = as.numeric(res_tree$accuracy),
  Random_Forest = as.numeric(res_rf$accuracy),
  SVM = as.numeric(res_svm$accuracy)
)
print(accuracies)

## Decision_Tree Random_Forest           SVM 
##     0.9325843     0.9325843     0.8651685

best_model_name <- names(which.max(accuracies))
cat("Model terbaik (berdasarkan akurasi validasi):", best_model_name, "\n")

## Model terbaik (berdasarkan akurasi validasi): Decision_Tree

• Confusion Matrix menunjukkan distribusi prediksi benar dan salah untuk tiap kelas (misalnya “Baik”, “Tercemar ringan”, “Tercemar berat”).

• Akurasi model Decision Tree dan Random Forest sama besar (≈93,26%), menunjukkan keduanya sangat baik dalam mengenali pola data.

• Nilai SVM lebih rendah (≈86,52%), artinya model ini sedikit kurang tepat dalam memisahkan kelas dibanding dua model lainnya.

• Karena akurasi tertinggi sama, kamu memilih Decision Tree sebagai model terbaik (bisa juga mempertimbangkan interpretabilitas dan kesederhanaan model).

2.7 Latih ulang model terbaik pada seluruh data training (train)

if(best_model_name == "Random_Forest"){
  final_clf <- randomForest(Status ~ pH + DO + BOD + TSS + Suhu, data = train, ntree = 500)
} else if(best_model_name == "Decision_Tree"){
  final_clf <- rpart(Status ~ pH + DO + BOD + TSS + Suhu, data = train, method = "class")
} else if(best_model_name == "SVM"){
  # perlu preproc dari seluruh train
  preproc_full <- preProcess(train[, predictors], method = c("center", "scale"))
  train_full_scaled <- predict(preproc_full, train[, predictors])
  final_clf <- svm(x = train_full_scaled, y = train$Status, probability = FALSE)
} else {
  stop("Tidak ada model terbaik yang dipilih.")
}

2.8 Prediksi status pada data testing

if("Status" %in% names(test)){
  test$Status <- NULL    # pastikan tidak ada
}
if(best_model_name == "SVM"){
  # gunakan preproc_full
  test_scaled <- predict(preproc_full, test[, predictors])
  test$Status_pred <- predict(final_clf, test_scaled)
} else {
  test$Status_pred <- predict(final_clf, test[, predictors], type = ifelse(best_model_name=="Decision_Tree","class","response"))
  # for rpart, type="class", for randomForest default predict returns factor.
}

cat("\nPrediksi Status (5 baris pertama pada test):\n")

## 
## Prediksi Status (5 baris pertama pada test):

print(head(test[, c("Lokasi", predictors, "Status_pred")]))

##   Lokasi     pH     DO    BOD     TSS     Suhu     Status_pred
## 1   S301 6.9977 5.0835 3.1813 50.9339 25.55730 tercemar_ringan
## 2   S302 7.3801 4.7482 3.3373 46.5210 27.09400 tercemar_ringan
## 3   S303 7.0195 6.5949 3.0039      NA 26.59540            baik
## 4   S304 7.3675     NA 3.4952 39.0091 26.65120 tercemar_ringan
## 5   S305 6.9268 6.2444 3.3449 47.1692 23.89283 tercemar_ringan
## 6   S306 6.9711 6.0028 3.4466 39.1027 27.70130 tercemar_ringan

————————-

SOAL 3: PREDIKSI VARIABEL DO

————————-

3.1 Bagi data (menggunakan bagian train_part/val_part dari sebelumnya)

3.2 Regresi Linear

formula_lm <- as.formula("DO ~ pH + BOD + TSS + Suhu")
model_lm <- lm(formula_lm, data = train_part)
summary(model_lm)  # koefisien, p-value, R2

## 
## Call:
## lm(formula = formula_lm, data = train_part)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.16202 -0.53333  0.01579  0.64500  2.32973 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.116912   1.378664   4.437 1.49e-05 ***
## pH          -0.007326   0.135300  -0.054    0.957    
## BOD          0.006608   0.082078   0.081    0.936    
## TSS          0.003463   0.007298   0.475    0.636    
## Suhu        -0.010762   0.031553  -0.341    0.733    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9458 on 206 degrees of freedom
## Multiple R-squared:  0.001641,   Adjusted R-squared:  -0.01774 
## F-statistic: 0.08468 on 4 and 206 DF,  p-value: 0.9871

Hasil analisis regresi linier pada gambar menunjukkan bahwa variabel dependen DO (Dissolved Oxygen) dipengaruhi oleh variabel independen pH, BOD, TSS, dan Suhu. Berdasarkan output model, nilai R-squared sebesar 0.001641 dan Adjusted R-squared sebesar -0.01774, yang berarti kemampuan variabel independen dalam menjelaskan variasi DO sangat rendah (hampir tidak ada pengaruh yang signifikan). Nilai p-value model (0.9871) juga jauh lebih besar dari 0.05, menandakan bahwa secara simultan, keempat variabel independen tidak berpengaruh signifikan terhadap DO. Selain itu, dari tabel koefisien terlihat bahwa seluruh variabel (pH, BOD, TSS, dan Suhu) memiliki nilai p-value > 0.05, sehingga tidak ada variabel yang berpengaruh signifikan secara parsial terhadap DO. Dengan demikian, model regresi ini tidak signifikan secara statistik dan tidak mampu menjelaskan hubungan yang berarti antara DO dan variabel-variabel bebasnya.

Prediksi di validation

pred_lm_val <- predict(model_lm, newdata = val_part)
mse_lm_val <- mse(val_part$DO, pred_lm_val)
rmse_lm_val <- rmse(val_part$DO, pred_lm_val)
r2_lm_val <- 1 - sum((val_part$DO - pred_lm_val)^2) / sum((val_part$DO - mean(val_part$DO))^2)

cat("\nLinear Regression (validasi): MSE =", mse_lm_val, "RMSE =", rmse_lm_val, "R2 =", r2_lm_val, "\n")

## 
## Linear Regression (validasi): MSE = 0.889373 RMSE = 0.9430657 R2 = -0.0073021

Hasil evaluasi model regresi linier pada data validasi menunjukkan bahwa nilai MSE = 0.889373, RMSE = 0.9430657, dan R² = -0.0073021. Nilai MSE dan RMSE yang cukup besar menunjukkan bahwa terdapat selisih atau galat yang relatif tinggi antara nilai prediksi dan nilai aktual DO. Selain itu, nilai R² negatif mengindikasikan bahwa model regresi tidak mampu menjelaskan variasi data pada tahap validasi — bahkan performanya lebih buruk dibandingkan model yang hanya menggunakan rata-rata nilai DO sebagai prediksi. # 3.3 Regresi Spline

model_spline <- lm(DO ~ bs(pH, df = 4) + bs(BOD, df = 4) + bs(TSS, df = 4) + bs(Suhu, df = 4),
                   data = train_part)
summary(model_spline)

## 
## Call:
## lm(formula = DO ~ bs(pH, df = 4) + bs(BOD, df = 4) + bs(TSS, 
##     df = 4) + bs(Suhu, df = 4), data = train_part)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.10897 -0.54157  0.02895  0.56943  2.49254 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.8544     0.9935   5.893 1.65e-08 ***
## bs(pH, df = 4)1    -0.4676     0.9002  -0.519    0.604    
## bs(pH, df = 4)2     0.4104     0.6324   0.649    0.517    
## bs(pH, df = 4)3    -0.6140     0.7742  -0.793    0.429    
## bs(pH, df = 4)4     0.1451     0.6630   0.219    0.827    
## bs(BOD, df = 4)1    1.0069     0.7615   1.322    0.188    
## bs(BOD, df = 4)2    0.9297     0.5908   1.574    0.117    
## bs(BOD, df = 4)3    0.3828     0.7220   0.530    0.597    
## bs(BOD, df = 4)4    1.0567     0.7101   1.488    0.138    
## bs(TSS, df = 4)1   -0.9967     0.7736  -1.289    0.199    
## bs(TSS, df = 4)2    0.7447     0.5805   1.283    0.201    
## bs(TSS, df = 4)3   -0.9819     0.7141  -1.375    0.171    
## bs(TSS, df = 4)4    0.3357     0.6340   0.529    0.597    
## bs(Suhu, df = 4)1  -0.2881     0.7130  -0.404    0.687    
## bs(Suhu, df = 4)2  -0.3953     0.5909  -0.669    0.504    
## bs(Suhu, df = 4)3  -0.1988     0.6925  -0.287    0.774    
## bs(Suhu, df = 4)4  -0.3270     0.6027  -0.542    0.588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.951 on 194 degrees of freedom
## Multiple R-squared:  0.04938,    Adjusted R-squared:  -0.02902 
## F-statistic: 0.6299 on 16 and 194 DF,  p-value: 0.8573

Hasil analisis regresi spline pada data pelatihan menunjukkan bahwa model yang menggunakan basis spline untuk variabel pH, BOD, TSS, dan Suhu belum memberikan hasil yang signifikan. Nilai R-squared sebesar 0.04938 dan Adjusted R-squared sebesar -0.02902 menandakan bahwa kemampuan model dalam menjelaskan variasi nilai DO sangat rendah. Nilai p-value F-statistic (0.8573) juga jauh di atas 0.05, sehingga model secara keseluruhan tidak signifikan. Selain itu, semua variabel spline memiliki nilai p-value > 0.05, yang berarti tidak ada pengaruh signifikan secara individual terhadap DO.

Prediksi di validation

pred_spline_val <- predict(model_spline, newdata = val_part)
mse_spline_val <- mse(val_part$DO, pred_spline_val)
rmse_spline_val <- rmse(val_part$DO, pred_spline_val)
r2_spline_val <- 1 - sum((val_part$DO - pred_spline_val)^2) / sum((val_part$DO - mean(val_part$DO))^2)

cat("\nSpline Regression (validasi): MSE =", mse_spline_val, "RMSE =", rmse_spline_val, "R2 =", r2_spline_val, "\n")

## 
## Spline Regression (validasi): MSE = 0.9157713 RMSE = 0.9569594 R2 = -0.03720082

Pada tahap validasi, hasil prediksi menunjukkan MSE = 0.9157713, RMSE = 0.9569594, dan R² = -0.0372. Nilai R² yang negatif kembali menunjukkan bahwa model spline tidak mampu memprediksi DO dengan baik — performanya bahkan lebih buruk daripada model yang hanya menggunakan rata-rata DO sebagai prediksi.

3.4 Pilih model DO terbaik (berdasarkan RMSE atau R2 pada validasi)

if(r2_spline_val > r2_lm_val){
  best_do_model <- "spline"
  final_do_model <- lm(DO ~ bs(pH, df = 4) + bs(BOD, df = 4) + bs(TSS, df = 4) + bs(Suhu, df = 4),
                       data = train)
} else {
  best_do_model <- "linear"
  final_do_model <- lm(DO ~ pH + BOD + TSS + Suhu, data = train)
}
cat("Model DO terbaik berdasarkan validasi:", best_do_model, "\n")

## Model DO terbaik berdasarkan validasi: linear

Berdasarkan hasil perbandingan model pada tahap validasi, model linear dipilih sebagai model terbaik karena memiliki performa yang lebih baik dibandingkan model spline, terutama dilihat dari nilai R² dan RMSE pada data validasi.

Dari hasil sebelumnya:

Model Linear: R² = -0.0073, RMSE = 0.9431

Model Spline: R² = -0.0372, RMSE = 0.9570

Nilai R² pada model linear sedikit lebih tinggi (kurang negatif), dan RMSE-nya juga sedikit lebih kecil dibandingkan model spline. Artinya, meskipun keduanya sama-sama memiliki performa yang kurang baik, model linear menghasilkan prediksi yang lebih mendekati nilai aktual DO dibandingkan model spline.

Dengan kata lain, penambahan kompleksitas melalui spline tidak meningkatkan akurasi model dan justru membuat performanya menurun. Oleh karena itu, secara objektif berdasarkan metrik evaluasi (R² dan RMSE) pada data validasi, model regresi linear sederhana merupakan model DO terbaik di antara kedua model yang diuji.

3.5 Visualisasi: prediksi vs aktual pada validation

plot_df <- data.frame(Aktual = val_part$DO,
                      Pred_LM = pred_lm_val,
                      Pred_Spline = pred_spline_val)

Plot linear

ggplot(plot_df, aes(x = Aktual, y = Pred_LM)) +
  geom_point() + geom_abline(intercept = 0, slope = 1) +
  ggtitle("DO: Prediksi Linear vs Aktual (Validation)") + xlab("Aktual DO") + ylab("Prediksi DO")

Visualisasi hasil prediksi terhadap nilai aktual menunjukkan bahwa baik model regresi linear maupun model regresi spline belum mampu memprediksi nilai DO dengan baik. Pada grafik regresi linear, titik-titik prediksi tersebar secara acak di sekitar garis diagonal, tetapi banyak titik yang jauh dari garis tersebut, menandakan bahwa perbedaan antara nilai aktual dan prediksi cukup besar. Hal ini sejalan dengan nilai R² yang sangat rendah dan RMSE yang tinggi pada tahap validasi.

Plot spline

ggplot(plot_df, aes(x = Aktual, y = Pred_Spline)) +
  geom_point() + geom_abline(intercept = 0, slope = 1) +
  ggtitle("DO: Prediksi Spline vs Aktual (Validation)") + xlab("Aktual DO") + ylab("Prediksi DO")

Grafik regresi spline menunjukkan pola sebaran yang serupa — titik-titik prediksi juga tidak mengikuti garis diagonal dengan baik, bahkan terlihat bahwa hasil prediksi cenderung menyempit pada rentang nilai tertentu, sehingga model kurang mampu menangkap variasi nilai DO aktual. Dengan demikian, kedua model belum memberikan hasil prediksi yang akurat, dan sebaran titik yang acak memperkuat kesimpulan bahwa hubungan antara DO dengan variabel pH, BOD, TSS, dan Suhu tidak dapat dijelaskan secara kuat oleh model yang digunakan.

3.6 Variabel yang paling memengaruhi DO

untuk linear

cat("\n--- Koefisien model linear (train_part) ---\n")

## 
## --- Koefisien model linear (train_part) ---

print(summary(model_lm)$coefficients)

##                 Estimate  Std. Error     t value     Pr(>|t|)
## (Intercept)  6.116911760 1.378664431  4.43683874 1.485821e-05
## pH          -0.007326187 0.135300002 -0.05414772 9.568699e-01
## BOD          0.006608475 0.082078132  0.08051444 9.359063e-01
## TSS          0.003463238 0.007297737  0.47456327 6.356010e-01
## Suhu        -0.010761736 0.031552509 -0.34107386 7.333956e-01

Pada model regresi linear, meskipun variabel seperti BOD (koefisien positif = 0.0066) dan TSS (koefisien positif = 0.0034) memiliki arah hubungan positif terhadap DO, serta pH (-0.0073) dan Suhu (-0.0107) memiliki hubungan negatif, namun tidak satupun yang signifikan secara statistik. Hal ini berarti perubahan nilai pada variabel-variabel tersebut tidak memberikan pengaruh yang berarti terhadap DO dalam model ini. # untuk Spline

cat("\n--- Summary spline (train_part) ---\n")

## 
## --- Summary spline (train_part) ---

print(summary(model_spline))

## 
## Call:
## lm(formula = DO ~ bs(pH, df = 4) + bs(BOD, df = 4) + bs(TSS, 
##     df = 4) + bs(Suhu, df = 4), data = train_part)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.10897 -0.54157  0.02895  0.56943  2.49254 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.8544     0.9935   5.893 1.65e-08 ***
## bs(pH, df = 4)1    -0.4676     0.9002  -0.519    0.604    
## bs(pH, df = 4)2     0.4104     0.6324   0.649    0.517    
## bs(pH, df = 4)3    -0.6140     0.7742  -0.793    0.429    
## bs(pH, df = 4)4     0.1451     0.6630   0.219    0.827    
## bs(BOD, df = 4)1    1.0069     0.7615   1.322    0.188    
## bs(BOD, df = 4)2    0.9297     0.5908   1.574    0.117    
## bs(BOD, df = 4)3    0.3828     0.7220   0.530    0.597    
## bs(BOD, df = 4)4    1.0567     0.7101   1.488    0.138    
## bs(TSS, df = 4)1   -0.9967     0.7736  -1.289    0.199    
## bs(TSS, df = 4)2    0.7447     0.5805   1.283    0.201    
## bs(TSS, df = 4)3   -0.9819     0.7141  -1.375    0.171    
## bs(TSS, df = 4)4    0.3357     0.6340   0.529    0.597    
## bs(Suhu, df = 4)1  -0.2881     0.7130  -0.404    0.687    
## bs(Suhu, df = 4)2  -0.3953     0.5909  -0.669    0.504    
## bs(Suhu, df = 4)3  -0.1988     0.6925  -0.287    0.774    
## bs(Suhu, df = 4)4  -0.3270     0.6027  -0.542    0.588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.951 on 194 degrees of freedom
## Multiple R-squared:  0.04938,    Adjusted R-squared:  -0.02902 
## F-statistic: 0.6299 on 16 and 194 DF,  p-value: 0.8573

pada model regresi spline, hasilnya juga menunjukkan bahwa seluruh komponen spline dari masing-masing variabel memiliki nilai p-value yang besar (semuanya > 0.1). Artinya, meskipun spline mencoba menangkap hubungan non-linier antara DO dan variabel bebas (pH, BOD, TSS, Suhu), model tetap gagal menemukan hubungan yang signifikan

Dengan demikian, dapat disimpulkan bahwa berdasarkan kedua pendekatan (linear dan spline), tidak ada variabel yang secara statistik paling memengaruhi DO dalam dataset ini. Hal ini menunjukkan bahwa faktor-faktor lain di luar pH, BOD, TSS, dan Suhu kemungkinan lebih dominan dalam menentukan kadar DO, atau hubungan antarvariabel bersifat kompleks dan tidak dapat dijelaskan dengan baik oleh model regresi sederhana maupun spline.

3.7 Prediksi DO pada data testing

if("DO" %in% names(test)){
  cat("Jika DO aktual ada di test, tampilkan MSE vs prediksi:\n")
  print(mse(test$DO, test$DO_pred))
}

## Jika DO aktual ada di test, tampilkan MSE vs prediksi:
## [1] NaN

————————-

SIMPAN HASIL

————————-

1) Simpan file csv hasil prediksi (Status_pred dan DO_pred)

output_file <- "hasil_prediksi_test.csv"
write.csv(test, output_file, row.names = FALSE)
cat("Hasil prediksi test disimpan ke:", output_file, "\n")

## Hasil prediksi test disimpan ke: hasil_prediksi_test.csv