3.1 Model dan penaksir

Pada regresi nonlinear, fungsi rata-rata nonlinear terhadap parameter (Fox & Weisberg, 2019):

\[ y_i=f(\mathbf{x}_i;\boldsymbol\theta)+\varepsilon_i,\qquad \boldsymbol\theta\in\mathbb{R}^{k}. \]

Perbedaannya dengan regresi polinomial: polinomial nonlinear terhadap \(\mathbf{x}\) tetapi linear terhadap \(\boldsymbol\beta\), sedangkan pada NLS turunan \(f\) terhadap sebagian parameter masih bergantung pada parameter itu sendiri (contoh: \(\theta_1 e^{\theta_2 z}\)). Penaksir NLS meminimumkan

\[ S(\boldsymbol\theta)=\sum_{i=1}^{n}\big(y_i-f(\mathbf{x}_i;\boldsymbol\theta)\big)^2 . \]

Dengan \(\mathbf{J}(\boldsymbol\theta)\) matriks Jacobian berukuran \(n\times k\) dengan elemen \(J_{ij}=\partial f(\mathbf{x}_i;\boldsymbol\theta)/\partial\theta_j\) dan \(\mathbf{r}(\boldsymbol\theta)=\mathbf{y}-f(\mathbf{X};\boldsymbol\theta)\), persamaan normal NLS adalah \(\mathbf{J}(\hat{\boldsymbol\theta})^\top\mathbf{r}(\hat{\boldsymbol\theta})=\mathbf{0}\). Persamaan ini umumnya tidak mempunyai solusi tertutup sehingga diselesaikan secara iteratif (Fox & Weisberg, 2019; Hansen, 2022). Bila galat \(\varepsilon_i\sim N(0,\sigma^2)\) saling bebas, penaksir NLS identik dengan MLE bagi \(\boldsymbol\theta\) (Hansen, 2022).

3.2 Algoritma iteratif: Gauss–Newton dan Levenberg–Marquardt

Metode Gauss–Newton memperbarui parameter dengan aproksimasi linear \(f(\boldsymbol\theta+\boldsymbol\delta)\approx f(\boldsymbol\theta)+\mathbf{J}\boldsymbol\delta\):

\[ \boldsymbol\theta^{(t+1)}=\boldsymbol\theta^{(t)}+\big(\mathbf{J}^\top\mathbf{J}\big)^{-1}\mathbf{J}^\top\mathbf{r}. \]

Algoritma Levenberg–Marquardt menstabilkan langkah tersebut dengan menambahkan peredam \(\mu\ge 0\) pada diagonal (Levenberg, 1944):

\[ \big(\mathbf{J}^\top\mathbf{J}+\mu\,\mathbf{D}\big)\boldsymbol\delta=\mathbf{J}^\top\mathbf{r},\qquad \boldsymbol\theta^{(t+1)}=\boldsymbol\theta^{(t)}+\boldsymbol\delta , \]

dengan \(\mathbf{D}\) matriks diagonal positif (matriks identitas pada formulasi Levenberg). Nilai \(\mu\) kecil membuat langkah mendekati Gauss–Newton, sedangkan \(\mu\) besar membuat langkah mendekati arah penurunan tercuram. Implementasi yang digunakan pada dokumen ini adalah subrutin MINPACK (Moré et al., 1980) melalui paket minpack.lm (Elzhov et al., 2023), yaitu algoritma yang sama dengan yang dipakai scipy.optimize.curve_fit pada notebook Python.

3.3 Nilai awal dan konvergensi

Algoritma iteratif NLS memerlukan nilai awal parameter, dan hasil dapat bergantung pada nilai awal tersebut (Fox & Weisberg, 2019). Karena itu pada penelitian ini nilai awal ditetapkan secara tetap dan sama untuk setiap fit. Bila algoritma tidak mencapai kriteria konvergensi, fit dinyatakan gagal dan nilai kriterianya dicatat sebagai NA.

5.1 AIC dan AICc

Pada galat normal, log-likelihood maksimum model regresi adalah \(\log L(\hat{\boldsymbol\theta},\hat\sigma^2)=-\tfrac{n}{2}\big[\log(2\pi)+\log(\mathrm{RSS}/n)+1\big]\). Karena \(\mathrm{AIC}=-2\log L+2k\), diperoleh bentuk yang berbeda dari AIC lengkap hanya oleh konstanta yang sama untuk semua model pada data yang sama (Hansen, 2022; Portet, 2020):

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

Untuk sampel berhingga digunakan koreksi (AICc):

\[ \mathrm{AICc}=\mathrm{AIC}+\frac{2k(k+1)}{n-k-1},\qquad n-k-1>0 . \]

Model dengan AICc terkecil dipilih (Portet, 2020). Catatan konvensi: pada implementasi penelitian ini \(k\) adalah jumlah parameter fungsi rata-rata (intercept + koefisien untuk Polynomial; \(\dim\boldsymbol\theta\) untuk NLS; \(\mathrm{edf}\) untuk GAM), tanpa menghitung \(\sigma^2\), mengikuti kode pada notebook Python. Bila \(n-k-1\le 0\), koreksi tidak ditambahkan.

5.2 K-fold Cross-Validation

Pada \(K\)-fold CV, data latih dibagi menjadi \(K\) lipatan yang saling lepas; secara bergantian satu lipatan menjadi data validasi dan sisanya menjadi data pelatihan (James et al., 2021; Arlot & Celisse, 2010). Untuk lipatan ke-\(m\) dengan himpunan indeks validasi \(\mathcal{V}_m\) berukuran \(n_m\):

\[ \mathrm{RMSE}_m=\sqrt{\frac{1}{n_m}\sum_{i\in\mathcal{V}_m}\big(y_i-\hat y_i^{(-m)}\big)^2}, \]

dengan \(\hat y_i^{(-m)}\) prediksi dari model yang dilatih tanpa lipatan ke-\(m\). Ringkasan lintas fold yang digunakan pada penelitian ini ada dua, sesuai notebook:

\[ \mathrm{CV\text{-}RMSE}^{\text{rms}}=\sqrt{\frac{1}{K}\sum_{m=1}^{K}\mathrm{RMSE}_m^{2}} \quad(\text{Polynomial dan NLS}),\qquad \mathrm{CV\text{-}RMSE}^{\text{mean}}=\frac{1}{K}\sum_{m=1}^{K}\mathrm{RMSE}_m \quad(\text{GAM}). \]

Kedua ringkasan sama bila seluruh \(\mathrm{RMSE}_m\) identik, dan keduanya berskala sama dengan RMSE. Model dengan CV-RMSE lebih kecil dinilai memiliki galat prediksi lebih rendah.

2.1 Eksplorasi Data dan Pemeriksaan Awal

Eksplorasi meliputi statistik deskriptif (count, mean, std, min, kuartil, median, modus, maks, IQR, skewness, kurtosis), pemeriksaan missing value dan duplikat, outlier aturan IQR, korelasi antarvariabel, VIF, serta pemeriksaan residual model linear awal (baseline) \(y=\beta_0+\sum_j\beta_jx_j+\varepsilon\) dengan uji Shapiro–Wilk dan Breusch–Pagan. Nilai VIF dihitung seperti pada statsmodels: untuk setiap kolom matriks desain (termasuk kolom konstanta) dilakukan regresi kolom tersebut terhadap kolom lainnya. Hanya VIF keempat prediktor yang dilaporkan dan diinterpretasikan; nilai untuk kolom konstanta (intercept) tidak bermakna dan bergantung pada versi statsmodels.

2.2 Pembagian Data dan 10-Fold Cross-Validation

Pada setiap sheet, data dibagi 80% data latih dan 20% data uji dengan indeks yang sama untuk semua sheet. Prosedurnya adalah permutasi acak penuh indeks \(1,\dots,n\) menggunakan set.seed(42); sebanyak \(n_{\text{test}}\) indeks pertama menjadi data uji dan sisanya data latih, dengan

\[ n_{\text{test}}=\lceil 0{,}20\,n\rceil,\qquad n_{\text{train}}=n-n_{\text{test}} . \]

Pada data latih, 10 fold dibentuk dengan permutasi acak penuh (set.seed(42)) lalu penempelan nomor fold secara round-robin (\(1,2,\dots,10,1,2,\dots\)), bukan blok berurutan. Setiap fold bergantian menjadi data validasi, dan RMSE tiap fold dihitung.

Catatan kesetaraan: generator bilangan acak R dan NumPy berbeda sehingga indeks acak yang persis sama secara numerik tidak dapat dicapai lintas bahasa yang disamakan adalah prosedur (bentuk permutasi, cara membagi, cara membentuk fold, rumus AICc, grid hyperparameter, dan cara agregasi RMSE).

2.3 Polynomial Regression

Model yang dipakai adalah polynomial multivariat berderajat total \(d\in\{1,2,3\}\) dengan suku interaksi, sehingga \(k=5,15,35\). Estimasi memakai OLS dengan solusi kuadrat terkecil berbasis SVD (setara LinearRegression pada scikit-learn yang memusatkan data lalu menyelesaikan kuadrat terkecil), karena kolom polynomial pada skala asli berkondisi buruk. Prediktor tidak distandarisasi pada Polynomial. Untuk setiap \(d\) dihitung \(\mathrm{AICc}\) pada seluruh data latih dan \(\mathrm{CV\text{-}RMSE}^{\text{rms}}\) dari 10 fold. Derajat terpilih rute AICc adalah \(d\) dengan AICc minimum; derajat terpilih rute CV adalah \(d\) dengan CV-RMSE minimum.

2.4 Nonlinear Least Squares (NLS)

Prediktor distandarisasi, \(z_{ij}=(x_{ij}-\bar x_j)/s_j\) dengan \(\bar x_j\) dan \(s_j\) (simpangan baku sampel, pembagi \(n-1\)) dihitung hanya dari data yang dipakai pada fit tersebut sehingga tidak ada kebocoran data (bila \(s_j=0\) atau tak terdefinisi digunakan 1). Standardisasi ini dimaksudkan agar skala prediktor setara bagi optimisasi numerik. Tiga bentuk fungsi nonlinear terhadap parameter dibandingkan, dengan indeks \(z_1,\dots,z_4\) berturut-turut untuk AT, V, AP, RH:

\[ \begin{aligned} \text{exp(AT)}:\quad & f=\theta_0+\theta_1e^{\theta_2 z_{1}}+\theta_3z_{2}+\theta_4z_{3}+\theta_5z_{4} && (k=6)\\ \text{exp(AT)+exp(V)}:\quad & f=\theta_0+\theta_1e^{\theta_2 z_{1}}+\theta_3e^{\theta_4 z_{2}}+\theta_5z_{3}+\theta_6z_{4} && (k=7)\\ \text{exp(linear)}:\quad & f=\exp\!\big(\theta_0+\theta_1z_{1}+\theta_2z_{2}+\theta_3z_{3}+\theta_4z_{4}\big) && (k=5) \end{aligned} \]

Bentuk fungsi eksponensial dipilih karena eksplorasi data menunjukkan hubungan AT dan V terhadap PE menurun dengan kemiringan yang melandai (pola peluruhan eksponensial), seperti dirancang pada notebook. Nilai awal ditetapkan tetap dan sama pada semua fit: koefisien linear \(b_j\) dari OLS pada \(\mathbf{Z}\), parameter pengali eksponen \(=1\), dan eksponen \(=-0{,}5\):

\[ \begin{aligned} \text{exp(AT)}:\quad & \boldsymbol\theta^{(0)}=(\bar y,\;1,\;-0{,}5,\;b_{2},\;b_{3},\;b_{4})\\ \text{exp(AT)+exp(V)}:\quad & \boldsymbol\theta^{(0)}=(\bar y,\;1,\;-0{,}5,\;1,\;-0{,}5,\;b_{3},\;b_{4})\\ \text{exp(linear)}:\quad & \boldsymbol\theta^{(0)}=(\log\bar y,\;0,\;0,\;0,\;0) \end{aligned} \]

Estimasi dengan Levenberg–Marquardt (batas evaluasi fungsi 20.000, toleransi relatif \(\sqrt{\epsilon_{\text{mesin}}}\approx1{,}49\times10^{-8}\)). Fit yang tidak mencapai kriteria konvergensi dicatat sebagai gagal (NA) dan tidak dipilih. Bentuk fungsi berperan sebagai hyperparameter yang dipilih lewat rute AICc dan rute CV, seperti derajat pada Polynomial; \(\mathrm{CV\text{-}RMSE}^{\text{rms}}\) dipakai seperti pada Polynomial.

2.5 Generalized Additive Model (GAM)

Model yang dipakai:

\[ y\sim s(\text{AT},q{=}20)+s(\text{V},q{=}20)+s(\text{AP},q{=}20)+s(\text{RH},q{=}20), \]

dengan basis P-spline. Dua rute pemilihan \(\lambda\):

• Rute AICc/GCV: \(\lambda\) tiap smooth dipilih otomatis dengan GCV (method = "GCV.Cp"). Nilai \(\mathrm{edf}\) diambil dari jumlah elemen diagonal matriks pengaruh, dan \(\mathrm{AICc}\) dihitung dengan \(k=\mathrm{edf}\). Konfigurasi dilaporkan sebagai rata-rata \(\lambda\) terpilih.
• Rute CV: \(\lambda\) yang sama diberlakukan pada keempat smooth dan dipilih dari grid \(\lambda\in 10^{\,\text{seq}(-3,3,\text{length}=10)}\); untuk tiap \(\lambda\) dihitung \(\mathrm{CV\text{-}RMSE}^{\text{mean}}\) dari 10 fold, lalu \(\lambda\) dengan nilai minimum dipilih.

Catatan kesetaraan: skala dan parameterisasi \(\lambda\) pada mgcv dan pyGAM tidak identik secara numerik; yang disamakan adalah desain model (4 smooth, \(q=20\), basis P-spline), desain grid, dan prosedur pemilihan.

2.6 Pemilihan Model per Sheet dan Model Akhir

Untuk tiap sheet, “model terbaik rute AICc” adalah model dengan AICc minimum di antara Polynomial (derajat terpilih AICc), NLS (bentuk terpilih AICc), dan GAM (rute GCV); “model terbaik rute CV” adalah model dengan CV-RMSE minimum di antara ketiganya (konfigurasi terpilih rute CV). Kesesuaian kedua rute dicatat pada kolom Sama.

Model akhir per sheet ditentukan dari tiga kandidat (satu per metode) dengan AICc terkecil tiap metode, CV-RMSE pada konfigurasi terpilih rute CV, serta metrik data uji dari rute CV. Kandidat bertanda AICc terbaik atau CV terbaik bila nilainya sama dengan minimum antar-kandidat (isclose, toleransi relatif \(10^{-5}\) dan mutlak \(10^{-8}\)). Jika terdapat kandidat yang terbaik pada kedua kriteria, dipilih yang ber-AICc terkecil di antaranya; jika tidak, dipilih kandidat ber-CV-RMSE terkecil.

2.7 Evaluasi pada Data Uji, Spearman, ANOVA, dan Kruskal–Wallis

Model terpilih dilatih ulang pada seluruh data latih dan dievaluasi pada data uji dengan RMSE dan \(R^2\). Kesesuaian peringkat AICc dan CV-RMSE antar tiga model diuji dengan korelasi Spearman per sheet (\(M=3\)). Perbedaan kinerja antar-sheet diuji dengan ANOVA satu arah dan Kruskal–Wallis pada 10 nilai RMSE fold dari konfigurasi terpilih rute CV, untuk masing-masing metode (\(G\) = jumlah sheet).

2.8 Tahapan Penelitian

1. membaca dan memvalidasi data seluruh sheet,
2. memeriksa kesamaan data antar-sheet,
3. melakukan eksplorasi data,
4. membagi data menjadi data latih dan data uji,
5. membentuk 10 fold pada data latih,
6. membangun Polynomial Regression, NLS, dan GAM pada setiap sheet,
7. memilih model berdasarkan AICc dan CV-RMSE,
8. mengevaluasi model pada data uji (RMSE dan \(R^2\)),
9. menganalisis kesesuaian peringkat (Spearman) dan perbedaan kinerja antar-sheet (ANOVA dan Kruskal–Wallis).

1.1 Struktur dan Karakteristik Data

Pemeriksaan dilakukan secara berurutan: pembacaan dan validasi seluruh sheet, pemeriksaan kesamaan antar-sheet, dan statistik deskriptif. Ubah fname sesuai lokasi berkas data pada komputer Anda.

fname <- "D:/FASTTRACK S2/REGRESI/PROJECT/Folds5x2_pp.xlsx"

read_ccpp_file <- function(path) {
  sheets <- readxl::excel_sheets(path)
  out <- lapply(sheets, function(sh) {
    d <- readxl::read_excel(path, sheet = sh)
    colnames(d) <- trimws(as.character(colnames(d)))
    as.data.frame(d)
  })
  names(out) <- sheets
  out
}

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)) {
    missing <- setdiff(need, colnames(data_sheets[[sh]]))
    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 x %d kolom\n", sh, nrow(data_sheets[[sh]]), ncol(data_sheets[[sh]])))

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

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

## 
## EDA memakai sheet pertama ('Sheet1') sebagai representasi.

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)
  sort_df <- function(d) { d <- d[, cols]; d <- d[do.call(order, unname(as.list(d))), ]; rownames(d) <- NULL; d }
  ref <- sort_df(data_sheets[[sheet_names[1]]])
  rows <- lapply(sheet_names, function(sh) {
    data.frame(Sheet = sh,
               n_baris = nrow(data_sheets[[sh]]),
               identik_dgn_sheet1_setelah_sort = isTRUE(all.equal(ref, sort_df(data_sheets[[sh]]),
                                                                  check.attributes = FALSE)))
  })
  do.call(rbind, rows)
}

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

Sheet	n_baris	identik_dgn_sheet1_setelah_sort
Sheet1	9568	TRUE
Sheet2	9568	TRUE
Sheet3	9568	TRUE
Sheet4	9568	TRUE
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 (definisi sama dengan pandas pada notebook)
calc_mode <- function(x) {                # modus; bila kembar, nilai terkecil (seperti pandas .mode().iloc[0])
  ux <- sort(unique(x))
  ux[which.max(tabulate(match(x, ux)))]
}
calc_skewness <- function(x) {            # G1 (terkoreksi bias), setara pandas .skew()
  n <- length(x); m <- mean(x); m2 <- mean((x - m)^2); m3 <- mean((x - m)^3)
  sqrt(n * (n - 1)) / (n - 2) * m3 / m2^1.5
}
calc_kurtosis <- function(x) {            # G2 (kurtosis berlebih terkoreksi bias), setara pandas .kurt()
  n <- length(x); m <- mean(x); m2 <- mean((x - m)^2); m4 <- mean((x - m)^4)
  g2 <- m4 / m2^2 - 3
  (n - 1) / ((n - 2) * (n - 3)) * ((n + 1) * g2 + 6)
}

desc_stats <- data.frame(
  Count    = sapply(df, function(x) sum(!is.na(x))),
  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),
  Q75      = apply(df, 2, quantile, probs = 0.75),
  Max      = apply(df, 2, max),
  Mode     = apply(df, 2, calc_mode),
  IQR      = apply(df, 2, IQR),
  Skewness = apply(df, 2, calc_skewness),
  Kurtosis = apply(df, 2, calc_kurtosis)
)
knitr::kable(desc_stats[, c("Count","Mean","Std","Min","Q25","Median","Q75","Max")], digits = 3, booktabs = TRUE)

	Count	Mean	Std	Min	Q25	Median	Q75	Max
AT	9568	19.651	7.452	1.81	13.510	20.345	25.72	37.11
V	9568	54.306	12.708	25.36	41.740	52.080	66.54	81.56
AP	9568	1013.259	5.939	992.89	1009.100	1012.940	1017.26	1033.30
RH	9568	73.309	14.600	25.56	63.328	74.975	84.83	100.16
PE	9568	454.365	17.067	420.26	439.750	451.550	468.43	495.76

knitr::kable(desc_stats[, c("Mode","IQR","Skewness","Kurtosis")], digits = 3, booktabs = TRUE)

	Mode	IQR	Skewness	Kurtosis
AT	25.21	12.210	-0.136	-1.038
V	41.17	24.800	0.199	-1.444
AP	1013.88	8.160	0.265	0.094
RH	100.09	21.502	-0.432	-0.445
PE	468.80	28.680	0.307	-1.049

Interpretasi : Seluruh sheet memiliki struktur variabel yang sama; hasil pemeriksaan kesamaan data setelah pengurutan: 5 dari 5 sheet identik dengan sheet pertama, masing-masing dengan 9568 baris. Pada sheet pertama, PE memiliki rata-rata 454,365 dengan simpangan baku 17,067. Nilai skewness berkisar antara -0,432 hingga 0,307 dan kurtosis berlebih antara -1,444 hingga 0,094.

1.2 Pemeriksaan Kualitas dan Distribusi Data

colSums(is.na(df))

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

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

## 
## Jumlah baris duplikat: 41

p_list <- lapply(colnames(df), function(col) {
  ggplot(df, aes(x = .data[[col]])) +
    geom_histogram(aes(y = after_stat(density)), fill = "steelblue", bins = 30, alpha = 0.7) +
    geom_density(color = "darkblue", linewidth = 0.8) +
    labs(title = paste("Distribusi", col))
})
grid.arrange(grobs = p_list, ncol = 3)

p_box <- lapply(colnames(df), function(col) {
  ggplot(df, aes(y = .data[[col]])) + geom_boxplot(fill = "lightcoral") + labs(title = col)
})
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
  sum(s < q1 - 1.5 * iqr | s > q3 + 1.5 * iqr)
}
out_counts <- sapply(df, iqr_outlier_count)
cat("Jumlah outlier (aturan IQR 1.5x) per kolom:\n"); print(out_counts)

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

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

Interpretasi : Hasil pemeriksaan menunjukkan 0 missing value dan 41 baris duplikat. Berdasarkan aturan IQR, outlier teridentifikasi pada AP (88 observasi), RH (12 observasi), sedangkan AT, V, PE tidak menunjukkan outlier. Histogram dan boxplot menunjukkan bentuk sebaran serta nilai ekstrem tiap variabel secara visual.

1.3 Hubungan Antarvariabel dan Pemeriksaan Asumsi Awal

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

p_reg <- lapply(features, function(col) {
  ggplot(df, aes(x = .data[[col]], y = .data[[target]])) +
    geom_point(alpha = 0.15, size = 1) +
    geom_smooth(method = "lm", color = "red", se = TRUE) +
    labs(title = paste(col, "vs", target))
})
grid.arrange(grobs = p_reg, ncol = 4)

# VIF seperti statsmodels: tiap kolom (termasuk konstanta) diregresikan pada kolom lain.
# R^2 dipusatkan bila regresor memuat konstanta, tak dipusatkan bila tidak (aturan statsmodels).
vif_statsmodels <- function(M) {
  sapply(seq_len(ncol(M)), function(i) {
    y <- M[, i]; Z <- M[, -i, drop = FALSE]
    b <- qr.solve(Z, y)
    res <- y - Z %*% b
    has_const <- any(apply(Z, 2, function(c) all(c == c[1]) && c[1] != 0))
    tss <- if (has_const) sum((y - mean(y))^2) else sum(y^2)
    1 / (1 - (1 - sum(res^2) / tss))
  })
}
X_check <- cbind(const = 1, as.matrix(df[, features]))
vif_all <- data.frame(variabel = colnames(X_check), VIF = vif_statsmodels(X_check))
vif_df  <- vif_all[vif_all$variabel %in% features, ]      # kolom konstanta tidak diinterpretasikan
tbl(vif_df)

variabel	VIF
AT	5.9776
V	3.9430
AP	1.4526
RH	1.7053

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

# Normalitas residual (sampel maks. 5000)
set.seed(RANDOM_STATE)
sample_resid <- sample(resid0, min(5000, length(resid0)))
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 = resid0), aes(x = resid)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "steelblue", alpha = 0.7) +
  geom_density(color = "darkblue") + labs(title = "Distribusi Residual")
p2 <- ggplot(data.frame(resid = resid0), aes(sample = resid)) +
  stat_qq() + stat_qq_line(color = "red") + labs(title = "QQ-Plot Residual")
grid.arrange(p1, p2, ncol = 2)

# Homoskedastisitas (Breusch-Pagan, studentized)
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)

Interpretasi : Korelasi terhadap PE: AT (r = -0,95), V (r = -0,87), AP (r = 0,52), RH (r = 0,39); korelasi antara AT dan V adalah 0,84. Nilai VIF: AT = 5,978, V = 3,943, AP = 1,453, RH = 1,705; AT memiliki VIF tertinggi, yaitu melebihi ambang 5 (James et al., 2021). Pada model linear awal, uji Shapiro–Wilk menghasilkan \(W=\) 0,9795 dengan p-value 2,73e-26 (residual menyimpang dari normalitas pada taraf 5%), dan uji Breusch–Pagan menghasilkan \(\mathrm{LM}=\) 33,9302 dengan p-value 7,7e-07 (terdapat indikasi heteroskedastisitas). Temuan ini dicatat sebagai bahan interpretasi dan menjadi pertimbangan untuk memakai metode yang mampu menangkap hubungan nonlinier.

1.4 Fungsi Bantu Pemodelan

Fungsi-fungsi berikut menyalin fungsi bantu pada notebook: kriteria informasi, metrik, pembagian data, pembentukan fold, standardisasi, dan pemasangan OLS berbasis SVD.

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)
  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)

# Permutasi penuh; n_test indeks pertama = uji, sisanya = latih
make_train_test_split <- function(n, test_size = 0.20, seed = 42) {
  set.seed(seed)
  shuffled <- sample(n)
  n_test <- ceiling(n * test_size)
  list(train = shuffled[(n_test + 1):n], test = shuffled[1:n_test])
}

# Permutasi penuh; nomor fold ditempel ROUND-ROBIN (rep(1:K, length.out = n))
make_folds <- function(n, K = 10, seed = 42) {
  set.seed(seed)
  shuffled <- sample(n)
  fold_id <- rep(1:K, length.out = n)
  lapply(1:K, function(m) sort(shuffled[fold_id == m]))
}

# Standardisasi: sd sampel (pembagi n-1); scale = 1 bila sd = 0 / NA
fit_scaler <- function(X) {
  X <- as.matrix(X)
  center <- colMeans(X)
  scale  <- apply(X, 2, sd)
  scale[scale == 0 | is.na(scale)] <- 1.0
  list(center = center, scale = scale)
}
apply_scaler <- function(X, center, scale) {
  sweep(sweep(as.matrix(X), 2, center, "-"), 2, scale, "/")
}

# OLS dengan intercept: pusatkan data lalu kuadrat terkecil via SVD (pseudo-inverse),
# setara sklearn.LinearRegression. Aman terhadap kolom polynomial yang berkondisi buruk.
ols_fit <- function(Xm, y) {
  Xm <- as.matrix(Xm)
  xm <- colMeans(Xm); ym <- mean(y)
  Xc <- sweep(Xm, 2, xm, "-")
  sv <- svd(Xc)
  tol <- .Machine$double.eps * max(dim(Xc)) * sv$d[1]
  pos <- sv$d > tol
  b <- sv$v[, pos, drop = FALSE] %*% ((t(sv$u[, pos, drop = FALSE]) %*% (y - ym)) / sv$d[pos])
  list(coef = as.numeric(b), intercept = ym - sum(xm * b))
}
ols_predict <- function(fit, Xm) as.numeric(as.matrix(Xm) %*% fit$coef + fit$intercept)

# Fitur polynomial derajat total d dengan suku interaksi, tanpa kolom konstanta
poly_design <- function(X, d) {
  as.matrix(do.call(polym, c(as.list(as.data.frame(X)), list(degree = d, raw = TRUE))))
}

2.1 Polynomial Regression

run_polynomial <- function(Xtr, ytr, Xte, yte, folds) {
  n <- length(ytr)
  rows <- list(); fold_scores <- list()
  Xp_all <- lapply(1:3, function(d) list(tr = poly_design(Xtr, d), te = poly_design(Xte, d)))

  for (degree in 1:3) {
    Xp  <- Xp_all[[degree]]$tr
    fit <- ols_fit(Xp, ytr)
    rss <- sum((ytr - ols_predict(fit, Xp))^2)
    k   <- ncol(Xp) + 1                                 # +1 intercept
    aicc <- compute_aicc(n, rss, k)

    fold_rmse <- numeric(length(folds))
    for (i in seq_along(folds)) {
      val_idx <- folds[[i]]
      tr_idx  <- setdiff(1:n, val_idx)
      f <- ols_fit(Xp[tr_idx, , drop = FALSE], ytr[tr_idx])
      fold_rmse[i] <- rmse(ytr[val_idx], ols_predict(f, Xp[val_idx, , drop = FALSE]))
    }
    fold_scores[[as.character(degree)]] <- fold_rmse
    rows[[degree]] <- data.frame(degree = degree, k = k, AICc = aicc,
                                 CV_RMSE = sqrt(mean(fold_rmse^2)))
  }
  res <- do.call(rbind, rows)
  degree_aic <- res$degree[which.min(res$AICc)]
  degree_cv  <- res$degree[which.min(res$CV_RMSE)]

  evaluate_route <- function(route, degree) {
    fit  <- ols_fit(Xp_all[[degree]]$tr, ytr)
    pred <- ols_predict(fit, Xp_all[[degree]]$te)
    row  <- res[res$degree == degree, ]
    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_R2 = r2_score(yte, pred))
  }
  list(table = res,
       final = list(AIC = evaluate_route("AIC", degree_aic),
                    CV  = evaluate_route("CV",  degree_cv)),
       cv_fold_rmse = fold_scores)
}

2.2 Nonlinear Least Squares (NLS)

NLS_FORMS <- c("exp(AT)", "exp(AT)+exp(V)", "exp(linear)")

# Fungsi rata-rata f(Z; theta) pada fitur terstandarisasi Z (kolom bernama AT, V, AP, RH)
nls_function <- function(name) {
  if (name == "exp(AT)")
    return(function(Z, b) b[1] + b[2] * exp(b[3] * Z[, "AT"]) + b[4] * Z[, "V"] +
                          b[5] * Z[, "AP"] + b[6] * Z[, "RH"])
  if (name == "exp(AT)+exp(V)")
    return(function(Z, b) b[1] + b[2] * exp(b[3] * Z[, "AT"]) + b[4] * exp(b[5] * Z[, "V"]) +
                          b[6] * Z[, "AP"] + b[7] * Z[, "RH"])
  function(Z, b) exp(b[1] + b[2] * Z[, "AT"] + b[3] * Z[, "V"] + b[4] * Z[, "AP"] + b[5] * Z[, "RH"])
}

# Nilai awal tetap: koefisien linear dari OLS, pengali eksponen = 1, eksponen = -0.5
nls_start_values <- function(name, Z, y) {
  b <- ols_fit(Z, y)$coef
  names(b) <- colnames(Z)
  if (name == "exp(AT)")
    return(c(mean(y), 1.0, -0.5, b[["V"]], b[["AP"]], b[["RH"]]))
  if (name == "exp(AT)+exp(V)")
    return(c(mean(y), 1.0, -0.5, 1.0, -0.5, b[["AP"]], b[["RH"]]))
  c(log(mean(y)), 0, 0, 0, 0)
}

# Estimasi NLS: Levenberg-Marquardt (MINPACK, seperti scipy curve_fit). Standardisasi hanya dari data fit.
# Gagal konvergen -> stop() (ditangkap pemanggil, setara RuntimeError pada Python).
nls_fit <- function(name, X_raw, y_raw) {
  y_raw <- as.numeric(y_raw)
  sc <- fit_scaler(X_raw)
  Z  <- apply_scaler(X_raw, sc$center, sc$scale)
  f  <- nls_function(name)
  fit <- minpack.lm::nls.lm(
    par = unname(nls_start_values(name, Z, y_raw)),
    fn  = function(p) y_raw - f(Z, p),
    control = minpack.lm::nls.lm.control(maxfev = 20000, maxiter = 1024)
  )
  if (!(fit$info %in% 1:4)) stop("NLS tidak konvergen")   # info 1-4 = konvergen
  theta <- fit$par
  list(theta = theta,
       predict = function(X_new) f(apply_scaler(X_new, sc$center, sc$scale), theta))
}

run_nls <- function(Xtr, ytr, Xte, yte, folds) {
  ytr <- as.numeric(ytr); yte <- as.numeric(yte)
  n <- length(ytr)
  rows <- list(); fold_scores <- list()

  for (name in NLS_FORMS) {
    full <- tryCatch({
      m <- nls_fit(name, Xtr, ytr)
      rss <- sum((ytr - m$predict(Xtr))^2)
      list(k = length(m$theta), aicc = compute_aicc(n, rss, length(m$theta)))
    }, error = function(e) {
      cat(sprintf("  [peringatan] NLS '%s' tidak konvergen pada data latih penuh\n", name))
      list(k = NA, aicc = NA)
    })

    fold_rmse <- numeric(length(folds))
    for (i in seq_along(folds)) {
      val_idx <- folds[[i]]
      tr_idx  <- setdiff(1:n, val_idx)
      fold_rmse[i] <- tryCatch({
        m <- nls_fit(name, Xtr[tr_idx, , drop = FALSE], ytr[tr_idx])
        rmse(ytr[val_idx], m$predict(Xtr[val_idx, , drop = FALSE]))
      }, error = function(e) {
        cat(sprintf("  [peringatan] NLS '%s' tidak konvergen pada fold %d\n", name, i))
        NA_real_
      })
    }
    fold_scores[[name]] <- fold_rmse
    rows[[name]] <- data.frame(form = name, k = full$k, AICc = full$aicc,
                               CV_RMSE = sqrt(mean(fold_rmse^2)))     # sama seperti Polynomial
  }
  res <- do.call(rbind, rows); rownames(res) <- NULL
  form_aic <- res$form[which.min(res$AICc)]
  form_cv  <- res$form[which.min(res$CV_RMSE)]

  evaluate_route <- function(route, name) {
    m <- nls_fit(name, Xtr, ytr)
    pred <- m$predict(Xte)
    row <- res[res$form == name, ]
    list(config = paste0("form=", name),
         AICc    = if (route == "AIC") row$AICc else NA,
         CV_RMSE = if (route == "CV")  row$CV_RMSE else NA,
         test_RMSE = rmse(yte, pred), test_R2 = r2_score(yte, pred))
  }
  list(table = res,
       final = list(AIC = evaluate_route("AIC", form_aic),
                    CV  = evaluate_route("CV",  form_cv)),
       cv_fold_rmse = fold_scores[[form_cv]])
}

2.3 Generalized Additive Model (GAM)

run_gam <- function(Xtr, ytr, Xte, yte, folds, lambda_grid = NULL, k_spline = 20) {
  if (is.null(lambda_grid)) lambda_grid <- 10^seq(-3, 3, length.out = 10)

  df_tr <- data.frame(y = as.numeric(ytr), Xtr)
  df_te <- data.frame(y = as.numeric(yte), Xte)
  n_tr  <- nrow(df_tr)
  gam_formula <- as.formula(paste("y ~", paste(
    sprintf("s(%s, bs = 'ps', k = %d)", colnames(Xtr), k_spline), collapse = " + ")))

  # Rute AICc: lambda dipilih GCV; AICc dengan k = edf
  gam_gcv <- gam(gam_formula, data = df_tr, method = "GCV.Cp")
  rss  <- sum(residuals(gam_gcv)^2)
  edof <- sum(gam_gcv$hat)
  aicc_gcv   <- compute_aicc(n_tr, rss, edof)
  lambda_gcv <- mean(gam_gcv$sp)

  # Rute CV: lambda sama untuk semua smooth, grid 10^seq(-3,3,length=10)
  cv_rows <- list(); fold_scores <- list()
  for (j in seq_along(lambda_grid)) {
    lam <- lambda_grid[j]
    sp_vec <- rep(lam, ncol(Xtr))
    fold_rmse <- numeric(length(folds))
    for (i in seq_along(folds)) {
      val_idx <- folds[[i]]
      tr_idx  <- setdiff(1:n_tr, val_idx)
      fit  <- gam(gam_formula, data = df_tr[tr_idx, ], sp = sp_vec)
      pred <- as.numeric(predict(fit, newdata = df_tr[val_idx, ]))
      fold_rmse[i] <- rmse(df_tr$y[val_idx], pred)
    }
    fold_scores[[j]] <- fold_rmse
    cv_rows[[j]] <- data.frame(lambda = lam, CV_RMSE = mean(fold_rmse))   # rata-rata RMSE fold
  }
  cv_df  <- do.call(rbind, cv_rows)
  best_j <- which.min(cv_df$CV_RMSE)
  best_cv <- cv_df[best_j, ]

  pred_test_gcv <- as.numeric(predict(gam_gcv, newdata = df_te))
  gam_cv <- gam(gam_formula, data = df_tr, sp = rep(best_cv$lambda, ncol(Xtr)))
  pred_test_cv <- as.numeric(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 = fold_scores[[best_j]],
       edof = edof, lambda_gcv = lambda_gcv)
}

2.4 Menjalankan Analisis pada Seluruh Sheet

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]]])
  for (sh in sheet_names)
    if (nrow(data_sheets[[sh]]) != n) stop("Semua sheet harus memiliki jumlah baris yang sama untuk eksperimen ini.")

  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, drop = FALSE],
                         X_test  = d[test_idx,  features, drop = FALSE],
                         y_train = d[train_idx, target],
                         y_test  = d[test_idx,  target])
  }

  poly_all <- list(); nls_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 NLS...\n", sh))
    nls_all[[sh]]  <- run_nls(sp$X_train, sp$y_train, sp$X_test, sp$y_test, folds)
    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, nls_all = nls_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 NLS...
## [Sheet1] menjalankan GAM...
## [Sheet2] menjalankan Polynomial...
## [Sheet2] menjalankan NLS...
## [Sheet2] menjalankan GAM...
## [Sheet3] menjalankan Polynomial...
## [Sheet3] menjalankan NLS...
## [Sheet3] menjalankan GAM...
## [Sheet4] menjalankan Polynomial...
## [Sheet4] menjalankan NLS...
## [Sheet4] menjalankan GAM...
## [Sheet5] menjalankan Polynomial...
## [Sheet5] menjalankan NLS...
## [Sheet5] menjalankan GAM...

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

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

n_train <- length(result$splits[[1]]$y_train); n_test <- length(result$splits[[1]]$y_test)

Data pada tiap sheet dibagi menjadi 7654 observasi latih dan 1914 observasi uji, dengan 10 fold pada data latih.

2.5 Tabel Hasil Pemilihan Hyperparameter (Sheet Pertama)

sh1 <- result$sheet_names[1]
cat("Polynomial -", sh1, "\n"); tbl(result$poly_all[[sh1]]$table)

## Polynomial - Sheet1

degree	k	AICc	CV_RMSE
1	5	23318.05	4.5880
2	15	22316.47	4.2963
3	35	21967.29	4.2016

cat("NLS -", sh1, "\n");        tbl(result$nls_all[[sh1]]$table)

## NLS - Sheet1

form	k	AICc	CV_RMSE
exp(AT)	6	22597.75	4.3776
exp(AT)+exp(V)	7	22595.91	4.3769
exp(linear)	5	23105.92	4.5249

cat("GAM (rute CV) -", sh1, "\n"); tbl(result$gam_all[[sh1]]$table)

## GAM (rute CV) - Sheet1

lambda	CV_RMSE
0.0010	4.1783
0.0046	4.1780
0.0215	4.1764
0.1000	4.1721
0.4642	4.1658
2.1544	4.1625
10.0000	4.1620
46.4159	4.1631
215.4435	4.1693
1000.0000	4.1832

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,
                   NLS        = result$nls_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,
                   NLS        = result$nls_all[[sh]]$final$CV$CV_RMSE,
                   GAM        = result$gam_all[[sh]]$final$CV$CV_RMSE)
    best_aic <- names(aicc_vals)[which.min(aicc_vals)]
    best_cv  <- names(cv_vals)[which.min(cv_vals)]
    rows[[sh]] <- data.frame(
      Sheet = sh,
      AICc_Polynomial = aicc_vals[["Polynomial"]], AICc_NLS = aicc_vals[["NLS"]], AICc_GAM = aicc_vals[["GAM"]],
      CV_RMSE_Polynomial = cv_vals[["Polynomial"]], CV_RMSE_NLS = cv_vals[["NLS"]], CV_RMSE_GAM = cv_vals[["GAM"]],
      Model_Terbaik_AIC = best_aic, Model_Terbaik_CV = best_cv, Sama = (best_aic == best_cv))
  }
  out <- do.call(rbind, rows); rownames(out) <- NULL
  out
}

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%

tbl(best_model_df[, c("Sheet", "AICc_Polynomial", "AICc_NLS", "AICc_GAM")])

tbl(best_model_df[, c("Sheet", "CV_RMSE_Polynomial", "CV_RMSE_NLS", "CV_RMSE_GAM")])

tbl(best_model_df[, c("Sheet", "Model_Terbaik_AIC", "Model_Terbaik_CV", "Sama")])

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_NLS, row$AICc_GAM)
    cv   <- c(row$CV_RMSE_Polynomial, row$CV_RMSE_NLS, row$CV_RMSE_GAM)
    res  <- tryCatch(cor.test(aicc, cv, method = "spearman"), error = function(e) NULL)
    rows[[i]] <- data.frame(Sheet = row$Sheet,
                            Spearman_rho = if (is.null(res)) NA_real_ else unname(res$estimate),
                            p_value      = if (is.null(res)) NA_real_ else res$p.value)
  }
  do.call(rbind, rows)
}
rank_corr_df <- build_rank_corr_df(best_model_df)
tbl(rank_corr_df)

3.1 Penentuan Model Akhir per Sheet

isclose <- function(a, b, rtol = 1e-5, atol = 1e-8) abs(a - b) <= (atol + rtol * abs(b))

build_overall_best <- function(result) {
  rows <- list()
  for (sh in result$sheet_names) {
    pt <- result$poly_all[[sh]]$table; nt <- result$nls_all[[sh]]$table
    fp <- result$poly_all[[sh]]$final; fn <- result$nls_all[[sh]]$final; fg <- result$gam_all[[sh]]$final
    deg_c  <- get_config_degree(fp$CV$config)
    form_c <- sub("^form=", "", fn$CV$config)
    cand <- data.frame(
      Model     = c("Polynomial", "NLS", "GAM"),
      Parameter = c(fp$CV$config, fn$CV$config, fg$CV$config),
      AICc      = c(min(pt$AICc, na.rm = TRUE), min(nt$AICc, na.rm = TRUE), fg$AIC$AICc),
      CV_RMSE   = c(fp$CV$CV_RMSE, fn$CV$CV_RMSE, fg$CV$CV_RMSE),
      Test_RMSE = c(fp$CV$test_RMSE, fn$CV$test_RMSE, fg$CV$test_RMSE),
      Test_R2   = c(fp$CV$test_R2, fn$CV$test_R2, fg$CV$test_R2),
      stringsAsFactors = FALSE)
    cand$AICc_Best <- isclose(cand$AICc,    min(cand$AICc,    na.rm = TRUE))
    cand$CV_Best   <- isclose(cand$CV_RMSE, min(cand$CV_RMSE, na.rm = TRUE))
    cand$Both_Best <- cand$AICc_Best & cand$CV_Best
    if (any(cand$Both_Best, na.rm = TRUE)) {
      b <- which(cand$Both_Best); best_idx <- b[which.min(cand$AICc[b])]
      status <- "Terbaik pada AICc dan CV_RMSE"
    } else {
      best_idx <- which.min(cand$CV_RMSE)
      status <- "AICc dan CV_RMSE tidak sepakat; mengikuti minimum CV_RMSE"
    }
    rows[[sh]] <- data.frame(Sheet = sh, Model_Akhir = cand$Model[best_idx], Konfigurasi = cand$Parameter[best_idx],
                             AICc = cand$AICc[best_idx], CV_RMSE = cand$CV_RMSE[best_idx],
                             Test_RMSE = cand$Test_RMSE[best_idx], Test_R2 = cand$Test_R2[best_idx],
                             Status = status, stringsAsFactors = FALSE)
  }
  out <- do.call(rbind, rows); rownames(out) <- NULL
  out
}
overall_df <- build_overall_best(result)
tbl(overall_df[, c("Sheet", "Model_Akhir", "Konfigurasi", "Status")])

Sheet	Model_Akhir	Konfigurasi	Status
Sheet1	GAM	lambda=10	Terbaik pada AICc dan CV_RMSE
Sheet2	GAM	lambda=10	Terbaik pada AICc dan CV_RMSE
Sheet3	GAM	lambda=10	Terbaik pada AICc dan CV_RMSE
Sheet4	GAM	lambda=10	Terbaik pada AICc dan CV_RMSE
Sheet5	GAM	lambda=10	Terbaik pada AICc dan CV_RMSE

tbl(overall_df[, c("Sheet", "AICc", "CV_RMSE", "Test_RMSE", "Test_R2")])

Sheet	AICc	CV_RMSE	Test_RMSE	Test_R2
Sheet1	21878.05	4.1620	4.0119	0.9443
Sheet2	21866.38	4.1657	4.0238	0.9462
Sheet3	21888.98	4.1727	3.9975	0.9453
Sheet4	21694.37	4.1275	4.2185	0.9372
Sheet5	21830.49	4.1547	4.0634	0.9435