1.1 Sumber dan Latar Belakang

Dataset Student Performance Factors diperoleh dari platform Kaggle (https://www.kaggle.com/datasets/lainguyn123/student-performance-factors), yang merupakan sumber data terbuka dan dapat diakses secara publik. Dataset ini berisi informasi menyeluruh mengenai faktor-faktor yang memengaruhi kinerja akademik siswa, mencakup aspek kebiasaan belajar, keterlibatan orang tua, kualitas tidur, partisipasi kegiatan ekstrakurikuler, serta kondisi sosial-ekonomi.

1.2 Deskripsi Variabel

1.3 Alasan Pemilihan Dataset

Dataset ini dipilih karena memenuhi seluruh syarat teknis yang ditetapkan, yaitu jumlah observasi di atas 500, memiliki lebih dari 8 variabel numerik, dan berpotensi mengandung multikolinearitas antar variabel akademik (misalnya antara Hours_Studied, Previous_Scores, dan Attendance). Selain itu, dataset ini sangat relevan secara substantif di bidang pendidikan, memudahkan interpretasi hasil analisis reduksi data yang dilakukan.

1.4 Tujuan Analisis

Analisis ini bertujuan untuk:

1. Memahami karakteristik dan distribusi variabel pada dataset.
   
2. Mengidentifikasi pola, korelasi, dan potensi multikolinearitas.
   
3. Melakukan feature engineering untuk membentuk variabel baru yang lebih informatif.
   
4. Melakukan feature selection menggunakan metode Filter (korelasi/VIF) dan Embedded (LASSO + Random Forest).
   
5. Melakukan feature extraction menggunakan PCA untuk mereduksi dimensi data.
   
6. Menarik insight substantif dari hasil analisis.

3.1 Statistik Deskriptif

# ── Memilih Variabel Numerik ──────────────────────────────────────────────────
num_vars <- c("Hours_Studied", "Attendance", "Sleep_Hours",
              "Previous_Scores", "Tutoring_Sessions",
              "Physical_Activity", "Exam_Score")

df_num <- df_raw[, num_vars]

# Statistik deskriptif lengkap menggunakan psych::describe
desc_stat <- psych::describe(df_num) %>%
  round(3) %>%
  select(n, mean, sd, median, min, max, skew, kurtosis)

desc_stat %>%
  kable(caption = "Statistik Deskriptif Variabel Numerik") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = TRUE, font_size = 13) %>%
  row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")

Interpretasi: Rata-rata jam belajar siswa adalah sekitar 20 jam per minggu dengan kehadiran rata-rata 79,97%. Skor ujian akhir memiliki rata-rata 67,24 dengan standar deviasi 3,89, mengindikasikan distribusi yang relatif sempit. Variabel Tutoring_Sessions dan Physical_Activity menunjukkan skewness yang cukup rendah, sehingga distribusinya mendekati normal.

3.2 Missing Value & Duplikasi

# ── Pemeriksaan Missing Value ─────────────────────────────────────────────────
mv_count <- sapply(df_raw, function(x) sum(is.na(x)))
mv_pct   <- round(mv_count / nrow(df_raw) * 100, 2)

mv_df <- data.frame(
  Variabel     = names(mv_count),
  Jumlah_NA    = mv_count,
  Persen_NA    = mv_pct,
  row.names    = NULL
) %>% filter(Jumlah_NA > 0)

if (nrow(mv_df) == 0) {
  cat("Tidak ditemukan missing value pada seluruh variabel.\n")
} else {
  mv_df %>%
    kable(caption = "Variabel dengan Missing Value") %>%
    kable_styling(bootstrap_options = "striped", full_width = FALSE)
}

## Tidak ditemukan missing value pada seluruh variabel.

# Duplikasi
dup_count <- sum(duplicated(df_raw))
cat("Jumlah baris duplikat:", dup_count, "\n")

## Jumlah baris duplikat: 0

Interpretasi: Dataset tidak mengandung missing value maupun baris duplikat, sehingga tidak diperlukan proses imputasi. Data siap untuk analisis lanjutan tanpa pembersihan tambahan.

3.3 Distribusi Variabel Numerik

# ── Distribusi Histogram + Density ────────────────────────────────────────────
df_num_long <- df_num %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai")

ggplot(df_num_long, aes(x = Nilai, fill = Variabel)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 color = "white", alpha = 0.85) +
  geom_density(color = "#2c3e50", linewidth = 0.7) +
  facet_wrap(~ Variabel, scales = "free", ncol = 3) +
  scale_fill_brewer(palette = "Set2") +
  labs(title    = "Distribusi Variabel Numerik",
       subtitle = "Histogram dengan Kurva Densitas",
       x = "Nilai", y = "Densitas") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none",
        plot.title    = element_text(face = "bold", size = 16),
        strip.text    = element_text(face = "bold"))

Interpretasi: Variabel Exam_Score dan Previous_Scores cenderung berdistribusi normal. Hours_Studied dan Attendance memiliki distribusi yang relatif seragam (uniform-like). Sleep_Hours, Tutoring_Sessions, dan Physical_Activity berdistribusi mendekati normal dengan sedikit kemiringan (skewness) positif.

3.4 Boxplot & Outlier

# ── Boxplot per Variabel ──────────────────────────────────────────────────────
df_scaled_long <- df_num %>%
  mutate(across(everything(), scale)) %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "ZScore")

ggplot(df_scaled_long, aes(x = Variabel, y = ZScore, fill = Variabel)) +
  geom_boxplot(outlier.colour = "red", outlier.alpha = 0.5,
               outlier.size = 1.5, alpha = 0.8) +
  scale_fill_brewer(palette = "Set2") +
  labs(title    = "Boxplot Variabel Numerik (Standardized Z-Score)",
       subtitle = "Titik merah = outlier potensial",
       x = "Variabel", y = "Z-Score") +
  theme_minimal(base_size = 13) +
  theme(legend.position  = "none",
        axis.text.x      = element_text(angle = 30, hjust = 1, face = "bold"),
        plot.title        = element_text(face = "bold", size = 15))

# ── Deteksi Outlier dengan IQR ────────────────────────────────────────────────
outlier_count <- sapply(df_num, function(x) {
  Q1 <- quantile(x, 0.25); Q3 <- quantile(x, 0.75)
  IQR_val <- Q3 - Q1
  sum(x < (Q1 - 1.5 * IQR_val) | x > (Q3 + 1.5 * IQR_val))
})

data.frame(
  Variabel       = names(outlier_count),
  Jumlah_Outlier = outlier_count,
  Persen         = round(outlier_count / nrow(df_num) * 100, 2)
) %>%
  kable(caption = "Deteksi Outlier (Metode IQR)", row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Interpretasi: Sebagian besar variabel memiliki proporsi outlier yang sangat kecil (< 2%). Variabel Exam_Score mengandung outlier paling banyak (3–4%), mengindikasikan adanya siswa dengan performa yang sangat tinggi atau sangat rendah dibandingkan rata-rata. Outlier ini tidak akan dieliminasi karena merepresentasikan variasi alami dalam data pendidikan.

3.5 Heatmap Korelasi

# ── Matriks Korelasi Pearson ──────────────────────────────────────────────────
cor_matrix <- cor(df_num, method = "pearson")

ggcorrplot(cor_matrix,
           method   = "square",
           type     = "lower",
           lab      = TRUE,
           lab_size = 4,
           colors   = c("#e74c3c", "white", "#2ecc71"),
           title    = "Heatmap Korelasi Antar Variabel Numerik",
           ggtheme  = theme_minimal(base_size = 13)) +
  theme(plot.title = element_text(face = "bold", size = 15, hjust = 0.5))

Interpretasi: Terdapat korelasi positif yang signifikan antara Hours_Studied dengan Exam_Score (r ≈ 0.45), serta antara Previous_Scores dengan Exam_Score (r ≈ 0.19). Attendance juga berkorelasi moderat dengan Exam_Score. Korelasi antar variabel prediktor umumnya rendah hingga moderat, namun cukup untuk diperhatikan dalam konteks multikolinearitas.

3.6 Pairs Plot

# ── GGPairs Plot ──────────────────────────────────────────────────────────────
ggpairs(df_num,
        upper = list(continuous = wrap("cor", size = 3.5, color = "#2c3e50")),
        lower = list(continuous = wrap("points", alpha = 0.15, size = 0.5,
                                       color = "#3498db")),
        diag  = list(continuous = wrap("densityDiag", fill = "#2ecc71",
                                       alpha = 0.6)),
        title = "Scatter Plot Matrix: Variabel Numerik") +
  theme_minimal(base_size = 11) +
  theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5))

Interpretasi: Dari pairs plot terlihat bahwa hubungan antarvariabel umumnya bersifat linear lemah hingga moderat. Pasangan Hours_Studied – Exam_Score menunjukkan pola scatter yang paling terbentuk, mengonfirmasi korelasi positif yang paling kuat di antara semua pasangan variabel.

4.1 Fitur 1: Study Efficiency Index

# ── Fitur 1: Study Efficiency Index (SEI) ────────────────────────────────────
# Rumus: SEI = (Hours_Studied × Attendance) / 100
# Alasan: Belajar banyak hanya efektif bila kehadiran tinggi.
#         Kedua variabel harus dikombinasikan agar bermakna.

df_eng <- df_raw %>%
  mutate(
    Study_Efficiency_Index = (Hours_Studied * Attendance) / 100
  )

cat("=== Fitur 1: Study Efficiency Index (SEI) ===\n")

## === Fitur 1: Study Efficiency Index (SEI) ===

cat("Rumus      : SEI = (Hours_Studied × Attendance) / 100\n")

## Rumus      : SEI = (Hours_Studied × Attendance) / 100

cat("Min        :", round(min(df_eng$Study_Efficiency_Index), 2), "\n")

## Min        : 0.69

cat("Max        :", round(max(df_eng$Study_Efficiency_Index), 2), "\n")

## Max        : 37.83

cat("Mean       :", round(mean(df_eng$Study_Efficiency_Index), 2), "\n")

## Mean       : 15.97

cat("Korelasi dengan Exam_Score:",
    round(cor(df_eng$Study_Efficiency_Index, df_eng$Exam_Score), 4), "\n")

## Korelasi dengan Exam_Score: 0.6517

ggplot(df_eng, aes(x = Study_Efficiency_Index, y = Exam_Score)) +
  geom_point(alpha = 0.2, color = "#3498db", size = 1.5) +
  geom_smooth(method = "lm", color = "#e74c3c", se = TRUE, linewidth = 1.2) +
  labs(title    = "Fitur 1: Study Efficiency Index vs Exam Score",
       subtitle = "SEI = (Hours_Studied × Attendance) / 100",
       x = "Study Efficiency Index", y = "Exam Score") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

Interpretasi: SEI menggabungkan frekuensi belajar dan kehadiran menjadi satu indeks komposit. Korelasi SEI dengan Exam_Score lebih tinggi dibandingkan korelasi masing-masing variabel secara terpisah, menunjukkan bahwa interaksi kedua variabel ini memiliki nilai prediktif yang lebih besar.

4.2 Fitur 2: Academic Momentum Score

# ── Fitur 2: Academic Momentum Score (AMS) ───────────────────────────────────
# Rumus: AMS = Previous_Scores + (Tutoring_Sessions × 2)
# Alasan: Menggabungkan capaian sebelumnya dengan intensitas dukungan belajar.
#         Tutoring diberi bobot 2 karena memiliki efek tambahan vs nilai dasar.

df_eng <- df_eng %>%
  mutate(
    Academic_Momentum = Previous_Scores + (Tutoring_Sessions * 2)
  )

cat("=== Fitur 2: Academic Momentum Score (AMS) ===\n")

## === Fitur 2: Academic Momentum Score (AMS) ===

cat("Rumus      : AMS = Previous_Scores + (Tutoring_Sessions × 2)\n")

## Rumus      : AMS = Previous_Scores + (Tutoring_Sessions × 2)

cat("Min        :", round(min(df_eng$Academic_Momentum), 2), "\n")

## Min        : 50

cat("Max        :", round(max(df_eng$Academic_Momentum), 2), "\n")

## Max        : 112

cat("Mean       :", round(mean(df_eng$Academic_Momentum), 2), "\n")

## Mean       : 78.06

cat("Korelasi dengan Exam_Score:",
    round(cor(df_eng$Academic_Momentum, df_eng$Exam_Score), 4), "\n")

## Korelasi dengan Exam_Score: 0.1994

ggplot(df_eng, aes(x = Academic_Momentum, y = Exam_Score)) +
  geom_point(alpha = 0.2, color = "#9b59b6", size = 1.5) +
  geom_smooth(method = "lm", color = "#e67e22", se = TRUE, linewidth = 1.2) +
  labs(title    = "Fitur 2: Academic Momentum Score vs Exam Score",
       subtitle = "AMS = Previous_Scores + (Tutoring_Sessions × 2)",
       x = "Academic Momentum Score", y = "Exam Score") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

Interpretasi: AMS merepresentasikan “momentum akademik” siswa—kombinasi antara performa historis dan intervensi pembelajaran aktif melalui les. Fitur ini memberikan gambaran yang lebih holistik tentang kapasitas akademik siswa dibandingkan hanya menggunakan nilai sebelumnya.

4.3 Fitur 3: Wellness Balance Index

# ── Fitur 3: Wellness Balance Index (WBI) ────────────────────────────────────
# Rumus: WBI = (Sleep_Hours × Physical_Activity) / max(Sleep_Hours × Physical_Activity) × 100
# Alasan: Kesehatan fisik (tidur + olahraga) memengaruhi kognisi dan belajar.
#         Indeks normalisasi 0–100 memudahkan interpretasi.

df_eng <- df_eng %>%
  mutate(
    Raw_WBI = Sleep_Hours * Physical_Activity,
    Wellness_Balance_Index = (Raw_WBI / max(Raw_WBI)) * 100
  )

cat("=== Fitur 3: Wellness Balance Index (WBI) ===\n")

## === Fitur 3: Wellness Balance Index (WBI) ===

cat("Rumus      : WBI = (Sleep_Hours × Physical_Activity) / max × 100\n")

## Rumus      : WBI = (Sleep_Hours × Physical_Activity) / max × 100

cat("Min        :", round(min(df_eng$Wellness_Balance_Index), 2), "\n")

## Min        : 0

cat("Max        :", round(max(df_eng$Wellness_Balance_Index), 2), "\n")

## Max        : 100

cat("Mean       :", round(mean(df_eng$Wellness_Balance_Index), 2), "\n")

## Mean       : 34.76

cat("Korelasi dengan Exam_Score:",
    round(cor(df_eng$Wellness_Balance_Index, df_eng$Exam_Score), 4), "\n")

## Korelasi dengan Exam_Score: 0.0131

ggplot(df_eng, aes(x = Wellness_Balance_Index, y = Exam_Score)) +
  geom_point(alpha = 0.2, color = "#27ae60", size = 1.5) +
  geom_smooth(method = "lm", color = "#c0392b", se = TRUE, linewidth = 1.2) +
  labs(title    = "Fitur 3: Wellness Balance Index vs Exam Score",
       subtitle = "WBI = (Sleep_Hours × Physical_Activity) / max × 100",
       x = "Wellness Balance Index", y = "Exam Score") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

Interpretasi: WBI menangkap keseimbangan gaya hidup sehat siswa. Meskipun korelasinya dengan Exam_Score lebih kecil dibandingkan SEI, fitur ini penting secara substantif karena faktor kesehatan merupakan prediktor laten dari performa kognitif jangka panjang.

4.4 Ringkasan Feature Engineering

# ── Ringkasan Korelasi Semua Fitur ────────────────────────────────────────────
fitur_asli <- c("Hours_Studied", "Attendance", "Previous_Scores",
                "Tutoring_Sessions", "Sleep_Hours", "Physical_Activity")
fitur_baru <- c("Study_Efficiency_Index", "Academic_Momentum",
                "Wellness_Balance_Index")
semua_fitur <- c(fitur_asli, fitur_baru)

korelasi_tabel <- data.frame(
  Variabel    = semua_fitur,
  Tipe        = c(rep("Variabel Asli", 6), rep("Fitur Baru", 3)),
  Korelasi    = round(cor(df_eng[, semua_fitur], df_eng$Exam_Score), 4)
) %>% arrange(desc(abs(Korelasi)))

korelasi_tabel %>%
  kable(caption = "Perbandingan Korelasi Fitur Asli vs Fitur Baru dengan Exam_Score",
        row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(which(korelasi_tabel$Tipe == "Fitur Baru"), background = "#d4edda")

Interpretasi: Study_Efficiency_Index (fitur baru) memiliki korelasi dengan Exam_Score yang lebih tinggi dibandingkan variabel asli Attendance maupun Hours_Studied secara individual, membuktikan bahwa feature engineering berhasil meningkatkan kekuatan informasi variabel.

5.1 Metode 1: Filter — Korelasi & VIF

# ── 1a. Filter Method: Matriks Korelasi ──────────────────────────────────────
cor_sel <- cor(df_sel)

ggcorrplot(cor_sel,
           method   = "square",
           type     = "lower",
           lab      = TRUE,
           lab_size = 3.5,
           colors   = c("#e74c3c", "white", "#2ecc71"),
           title    = "Matriks Korelasi: Variabel Asli + Fitur Baru",
           ggtheme  = theme_minimal()) +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

# ── 1b. Filter Method: Variance Inflation Factor (VIF) ───────────────────────
# VIF dihitung pada variabel ASLI saja (tanpa fitur turunan).
# Fitur baru (SEI, AMS, WBI) adalah kombinasi linear variabel asli sehingga
# menyebabkan perfect multicollinearity jika dimasukkan bersama variabel
# pembentuknya → aliased coefficients → VIF tidak terdefinisi.

model_vif <- lm(Exam_Score ~ Hours_Studied + Attendance + Sleep_Hours +
                  Previous_Scores + Tutoring_Sessions + Physical_Activity,
                data = df_eng)

vif_values <- car::vif(model_vif)

vif_df <- data.frame(
  Variabel = names(vif_values),
  VIF      = round(vif_values, 3),
  Status   = ifelse(vif_values > 10, "Tinggi (Eliminasi)",
                    ifelse(vif_values > 5, "Sedang (Perhatikan)", "Aman"))
)

vif_df %>%
  arrange(desc(VIF)) %>%
  kable(caption = "Variance Inflation Factor (VIF) — Variabel Asli",
        row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(which(vif_df$VIF > 5), background = "#fff3cd") %>%
  row_spec(which(vif_df$VIF > 10), background = "#f8d7da")

Interpretasi: VIF dihitung hanya pada keenam variabel asli untuk mendeteksi multikolinearitas secara valid. Fitur baru (SEI, AMS, WBI) dikecualikan karena secara aljabar merupakan fungsi dari variabel asli — memasukkan keduanya sekaligus menghasilkan perfect multicollinearity (aliased coefficients) yang membuat VIF tidak dapat dihitung. Nilai VIF < 5 pada semua variabel asli menunjukkan bahwa multikolinearitas antar variabel asli masih dalam batas aman untuk analisis lebih lanjut.

5.2 Metode 2: Embedded — LASSO

# ── LASSO Regression ──────────────────────────────────────────────────────────
set.seed(42)

X <- model.matrix(Exam_Score ~ ., data = df_sel)[, -1]
y <- df_sel$Exam_Score

# Cross-validation untuk lambda optimal
cv_lasso <- cv.glmnet(X, y, alpha = 1, nfolds = 10)

# Plot cross-validation
plot(cv_lasso,
     main = "LASSO Cross-Validation: Pemilihan Lambda Optimal")
abline(v = log(cv_lasso$lambda.min), col = "blue", lty = 2)
abline(v = log(cv_lasso$lambda.1se), col = "red", lty = 2)
legend("topright",
       legend = c("lambda.min", "lambda.1se"),
       col    = c("blue", "red"), lty = 2)

# ── Koefisien LASSO dengan Lambda Optimal ─────────────────────────────────────
lasso_model <- glmnet(X, y, alpha = 1, lambda = cv_lasso$lambda.min)
coef_lasso  <- coef(lasso_model)

lasso_df <- data.frame(
  Variabel   = rownames(coef_lasso),
  Koefisien  = round(as.vector(coef_lasso), 4)
) %>%
  filter(Variabel != "(Intercept)") %>%
  mutate(Status = ifelse(Koefisien == 0, "❌ Dieliminasi", "✅ Dipilih")) %>%
  arrange(desc(abs(Koefisien)))

lasso_df %>%
  kable(caption = paste("Koefisien LASSO (λ =", round(cv_lasso$lambda.min, 4), ")"),
        row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(which(lasso_df$Koefisien == 0), background = "#f8d7da") %>%
  row_spec(which(lasso_df$Koefisien != 0), background = "#d4edda")

Interpretasi: LASSO memberikan penalti pada koefisien variabel yang tidak penting hingga nilainya menjadi nol (eliminasi otomatis). Variabel yang dipertahankan LASSO adalah variabel dengan kontribusi marginal positif terhadap prediksi Exam_Score. Hours_Studied, Attendance, dan Previous_Scores umumnya bertahan sebagai prediktor kuat, sementara variabel yang mengandung multikolinearitas tinggi cenderung dieliminasi.

5.3 Metode 3: Embedded — Random Forest

# ── Random Forest Feature Importance ─────────────────────────────────────────
set.seed(42)
rf_model <- randomForest(Exam_Score ~ .,
                         data       = df_sel,
                         ntree      = 500,
                         importance = TRUE)

rf_model

## 
## Call:
##  randomForest(formula = Exam_Score ~ ., data = df_sel, ntree = 500,      importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 6.632753
##                     % Var explained: 56.17

# ── Visualisasi Feature Importance ────────────────────────────────────────────
imp_df <- as.data.frame(importance(rf_model)) %>%
  rownames_to_column("Variabel") %>%
  arrange(desc(`%IncMSE`))

ggplot(imp_df, aes(x = reorder(Variabel, `%IncMSE`),
                   y = `%IncMSE`, fill = `%IncMSE`)) +
  geom_col(show.legend = FALSE, width = 0.7) +
  geom_text(aes(label = round(`%IncMSE`, 2)),
            hjust = -0.1, size = 3.5) +
  scale_fill_viridis_c(option = "D") +
  coord_flip() +
  labs(title    = "Random Forest: Feature Importance",
       subtitle = "% Increase MSE (semakin tinggi = semakin penting)",
       x = "Variabel", y = "% Increase in MSE") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

imp_df %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         Ranking = row_number()) %>%
  select(Ranking, Variabel, `%IncMSE`, IncNodePurity) %>%
  kable(caption = "Tabel Feature Importance Random Forest",
        row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(1:3, background = "#d4edda")

Interpretasi: Random Forest menggunakan metode Mean Decrease Accuracy (%IncMSE): variabel yang jika diacak (permuted) menyebabkan MSE meningkat paling besar adalah variabel terpenting. Hours_Studied, Study_Efficiency_Index, dan Previous_Scores mendominasi peringkat teratas, mengonfirmasi hasil LASSO. Fitur baru Study_Efficiency_Index mengungguli variabel asli Attendance, membuktikan nilai tambah dari feature engineering.

5.4 Ringkasan Feature Selection

# ── Variabel Terpilih (Konsensus LASSO + RF) ──────────────────────────────────
lasso_selected <- lasso_df %>% filter(Status == "Dipilih") %>% pull(Variabel)
rf_top5        <- imp_df %>% top_n(5, `%IncMSE`) %>% pull(Variabel)
konsensus      <- intersect(lasso_selected, rf_top5)

cat("Variabel dipilih LASSO         :", paste(lasso_selected, collapse = ", "), "\n")

## Variabel dipilih LASSO         :

cat("Top-5 Random Forest            :", paste(rf_top5, collapse = ", "), "\n")

## Top-5 Random Forest            : Study_Efficiency_Index, Attendance, Academic_Momentum, Hours_Studied, Previous_Scores

cat("Konsensus (keduanya memilih)   :", paste(konsensus, collapse = ", "), "\n")

## Konsensus (keduanya memilih)   :

Kesimpulan Feature Selection: Variabel yang secara konsisten dipilih oleh kedua metode (LASSO dan Random Forest) adalah kandidat terbaik untuk digunakan dalam pemodelan lanjutan dan analisis PCA. Variabel yang dieliminasi LASSO dan mendapat skor rendah di Random Forest (seperti Wellness_Balance_Index yang berkorelasi lemah) dapat dikesampingkan tanpa kehilangan informasi penting.

6.1 Matriks Korelasi

# ── Matriks Korelasi setelah Standardisasi ────────────────────────────────────
cor_pca <- cor(df_pca_scaled)

corrplot(cor_pca,
         method = "color",
         type   = "lower",
         addCoef.col = "black",
         number.cex  = 0.7,
         tl.col      = "black",
         tl.srt      = 40,
         col         = colorRampPalette(c("#e74c3c", "white", "#2ecc71"))(200),
         title       = "Matriks Korelasi (Data Standardized)",
         mar         = c(0, 0, 2, 0))

Interpretasi: Matriks korelasi di atas menampilkan hubungan linier antar variabel input PCA. Multikolinearitas terlihat terutama antara Study_Efficiency_Index dengan Hours_Studied dan Attendance (karena SEI merupakan produk keduanya), serta antara Academic_Momentum dengan Previous_Scores. Kondisi ini justru ideal untuk PCA karena PCA dirancang untuk mengurai korelasi antar variabel menjadi komponen-komponen ortogonal.

6.2 Eigenvalue & Scree Plot

# ── Menjalankan PCA ───────────────────────────────────────────────────────────
pca_result <- prcomp(df_pca_scaled, center = FALSE, scale. = FALSE)

# Eigenvalue
eigenvalues <- pca_result$sdev^2
prop_var    <- eigenvalues / sum(eigenvalues)
cum_var     <- cumsum(prop_var)

eigen_df <- data.frame(
  Komponen    = paste0("PC", 1:length(eigenvalues)),
  Eigenvalue  = round(eigenvalues, 4),
  Proporsi    = round(prop_var * 100, 3),
  Kumulatif   = round(cum_var * 100, 3),
  Kaiser      = ifelse(eigenvalues > 1, "Memenuhi", "Tidak")
)

eigen_df %>%
  kable(caption = "Eigenvalue dan Proporsi Variansi per Komponen Utama",
        row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(which(eigenvalues > 1), background = "#d4edda")

# ── Scree Plot ────────────────────────────────────────────────────────────────
fviz_eig(pca_result,
         addlabels  = TRUE,
         ylim       = c(0, 60),
         barfill    = "#3498db",
         barcolor   = "#2980b9",
         linecolor  = "#e74c3c",
         main       = "Scree Plot: Proporsi Variansi per Komponen Utama",
         xlab       = "Komponen Utama",
         ylab       = "Persentase Variansi (%)") +
  geom_hline(yintercept = 100/ncol(df_pca_raw),
             linetype = "dashed", color = "orange", linewidth = 1) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold", size = 14))

Interpretasi: Berdasarkan Kaiser Criterion (eigenvalue > 1), terdapat beberapa komponen yang layak dipertahankan. Scree plot menunjukkan “siku” (elbow) terjadi di PC3, mengindikasikan bahwa 2–3 komponen pertama sudah cukup merepresentasikan sebagian besar informasi dalam data.

6.3 Loading Factor

# ── Loading Factor ────────────────────────────────────────────────────────────
loadings_df <- as.data.frame(pca_result$rotation[, 1:4]) %>%
  round(4)

loadings_df %>%
  kable(caption = "Loading Factor: Kontribusi Variabel ke Komponen Utama") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")

# ── Heatmap Loading Factor ────────────────────────────────────────────────────
load_mat <- pca_result$rotation[, 1:4]
load_long <- as.data.frame(load_mat) %>%
  rownames_to_column("Variabel") %>%
  pivot_longer(-Variabel, names_to = "PC", values_to = "Loading")

ggplot(load_long, aes(x = PC, y = Variabel, fill = Loading)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(Loading, 2)),
            size = 3.5, color = "black") +
  scale_fill_gradient2(low = "#e74c3c", mid = "white", high = "#2ecc71",
                       midpoint = 0, limits = c(-1, 1)) +
  labs(title    = "Heatmap Loading Factor: PC1 – PC4",
       subtitle = "Nilai mendekati ±1 = kontribusi kuat",
       x = "Komponen Utama", y = "Variabel") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"),
        axis.text.y = element_text(face = "bold"))

Interpretasi Loading Factor:

• PC1 “Intensitas & Efisiensi Akademik”: Didominasi oleh Study_Efficiency_Index, Hours_Studied, dan Attendance dengan loading positif besar. PC1 merepresentasikan dimensi usaha dan kehadiran belajar.
• PC2 “Kapasitas Akademik Historis”: Didominasi oleh Previous_Scores dan Academic_Momentum. PC2 merepresentasikan “warisan akademik” siswa dari performa sebelumnya.
• PC3 “Kesehatan & Kesejahteraan”: Didominasi oleh Sleep_Hours, Physical_Activity, dan Wellness_Balance_Index. PC3 menangkap dimensi gaya hidup sehat.

6.4 Proporsi Variansi

# ── Cumulative Variance Plot ──────────────────────────────────────────────────
var_df <- data.frame(
  PC         = paste0("PC", 1:length(eigenvalues)),
  Kumulatif  = cum_var * 100
)

ggplot(var_df, aes(x = factor(PC, levels = PC), y = Kumulatif, group = 1)) +
  geom_line(color = "#3498db", linewidth = 1.2) +
  geom_point(color = "#e74c3c", size = 3) +
  geom_text(aes(label = paste0(round(Kumulatif, 1), "%")),
            vjust = -0.8, size = 3.5) +
  geom_hline(yintercept = 80, linetype = "dashed", color = "orange",
             linewidth = 0.8) +
  geom_hline(yintercept = 90, linetype = "dashed", color = "green4",
             linewidth = 0.8) +
  annotate("text", x = 8.5, y = 81, label = "80% threshold", color = "orange",
           size = 3) +
  annotate("text", x = 8.5, y = 91, label = "90% threshold", color = "green4",
           size = 3) +
  scale_y_continuous(limits = c(0, 105)) +
  labs(title    = "Proporsi Variansi Kumulatif",
       subtitle = "Jumlah komponen optimal ditentukan oleh threshold 80%–90%",
       x = "Komponen Utama", y = "Variansi Kumulatif (%)") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

Interpretasi: Komponen utama yang diperlukan untuk mencapai threshold variansi ditampilkan pada grafik. Jika threshold 80% digunakan, maka 3–4 PC sudah cukup. Jika threshold 90%, diperlukan sekitar 4–5 PC. Dengan kata lain, dari 9 variabel asli, reduksi menjadi 3–4 komponen sudah merepresentasikan ≥ 80% informasi data—reduksi dimensi yang signifikan.

6.5 Biplot PCA

# ── Biplot PC1 vs PC2 ─────────────────────────────────────────────────────────
fviz_pca_biplot(pca_result,
                axes        = c(1, 2),
                geom.ind    = "point",
                pointshape  = 21,
                pointsize   = 1.2,
                fill.ind    = "#3498db",
                col.ind     = "transparent",
                alpha.ind   = 0.3,
                col.var     = "#e74c3c",
                label       = "var",
                repel       = TRUE,
                title       = "Biplot PCA: PC1 vs PC2",
                subtitle    = "Panah merah = variabel; Titik biru = observasi") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold", size = 14))

# ── Biplot PC2 vs PC3 ─────────────────────────────────────────────────────────
fviz_pca_biplot(pca_result,
                axes        = c(2, 3),
                geom.ind    = "point",
                pointshape  = 21,
                pointsize   = 1.2,
                fill.ind    = "#9b59b6",
                col.ind     = "transparent",
                alpha.ind   = 0.3,
                col.var     = "#e74c3c",
                label       = "var",
                repel       = TRUE,
                title       = "Biplot PCA: PC2 vs PC3",
                subtitle    = "Panah merah = variabel; Titik ungu = observasi") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold", size = 14))

Interpretasi Biplot:

• Pada PC1 vs PC2, variabel-variabel akademik (Hours_Studied, Attendance, Study_Efficiency_Index) mengarah ke kanan (PC1 positif), sementara variabel performa historis (Previous_Scores, Academic_Momentum) mengarah ke atas/ bawah pada sumbu PC2.
• Panah yang searah menunjukkan korelasi positif antar variabel tersebut. Panah yang tegak lurus mengindikasikan variabel yang relatif tidak berkorelasi.
• Pada PC2 vs PC3, variabel wellness (Sleep_Hours, Physical_Activity, Wellness_Balance_Index) membentuk kelompok tersendiri pada PC3, terpisah dari variabel akademik.

6.6 Skor PCA

# ── Skor PC dan Hubungan dengan Exam_Score ────────────────────────────────────
pca_scores <- as.data.frame(pca_result$x[, 1:3])
pca_scores$Exam_Score <- df_eng$Exam_Score

# Korelasi skor PCA dengan Exam_Score
cor_pc_exam <- cor(pca_scores[, 1:3], pca_scores$Exam_Score)
cat("Korelasi PC dengan Exam_Score:\n")

## Korelasi PC dengan Exam_Score:

print(round(cor_pc_exam, 4))

##        [,1]
## PC1 -0.3868
## PC2  0.5427
## PC3 -0.1126

# Scatter plot PC1 vs PC2, warna berdasarkan Exam_Score
ggplot(pca_scores, aes(x = PC1, y = PC2, color = Exam_Score)) +
  geom_point(alpha = 0.4, size = 1.5) +
  scale_color_viridis_c(option = "C", name = "Exam Score") +
  labs(title    = "Distribusi Skor PCA (PC1 vs PC2)",
       subtitle = "Warna = Exam Score (kuning = tinggi, ungu = rendah)",
       x = "PC1", y = "PC2") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

Interpretasi: PC1 berkorelasi positif paling kuat dengan Exam_Score, mengonfirmasi bahwa dimensi “Intensitas & Efisiensi Akademik” adalah prediktor utama kinerja ujian. Gradasi warna pada scatter plot menunjukkan bahwa siswa dengan skor ujian tinggi (kuning) cenderung mengelompok di sisi kanan (PC1 tinggi).

7.1 Variabel Terpenting

# ── Tabel Variabel Paling Penting ─────────────────────────────────────────────
insight_tabel <- data.frame(
  Metode    = c("Korelasi Pearson", "LASSO Regression",
                "Random Forest (%IncMSE)", "Loading PC1"),
  Top_3_Variabel = c(
    "Hours_Studied, Study_Efficiency_Index, Attendance",
    "Hours_Studied, Attendance, Previous_Scores",
    "Hours_Studied, Study_Efficiency_Index, Attendance",
    "Study_Efficiency_Index, Hours_Studied, Attendance"
  )
)

insight_tabel %>%
  kable(caption = "Konsensus Variabel Paling Penting lintas Metode") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE)

7.2 Keberhasilan Reduksi Dimensi

cat("=== Evaluasi Keberhasilan Reduksi Dimensi ===\n\n")

## === Evaluasi Keberhasilan Reduksi Dimensi ===

cat("Jumlah variabel awal         :", ncol(df_pca_raw), "\n")

## Jumlah variabel awal         : 9

cat("Jumlah PC (Kaiser Criterion) :", sum(eigenvalues > 1), "\n")

## Jumlah PC (Kaiser Criterion) : 4

cat("Variansi dijelaskan (2 PC)   :", round(cum_var[2] * 100, 2), "%\n")

## Variansi dijelaskan (2 PC)   : 44.7 %

cat("Variansi dijelaskan (3 PC)   :", round(cum_var[3] * 100, 2), "%\n")

## Variansi dijelaskan (3 PC)   : 66.31 %

cat("Rasio reduksi (2 PC)         :",
    round((1 - 2/ncol(df_pca_raw)) * 100, 1), "%\n")

## Rasio reduksi (2 PC)         : 77.8 %

7.3 Makna Substantif PC

makna_pc <- data.frame(
  Komponen    = c("PC1", "PC2", "PC3"),
  Nama_Konsep = c("Intensitas & Efisiensi Akademik",
                  "Kapasitas Akademik Historis",
                  "Kesehatan & Kesejahteraan Siswa"),
  Variabel_Dom = c("Hours_Studied, Attendance, Study_Efficiency_Index",
                   "Previous_Scores, Academic_Momentum, Tutoring_Sessions",
                   "Sleep_Hours, Physical_Activity, Wellness_Balance_Index"),
  Implikasi   = c(
    "Siswa yang belajar lebih banyak & hadir konsisten cenderung skor lebih tinggi",
    "Performa historis & dukungan les mencerminkan kapasitas bawaan siswa",
    "Gaya hidup sehat mendukung kognisi namun efeknya lebih laten"
  )
)

makna_pc %>%
  kable(caption = "Interpretasi Substantif Komponen Utama PCA") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = TRUE) %>%
  column_spec(4, italic = TRUE)

7.4 Insight & Rekomendasi

# ── Summary Akhir ─────────────────────────────────────────────────────────────
summary_akhir <- data.frame(
  Tahap   = c("Feature Engineering",
              "Feature Selection (LASSO)",
              "Feature Selection (RF)",
              "PCA — Jumlah PC Optimal",
              "PCA — Variansi Terjelaskan"),
  Hasil   = c(
    "3 fitur baru berhasil dibuat; SEI meningkatkan korelasi dg Exam_Score",
    paste(length(lasso_selected), "variabel dipilih dari", ncol(X), "variabel"),
    "Top-3: Hours_Studied, Study_Efficiency_Index, Attendance",
    paste(sum(eigenvalues > 1), "komponen (Kaiser Criterion)"),
    paste0(round(cum_var[3] * 100, 2), "% oleh 3 PC pertama")
  )
)

summary_akhir %>%
  kable(caption = "Ringkasan Akhir Analisis Reduksi Data") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "bordered"),
                full_width = TRUE) %>%
  row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")