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

Variabel PE berperan sebagai variabel respons (Y), sedangkan AT, V, AP, dan RH berperan sebagai prediktor (X). 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.

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.

1. Hasil Eksplorasi dan Pemeriksaan Data

Tahap eksplorasi dan pemeriksaan data dilakukan untuk mengetahui struktur, karakteristik, serta kondisi awal data sebelum proses pemodelan. Analisis dilakukan secara berurutan meliputi pembacaan dan validasi data, pemeriksaan kesamaan antar-sheet, statistik deskriptif, pemeriksaan missing value dan duplikasi, identifikasi outlier, analisis hubungan antarvariabel, pemeriksaan multikolinearitas, serta pemeriksaan residual pada model regresi awal.

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

sheets <- readxl::excel_sheets(fname)

data_sheets <- lapply(sheets, function(sh) {
  d <- readxl::read_excel(fname, sheet = sh)
  colnames(d) <- trimws(as.character(colnames(d)))
  as.data.frame(d)
})

names(data_sheets) <- sheets
sheet_names <- names(data_sheets)

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

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)

df <- data_sheets[[1]]

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

# 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

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

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)

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

2. Hasil Pemodelan

Setelah seluruh model dijalankan, hasil pemodelan diperoleh untuk setiap sheet dengan menggunakan tiga metode, yaitu Polynomial Regression, Elastic Net, dan Generalized Additive Model (GAM). Proses pemodelan dilakukan menggunakan pembagian data dan skema 10-fold Cross-Validation yang sama sehingga hasil antar-metode dapat dibandingkan.

Hasil pemodelan menunjukkan bahwa ketiga metode menghasilkan nilai kesalahan prediksi yang berbeda. Nilai CV_RMSE Polynomial berada pada kisaran sekitar 4,23–4,28, sedangkan Elastic Net berada pada kisaran sekitar 4,57–4,58. GAM menghasilkan CV_RMSE sekitar 4,12–4,16. Berdasarkan nilai tersebut, GAM menghasilkan nilai CV_RMSE yang lebih rendah dibandingkan Polynomial Regression dan Elastic Net pada seluruh sheet.

Untuk melihat konfigurasi dan nilai kriteria masing-masing metode secara lebih lengkap, hasil pemodelan selanjutnya disusun berdasarkan AICc dan CV_RMSE. Hasil tersebut digunakan pada tahap perbandingan dan pemilihan model pada subbab berikutnya.

Syntax yang digunakan:

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

3. Perbandingan dan Pemilihan Model

Perbandingan model dilakukan berdasarkan nilai AICc dan CV_RMSE yang diperoleh dari proses pemodelan. AICc digunakan untuk mempertimbangkan kesesuaian model dan kompleksitasnya, sedangkan CV_RMSE digunakan untuk mengukur kesalahan prediksi melalui 10-fold Cross-Validation.

Hasil pemodelan menunjukkan bahwa GAM dipilih berdasarkan AICc maupun CV_RMSE pada seluruh sheet. Polynomial Regression memiliki nilai CV_RMSE yang lebih tinggi dibandingkan GAM, sedangkan Elastic Net menghasilkan nilai CV_RMSE yang paling tinggi. Dengan demikian, kedua kriteria pemilihan model memberikan hasil yang konsisten dalam menentukan model yang digunakan.

Kesesuaian peringkat model berdasarkan AICc dan CV_RMSE selanjutnya diuji menggunakan korelasi peringkat Spearman. Hasil pengujian menunjukkan koefisien Spearman sebesar 1 pada seluruh sheet dengan p-value sebesar 0,3333. Nilai koefisien tersebut menunjukkan bahwa urutan ketiga model berdasarkan AICc dan CV_RMSE sama pada setiap sheet. Namun, p-value lebih besar dari 0,05 sehingga hubungan tersebut tidak menunjukkan signifikansi statistik pada taraf 5%.

Secara keseluruhan, hasil perbandingan menunjukkan bahwa AICc dan 10-fold Cross-Validation memberikan hasil pemilihan model yang konsisten pada seluruh sheet. Hasil tersebut selanjutnya digunakan sebagai dasar dalam evaluasi model pada data uji.

Syntax yang digunakan:

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

4. Evaluasi Model pada Data Uji

Evaluasi model pada data uji dilakukan untuk mengetahui kemampuan model dalam menghasilkan prediksi pada observasi yang tidak digunakan dalam proses pelatihan. Kinerja model dievaluasi menggunakan RMSE dan R².

Hasil evaluasi menunjukkan bahwa model yang diperoleh mampu memberikan nilai kesalahan prediksi yang relatif rendah pada data uji. Nilai RMSE berada pada kisaran sekitar 4,0, sedangkan nilai R² berada pada kisaran sekitar 0,94. Hasil tersebut menunjukkan bahwa model mampu menjelaskan sebagian besar variasi PE pada data uji.

Evaluasi data uji juga digunakan untuk melihat kesesuaian hasil dengan pemilihan model berdasarkan AICc dan CV_RMSE. Hasil ini menjadi dasar untuk menilai kinerja model sebelum dilakukan analisis perbedaan kinerja antar-sheet.

Syntax yang digunakan:

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

5. Perbedaan Kinerja Antar Randomisasi

Perbedaan kinerja antar-randomisasi dianalisis menggunakan ANOVA satu arah dan uji Kruskal-Wallis berdasarkan nilai CV-fold RMSE. Kedua pengujian digunakan untuk mengetahui apakah terdapat perbedaan kinerja model antar-sheet.

Hasil ANOVA menunjukkan p-value sebesar 0,9961 untuk Polynomial Regression, 0,99995 untuk Elastic Net, dan 0,9967 untuk GAM. Sementara itu, uji Kruskal-Wallis menghasilkan p-value sebesar 0,8920, 0,9328, dan 0,9067 untuk ketiga metode secara berturut-turut.

Seluruh p-value lebih besar dari 0,05. Dengan demikian, tidak terdapat bukti statistik yang menunjukkan adanya perbedaan kinerja antar-randomisasi berdasarkan nilai CV-fold RMSE, baik menurut ANOVA maupun Kruskal-Wallis.

Hasil ini diperkuat oleh visualisasi boxplot CV-fold RMSE yang menunjukkan pola sebaran yang relatif serupa antar-sheet. Dengan demikian, variasi sheet tidak menghasilkan perbedaan kinerja model yang berarti pada proses Cross-Validation.

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

6. Pembahasan

Hasil eksplorasi menunjukkan bahwa data memiliki hubungan yang kuat antara beberapa variabel prediktor dengan PE, terutama AT dan V. Selain itu, ditemukan adanya indikasi multikolinearitas serta penyimpangan normalitas dan heteroskedastisitas pada residual model linear awal. Kondisi tersebut menunjukkan bahwa hubungan antara prediktor dan PE tidak sepenuhnya sederhana dan linear, sehingga penggunaan metode yang mampu menangkap pola nonlinier menjadi relevan.

Hasil pemodelan menunjukkan bahwa GAM menghasilkan nilai CV_RMSE yang lebih rendah dibandingkan Polynomial Regression dan Elastic Net. Hasil tersebut menunjukkan bahwa pendekatan fungsi smooth pada GAM mampu menangkap pola hubungan nonlinier antara prediktor dan PE dengan lebih baik pada data penelitian. Polynomial Regression juga mampu memodelkan hubungan nonlinier, tetapi bentuk hubungan dibatasi oleh derajat polynomial yang digunakan. Sementara itu, Elastic Net lebih berfokus pada pengendalian kompleksitas dan multikolinearitas melalui penalti sehingga fleksibilitasnya dalam menangkap pola nonlinier lebih terbatas.

Konsistensi hasil pemilihan model berdasarkan AICc dan 10-fold Cross-Validation menunjukkan bahwa kedua kriteria memberikan urutan model yang sama pada setiap randomisasi. Hasil korelasi Spearman yang bernilai 1 juga memperkuat kesamaan peringkat tersebut. Hal ini menunjukkan bahwa keputusan pemilihan model relatif konsisten terhadap kriteria evaluasi yang digunakan.

Hasil pengujian antar-randomisasi menunjukkan tidak terdapat perbedaan kinerja yang signifikan berdasarkan CV-fold RMSE. Kondisi ini dapat dikaitkan dengan kesamaan karakteristik data yang digunakan pada setiap randomisasi. Dengan demikian, perubahan randomisasi tidak menghasilkan perubahan kinerja model yang berarti.

Secara keseluruhan, hasil penelitian menunjukkan bahwa GAM memberikan kemampuan pemodelan yang paling sesuai untuk menangkap hubungan antara AT, V, AP, RH, dan PE pada data penelitian. Hasil evaluasi pada data uji yang menunjukkan nilai RMSE sekitar 4 dan R² sekitar 0,94 juga mendukung kemampuan model dalam menghasilkan prediksi yang baik pada observasi yang tidak digunakan dalam proses pelatihan.