ANALISIS PEMILIHAN MODEL DENGAN AIC DAN CROSS-VALIDATION
UNTUK PREDIKSI OUTPUT LISTRIK CCPP

Pemilihan model, validasi silang, dan stabilitas hasil pada lima pengacakan

Dosen Pembimbing:

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

Disusun Oleh:

Rasendriya Nandana Kurniawan (140720260011)
Yuda Taufiqurahman Wenske (140720260017)

PROGRAM STUDI MAGISTER STATISTIKA TERAPAN
DEPARTEMEN STATISTIKA
FAKULTAS MATEMATIKA DAN ILMU PENGETAHUAN ALAM
UNIVERSITAS PADJADJARAN
TAHUN AJARAN 2026/2027 

knitr::opts_chunk$set(warning = FALSE, message = FALSE, fig.width = 9, fig.height = 4.5)

library(readxl)
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(tidyr)
library(ggplot2)
library(mgcv)      # GAM
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
library(ranger)    # Random Forest (benchmark)
library(car)       # VIF
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(lmtest)    # Breusch-Pagan, coeftest
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich)  # SE robust (HC3)
library(knitr)

theme_set(theme_minimal(base_size = 11))

tabel <- function(x, caption = NULL, format = "html",
                  table.attr = 'class="report-table"', ...) {
  out <- as.character(knitr::kable(x, format = "html", caption = caption,
                                   table.attr = table.attr, ...))
  out <- sub("<table",
             '<table style="border-collapse:collapse;font-size:13px;line-height:1.45;margin:0;"',
             out, fixed = TRUE)
  out <- sub("<caption>",
             '<caption style="caption-side:top;text-align:left;font-weight:700;color:#1f2937;padding:0 0 8px;">',
             out, fixed = TRUE)
  out <- gsub("<th>", '<th style=""', out, fixed = TRUE)
  out <- gsub("<td>", '<td style=""', out, fixed = TRUE)
  out <- gsub('<th style="',
              '<th style="background:#e5edf6;color:#1f2937;border:1px solid #9eacbb;padding:7px 10px;font-weight:700;vertical-align:middle;white-space:nowrap;',
              out, fixed = TRUE)
  out <- gsub('<td style="',
              '<td style="border:1px solid #c3ccd6;padding:7px 10px;vertical-align:middle;',
              out, fixed = TRUE)
  knitr::asis_output(paste0('<div style="overflow-x:auto;margin:12px 0 22px;">', out, '</div>'))
}

RANDOM_STATE <- 42
K_FOLDS      <- 10
TEST_SIZE    <- 0.20
B_BOOT       <- 1000
RF_TREES_CV  <- 200
RF_TREES_FINAL <- 500
DATA_FILE    <- "DATASET REGRESI.xlsx"
target       <- "PE"
features     <- c("AT", "V", "AP", "RH")
set.seed(RANDOM_STATE)

1 Pertanyaan penelitian

Pertanyaan utama: seberapa kompleks bentuk hubungan antara kondisi lingkungan (AT, V, AP, RH) dan output listrik (PE) yang didukung oleh data, dan apakah AIC dan 10-fold cross-validation (CV) memilih model yang sama?

Untuk menjawabnya, disusun satu himpunan kandidat model dengan kompleksitas yang meningkat secara bertahap. Setiap model punya peran yang jelas:

Peran Model Dinilai dengan
Baseline M0: regresi linear (OLS) AICc + CV + data uji
Kandidat M1: polinomial aditif derajat 2 AICc + CV + data uji
Kandidat M2: polinomial aditif derajat 3 AICc + CV + data uji
Kandidat M3: polinomial derajat 2 + interaksi AICc + CV + data uji
Kandidat M4: GAM (spline berpenalti, REML) AICc (edf) + CV + data uji
Benchmark Random Forest CV + data uji saja (tanpa likelihood → tanpa AIC)

Semua model statistik (M0–M4) diestimasi dengan likelihood Gaussian pada respons dan data latih yang sama, sehingga AICc-nya dapat dibandingkan. Random Forest tidak memiliki likelihood, sehingga hanya dibandingkan dari sisi akurasi prediksi.

2 Data

sheets <- excel_sheets(DATA_FILE)
data_sheets <- lapply(sheets, function(sh) {
  d <- as.data.frame(read_excel(DATA_FILE, sheet = sh))
  names(d) <- trimws(names(d))
  miss <- setdiff(c(features, target), names(d))
  if (length(miss) > 0) stop(sh, ": kolom tidak ditemukan -> ", paste(miss, collapse = ", "))
  d[, c(features, target)]
})
names(data_sheets) <- sheets

# Apakah setiap sheet berisi data yang sama (hanya urutannya diacak)?
sort_df <- function(d) { d <- d[do.call(order, d), ]; rownames(d) <- NULL; d }
ref <- sort_df(data_sheets[[1]])
verify_df <- data.frame(
  Sheet = sheets,
  n = sapply(data_sheets, nrow),
  Identik_dgn_Sheet1_setelah_sort = sapply(data_sheets, function(d) identical(sort_df(d), ref))
)
tabel(verify_df, caption = "Verifikasi isi tiap sheet dibanding Sheet1 (setelah diurutkan)")
Verifikasi isi tiap sheet dibanding Sheet1 (setelah diurutkan)
Sheet n 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 
df <- data_sheets[[1]]
n_total <- nrow(df)

Unit observasi: satu baris adalah rata-rata per jam kondisi lingkungan dan output listrik bersih CCPP saat beroperasi pada beban penuh (2006–2011; Tüfekci, 2014). Kelima sheet berisi 9568 observasi yang sama, dengan urutan baris berbeda. Karena itu, kelima sheet bukan lima sampel independen. Di analisis ini, setiap sheet digunakan untuk satu pengacakan pembagian latih/uji dan fold CV yang berbeda, sehingga kami bisa mengukur sensitivitas hasil terhadap pembagian data. Ini bisa menunjukkan bias prediksi sistematis pada data uji masing-masing pengacakan, tetapi tidak membuktikan ada atau tidaknya bias pada populasi pembangkit secara umum.

data_dict <- data.frame(
  Variabel = c("AT", "V", "AP", "RH", "PE"),
  Nama     = c("Ambient Temperature", "Exhaust Vacuum", "Ambient Pressure",
               "Relative Humidity", "Net hourly electrical energy output"),
  Satuan   = c("°C", "cm Hg", "mbar", "%", "MW"),
  Peran    = c("Prediktor", "Prediktor", "Prediktor", "Prediktor", "Respons"),
  Rentang  = sapply(df[, c(features, target)], function(x) sprintf("%.2f – %.2f", min(x), max(x)))
)
tabel(data_dict, caption = "Kamus data")
Kamus data
Variabel Nama Satuan Peran Rentang
AT AT Ambient Temperature °C Prediktor 1.81 – 37.11
V V Exhaust Vacuum cm Hg Prediktor 25.36 – 81.56
AP AP Ambient Pressure mbar Prediktor 992.89 – 1033.30
RH RH Relative Humidity % Prediktor 25.56 – 100.16
PE PE Net hourly electrical energy output MW Respons 420.26 – 495.76

Alasan pemilihan prediktor dan pencegahan leakage: keempat prediktor adalah kondisi lingkungan yang relevan untuk menjelaskan output. Tidak ada variabel turunan dari PE yang dimasukkan sebagai prediktor. Analisis ini mengestimasi PE ketika AT, V, AP, dan RH tersedia; untuk klaim peramalan operasional ke depan, ketersediaan nilai prediktor pada saat prediksi tetap harus dipastikan.

3 Eksplorasi data terarah

desc_stats <- data.frame(
  Mean = sapply(df, mean), SD = sapply(df, sd), Min = sapply(df, min),
  Q1 = sapply(df, quantile, 0.25), Median = sapply(df, median),
  Q3 = sapply(df, quantile, 0.75), Max = sapply(df, max)
)
desc_stats_table <- tibble::rownames_to_column(round(desc_stats, 3), var = "Variabel")
tabel(desc_stats_table, format = "html", caption = "Statistik deskriptif setiap variabel",
             table.attr = 'class="report-table"')
Statistik deskriptif setiap variabel
Variabel Mean SD Min Q1 Median Q3 Max
AT 19.651 7.452 1.81 13.510 20.345 25.72 37.11
V 54.306 12.708 25.36 41.740 52.080 66.54 81.56
AP 1013.259 5.939 992.89 1009.100 1012.940 1017.26 1033.30
RH 73.309 14.600 25.56 63.328 74.975 84.83 100.16
PE 454.365 17.067 420.26 439.750 451.550 468.43 495.76 
n_missing <- sum(is.na(df))
n_dup     <- sum(duplicated(df))
iqr_out   <- sapply(df, function(s) {
  q <- quantile(s, c(.25, .75)); i <- diff(q)
  sum(s < q[1] - 1.5 * i | s > q[2] + 1.5 * i)
})
iqr_table <- data.frame(
  Variabel = names(iqr_out),
  Jumlah_outlier_IQR = as.integer(iqr_out),
  row.names = NULL
)
tabel(iqr_table, format = "html", caption = "Jumlah observasi di luar batas IQR",
             table.attr = 'class="report-table"')
Jumlah observasi di luar batas IQR
Variabel Jumlah_outlier_IQR
AT 0
V 0
AP 88
RH 12
PE 0

Tidak ada nilai hilang (0). Ditemukan 41 baris duplikat. Duplikat dipertahankan karena data berupa rata-rata per jam dengan presisi terbatas, sehingga kombinasi nilai yang sama dapat terjadi secara wajar pada jam berbeda. Jumlahnya juga sangat kecil (0.43% dari data). Outlier IQR pada AP dan RH tetap dipertahankan karena nilainya masih dalam rentang fisik yang wajar. Pengaruhnya diperiksa lewat Cook’s distance (bagian Diagnostik model terpilih).

df_long <- df %>% pivot_longer(all_of(features), names_to = "Prediktor", values_to = "Nilai")
ggplot(df_long, aes(Nilai, PE)) +
  geom_point(alpha = 0.08, size = 0.6) +
  geom_smooth(method = "lm", se = FALSE, colour = "grey40", linetype = 2, linewidth = 0.7) +
  geom_smooth(method = "gam", formula = y ~ s(x), se = FALSE, colour = "#D55E00", linewidth = 0.9) +
  facet_wrap(~ Prediktor, scales = "free_x", nrow = 1) +
  labs(title = "PE terhadap tiap prediktor",
       subtitle = "Garis putus-putus = linear; garis oranye = smooth. Selisih keduanya = indikasi non-linearitas",
       x = NULL)

corr_table <- as.data.frame(round(cor(df), 3))
corr_table <- tibble::rownames_to_column(corr_table, var = "Variabel")
tabel(corr_table, format = "html", caption = "Matriks korelasi antarvariabel",
             table.attr = 'class="report-table"')
Matriks korelasi antarvariabel
Variabel AT V AP RH PE
AT 1.000 0.844 -0.508 -0.543 -0.948
V 0.844 1.000 -0.414 -0.312 -0.870
AP -0.508 -0.414 1.000 0.100 0.518
RH -0.543 -0.312 0.100 1.000 0.390
PE -0.948 -0.870 0.518 0.390 1.000 
vif_df <- data.frame(Variabel = names(car::vif(lm(PE ~ AT + V + AP + RH, data = df))),
                     VIF = as.numeric(car::vif(lm(PE ~ AT + V + AP + RH, data = df))))
tabel(vif_df %>% mutate(VIF = round(VIF, 3)), format = "html",
             caption = "Variance Inflation Factor (VIF)",
             table.attr = 'class="report-table"')
Variance Inflation Factor (VIF)
Variabel VIF
AT 5.978
V 3.943
AP 1.453
RH 1.705

AT dan V berkorelasi kuat negatif dengan PE dan juga saling berkorelasi (VIF AT = NA). Multikolinearitas ini moderat (VIF < 10), tetapi membuat koefisien individual pada model linear kurang stabil untuk diinterpretasi. Hal ini dibahas di bagian Inferensi vs prediksi.

4 Metode

4.1 Pembagian data

Setiap sheet sudah berupa pengacakan baris yang berbeda. Untuk memanfaatkan pengacakan itu, 80% baris awal pada setiap sheet menjadi data latih dan 20% sisanya menjadi data uji. Data latih dipakai untuk menghitung AICc dan menjalankan 10-fold CV; data uji hanya dipakai sekali untuk evaluasi akhir. Fold CV yang sama dipakai untuk semua model dalam sheet tersebut, sehingga perbandingan CV berpasangan. Dokumentasi UCI menyebut lima pengacakan tersebut disediakan untuk 5×2-fold CV; jadi rancangan saat ini bukan implementasi persis prosedur validasi asli UCI.

make_split <- function(n, test_size = TEST_SIZE) {
  # Setiap worksheet sudah diacak oleh penyedia dataset. Ambil 80% baris awal
  # sebagai data latih dan 20% sisanya sebagai data uji; antar-sheet, observasi
  # yang masuk ke data uji berbeda karena urutan barisnya berbeda.
  n_train <- floor(n * (1 - test_size))
  list(train = seq_len(n_train), test = seq.int(n_train + 1L, n))
}

make_folds <- function(n, K = K_FOLDS, seed = RANDOM_STATE) {
  if (K > n) stop("Jumlah fold tidak boleh melebihi jumlah observasi latih.")
  set.seed(seed)
  idx <- sample.int(n)
  split(idx, rep(seq_len(K), length.out = n))
}

4.2 Kandidat model

model_specs <- list(
  "M0 Linear (baseline)"     = list(type = "lm",
      f = PE ~ AT + V + AP + RH),
  "M1 Polinomial-2 aditif"   = list(type = "lm",
      f = PE ~ poly(AT, 2) + poly(V, 2) + poly(AP, 2) + poly(RH, 2)),
  "M2 Polinomial-3 aditif"   = list(type = "lm",
      f = PE ~ poly(AT, 3) + poly(V, 3) + poly(AP, 3) + poly(RH, 3)),
  "M3 Polinomial-2 + interaksi" = list(type = "lm",
      f = PE ~ (AT + V + AP + RH)^2 + I(AT^2) + I(V^2) + I(AP^2) + I(RH^2)),
  "M4 GAM"                   = list(type = "gam",
      f = PE ~ s(AT, k = 20, bs = "cr") + s(V, k = 20, bs = "cr") +
          s(AP, k = 20, bs = "cr") + s(RH, k = 20, bs = "cr")),
  "Random Forest (benchmark)" = list(type = "rf", f = PE ~ AT + V + AP + RH)
)
stat_models <- names(model_specs)[sapply(model_specs, `[[`, "type") != "rf"]
rf_name     <- "Random Forest (benchmark)"

M1 dan M2 memakai basis polinomial ortogonal (poly()) agar stabil secara numerik. M3 memuat efek utama, seluruh interaksi dua arah, dan kuadrat masing-masing prediktor sehingga menjadi kandidat polinomial derajat dua yang lebih fleksibel. GAM memakai cubic regression spline (k = 20 per prediktor) dengan smoothing parameter yang diestimasi lewat REML.

4.3 Kriteria

  • AICc dihitung dengan rumus yang sama untuk semua model statistik: \(\text{AIC} = -2\log L + 2k\) dan \(\text{AICc} = \text{AIC} + \frac{2k(k+1)}{n-k-1}\). Untuk OLS, \(k\) = jumlah koefisien + 1 (varians error). Untuk GAM, \(k\) = derajat bebas efektif (edf) + 1, yang diambil dari logLik(). Dari AICc dihitung juga \(\Delta\)AICc dan bobot Akaike.
  • CV_RMSE \(= \sqrt{\frac{1}{K}\sum_k \text{MSE}_k}\), dengan rumus yang sama untuk semua model. Ketidakpastiannya dilaporkan sebagai SE dari RMSE per fold. Aturan 1-SE memilih model paling sederhana yang CV_RMSE-nya masih dalam 1 SE dari CV_RMSE terkecil.
  • Data uji: RMSE, MAE, dan R², dengan interval kepercayaan 95% bootstrap (B = 1000) untuk RMSE.
rmse <- function(y, p) sqrt(mean((y - p)^2))
mae  <- function(y, p) mean(abs(y - p))
r2   <- function(y, p) 1 - sum((y - p)^2) / sum((y - mean(y))^2)

# Tuning mtry memakai OOB error hanya pada data analisis yang tersedia.
# Fungsi ini dipanggil ulang pada setiap fold CV, bukan sekali memakai seluruh data.
tune_rf_mtry <- function(train, num.trees = RF_TREES_CV, seed = RANDOM_STATE) {
  scores <- vapply(seq_along(features), function(m) {
    fit <- ranger(
      PE ~ ., data = train, num.trees = num.trees, mtry = m,
      seed = seed + m, importance = "none"
    )
    fit$prediction.error
  }, numeric(1))
  which.min(scores)
}

fit_model <- function(spec, data, mtry = NULL, importance = "none",
                      seed = RANDOM_STATE, num.trees = RF_TREES_FINAL) {
  switch(spec$type,
    lm  = lm(spec$f, data = data),
    gam = gam(spec$f, data = data, method = "REML"),
    rf  = ranger(spec$f, data = data, num.trees = num.trees,
                 mtry = mtry, seed = seed, importance = importance)
  )
}

predict_model <- function(m, newdata) {
  if (inherits(m, "ranger")) predict(m, data = newdata)$predictions
  else as.numeric(predict(m, newdata = newdata))
}

aicc_of <- function(m) {
  if (inherits(m, "ranger")) return(c(k = NA, logLik = NA, AIC = NA, AICc = NA))
  ll <- logLik(m)
  k  <- attr(ll, "df")
  n  <- nobs(m)
  aic <- -2 * as.numeric(ll) + 2 * k
  c(k = k, logLik = as.numeric(ll), AIC = aic,
    AICc = aic + 2 * k * (k + 1) / (n - k - 1))
}

# Hyperparameter RF dituning di dalam setiap fold menggunakan OOB pada analysis fold.
cv_fold_eval <- function(spec, train, folds, seed = RANDOM_STATE) {
  fold_rmse <- numeric(length(folds))
  fold_mtry <- rep(NA_integer_, length(folds))

  for (j in seq_along(folds)) {
    val_idx <- folds[[j]]
    analysis <- train[-val_idx, , drop = FALSE]
    validation <- train[val_idx, , drop = FALSE]

    if (spec$type == "rf") {
      mtry_j <- tune_rf_mtry(analysis, num.trees = RF_TREES_CV,
                             seed = seed + 100L * j)
      m <- fit_model(spec, analysis, mtry = mtry_j, seed = seed + 200L * j,
                     num.trees = RF_TREES_CV)
      fold_mtry[j] <- mtry_j
    } else {
      m <- fit_model(spec, analysis, seed = seed + 200L * j)
    }

    pred <- predict_model(m, validation)
    fold_rmse[j] <- rmse(validation$PE, pred)
  }

  list(rmse = fold_rmse, mtry = fold_mtry)
}

boot_rmse_ci <- function(err, B = B_BOOT, seed = RANDOM_STATE) {
  set.seed(seed)
  bs <- replicate(B, sqrt(mean(sample(err, size = length(err), replace = TRUE)^2)))
  unname(quantile(bs, c(0.025, 0.975)))
}

boot_mean_ci <- function(x, B = B_BOOT, seed = RANDOM_STATE) {
  set.seed(seed)
  bs <- replicate(B, mean(sample(x, size = length(x), replace = TRUE)))
  unname(quantile(bs, c(0.025, 0.975)))
}

# Satu model set di-fit pada satu sheet/permutasi; seluruh kandidat memakai split dan fold sama.
run_sheet <- function(d, sheet_name, model_seed, fold_seed, keep_models = FALSE) {
  split_idx <- make_split(nrow(d), test_size = TEST_SIZE)
  train <- d[split_idx$train, , drop = FALSE]
  test  <- d[split_idx$test, , drop = FALSE]
  folds <- make_folds(nrow(train), K = K_FOLDS, seed = fold_seed)

  rows <- list(); fold_list <- list(); models <- list(); test_err <- list()
  rf_tuning <- data.frame()

  for (i in seq_along(model_specs)) {
    nm <- names(model_specs)[i]
    spec <- model_specs[[nm]]
    cv <- cv_fold_eval(spec, train, folds, seed = fold_seed + 1000L * i)
    fr <- cv$rmse

    if (spec$type == "rf") {
      mtry_final <- tune_rf_mtry(train, num.trees = RF_TREES_CV,
                                 seed = model_seed + 50000L)
      m <- fit_model(spec, train, mtry = mtry_final,
                     importance = if (keep_models) "permutation" else "none",
                     seed = model_seed + 60000L, num.trees = RF_TREES_FINAL)
      oob_rmse <- sqrt(m$prediction.error)
      rf_tuning <- data.frame(
        Sheet = sheet_name, Fold = seq_along(cv$mtry),
        Selected_mtry = cv$mtry, CV_fold_RMSE = cv$rmse
      )
    } else {
      mtry_final <- NA_integer_
      m <- fit_model(spec, train, seed = model_seed + 60000L + i)
      oob_rmse <- NA_real_
    }

    pred_test <- predict_model(m, test)
    err_test <- test$PE - pred_test  # positif = model cenderung underpredict
    pred_train <- predict_model(m, train)
    ci_rmse <- boot_rmse_ci(err_test, seed = RANDOM_STATE + 10L * i)
    ci_bias  <- boot_mean_ci(err_test, seed = RANDOM_STATE + 10L * i)
    a <- aicc_of(m)

    rows[[nm]] <- data.frame(
      Model = nm, k = unname(a[["k"]]), AICc = unname(a[["AICc"]]),
      CV_RMSE = sqrt(mean(fr^2)), CV_SE = sd(fr) / sqrt(length(fr)),
      Train_RMSE = rmse(train$PE, pred_train),
      Test_RMSE = rmse(test$PE, pred_test),
      Test_RMSE_lo = ci_rmse[1], Test_RMSE_hi = ci_rmse[2],
      Test_MAE = mae(test$PE, pred_test), Test_R2 = r2(test$PE, pred_test),
      Test_Bias = mean(err_test), Test_Bias_lo = ci_bias[1], Test_Bias_hi = ci_bias[2],
      Test_Bias_rel_pct = 100 * mean(err_test) / mean(test$PE),
      RF_mtry_final = mtry_final, RF_OOB_RMSE = oob_rmse
    )

    fold_list[[nm]] <- fr
    test_err[[nm]] <- err_test
    if (keep_models) models[[nm]] <- m
  }

  tab <- do.call(rbind, rows)
  rownames(tab) <- NULL

  # AICc dan bobot Akaike hanya untuk model statistik M0-M4.
  is_stat <- tab$Model %in% stat_models
  tab$dAICc <- NA_real_; tab$w_Akaike <- NA_real_
  tab$dAICc[is_stat] <- tab$AICc[is_stat] - min(tab$AICc[is_stat])
  tab$w_Akaike[is_stat] <- exp(-tab$dAICc[is_stat] / 2) /
    sum(exp(-tab$dAICc[is_stat] / 2))

  st <- tab[is_stat, , drop = FALSE]
  best_cv_idx <- which.min(st$CV_RMSE)
  threshold <- min(st$CV_RMSE) + st$CV_SE[best_cv_idx]
  candidates_1se <- st[st$CV_RMSE <= threshold, , drop = FALSE]

  list(
    sheet = sheet_name, table = tab,
    pick = c(AICc = st$Model[which.min(st$AICc)],
             CV = st$Model[best_cv_idx],
             OneSE = candidates_1se$Model[which.min(candidates_1se$k)]),
    folds = fold_list, test_err = test_err, models = models,
    train = train, test = test, split_idx = split_idx,
    rf_tuning = rf_tuning
  )
}

5 Hasil model selection pada lima pengacakan

Setiap sheet dijalankan sebagai satu replikasi evaluasi terpisah. Karena susunan baris tiap sheet berbeda, holdout 20% terakhir menghasilkan bagian uji berbeda; fold CV menggunakan seed berbeda antar-sheet. Semua model dalam satu sheet memakai split dan fold yang sama. Model M0–M4 dievaluasi dengan AICc, CV, dan data uji; Random Forest hanya dibandingkan dengan metrik prediksi. UCI memang menyediakan lima pengacakan untuk 5×2-fold CV, tetapi di sini kami mempertahankan desain 80/20 + 10-fold CV, sehingga hasil lima sheet dibaca sebagai uji sensitivitas terhadap pembagian data, bukan sebagai lima dataset independen.

res_all <- lapply(seq_along(sheets), function(j) {
  run_sheet(
    d = data_sheets[[sheets[j]]], sheet_name = sheets[j],
    model_seed = RANDOM_STATE + 100L * j,
    fold_seed = RANDOM_STATE + 1000L * j,
    keep_models = (j == 1)
  )
})
names(res_all) <- sheets
res1 <- res_all[[1]]
m_list <- res1$models
n_train <- nrow(res1$train)

stab_tab <- dplyr::bind_rows(lapply(sheets, function(sh) {
  dplyr::mutate(res_all[[sh]]$table, Sheet = sh, .before = 1)
}))

picks <- dplyr::bind_rows(lapply(sheets, function(sh) {
  p <- res_all[[sh]]$pick
  data.frame(Sheet = sh, AICc = unname(p[["AICc"]]),
             CV = unname(p[["CV"]]), OneSE = unname(p[["OneSE"]]))
}))

pick_long <- picks %>% tidyr::pivot_longer(
  cols = c(AICc, CV, OneSE), names_to = "Kriteria", values_to = "Model"
)
pick_frequency <- pick_long %>%
  count(Kriteria, Model, name = "Jumlah_sheet") %>%
  arrange(Kriteria, desc(Jumlah_sheet), Model)

rf_cv_tuning <- dplyr::bind_rows(lapply(sheets, function(sh) res_all[[sh]]$rf_tuning))
rf_tuning_frequency <- rf_cv_tuning %>%
  count(Sheet, Selected_mtry, name = "Jumlah_fold") %>%
  arrange(Sheet, Selected_mtry)

stab_summary <- stab_tab %>% group_by(Model) %>%
  summarise(
    CV_RMSE_mean = mean(CV_RMSE), CV_RMSE_sd = sd(CV_RMSE),
    Test_RMSE_mean = mean(Test_RMSE), Test_RMSE_sd = sd(Test_RMSE),
    Test_MAE_mean = mean(Test_MAE), Test_R2_mean = mean(Test_R2),
    Test_Bias_mean = mean(Test_Bias), Test_Bias_sd = sd(Test_Bias),
    Mean_abs_Test_Bias = mean(abs(Test_Bias)),
    Train_RMSE_mean = mean(Train_RMSE),
    CV_to_Test_gap_mean = mean(Test_RMSE - CV_RMSE),
    .groups = "drop"
  ) %>% arrange(CV_RMSE_mean)

bias_by_sheet <- stab_tab %>%
  select(Sheet, Model, Test_Bias, Test_Bias_lo, Test_Bias_hi, Test_Bias_rel_pct)

rank_agreement <- dplyr::bind_rows(lapply(sheets, function(sh) {
  st <- res_all[[sh]]$table %>% filter(Model %in% stat_models)
  data.frame(
    Sheet = sh,
    Spearman_rank_AICc_CV = cor(st$AICc, st$CV_RMSE, method = "spearman"),
    Pemenang_AICc = st$Model[which.min(st$AICc)],
    Pemenang_CV = st$Model[which.min(st$CV_RMSE)],
    Pemenang_sama = st$Model[which.min(st$AICc)] == st$Model[which.min(st$CV_RMSE)]
  )
}))

# Model pilihan pada sheet pertama digunakan untuk diagnostik rinci dan inferensi ilustratif.
ref_model <- unname(res1$pick[["CV"]])
mtry_rf <- res1$table$RF_mtry_final[res1$table$Model == rf_name]

# Tabel lengkap per sheet; nilai AICc, CV, dan metrik uji tidak dirata-ratakan dahulu.
comparison_by_sheet <- stab_tab %>%
  select(Sheet, Model, k, AICc, dAICc, w_Akaike, CV_RMSE, CV_SE,
         Train_RMSE, Test_RMSE, Test_MAE, Test_R2,
         Test_Bias, Test_Bias_lo, Test_Bias_hi, Test_Bias_rel_pct,
         RF_mtry_final, RF_OOB_RMSE)
tabel(comparison_by_sheet %>% mutate(across(where(is.numeric), ~ round(.x, 4))),
             format = "html", caption = "Perbandingan setiap model pada masing-masing pengacakan",
             table.attr = 'class="report-table wide-table"')
Perbandingan setiap model pada masing-masing pengacakan
Sheet Model k AICc dAICc w_Akaike CV_RMSE CV_SE Train_RMSE Test_RMSE Test_MAE Test_R2 Test_Bias Test_Bias_lo Test_Bias_hi Test_Bias_rel_pct RF_mtry_final RF_OOB_RMSE
Sheet1 M0 Linear (baseline) 6.0000 44862.48 1476.4069 0 4.5344 0.0570 4.5310 4.6617 3.6594 0.9247 -0.0486 -0.2693 0.1526 -0.0107 NA NA
Sheet1 M1 Polinomial-2 aditif 10.0000 43960.43 574.3548 0 4.2758 0.0617 4.2694 4.3732 3.4256 0.9338 -0.0597 -0.2559 0.1343 -0.0132 NA NA
Sheet1 M2 Polinomial-3 aditif 14.0000 43783.76 397.6815 0 4.2254 0.0655 4.2182 4.3267 3.3602 0.9352 -0.0667 -0.2593 0.1237 -0.0147 NA NA
Sheet1 M3 Polinomial-2 + interaksi 16.0000 43819.26 433.1831 0 4.2362 0.0612 4.2269 4.3686 3.3995 0.9339 -0.0737 -0.2724 0.1114 -0.0162 NA NA
Sheet1 M4 GAM 46.7984 43386.07 0.0000 1 4.1170 0.0666 4.0923 4.2147 3.2057 0.9385 -0.0459 -0.2360 0.1472 -0.0101 NA NA
Sheet1 Random Forest (benchmark) NA NA NA NA 3.2459 0.0800 1.4689 3.3276 2.2912 0.9617 -0.0540 -0.1996 0.0986 -0.0119 2 3.2026
Sheet2 M0 Linear (baseline) 6.0000 44853.70 1451.6997 0 4.5312 0.0494 4.5284 4.6720 3.6387 0.9260 0.0249 -0.1849 0.2285 0.0055 NA NA
Sheet2 M1 Polinomial-2 aditif 10.0000 43966.63 564.6257 0 4.2777 0.0536 4.2712 4.3654 3.3926 0.9354 0.0407 -0.1732 0.2181 0.0089 NA NA
Sheet2 M2 Polinomial-3 aditif 14.0000 43783.26 381.2522 0 4.2260 0.0530 4.2181 4.3272 3.3444 0.9365 0.0726 -0.1145 0.2836 0.0160 NA NA
Sheet2 M3 Polinomial-2 + interaksi 16.0000 43838.50 436.4961 0 4.2407 0.0534 4.2322 4.3486 3.3585 0.9359 0.0483 -0.1461 0.2540 0.0106 NA NA
Sheet2 M4 GAM 47.4066 43402.01 0.0000 1 4.1228 0.0571 4.0963 4.1936 3.1788 0.9404 0.0471 -0.1322 0.2298 0.0104 NA NA
Sheet2 Random Forest (benchmark) NA NA NA NA 3.2519 0.0779 1.4658 3.3183 2.2875 0.9627 -0.0137 -0.1577 0.1403 -0.0030 2 3.2028
Sheet3 M0 Linear (baseline) 6.0000 44918.21 1438.2949 0 4.5522 0.0514 4.5475 4.5971 3.7153 0.9270 -0.2428 -0.4389 -0.0378 -0.0535 NA NA
Sheet3 M1 Polinomial-2 aditif 10.0000 44042.31 562.3896 0 4.2992 0.0507 4.2923 4.2838 3.4650 0.9366 -0.2479 -0.4445 -0.0642 -0.0546 NA NA
Sheet3 M2 Polinomial-3 aditif 14.0000 43859.71 379.7940 0 4.2475 0.0506 4.2392 4.2460 3.4039 0.9377 -0.2552 -0.4483 -0.0729 -0.0562 NA NA
Sheet3 M3 Polinomial-2 + interaksi 16.0000 43919.88 439.9579 0 4.2654 0.0518 4.2548 4.2589 3.4240 0.9374 -0.2558 -0.4389 -0.0698 -0.0563 NA NA
Sheet3 M4 GAM 49.6225 43479.92 0.0000 1 4.1417 0.0559 4.1160 4.1138 3.2566 0.9416 -0.1994 -0.3783 -0.0224 -0.0439 NA NA
Sheet3 Random Forest (benchmark) NA NA NA NA 3.3036 0.0605 1.5064 3.2245 2.3538 0.9641 -0.1383 -0.2988 0.0075 -0.0305 2 3.2763
Sheet4 M0 Linear (baseline) 6.0000 45074.30 1499.3439 0 4.5975 0.0709 4.5941 4.4084 3.5753 0.9311 -0.0265 -0.2212 0.1656 -0.0058 NA NA
Sheet4 M1 Polinomial-2 aditif 10.0000 44144.84 569.8923 0 4.3288 0.0799 4.3212 4.1665 3.3722 0.9384 -0.0318 -0.2103 0.1534 -0.0070 NA NA
Sheet4 M2 Polinomial-3 aditif 14.0000 43979.19 404.2408 0 4.2810 0.0787 4.2724 4.1095 3.2997 0.9401 -0.0480 -0.2396 0.1272 -0.0106 NA NA
Sheet4 M3 Polinomial-2 + interaksi 16.0000 44042.92 467.9645 0 4.2999 0.0778 4.2891 4.1189 3.3098 0.9398 -0.0484 -0.2328 0.1467 -0.0107 NA NA
Sheet4 M4 GAM 49.1338 43574.95 0.0000 1 4.1716 0.0784 4.1419 4.0058 3.1614 0.9431 -0.1064 -0.2671 0.0682 -0.0234 NA NA
Sheet4 Random Forest (benchmark) NA NA NA NA 3.3202 0.0945 1.5118 3.0227 2.2134 0.9676 -0.0966 -0.2334 0.0307 -0.0213 2 3.2678
Sheet5 M0 Linear (baseline) 6.0000 45021.47 1486.1306 0 4.5818 0.0672 4.5783 4.4731 3.5663 0.9335 0.2245 0.0248 0.4103 0.0493 NA NA
Sheet5 M1 Polinomial-2 aditif 10.0000 44077.47 542.1312 0 4.3081 0.0695 4.3022 4.2422 3.3854 0.9402 0.1306 -0.0478 0.3134 0.0287 NA NA
Sheet5 M2 Polinomial-3 aditif 14.0000 43910.44 375.1044 0 4.2615 0.0675 4.2533 4.1880 3.3216 0.9417 0.1296 -0.0606 0.3152 0.0285 NA NA
Sheet5 M3 Polinomial-2 + interaksi 16.0000 43972.92 437.5765 0 4.2779 0.0704 4.2696 4.1977 3.3433 0.9415 0.1477 -0.0398 0.3262 0.0325 NA NA
Sheet5 M4 GAM 48.3565 43535.34 0.0000 1 4.1600 0.0718 4.1316 4.0567 3.1595 0.9453 0.1281 -0.0375 0.3197 0.0282 NA NA
Sheet5 Random Forest (benchmark) NA NA NA NA 3.3213 0.0802 1.5028 3.1430 2.2798 0.9672 0.1371 -0.0074 0.2763 0.0301 2 3.2536 
tabel(picks, format = "html", caption = "Model terpilih menurut AICc, CV_RMSE, dan aturan 1-SE",
             table.attr = 'class="report-table"')
Model terpilih menurut AICc, CV_RMSE, dan aturan 1-SE
Sheet AICc CV OneSE
Sheet1 M4 GAM M4 GAM M4 GAM
Sheet2 M4 GAM M4 GAM M4 GAM
Sheet3 M4 GAM M4 GAM M4 GAM
Sheet4 M4 GAM M4 GAM M4 GAM
Sheet5 M4 GAM M4 GAM M4 GAM 
tabel(pick_frequency, format = "html", caption = "Frekuensi model terpilih di lima pengacakan",
             table.attr = 'class="report-table"')
Frekuensi model terpilih di lima pengacakan
Kriteria Model Jumlah_sheet
AICc M4 GAM 5
CV M4 GAM 5
OneSE M4 GAM 5 
tabel(rf_tuning_frequency, format = "html", caption = "Frekuensi mtry terpilih pada fold CV Random Forest",
             table.attr = 'class="report-table"')
Frekuensi mtry terpilih pada fold CV Random Forest
Sheet Selected_mtry Jumlah_fold
Sheet1 2 10
Sheet2 2 10
Sheet3 2 10
Sheet4 2 10
Sheet5 2 10 
tabel(stab_summary %>% mutate(across(where(is.numeric), ~ round(.x, 4))),
             format = "html", caption = "Rata-rata dan simpangan baku metrik pada lima pengacakan",
             table.attr = 'class="report-table wide-table"')
Rata-rata dan simpangan baku metrik pada lima pengacakan
Model CV_RMSE_mean CV_RMSE_sd Test_RMSE_mean Test_RMSE_sd Test_MAE_mean Test_R2_mean Test_Bias_mean Test_Bias_sd Mean_abs_Test_Bias Train_RMSE_mean CV_to_Test_gap_mean
Random Forest (benchmark) 3.2886 0.0369 3.2072 0.1278 2.2851 0.9646 -0.0331 0.1059 0.0879 1.4911 -0.0814
M4 GAM 4.1426 0.0234 4.1169 0.0886 3.1924 0.9418 -0.0353 0.1281 0.1054 4.1156 -0.0257
M2 Polinomial-3 aditif 4.2483 0.0238 4.2395 0.0934 3.3459 0.9383 -0.0335 0.1486 0.1144 4.2403 -0.0088
M3 Polinomial-2 + interaksi 4.2640 0.0265 4.2586 0.1042 3.3670 0.9377 -0.0364 0.1506 0.1148 4.2545 -0.0055
M1 Polinomial-2 aditif 4.2979 0.0221 4.2862 0.0868 3.4082 0.9369 -0.0336 0.1406 0.1021 4.2913 -0.0117
M0 Linear (baseline) 4.5594 0.0293 4.5625 0.1170 3.6310 0.9285 -0.0137 0.1674 0.1135 4.5558 0.0030 
tabel(rank_agreement, format = "html", caption = "Kesepakatan peringkat/pemenang AICc dan CV",
             table.attr = 'class="report-table"')
Kesepakatan peringkat/pemenang AICc dan CV
Sheet Spearman_rank_AICc_CV Pemenang_AICc Pemenang_CV Pemenang_sama
Sheet1 1 M4 GAM M4 GAM TRUE
Sheet2 1 M4 GAM M4 GAM TRUE
Sheet3 1 M4 GAM M4 GAM TRUE
Sheet4 1 M4 GAM M4 GAM TRUE
Sheet5 1 M4 GAM M4 GAM TRUE

Cara membaca bias prediksi: Test_Bias didefinisikan sebagai rata-rata aktual - prediksi. Nilai positif berarti model cenderung underpredict; nilai negatif berarti overpredict. Interval bootstrap 95% yang mencakup nol berarti tidak ada bukti kuat tentang bias rata-rata pada subset uji tersebut, tetapi bukan bukti bahwa model pasti bebas bias. Test_Bias_rel_pct menyatakan bias rata-rata sebagai persentase dari rata-rata PE aktual pada subset uji.

Karena kelima sheet berasal dari observasi yang sama dalam urutan berbeda, simpangan baku pada stab_summary mengukur sensitivitas terhadap pengacakan/pembagian data, bukan variasi antar-populasi. Konsistensi pemenang di pick_frequency adalah indikator stabilitas pemilihan model. Train_RMSE merupakan error pada data yang digunakan untuk fit; untuk menilai generalisasi Random Forest, prioritaskan OOB RMSE, CV_RMSE, dan Test_RMSE.

ggplot(stab_tab, aes(Test_RMSE, factor(Model, levels = rev(names(model_specs))), colour = Sheet)) +
  geom_point(position = position_jitter(height = 0.10, width = 0, seed = 1), size = 2) +
  labs(title = "Test RMSE menurut model dan pengacakan", x = "Test RMSE (MW)", y = NULL) +
  theme(legend.position = "bottom")

ggplot(stab_tab, aes(Test_Bias, factor(Model, levels = rev(names(model_specs))), colour = Sheet)) +
  geom_vline(xintercept = 0, linetype = 2, colour = "grey35") +
  geom_errorbarh(aes(xmin = Test_Bias_lo, xmax = Test_Bias_hi), height = 0.18,
                 position = position_dodge(width = 0.45)) +
  geom_point(position = position_dodge(width = 0.45), size = 2) +
  labs(title = "Bias prediksi pada data uji (aktual − prediksi)",
       subtitle = "Interval bootstrap 95%; nilai positif berarti underprediction",
       x = "Mean Error / Bias (MW)", y = NULL) +
  theme(legend.position = "bottom")

5.1 Rincian satu pengacakan: Sheet pertama

Sheet pertama dipakai untuk tabel rinci, perbandingan berpasangan di data uji, diagnostik, dan pembahasan inferensi. Keputusan akhir tetap dibaca bersama hasil lima pengacakan, bukan dari Sheet pertama saja.

tab1 <- res1$table %>%
  mutate(across(c(AICc, dAICc), ~ round(.x, 2)),
         across(c(k, CV_RMSE, CV_SE, Train_RMSE, Test_RMSE, Test_MAE,
                 Test_Bias, Test_Bias_lo, Test_Bias_hi), ~ round(.x, 3)),
         w_Akaike = signif(w_Akaike, 3), Test_R2 = round(Test_R2, 4),
         Test_RMSE_CI95 = sprintf("[%.3f, %.3f]", Test_RMSE_lo, Test_RMSE_hi),
         Test_Bias_CI95 = sprintf("[%.3f, %.3f]", Test_Bias_lo, Test_Bias_hi)) %>%
  select(Model, k, AICc, dAICc, w_Akaike, CV_RMSE, CV_SE, Train_RMSE,
         Test_RMSE, Test_RMSE_CI95, Test_MAE, Test_R2,
         Test_Bias, Test_Bias_CI95, Test_Bias_rel_pct, RF_mtry_final, RF_OOB_RMSE)
tabel(tab1, format = "html", caption = paste("Hasil lengkap model pada", sheets[1]),
             table.attr = 'class="report-table wide-table"')
Hasil lengkap model pada Sheet1
Model k AICc dAICc w_Akaike CV_RMSE CV_SE Train_RMSE Test_RMSE Test_RMSE_CI95 Test_MAE Test_R2 Test_Bias Test_Bias_CI95 Test_Bias_rel_pct RF_mtry_final RF_OOB_RMSE
M0 Linear (baseline) 6.000 44862.48 1476.41 0 4.534 0.057 4.531 4.662 [4.434, 4.962] 3.659 0.9247 -0.049 [-0.269, 0.153] -0.0107105 NA NA
M1 Polinomial-2 aditif 10.000 43960.43 574.35 0 4.276 0.062 4.269 4.373 [4.110, 4.693] 3.426 0.9338 -0.060 [-0.256, 0.134] -0.0131562 NA NA
M2 Polinomial-3 aditif 14.000 43783.76 397.68 0 4.225 0.065 4.218 4.327 [4.066, 4.663] 3.360 0.9352 -0.067 [-0.259, 0.124] -0.0146973 NA NA
M3 Polinomial-2 + interaksi 16.000 43819.26 433.18 0 4.236 0.061 4.227 4.369 [4.091, 4.740] 3.399 0.9339 -0.074 [-0.272, 0.111] -0.0162419 NA NA
M4 GAM 46.798 43386.07 0.00 1 4.117 0.067 4.092 4.215 [3.932, 4.532] 3.206 0.9385 -0.046 [-0.236, 0.147] -0.0101011 NA NA
Random Forest (benchmark) NA NA NA NA 3.246 0.080 1.469 3.328 [2.987, 3.718] 2.291 0.9617 -0.054 [-0.200, 0.099] -0.0118981 2 3.202604 
picked_table <- data.frame(Kriteria = names(res1$pick), Model_terpilih = unname(res1$pick))
tabel(picked_table, format = "html", caption = "Ringkasan model terpilih pada pengacakan pertama",
             table.attr = 'class="report-table"')
Ringkasan model terpilih pada pengacakan pertama
Kriteria Model_terpilih
AICc M4 GAM
CV M4 GAM
OneSE M4 GAM 
pd <- res1$table %>% mutate(Model = factor(Model, levels = rev(names(model_specs))))
ggplot(pd, aes(CV_RMSE, Model)) +
  geom_errorbarh(aes(xmin = CV_RMSE - CV_SE, xmax = CV_RMSE + CV_SE), height = 0.25) +
  geom_point(size = 2.5, aes(colour = Model == rf_name), show.legend = FALSE) +
  scale_colour_manual(values = c(`FALSE` = "black", `TRUE` = "#0072B2")) +
  labs(title = "CV_RMSE ± 1 SE pada data latih pengacakan pertama",
       x = "CV_RMSE (MW)", y = NULL)

5.2 Apakah perbedaan akurasi data uji meyakinkan?

Selisih RMSE uji setiap model terhadap model statistik terpilih menurut CV pada sheet pertama dihitung dengan bootstrap berpasangan (observasi uji yang sama diresampling untuk kedua model):

paired_diff <- function(e_a, e_b, B = B_BOOT, seed = RANDOM_STATE) {
  set.seed(seed)
  n <- length(e_a)
  d <- replicate(B, {
    i <- sample.int(n, size = n, replace = TRUE)
    sqrt(mean(e_a[i]^2)) - sqrt(mean(e_b[i]^2))
  })
  c(est = sqrt(mean(e_a^2)) - sqrt(mean(e_b^2)),
    quantile(d, c(.025, .975)))
}

diff_df <- do.call(rbind, lapply(setdiff(names(model_specs), ref_model), function(nm) {
  x <- paired_diff(res1$test_err[[nm]], res1$test_err[[ref_model]])
  data.frame(Model = nm, Selisih_RMSE = x[["est"]],
             CI_lo = x[["2.5%"]], CI_hi = x[["97.5%"]],
             Signifikan = x[["2.5%"]] > 0 | x[["97.5%"]] < 0)
}))
diff_df$Pembanding <- ref_model
tabel(diff_df %>% mutate(across(c(Selisih_RMSE, CI_lo, CI_hi), ~ round(.x, 3))),
             format = "html", caption = paste("Selisih RMSE uji relatif terhadap", ref_model),
             table.attr = 'class="report-table"')
Selisih RMSE uji relatif terhadap M4 GAM
Model Selisih_RMSE CI_lo CI_hi Signifikan Pembanding
M0 Linear (baseline) 0.447 0.344 0.554 TRUE M4 GAM
M1 Polinomial-2 aditif 0.158 0.114 0.205 TRUE M4 GAM
M2 Polinomial-3 aditif 0.112 0.075 0.152 TRUE M4 GAM
M3 Polinomial-2 + interaksi 0.154 0.106 0.202 TRUE M4 GAM
Random Forest (benchmark) -0.887 -1.019 -0.758 TRUE M4 GAM

Selisih positif berarti model pembanding memiliki RMSE lebih tinggi daripada M4 GAM. Jika interval mencakup nol, perbedaan akurasi pada subset uji pertama belum meyakinkan. Hasil ini merupakan pelengkap; stabilitas lintas pengacakan tetap dilihat di tabel sebelumnya.

6 Diagnostik model terpilih

Diagnostik dilakukan pada model terpilih menurut CV (M4 GAM) dan pada baseline linear sebagai pembanding.

m_sel  <- res1$models[[ref_model]]
m_base <- res1$models[["M0 Linear (baseline)"]]
train1 <- res1$train

diag_data <- function(m, label) {
  data.frame(Model = label, fitted = as.numeric(fitted(m)),
             resid = as.numeric(residuals(m)), cook = as.numeric(cooks.distance(m)))
}
dd <- rbind(diag_data(m_base, "M0 Linear (baseline)"), diag_data(m_sel, ref_model))
ggplot(dd, aes(fitted, resid)) +
  geom_point(alpha = 0.1, size = 0.6) +
  geom_hline(yintercept = 0, colour = "red") +
  geom_smooth(method = "gam", formula = y ~ s(x), se = FALSE, colour = "#D55E00") +
  facet_wrap(~ Model) + labs(title = "Residual vs fitted", x = "Fitted PE", y = "Residual")

ggplot(dd, aes(sample = resid)) + stat_qq(alpha = 0.2, size = 0.6) + stat_qq_line(colour = "red") +
  facet_wrap(~ Model) + labs(title = "QQ-plot residual")

# Breusch-Pagan versi umum (nR^2 dari regresi residual^2 ke prediktor) agar berlaku untuk lm dan GAM
bp_general <- function(resid, X) {
  aux <- lm(resid^2 ~ ., data = X); n <- length(resid)
  stat <- n * summary(aux)$r.squared
  c(BP = stat, p = pchisq(stat, df = ncol(X), lower.tail = FALSE))
}
set.seed(RANDOM_STATE)
sw <- function(r) shapiro.test(sample(r, min(5000, length(r))))$p.value

assump <- rbind(
  data.frame(Model = "M0 Linear (baseline)", t(bp_general(residuals(m_base), train1[, features])),
             Shapiro_p = sw(residuals(m_base))),
  data.frame(Model = ref_model, t(bp_general(residuals(m_sel), train1[, features])),
             Shapiro_p = sw(residuals(m_sel))))
tabel(assump %>% mutate(BP = round(BP, 2), across(c(p, Shapiro_p), ~ signif(.x, 3))),
      caption = "Uji asumsi residual: Breusch-Pagan versi umum dan Shapiro-Wilk")
Uji asumsi residual: Breusch-Pagan versi umum dan Shapiro-Wilk
Model BP p Shapiro_p
M0 Linear (baseline) 35.16 4e-07 0
M4 GAM 51.24 0e+00 0

Dengan n ≈ 7654, uji formal hampir selalu menolak H0 meskipun penyimpangannya kecil. Karena itu, keputusan didasarkan pada plot. Penyimpangan normalitas dan heteroskedastisitas tidak mengganggu prediksi titik maupun CV, tetapi memengaruhi SE/p-value klasik. Inilah alasan inferensi baseline di bagian Inferensi memakai SE robust HC3.

thr_cook <- 4 / n_train
cook_sum <- dd %>% group_by(Model) %>%
  summarise(n_cook_gt_4n = sum(cook > thr_cook), max_cook = max(cook), .groups = "drop")
tabel(cook_sum, caption = "Ringkasan Cook's distance per model")
Ringkasan Cook’s distance per model
Model n_cook_gt_4n max_cook
M0 Linear (baseline) 343 0.0185220
M4 GAM 272 0.0190888 
ggplot(dd %>% group_by(Model) %>% mutate(idx = row_number()), aes(idx, cook)) +
  geom_segment(aes(xend = idx, yend = 0), alpha = 0.4) +
  geom_hline(yintercept = thr_cook, colour = "red", linetype = 2) +
  facet_wrap(~ Model, scales = "free_y") + labs(title = "Cook's distance (garis merah = 4/n)", x = "Observasi", y = "Cook's D")

# Sensitivitas: refit model terpilih tanpa observasi Cook's D > 4/n, lalu bandingkan RMSE uji
infl <- which(cooks.distance(m_sel) > thr_cook)
train_noinfl <- if (length(infl) > 0) train1[-infl, , drop = FALSE] else train1
m_sel_noinfl <- fit_model(model_specs[[ref_model]], train_noinfl,
                          seed = RANDOM_STATE + 91001)
sens_infl <- data.frame(
  Skenario = c("Semua data latih", sprintf("Tanpa %d observasi berpengaruh", length(infl))),
  Test_RMSE = c(rmse(res1$test$PE, predict_model(m_sel, res1$test)),
                rmse(res1$test$PE, predict_model(m_sel_noinfl, res1$test))))
tabel(sens_infl, caption = "Sensitivitas: RMSE uji dengan dan tanpa observasi berpengaruh")
Sensitivitas: RMSE uji dengan dan tanpa observasi berpengaruh
Skenario Test_RMSE
Semua data latih 4.214690
Tanpa 272 observasi berpengaruh 4.223007 
if (inherits(m_sel, "gam")) {
  kcheck_table <- as.data.frame(k.check(m_sel))
  kcheck_table <- tibble::rownames_to_column(kcheck_table, var = "Smooth")
  tabel(kcheck_table, format = "html",
               caption = "Pemeriksaan dimensi basis GAM (k-index < 1 dan p kecil perlu ditinjau)",
               table.attr = 'class="report-table"')

  concurvity_table <- as.data.frame(round(concurvity(m_sel, full = TRUE), 3))
  concurvity_table <- tibble::rownames_to_column(concurvity_table, var = "Komponen")
  tabel(concurvity_table, format = "html",
               caption = "Concurvity GAM (nilai tinggi perlu diwaspadai)",
               table.attr = 'class="report-table"')
}
Concurvity GAM (nilai tinggi perlu diwaspadai)
Komponen para s(AT) s(V) s(AP) s(RH)
worst 0 0.886 0.830 0.454 0.534
observed 0 0.848 0.783 0.196 0.506
estimate 0 0.108 0.120 0.062 0.051

7 Inferensi vs prediksi

7.1 Inferensi: apa yang bisa dikatakan tentang hubungan

Baseline linear dengan SE robust (HC3). Koefisiennya mudah ditafsirkan sebagai perubahan PE per satuan prediktor dengan prediktor lain tetap, tetapi hanya valid jika hubungannya memang linear.

ct <- coeftest(m_base, vcov. = vcovHC(m_base, type = "HC3"))
ci <- coefci(m_base, vcov. = vcovHC(m_base, type = "HC3"))
tabel(data.frame(Estimasi = ct[, 1], SE_HC3 = ct[, 2], CI_lo = ci[, 1], CI_hi = ci[, 2], p = ct[, 4]) %>%
  round(4),
      caption = "Koefisien baseline linear dengan SE robust HC3")
Koefisien baseline linear dengan SE robust HC3
Estimasi SE_HC3 CI_lo CI_hi p
(Intercept) 454.9770 10.5968 434.2045 475.7496 0
AT -1.9923 0.0179 -2.0273 -1.9573 0
V -0.2271 0.0081 -0.2430 -0.2113 0
AP 0.0618 0.0103 0.0417 0.0820 0
RH -0.1606 0.0046 -0.1696 -0.1516 0

Uji F bersarang. Apakah penambahan kelengkungan dan interaksi didukung data?

# Model bersarang: M0, M1, dan M2
m0 <- m_list[["M0 Linear (baseline)"]]
m1 <- m_list[["M1 Polinomial-2 aditif"]]
m2 <- m_list[["M2 Polinomial-3 aditif"]]
m3 <- m_list[["M3 Polinomial-2 + interaksi"]]

anova_012 <- as.data.frame(anova(m0, m1, m2))
anova_012 <- tibble::rownames_to_column(
  anova_012, var = "Model"
)

tabel(
  anova_012,
  format = "html",
  caption = "Perbandingan model bersarang: M0, M1, dan M2",
  table.attr = 'class="report-table"'
)
Perbandingan model bersarang: M0, M1, dan M2
Model Res.Df RSS Df Sum of Sq F Pr(>F)
1 7649 157133.2 NA NA NA NA
2 7645 139517.9 4 17615.341 247.07647 0
3 7641 136191.5 4 3326.423 46.65711 0 
# Perbandingan M1 dan M3
anova_13 <- as.data.frame(anova(m1, m3))
anova_13 <- tibble::rownames_to_column(
  anova_13, var = "Model"
)

tabel(
  anova_13,
  format = "html",
  caption = "Perbandingan model bersarang: M1 dan M3",
  table.attr = 'class="report-table"'
)
Perbandingan model bersarang: M1 dan M3
Model Res.Df RSS Df Sum of Sq F Pr(>F)
1 7645 139517.9 NA NA NA NA
2 7639 136752.8 6 2765.036 25.74244 0

GAM: signifikansi dan bentuk tiap fungsi smooth. edf ≈ 1 berarti hubungannya praktis linear, sedangkan edf yang besar berarti hubungannya melengkung.

m_gam <- m_list[["M4 GAM"]]
gam_s <- as.data.frame(round(summary(m_gam)$s.table, 4))
gam_s <- tibble::rownames_to_column(gam_s, var = "Smooth")
tabel(gam_s, format = "html", caption = "Signifikansi dan derajat bebas efektif smooth GAM",
             table.attr = 'class="report-table"')
Signifikansi dan derajat bebas efektif smooth GAM
Smooth edf Ref.df F p-value
s(AT) 10.5247 12.8398 1030.1359 0
s(V) 16.3571 17.8190 85.0548 0
s(AP) 9.7080 11.7820 22.2681 0
s(RH) 4.8852 6.1359 102.6389 0 
plot(m_gam, pages = 1, shade = TRUE, seWithMean = TRUE, scale = 0, residuals = FALSE)

7.2 Prediksi: apa yang bisa dan tidak bisa dikatakan dari benchmark Random Forest

m_rf <- m_list[[rf_name]]
imp <- sort(importance(m_rf), decreasing = TRUE)
importance_table <- data.frame(Variabel = names(imp),
                               Permutation_importance = round(as.numeric(imp), 2))
tabel(importance_table, format = "html", caption = "Permutation importance Random Forest",
             table.attr = 'class="report-table"')
Permutation importance Random Forest
Variabel Permutation_importance
AT 313.35
V 82.87
AP 14.12
RH 9.50 
# Pilihan mtry akhir untuk pengacakan pertama dan OOB RMSE-nya
res1$table %>% filter(Model == rf_name) %>%
  select(Model, RF_mtry_final, RF_OOB_RMSE, CV_RMSE, Test_RMSE) %>%
  mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
  tabel(format = "html", caption = "Kinerja Random Forest pada pengacakan pertama",
               table.attr = 'class="report-table"')
Kinerja Random Forest pada pengacakan pertama
Model RF_mtry_final RF_OOB_RMSE CV_RMSE Test_RMSE
Random Forest (benchmark) 2 3.203 3.246 3.328 
# Partial dependence AT: RF vs GAM, dihitung dengan cara yang sama agar sebanding
set.seed(RANDOM_STATE)
pdp_base <- train1[sample(nrow(train1), 1000), ]
grid_at  <- seq(quantile(train1$AT, .01), quantile(train1$AT, .99), length.out = 40)
pdp <- do.call(rbind, lapply(grid_at, function(g) {
  nd <- pdp_base; nd$AT <- g
  data.frame(AT = g, RF = mean(predict_model(m_rf, nd)), GAM = mean(predict_model(m_gam, nd)))
})) %>% pivot_longer(-AT, names_to = "Model", values_to = "PE")
ggplot(pdp, aes(AT, PE, colour = Model)) + geom_line(linewidth = 1) +
  scale_colour_manual(values = c(GAM = "#D55E00", RF = "#0072B2")) +
  labs(title = "Partial dependence PE terhadap AT", y = "Rata-rata prediksi PE (MW)")

Perbedaan inferensi statistik dan predictive accuracy:

  • Model M0–M4 menjawab bagaimana PE berubah terhadap tiap prediktor. Jawabannya datang dengan ketidakpastian yang terukur (CI, uji F, signifikansi smooth) dan dapat dibandingkan lewat AICc karena semuanya punya likelihood.
  • Random Forest hanya menjawab seberapa akurat PE dapat diprediksi. Jika RMSE-nya lebih kecil, itu berarti ada struktur (misalnya interaksi tingkat tinggi) yang belum ditangkap kandidat model statistik. Namun itu tidak memberi koefisien, CI, atau uji hipotesis.
  • Yang tidak boleh disimpulkan dari benchmark: (1) permutation importance bukan efek kausal dan bisa terdistorsi oleh korelasi AT–V; (2) RMSE yang lebih kecil tidak berarti model “benar”, melainkan hanya lebih akurat pada distribusi data ini; (3) RF tidak dapat diekstrapolasi di luar rentang data latih (prediksinya mendatar); (4) RF tidak bisa diikutkan dalam perbandingan AICc.

8 Pembahasan

cv_mean <- setNames(stab_summary$CV_RMSE_mean, stab_summary$Model)
te_mean <- setNames(stab_summary$Test_RMSE_mean, stab_summary$Model)
n_agree <- sum(rank_agreement$Pemenang_sama)
n_sheet <- nrow(rank_agreement)
pair_count <- function(metric) sum(sapply(sheets, function(sh) {
  tb <- res_all[[sh]]$table
  v  <- tb[[metric]]
  v[tb$Model == "M2 Polinomial-3 aditif"] < v[tb$Model == "M3 Polinomial-2 + interaksi"]
}))
m2_lebih_cv    <- pair_count("CV_RMSE")
m2_lebih_aicc <- pair_count("AICc")
m2_lebih_test <- pair_count("Test_RMSE")
winner_cv <- paste(unique(rank_agreement$Pemenang_CV), collapse = ", ")

Model statistik terbaik. Pada 5 dari 5 pengacakan, model terpilih oleh AICc sama dengan model terpilih oleh CV. Model yang menang pada kedua kriteria adalah M4 GAM. Rata-rata CV_RMSE-nya 4.143 MW, sedangkan polinomial terbaik (M2) 4.248 MW dan baseline linear (M0) 4.559 MW. Ini menunjukkan bahwa hubungan PE dengan kondisi lingkungan memang tidak linear, dan fleksibilitas yang lebih besar dari GAM menangkap bentuk tersebut lebih baik daripada polinomial dengan derajat tetap.

Peran polinomial dan interaksi. Uji F bersarang antara M1 dan M3 signifikan (F = 25.74, p < 0,001), sehingga interaksi memang menurunkan RSS pada data latih. Namun, M2 (aditif derajat 3) mengungguli M3 (derajat 2 dengan interaksi) pada AICc di 5 dari 5 pengacakan, pada CV_RMSE di 5 dari 5, dan pada RMSE data uji di 5 dari 5. Dengan kata lain, kelengkungan yang lebih fleksibel lebih berguna untuk prediksi daripada interaksi dua arah, meskipun interaksi tetap memberi perbaikan dalam sampel.

Random Forest sebagai benchmark. Random Forest memiliki CV_RMSE rata-rata 3.289 MW, lebih rendah daripada GAM. Selisih ini menandakan bahwa ada struktur yang belum sepenuhnya ditangkap oleh kandidat statistik, misalnya interaksi tingkat tinggi. Namun, Random Forest tidak memberi koefisien, interval kepercayaan, atau uji hipotesis, sehingga ia tidak menggantikan model statistik untuk tujuan inferensi.

9 Keterbatasan

  1. Data hanya mencakup kondisi beban penuh pada satu pembangkit, sehingga model tidak berlaku untuk beban parsial atau pembangkit lain.
  2. Tidak ada informasi waktu. Jika observasi berurutan per jam saling berkorelasi, CV acak dapat sedikit terlalu optimis.
  3. AICc hanya membandingkan model di dalam himpunan kandidat. Model terbaik menurut AICc belum tentu model yang benar.
  4. Hasil AICc untuk GAM bergantung pada pendekatan edf. Smoothing parameter yang diestimasi menambah ketidakpastian yang tidak sepenuhnya tercermin di \(k\).
  5. Residual tidak normal dan heteroskedastik, sehingga inferensi klasik pada model polinomial perlu dibaca hati-hati (SE robust dipakai pada baseline).
  6. Random Forest adalah benchmark prediksi, bukan model inferensi yang dibandingkan dengan AICc. Tuning mtry dilakukan hanya pada data analisis yang tersedia: menggunakan OOB error di dalam setiap fold CV dan di data latih untuk fit final. Pendekatan ini mengurangi risiko tuning memakai fold validasi, meskipun performa tetap harus ditafsirkan sebagai hasil pada distribusi data CCPP ini.

10 Kesimpulan

  1. Hubungan PE dengan AT, V, AP, dan RH tidak sepenuhnya linear. Model GAM, yang tidak membatasi bentuk kurva, menjadi kandidat statistik terbaik pada kelima pengacakan.
  2. AICc dan CV memilih model yang sama pada 5 dari 5 pengacakan, sehingga kedua kriteria cukup konsisten untuk himpunan kandidat ini.
  3. Interaksi dua arah (M3) memperbaiki M1 secara signifikan menurut uji F, tetapi M2 (aditif derajat 3) tetap mengungguli M3 di kelima pengacakan pada AICc, CV, dan data uji. Dalam data ini, kelengkungan lebih penting daripada interaksi.
  4. Random Forest memberikan akurasi prediksi tertinggi, sehingga masih ada struktur yang belum tertangkap oleh model statistik. Untuk tujuan menjelaskan hubungan dan menguji hipotesis, model statistik (terutama GAM) tetap menjadi pilihan utama; untuk prediksi semata, Random Forest layak dipertimbangkan.