div style=“text-align:center; font-size:20px; font-weight:bold;”>

Dosen Pembimbing:

Prof. Dr. I Gede Nyoman Mindra Jaya, M.Si.

Disusun Oleh:

Rasendriya Nandana Kurniawan - 140720260011 Yuda Taufiqurahman Wenske - 1407202600

PROGRAM STUDI MAGISTER STATISTIKA TERAPAN

DEPARTEMEN STATISTIKA

FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM

UNIVERSITAS PADJADJARAN

TAHUN AJARAN 2025/2026

Bab 1. Pendahuluan

1.1. Latar Belakang

Energi listrik merupakan kebutuhan dasar yang terus meningkat seiring pertumbuhan penduduk, industrialisasi, dan digitalisasi. Dalam sistem pembangkit listrik, khususnya Combined Cycle Power Plant (CCPP), kemampuan memprediksi daya keluaran listrik (power output) secara akurat menjadi faktor penting untuk meningkatkan efisiensi operasional, menekan biaya bahan bakar, dan menjaga keandalan pasokan. CCPP bekerja dengan menggabungkan turbin gas dan turbin uap, sehingga kinerjanya dipengaruhi oleh berbagai kondisi lingkungan seperti suhu udara, kelembapan, tekanan udara, dan vakum buang.

Dataset yang digunakan dalam penelitian ini memuat lima variabel utama, yaitu AT (suhu lingkungan), V (vakum buang), AP (tekanan udara ambien), RH (kelembapan relatif), dan PE (daya keluaran listrik). Variabel AT, V, AP, dan RH berperan sebagai prediktor, sedangkan PE merupakan variabel respons yang akan diprediksi. Hubungan antara prediktor dan respons tidak selalu linear, sehingga pendekatan regresi linear sederhana berpotensi kurang mampu menangkap pola yang kompleks.

Perkembangan metode statistika dan pembelajaran mesin memberikan alternatif untuk memodelkan hubungan non-linear dan multikolinearitas. Polynomial Regression dapat menangkap efek non-linear melalui suku pangkat, Elastic Net mampu menangani multikolinearitas dan melakukan seleksi variabel melalui regularisasi, sedangkan Generalized Additive Model (GAM) memodelkan hubungan non-linear secara fleksibel melalui fungsi halus (smooth functions). Ketiga metode tersebut memiliki karakteristik yang berbeda, sehingga perlu dibandingkan secara objektif.

Pemilihan model yang tepat memerlukan kriteria yang dapat menyeimbangkan kecocokan model dan kompleksitas. AICc (Corrected Akaike Information Criterion) digunakan untuk pemilihan model pada data berukuran kecil hingga sedang, sedangkan 10-fold Cross-Validation digunakan untuk mengukur kemampuan generalisasi model terhadap data yang tidak digunakan dalam pelatihan. Selain itu, dataset dalam penelitian ini tersedia dalam beberapa sheet, sehingga perlu dianalisis apakah kinerja model berbeda secara signifikan antar sheet.

Berdasarkan uraian tersebut, penelitian ini dilakukan untuk membandingkan Polynomial Regression, Elastic Net, dan GAM dalam memprediksi PE pada data CCPP, memilih model terbaik berdasarkan AICc dan 10-fold Cross-Validation, serta menganalisis perbedaan kinerja model antar sheet menggunakan ANOVA dan Kruskal-Wallis. Hasil penelitian diharapkan dapat memberikan gambaran mengenai metode yang paling sesuai dan memberikan rekomendasi bagi pengelolaan data energi.

1.2. Identifikasi Masalah

Berdasarkan latar belakang, masalah dalam penelitian ini dapat diidentifikasi sebagai berikut.

1. Prediksi PE pada CCPP penting, tetapi hubungan antara AT, V, AP, RH dan PE tidak selalu linear.
2. Terdapat potensi multikolinearitas dan interaksi antar prediktor yang dapat memengaruhi kinerja model regresi.
3. Metode Polynomial Regression, Elastic Net, dan GAM memiliki kelebihan dan keterbatasan yang berbeda, sehingga perlu dibandingkan secara empiris.
4. Pemilihan model terbaik tidak cukup hanya berdasarkan satu kriteria, sehingga diperlukan perbandingan AICc dan 10-fold Cross-Validation.
5. Dataset terdiri atas beberapa sheet, sehingga perlu diketahui apakah kinerja model berbeda nyata antar sheet.
6. Belum banyak penelitian yang secara khusus membandingkan ketiga metode tersebut pada dataset CCPP dengan pendekatan pemilihan model AICc dan cross-validation.

1.3. Tujuan

Penelitian ini dimaksudkan untuk menerapkan dan membandingkan metode Polynomial Regression, Elastic Net, dan Generalized Additive Model (GAM) dalam memprediksi PE berdasarkan variabel AT, V, AP, dan RH pada dataset CCPP yang tersedia dalam beberapa sheet.

Adapun tujuan penelitian ini adalah: 1. Mendeskripsikan karakteristik data AT, V, AP, RH, dan PE pada setiap sheet. 2. Membandingkan kinerja Polynomial Regression, Elastic Net, dan GAM dalam memprediksi PE. 3. Memilih model terbaik berdasarkan AICc dan 10-fold Cross-Validation. 4. Menganalisis perbedaan kinerja model antar sheet menggunakan ANOVA dan Kruskal-Wallis. 5. Memberikan rekomendasi model prediksi PE yang paling sesuai berdasarkan hasil analisis.

1.4. Manfaat Penelitian

Secara praktis, hasil penelitian ini diharapkan dapat memberikan informasi mengenai model prediksi PE yang paling akurat dan stabil. Informasi tersebut dapat dimanfaatkan oleh operator CCPP atau pengelola data energi sebagai dasar pengambilan keputusan, optimasi operasional, dan evaluasi kinerja pembangkit.

Secara akademik, penelitian ini diharapkan memberikan kontribusi pada penerapan metode regresi dan pembelajaran mesin, khususnya perbandingan Polynomial Regression, Elastic Net, dan GAM dengan pemilihan model AICc dan 10-fold Cross-Validation. Penelitian ini juga dapat menjadi referensi untuk studi lanjutan pada dataset energi atau data dengan struktur multi-sheet.

1.5. Batasan Penelitian

Agar ruang lingkup penelitian lebih terarah, ditetapkan batasan sebagai berikut.

1. Data yang digunakan adalah dataset regresi dengan variabel AT, V, AP, RH, dan PE.
2. Dataset terdiri atas beberapa sheet (Sheet1–Sheet4) yang dianalisis secara terpisah.
3. Metode yang dibandingkan adalah Polynomial Regression (derajat 1–3), Elastic Net, dan GAM.
4. Pemilihan model menggunakan AICc dan 10-fold Cross-Validation.
5. Evaluasi model uji menggunakan RMSE dan R².
6. Uji beda kinerja antar sheet menggunakan ANOVA dan Kruskal-Wallis.
7. Penelitian tidak membahas metode deep learning, neural network, atau pemodelan deret waktu.

Bab 2. Tinjauan Pustaka

2.1. Pendahuluan

Penelitian ini menggunakan pendekatan regresi untuk memprediksi PE berdasarkan AT, V, AP, dan RH. Pendekatan yang digunakan mencakup Polynomial Regression, Elastic Net, dan Generalized Additive Model (GAM). Bab ini membahas konsep dasar data mining, machine learning, regresi, regularisasi, GAM, pemilihan model, evaluasi kinerja, serta penelitian terdahulu pada dataset CCPP. Sesuai ketentuan, bab ini disajikan tanpa penulisan rumus.

2.2. Data Mining dan Machine Learning

Data mining adalah proses menggali pola, hubungan, atau pengetahuan baru dari data berukuran besar menggunakan pendekatan komputasional. Machine learning merupakan cabang kecerdasan buatan yang memungkinkan sistem belajar dari data. Dalam konteks prediksi, machine learning digunakan untuk membangun model yang mampu memetakan input ke output. Metode yang digunakan dalam penelitian ini termasuk dalam pendekatan supervised learning karena data memiliki variabel respons PE.

2.3. Regresi Linear dan Polynomial Regression

Regresi linear memodelkan hubungan antara variabel respons dan prediktor secara linear. Namun, hubungan antara kondisi lingkungan dan daya keluaran listrik pada CCPP sering tidak linear. Polynomial Regression memperluas regresi linear dengan menambahkan suku pangkat dan interaksi, sehingga mampu menangkap pola non-linear. Derajat polinomial menentukan fleksibilitas model; derajat terlalu rendah dapat menyebabkan underfitting, sedangkan derajat terlalu tinggi dapat menyebabkan overfitting.

2.4. Elastic Net

Elastic Net merupakan metode regularisasi yang menggabungkan penalti L1 (lasso) dan L2 (ridge). Penalti L1 mendorong seleksi variabel dengan menyusutkan koefisien tertentu menjadi nol, sedangkan penalti L2 menangani multikolinearitas dan menjaga stabilitas estimasi. Elastic Net sangat berguna ketika jumlah prediktor banyak atau terdapat korelasi antar prediktor. Parameter lambda mengatur besarnya penalti, sedangkan alpha mengatur proporsi antara L1 dan L2.

2.5. Generalized Additive Model (GAM)

Generalized Additive Model (GAM) merupakan perluasan model linear yang memodelkan hubungan antara respons dan prediktor melalui fungsi halus (smooth functions). GAM tidak mengasumsikan bentuk hubungan yang kaku, sehingga sangat fleksibel untuk menangkap non-linearitas. Dalam penelitian ini, GAM digunakan dengan fungsi spline dan parameter penghalus (smoothing parameter) yang dipilih melalui GCV atau grid search cross-validation.

2.6. Pemilihan Model: AICc dan Cross-Validation

Pemilihan model perlu mempertimbangkan keseimbangan antara kecocokan dan kompleksitas. AIC (Akaike Information Criterion) mengukur kecocokan model dengan memberikan penalti pada jumlah parameter. AICc merupakan koreksi AIC untuk ukuran sampel kecil. Semakin kecil nilai AICc, semakin baik model.

Cross-Validation (CV) mengukur kemampuan generalisasi model dengan membagi data menjadi beberapa fold. Pada 10-fold CV, data dibagi menjadi 10 bagian; model dilatih pada 9 bagian dan diuji pada 1 bagian secara bergantian. RMSE dari setiap fold dirata-ratakan untuk memperoleh ukuran kinerja yang lebih stabil.

2.7. Evaluasi Kinerja Model

Kinerja model prediksi dievaluasi menggunakan RMSE dan R². RMSE mengukur besarnya kesalahan prediksi dalam satuan variabel respons; semakin kecil RMSE, semakin baik. R² mengukur proporsi keragaman respons yang dapat dijelaskan oleh model; semakin mendekati 1, semakin baik.

2.8. Uji Beda Kinerja antar Sheet

Untuk mengetahui apakah kinerja model berbeda secara signifikan antar sheet, digunakan ANOVA jika asumsi normalitas dan homogenitas varians terpenuhi, atau Kruskal-Wallis sebagai alternatif non-parametrik. Uji ini diterapkan pada nilai CV-fold RMSE dari setiap sheet.

2.9. Penelitian Terdahulu pada Dataset CCPP

Dataset CCPP telah banyak digunakan dalam penelitian machine learning. Tüfekci (2014) membandingkan berbagai metode seperti regresi linear, neural network, dan support vector regression untuk memprediksi daya keluaran. Kaya et al. (2012) menggunakan metode local learning dan global learning. Penelitian-penelitian tersebut menunjukkan bahwa pendekatan non-linear dan regularisasi dapat meningkatkan akurasi prediksi. Namun, perbandingan khusus antara Polynomial Regression, Elastic Net, dan GAM dengan pemilihan model AICc dan 10-fold CV pada data multi-sheet masih terbatas.

2.10. Kerangka Pemikiran

Berdasarkan tinjauan pustaka, penelitian ini berangkat dari kebutuhan prediksi PE yang akurat. Tiga metode regresi dibandingkan pada setiap sheet. Model terbaik dipilih berdasarkan AICc dan 10-fold CV. Kinerja antar sheet dianalisis menggunakan ANOVA dan Kruskal-Wallis. Hasilnya digunakan untuk merekomendasikan model prediksi yang paling sesuai.

Bab 3. Metodologi Penelitian

3.1. Data Penelitian

Penelitian ini menggunakan data sekunder dari Combined Cycle Power Plant (CCPP) Data Set yang tersedia secara publik di UCI Machine Learning Repository. Dataset ini didonasikan pada 25 Maret 2014 dan pertama kali digunakan oleh Tüfekci (2014) untuk memprediksi daya keluaran listrik pada pembangkit siklus gabungan (Combined Cycle Power Plant). Dataset asli memuat 9.568 titik data yang dikumpulkan selama 6 tahun (2006–2011) pada kondisi pembangkit beroperasi dengan beban penuh (full load)

** Tabel 3.1. Variabel penelitian **

\[ \begin{array}{|l|l|l|} \hline \textbf{Variabel} & \textbf{Keterangan} & \textbf{Skala Data} \\ \hline AT & \text{Suhu lingkungan} & \text{Interval} \\ V & \text{Vakum buang} & \text{Interval} \\ AP & \text{Tekanan udara ambien} & \text{Interval} \\ RH & \text{Kelembapan relatif} & \text{Interval} \\ PE & \text{Daya keluaran listrik} & \text{Interval} \\ \hline \end{array} \] Variabel PE berperan sebagai variabel respons, sedangkan AT, V, AP, dan RH berperan sebagai prediktor. Data yang digunakan terdiri atas 5 dataset yang merupakan hasil partisi dari dataset CCPP asli ke dalam beberapa bagian untuk keperluan analisis multi-sheet dalam penelitian ini.

3.2 Metode Analisis

Penelitian ini menggunakan pendekatan kuantitatif untuk memodelkan Net Electrical Output (PE) berdasarkan variabel Ambient Temperature (AT), Exhaust Vacuum (V), Ambient Pressure (AP), dan Relative Humidity (RH). Analisis dilakukan pada setiap sheet menggunakan Polynomial Regression, Elastic Net, dan Generalized Additive Model (GAM). Pemilihan model dilakukan berdasarkan AICc dan 10-fold Cross-Validation, kemudian kinerja model dievaluasi menggunakan data uji.

3.2.1 Eksplorasi Data dan Pemeriksaan Awal

Eksplorasi data dilakukan untuk mengetahui karakteristik dan kualitas data. Analisis meliputi statistik deskriptif, pemeriksaan missing values, data duplikat, outlier, korelasi antarvariabel, dan multikolinearitas. Statistik deskriptif yang digunakan meliputi mean, standar deviasi, minimum, kuartil, median, maksimum, IQR, skewness, dan kurtosis.

Identifikasi outlier dilakukan menggunakan aturan Interquartile Range (IQR). Nilai IQR dihitung sebagai:

\[ IQR=Q_3-Q_1 \]

Batas bawah dan batas atas outlier ditentukan dengan:

\[ Batas\ bawah=Q_1-1.5(IQR) \]

\[ Batas\ atas=Q_3+1.5(IQR) \]

Skewness dan kurtosis digunakan untuk menggambarkan karakteristik distribusi data. Skewness dihitung dengan:

\[ Skewness= \frac{\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^3} {\left(\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2\right)^{3/2}} \]

Kurtosis dihitung dengan:

\[ Kurtosis= \frac{\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^4} {\left(\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2\right)^2}-3 \]

Multikolinearitas antarvariabel prediktor diperiksa menggunakan Variance Inflation Factor (VIF) dengan:

\[ VIF_j=\frac{1}{1-R_j^2} \]

dengan \(R_j^2\) merupakan koefisien determinasi dari regresi prediktor ke-\(j\) terhadap prediktor lainnya.

3.2.2 Pembagian Data dan 10-Fold Cross-Validation

Data pada setiap sheet dibagi menjadi 80% data latih dan 20% data uji secara acak menggunakan set.seed(42). Jika \(n\) merupakan jumlah observasi, jumlah data uji ditentukan dengan:

\[ n_{test}=\lceil0.20n\rceil \]

sedangkan jumlah data latih adalah:

\[ n_{train}=n-n_{test} \]

Data latih digunakan untuk membangun model dan melakukan 10-fold Cross-Validation, sedangkan data uji digunakan untuk mengevaluasi kinerja akhir model. Pada setiap iterasi Cross-Validation, sembilan fold digunakan sebagai data pelatihan dan satu fold sebagai data validasi.

3.2.3 Polynomial Regression

Polynomial Regression digunakan untuk menangkap hubungan nonlinier antara variabel prediktor dan respons. Model polynomial dengan \(p\) prediktor dan derajat \(d\) dinyatakan sebagai:

\[ y_i=\beta_0+ \sum_{j=1}^{p}\sum_{r=1}^{d} \beta_{jr}x_{ij}^{r}+\varepsilon_i \]

Dalam penelitian ini digunakan derajat polynomial 1, 2, dan 3. Estimasi parameter dilakukan menggunakan Ordinary Least Squares (OLS) dengan:

\[ \hat{\beta}=(X^TX)^{-1}X^Ty \]

Residual Sum of Squares (RSS) digunakan sebagai ukuran kesalahan model:

\[ RSS=\sum_{i=1}^{n}(y_i-\hat{y}_i)^2 \]

3.2.4 Elastic Net

Elastic Net digunakan untuk mengatasi kompleksitas model dan potensi multikolinearitas dengan menggabungkan penalti L1 dan L2. Fungsi objektifnya adalah:

\[ \min_{\beta_0,\beta} \left\{ \frac{1}{2n}\sum_{i=1}^{n} (y_i-\beta_0-x_i^T\beta)^2 + \lambda \left[ \frac{1-\alpha}{2}\|\beta\|_2^2+ \alpha\|\beta\|_1 \right] \right\} \]

dengan \(\lambda\) sebagai parameter penalti dan \(\alpha\) sebagai proporsi penalti L1. Nilai \(\alpha=1\) menghasilkan Lasso, sedangkan \(\alpha=0\) menghasilkan Ridge.

Sebelum pemodelan, variabel prediktor distandarisasi berdasarkan data latih dengan:

\[ x_{scaled}=\frac{x-\bar{x}}{s} \]

Parameter \(\lambda\) dan \(\alpha\) ditentukan melalui kombinasi nilai yang diuji pada proses pemodelan.

3.2.5 Generalized Additive Model (GAM)

Generalized Additive Model (GAM) digunakan untuk menangkap hubungan nonlinier melalui fungsi smooth. Bentuk umum model GAM adalah:

\[ y_i=\beta_0+\sum_{j=1}^{p}f_j(x_{ij})+\varepsilon_i \]

Dalam penelitian ini, fungsi smooth digunakan pada masing-masing prediktor dengan basis spline:

\[ y\sim s(AT,k=20)+s(V,k=20)+s(AP,k=20)+s(RH,k=20) \]

Fungsi smooth dibentuk dari kombinasi fungsi basis:

\[ f_j(x)=\sum_{m=1}^{k}\beta_{jm}b_{jm}(x) \]

Kehalusan fungsi dikendalikan menggunakan parameter penghalus \(\lambda\). Estimasi penghalusan juga mempertimbangkan Generalized Cross-Validation (GCV):

\[ GCV= \frac{RSS/n}{(1-edf/n)^2} \]

dengan \(edf\) merupakan effective degrees of freedom.

3.2.6 Pemilihan Model dengan AICc

Pemilihan model berdasarkan Akaike Information Criterion corrected (AICc) dilakukan dengan mempertimbangkan kesesuaian dan kompleksitas model. Nilai AIC dihitung dengan:

\[ AIC=n\log\left(\frac{RSS}{n}\right)+2k \]

Kemudian AICc dihitung sebagai:

\[ AICc=AIC+\frac{2k(k+1)}{n-k-1} \]

dengan \(n\) merupakan jumlah observasi dan \(k\) merupakan jumlah parameter model. Model dengan nilai AICc yang lebih kecil dipilih berdasarkan kriteria AICc.

3.2.7 Evaluasi 10-Fold Cross-Validation

Kinerja model dalam proses Cross-Validation dievaluasi menggunakan Root Mean Squared Error (RMSE). RMSE pada fold ke-\(k\) dihitung dengan:

\[ RMSE_k= \sqrt{ \frac{1}{n_k} \sum_{i\in fold_k} (y_i-\hat{y}_i)^2 } \]

Nilai CV_RMSE diperoleh dari akar rata-rata MSE seluruh fold:

\[ CV\_RMSE= \sqrt{ \frac{1}{K} \sum_{k=1}^{K}MSE_k } \]

dengan \(K=10\). Model dengan nilai CV_RMSE yang lebih kecil menunjukkan kesalahan prediksi yang lebih rendah.

3.2.8 Evaluasi Model pada Data Uji

Model terpilih selanjutnya dievaluasi menggunakan data uji berdasarkan RMSE dan R². RMSE dihitung dengan:

\[ RMSE= \sqrt{ \frac{1}{n_{test}} \sum_{i=1}^{n_{test}} (y_i-\hat{y}_i)^2 } \]

Sedangkan koefisien determinasi dihitung dengan:

\[ R^2= 1- \frac{ \sum_{i=1}^{n_{test}}(y_i-\hat{y}_i)^2 }{ \sum_{i=1}^{n_{test}}(y_i-\bar{y})^2 } \]

RMSE menunjukkan besarnya kesalahan prediksi, sedangkan R² menunjukkan proporsi variasi PE yang dapat dijelaskan oleh model.

3.2.9 Korelasi Peringkat Spearman

Korelasi peringkat Spearman digunakan untuk mengetahui kesesuaian peringkat model berdasarkan AICc dan CV_RMSE. Koefisien Spearman dihitung dengan:

\[ \rho= 1-\frac{6\sum_{i=1}^{n}d_i^2} {n(n^2-1)} \]

dengan \(d_i\) merupakan selisih peringkat antara kedua kriteria dan \(n\) merupakan jumlah model yang dibandingkan.

3.2.10 ANOVA dan Kruskal-Wallis

Perbedaan kinerja model antar-sheet dianalisis menggunakan ANOVA satu arah dan Kruskal-Wallis berdasarkan nilai CV-fold RMSE.

Statistik ANOVA dihitung dengan:

\[ F= \frac{SSB/(k-1)} {SSW/(N-k)} \]

dengan \(SSB\) merupakan jumlah kuadrat antar kelompok dan \(SSW\) merupakan jumlah kuadrat dalam kelompok.

Apabila diperlukan pendekatan nonparametrik, digunakan uji Kruskal-Wallis dengan statistik:

\[ H= \frac{12}{N(N+1)} \sum_{i=1}^{k}\frac{R_i^2}{n_i} -3(N+1) \]

Kedua pengujian digunakan untuk mengetahui apakah terdapat perbedaan kinerja model antar-sheet.

3.2.11 Uji Normalitas dan Heteroskedastisitas

Pemeriksaan asumsi dilakukan menggunakan uji Shapiro-Wilk untuk normalitas residual dan uji Breusch-Pagan untuk heteroskedastisitas.

Statistik Shapiro-Wilk dinyatakan sebagai:

\[ W= \frac{ \left(\sum_{i=1}^{n}a_ix_{(i)}\right)^2 }{ \sum_{i=1}^{n}(x_i-\bar{x})^2 } \]

Sementara itu, statistik Breusch-Pagan dihitung dengan:

\[ LM=nR^2 \]

Hasil kedua pengujian digunakan sebagai informasi pendukung dalam interpretasi model.

3.2.12 Tahapan Penelitian

Tahapan penelitian dilakukan secara sistematis melalui proses

1. membaca dan memvalidasi data,
2. memeriksa kesamaan data antar-sheet,
3. melakukan eksplorasi data,
4. membagi data menjadi data latih dan data uji,
5. melakukan 10-fold Cross-Validation,
6. membangun Polynomial Regression, Elastic Net, dan GAM,
7. melakukan pemilihan model berdasarkan AICc dan CV_RMSE,
8. mengevaluasi model menggunakan RMSE dan R² pada data uji,
9. menganalisis kesesuaian peringkat dan perbedaan kinerja antar-sheet menggunakan korelasi Spearman, ANOVA, dan Kruskal-Wallis.

Bab 4. Hasil Penelitian

cat("\014")

rm(list=ls())

# Tentukan nama file dataset
fname <- "DATASET REGRESI.xlsx"

read_ccpp_file <- function(path) {
  ext <- tolower(tools::file_ext(path))
  if (ext %in% c("xlsx", "xls")) {
    sheets <- readxl::excel_sheets(path)
    data_sheets <- lapply(sheets, function(sh) readxl::read_excel(path, sheet = sh))
    names(data_sheets) <- sheets
  } else if (ext == "csv") {
    data_sheets <- list(Sheet1 = read.csv(path))
  } else {
    stop("Format file harus .xlsx, .xls, atau .csv")
  }
  
  cleaned <- lapply(data_sheets, function(d) {
    colnames(d) <- trimws(as.character(colnames(d)))
    return(as.data.frame(d))
  })
  return(cleaned)
}

target <- "PE"
features <- c("AT", "V", "AP", "RH")

data_sheets <- read_ccpp_file(fname)
sheet_names <- names(data_sheets)

validate_ccpp_data <- function(data_sheets, target, features) {
  if (length(data_sheets) == 0) stop("Tidak ada sheet yang terbaca.")
  need <- c(features, target)
  problems <- c()
  for (sh in names(data_sheets)) {
    d <- data_sheets[[sh]]
    missing <- setdiff(need, colnames(d))
    if (length(missing) > 0) {
      problems <- c(problems, paste0(sh, ": ", paste(missing, collapse = ", ")))
    }
  }
  if (length(problems) > 0) {
    stop(paste("Variabel tidak ditemukan ->", paste(problems, collapse = " | ")))
  }
}

validate_ccpp_data(data_sheets, target, features)

cat(sprintf("Jumlah sheet terbaca: %d -> %s\n", length(sheet_names), paste(sheet_names, collapse = ", ")))

## Jumlah sheet terbaca: 5 -> Sheet1, Sheet2, Sheet3, Sheet4, Sheet5

for (sh in sheet_names) {
  cat(sprintf("  - %s: %d baris, %d kolom\n", sh, nrow(data_sheets[[sh]]), ncol(data_sheets[[sh]])))
}

##   - Sheet1: 9568 baris, 5 kolom
##   - Sheet2: 9568 baris, 5 kolom
##   - Sheet3: 9568 baris, 5 kolom
##   - Sheet4: 9568 baris, 5 kolom
##   - Sheet5: 9568 baris, 5 kolom

# df dipakai untuk EDA di sel-sel berikutnya
df <- data_sheets[[sheet_names[1]]]
cat(sprintf("\nEDA di sel-sel berikutnya memakai sheet pertama ('%s') sebagai representasi.\n", sheet_names[1]))

## 
## EDA di sel-sel berikutnya memakai sheet pertama ('Sheet1') sebagai representasi.

cat(sprintf("(%d, %d)\n", nrow(df), ncol(df)))

## (9568, 5)

head(df)

##      AT     V      AP    RH     PE
## 1 14.96 41.76 1024.07 73.17 463.26
## 2 25.18 62.96 1020.04 59.08 444.37
## 3  5.11 39.40 1012.16 92.14 488.56
## 4 20.86 57.32 1010.24 76.64 446.48
## 5 10.82 37.50 1009.23 96.62 473.90
## 6 26.27 59.44 1012.23 58.77 443.67

verify_identical_sheets <- function(data_sheets, features, target) {
  cols <- c(features, target)
  sheet_names <- names(data_sheets)
  
  ref <- data_sheets[[sheet_names[1]]][, cols]
  ref <- ref[do.call(order, ref), ]
  rownames(ref) <- NULL
  
  rows <- list()
  for (sh in sheet_names) {
    d_sorted <- data_sheets[[sh]][, cols]
    d_sorted <- d_sorted[do.call(order, d_sorted), ]
    rownames(d_sorted) <- NULL
    
    identik <- identical(ref, d_sorted)
    rows[[sh]] <- data.frame(
      Sheet = sh,
      n_baris = nrow(data_sheets[[sh]]),
      identik_dgn_sheet1_setelah_sort = identik
    )
  }
  return(do.call(rbind, rows))
}

verify_df <- verify_identical_sheets(data_sheets, features, target)
verify_df

##         Sheet n_baris identik_dgn_sheet1_setelah_sort
## Sheet1 Sheet1    9568                            TRUE
## Sheet2 Sheet2    9568                            TRUE
## Sheet3 Sheet3    9568                            TRUE
## Sheet4 Sheet4    9568                            TRUE
## Sheet5 Sheet5    9568                            TRUE

str(df)

## 'data.frame':    9568 obs. of  5 variables:
##  $ AT: num  14.96 25.18 5.11 20.86 10.82 ...
##  $ V : num  41.8 63 39.4 57.3 37.5 ...
##  $ AP: num  1024 1020 1012 1010 1009 ...
##  $ RH: num  73.2 59.1 92.1 76.6 96.6 ...
##  $ PE: num  463 444 489 446 474 ...

# Ukuran pusat & sebaran lengkap
calc_mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

calc_skewness <- function(x) {
  n <- length(x)
  m3 <- sum((x - mean(x))^3) / n
  s3 <- (sum((x - mean(x))^2) / n)^(3/2)
  return(m3 / s3)
}

calc_kurtosis <- function(x) {
  n <- length(x)
  m4 <- sum((x - mean(x))^4) / n
  s4 <- (sum((x - mean(x))^2) / n)^2
  return((m4 / s4) - 3)
}

desc_stats <- data.frame(
  Mean = colMeans(df),
  Std = apply(df, 2, sd),
  Min = apply(df, 2, min),
  Q25 = apply(df, 2, quantile, probs = 0.25),
  Median = apply(df, 2, median),
  Mode = apply(df, 2, calc_mode),
  Q75 = apply(df, 2, quantile, probs = 0.75),
  Max = apply(df, 2, max),
  IQR = apply(df, 2, IQR),
  Skewness = apply(df, 2, calc_skewness),
  Kurtosis = apply(df, 2, calc_kurtosis)
)
desc_stats

##          Mean       Std    Min       Q25   Median    Mode     Q75     Max
## AT   19.65123  7.452473   1.81   13.5100   20.345   25.21   25.72   37.11
## V    54.30580 12.707893  25.36   41.7400   52.080   41.17   66.54   81.56
## AP 1013.25908  5.938784 992.89 1009.1000 1012.940 1013.88 1017.26 1033.30
## RH   73.30898 14.600269  25.56   63.3275   74.975  100.09   84.83  100.16
## PE  454.36501 17.066995 420.26  439.7500  451.550  468.80  468.43  495.76
##        IQR   Skewness    Kurtosis
## AT 12.2100 -0.1363717 -1.03763410
## V  24.8000  0.1984899 -1.44420904
## AP  8.1600  0.2654031  0.09356093
## RH 21.5025 -0.4317710 -0.44492113
## PE 28.6800  0.3064614 -1.04860014

cat("Missing values per kolom:\n")

## Missing values per kolom:

print(colSums(is.na(df)))

## AT  V AP RH PE 
##  0  0  0  0  0

cat("\nJumlah baris duplikat:", sum(duplicated(df)), "\n")

## 
## Jumlah baris duplikat: 41

p_list <- lapply(colnames(df), function(col) {
  ggplot(df, aes_string(x = col)) +
    geom_histogram(aes(y = ..density..), fill = "steelblue", bins = 30, alpha = 0.7) +
    geom_density(color = "darkblue", size = 0.8) +
    labs(title = paste("Distribusi", col)) +
    theme_minimal()
})

grid.arrange(grobs = p_list, ncol = 3)

p_box <- lapply(colnames(df), function(col) {
  ggplot(df, aes_string(y = col)) +
    geom_boxplot(fill = "lightcoral") +
    labs(title = col) +
    theme_minimal()
})

grid.arrange(grobs = p_box, ncol = 5)

iqr_outlier_count <- function(s) {
  q1 <- quantile(s, 0.25)
  q3 <- quantile(s, 0.75)
  iqr <- q3 - q1
  lb <- q1 - 1.5 * iqr
  ub <- q3 + 1.5 * iqr
  sum(s < lb | s > ub)
}

cat("Jumlah outlier (aturan IQR 1.5x) per kolom:\n")

## Jumlah outlier (aturan IQR 1.5x) per kolom:

print(sapply(df, iqr_outlier_count))

## AT  V AP RH PE 
##  0  0 88 12  0

cor_matrix <- cor(df)
corrplot(cor_matrix, method = "color", addCoef.col = "black", 
         tl.col = "black", number.digits = 2, cl.lim = c(-1, 1),
         title = "Matriks Korelasi", mar = c(0,0,2,0))

p_reg <- lapply(features, function(col) {
  ggplot(df, aes_string(x = col, y = target)) +
    geom_point(alpha = 0.15, size = 1) +
    geom_smooth(method = "lm", color = "red", se = FALSE) +
    labs(title = paste(col, "vs", target)) +
    theme_minimal()
})

grid.arrange(grobs = p_reg, ncol = 4)

X_check <- lm(as.formula(paste(target, "~", paste(features, collapse = " + "))), data = df)
vif_vals <- car::vif(X_check)
vif_df <- data.frame(
  variabel = names(vif_vals),
  VIF = as.numeric(vif_vals)
)
vif_df

##   variabel      VIF
## 1       AT 5.977602
## 2        V 3.943003
## 3       AP 1.452639
## 4       RH 1.705290

model0 <- lm(as.formula(paste(target, "~", paste(features, collapse = " + "))), data = df)
resid <- residuals(model0)

set.seed(RANDOM_STATE)
sample_resid <- sample(resid, min(5000, length(resid)))
shapiro_res <- shapiro.test(sample_resid)
cat(sprintf("Shapiro-Wilk: statistic=%.4f, p-value=%.4g\n", shapiro_res$statistic, shapiro_res$p.value))

## Shapiro-Wilk: statistic=0.9795, p-value=2.728e-26

p1 <- ggplot(data.frame(resid = resid), aes(x = resid)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "steelblue", alpha = 0.7) +
  geom_density(color = "darkblue") +
  labs(title = "Distribusi Residual") + theme_minimal()

p2 <- ggplot(data.frame(resid = resid), aes(sample = resid)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "QQ-Plot Residual") + theme_minimal()

grid.arrange(p1, p2, ncol = 2)

bp_res <- lmtest::bptest(model0)
cat(sprintf("Breusch-Pagan: statistic=%.4f, p-value=%.4g\n", bp_res$statistic, bp_res$p.value))

## Breusch-Pagan: statistic=33.9302, p-value=7.702e-07

cat("(p-value kecil -> indikasi heteroskedastisitas, perlu dicatat sebagai catatan interpretasi, bukan alasan transformasi otomatis)\n")

## (p-value kecil -> indikasi heteroskedastisitas, perlu dicatat sebagai catatan interpretasi, bukan alasan transformasi otomatis)

Model Selection: AICc vs 10-Fold CV — Semua Sheet (replikasi penuh run_full_analysis())

Bagian ini mereplikasi persis alur analisis: - Polynomial, Elastic Net, dan GAM dijalankan untuk setiap sheet. - Model terbaik (AICc vs CV) dibandingkan per sheet. - Korelasi Spearman AICc vs CV per sheet (rank_corr_df). - ANOVA dan Kruskal-Wallis pada CV-fold RMSE antar sheet.

compute_aic <- function(n, rss, k) {
  n * log(rss / n) + 2 * k
}

compute_aicc <- function(n, rss, k) {
  aic <- compute_aic(n, rss, k)
  if (n - k - 1 > 0) {
    aic <- aic + (2 * k * (k + 1)) / (n - k - 1)
  }
  return(aic)
}

rmse <- function(y_true, y_pred) {
  sqrt(mean((y_true - y_pred)^2))
}

r2_score <- function(y_true, y_pred) {
  1 - sum((y_true - y_pred)^2) / sum((y_true - mean(y_true))^2)
}

make_train_test_split <- function(n, test_size = 0.20, seed = 42) {
  set.seed(seed)
  shuffled <- sample(n)
  n_test <- ceiling(n * test_size)
  test_idx <- shuffled[1:n_test]
  train_idx <- shuffled[(n_test + 1):n]
  return(list(train = train_idx, test = test_idx))
}

make_folds <- function(n, K = 10, seed = 42) {
  set.seed(seed)
  shuffled <- sample(n)
  fold_id <- rep(1:K, length.out = n)
  folds <- lapply(1:K, function(k) sort(shuffled[fold_id == k]))
  return(folds)
}

fit_scaler <- function(X) {
  center <- colMeans(X)
  scale <- apply(X, 2, sd)
  scale[scale == 0 | is.na(scale)] <- 1.0
  return(list(center = center, scale = scale))
}

apply_scaler <- function(X, center, scale) {
  scale(X, center = center, scale = scale)
}

run_polynomial <- function(Xtr, ytr, Xte, yte, folds) {
  n <- length(ytr)
  rows <- list()
  fold_scores <- list()
  # ============================================================
  # Polynomial Regression degree 1-3
  # ============================================================
  for (degree in 1:3) {
    # Membuat polynomial features secara eksplisit
    Xtr_df <- as.data.frame(Xtr)
    Xte_df <- as.data.frame(Xte)
    # Formula polynomial
    poly_terms <- unlist(
      lapply(colnames(Xtr_df), function(v) {
        paste0(
          "poly(", v, ", degree = ", degree,
          ", raw = TRUE)"
        )
      })
    )
    
    formula_poly <- as.formula(
      paste("~", paste(poly_terms, collapse = " + "))
    )
    
    # Design matrix training dan testing
    Xp_tr <- model.matrix(formula_poly, data = Xtr_df)
    Xp_te <- model.matrix(formula_poly, data = Xte_df)
    
    # ------------------------------------------------------------
    # Model pada seluruh training data
    # ------------------------------------------------------------
    lr <- lm(ytr ~ Xp_tr - 1)
    rss <- sum(residuals(lr)^2)
    
    # Jumlah parameter
    k <- ncol(Xp_tr) + 1
    
    aicc <- compute_aicc(
      n = n,
      rss = rss,
      k = k
    )
    
    # ------------------------------------------------------------
    # 10-fold Cross Validation
    # ------------------------------------------------------------
    
    fold_rmse <- numeric(length(folds))
    
    for (i in seq_along(folds)) {
      
      val_idx <- folds[[i]]
      tr_idx <- setdiff(1:n, val_idx)
      
      fold_fit <- lm(
        ytr[tr_idx] ~ Xp_tr[tr_idx, , drop = FALSE] - 1
      )
      
      pred_val <- Xp_tr[val_idx, , drop = FALSE] %*%
        coef(fold_fit)
      
      pred_val <- as.numeric(pred_val)
      
      # Penanganan koefisien NA
      pred_val[is.na(pred_val)] <- mean(ytr[tr_idx])
      
      fold_rmse[i] <- rmse(
        ytr[val_idx],
        pred_val
      )
    }
    
    # Simpan hasil CV setiap degree
    fold_scores[[as.character(degree)]] <- fold_rmse
    
    rows[[degree]] <- data.frame(
      degree = degree,
      k = k,
      AICc = aicc,
      CV_RMSE = sqrt(mean(fold_rmse^2))
    )
  }
  
  # ============================================================
  # Tabel hasil pemilihan degree
  # ============================================================
  
  res <- do.call(rbind, rows)
  
  degree_aic <- res$degree[
    which.min(res$AICc)
  ]
  
  degree_cv <- res$degree[
    which.min(res$CV_RMSE)
  ]
  
  # ============================================================
  # Evaluasi model terpilih
  # ============================================================
  
  evaluate_route <- function(route, degree) {
    
    Xtr_df <- as.data.frame(Xtr)
    Xte_df <- as.data.frame(Xte)
    
    poly_terms <- unlist(
      lapply(colnames(Xtr_df), function(v) {
        paste0(
          "poly(", v, ", degree = ", degree,
          ", raw = TRUE)"
        )
      })
    )
    
    formula_poly <- as.formula(
      paste("~", paste(poly_terms, collapse = " + "))
    )
    
    Xp_tr <- model.matrix(
      formula_poly,
      data = Xtr_df
    )
    
    Xp_te <- model.matrix(
      formula_poly,
      data = Xte_df
    )
    
    fit <- lm(
      ytr ~ Xp_tr - 1
    )
    
    pred_test <- Xp_te %*%
      coef(fit)
    
    pred_test <- as.numeric(pred_test)
    
    pred_test[is.na(pred_test)] <- mean(ytr)
    
    row <- res[
      res$degree == degree,
      ,
      drop = FALSE
    ]
    
    list(
      config = paste0(
        "degree=", degree
      ),
      
      AICc = if (
        route == "AIC"
      ) row$AICc else NA,
      
      CV_RMSE = if (
        route == "CV"
      ) row$CV_RMSE else NA,
      
      test_RMSE = rmse(
        yte,
        pred_test
      ),
      
      test_R2 = r2_score(
        yte,
        pred_test
      )
    )
  }
  # ============================================================
  # Output
  # ============================================================
  list(
    table = res,
    
    final = list(
      AIC = evaluate_route(
        "AIC",
        degree_aic
      ),
      
      CV = evaluate_route(
        "CV",
        degree_cv
      )
    ),
    cv_fold_rmse = fold_scores
  )
}

run_elasticnet <- function(Xtr, ytr, Xte, yte, folds, 
                            lambda_grid = NULL, l1_ratio_grid = c(0.1, 0.5, 0.9, 1.0), seed = 42) {
  if (is.null(lambda_grid)) {
    lambda_grid <- sort(10^seq(-4, 1, length.out = 20), decreasing = TRUE)
  }
  
  scaler <- fit_scaler(Xtr)
  Xtr_s <- apply_scaler(Xtr, scaler$center, scaler$scale)
  Xte_s <- apply_scaler(Xte, scaler$center, scaler$scale)
  n <- nrow(Xtr_s)
  
  # 1. Grid Search AICc
  aic_rows <- list()
  for (l1r in l1_ratio_grid) {
    for (lam in lambda_grid) {
      m <- glmnet(Xtr_s, ytr, alpha = l1r, lambda = lam, standardise = FALSE)
      pred <- predict(m, newx = Xtr_s)
      rss <- sum((ytr - pred)^2)
      k <- sum(as.vector(coef(m)) != 0)
      
      aic_rows[[length(aic_rows) + 1]] <- data.frame(
        l1_ratio = l1r, lambda = lam, k = k, AICc = compute_aicc(n, rss, k)
      )
    }
  }
  en_aic <- do.call(rbind, aic_rows)
  best_aicc <- en_aic[which.min(en_aic$AICc), ]
  
  # 2. Grid Search Cross-Validation
  cv_rows <- list()
  for (l1r in l1_ratio_grid) {
    for (lam in lambda_grid) {
      fold_mse <- numeric(length(folds))
      for (i in seq_along(folds)) {
        val_idx <- folds[[i]]
        tr_idx <- setdiff(1:n, val_idx)
        
        m <- glmnet(Xtr_s[tr_idx, ], ytr[tr_idx], alpha = l1r, lambda = lam, standardise = FALSE)
        pred <- predict(m, newx = Xtr_s[val_idx, ])
        fold_mse[i] <- mean((ytr[val_idx] - pred)^2)
      }
      cv_rows[[length(cv_rows) + 1]] <- data.frame(
        l1_ratio = l1r, lambda = lam, CV_RMSE = sqrt(mean(fold_mse))
      )
    }
  }
  cv_grid <- do.call(rbind, cv_rows)
  best_cv <- cv_grid[which.min(cv_grid$CV_RMSE), ]
  
  # Hitung RMSE per fold untuk model CV terbaik
  cv_fold_rmse <- numeric(length(folds))
  for (i in seq_along(folds)) {
    val_idx <- folds[[i]]
    tr_idx <- setdiff(1:n, val_idx)
    m <- glmnet(Xtr_s[tr_idx, ], ytr[tr_idx], alpha = best_cv$l1_ratio, lambda = best_cv$lambda, standardise = FALSE)
    pred <- predict(m, newx = Xtr_s[val_idx, ])
    cv_fold_rmse[i] <- rmse(ytr[val_idx], pred)
  }
  
  # 3. Evaluasi Test Set
  fit_aic <- glmnet(Xtr_s, ytr, alpha = best_aicc$l1_ratio, lambda = best_aicc$lambda, standardise = FALSE)
  pred_aic <- predict(fit_aic, newx = Xte_s)
  
  fit_cv <- glmnet(Xtr_s, ytr, alpha = best_cv$l1_ratio, lambda = best_cv$lambda, standardise = FALSE)
  pred_cv <- predict(fit_cv, newx = Xte_s)
  
  list(
    table = en_aic,
    cv_grid = cv_grid,
    final = list(
      AIC = list(
        config = sprintf("lambda=%.4f, l1=%.1f", best_aicc$lambda, best_aicc$l1_ratio),
        AICc = best_aicc$AICc, CV_RMSE = NA,
        test_RMSE = rmse(yte, pred_aic), test_R2 = r2_score(yte, pred_aic)
      ),
      CV = list(
        config = sprintf("lambda=%.4f, l1=%.1f", best_cv$lambda, best_cv$l1_ratio),
        AICc = NA, CV_RMSE = mean(cv_fold_rmse),
        test_RMSE = rmse(yte, pred_cv), test_R2 = r2_score(yte, pred_cv)
      )
    ),
    cv_fold_rmse = cv_fold_rmse,
    scaler = scaler
  )
}

run_gam <- function(Xtr, ytr, Xte, yte, folds, 
                    lambda_grid = NULL, lambda_grid_gcv = NULL, k_spline = 20) {
  if (is.null(lambda_grid)) lambda_grid <- 10^seq(-3, 3, length.out = 10)
  
  df_tr <- data.frame(y = ytr, Xtr)
  df_te <- data.frame(y = yte, Xte)
  
  gam_formula <- as.formula(paste("y ~", paste(sprintf("s(%s, k=%d)", colnames(Xtr), k_spline), collapse = " + ")))
  
  gam_gcv <- gam(gam_formula, data = df_tr, method = "GCV.Cp")
  rss <- sum(residuals(gam_gcv)^2)
  edof <- sum(influence(gam_gcv))
  aicc_gcv <- compute_aicc(nrow(df_tr), rss, edof)
  lambda_gcv <- mean(gam_gcv$sp)
  
  cv_rows <- list()
  fold_scores <- list()
  for (lam in lambda_grid) {
    fold_rmse <- numeric(length(folds))
    sp_vec <- rep(lam, length(colnames(Xtr)))
    for (i in seq_along(folds)) {
      val_idx <- folds[[i]]
      tr_idx <- setdiff(1:nrow(df_tr), val_idx)
      
      fit <- gam(gam_formula, data = df_tr[tr_idx, ], sp = sp_vec)
      pred <- predict(fit, newdata = df_tr[val_idx, ])
      fold_rmse[i] <- rmse(df_tr$y[val_idx], pred)
    }
    fold_scores[[as.character(lam)]] <- fold_rmse
    cv_rows[[length(cv_rows) + 1]] <- data.frame(lambda = lam, CV_RMSE = mean(fold_rmse))
  }
  
  cv_df <- do.call(rbind, cv_rows)
  best_cv <- cv_df[which.min(cv_df$CV_RMSE), ]
  best_cv_fold_rmse <- fold_scores[[as.character(best_cv$lambda)]]
  
  pred_test_gcv <- predict(gam_gcv, newdata = df_te)
  
  gam_cv <- gam(gam_formula, data = df_tr, sp = rep(best_cv$lambda, length(colnames(Xtr))))
  pred_test_cv <- predict(gam_cv, newdata = df_te)
  
  list(
    table = cv_df,
    final = list(
      AIC = list(config = sprintf("lambda(GCV)~=%.4g", lambda_gcv), AICc = aicc_gcv, CV_RMSE = NA,
                 test_RMSE = rmse(yte, pred_test_gcv), test_R2 = r2_score(yte, pred_test_gcv)),
      CV = list(config = sprintf("lambda=%.4g", best_cv$lambda), AICc = NA, CV_RMSE = best_cv$CV_RMSE,
                test_RMSE = rmse(yte, pred_test_cv), test_R2 = r2_score(yte, pred_test_cv))
    ),
    cv_fold_rmse = best_cv_fold_rmse,
    gcv_model = gam_gcv,
    cv_model = gam_cv
  )
}

run_full_analysis <- function(data_sheets, target = "PE", features = c("AT", "V", "AP", "RH"), K = 10, seed = 42) {
  sheet_names <- names(data_sheets)
  n <- nrow(data_sheets[[sheet_names[1]]])
  
  split_res <- make_train_test_split(n, test_size = 0.20, seed = seed)
  train_idx <- split_res$train
  test_idx <- split_res$test
  folds <- make_folds(length(train_idx), K = K, seed = seed)
  
  splits <- list()
  for (sh in sheet_names) {
    d <- data_sheets[[sh]]
    splits[[sh]] <- list(
      X_train = d[train_idx, features],
      X_test = d[test_idx, features],
      y_train = d[train_idx, target],
      y_test = d[test_idx, target]
    )
  }
  
  poly_all <- list()
  en_all <- list()
  gam_all <- list()
  
  for (sh in sheet_names) {
    sp <- splits[[sh]]
    cat(sprintf("[%s] menjalankan Polynomial...\n", sh))
    poly_all[[sh]] <- run_polynomial(sp$X_train, sp$y_train, sp$X_test, sp$y_test, folds)
    cat(sprintf("[%s] menjalankan Elastic Net...\n", sh))
    en_all[[sh]] <- run_elasticnet(sp$X_train, sp$y_train, sp$X_test, sp$y_test, folds, seed = seed)
    cat(sprintf("[%s] menjalankan GAM...\n", sh))
    gam_all[[sh]] <- run_gam(sp$X_train, sp$y_train, sp$X_test, sp$y_test, folds)
  }
  
  list(sheet_names = sheet_names, splits = splits, folds = folds,
       poly_all = poly_all, en_all = en_all, gam_all = gam_all)
}

result <- run_full_analysis(data_sheets, target = target, features = features, K = 10, seed = RANDOM_STATE)

## [Sheet1] menjalankan Polynomial...
## [Sheet1] menjalankan Elastic Net...
## [Sheet1] menjalankan GAM...
## [Sheet2] menjalankan Polynomial...
## [Sheet2] menjalankan Elastic Net...
## [Sheet2] menjalankan GAM...
## [Sheet3] menjalankan Polynomial...
## [Sheet3] menjalankan Elastic Net...
## [Sheet3] menjalankan GAM...
## [Sheet4] menjalankan Polynomial...
## [Sheet4] menjalankan Elastic Net...
## [Sheet4] menjalankan GAM...
## [Sheet5] menjalankan Polynomial...
## [Sheet5] menjalankan Elastic Net...
## [Sheet5] menjalankan GAM...

cat("\nSelesai menjalankan Polynomial/Elastic Net/GAM untuk seluruh sheet:", paste(result$sheet_names, collapse = ", "), "\n")

## 
## Selesai menjalankan Polynomial/Elastic Net/GAM untuk seluruh sheet: Sheet1, Sheet2, Sheet3, Sheet4, Sheet5

get_config_degree <- function(config) {
  as.integer(sub(".*=", "", config))
}

build_best_model_df <- function(result) {
  rows <- list()
  for (sh in result$sheet_names) {
    aicc_vals <- c(
      Polynomial = result$poly_all[[sh]]$final$AIC$AICc,
      "Elastic Net" = result$en_all[[sh]]$final$AIC$AICc,
      GAM = result$gam_all[[sh]]$final$AIC$AICc
    )
    cv_vals <- c(
      Polynomial = result$poly_all[[sh]]$final$CV$CV_RMSE,
      "Elastic Net" = result$en_all[[sh]]$final$CV$CV_RMSE,
      GAM = result$gam_all[[sh]]$final$CV$CV_RMSE
    )
    model_best_aic <- names(aicc_vals)[which.min(aicc_vals)]
    model_best_cv <- names(cv_vals)[which.min(cv_vals)]
    
    rows[[sh]] <- data.frame(
      Sheet = sh,
      AICc_Polynomial = aicc_vals["Polynomial"],
      AICc_ElasticNet = aicc_vals["Elastic Net"],
      AICc_GAM = aicc_vals["GAM"],
      CV_RMSE_Polynomial = cv_vals["Polynomial"],
      CV_RMSE_ElasticNet = cv_vals["Elastic Net"],
      CV_RMSE_GAM = cv_vals["GAM"],
      Model_Terbaik_AIC = model_best_aic,
      Model_Terbaik_CV = model_best_cv,
      Sama = (model_best_aic == model_best_cv)
    )
  }
  do.call(rbind, rows)
}

best_model_df <- build_best_model_df(result)
agreement_rate <- mean(best_model_df$Sama)
cat(sprintf("Tingkat kesepakatan model terbaik (AICc vs CV) antar sheet: %.2f%%\n", agreement_rate * 100))

## Tingkat kesepakatan model terbaik (AICc vs CV) antar sheet: 100.00%

best_model_df

##         Sheet AICc_Polynomial AICc_ElasticNet AICc_GAM CV_RMSE_Polynomial
## Sheet1 Sheet1        22252.78        23318.06 21837.24           4.279636
## Sheet2 Sheet2        22241.77        23265.64 21818.80           4.273885
## Sheet3 Sheet3        22249.39        23304.35 21848.26           4.276527
## Sheet4 Sheet4        22069.54        23250.72 21655.34           4.227976
## Sheet5 Sheet5        22203.59        23306.63 21783.32           4.264502
##        CV_RMSE_ElasticNet CV_RMSE_GAM Model_Terbaik_AIC Model_Terbaik_CV Sama
## Sheet1           4.574881    4.151538               GAM              GAM TRUE
## Sheet2           4.565893    4.153729               GAM              GAM TRUE
## Sheet3           4.577456    4.163720               GAM              GAM TRUE
## Sheet4           4.567286    4.117419               GAM              GAM TRUE
## Sheet5           4.577479    4.142578               GAM              GAM TRUE

build_rank_corr_df <- function(best_model_df) {
  rows <- list()
  for (i in 1:nrow(best_model_df)) {
    row <- best_model_df[i, ]
    aicc <- c(row$AICc_Polynomial, row$AICc_ElasticNet, row$AICc_GAM)
    cv <- c(row$CV_RMSE_Polynomial, row$CV_RMSE_ElasticNet, row$CV_RMSE_GAM)
    res <- cor.test(aicc, cv, method = "spearman")
    rows[[i]] <- data.frame(
      Sheet = row$Sheet,
      Spearman_rho = unname(res$estimate),
      p_value = res$p.value
    )
  }
  do.call(rbind, rows)
}

rank_corr_df <- build_rank_corr_df(best_model_df)
rank_corr_df

##    Sheet Spearman_rho   p_value
## 1 Sheet1            1 0.3333333
## 2 Sheet2            1 0.3333333
## 3 Sheet3            1 0.3333333
## 4 Sheet4            1 0.3333333
## 5 Sheet5            1 0.3333333

cv_fold_store <- list(Polynomial = list(), "Elastic Net" = list(), GAM = list())

for (sh in result$sheet_names) {
  deg <- get_config_degree(result$poly_all[[sh]]$final$CV$config)
  cv_fold_store$Polynomial[[sh]] <- result$poly_all[[sh]]$cv_fold_rmse[[as.character(deg)]]
  cv_fold_store$`Elastic Net`[[sh]] <- result$en_all[[sh]]$cv_fold_rmse
  cv_fold_store$GAM[[sh]] <- result$gam_all[[sh]]$cv_fold_rmse
}

anova_rows <- list()
for (model_name in names(cv_fold_store)) {
  per_sheet <- cv_fold_store[[model_name]]
  df_anova <- do.call(rbind, lapply(names(per_sheet), function(sh) {
    data.frame(Sheet = sh, RMSE = per_sheet[[sh]])
  }))
  
  aov_res <- summary(aov(RMSE ~ Sheet, data = df_anova))[[1]]
  f_stat <- aov_res["Sheet", "F value"]
  p_anova <- aov_res["Sheet", "Pr(>F)"]
  
  kw_res <- kruskal.test(RMSE ~ Sheet, data = df_anova)
  h_stat <- unname(kw_res$statistic)
  p_kw <- kw_res$p.value
  
  anova_rows[[model_name]] <- data.frame(
    Model = model_name,
    ANOVA_F = f_stat, ANOVA_p = p_anova,
    KruskalWallis_H = h_stat, KruskalWallis_p = p_kw,
    Signifikan_p_lt_0_05 = (p_anova < 0.05)
  )
}

anova_df <- do.call(rbind, anova_rows)
anova_df

##                   Model     ANOVA_F   ANOVA_p KruskalWallis_H KruskalWallis_p
## Polynomial   Polynomial 0.044540876 0.9961140       1.1143529       0.8919883
## Elastic Net Elastic Net 0.004830853 0.9999516       0.8414118       0.9328116
## GAM                 GAM 0.040660400 0.9967435       1.0202353       0.9067123
##             Signifikan_p_lt_0_05
## Polynomial                 FALSE
## Elastic Net                FALSE
## GAM                        FALSE

plot_data_list <- list()
for (model_name in names(cv_fold_store)) {
  per_sheet <- cv_fold_store[[model_name]]
  for (sh in names(per_sheet)) {
    plot_data_list[[length(plot_data_list) + 1]] <- data.frame(
      Model = model_name, Sheet = sh, CV_Fold_RMSE = per_sheet[[sh]]
    )
  }
}
plot_df <- do.call(rbind, plot_data_list)

ggplot(plot_df, aes(x = Sheet, y = CV_Fold_RMSE)) +
  geom_boxplot(fill = "lightsteelblue") +
  facet_wrap(~ Model, scales = "free_y") +
  labs(title = "Sebaran CV-Fold RMSE per Sheet (dasar uji ANOVA/Kruskal-Wallis)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

detail_rows <- list()
for (sh in result$sheet_names) {
  for (model_name in c("Polynomial", "Elastic Net", "GAM")) {
    all_dict <- switch(model_name,
                       "Polynomial" = result$poly_all,
                       "Elastic Net" = result$en_all,
                       "GAM" = result$gam_all)
    for (route in c("AIC", "CV")) {
      f <- all_dict[[sh]]$final[[route]]
      detail_rows[[length(detail_rows) + 1]] <- data.frame(
        Sheet = sh, Model = model_name, Rute = route,
        Konfigurasi = f$config, AICc = f$AICc, CV_RMSE = f$CV_RMSE,
        Test_RMSE = f$test_RMSE, Test_R2 = f$test_R2
      )
    }
  }
}

detail_df <- do.call(rbind, detail_rows)
detail_df

##     Sheet       Model Rute           Konfigurasi     AICc  CV_RMSE Test_RMSE
## 1  Sheet1  Polynomial  AIC              degree=3 22252.78       NA  4.116409
## 2  Sheet1  Polynomial   CV              degree=3       NA 4.279636  4.116409
## 3  Sheet1 Elastic Net  AIC lambda=0.0038, l1=0.1 23318.06       NA  4.449722
## 4  Sheet1 Elastic Net   CV lambda=0.0038, l1=0.1       NA 4.574881  4.449722
## 5  Sheet1         GAM  AIC  lambda(GCV)~=0.05372 21837.24       NA  3.992165
## 6  Sheet1         GAM   CV       lambda=0.004642       NA 4.151538  3.998908
## 7  Sheet2  Polynomial  AIC              degree=3 22241.77       NA  4.129567
## 8  Sheet2  Polynomial   CV              degree=3       NA 4.273885  4.129567
## 9  Sheet2 Elastic Net  AIC lambda=0.0038, l1=0.1 23265.64       NA  4.512531
## 10 Sheet2 Elastic Net   CV lambda=0.0038, l1=0.1       NA 4.565893  4.512531
## 11 Sheet2         GAM  AIC  lambda(GCV)~=0.05199 21818.80       NA  4.017704
## 12 Sheet2         GAM   CV       lambda=0.004642       NA 4.153729  4.020662
## 13 Sheet3  Polynomial  AIC              degree=3 22249.39       NA  4.120109
## 14 Sheet3  Polynomial   CV              degree=3       NA 4.276527  4.120109
## 15 Sheet3 Elastic Net  AIC lambda=0.0038, l1=0.1 23304.35       NA  4.465787
## 16 Sheet3 Elastic Net   CV lambda=0.0038, l1=0.1       NA 4.577456  4.465787
## 17 Sheet3         GAM  AIC   lambda(GCV)~=0.1037 21848.26       NA  3.983542
## 18 Sheet3         GAM   CV       lambda=0.004642       NA 4.163720  3.983234
## 19 Sheet4  Polynomial  AIC              degree=3 22069.54       NA  4.322818
## 20 Sheet4  Polynomial   CV              degree=3       NA 4.227976  4.322818
## 21 Sheet4 Elastic Net  AIC lambda=0.0038, l1=0.1 23250.72       NA  4.531404
## 22 Sheet4 Elastic Net   CV lambda=0.0038, l1=0.1       NA 4.567286  4.531404
## 23 Sheet4         GAM  AIC    lambda(GCV)~=0.047 21655.34       NA  4.195046
## 24 Sheet4         GAM   CV        lambda=0.02154       NA 4.117419  4.208255
## 25 Sheet5  Polynomial  AIC              degree=3 22203.59       NA  4.171133
## 26 Sheet5  Polynomial   CV              degree=3       NA 4.264502  4.171133
## 27 Sheet5 Elastic Net  AIC lambda=0.0038, l1=0.1 23306.63       NA  4.462542
## 28 Sheet5 Elastic Net   CV lambda=0.0038, l1=0.1       NA 4.577479  4.462542
## 29 Sheet5         GAM  AIC  lambda(GCV)~=0.04472 21783.32       NA  4.050962
## 30 Sheet5         GAM   CV       lambda=0.004642       NA 4.142578  4.059058
##      Test_R2
## 1  0.9414030
## 2  0.9414030
## 3  0.9315294
## 4  0.9315294
## 5  0.9448868
## 6  0.9447005
## 7  0.9433802
## 8  0.9433802
## 9  0.9323917
## 10 0.9323917
## 11 0.9464061
## 12 0.9463272
## 13 0.9419333
## 14 0.9419333
## 15 0.9317809
## 16 0.9317809
## 17 0.9457189
## 18 0.9457273
## 19 0.9340119
## 20 0.9340119
## 21 0.9274901
## 22 0.9274901
## 23 0.9378551
## 24 0.9374632
## 25 0.9404738
## 26 0.9404738
## 27 0.9318659
## 28 0.9318659
## 29 0.9438543
## 30 0.9436296