1. Paket yang dibutuhkan

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(stringr)
library(ggplot2)
library(readxl)
library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(caret)
## Loading required package: lattice
library(e1071)
## 
## Attaching package: 'e1071'
## The following object is masked from 'package:ggplot2':
## 
##     element
library(rpart)
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:psych':
## 
##     outlier
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:dplyr':
## 
##     combine
library(splines)
library(tidyr)

2. Import Data & Pemeriksaan awal

## Import 
data_all <- read_excel("D:/s5/Komputasi Awan/kalitasair.xlsx")  # ganti path bila perlu

## Tampilkan struktur dan ringkasan awal
str(data_all)
## tibble [300 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Lokasi: chr [1:300] "S1" "S2" "S3" "S4" ...
##  $ pH    : num [1:300] 7.69 6.72 7.18 7.32 7.2 ...
##  $ DO    : num [1:300] NA 5.72 4.89 6.13 7.79 ...
##  $ BOD   : num [1:300] 1.71 1.44 2.73 3.14 1.18 ...
##  $ TSS   : num [1:300] 43.1 44.3 NA 41 48.1 ...
##  $ Suhu  : num [1:300] 26.8 27.7 26 29.7 26.4 ...
##  $ Status: chr [1:300] "Tercemar ringan" "Tercemar ringan" "Tercemar ringan" "Tercemar ringan" ...
summary(data_all)
##     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
head(data_all)
## # A tibble: 6 × 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           
## 6 S6      6.95  8.42  3.23  48.6  28.7 Tercemar ringan

Pada tahap awal saya memuat dataset dan melakukan pemeriksaan struktur dan ringkasan statistik untuk memahami variabel serta mendeteksi masalah awal seperti missing value atau penulisan kategori yang tidak konsisten.

3. Penanganan missing value

colSums(is.na(data_all))
## Lokasi     pH     DO    BOD    TSS   Suhu Status 
##      0      0     23     22     24      0      0
data <- data_all %>%
mutate(across(c(pH, DO, BOD, TSS, Suhu),
~ ifelse(is.na(.), median(., na.rm = TRUE), .)))

menggunakan imputasi median untuk variabel numerik karena median lebih tahan terhadap outlier dibanding mean. Jika missing banyak, pendekatan lain (mis. penghapusan baris) bisa dipertimbangkan dan harus didokumentasikan.

4. Deteksi dan Penangan Outlier

boxplot(data %>% select(pH, DO, BOD, TSS, Suhu),
main = "Boxplot Deteksi Outlier", las = 2)

remove_outliers <- function(df, cols) {
for (col in cols) {
Q1 <- quantile(df[[col]], 0.25, na.rm = TRUE)
Q3 <- quantile(df[[col]], 0.75, na.rm = TRUE)
IQRv <- Q3 - Q1
lower <- Q1 - 1.5 * IQRv
upper <- Q3 + 1.5 * IQRv
df <- df %>% filter(.data[[col]] >= lower & .data[[col]] <= upper)
}
df
}

data_clean <- remove_outliers(data, c("pH","DO","BOD","TSS","Suhu"))

Hasil boxplot menunjukkan bahwa variabel TSS memiliki sebaran data paling lebar dan banyak outlier, menandakan variasi kadar padatan tersuspensi yang tinggi antar sampel. Variabel pH, DO, BOD, dan Suhu relatif stabil namun masih terdapat beberapa nilai ekstrem di luar batas normal. Kehadiran outlier ini menunjukkan adanya perbedaan kondisi lingkungan pada beberapa titik pengamatan yang dapat memengaruhi hasil analisis selanjutnya.

5. Standarisasi Kategori Status

data_clean <- data_clean %>%
mutate(Status = tolower(str_trim(as.character(Status))),
Status = case_when(
Status %in% c("baik", "1", "baik=1", "baik = 1") ~ "1",
Status %in% c("tercemar ringan", "tercemar_ringan", "2", "tercemar ringan=2") ~ "2",
Status %in% c("tercemar berat", "tercemar_berat", "3", "tercemar berat=3") ~ "3",
TRUE ~ NA_character_
)) %>%
mutate(Status = factor(Status, levels = c("1","2","3"),
labels = c("Baik","Tercemar_ringan","Tercemar_berat")))




table(data_clean$Status, useNA = "ifany")
## 
##            Baik Tercemar_ringan  Tercemar_berat 
##              71             207               3

Standarisasi teks ke kategori faktorial memastikan konsistensi label (Baik, Tercemar ringan, Tercemar berat). Baris yang gagal dipetakan menjadi NA perlu diperiksa manual.

6. Statistika Deskriptif

describe(data_clean %>% select(pH, DO, BOD, TSS, Suhu))
##      vars   n  mean   sd median trimmed  mad   min   max range  skew kurtosis
## pH      1 281  6.99 0.48   6.99    6.99 0.49  5.78  8.23  2.45  0.05    -0.39
## DO      2 281  5.96 0.91   5.99    5.98 0.88  3.81  8.24  4.43 -0.16    -0.29
## BOD     3 281  2.99 0.76   3.07    2.99 0.80  0.90  5.10  4.20 -0.08    -0.12
## TSS     4 281 49.76 8.76  49.52   49.78 8.22 27.73 72.41 44.69  0.01    -0.18
## Suhu    5 281 28.12 2.00  28.02   28.12 2.07 22.77 33.25 10.47  0.00    -0.25
##        se
## pH   0.03
## DO   0.05
## BOD  0.05
## TSS  0.52
## Suhu 0.12
prop.table(table(data_clean$Status))
## 
##            Baik Tercemar_ringan  Tercemar_berat 
##      0.25266904      0.73665480      0.01067616

Ringkasan deskriptif (mean, sd, median) memberi gambaran distribusi tiap variabe, proporsi kelas Status penting untuk melihat apakah data imbalance yang mempengaruhi pemilihan metrik evaluasi.

7. Klasifikasi Status Kualitas Air

set.seed(123)
df_model <- data_clean %>%
filter(!is.na(Status)) %>%
select(Status, pH, DO, BOD, TSS, Suhu)

train_idx <- createDataPartition(df_model$Status, p = 0.8, list = FALSE)
train_data <- df_model[train_idx, ]
test_data <- df_model[-train_idx, ]

train_data$Status <- as.factor(train_data$Status)
test_data$Status <- as.factor(test_data$Status)
## SVM
svm_model <- svm(Status ~ pH + DO + BOD + TSS + Suhu,
data = train_data, kernel = "radial", probability = TRUE)
svm_pred <- predict(svm_model, test_data)
svm_cm <- confusionMatrix(svm_pred, test_data$Status)

## Decision Tree 
tree_model <- rpart(Status ~ pH + DO + BOD + TSS + Suhu,
data = train_data, method = "class",
control = rpart.control(cp = 0.01))
tree_pred <- predict(tree_model, test_data, type = "class")
tree_cm <- confusionMatrix(tree_pred, test_data$Status)

## Random Forest 
set.seed(123)
rf_model <- randomForest(Status ~ pH + DO + BOD + TSS + Suhu,
data = train_data, ntree = 500, importance = TRUE)
rf_pred <- predict(rf_model, test_data)
rf_cm <- confusionMatrix(rf_pred, test_data$Status)
akurasi_df <- data.frame(
Model = c("SVM", "Decision Tree", "Random Forest"),
Accuracy = c(svm_cm$overall["Accuracy"], tree_cm$overall["Accuracy"], rf_cm$overall["Accuracy"]),
Kappa = c(svm_cm$overall["Kappa"], tree_cm$overall["Kappa"], rf_cm$overall["Kappa"])
)
print(akurasi_df)
##           Model  Accuracy     Kappa
## 1           SVM 0.8727273 0.6199408
## 2 Decision Tree 0.9636364 0.8994516
## 3 Random Forest 0.9636364 0.8994516

Hasil evaluasi menunjukkan bahwa model Decision Tree dan Random Forest memiliki akurasi tertinggi sebesar 96,36% dengan nilai Kappa 0,899, yang menandakan tingkat kesesuaian prediksi yang sangat baik. Sementara itu, model SVM memiliki akurasi lebih rendah yaitu 87,27% dengan Kappa 0,62, yang menunjukkan performa klasifikasi masih cukup baik namun kurang optimal dibanding dua model lainnya. Secara keseluruhan, Random Forest dan Decision Tree merupakan model terbaik untuk memprediksi status kualitas air pada data ini

8. Prediksi Model DO (Regresi)

df_do <- data_clean %>% select(DO, pH, BOD, TSS, Suhu) %>% na.omit()
set.seed(123)
train_idx2 <- createDataPartition(df_do$DO, p = 0.8, list = FALSE)
train_do <- df_do[train_idx2, ]
test_do <- df_do[-train_idx2, ]

## Linear Regression
lm_model <- lm(DO ~ pH + BOD + TSS + Suhu, data = train_do)
lm_pred <- predict(lm_model, newdata = test_do)
lm_mse <- mean((lm_pred - test_do$DO)^2)
lm_rmse <- sqrt(lm_mse)
lm_r2 <- summary(lm_model)$r.squared

Regresi linear pertama kali dipilih sebagai baseline untuk menilai hubungan linier antar variabel prediktor dan DO.

spline_model <- lm(DO ~
bs(pH, degree = 3, Boundary.knots = range(pH)) +
bs(BOD, degree = 3, Boundary.knots = range(BOD)) +
bs(TSS, degree = 3, Boundary.knots = range(TSS)) +
bs(Suhu, degree = 3, Boundary.knots = range(Suhu)),
data = train_do)

spline_pred <- predict(spline_model, newdata = test_do)
## Warning in bs(TSS, degree = 3L, knots = numeric(0), Boundary.knots = c(27.7275,
## : some 'x' values beyond boundary knots may cause ill-conditioned bases
spline_mse <- mean((spline_pred - test_do$DO)^2)
spline_rmse <- sqrt(spline_mse)
spline_r2 <- summary(spline_model)$r.squared
eval_df <- rbind(
data.frame(Model = "Linear", R2 = lm_r2, MSE = lm_mse, RMSE = lm_rmse),
data.frame(Model = "Spline", R2 = spline_r2, MSE = spline_mse, RMSE = spline_rmse)
)
print(eval_df)
##    Model         R2      MSE     RMSE
## 1 Linear 0.01325826 1.045882 1.022684
## 2 Spline 0.06138635 1.087685 1.042921

Hasil evaluasi model regresi menunjukkan bahwa model Spline memiliki nilai R² sebesar 0,061 yang sedikit lebih tinggi dibandingkan model Linear dengan R² sebesar 0,013, menandakan kemampuan penjelasan variabilitas DO yang sedikit lebih baik. Namun, nilai MSE dan RMSE keduanya relatif besar dan tidak berbeda jauh antar model, sehingga baik model Linear maupun Spline sama-sama belum mampu memprediksi kadar DO dengan akurasi tinggi. Secara keseluruhan, Spline memberikan performa sedikit lebih baik, tetapi peningkatannya masih terbatas.

9. Visualisasi Prediksi dan Aktual

plot_df <- data.frame(
Aktual = test_do$DO,
Pred_LM = lm_pred,
Pred_Spline = spline_pred
)

ggplot(plot_df, aes(x = Aktual, y = Pred_LM)) +
geom_point(alpha = 0.6) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
labs(title = "Prediksi vs Aktual DO (Linear)", x = "Aktual DO", y = "Prediksi DO") +
theme_minimal()

ggplot(plot_df, aes(x = Aktual, y = Pred_Spline)) +
geom_point(alpha = 0.6) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
labs(title = "Prediksi vs Aktual DO (Spline)", x = "Aktual DO", y = "Prediksi DO") +
theme_minimal()

Gambar scatter plot Prediksi vs Aktual DO (Spline) menunjukkan bahwa titik-titik data menyebar cukup jauh dari garis diagonal putus-putus (yang merepresentasikan prediksi sempurna). Pola sebaran ini menandakan bahwa model spline belum mampu memprediksi nilai DO dengan akurasi tinggi, karena hasil prediksi cenderung berpusat pada rentang sempit (sekitar 5,5–6,2), sementara nilai aktual DO memiliki variasi lebih luas. Dengan kata lain, model cenderung underfitting, belum menangkap pola hubungan nonlinier antar variabel secara optimal.Gambar scatter plot Prediksi vs Aktual DO (Linear) menunjukkan Model regresi linear yang digunakan belum mampu menggambarkan hubungan yang kuat antara nilai aktual dan prediksi DO. Hal ini terlihat dari sebaran titik yang cukup menyebar dan tidak mengikuti pola garis diagonal (yang idealnya menunjukkan prediksi mendekati nilai aktual). Garis putus-putus menunjukkan tren prediksi model, namun masih terdapat penyimpangan cukup besar di beberapa titik data. # 10. Variabel Importance untuk DO

set.seed(123)
rf_reg <- randomForest(DO ~ pH + BOD + TSS + Suhu,
data = train_do, ntree = 500, importance = TRUE)
importance(rf_reg)
##        %IncMSE IncNodePurity
## pH    1.929466      37.12633
## BOD  -1.677709      40.95017
## TSS  -1.836042      37.36876
## Suhu -1.860369      37.89069
varImpPlot(rf_reg)

Secara umum, model Random Forest belum memberikan peningkatan akurasi yang signifikan, namun hasil IncNodePurity mengindikasikan bahwa BOD merupakan faktor yang paling berpengaruh terhadap kadar DO di antara variabel lain yang diuji # 11. Prediksi data testing

## Baca data tanpa label Status (misal file test)

test_nolabel <- read_excel("D:/s5/Statistika Lingkungan/training.xlsx")

## Pastikan preprocessing sama seperti data training

test_nolabel <- test_nolabel %>%
mutate(across(c(pH, DO, BOD, TSS, Suhu),
~ ifelse(is.na(.), median(data_clean[[cur_column()]], na.rm = TRUE), .)))

## Prediksi Status

pred_status_test75 <- predict(rf_model, newdata = test_nolabel)

## Gabungkan hasil prediksi ke data asli

out_df <- test_nolabel %>%
mutate(Predicted_Status = pred_status_test75)

## Simpan hasil ke file CSV

write_csv(out_df, "D:/s5/Statistika Lingkungan/prediksi_DO_test75.csv")

## TAMPILKAN 10 BARIS PERTAMA DI OUTPUT HTML

library(knitr)
kable(head(out_df, 10), caption = "Tabel 1. Contoh 10 Baris Pertama Prediksi Status Data 75 Baris")
Tabel 1. Contoh 10 Baris Pertama Prediksi Status Data 75 Baris
Lokasi pH DO BOD TSS Suhu Predicted_Status
S301 6.9977 5.0835 3.1813 50.93390 25.5573 Tercemar_ringan
S302 7.3801 4.7482 3.3373 46.52100 27.0940 Tercemar_ringan
S303 7.0195 6.5949 3.0039 49.52205 26.5954 Baik
S304 7.3675 5.9909 3.4952 39.00910 26.6512 Tercemar_ringan
S305 6.9268 6.2444 3.3449 47.16920 23.3940 Tercemar_ringan
S306 6.9711 6.0028 3.4466 39.10270 27.7013 Tercemar_ringan
S307 7.2412 4.6718 3.3968 45.94330 29.8974 Tercemar_ringan
S308 7.4965 7.1797 4.3295 55.26910 30.3127 Tercemar_ringan
S309 6.3768 5.4072 2.1560 52.40420 32.1487 Tercemar_ringan
S310 6.9833 7.2000 4.2067 58.98360 28.5306 Tercemar_ringan