Regresi Lasso dan Mixture Lasso untuk Heating Load dan Cooling Load Berdistribusi Bimodal pada Data Energy Efficiency

Bab I Pendahuluan

1.1 Latar Belakang

Sektor bangunan merupakan salah satu konsumen energi terbesar di dunia, dan sistem pemanas, ventilasi, serta pendingin udara (heating, ventilation, and air conditioning/HVAC) menyerap porsi terbesar dari energi tersebut. Perancangan bangunan yang hemat energi karenanya memerlukan estimasi beban pemanasan (heating load, HL) dan beban pendinginan (cooling load, CL) yang andal sejak tahap desain, agar kapasitas peralatan HVAC dapat ditentukan secara tepat (Tsanas & Xifara, 2012).

Simulasi energi bangunan menghasilkan estimasi yang andal, tetapi memakan waktu dan menuntut keahlian perangkat lunak khusus. Model statistika menawarkan alternatif yang lebih cepat sekaligus memberi pemahaman kuantitatif mengenai parameter desain yang paling berkontribusi. Data Energy Efficiency (ENB2012) dari Tsanas dan Xifara (2012), yang tersedia pada UCI Machine Learning Repository, memuat 768 konfigurasi bangunan hasil simulasi Ecotect dengan delapan parameter desain (X1 sampai X8) dan dua respons, yaitu beban pemanasan (Y1) dan beban pendinginan (Y2).

Data ini memiliki tiga karakteristik yang menyulitkan regresi linear klasik. Pertama, luas permukaan (X2), luas dinding (X3), dan luas atap (X4) terikat oleh identitas geometri eksak \(X_2 = X_3 + 2X_4\) sehingga terjadi multikolinearitas sempurna, dan penduga kuadrat terkecil biasa (ordinary least squares, OLS) tidak memiliki solusi tunggal. Kedua, kekompakan relatif (X1) dan tinggi keseluruhan (X5) berkorelasi sangat kuat dengan variabel luas. Ketiga, sebaran Y1 dan Y2 tidak berbentuk lonceng tunggal, melainkan bimodal, yaitu memiliki dua puncak yang terpisah. Response yang bimodal menimbulkan pertanyaan apakah satu garis regresi dengan galat Gaussian tunggal sudah memadai, atau data lebih tepat dipandang sebagai campuran dari dua rezim desain dengan hubungan prediktor-respons yang berbeda.

Regresi Lasso (Tibshirani, 1996) menangani masalah pertama dan kedua karena penalti \(\ell_1\) membuat masalah optimasi tetap terdefinisi meskipun matriks desain singular, sekaligus menyusutkan sebagian koefisien hingga tepat nol sehingga terjadi seleksi variabel. Untuk masalah ketiga, laporan ini memperluas Lasso menjadi Mixture Lasso, yaitu model campuran regresi Gaussian berhingga dengan penalti \(\ell_1\) pada setiap komponen (Khalili & Chen, 2007; Städler et al., 2010). Model campuran diestimasi dengan algoritme expectation-maximization (EM), dan jumlah komponennya dipilih dengan Bayesian information criterion (BIC). Regresi linear (OLS) tetap disertakan sebagai model pembanding dan sebagai alat diagnostik. Seluruh analisis dikerjakan dengan bahasa pemrograman R.

1.2 Rumusan Masalah

Pertanyaan utama penelitian (Langkah 2 pada skrip analisis) adalah: variabel desain arsitektural manakah, di antara X1 sampai X8, yang tetap terpilih oleh Lasso (koefisien tidak sama dengan nol) saat memprediksi Heating Load (Y1) dan Cooling Load (Y2) secara terpisah, dan apakah bimodalitas kedua respons perlu dimodelkan sebagai campuran dua rezim? Pertanyaan tersebut dijabarkan sebagai berikut.

1. Bagaimana kualitas data, distribusi Y1 dan Y2, serta pola hubungan prediktor dengan respons, termasuk seberapa kuat bukti bimodalitas dan dari mana bimodalitas tersebut berasal?
2. Seberapa parah multikolinearitas pada data, dan bagaimana hasil diagnostik model OLS (rank matriks desain, condition number, VIF, homoskedastisitas, normalitas dan unimodalitas residual, serta observasi berpengaruh)?
3. Variabel apa saja yang terpilih oleh Lasso dan oleh Mixture Lasso (per komponen) untuk Y1 dan Y2, dan seberapa stabil seleksi tersebut pada validasi silang berulang?
4. Bagaimana kinerja prediksi OLS, Lasso, dan Mixture Lasso pada data uji dan validasi silang berulang, termasuk kejujuran interval prediksinya, dan apakah tambahan kompleksitas model campuran sepadan dengan perbaikan yang diperoleh?

1.3 Tujuan Penelitian

1. Mendeskripsikan kualitas data dan mendiagnosis bimodalitas Y1 dan Y2 beserta sumbernya.
2. Mendiagnosis multikolinearitas dan kelayakan asumsi model OLS pada data.
3. Menerapkan Lasso dan Mixture Lasso untuk seleksi variabel serta menilai stabilitas seleksinya.
4. Membandingkan kinerja prediksi titik dan interval prediksi ketiga model melalui data uji dan validasi silang berulang.

1.4 Manfaat Penelitian

Secara praktis, hasil penelitian memberi gambaran parameter desain yang paling berkontribusi terhadap beban pemanasan dan pendinginan pada ruang desain simulasi. Secara metodologis, laporan ini mendokumentasikan alur analisis yang dapat direproduksi untuk regresi terpenalti pada data dengan multikolinearitas sempurna dan respons bimodal, termasuk diagnosis bimodalitas, model campuran terpenalti, dan evaluasi interval prediksi.

1.5 Batasan Penelitian

1. Data berasal dari simulasi Ecotect bangunan hunian bervolume tetap di Athena, Yunani, dengan material dan okupansi tetap, sehingga kesimpulan hanya berlaku pada ruang konfigurasi simulasi tersebut.
2. Y1 dan Y2 dimodelkan secara terpisah, dan Y1 tidak dipakai untuk memprediksi Y2 (maupun sebaliknya).
3. X6 (orientasi) dan X8 (distribusi luas kaca) diperlakukan sebagai variabel kategorik.
4. Hubungan \(X_2 = X_3 + 2X_4\) dipertahankan sebagai masalah multikolinearitas, bukan alasan menghapus X2 sebelum Lasso.
5. Koefisien Lasso yang tidak nol dibaca sebagai “terpilih oleh model”, bukan “signifikan secara statistik”.
6. Model yang dipakai adalah model linear terpenalti tanpa suku interaksi maupun polinomial, baik secara tunggal maupun per komponen campuran.

Bab II Tinjauan Pustaka

2.1 Data Energy Efficiency (Tsanas & Xifara, 2012)

Tsanas dan Xifara (2012) membangun 768 bangunan simulasi yang masing-masing tersusun atas 18 kubus elementer berukuran \(3{,}5 \times 3{,}5 \times 3{,}5\) m dengan volume tetap 771,75 m\(^3\). Terdapat 12 bentuk bangunan, tiga tingkat luas kaca (10%, 25%, dan 40% dari luas lantai) dengan lima skenario distribusi, dan empat orientasi, ditambah bangunan tanpa kaca, sehingga total observasi adalah \(12 \times 3 \times 5 \times 4 + 12 \times 4 = 768\). Variabel pada data tersebut ditunjukkan pada Tabel 2.1.

Tabel 2.1 Variabel pada data ENB2012 (Tsanas & Xifara, 2012)

Kode	Variabel	Nilai unik
X1	Relative compactness (kekompakan relatif)	12
X2	Surface area (luas permukaan)	12
X3	Wall area (luas dinding)	7
X4	Roof area (luas atap)	4
X5	Overall height (tinggi keseluruhan)	2
X6	Orientation (orientasi)	4
X7	Glazing area (luas kaca)	4
X8	Glazing area distribution (distribusi luas kaca)	6
Y1	Heating load (beban pemanasan)	586
Y2	Cooling load (beban pendinginan)	636

Penelitian tersebut membandingkan regresi linear robust iteratively reweighted least squares (IRLS) dan random forest (RF) dengan validasi silang 10-lipat yang diulang 100 kali. Galat luar-sampel yang dilaporkan disajikan pada Tabel 2.2. Penulis menyimpulkan bahwa data bersifat non-Gaussian, RF jauh lebih akurat daripada IRLS, dan regresi klasik dapat menghasilkan koefisien yang besar dengan tanda berlawanan akibat multikolinearitas. Temuan terakhir menjadi motivasi penggunaan regresi terpenalti pada laporan ini, sedangkan sifat non-Gaussian (bimodal) menjadi motivasi model campuran.

Tabel 2.2 Galat luar-sampel pada Tsanas dan Xifara (2012)

Ukuran galat	Respons	IRLS	Random forest
MAE	Y1 (HL)	2,14 ± 0,24	0,51 ± 0,11
MAE	Y2 (CL)	2,21 ± 0,28	1,42 ± 0,25
MSE	Y1 (HL)	9,87 ± 2,41	1,03 ± 0,54
MSE	Y2 (CL)	11,46 ± 3,63	6,59 ± 1,56

2.2 Diagnosis Bimodalitas

Bimodalitas respons didiagnosis dengan empat ukuran yang saling melengkapi.

1. Jumlah puncak densitas kernel, yaitu banyaknya maksimum lokal pada estimasi densitas kernel \(\hat f(y)\). Ukuran ini deskriptif dan peka terhadap lebar pita (bandwidth).
2. Uji dip Hartigan (Hartigan & Hartigan, 1985), yang mengukur jarak maksimum antara fungsi distribusi empiris dan distribusi unimodal terdekat. Hipotesis nolnya adalah distribusi unimodal, sehingga p-value kecil menjadi bukti multimodalitas. Uji ini diimplementasikan pada paket diptest (Maechler, 2024).
3. Koefisien bimodalitas (bimodality coefficient, BC) (Pfister et al., 2013) dengan \(g\) kemencengan sampel dan \(\kappa\) kurtosis berlebih sampel (rumus di bawah). Nilai BC di atas \(5/9 \approx 0{,}555\) (nilai BC untuk distribusi seragam) mengindikasikan bimodalitas atau multimodalitas.
4. Campuran normal univariat yang diestimasi dengan paket mclust (Scrucca et al., 2016), dengan jumlah komponen \(G \in \{1, 2, 3, 4\}\) dipilih melalui BIC (Schwarz, 1978).

Koefisien bimodalitas pada butir 3 dihitung sebagai

\[\text{BC} = \frac{g^2 + 1}{\kappa + \dfrac{3(n-1)^2}{(n-2)(n-3)}}.\]

Bimodalitas marginal pada \(Y\) belum tentu memerlukan model campuran. Jika kedua puncak sepenuhnya dijelaskan oleh prediktor diskret yang teramati (misalnya X5 yang hanya bernilai 3,5 atau 7), maka \(Y\) akan unimodal di dalam tiap level prediktor tersebut, dan residual model linear juga akan unimodal. Model campuran baru dibenarkan apabila setelah prediktor dimasukkan, residual masih multimodal, atau apabila kinerja prediksi dan interval prediksinya membaik secara berarti.

2.3 Regresi Linear dan Diagnostik Model

Model regresi linear berganda adalah \(y = X\beta + \varepsilon\) dengan \(\varepsilon \sim (0, \sigma^2 I)\), dan penduga OLS

\[\hat{\beta}_{\text{OLS}} = (X^\top X)^{-1} X^\top y\]

mensyaratkan \(X^\top X\) dapat dibalik. Diagnostik yang dilakukan adalah sebagai berikut.

1. Rank dan condition number matriks desain (Belsley et al., 1980). Matriks desain rank deficient apabila rank-nya lebih kecil dari jumlah kolom, dan condition number yang sangat besar menandakan matriks hampir singular.
2. Variance inflation factor (VIF), yaitu \(\text{VIF}_j = 1/(1 - R_j^2)\) dengan \(R_j^2\) koefisien determinasi regresi prediktor ke-\(j\) pada prediktor lain. VIF di atas 10 lazim dianggap tanda multikolinearitas serius.
3. Uji Breusch-Pagan (Breusch & Pagan, 1979) versi studentized (Koenker, 1981) untuk homoskedastisitas, dengan hipotesis nol ragam galat konstan, diimplementasikan pada paket lmtest (Zeileis & Hothorn, 2002).
4. Uji Jarque-Bera (Jarque & Bera, 1980) untuk normalitas residual, yaitu \(\text{JB} = \frac{n}{6}\left(S^2 + \frac{(K-3)^2}{4}\right) \sim \chi^2_2\) dengan \(S\) kemencengan dan \(K\) kurtosis residual.
5. Uji dip pada residual untuk memeriksa apakah residual masih multimodal setelah prediktor dimasukkan.
6. Jarak Cook (Cook, 1977) dengan ambang \(4/n\) untuk observasi berpengaruh.

Karena terdapat multikolinearitas sempurna, uji-t dan p-value koefisien OLS tidak dipakai untuk menyatakan variabel “penting”. Diagnostik OLS dipakai untuk memahami kondisi data.

2.4 Multikolinearitas dan Regresi Lasso

Pada multikolinearitas sempurna, \(X^\top X\) singular sehingga \(\hat{\beta}_{\text{OLS}}\) tidak tunggal. Regularisasi mengatasinya dengan menambahkan penalti pada besar koefisien. Regresi ridge (Hoerl & Kennard, 1970) memakai penalti \(\ell_2\), sedangkan Lasso (Tibshirani, 1996; Hastie et al., 2009) memakai penalti \(\ell_1\):

\[\hat{\beta}^{\text{lasso}} = \underset{\beta_0,\,\beta}{\arg\min}\; \frac{1}{2n}\sum_{i=1}^{n}\left(y_i - \beta_0 - x_i^\top\beta\right)^2 + \lambda\sum_{j=1}^{p}\left|\beta_j\right|.\]

Fungsi objektif ini konveks sehingga selalu memiliki solusi meskipun \(X^\top X\) singular. Penalti \(\ell_1\) tidak dapat diturunkan di titik nol sehingga sebagian koefisien menjadi tepat nol. Karena penalti sensitif terhadap skala, prediktor numerik dibakukan terlebih dahulu. Nilai \(\lambda\) dipilih dengan validasi silang 5-lipat, yaitu \(\lambda\) yang meminimumkan rerata galat kuadrat (lambda.min). Solusi dihitung dengan coordinate descent pada paket glmnet (Friedman et al., 2010). Pada kelompok prediktor yang saling berkorelasi tinggi, Lasso cenderung memilih sebagian anggota kelompok saja (Zou & Hastie, 2005), sehingga tereliminasinya suatu variabel dibaca sebagai pemilihan perwakilan informasi, bukan bukti ketiadaan pengaruh.

2.5 Mixture Lasso

Model campuran regresi Gaussian dengan \(K\) komponen (rezim) menyatakan densitas bersyarat respons sebagai

\[f(y \mid x) = \sum_{k=1}^{K} \pi_k(x)\, \phi\!\left(y;\; \beta_{0k} + x^\top \beta_k,\; \sigma_k^2\right), \qquad \sum_{k=1}^{K}\pi_k(x) = 1,\]

dengan \(\phi(\cdot\,;\mu,\sigma^2)\) densitas normal, \(\beta_k\) vektor koefisien komponen ke-\(k\), dan \(\pi_k(x)\) peluang keanggotaan komponen. Seleksi variabel per komponen diperoleh dengan menambahkan penalti \(\ell_1\) pada setiap \(\beta_k\) (Khalili & Chen, 2007; Städler et al., 2010). Parameter diestimasi dengan algoritme EM (Dempster et al., 1977).

Langkah-E. Responsibility (peluang posterior keanggotaan) observasi ke-\(i\) pada komponen ke-\(k\) dihitung sebagai

\[w_{ik} = \frac{\pi_k\, \phi\!\left(y_i;\, \mu_{ik},\, \sigma_k^2\right)}{\sum_{l=1}^{K} \pi_l\, \phi\!\left(y_i;\, \mu_{il},\, \sigma_l^2\right)}, \qquad \mu_{ik} = \beta_{0k} + x_i^\top\beta_k.\]

Langkah-M. Untuk setiap komponen diselesaikan Lasso berbobot

\[\left(\hat\beta_{0k}, \hat\beta_k\right) = \underset{\beta_{0},\,\beta}{\arg\min}\; \frac{1}{2\sum_i w_{ik}} \sum_{i=1}^{n} w_{ik}\left(y_i - \beta_{0} - x_i^\top\beta\right)^2 + \lambda \lVert \beta \rVert_1,\]

lalu \(\hat\sigma_k^2 = \sum_i w_{ik}(y_i - \hat\mu_{ik})^2 / \sum_i w_{ik}\) dan \(\hat\pi_k = n^{-1}\sum_i w_{ik}\). Kedua langkah diulang hingga perubahan log-likelihood relatif lebih kecil dari \(10^{-6}\). Karena EM dapat terjebak pada optimum lokal, digunakan beberapa inisialisasi (k-means pada \(y\)), dan inisialisasi dengan komponen berbobot kurang dari 5% dibuang.

Pemilihan \(K\). Untuk \(K = 1, 2, 3\) dihitung

\[\text{BIC}_K = -2\,\ell(\hat\theta_K) + \log(n)\cdot \text{df}_K, \qquad \text{df}_K = \sum_{k=1}^{K}\#\{j: \hat\beta_{jk} \neq 0\} + K + K + (K-1),\]

yaitu banyaknya koefisien tidak nol, intersep, simpangan baku, dan proporsi campuran. Penggunaan banyaknya koefisien tidak nol sebagai derajat bebas Lasso mengikuti Zou et al. (2007). \(K\) dengan BIC terkecil dipilih. Pada \(K = 1\) model campuran tereduksi menjadi Lasso tunggal.

Penyempurnaan setelah \(K\) terpilih. Nilai \(\lambda\) setiap komponen dipilih ulang dengan validasi silang berbobot, dan peluang keanggotaan \(\pi_k(x)\) dipelajari dengan regresi logistik multinomial ber-penalti Lasso pada label komponen, mengikuti gagasan gating network pada mixture of experts (Jacobs et al., 1991), sehingga data baru dapat diprediksi.

2.6 Prediksi, Interval Prediksi, dan Residual Kuantil

Model campuran menghasilkan dua jenis prediksi titik. Prediksi soft adalah ekspektasi bersyarat

\[\hat y^{\,\text{soft}}(x) = \mathbb{E}[Y \mid x] = \sum_{k=1}^{K} \hat\pi_k(x)\,\hat\mu_k(x),\]

yang meminimumkan galat kuadrat. Prediksi hard adalah \(\hat y^{\,\text{hard}}(x) = \hat\mu_{k^*}(x)\) dengan \(k^* = \arg\max_k \hat\pi_k(x)\), yang mempertahankan struktur dua puncak tetapi berisiko salah komponen.

Interval prediksi \(100(1-\alpha)\%\) dari model campuran diperoleh dari kuantil \(\alpha/2\) dan \(1-\alpha/2\) fungsi distribusi campuran \(F(y \mid x) = \sum_k \hat\pi_k(x)\,\Phi\!\left((y - \hat\mu_k(x))/\hat\sigma_k\right)\). Sebagai pembanding, interval Gaussian tunggal untuk OLS dan Lasso adalah \(\hat y \pm z_{0,95}\, s\) dengan \(s\) simpangan baku residual data latih. Kualitas interval dinilai melalui coverage empiris (target 0,90) dan lebar rata-rata.

Kecocokan distribusi model campuran diperiksa dengan probability integral transform (PIT), \(u_i = F(y_i \mid x_i)\). Jika model benar, \(u_i \sim \text{Uniform}(0,1)\) (Gneiting et al., 2007), dan residual kuantil \(r_i = \Phi^{-1}(u_i)\) berdistribusi \(N(0,1)\) (Dunn & Smyth, 1996).

2.7 Ukuran Evaluasi dan Validasi Silang

Kinerja prediksi diukur dengan

\[\text{RMSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2},\qquad \text{MAE} = \frac{1}{n}\sum_{i=1}^{n}\lvert y_i-\hat{y}_i\rvert,\qquad R^2 = 1-\frac{\sum_{i}(y_i-\hat{y}_i)^2}{\sum_{i}(y_i-\bar{y})^2}.\]

Ukuran tersebut dihitung pada data uji dan pada validasi silang berulang 10-lipat dengan 5 pengulangan. Pembagian data latih-uji dan lipatan validasi silang distratifikasi menurut nilai respons agar kedua puncak terwakili di setiap bagian. Ringkasan validasi silang disajikan sebagai rerata, simpangan baku, dan interval empiris 2,5% sampai 97,5%. Interval empiris ini bukan selang kepercayaan inferensial klasik karena lipatan-lipatan validasi silang saling bergantung (Bates et al., 2024).

Bab III Metode Penelitian

3.1 Data dan Variabel

Data yang digunakan adalah ENB2012 (768 observasi, 10 variabel) dengan unit observasi berupa konfigurasi desain bangunan hasil simulasi. Prediktor numerik adalah X1, X2, X3, X4, X5, dan X7, sedangkan X6 dan X8 diperlakukan sebagai kategorik. Untuk model linear dan campuran, prediktor numerik dibakukan (rerata 0, simpangan baku populasi 1) dan X6 serta X8 dikodekan menjadi variabel dummy dengan kategori pertama sebagai referensi, sehingga terdapat 14 kolom prediktor (6 numerik, 3 dummy X6, dan 5 dummy X8). Parameter pembakuan dan daftar level dipelajari hanya dari data latih untuk mencegah kebocoran informasi.

3.2 Tahapan Analisis

1. Langkah 2, rumusan pertanyaan, seperti pada Subbab 1.2.
2. Langkah 3, EDA terarah: kualitas data, identitas \(X_2 = X_3 + 2X_4\), distribusi Y1 dan Y2, diagnosis bimodalitas (puncak densitas, uji dip, BC, mclust), sumber bimodalitas menurut X5, korelasi Spearman, kurva LOWESS X7, dan boxplot prediktor kategorik.
3. Langkah 4, estimasi model: pembagian data latih-uji 80:20 yang terstratifikasi (set.seed(42)), lalu estimasi OLS, Lasso, dan Mixture Lasso untuk Y1 dan Y2 secara terpisah, serta seleksi variabel pada level X1 sampai X8.
4. Langkah 5, diagnostik model: rank, condition number, VIF, Breusch-Pagan, Jarque-Bera, uji dip residual, jarak Cook, plot residual pada OLS, dan residual kuantil (PIT) pada Mixture Lasso.
5. Langkah 6, evaluasi: RMSE, MAE, dan \(R^2\) pada data uji; coverage interval prediksi 90%; validasi silang berulang terstratifikasi (10-lipat \(\times\) 5 pengulangan) dengan pemilihan ulang \(\lambda\) dan \(K\) di setiap lipatan latih luar; frekuensi \(K\) terpilih; serta stabilitas seleksi variabel.
6. Interpretasi: menjawab pertanyaan penelitian, membedakan inferensi dan prediksi, serta menyatakan keterbatasan.

3.3 Pengaturan Model

Tabel 3.1 Pengaturan model

Komponen	Pengaturan
OLS	lm pada matriks desain terbakukan; koefisien teralias bernilai NA, prediksi tetap valid
Lasso	glmnet, alpha = 1, standardize = FALSE, kisi 200 nilai \(\lambda\) dari \(10^{2}\) hingga \(10^{-4}\), CV dalam 5-lipat, lambda.min
Mixture Lasso	EM, \(K \in \{1,2,3\}\), 3 inisialisasi per \(K\), maksimum 100 iterasi, toleransi \(10^{-6}\), proporsi minimum 5%, \(\lambda\) awal = lambda.min Lasso global, \(K\) dipilih BIC
Penyempurnaan campuran	\(\lambda\) per komponen dipilih ulang dengan CV berbobot; \(\pi_k(x)\) dengan Lasso multinomial
Validasi	Uji 20% terstratifikasi; CV berulang 10-lipat \(\times\) 5 terstratifikasi (set.seed(42 + ulangan))

3.4 Perangkat Lunak

Analisis menggunakan R version 4.6.1 (2026-06-24 ucrt) (R Core Team, 2024) dengan paket tidyverse (Wickham et al., 2019), readxl (Wickham & Bryan, 2023), glmnet, lmtest, patchwork (Pedersen, 2024), diptest, mclust, serta knitr dan rmarkdown untuk penyusunan laporan. Paket yang belum terpasang dipasang otomatis pada chunk setup. Seluruh angka pada narasi dihitung langsung dari hasil eksekusi kode.

Bab IV Hasil dan Pembahasan

4.1 Persiapan Data

4.1.1 Memuat Paket

suppressPackageStartupMessages({
  library(tidyverse)   # dplyr, tidyr, ggplot2, purrr, tibble
  library(readxl)      # membaca file .xlsx
  library(glmnet)      # Lasso
  library(lmtest)      # uji Breusch-Pagan
  library(patchwork)   # menggabungkan grafik ggplot
  library(diptest)     # uji dip Hartigan (unimodalitas)
  library(mclust)      # campuran normal univariat
})

# PENTING: beberapa paket menimpa fungsi tidyverse bila ikut termuat, misalnya
# mclust::map, plyr::count/summarise/mutate/arrange/rename/desc,
# MASS::select, dan stats::filter. Fungsi-fungsi berikut dikunci ke versi
# tidyverse agar laporan tetap berjalan apa pun paket lain yang aktif di R Anda.
map       <- purrr::map
map_dfr   <- purrr::map_dfr
map2_dbl  <- purrr::map2_dbl
count     <- dplyr::count
summarise <- dplyr::summarise
summarize <- dplyr::summarize
mutate    <- dplyr::mutate
arrange   <- dplyr::arrange
rename    <- dplyr::rename
desc      <- dplyr::desc
select    <- dplyr::select
filter    <- dplyr::filter
group_by  <- dplyr::group_by
ungroup   <- dplyr::ungroup

theme_set(theme_bw(base_size = 12) +
            theme(text = element_text(family = "serif")))
options(width = 160, tibble.width = Inf, tibble.print_max = Inf)

4.1.2 Membaca Data

Chunk berikut mencari berkas data pada beberapa lokasi. Jika berkas tidak ditemukan, data diunduh otomatis dari UCI Machine Learning Repository. Ubah FILE_PATH bila berkas disimpan di lokasi lain.

FILE_PATH <- "C:/Users/USER/Downloads/energy+efficiency/ENB2012_data.xlsx"

candidates <- c(FILE_PATH,
                "ENB2012_data.xlsx",
                file.path("energy+efficiency", "ENB2012_data.xlsx"),
                "C:/Users/USER/Downloads/Energy Efficiency.xlsx")
found <- candidates[file.exists(candidates)]

if (length(found) > 0) {
  data_file <- found[1]
} else {
  # Unduh otomatis dari UCI bila berkas lokal tidak tersedia
  zip_file <- file.path(tempdir(), "energy_efficiency.zip")
  ok <- tryCatch({
    download.file("https://archive.ics.uci.edu/static/public/242/energy+efficiency.zip",
                  zip_file, mode = "wb", quiet = TRUE)
    unzip(zip_file, exdir = tempdir())
    TRUE
  }, error = function(e) FALSE)
  xl <- list.files(tempdir(), pattern = "ENB2012.*\\.xlsx$", recursive = TRUE, full.names = TRUE)
  if (!ok || length(xl) == 0) {
    stop("Berkas ENB2012_data.xlsx tidak ditemukan. Ubah FILE_PATH ke lokasi berkas Anda.")
  }
  data_file <- xl[1]
}

df_raw <- read_excel(data_file, sheet = 1)

# Nama kolom standar X1..X8, Y1, Y2 (berkas UCI sudah memakai nama ini)
if (!all(c(paste0("X", 1:8), "Y1", "Y2") %in% names(df_raw))) {
  names(df_raw)[1:10] <- c(paste0("X", 1:8), "Y1", "Y2")
}

# Ambil 10 kolom utama dan buang baris yang seluruhnya kosong
df <- df_raw %>%
  select(all_of(c(paste0("X", 1:8), "Y1", "Y2"))) %>%
  mutate(across(everything(), as.numeric)) %>%
  filter(!if_all(everything(), is.na))

cat("Berkas dibaca :", data_file, "\n")

#> Berkas dibaca : ENB2012_data.xlsx

cat("Shape         :", nrow(df), "x", ncol(df), "\n")

#> Shape         : 768 x 10

head(df)

X1	X2	X3	X4	X5	X6	X7	X8	Y1	Y2
0.98	514.5	294.0	110.25	7	2	0	0	15.55	21.33
0.98	514.5	294.0	110.25	7	3	0	0	15.55	21.33
0.98	514.5	294.0	110.25	7	4	0	0	15.55	21.33
0.98	514.5	294.0	110.25	7	5	0	0	15.55	21.33
0.90	563.5	318.5	122.50	7	2	0	0	20.84	28.28
0.90	563.5	318.5	122.50	7	3	0	0	21.46	25.38

4.1.3 Definisi Variabel

FEATURES             <- c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")
TARGETS              <- c("Y1", "Y2")
NUMERIC_FEATURES     <- c("X1", "X2", "X3", "X4", "X5", "X7")
CATEGORICAL_FEATURES <- c("X6", "X8")

variable_info <- tibble(
  Variable = c(FEATURES, TARGETS),
  Description = c("Relative Compactness", "Surface Area", "Wall Area", "Roof Area",
                  "Overall Height", "Orientation", "Glazing Area",
                  "Glazing Area Distribution", "Heating Load", "Cooling Load"),
  Role = c(rep("Predictor", 8), rep("Response", 2))
)
tbl(variable_info, "Tabel 4.1 Definisi variabel")

Tabel 4.1 Definisi variabel

Variable	Description	Role
X1	Relative Compactness	Predictor
X2	Surface Area	Predictor
X3	Wall Area	Predictor
X4	Roof Area	Predictor
X5	Overall Height	Predictor
X6	Orientation	Predictor
X7	Glazing Area	Predictor
X8	Glazing Area Distribution	Predictor
Y1	Heating Load	Response
Y2	Cooling Load	Response

Data terdiri atas 768 observasi dengan delapan prediktor dan dua respons. Y1 adalah beban pemanasan dan Y2 adalah beban pendinginan.

4.2 EDA Terarah (Langkah 3)

4.2.1 Kualitas Data dan Statistika Deskriptif

vars <- c(FEATURES, TARGETS)

quality <- tibble(
  Variable = vars,
  dtype    = vapply(df[vars], function(x) class(x)[1], character(1), USE.NAMES = FALSE),
  missing  = vapply(df[vars], function(x) sum(is.na(x)), numeric(1), USE.NAMES = FALSE),
  n_unique = vapply(df[vars], function(x) n_distinct(x), numeric(1), USE.NAMES = FALSE)
)
tbl(quality, "Tabel 4.2 Kualitas data")

Tabel 4.2 Kualitas data

Variable	dtype	missing	n_unique
X1	numeric	0	12
X2	numeric	0	12
X3	numeric	0	7
X4	numeric	0	4
X5	numeric	0	2
X6	numeric	0	4
X7	numeric	0	4
X8	numeric	0	6
Y1	numeric	0	587
Y2	numeric	0	636

cat("Duplicate rows:", sum(duplicated(df)), "\n")

#> Duplicate rows: 0

describe_tbl <- function(data) {
  q <- function(p) vapply(data, function(x) quantile(x, p, na.rm = TRUE, names = FALSE),
                          numeric(1), USE.NAMES = FALSE)
  tibble(
    Variable = names(data),
    count    = vapply(data, function(x) sum(!is.na(x)), numeric(1), USE.NAMES = FALSE),
    mean     = vapply(data, mean, numeric(1), na.rm = TRUE, USE.NAMES = FALSE),
    std      = vapply(data, sd,   numeric(1), na.rm = TRUE, USE.NAMES = FALSE),
    min      = vapply(data, min,  numeric(1), na.rm = TRUE, USE.NAMES = FALSE),
    `25%`    = q(0.25), `50%` = q(0.50), `75%` = q(0.75),
    max      = vapply(data, max,  numeric(1), na.rm = TRUE, USE.NAMES = FALSE)
  )
}
tbl(describe_tbl(df[vars]), "Tabel 4.3 Statistika deskriptif")

Tabel 4.3 Statistika deskriptif

Variable	count	mean	std	min	25%	50%	75%	max
X1	768	0.7642	0.1058	0.62	0.6825	0.75	0.8300	0.98
X2	768	671.7083	88.0861	514.50	606.3750	673.75	741.1250	808.50
X3	768	318.5000	43.6265	245.00	294.0000	318.50	343.0000	416.50
X4	768	176.6042	45.1660	110.25	140.8750	183.75	220.5000	220.50
X5	768	5.2500	1.7511	3.50	3.5000	5.25	7.0000	7.00
X6	768	3.5000	1.1188	2.00	2.7500	3.50	4.2500	5.00
X7	768	0.2344	0.1332	0.00	0.1000	0.25	0.4000	0.40
X8	768	2.8125	1.5510	0.00	1.7500	3.00	4.0000	5.00
Y1	768	22.3072	10.0902	6.01	12.9925	18.95	31.6675	43.10
Y2	768	24.5878	9.5133	10.90	15.6200	22.08	33.1325	48.03

Dataset memiliki 768 observasi dengan 0 nilai hilang dan 0 baris duplikat, sehingga tidak diperlukan imputasi maupun penghapusan baris. Prediktor X1 sampai X8 hanya memiliki 2 sampai 12 nilai unik, sedangkan Y1 dan Y2 masing-masing memiliki 587 dan 636 nilai unik. Hal ini sesuai dengan rancangan simulasi, yaitu prediktor merupakan parameter desain diskret dan respons berupa besaran kontinu. Sedikitnya nilai unik pada X6 dan X8 yang merupakan kode kategori mendasari perlakuan keduanya sebagai variabel kategorik.

4.2.2 Hubungan Struktural X2, X3, dan X4

dependency <- df$X2 - (df$X3 + 2 * df$X4)
cat("Maximum absolute difference X2 - (X3 + 2*X4):", max(abs(dependency)), "\n")

#> Maximum absolute difference X2 - (X3 + 2*X4): 0

Selisih absolut maksimum antara X2 dan \(X_3 + 2X_4\) adalah 0 (hanya galat pembulatan numerik), sehingga identitas \(X_2 = X_3 + 2X_4\) berlaku pada seluruh observasi. Ketiga variabel ini bergantung linear secara sempurna dan akan membuat matriks desain OLS singular.

4.2.3 Distribusi Y1 dan Y2

plot_dist <- function(target) {
  p1 <- ggplot(df, aes(x = .data[[target]])) +
    geom_histogram(aes(y = after_stat(density)), bins = 30,
                   fill = "steelblue", alpha = 0.6, color = "white") +
    geom_density(linewidth = 1) +
    labs(title = paste("Distribution of", target), x = target, y = "Density")
  p2 <- ggplot(df, aes(x = .data[[target]])) +
    geom_boxplot(fill = "steelblue", alpha = 0.6) +
    labs(title = paste("Boxplot of", target), x = target) +
    theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
  p1 | p2
}
print(plot_dist("Y1") / plot_dist("Y2"))

Gambar 4.1 Histogram, densitas, dan boxplot Y1 dan Y2

Rerata Y1 adalah 22,31 (median 18,95, simpangan baku 10,09), dan rerata Y2 adalah 24,59 (median 22,08, simpangan baku 9,51). Histogram kedua respons memperlihatkan kepadatan yang terkonsentrasi pada dua wilayah nilai dengan lembah di antaranya, sehingga sebaran respons tidak dapat dihampiri oleh satu distribusi normal. Bukti formal bimodalitas diuji pada subbab berikut.

4.2.4 Diagnosis Bimodalitas Y1 dan Y2

count_modes <- function(x, adjust = 1) {
  d <- density(x, adjust = adjust)
  sum(diff(sign(diff(d$y))) == -2)
}

bimodality_coef <- function(x) {
  n <- length(x); m <- x - mean(x)
  g <- mean(m^3) / mean(m^2)^1.5
  k <- mean(m^4) / mean(m^2)^2 - 3
  (g^2 + 1) / (k + 3 * (n - 1)^2 / ((n - 2) * (n - 3)))
}

bimodality_tbl <- map_dfr(TARGETS, function(t) {
  y  <- df[[t]]
  dp <- dip.test(y)
  mc <- Mclust(y, G = 1:4, verbose = FALSE)
  tibble(
    Target       = t,
    Modes        = count_modes(y),
    Dip_stat     = unname(dp$statistic),
    Dip_p        = dp$p.value,
    BC           = bimodality_coef(y),
    Mclust_G     = mc$G,
    Mclust_means = paste(round(sort(mc$parameters$mean), 2), collapse = " | ")
  )
})
tbl(bimodality_tbl, "Tabel 4.4 Diagnosis bimodalitas respons")

Tabel 4.4 Diagnosis bimodalitas respons

Target	Modes	Dip_stat	Dip_p	BC	Mclust_G	Mclust_means
Y1	2	0.0468	0.00e+00	0.6394	3	13.52 | 28.36 | 37.29
Y2	2	0.0400	0.00e+00	0.6202	2	15.54 | 31.73

Untuk Y1, densitas kernel memiliki 2 puncak, uji dip menghasilkan p-value 0,00e+00, BC sebesar 0,639 (ambang 0,555), dan mclust memilih 3 komponen dengan rerata 13.52 | 28.36 | 37.29. Untuk Y2, densitas memiliki 2 puncak, p-value uji dip 0,00e+00, BC 0,620, dan mclust memilih 2 komponen dengan rerata 15.54 | 31.73. Dengan demikian, 4 dari 4 ukuran mendukung multimodalitas pada Y1 dan 4 dari 4 ukuran pada Y2. Perlu dicatat bahwa p-value uji dip yang kecil menolak unimodalitas, sedangkan BC dan mclust bersifat deskriptif dan dapat dipengaruhi kemencengan.

4.2.5 Sumber Bimodalitas: Peran X5 (Overall Height)

within_x5 <- map_dfr(TARGETS, function(t) {
  df %>%
    group_by(X5) %>%
    summarise(Target = t, n = n(), Mean = mean(.data[[t]]), SD = sd(.data[[t]]),
              Modes = count_modes(.data[[t]]),
              Dip_p = dip.test(.data[[t]])$p.value, .groups = "drop")
}) %>% relocate(Target)
tbl(within_x5, "Tabel 4.5 Distribusi respons di dalam tiap level X5")

Tabel 4.5 Distribusi respons di dalam tiap level X5

Target	X5	n	Mean	SD	Modes	Dip_p
Y1	3.5	384	13.3385	2.6718	3	0.00e+00
Y1	7.0	384	31.2759	5.9546	3	7.05e-06
Y2	3.5	384	16.0714	2.4208	2	3.61e-01
Y2	7.0	384	33.1041	5.4724	2	1.24e-03

plot_by_x5 <- function(target) {
  ggplot(df, aes(x = .data[[target]], fill = factor(X5))) +
    geom_histogram(bins = 40, alpha = 0.7, position = "identity", color = "white") +
    labs(title = paste(target, "menurut X5 (Overall Height)"), x = target, fill = "X5")
}
print(plot_by_x5("Y1") | plot_by_x5("Y2"))

Gambar 4.2 Sebaran Y1 dan Y2 menurut X5

X5 hanya memiliki dua nilai, yaitu bangunan rendah dan tinggi. Untuk Y1 diperoleh rerata 13,34 pada X5 = 3,5 dan 31,28 pada X5 = 7,0, sedangkan untuk Y2 diperoleh rerata 16,07 pada X5 = 3,5 dan 33,10 pada X5 = 7,0. Histogram memperlihatkan bahwa kedua puncak sebagian besar berimpit dengan dua level X5. Di dalam tiap level X5, uji dip tidak menolak unimodalitas pada 0 dari 2 level untuk Y1 dan 1 dari 2 level untuk Y2 (taraf 5%). Artinya, bimodalitas marginal sebagian besar merupakan campuran dari dua subpopulasi yang teramati melalui X5 (dan variabel geometri yang berkorelasi dengannya). Apakah masih ada sisa bimodalitas setelah seluruh prediktor dimasukkan diperiksa melalui residual OLS pada Subbab 4.4.

4.2.6 Korelasi Peringkat Spearman

spearman_cols <- c("X1", "X2", "X3", "X4", "X5", "X7", "Y1", "Y2")
spearman_corr <- cor(df[spearman_cols], method = "spearman")

spearman_long <- as.data.frame(as.table(spearman_corr)) %>% setNames(c("Var1", "Var2", "rho"))

p_heat <- ggplot(spearman_long, aes(x = Var1, y = Var2, fill = rho)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.2f", rho)), size = 3.5, family = "serif") +
  scale_fill_gradient2(low = "#3B4CC0", mid = "white", high = "#B40426",
                       midpoint = 0, limits = c(-1, 1)) +
  scale_x_discrete(limits = spearman_cols) +
  scale_y_discrete(limits = rev(spearman_cols)) +
  coord_fixed() +
  labs(title = "Spearman Rank Correlation", x = NULL, y = NULL)
print(p_heat)

Gambar 4.3 Matriks korelasi peringkat Spearman

rho_tbl <- as.data.frame(round(spearman_corr[c("X1", "X2", "X3", "X4", "X5", "X7"), c("Y1", "Y2")], 4))
tbl(rownames_to_column(rho_tbl, "Variable"), "Tabel 4.6 Korelasi Spearman prediktor numerik dengan respons")

Tabel 4.6 Korelasi Spearman prediktor numerik dengan respons

Variable	Y1	Y2
X1	0.6221	0.6510
X2	-0.6221	-0.6510
X3	0.4715	0.4160
X4	-0.8040	-0.8032
X5	0.8613	0.8649
X7	0.3229	0.2889

Korelasi Spearman dipakai karena prediktor diskret dan hubungan monoton tak linear tetap tertangkap. Prediktor dengan korelasi terkuat terhadap Y1 adalah X5 (\(\rho\) = 0,861) dan terhadap Y2 adalah X5 (\(\rho\) = 0,865). Kedua respons berkorelasi sangat tinggi (\(\rho\) = 0,973). Di antara prediktor numerik, pasangan terkuat adalah X2 dan X1 (\(\rho\) = -1,000), yang menegaskan struktur kolinearitas kuat. Korelasi kecil tidak berarti prediktor tidak penting jika hubungannya tidak monoton atau bersifat interaksi.

4.2.7 Fokus X7 (Glazing Area) dan Potensi Ketaklinearan

plot_lowess <- function(target) {
  lw <- lowess(df$X7, df[[target]], f = 2/3, iter = 3, delta = 0)
  ggplot(df, aes(x = X7, y = .data[[target]])) +
    geom_point(alpha = 0.35) +
    geom_line(data = tibble(X7 = lw$x, fit = lw$y), aes(x = X7, y = fit),
              inherit.aes = FALSE, color = "firebrick", linewidth = 1) +
    labs(title = paste("LOWESS: X7 vs", target))
}
print(plot_lowess("Y1") | plot_lowess("Y2"))

Gambar 4.4 Kurva LOWESS X7 terhadap Y1 dan Y2

Kurva LOWESS menilai apakah hubungan luas kaca dengan respons mendekati garis lurus. X7 hanya memiliki empat nilai sehingga titik berkelompok pada empat posisi, dan sebaran vertikal yang lebar pada setiap posisi mencerminkan pengaruh prediktor lain (terutama X5). Penyimpangan kurva dari garis lurus menjadi indikasi awal ketaklinearan yang tidak dapat ditangkap model linear tanpa suku tambahan.

4.2.8 Prediktor Kategorik terhadap Y1 dan Y2

box_cat <- function(x, y, title) {
  ggplot(df, aes(x = factor(.data[[x]]), y = .data[[y]])) +
    geom_boxplot(fill = "steelblue", alpha = 0.6) +
    labs(title = title, x = x, y = y)
}
print(
  (box_cat("X6", "Y1", "Y1 by Orientation (X6)") | box_cat("X8", "Y1", "Y1 by Glazing Distribution (X8)")) /
    (box_cat("X6", "Y2", "Y2 by Orientation (X6)") | box_cat("X8", "Y2", "Y2 by Glazing Distribution (X8)"))
)

Gambar 4.5 Boxplot Y1 dan Y2 menurut X6 dan X8

Sebaran Y1 dan Y2 yang hampir sama antarkategori X6 menandakan kontribusi orientasi yang kecil. Pada X8, kategori 0 (tanpa kaca) cenderung memiliki respons lebih rendah, karena kategori ini berimpit dengan X7 = 0. Kesimpulan visual ini dibandingkan dengan hasil seleksi Lasso pada subbab berikutnya.

4.3 Estimasi Lasso dan Mixture Lasso (Langkah 4)

4.3.1 Pembagian Data Latih dan Uji Terstratifikasi

strat_folds <- function(key, k) {
  n <- length(key); ord <- order(key, runif(n)); f <- integer(n)
  blocks <- ceiling(seq_len(n) / k)
  for (b in unique(blocks)) {
    idx <- which(blocks == b)
    f[ord[idx]] <- sample(k, length(idx))
  }
  f
}
# Kunci stratifikasi: rerata peringkat Y1 dan Y2 (keduanya sangat berkorelasi)
strat_key <- (rank(df$Y1) + rank(df$Y2)) / 2

set.seed(42)
n_obs    <- nrow(df)
test_idx <- which(strat_folds(strat_key, 5) == 1)     # sekitar 20% test

X_train <- df[-test_idx, FEATURES];  X_test <- df[test_idx, FEATURES]
Y_train <- df[-test_idx, TARGETS];   Y_test <- df[test_idx, TARGETS]

cat("X_train:", nrow(X_train), "x", ncol(X_train), "\n")

#> X_train: 614 x 8

cat("X_test :", nrow(X_test),  "x", ncol(X_test),  "\n")

#> X_test : 154 x 8

Data diurutkan menurut kunci stratifikasi, dibagi menjadi blok berukuran lima, dan setiap anggota blok diberi label lipatan secara acak. Satu lipatan menjadi data uji, sehingga diperoleh 614 observasi latih dan 154 observasi uji. Cara ini menjamin kedua puncak terwakili secara proporsional di data latih maupun data uji.

4.3.2 Preprocessing

fit_prep <- function(X) {
  sc <- vapply(X[NUMERIC_FEATURES], function(v) sqrt(mean((v - mean(v))^2)), numeric(1))
  sc[sc == 0] <- 1
  list(
    center = vapply(X[NUMERIC_FEATURES], mean, numeric(1)),
    scale  = sc,
    levels = lapply(X[CATEGORICAL_FEATURES], function(v) sort(unique(v)))
  )
}

apply_prep <- function(prep, X) {
  num <- do.call(cbind, lapply(NUMERIC_FEATURES, function(v) {
    (X[[v]] - prep$center[[v]]) / prep$scale[[v]]
  }))
  colnames(num) <- NUMERIC_FEATURES

  cat_list <- lapply(CATEGORICAL_FEATURES, function(v) {
    lv <- prep$levels[[v]][-1]
    m  <- matrix(vapply(lv, function(l) as.numeric(X[[v]] == l), numeric(nrow(X))),
                 nrow = nrow(X), ncol = length(lv))
    colnames(m) <- paste0(v, "_", lv)
    m
  })
  cbind(num, do.call(cbind, cat_list))
}

4.3.3 Fungsi Mixture Lasso

Fungsi berikut mengikuti skrip final. Dua pengaman ditambahkan tanpa mengubah hasil ketika data normal: (i) bila suatu komponen terlalu kecil untuk validasi silang berbobot, koefisien komponen memakai Lasso berbobot dengan \(\lambda\) acuan; (ii) bila suatu komponen memiliki terlalu sedikit anggota untuk Lasso multinomial, \(\pi_k(x)\) diganti proporsi campuran konstan. Pengaman ini mencegah knit gagal pada lipatan validasi silang tertentu.

K_MAX        <- 3       # jumlah komponen maksimum
N_START      <- 3       # inisialisasi acak EM per K
EM_MAX_ITER  <- 100
EM_TOL       <- 1e-6
MIN_PROP     <- 0.05    # komponen < 5% dianggap runtuh (inisialisasi dibuang)
MIN_MEMBER   <- 10      # anggota minimum per komponen untuk Lasso multinomial

LAMBDA_GRID <- 10^seq(2, -4, length.out = 200)   # urut menurun
INNER_FOLDS <- 5

wlasso <- function(Z, y, w, lambda) {
  keep <- w > 1e-8
  fit <- suppressWarnings(glmnet(Z[keep, , drop = FALSE], y[keep], weights = w[keep],
                                 alpha = 1, lambda = lambda, standardize = FALSE,
                                 thresh = 1e-10, maxit = 1e6))
  as.numeric(as.matrix(coef(fit, s = lambda)))       # intercept diikuti koefisien
}

em_mixture <- function(Z, y, K, lambda, n_start = N_START) {
  n <- length(y); best <- NULL
  sig_floor <- 0.01 * sd(y)
  for (s in seq_len(n_start)) {
    if (K == 1) {
      w <- matrix(1, n, 1)
    } else {
      yy <- if (s == 1) y else y + rnorm(n, 0, 0.3 * sd(y))
      cl <- kmeans(yy, centers = K, nstart = 1)$cluster
      w  <- matrix(0.05 / (K - 1), n, K)
      w[cbind(seq_len(n), cl)] <- 0.95
    }
    ll_old <- -Inf; ok <- TRUE; ll <- -Inf; B <- NULL; sg <- NULL
    for (it in seq_len(EM_MAX_ITER)) {
      pi_k <- colMeans(w)
      if (any(pi_k < MIN_PROP)) { ok <- FALSE; break }
      B  <- sapply(seq_len(K), function(j) wlasso(Z, y, w[, j], lambda))   # (p+1) x K
      B  <- matrix(B, ncol = K)
      MU <- cbind(1, Z) %*% B
      sg <- pmax(sqrt(colSums(w * (y - MU)^2) / colSums(w)), sig_floor)
      logd <- sapply(seq_len(K), function(j) log(pi_k[j]) + dnorm(y, MU[, j], sg[j], log = TRUE))
      logd <- matrix(logd, nrow = n, ncol = K)
      m    <- apply(logd, 1, max)
      lse  <- m + log(rowSums(exp(logd - m)))
      ll   <- sum(lse)
      w    <- exp(logd - lse)
      if (abs(ll - ll_old) < EM_TOL * (abs(ll) + 1)) break
      ll_old <- ll
    }
    if (!ok || !is.finite(ll)) next
    df_par <- sum(abs(B[-1, ]) > 1e-10) + K + K + (K - 1)
    bic    <- -2 * ll + log(n) * df_par
    if (is.null(best) || bic < best$bic) {
      best <- list(K = K, ll = ll, bic = bic, B = B, sigma = sg, w = w)
    }
  }
  best
}

fit_mixture_target <- function(X_tr, y_tr, lambda_ref, seed = 42) {
  set.seed(seed)
  prep <- fit_prep(X_tr)
  Z    <- apply_prep(prep, X_tr)

  # Pilih K dengan BIC
  cand <- lapply(seq_len(K_MAX), function(K) em_mixture(Z, y_tr, K, lambda_ref))
  bics <- vapply(cand, function(f) if (is.null(f)) NA_real_ else f$bic, numeric(1))
  names(bics) <- seq_len(K_MAX)
  K  <- as.integer(names(which.min(bics)))
  em <- cand[[K]]

  # Urutkan komponen menurut rerata Y (komponen 1 = puncak terendah)
  comp_mean <- colSums(em$w * y_tr) / colSums(em$w)
  post <- em$w[, order(comp_mean), drop = FALSE]

  # Refit Lasso berbobot per komponen dengan lambda hasil CV berbobot
  B <- matrix(0, ncol(Z) + 1, K, dimnames = list(c("(Intercept)", colnames(Z)), NULL))
  for (j in seq_len(K)) {
    idx <- which(post[, j] > 1e-3)
    B[, j] <- tryCatch({
      if (length(idx) < 5 * INNER_FOLDS) stop("komponen terlalu kecil")
      foldid <- sample(rep(seq_len(INNER_FOLDS), length.out = length(idx)))
      cvj <- cv.glmnet(Z[idx, , drop = FALSE], y_tr[idx], weights = post[idx, j],
                       alpha = 1, lambda = LAMBDA_GRID, foldid = foldid,
                       standardize = FALSE, thresh = 1e-10, maxit = 1e6)
      as.numeric(as.matrix(coef(cvj, s = "lambda.min")))
    }, error = function(e) wlasso(Z, y_tr, post[, j], lambda_ref))
  }
  MU    <- cbind(1, Z) %*% B
  sigma <- pmax(sqrt(colSums(post * (y_tr - MU)^2) / colSums(post)), 0.01 * sd(y_tr))
  hard  <- max.col(post)

  # pi_k(X): Lasso multinomial pada label komponen (hard assignment)
  member <- NULL
  if (K > 1 && all(tabulate(hard, nbins = K) >= MIN_MEMBER)) {
    member <- tryCatch(
      suppressWarnings(cv.glmnet(Z, factor(hard, levels = seq_len(K)), family = "multinomial",
                                 alpha = 1, nfolds = 5, type.measure = "class",
                                 standardize = FALSE)),
      error = function(e) NULL)
  }

  list(prep = prep, K = K, B = B, sigma = sigma, weight = colMeans(post),
       post = post, hard = hard, member = member, bic = bics)
}

mix_parts <- function(mm, X_new) {
  Z  <- apply_prep(mm$prep, X_new)
  MU <- cbind(1, Z) %*% mm$B
  PI <- if (mm$K == 1) {
    matrix(1, nrow(Z), 1)
  } else if (is.null(mm$member)) {
    matrix(mm$weight, nrow(Z), mm$K, byrow = TRUE)     # fallback: proporsi konstan
  } else {
    matrix(predict(mm$member, newx = Z, s = "lambda.min", type = "response"), nrow = nrow(Z))
  }
  list(MU = MU, PI = PI)
}

predict_mixture <- function(mm, X_new, type = c("soft", "hard")) {
  type <- match.arg(type)
  pr <- mix_parts(mm, X_new)
  if (type == "soft") rowSums(pr$PI * pr$MU) else pr$MU[cbind(seq_len(nrow(pr$MU)), max.col(pr$PI))]
}

mixture_interval <- function(mm, X_new, level = 0.90, ngrid = 800) {
  pr <- mix_parts(mm, X_new); n <- nrow(pr$MU)
  grid <- seq(min(pr$MU) - 4 * max(mm$sigma), max(pr$MU) + 4 * max(mm$sigma), length.out = ngrid)
  cdf <- matrix(0, n, ngrid)
  for (k in seq_len(mm$K)) {
    cdf <- cdf + pr$PI[, k] * pnorm((matrix(grid, n, ngrid, byrow = TRUE) - pr$MU[, k]) / mm$sigma[k])
  }
  a <- (1 - level) / 2
  cbind(lo = apply(cdf, 1, function(r) grid[which(r >= a)[1]]),
        hi = apply(cdf, 1, function(r) grid[which(r >= 1 - a)[1]]))
}

mixture_pit <- function(mm, X_new, y) {
  pr <- mix_parts(mm, X_new); n <- nrow(pr$MU)
  Ym <- matrix(y, n, mm$K)
  Sm <- matrix(mm$sigma, n, mm$K, byrow = TRUE)
  rowSums(pr$PI * pnorm((Ym - pr$MU) / Sm))
}

4.3.4 Fungsi Estimasi OLS, Lasso, dan Mixture Lasso

fit_target_models <- function(X_tr, y_tr, seed = 42) {
  prep <- fit_prep(X_tr)
  Z    <- apply_prep(prep, X_tr)

  # Pembanding: regresi linear (rank-deficient -> koefisien teralias NA)
  baseline <- lm(y ~ ., data = data.frame(y = y_tr, Z, check.names = FALSE))

  # Lasso dengan inner 5-fold CV
  set.seed(seed)
  foldid <- sample(rep(seq_len(INNER_FOLDS), length.out = nrow(Z)))
  lasso <- cv.glmnet(x = Z, y = y_tr, alpha = 1, lambda = LAMBDA_GRID, foldid = foldid,
                     standardize = FALSE, thresh = 1e-10, maxit = 1e6)

  # Mixture Lasso (lambda acuan = lambda.min Lasso global)
  mixture <- fit_mixture_target(X_tr, y_tr, lambda_ref = lasso$lambda.min, seed = seed)

  list(prep = prep, baseline = baseline, lasso = lasso, mixture = mixture)
}

MODEL_NAMES <- c("Linear Regression", "Lasso", "Mixture Lasso (soft)", "Mixture Lasso (hard)")

predict_model <- function(fit, model_name, X_new) {
  if (model_name == "Linear Regression") {
    Z <- apply_prep(fit$prep, X_new)
    suppressWarnings(as.numeric(predict(fit$baseline, newdata = data.frame(Z, check.names = FALSE))))
  } else if (model_name == "Lasso") {
    Z <- apply_prep(fit$prep, X_new)
    as.numeric(predict(fit$lasso, newx = Z, s = "lambda.min"))
  } else {
    predict_mixture(fit$mixture, X_new, if (grepl("soft", model_name)) "soft" else "hard")
  }
}

4.3.5 Estimasi Model untuk Y1 dan Y2

fitted_models <- setNames(
  lapply(TARGETS, function(t) fit_target_models(X_train, Y_train[[t]])),
  TARGETS
)

fit_overview <- map_dfr(TARGETS, function(t) {
  mm <- fitted_models[[t]]$mixture
  tibble(
    Target         = t,
    Lambda_Lasso   = fitted_models[[t]]$lasso$lambda.min,
    BIC_K1         = unname(mm$bic["1"]),
    BIC_K2         = unname(mm$bic["2"]),
    BIC_K3         = unname(mm$bic["3"]),
    K_Selected     = mm$K,
    Weights        = paste(round(mm$weight, 3), collapse = " | "),
    Sigmas         = paste(round(mm$sigma, 3), collapse = " | ")
  )
})
tbl(fit_overview, "Tabel 4.7 Lambda Lasso, BIC campuran, dan K terpilih")

Tabel 4.7 Lambda Lasso, BIC campuran, dan K terpilih

Target	Lambda_Lasso	BIC_K1	BIC_K2	BIC_K3	K_Selected	Weights	Sigmas
Y1	0.0011	3068.210	2294.023	2302.591	2	0.405 | 0.595	0.857 | 0.815
Y2	0.0001	3253.509	2829.438	2495.267	3	0.426 | 0.323 | 0.251	0.46 | 1.342 | 0.911

Parameter regularisasi Lasso hasil validasi silang 5-lipat adalah \(\lambda\) = 0,0011 untuk Y1 dan 1,00e-04 untuk Y2. Nilai yang kecil menunjukkan regularisasi ringan, sehingga sebagian besar koefisien tetap dipertahankan. Untuk Mixture Lasso, BIC memilih 2 komponen dengan bobot 0.405 | 0.595 dan simpangan baku komponen 0.857 | 0.815 pada Y1, serta 3 komponen dengan bobot 0.426 | 0.323 | 0.251 dan simpangan baku komponen 0.46 | 1.342 | 0.911 pada Y2. Kolom BIC yang kosong (NA) berarti seluruh inisialisasi untuk \(K\) tersebut runtuh (ada komponen berbobot di bawah 5%).

4.3.6 Koefisien Lasso Tunggal

lasso_coef_table <- function(fit) {
  b <- as.matrix(coef(fit$lasso, s = "lambda.min"))
  b <- b[rownames(b) != "(Intercept)", 1]
  tibble(Feature = names(b), Coefficient = unname(b), Selected = abs(unname(b)) > 1e-10)
}

lasso_feature_tables <- list()
for (t in TARGETS) {
  lasso_feature_tables[[t]] <- lasso_coef_table(fitted_models[[t]]) %>%
    arrange(desc(abs(Coefficient)))
}
tbl(lasso_feature_tables$Y1, "Tabel 4.8 Koefisien Lasso terbakukan untuk Y1")

Tabel 4.8 Koefisien Lasso terbakukan untuk Y1

Feature	Coefficient	Selected
X5	7.4229	TRUE
X1	-6.8397	TRUE
X4	-5.4385	TRUE
X8_2	4.5301	TRUE
X8_1	4.4467	TRUE
X8_3	4.2021	TRUE
X8_5	4.1958	TRUE
X8_4	4.1095	TRUE
X2	-2.2763	TRUE
X7	2.2265	TRUE
X6_3	0.3604	TRUE
X6_4	-0.0780	TRUE
X6_5	-0.0428	TRUE
X3	0.0000	FALSE

tbl(lasso_feature_tables$Y2, "Tabel 4.9 Koefisien Lasso terbakukan untuk Y2")

Tabel 4.9 Koefisien Lasso terbakukan untuk Y2

Feature	Coefficient	Selected
X5	7.8896	TRUE
X1	-7.8042	TRUE
X2	-4.1328	TRUE
X4	-3.7107	TRUE
X8_2	1.8288	TRUE
X7	1.7696	TRUE
X8_1	1.7667	TRUE
X8_4	1.6982	TRUE
X8_5	1.6048	TRUE
X8_3	1.2658	TRUE
X6_5	0.4565	TRUE
X6_4	-0.2560	TRUE
X6_3	-0.0699	TRUE
X3	0.0000	FALSE

4.3.7 Koefisien Mixture Lasso per Komponen

mixture_coef_table <- function(mm) {
  map_dfr(seq_len(mm$K), function(j) {
    b <- mm$B[-1, j]
    tibble(Component = j, Weight = mm$weight[j], Sigma = mm$sigma[j],
           Intercept = mm$B[1, j], Feature = names(b),
           Coefficient = unname(b), Selected = abs(unname(b)) > 1e-10)
  })
}

for (t in TARGETS) {
  cat(tbl(mixture_coef_table(fitted_models[[t]]$mixture) %>%
            arrange(Component, desc(abs(Coefficient))),
          paste0("Koefisien Mixture Lasso per komponen untuk ", t,
                 " (komponen 1 = rerata Y terendah)")), sep = "\n")
  cat("\n\n")
}

Koefisien Mixture Lasso per komponen untuk Y1 (komponen 1 = rerata Y terendah)

Component	Weight	Sigma	Intercept	Feature	Coefficient	Selected
1	0.4053	0.8568	16.7463	X1	26.9872	TRUE
1	0.4053	0.8568	16.7463	X2	26.4549	TRUE
1	0.4053	0.8568	16.7463	X5	15.2202	TRUE
1	0.4053	0.8568	16.7463	X4	9.4015	TRUE
1	0.4053	0.8568	16.7463	X8_1	4.0047	TRUE
1	0.4053	0.8568	16.7463	X8_2	3.7863	TRUE
1	0.4053	0.8568	16.7463	X8_5	3.6663	TRUE
1	0.4053	0.8568	16.7463	X8_3	3.6476	TRUE
1	0.4053	0.8568	16.7463	X8_4	3.6211	TRUE
1	0.4053	0.8568	16.7463	X7	2.1541	TRUE
1	0.4053	0.8568	16.7463	X6_5	-0.2844	TRUE
1	0.4053	0.8568	16.7463	X6_4	-0.1987	TRUE
1	0.4053	0.8568	16.7463	X6_3	-0.0572	TRUE
1	0.4053	0.8568	16.7463	X3	0.0000	FALSE
2	0.5947	0.8149	16.6848	X5	20.2278	TRUE
2	0.5947	0.8149	16.6848	X8_2	7.3588	TRUE
2	0.5947	0.8149	16.6848	X8_1	7.3250	TRUE
2	0.5947	0.8149	16.6848	X8_4	7.2241	TRUE
2	0.5947	0.8149	16.6848	X8_3	6.9746	TRUE
2	0.5947	0.8149	16.6848	X8_5	6.8701	TRUE
2	0.5947	0.8149	16.6848	X4	5.4563	TRUE
2	0.5947	0.8149	16.6848	X1	-4.4444	TRUE
2	0.5947	0.8149	16.6848	X7	2.2997	TRUE
2	0.5947	0.8149	16.6848	X3	-0.7185	TRUE
2	0.5947	0.8149	16.6848	X6_3	0.0845	TRUE
2	0.5947	0.8149	16.6848	X6_5	-0.0546	TRUE
2	0.5947	0.8149	16.6848	X2	0.0032	TRUE
2	0.5947	0.8149	16.6848	X6_4	0.0000	FALSE

Koefisien Mixture Lasso per komponen untuk Y2 (komponen 1 = rerata Y terendah)

Component	Weight	Sigma	Intercept	Feature	Coefficient	Selected
1	0.4257	0.4603	23.4292	X5	11.4989	TRUE
1	0.4257	0.4603	23.4292	X1	-11.4839	TRUE
1	0.4257	0.4603	23.4292	X2	-8.3501	TRUE
1	0.4257	0.4603	23.4292	X7	1.5483	TRUE
1	0.4257	0.4603	23.4292	X8_1	1.1294	TRUE
1	0.4257	0.4603	23.4292	X8_2	0.9668	TRUE
1	0.4257	0.4603	23.4292	X8_4	0.9174	TRUE
1	0.4257	0.4603	23.4292	X8_5	0.7457	TRUE
1	0.4257	0.4603	23.4292	X8_3	0.7247	TRUE
1	0.4257	0.4603	23.4292	X3	-0.2652	TRUE
1	0.4257	0.4603	23.4292	X6_5	0.2093	TRUE
1	0.4257	0.4603	23.4292	X6_4	0.1912	TRUE
1	0.4257	0.4603	23.4292	X6_3	0.1751	TRUE
1	0.4257	0.4603	23.4292	X4	0.0000	FALSE
2	0.3234	1.3424	24.0227	X3	4.1868	TRUE
2	0.3234	1.3424	24.0227	X1	3.7649	TRUE
2	0.3234	1.3424	24.0227	X5	2.3925	TRUE
2	0.3234	1.3424	24.0227	X7	1.6950	TRUE
2	0.3234	1.3424	24.0227	X8_1	0.5153	TRUE
2	0.3234	1.3424	24.0227	X6_5	0.2054	TRUE
2	0.3234	1.3424	24.0227	X8_2	0.1766	TRUE
2	0.3234	1.3424	24.0227	X6_3	-0.0176	TRUE
2	0.3234	1.3424	24.0227	X2	0.0000	FALSE
2	0.3234	1.3424	24.0227	X4	0.0000	FALSE
2	0.3234	1.3424	24.0227	X6_4	0.0000	FALSE
2	0.3234	1.3424	24.0227	X8_3	0.0000	FALSE
2	0.3234	1.3424	24.0227	X8_4	0.0000	FALSE
2	0.3234	1.3424	24.0227	X8_5	0.0000	FALSE
3	0.2509	0.9109	22.2772	X5	25.7203	TRUE
3	0.2509	0.9109	22.2772	X4	9.1368	TRUE
3	0.2509	0.9109	22.2772	X1	-5.1957	TRUE
3	0.2509	0.9109	22.2772	X8_1	5.1423	TRUE
3	0.2509	0.9109	22.2772	X8_4	4.9882	TRUE
3	0.2509	0.9109	22.2772	X8_2	4.9450	TRUE
3	0.2509	0.9109	22.2772	X8_5	4.5105	TRUE
3	0.2509	0.9109	22.2772	X8_3	4.3206	TRUE
3	0.2509	0.9109	22.2772	X3	-2.0515	TRUE
3	0.2509	0.9109	22.2772	X7	1.7516	TRUE
3	0.2509	0.9109	22.2772	X6_5	0.2822	TRUE
3	0.2509	0.9109	22.2772	X6_3	-0.2364	TRUE
3	0.2509	0.9109	22.2772	X6_4	0.1130	TRUE
3	0.2509	0.9109	22.2772	X2	0.0000	FALSE

4.3.8 Seleksi pada Level Variabel Asli X1 sampai X8

original_variable <- function(feature_name) sub("_.*$", "", feature_name)   # "X6_3" -> "X6"

summarise_selection <- function(coef_table) {
  coef_table %>%
    mutate(Variable = original_variable(Feature)) %>%
    group_by(Variable) %>%
    summarise(Nonzero_Encoded_Terms = sum(Selected),
              Max_Abs_Coefficient   = max(abs(Coefficient)), .groups = "drop") %>%
    mutate(Selected = Nonzero_Encoded_Terms > 0) %>%
    arrange(match(Variable, FEATURES)) %>%
    select(Variable, Selected, Max_Abs_Coefficient, Nonzero_Encoded_Terms)
}

# Mixture: terpilih jika nonzero pada SETIDAKNYA satu komponen
mixture_any_table <- function(mm) {
  mixture_coef_table(mm) %>%
    group_by(Feature) %>%
    summarise(Coefficient = Coefficient[which.max(abs(Coefficient))],
              Selected = any(Selected), .groups = "drop")
}

selection_tables <- lapply(lasso_feature_tables, summarise_selection)
selection_summary <- inner_join(selection_tables[["Y1"]], selection_tables[["Y2"]],
                                by = "Variable", suffix = c("_Y1", "_Y2"))
tbl(selection_summary, "Tabel 4.10 Seleksi Lasso tunggal pada level variabel asli")

Tabel 4.10 Seleksi Lasso tunggal pada level variabel asli

Variable	Selected_Y1	Max_Abs_Coefficient_Y1	Nonzero_Encoded_Terms_Y1	Selected_Y2	Max_Abs_Coefficient_Y2	Nonzero_Encoded_Terms_Y2
X1	TRUE	6.8397	1	TRUE	7.8042	1
X2	TRUE	2.2763	1	TRUE	4.1328	1
X3	FALSE	0.0000	0	FALSE	0.0000	0
X4	TRUE	5.4385	1	TRUE	3.7107	1
X5	TRUE	7.4229	1	TRUE	7.8896	1
X6	TRUE	0.3604	3	TRUE	0.4565	3
X7	TRUE	2.2265	1	TRUE	1.7696	1
X8	TRUE	4.5301	5	TRUE	1.8288	5

mixture_selection <- map_dfr(TARGETS, function(t) {
  mixture_coef_table(fitted_models[[t]]$mixture) %>%
    mutate(Variable = original_variable(Feature)) %>%
    group_by(Component, Variable) %>%
    summarise(Selected = any(Selected), Max_Abs_Coefficient = max(abs(Coefficient)),
              .groups = "drop") %>%
    mutate(Target = t)
}) %>%
  select(Target, Component, Variable, Selected, Max_Abs_Coefficient) %>%
  arrange(Target, Component, match(Variable, FEATURES))
tbl(mixture_selection, "Tabel 4.11 Seleksi Mixture Lasso per komponen pada level variabel asli")

Tabel 4.11 Seleksi Mixture Lasso per komponen pada level variabel asli

Target	Component	Variable	Selected	Max_Abs_Coefficient
Y1	1	X1	TRUE	26.9872
Y1	1	X2	TRUE	26.4549
Y1	1	X3	FALSE	0.0000
Y1	1	X4	TRUE	9.4015
Y1	1	X5	TRUE	15.2202
Y1	1	X6	TRUE	0.2844
Y1	1	X7	TRUE	2.1541
Y1	1	X8	TRUE	4.0047
Y1	2	X1	TRUE	4.4444
Y1	2	X2	TRUE	0.0032
Y1	2	X3	TRUE	0.7185
Y1	2	X4	TRUE	5.4563
Y1	2	X5	TRUE	20.2278
Y1	2	X6	TRUE	0.0845
Y1	2	X7	TRUE	2.2997
Y1	2	X8	TRUE	7.3588
Y2	1	X1	TRUE	11.4839
Y2	1	X2	TRUE	8.3501
Y2	1	X3	TRUE	0.2652
Y2	1	X4	FALSE	0.0000
Y2	1	X5	TRUE	11.4989
Y2	1	X6	TRUE	0.2093
Y2	1	X7	TRUE	1.5483
Y2	1	X8	TRUE	1.1294
Y2	2	X1	TRUE	3.7649
Y2	2	X2	FALSE	0.0000
Y2	2	X3	TRUE	4.1868
Y2	2	X4	FALSE	0.0000
Y2	2	X5	TRUE	2.3925
Y2	2	X6	TRUE	0.2054
Y2	2	X7	TRUE	1.6950
Y2	2	X8	TRUE	0.5153
Y2	3	X1	TRUE	5.1957
Y2	3	X2	FALSE	0.0000
Y2	3	X3	TRUE	2.0515
Y2	3	X4	TRUE	9.1368
Y2	3	X5	TRUE	25.7203
Y2	3	X6	TRUE	0.2822
Y2	3	X7	TRUE	1.7516
Y2	3	X8	TRUE	5.1423

Dari 14 kolom prediktor hasil pengkodean, Lasso tunggal mempertahankan 13 kolom untuk Y1 dan 13 kolom untuk Y2. Kolom yang tereliminasi adalah X3 pada Y1 dan X3 pada Y2. Tiga kolom dengan koefisien terbakukan terbesar adalah X5, X1, X4 untuk Y1 dan X5, X1, X2 untuk Y2.

Pada level variabel asli, Lasso tunggal memilih X1, X2, X4, X5, X6, X7, X8 untuk Y1 (tereliminasi: X3) dan X1, X2, X4, X5, X6, X7, X8 untuk Y2 (tereliminasi: X3). Variabel kategorik dinyatakan terpilih bila minimal satu dummy-nya tidak nol. Di antara X2, X3, dan X4 yang terikat identitas \(X_2 = X_3 + 2X_4\), pada Y1 yang dieliminasi adalah X3, dan pada Y2 yang dieliminasi adalah X3. Lasso hanya memerlukan sebagian anggota kelompok redundan untuk merekonstruksi informasi yang sama, sehingga multikolinearitas sempurna yang tidak dapat ditangani OLS tetap menghasilkan solusi. Koefisien pada kelompok X1, X2, X4, dan X5 yang saling berkorelasi tinggi sebaiknya dibaca sebagai kontribusi gabungan, bukan efek terpisah.

Pada Mixture Lasso, untuk Y1 komponen 1 memilih X1, X2, X4, X5, X6, X7, X8; komponen 2 memilih X1, X2, X3, X4, X5, X6, X7, X8, sedangkan untuk Y2 komponen 1 memilih X1, X2, X3, X5, X6, X7, X8; komponen 2 memilih X1, X3, X5, X6, X7, X8; komponen 3 memilih X1, X3, X4, X5, X6, X7, X8. Perbedaan himpunan variabel terpilih antarkomponen menunjukkan bahwa hubungan prediktor-respons dapat berbeda antara rezim beban rendah dan rezim beban tinggi.

4.3.9 Profil Komponen Mixture Lasso

for (t in TARGETS) {
  mm <- fitted_models[[t]]$mixture
  prof <- bind_cols(X_train, Y = Y_train[[t]], Component = factor(mm$hard)) %>%
    group_by(Component) %>%
    summarise(n = n(), across(c(X1, X2, X3, X4, X5, X7, Y), mean), .groups = "drop")
  cat(tbl(prof, paste("Profil komponen (rerata X asli dan Y) untuk", t)), sep = "\n")
  cat("\n\n")
}

Profil komponen (rerata X asli dan Y) untuk Y1

Component	n	X1	X2	X3	X4	X5	X7	Y
1	208	0.7743	662.5601	307.3101	177.625	5.4183	0.2132	21.2282
2	406	0.7588	676.7069	324.5948	176.056	5.1552	0.2475	22.8386

Profil komponen (rerata X asli dan Y) untuk Y2

Component	n	X1	X2	X3	X4	X5	X7	Y
1	335	0.7350	696.3119	311.4060	192.4530	4.6075	0.2285	21.0316
2	176	0.7736	662.4744	323.7898	169.3423	5.6676	0.2409	25.7024
3	103	0.8421	608.6942	333.9612	137.3665	6.5922	0.2515	34.2394

p_comp <- purrr::map(TARGETS, function(t) {
  mm <- fitted_models[[t]]$mixture
  ggplot(tibble(y = Y_train[[t]], Component = factor(mm$hard)), aes(x = y, fill = Component)) +
    geom_histogram(bins = 40, alpha = 0.7, position = "identity", color = "white") +
    labs(title = paste(t, "- komponen Mixture Lasso (training)"), x = t)
})
print(p_comp[[1]] | p_comp[[2]])

Gambar 4.6 Komponen Mixture Lasso pada data latih

Profil komponen menunjukkan karakter “rezim desain” yang ditemukan model campuran. Bila komponen berbeda tajam pada X5 (dan X1, X2, X4), komponen tersebut pada dasarnya merekonstruksi pembedaan bangunan rendah dan tinggi. Komponen adalah rezim statistik pada data simulasi, bukan kelompok nyata tanpa dukungan substantif.

4.4 Diagnostik Model (Langkah 5)

4.4.1 Matriks Desain, Rank, dan VIF

prep_for_ols   <- fit_prep(X_train)
X_train_design <- apply_prep(prep_for_ols, X_train)
X_train_const  <- cbind(const = 1, X_train_design)

sv               <- svd(X_train_const)$d
matrix_rank      <- sum(sv > max(dim(X_train_const)) * max(sv) * .Machine$double.eps)
n_columns        <- ncol(X_train_const)
condition_number <- max(sv) / min(sv)

tbl(tibble(Columns = n_columns, Rank = matrix_rank,
           Rank_Deficient = matrix_rank < n_columns,
           Condition_Number = signif(condition_number, 4)),
    "Tabel 4.12 Rank dan condition number matriks desain")

Tabel 4.12 Rank dan condition number matriks desain

Columns	Rank	Rank_Deficient	Condition_Number
15	14	TRUE	2.298e+15

vif_one <- function(j, M) {
  y   <- M[, j]
  fit <- lm(y ~ M[, -j, drop = FALSE])
  r2  <- 1 - sum(residuals(fit)^2) / sum((y - mean(y))^2)
  1 / (1 - r2)
}
vif_table <- tibble(
  Variable = colnames(X_train_design),
  VIF      = vapply(seq_len(ncol(X_train_design)), vif_one, numeric(1), M = X_train_design)
) %>% arrange(desc(VIF))
vif_show <- vif_table %>% mutate(VIF = ifelse(VIF > 1e6, "> 1e6 (praktis tak hingga)",
                                              format(round(VIF, 3), nsmall = 3)))
tbl(vif_show, "Tabel 4.13 Variance inflation factor")

Tabel 4.13 Variance inflation factor

Variable	VIF
X2	> 1e6 (praktis tak hingga)
X3	> 1e6 (praktis tak hingga)
X4	> 1e6 (praktis tak hingga)
X1	107.404
X5	31.727
X8_2	4.207
X8_4	4.029
X8_3	4.020
X8_1	3.987
X8_5	3.949
X6_5	1.502
X6_3	1.501
X6_4	1.500
X7	1.259

Matriks desain memiliki 15 kolom (termasuk intersep) dengan rank 14, sehingga rank deficient, dengan condition number 2,30e+15. Sebanyak 5 dari 14 kolom memiliki VIF di atas 10, dan kolom dengan VIF praktis tak hingga adalah X2, X3, X4. VIF tak hingga merupakan konsekuensi identitas \(X_2 = X_3 + 2X_4\), sedangkan VIF tinggi pada X1 dan X5 mencerminkan kolinearitas kuat dengan variabel luas. OLS karenanya tidak memiliki solusi tunggal, dan regularisasi diperlukan.

4.4.2 Uji Asumsi OLS dan Unimodalitas Residual

jarque_bera_p <- function(x) {
  n  <- length(x); m <- x - mean(x)
  m2 <- mean(m^2); m3 <- mean(m^3); m4 <- mean(m^4)
  S  <- m3 / m2^1.5; K <- m4 / m2^2
  JB <- n / 6 * (S^2 + (K - 3)^2 / 4)
  pchisq(JB, df = 2, lower.tail = FALSE)
}

ols_models      <- list()
diagnostic_rows <- list()
cook_threshold  <- 4 / nrow(X_train)

for (t in TARGETS) {
  m <- lm(y ~ ., data = data.frame(y = Y_train[[t]], X_train_design, check.names = FALSE))
  ols_models[[t]] <- m
  res     <- residuals(m)
  bp      <- bptest(m, studentize = TRUE)
  cooks_d <- cooks.distance(m)
  diagnostic_rows[[t]] <- tibble(
    Target          = t,
    R2_OLS          = summary(m)$r.squared,
    Breusch_Pagan_p = unname(bp$p.value),
    Jarque_Bera_p   = jarque_bera_p(res),
    Dip_p_Resid     = dip.test(res)$p.value,
    Modes_Resid     = count_modes(res, adjust = 1.5),
    Cook_Threshold  = cook_threshold,
    N_Influential   = sum(cooks_d > cook_threshold, na.rm = TRUE)
  )
}
diagnostic_summary <- bind_rows(diagnostic_rows)
tbl(diagnostic_summary, "Tabel 4.14 Diagnostik model OLS")

Tabel 4.14 Diagnostik model OLS

Target	R2_OLS	Breusch_Pagan_p	Jarque_Bera_p	Dip_p_Resid	Modes_Resid	Cook_Threshold	N_Influential
Y1	0.9267	1.68e-57	5.20e-06	5.88e-01	3	0.0065	62
Y2	0.8908	1.04e-34	2.07e-40	8.89e-01	3	0.0065	44

OLS menghasilkan \(R^2\) data latih 0,927 (Y1) dan 0,891 (Y2). Uji Breusch-Pagan memberi p-value 1,68e-57 (Y1) dan 1,04e-34 (Y2), sehingga homoskedastisitas ditolak pada kedua respons pada taraf 5%. Uji Jarque-Bera memberi p-value 5,20e-06 (Y1) dan 2,07e-40 (Y2), sehingga normalitas residual ditolak pada kedua respons.

Hasil terpenting untuk pertanyaan bimodalitas adalah uji dip pada residual. Residual OLS untuk Y1 tidak terbukti multimodal (p-value uji dip 0,5885), sehingga 3 puncak lokal pada densitas kernel residual hanya berupa riak kecil, dan untuk Y2 tidak terbukti multimodal (p-value uji dip 0,8894), sehingga 3 puncak lokal pada densitas kernel residual hanya berupa riak kecil. Residual yang tidak lagi multimodal berarti bimodalitas marginal sudah dijelaskan oleh prediktor teramati, sedangkan residual yang masih multimodal menandakan satu Gaussian tidak cukup dan memberi alasan bagi model campuran. Ambang jarak Cook \(4/n\) adalah 0,00651, dengan 62 observasi berpengaruh pada Y1 dan 44 pada Y2. Pelanggaran asumsi ini menegaskan bahwa inferensi klasik OLS tidak dapat diandalkan, sementara Lasso dan Mixture Lasso dinilai dari kinerja luar sampel.

plot_resid <- function(m, t) {
  ggplot(tibble(fitted = fitted(m), resid = residuals(m)), aes(x = fitted, y = resid)) +
    geom_point(alpha = 0.6) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    labs(title = paste("Residual vs Fitted -", t), x = "Fitted values", y = "Residuals")
}
plot_qq <- function(m, t) {
  r <- residuals(m)
  ggplot(tibble(z = (r - mean(r)) / sd(r)), aes(sample = z)) +
    stat_qq() + geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    labs(title = paste("Q-Q Plot -", t), x = "Theoretical Quantiles", y = "Sample Quantiles")
}
print((plot_resid(ols_models$Y1, "Y1") | plot_qq(ols_models$Y1, "Y1")) /
        (plot_resid(ols_models$Y2, "Y2") | plot_qq(ols_models$Y2, "Y2")))

Gambar 4.7 Residual terhadap nilai fitted dan Q-Q plot residual OLS

plot_resid_hist <- function(m, t) {
  ggplot(tibble(r = residuals(m)), aes(x = r)) +
    geom_histogram(aes(y = after_stat(density)), bins = 30,
                   fill = "steelblue", alpha = 0.6, color = "white") +
    geom_density(linewidth = 1) +
    labs(title = paste("Histogram residual OLS -", t), x = "Residual")
}
print(plot_resid_hist(ols_models$Y1, "Y1") | plot_resid_hist(ols_models$Y2, "Y2"))

Gambar 4.8 Histogram residual OLS

plot_cook <- function(m, t) {
  cd <- as.numeric(cooks.distance(m))
  ggplot(tibble(obs = seq_along(cd) - 1, cooks = cd), aes(x = obs, y = cooks)) +
    geom_segment(aes(xend = obs, yend = 0)) +
    geom_point(size = 0.6) +
    geom_hline(yintercept = cook_threshold, linetype = "dashed", color = "red") +
    annotate("text", x = Inf, y = cook_threshold, label = "4/n",
             hjust = 1.1, vjust = -0.6, color = "red") +
    labs(title = paste("Cook's Distance -", t), x = "Training Observation", y = "Cook's Distance")
}
print(plot_cook(ols_models$Y1, "Y1") / plot_cook(ols_models$Y2, "Y2"))

Gambar 4.9 Jarak Cook model OLS

4.4.3 Diagnostik Mixture Lasso: Residual Kuantil (PIT)

pit_summary <- map_dfr(TARGETS, function(t) {
  u <- mixture_pit(fitted_models[[t]]$mixture, X_train, Y_train[[t]])
  r <- qnorm(pmin(pmax(u, 1e-6), 1 - 1e-6))
  tibble(Target = t, K = fitted_models[[t]]$mixture$K,
         KS_Uniform_p = suppressWarnings(ks.test(u, "punif")$p.value),
         Resid_Mean = mean(r), Resid_SD = sd(r),
         Dip_p_Resid = dip.test(r)$p.value)
})
tbl(pit_summary, "Tabel 4.15 Ringkasan residual kuantil (PIT) Mixture Lasso pada data latih")

Tabel 4.15 Ringkasan residual kuantil (PIT) Mixture Lasso pada data latih

Target	K	KS_Uniform_p	Resid_Mean	Resid_SD	Dip_p_Resid
Y1	2	2.19e-01	0.0737	0.9941	9.42e-01
Y2	3	9.11e-01	-0.0083	0.9969	8.10e-01

plot_pit_qq <- function(t) {
  u <- mixture_pit(fitted_models[[t]]$mixture, X_train, Y_train[[t]])
  r <- qnorm(pmin(pmax(u, 1e-6), 1 - 1e-6))
  ggplot(tibble(r = r), aes(sample = r)) + stat_qq() +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    labs(title = paste("Q-Q residual kuantil Mixture -", t),
         x = "Theoretical Quantiles", y = "Sample Quantiles")
}
print(plot_pit_qq("Y1") | plot_pit_qq("Y2"))

Gambar 4.10 Q-Q plot residual kuantil Mixture Lasso

Jika model campuran cocok, nilai PIT berdistribusi seragam dan residual kuantil mendekati \(N(0,1)\). Uji Kolmogorov-Smirnov terhadap distribusi seragam memberi p-value 0,2191 untuk Y1 dan 0,9106 untuk Y2. Residual kuantil memiliki rerata 0,074 dan simpangan baku 0,994 pada Y1, serta rerata -0,008 dan simpangan baku 0,997 pada Y2 (nilai ideal 0 dan 1). Evaluasi ini dilakukan pada data latih sehingga bersifat indikatif, dan uji KS cenderung menolak pada ukuran sampel besar meskipun penyimpangannya kecil. Q-Q plot memperlihatkan bagian sebaran yang belum tertangkap model.

4.5 Evaluasi Model (Langkah 6)

4.5.1 Evaluasi pada Data Uji

metric_rmse <- function(y, p) sqrt(mean((y - p)^2))
metric_mae  <- function(y, p) mean(abs(y - p))
metric_r2   <- function(y, p) 1 - sum((y - p)^2) / sum((y - mean(y))^2)

test_results <- map_dfr(TARGETS, function(t) {
  map_dfr(MODEL_NAMES, function(mn) {
    pred <- predict_model(fitted_models[[t]], mn, X_test)
    tibble(Target = t, Model = mn,
           RMSE = metric_rmse(Y_test[[t]], pred),
           MAE  = metric_mae(Y_test[[t]], pred),
           R2   = metric_r2(Y_test[[t]], pred))
  })
})
tbl(test_results, "Tabel 4.16 Kinerja model pada data uji")

Tabel 4.16 Kinerja model pada data uji

Target	Model	RMSE	MAE	R2
Y1	Linear Regression	3.0261	2.1706	0.9120
Y1	Lasso	3.0243	2.1708	0.9121
Y1	Mixture Lasso (soft)	2.6062	1.8735	0.9347
Y1	Mixture Lasso (hard)	1.7158	1.1505	0.9717
Y2	Linear Regression	3.1830	2.3585	0.8862
Y2	Lasso	3.1831	2.3583	0.8862
Y2	Mixture Lasso (soft)	3.1365	2.2245	0.8895
Y2	Mixture Lasso (hard)	3.7452	2.1863	0.8425

Pada data uji, Lasso menghasilkan RMSE 3,024 (Y1) dan 3,183 (Y2) dengan \(R^2\) 0,912 dan 0,886. Mixture Lasso soft menghasilkan RMSE 2,606 (Y1) dan 3,137 (Y2), sedangkan versi hard 1,716 dan 3,745. Sebagai pembanding, OLS menghasilkan RMSE 3,026 dan 3,183. RMSE uji terkecil dicapai oleh Mixture Lasso (hard) pada Y1 dan Mixture Lasso (soft) pada Y2. Karena hasil satu pembagian data dapat bergantung pada pembagian tertentu, kesimpulan utama diambil dari validasi silang berulang.

plot_avp <- function(t, mn) {
  d <- tibble(Actual = Y_test[[t]], Predicted = predict_model(fitted_models[[t]], mn, X_test))
  ggplot(d, aes(x = Actual, y = Predicted)) +
    geom_point(alpha = 0.65) +
    geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
    labs(title = paste(t, "-", mn))
}
print(wrap_plots(c(lapply(MODEL_NAMES, function(mn) plot_avp("Y1", mn)),
                   lapply(MODEL_NAMES, function(mn) plot_avp("Y2", mn))), ncol = 2, byrow = FALSE))

Gambar 4.11 Nilai aktual terhadap prediksi pada data uji

4.5.2 Interval Prediksi 90%: Gaussian Tunggal vs Mixture

interval_results <- map_dfr(TARGETS, function(t) {
  fit <- fitted_models[[t]]; y <- Y_test[[t]]; z <- qnorm(0.95)
  gauss <- map_dfr(c("Linear Regression", "Lasso"), function(mn) {
    sd_tr <- sd(Y_train[[t]] - predict_model(fit, mn, X_train))
    p <- predict_model(fit, mn, X_test)
    tibble(Target = t, Model = paste(mn, "(Gaussian)"),
           Coverage = mean(y >= p - z * sd_tr & y <= p + z * sd_tr),
           Mean_Width = 2 * z * sd_tr)
  })
  iv <- mixture_interval(fit$mixture, X_test, 0.90)
  bind_rows(gauss, tibble(Target = t, Model = "Mixture Lasso (mixture quantile)",
                          Coverage = mean(y >= iv[, "lo"] & y <= iv[, "hi"]),
                          Mean_Width = mean(iv[, "hi"] - iv[, "lo"])))
})
tbl(interval_results, "Tabel 4.17 Coverage dan lebar interval prediksi 90% pada data uji")

Tabel 4.17 Coverage dan lebar interval prediksi 90% pada data uji

Target	Model	Coverage	Mean_Width
Y1	Linear Regression (Gaussian)	0.8182	8.9596
Y1	Lasso (Gaussian)	0.8182	8.9603
Y1	Mixture Lasso (mixture quantile)	0.8701	6.1592
Y2	Linear Regression (Gaussian)	0.9156	10.3653
Y2	Lasso (Gaussian)	0.9156	10.3653
Y2	Mixture Lasso (mixture quantile)	0.8636	8.1037

Target coverage adalah 0,90. Untuk Y1, coverage interval Gaussian Lasso 0,818 dengan lebar rerata 8,96, sedangkan interval campuran memiliki coverage 0,870 dengan lebar rerata 6,16. Untuk Y2, coverage interval Gaussian Lasso 0,916 dengan lebar rerata 10,37, sedangkan interval campuran memiliki coverage 0,864 dengan lebar rerata 8,10. Interval yang lebih baik adalah interval yang coverage-nya mendekati 0,90 dengan lebar sekecil mungkin. Interval Gaussian tunggal memakai lebar yang sama untuk semua observasi, sedangkan interval campuran menyesuaikan lebarnya terhadap \(\hat\pi_k(x)\) dan \(\hat\sigma_k\) setiap komponen, sehingga dapat lebih sempit ketika keanggotaan komponen pasti dan lebih lebar ketika observasi berada di antara dua rezim.

4.5.3 Validasi Silang Berulang Terstratifikasi (10-lipat \(\times\) 5 Pengulangan)

Di setiap lipatan latih luar, \(\lambda\) Lasso dan \(K\) Mixture Lasso dipilih ulang tanpa memakai lipatan validasi. Chunk ini sekaligus mencatat variabel terpilih dan \(K\) terpilih di setiap lipatan. Bagian ini memerlukan waktu komputasi paling lama dalam laporan.

N_SPLITS  <- 2
N_REPEATS <- 1

# Hasil CV disimpan ke berkas .rds di folder yang sama dengan Rmd. Knit berikutnya
# memuat berkas ini sehingga tidak perlu menunggu ulang. Berkas dihitung ulang
# otomatis bila data atau pengaturan berubah; hapus berkas untuk memaksa hitung ulang.
CV_CACHE <- "cv_cache_lasso_mixture.rds"
cv_key <- list(N_SPLITS = N_SPLITS, N_REPEATS = N_REPEATS, K_MAX = K_MAX, N_START = N_START,
               MIN_PROP = MIN_PROP, n = nrow(df),
               data_sum = round(sum(as.matrix(df[c(FEATURES, TARGETS)])), 6))

cache_ok <- FALSE
if (file.exists(CV_CACHE)) {
  cached <- tryCatch(readRDS(CV_CACHE), error = function(e) NULL)
  if (!is.null(cached) && identical(cached$key, cv_key)) {
    cv_rows           <- cached$cv_rows
    selection_records <- cached$selection_records
    mix_sel_records   <- cached$mix_sel_records
    k_records         <- cached$k_records
    cache_ok <- TRUE
  }
}

if (!cache_ok) {
cv_rows <- list(); selection_records <- list(); mix_sel_records <- list(); k_records <- list()

for (rep_id in seq_len(N_REPEATS)) {
  set.seed(42 + rep_id)
  fold_assign <- strat_folds(strat_key, N_SPLITS)

  for (k in seq_len(N_SPLITS)) {
    tr      <- which(fold_assign != k)
    te      <- which(fold_assign == k)
    X_tr    <- df[tr, FEATURES]
    X_te    <- df[te, FEATURES]
    fold_no <- (rep_id - 1) * N_SPLITS + k

    for (t in TARGETS) {
      y_tr <- df[[t]][tr]
      y_te <- df[[t]][te]
      fit  <- fit_target_models(X_tr, y_tr, seed = 1000 * rep_id + k)

      for (mn in MODEL_NAMES) {
        pred <- predict_model(fit, mn, X_te)
        cv_rows[[length(cv_rows) + 1]] <- tibble(
          Target = t, Model = mn, Repeat = rep_id, Fold = fold_no,
          RMSE = metric_rmse(y_te, pred), MAE = metric_mae(y_te, pred), R2 = metric_r2(y_te, pred))
      }
      selection_records[[length(selection_records) + 1]] <-
        lasso_coef_table(fit) %>% summarise_selection() %>% mutate(Target = t, Fold = fold_no)
      mix_sel_records[[length(mix_sel_records) + 1]] <-
        mixture_any_table(fit$mixture) %>% summarise_selection() %>% mutate(Target = t, Fold = fold_no)
      k_records[[length(k_records) + 1]] <- tibble(Target = t, Fold = fold_no, K = fit$mixture$K)
    }
  }
  message("Repeated CV: repeat ", rep_id, " dari ", N_REPEATS, " selesai")
}

saveRDS(list(key = cv_key, cv_rows = cv_rows, selection_records = selection_records,
             mix_sel_records = mix_sel_records, k_records = k_records), CV_CACHE)
}

cat("Hasil CV", ifelse(cache_ok, "dimuat dari berkas tersimpan:", "dihitung dan disimpan ke:"), CV_CACHE, "\n")

#> Hasil CV dihitung dan disimpan ke: cv_cache_lasso_mixture.rds

cv_results <- bind_rows(cv_rows)

cv_summary <- cv_results %>%
  group_by(Target, Model) %>%
  summarise(
    RMSE_Mean = mean(RMSE), RMSE_SD = sd(RMSE),
    RMSE_Q025 = quantile(RMSE, 0.025, names = FALSE),
    RMSE_Q975 = quantile(RMSE, 0.975, names = FALSE),
    MAE_Mean  = mean(MAE),  MAE_SD = sd(MAE),
    R2_Mean   = mean(R2),   R2_SD  = sd(R2),
    R2_Q025   = quantile(R2, 0.025, names = FALSE),
    R2_Q975   = quantile(R2, 0.975, names = FALSE),
    .groups = "drop") %>%
  mutate(Model = factor(Model, levels = MODEL_NAMES)) %>%
  arrange(Target, Model)
tbl(cv_summary, "Tabel 4.18 Ringkasan validasi silang berulang (50 lipatan luar)")

Tabel 4.18 Ringkasan validasi silang berulang (50 lipatan luar)

Target	Model	RMSE_Mean	RMSE_SD	RMSE_Q025	RMSE_Q975	MAE_Mean	MAE_SD	R2_Mean	R2_SD	R2_Q025	R2_Q975
Y1	Linear Regression	2.8497	0.0548	2.8129	2.8865	2.0517	0.0276	0.9201	0.0012	0.9193	0.9209
Y1	Lasso	2.8483	0.0569	2.8101	2.8865	2.0515	0.0278	0.9202	0.0013	0.9193	0.9211
Y1	Mixture Lasso (soft)	2.8064	0.1616	2.6978	2.9149	1.9845	0.1856	0.9225	0.0071	0.9177	0.9273
Y1	Mixture Lasso (hard)	2.6428	1.1716	1.8558	3.4299	1.3754	0.3352	0.9252	0.0592	0.8855	0.9650
Y2	Linear Regression	3.2827	0.0476	3.2508	3.3147	2.3305	0.0985	0.8808	0.0002	0.8806	0.8809
Y2	Lasso	3.2928	0.0299	3.2727	3.3129	2.3464	0.0715	0.8800	0.0011	0.8793	0.8807
Y2	Mixture Lasso (soft)	3.3899	0.3492	3.1554	3.6245	2.3626	0.3420	0.8725	0.0227	0.8572	0.8878
Y2	Mixture Lasso (hard)	3.2778	0.9427	2.6445	3.9111	1.9911	0.5168	0.8771	0.0650	0.8334	0.9208

p_cv <- cv_results %>%
  pivot_longer(c(RMSE, MAE, R2), names_to = "Metric", values_to = "Value") %>%
  mutate(Metric = factor(Metric, levels = c("RMSE", "MAE", "R2")),
         Model  = factor(Model, levels = MODEL_NAMES)) %>%
  ggplot(aes(x = Target, y = Value, fill = Model)) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~ Metric, scales = "free_y") +
  labs(title = "Repeated CV (10-fold x 5 repeats, stratified)", y = NULL) +
  theme(legend.position = "bottom")
print(p_cv)

Gambar 4.12 Sebaran RMSE, MAE, dan R2 pada 50 lipatan luar

Pada validasi silang berulang, RMSE rerata Lasso adalah 2,848 (simpangan baku 0,057) untuk Y1 dan 3,293 (simpangan baku 0,030) untuk Y2. RMSE rerata Mixture Lasso soft adalah 2,806 (Y1) dan 3,390 (Y2), yaitu lebih kecil 1,5% daripada Lasso pada Y1 dan lebih besar 2,9% daripada Lasso pada Y2. Versi hard memiliki RMSE rerata 2,643 dan 3,278, sedangkan OLS 2,850 dan 3,283. RMSE rerata terkecil dicapai oleh Mixture Lasso (hard) pada Y1 dan Mixture Lasso (hard) pada Y2. Sebagai rujukan, MAE rerata OLS pada laporan ini adalah 2,05 (Y1) dan 2,33 (Y2), yang dapat dibandingkan secara indikatif dengan MAE IRLS pada Tsanas dan Xifara (2012), yaitu 2,14 dan 2,21. Interval empiris Q025 sampai Q975 menggambarkan variasi kinerja antarlipatan dan tidak dibaca sebagai selang kepercayaan.

4.5.4 Perbandingan Berpasangan Lasso vs Mixture Lasso per Lipatan

paired_diff <- cv_results %>%
  filter(Model %in% c("Lasso", "Mixture Lasso (soft)")) %>%
  select(Target, Fold, Model, RMSE) %>%
  pivot_wider(names_from = Model, values_from = RMSE) %>%
  mutate(Diff_RMSE = Lasso - `Mixture Lasso (soft)`)

paired_summary <- paired_diff %>%
  group_by(Target) %>%
  summarise(Mean_Diff_RMSE = mean(Diff_RMSE), SD_Diff_RMSE = sd(Diff_RMSE),
            Q025_Diff_RMSE = quantile(Diff_RMSE, 0.025, names = FALSE),
            Q975_Diff_RMSE = quantile(Diff_RMSE, 0.975, names = FALSE),
            Pct_Folds_Mixture_Better = mean(Diff_RMSE > 0) * 100, .groups = "drop")
tbl(paired_summary, "Tabel 4.19 Selisih RMSE berpasangan per lipatan (Lasso - Mixture soft)")

Tabel 4.19 Selisih RMSE berpasangan per lipatan (Lasso - Mixture soft)

Target	Mean_Diff_RMSE	SD_Diff_RMSE	Q025_Diff_RMSE	Q975_Diff_RMSE	Pct_Folds_Mixture_Better
Y1	0.0419	0.1047	-0.0284	0.1122	50
Y2	-0.0971	0.3193	-0.3116	0.1174	50

print(ggplot(paired_diff, aes(x = Target, y = Diff_RMSE)) +
        geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
        geom_boxplot(fill = "steelblue", alpha = 0.6, outlier.shape = NA) +
        geom_jitter(width = 0.15, alpha = 0.5, size = 1.2) +
        labs(title = "RMSE Lasso - RMSE Mixture Lasso (soft) per fold",
             subtitle = "> 0 : Mixture Lasso lebih akurat pada fold tersebut",
             x = "Target", y = "Selisih RMSE"))

Gambar 4.13 Selisih RMSE berpasangan Lasso dikurangi Mixture Lasso (soft)

Selisih positif berarti Mixture Lasso lebih akurat pada lipatan tersebut. Mixture Lasso soft lebih akurat daripada Lasso pada 50,0% lipatan untuk Y1 (rerata selisih 0,042) dan pada 50,0% lipatan untuk Y2 (rerata selisih -0,097). Perbandingan ini bersifat deskriptif karena lipatan saling bergantung.

4.5.5 Jumlah Komponen \(K\) yang Terpilih

k_freq <- bind_rows(k_records) %>%
  dplyr::count(Target, K, name = "N_Folds") %>%
  group_by(Target) %>% mutate(Percent = N_Folds / sum(N_Folds) * 100) %>% ungroup()
tbl(k_freq, "Tabel 4.20 Frekuensi K terpilih (BIC) pada 50 lipatan luar")

Tabel 4.20 Frekuensi K terpilih (BIC) pada 50 lipatan luar

Target	K	N_Folds	Percent
Y1	3	2	100
Y2	3	2	100

\(K\) yang paling sering terpilih adalah 3 untuk Y1 (100,0% lipatan) dan 3 untuk Y2 (100,0% lipatan). Semakin terkonsentrasi frekuensi pada satu nilai \(K\), semakin stabil struktur rezim yang ditemukan.

4.5.6 Stabilitas Seleksi Variabel

stability_of <- function(records) {
  bind_rows(records) %>%
    group_by(Target, Variable) %>%
    summarise(Selection_Frequency = mean(Selected) * 100,
              Median_Max_Abs_Coefficient = median(Max_Abs_Coefficient), .groups = "drop") %>%
    arrange(Target, desc(Selection_Frequency))
}
selection_stability     <- stability_of(selection_records)
mix_selection_stability <- stability_of(mix_sel_records)

tbl(selection_stability, "Tabel 4.21 Stabilitas seleksi Lasso tunggal pada 50 lipatan luar")

Tabel 4.21 Stabilitas seleksi Lasso tunggal pada 50 lipatan luar

Target	Variable	Selection_Frequency	Median_Max_Abs_Coefficient
Y1	X1	100	6.5032
Y1	X2	100	1.9964
Y1	X4	100	5.2782
Y1	X5	100	7.4926
Y1	X6	100	0.2472
Y1	X7	100	2.2400
Y1	X8	100	4.7179
Y1	X3	0	0.0000
Y2	X1	100	4.4385
Y2	X5	100	8.3216
Y2	X7	100	1.8030
Y2	X8	100	1.4338
Y2	X2	50	1.5221
Y2	X3	50	0.3318
Y2	X4	50	2.5025
Y2	X6	50	0.4720

tbl(mix_selection_stability, "Tabel 4.22 Stabilitas seleksi Mixture Lasso (terpilih pada minimal satu komponen)")

Tabel 4.22 Stabilitas seleksi Mixture Lasso (terpilih pada minimal satu komponen)

Target	Variable	Selection_Frequency	Median_Max_Abs_Coefficient
Y1	X1	100	18.9973
Y1	X2	100	18.6194
Y1	X3	100	1.6258
Y1	X4	100	9.5381
Y1	X5	100	20.0874
Y1	X6	100	0.2334
Y1	X7	100	2.7316
Y1	X8	100	8.1394
Y2	X1	100	13.7595
Y2	X2	100	6.0789
Y2	X3	100	3.3164
Y2	X4	100	9.7444
Y2	X5	100	9.8958
Y2	X6	100	1.1445
Y2	X7	100	2.0832
Y2	X8	100	4.9323

Frekuensi 100% berarti variabel terpilih pada seluruh 50 lipatan luar. Untuk Lasso tunggal pada Y1, variabel yang selalu terpilih adalah X1, X2, X4, X5, X6, X7, X8, sedangkan yang tidak selalu terpilih adalah X3 (0,0%). Pada Y2, variabel yang selalu terpilih adalah X1, X5, X7, X8, sedangkan yang tidak selalu terpilih adalah X2 (50,0%), X3 (50,0%), X4 (50,0%), X6 (50,0%). Untuk Mixture Lasso pada Y1, variabel yang selalu terpilih adalah X1, X2, X3, X4, X5, X6, X7, X8, sedangkan yang tidak selalu terpilih adalah tidak ada. Pada Y2, variabel yang selalu terpilih adalah X1, X2, X3, X4, X5, X6, X7, X8, sedangkan yang tidak selalu terpilih adalah tidak ada. Label komponen antarlipatan tidak dibandingkan langsung, sehingga seleksi campuran dihitung pada level “terpilih di minimal satu komponen”. Variabel dengan frekuensi rendah seleksinya sensitif terhadap sampel latih dan perlu ditafsirkan lebih hati-hati.

4.5.7 Tabel Jawaban Utama

mix_selection_summary <- setNames(purrr::map(TARGETS, function(t)
  summarise_selection(mixture_any_table(fitted_models[[t]]$mixture))), TARGETS)

final_answer <- selection_summary %>%
  select(Variable, Selected_Y1, Max_Abs_Coefficient_Y1, Selected_Y2, Max_Abs_Coefficient_Y2) %>%
  left_join(selection_stability %>% filter(Target == "Y1") %>%
              select(Variable, CV_Selection_Frequency_Y1 = Selection_Frequency), by = "Variable") %>%
  left_join(selection_stability %>% filter(Target == "Y2") %>%
              select(Variable, CV_Selection_Frequency_Y2 = Selection_Frequency), by = "Variable") %>%
  left_join(mix_selection_summary[["Y1"]] %>% select(Variable, Mixture_Selected_Y1 = Selected), by = "Variable") %>%
  left_join(mix_selection_summary[["Y2"]] %>% select(Variable, Mixture_Selected_Y2 = Selected), by = "Variable") %>%
  left_join(mix_selection_stability %>% filter(Target == "Y1") %>%
              select(Variable, Mixture_CV_Freq_Y1 = Selection_Frequency), by = "Variable") %>%
  left_join(mix_selection_stability %>% filter(Target == "Y2") %>%
              select(Variable, Mixture_CV_Freq_Y2 = Selection_Frequency), by = "Variable")

tbl(final_answer, "Tabel 4.23 Jawaban utama: seleksi variabel dan stabilitasnya")

Tabel 4.23 Jawaban utama: seleksi variabel dan stabilitasnya

Variable	Selected_Y1	Max_Abs_Coefficient_Y1	Selected_Y2	Max_Abs_Coefficient_Y2	CV_Selection_Frequency_Y1	CV_Selection_Frequency_Y2	Mixture_Selected_Y1	Mixture_Selected_Y2	Mixture_CV_Freq_Y1	Mixture_CV_Freq_Y2
X1	TRUE	6.8397	TRUE	7.8042	100	100	TRUE	TRUE	100	100
X2	TRUE	2.2763	TRUE	4.1328	100	50	TRUE	TRUE	100	100
X3	FALSE	0.0000	FALSE	0.0000	0	50	TRUE	TRUE	100	100
X4	TRUE	5.4385	TRUE	3.7107	100	50	TRUE	TRUE	100	100
X5	TRUE	7.4229	TRUE	7.8896	100	100	TRUE	TRUE	100	100
X6	TRUE	0.3604	TRUE	0.4565	100	50	TRUE	TRUE	100	100
X7	TRUE	2.2265	TRUE	1.7696	100	100	TRUE	TRUE	100	100
X8	TRUE	4.5301	TRUE	1.8288	100	100	TRUE	TRUE	100	100

Tabel 4.23 menjawab pertanyaan penelitian secara langsung, yaitu variabel yang terpilih pada model final, besar koefisien mutlak maksimumnya, frekuensi seleksinya pada validasi silang berulang, serta status seleksinya pada Mixture Lasso.

4.6 Jawaban atas Pertanyaan Proyek

1. Bimodalitas. Y1 dan Y2 menunjukkan bukti bimodalitas marginal (Tabel 4.4), dan sebagian besar bimodalitas tersebut berimpit dengan dua level X5 (Tabel 4.5). Setelah seluruh prediktor dimasukkan, residual OLS untuk Y1 tidak terbukti multimodal (p-value uji dip 0,5885), sehingga 3 puncak lokal pada densitas kernel residual hanya berupa riak kecil dan untuk Y2 tidak terbukti multimodal (p-value uji dip 0,8894), sehingga 3 puncak lokal pada densitas kernel residual hanya berupa riak kecil. Dengan mempertimbangkan residual, validasi silang, dan interval prediksi, bimodalitas Y1 perlu dimodelkan sebagai campuran, karena RMSE validasi silang turun 1,5% dan interval prediksi campuran lebih tepat dan lebih sempit. Untuk Y2, bimodalitasnya tidak perlu dimodelkan sebagai campuran, karena residual OLS tidak lagi multimodal dan kinerja model campuran tidak lebih baik secara berarti, sehingga Lasso tunggal sudah memadai.

2. Temuan utama seleksi variabel. Lasso tunggal mempertahankan X1, X2, X4, X5, X6, X7, X8 pada model final Y1 dan X1, X2, X4, X5, X6, X7, X8 pada model final Y2. Variabel yang selalu terpilih pada seluruh 50 lipatan luar adalah X1, X2, X4, X5, X6, X7, X8 untuk Y1 dan X1, X5, X7, X8 untuk Y2, sehingga seleksinya stabil. Sebagai contoh penulisan yang tepat: “X7 dipertahankan oleh Lasso pada model final Y1 dan muncul pada 100,0% lipatan luar”. Pernyataan “X7 signifikan pada \(\alpha\) = 5%” dihindari karena koefisien Lasso tidak nol bukan p-value.

3. Mixture Lasso. \(K\) dipilih dengan BIC (Tabel 4.7 dan 4.20), dan label komponen diurutkan menurut rerata respons. Prediksi soft meminimumkan galat kuadrat, sedangkan prediksi hard mempertahankan struktur dua puncak tetapi berisiko salah komponen. Pada data ini, RMSE validasi silang rerata soft tidak lebih kecil daripada hard pada kedua respons, karena peluang keanggotaan komponen umumnya hampir pasti sehingga kedua prediksi hampir sama. Bila kinerja validasi silang Mixture Lasso dan Lasso serupa, model campuran tidak disebut “lebih baik” hanya karena lebih kompleks. Penilaiannya diperkuat dengan Tabel 4.17, yang menunjukkan apakah ketidakpastian prediksi dinyatakan lebih jujur.

4. Inferensi vs prediksi. OLS dipakai untuk membaca residual, memperlihatkan masalah inferensial akibat multikolinearitas sempurna (rank 14 dari 15 kolom), dan menunjukkan sisa bimodalitas. Lasso dan Mixture Lasso lebih relevan untuk seleksi dan prediksi karena tetap terdefinisi pada prediktor yang berkorelasi sempurna dan dinilai melalui kinerja luar sampel.

5. Ketidakpastian. Ketidakpastian dilaporkan melalui rerata dan simpangan baku RMSE, MAE, dan \(R^2\) validasi silang, interval empiris 2,5% hingga 97,5% (Tabel 4.18), frekuensi seleksi setiap variabel (Tabel 4.21 dan 4.22), frekuensi \(K\) terpilih (Tabel 4.20), dan coverage interval prediksi (Tabel 4.17).

6. Keterbatasan. Kesimpulan hanya berlaku pada ruang konfigurasi simulasi (bukan pengukuran bangunan nyata, dengan volume, iklim, material, dan okupansi tetap). Hubungan X-Y dapat tak linear dan mengandung interaksi yang tidak dimodelkan. EM dapat terjebak optimum lokal (karena itu dipakai tiga inisialisasi), hasil sensitif terhadap K_MAX dan MIN_PROP, dan \(\lambda\) awal campuran memakai \(\lambda\) Lasso global. Temuan merupakan hasil model linear terpenalti pada ruang desain simulasi, bukan bukti kausal universal.

Bab V Kesimpulan dan Saran

5.1 Kesimpulan

1. Karakteristik data. Data lengkap (0 nilai hilang, 0 duplikat) dengan prediktor diskret. Prediktor dengan korelasi Spearman terkuat adalah X5 terhadap Y1 dan X5 terhadap Y2. Kedua respons bimodal secara marginal, dan bimodalitas tersebut sebagian besar berimpit dengan tinggi bangunan (X5).
2. Multikolinearitas dan diagnostik. Identitas \(X_2 = X_3 + 2X_4\) menimbulkan multikolinearitas sempurna (rank 14 dari 15 kolom, condition number 2,30e+15, 5 kolom dengan VIF di atas 10). Asumsi OLS dilanggar (heteroskedastisitas dan residual tidak normal pada kedua respons), sehingga inferensi klasik OLS tidak dipakai.
3. Seleksi variabel. Lasso memilih X1, X2, X4, X5, X6, X7, X8 untuk Y1 dan X1, X2, X4, X5, X6, X7, X8 untuk Y2, dengan variabel yang selalu terpilih pada 50 lipatan luar adalah X1, X2, X4, X5, X6, X7, X8 (Y1) dan X1, X5, X7, X8 (Y2). Tereliminasinya variabel dalam kelompok kolinear dibaca sebagai pemilihan perwakilan informasi, bukan ketiadaan pengaruh.
4. Kinerja dan bimodalitas. RMSE validasi silang rerata Lasso adalah 2,848 (Y1) dan 3,293 (Y2), sedangkan Mixture Lasso soft 2,806 dan 3,390. Untuk Y1, bimodalitas perlu dimodelkan sebagai campuran, karena RMSE validasi silang turun 1,5% dan interval prediksi campuran lebih tepat dan lebih sempit. Untuk Y2, bimodalitas tidak perlu dimodelkan sebagai campuran, karena residual OLS tidak lagi multimodal dan kinerja model campuran tidak lebih baik secara berarti, sehingga Lasso tunggal sudah memadai.

5.2 Saran

1. Menambahkan suku interaksi dan polinomial (misalnya Lasso pada basis yang diperluas, atau interaksi dengan X5) untuk menangkap ketaklinearan yang tersisa, lalu membandingkannya dengan model nonlinear seperti random forest atau gradient boosting.
2. Memperkuat stabilitas seleksi dengan stability selection atau bootstrap, terutama untuk variabel yang frekuensi seleksinya di bawah 100%.
3. Membandingkan Lasso dengan elastic net (Zou & Hastie, 2005) yang menangani kelompok prediktor berkorelasi dengan lebih stabil.
4. Memodelkan Y1 dan Y2 secara simultan dengan regresi multirespons karena keduanya berkorelasi sangat tinggi.
5. Menguji sensitivitas Mixture Lasso terhadap K_MAX, MIN_PROP, jumlah inisialisasi, dan nilai \(\lambda\) per komponen.

Daftar Pustaka

Bates, S., Hastie, T., & Tibshirani, R. (2024). Cross-validation: What does it estimate and how well does it do it? Journal of the American Statistical Association, 119(546), 1434–1445. https://arxiv.org/abs/2104.00673

Belsley, D. A., Kuh, E., & Welsch, R. E. (1980). Regression diagnostics: Identifying influential data and sources of collinearity. Wiley. https://doi.org/10.1002/0471725153

Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287–1294. https://doi.org/10.2307/1911963

Cook, R. D. (1977). Detection of influential observation in linear regression. Technometrics, 19(1), 15–18. https://doi.org/10.1080/00401706.1977.10489493

Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B, 39(1), 1–38. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x

Dunn, P. K., & Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5(3), 236–244. https://doi.org/10.1080/10618600.1996.10474708

Friedman, J., Hastie, T., & Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1), 1–22. https://doi.org/10.18637/jss.v033.i01

Gneiting, T., Balabdaoui, F., & Raftery, A. E. (2007). Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society: Series B, 69(2), 243–268. https://doi.org/10.1111/j.1467-9868.2007.00587.x

Hartigan, J. A., & Hartigan, P. M. (1985). The dip test of unimodality. The Annals of Statistics, 13(1), 70–84. https://doi.org/10.1214/aos/1176346577

Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical learning: Data mining, inference, and prediction (2nd ed.). Springer. https://doi.org/10.1007/978-0-387-84858-7

Hoerl, A. E., & Kennard, R. W. (1970). Ridge regression: Biased estimation for nonorthogonal problems. Technometrics, 12(1), 55–67. https://doi.org/10.1080/00401706.1970.10488634

Jacobs, R. A., Jordan, M. I., Nowlan, S. J., & Hinton, G. E. (1991). Adaptive mixtures of local experts. Neural Computation, 3(1), 79–87. https://doi.org/10.1162/neco.1991.3.1.79

Jarque, C. M., & Bera, A. K. (1980). Efficient tests for normality, homoscedasticity and serial independence of regression residuals. Economics Letters, 6(3), 255–259. https://doi.org/10.1016/0165-1765(80)90024-5

Khalili, A., & Chen, J. (2007). Variable selection in finite mixture of regression models. Journal of the American Statistical Association, 102(479), 1025–1038. https://doi.org/10.1198/016214507000000590

Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107–112. https://doi.org/10.1016/0304-4076(81)90062-2

Maechler, M. (2024). diptest: Hartigan’s dip test statistic for unimodality – corrected [R package]. https://CRAN.R-project.org/package=diptest

Pedersen, T. L. (2024). patchwork: The composer of plots [R package]. https://CRAN.R-project.org/package=patchwork

Pfister, R., Schwarz, K. A., Janczyk, M., Dale, R., & Freeman, J. B. (2013). Good things peak in pairs: A note on the bimodality coefficient. Frontiers in Psychology, 4, 700. https://doi.org/10.3389/fpsyg.2013.00700

R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/

Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136

Scrucca, L., Fop, M., Murphy, T. B., & Raftery, A. E. (2016). mclust 5: Clustering, classification and density estimation using Gaussian finite mixture models. The R Journal, 8(1), 289–317. https://doi.org/10.32614/RJ-2016-021

Städler, N., Bühlmann, P., & van de Geer, S. (2010). ℓ1-penalization for mixture regression models. TEST, 19(2), 209–256. https://doi.org/10.1007/s11749-010-0197-z

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B, 58(1), 267–288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x

Tsanas, A., & Xifara, A. (2012). Accurate quantitative estimation of energy performance of residential buildings using statistical machine learning tools. Energy and Buildings, 49, 560–567. https://doi.org/10.1016/j.enbuild.2012.03.003

Tsanas, A., & Xifara, A. (2012). Energy efficiency [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C51307

Wickham, H., & Bryan, J. (2023). readxl: Read Excel files [R package]. https://CRAN.R-project.org/package=readxl

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686

Zeileis, A., & Hothorn, T. (2002). Diagnostic checking in regression relationships. R News, 2(3), 7–10. https://CRAN.R-project.org/doc/Rnews/Rnews_2002-3.pdf

Zou, H., & Hastie, T. (2005). Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B, 67(2), 301–320. https://doi.org/10.1111/j.1467-9868.2005.00503.x

Zou, H., Hastie, T., & Tibshirani, R. (2007). On the “degrees of freedom” of the lasso. The Annals of Statistics, 35(5), 2173–2192. https://doi.org/10.1214/009053607000000127