Prakata

Ebook ini ditulis ulang untuk kebutuhan mahasiswa S2 Statistika. Karena sasaran pembacanya bukan lagi pengantar umum machine learning, penekanan utama buku ini diletakkan pada lima hal. Pertama, setiap bab harus menjelaskan objek statistik yang sedang dimodelkan: apakah targetnya mean bersyarat, probabilitas posterior, struktur laten, fungsi keputusan, atau risk function. Kedua, setiap metode harus dituliskan dalam notasi matematis yang eksplisit agar mahasiswa memahami struktur parameter, asumsi, dan ruang optimisasi yang sedang dihadapi. Ketiga, setiap algoritma harus dijelaskan bukan sekadar sebagai resep komputasi, tetapi sebagai prosedur estimasi parameter atau prosedur pencarian fungsi prediksi yang mengoptimalkan fungsi objektif tertentu. Keempat, implementasi dalam R diposisikan sebagai jembatan antara teori statistika dan komputasi, bukan sebagai pengganti teori. Kelima, perbandingan antar-metode dilakukan dari sudut pandang statistika modern: bias-variance trade-off, regularisasi, generalisasi, identifikasi, stabilitas, dan interpretabilitas.

Dalam tradisi statistika, pertanyaan paling penting tidak berhenti pada “algoritma apa yang paling akurat”, tetapi meluas ke “objek apa yang sedang diestimasi”, “bagaimana loss function didefinisikan”, “apakah estimator teridentifikasi”, “bagaimana perilaku estimator saat dimensi membesar”, dan “bagaimana hasil dipertanggungjawabkan secara ilmiah”. Karena itu, buku ini menempatkan machine learning sebagai kelanjutan dari statistika, bukan sebagai dunia terpisah. Regresi penalized, support vector machine, neural network, clustering, dan gradient boosting semuanya dapat dibaca sebagai variasi masalah optimisasi statistik dengan struktur asumsi, penalti, dan prosedur validasi yang berbeda.

Setiap bab memiliki urutan yang relatif konsisten:

posisi metode dalam statistik dan machine learning,
definisi konseptual,
formulasi matematis,
fungsi objektif atau risk function,
pendekatan estimasi parameter,
isu inferensial dan diagnostik,
implementasi R,
latihan dan tugas.

Pendekatan ini dipilih agar mahasiswa dapat bergerak dari intuisi ke model, lalu dari model ke algoritma, dan akhirnya dari algoritma ke implementasi empiris.

Cara Membaca Buku Ini

0.1 Sasaran

Ebook ini cocok untuk:

mahasiswa S2 Statistika,
dosen yang mengajar machine learning dari perspektif statistik,
peneliti yang ingin memahami koneksi antara inferensi statistik, regularisasi, dan prediksi.

0.2 Notasi Umum

Notasi	Arti
\(n\)	jumlah observasi
\(p\)	jumlah prediktor
\(\mathbf{X}\)	matriks desain berukuran \(n \times p\)
\(\mathbf{y}\)	vektor target
\(f(x)\)	fungsi prediktor
\(\theta\)	parameter umum model
\(\beta\)	parameter regresi
\(L(y, \hat{y})\)	loss function
\(R(f)\)	population risk
\(\widehat{R}(f)\)	empirical risk
\(\lambda\)	parameter regularisasi
\(K\)	jumlah fold atau jumlah cluster, tergantung konteks

0.3 Paket R yang Sering Dipakai

paket <- c(
  "tidyverse", "glmnet", "cluster", "e1071", "class",
  "rpart", "rpart.plot", "randomForest", "gbm", "nnet"
)

install.packages(setdiff(paket, rownames(installed.packages())))

0.4 Peta Buku

Bab	Topik	Fokus Statistik
1	Perkenalan Machine Learning	risk minimization dan generalisasi
2	Regresi dan Gradient Descent	OLS, convex optimization, iterasi gradien
3	Ridge, Lasso, Elastic Net	penalized estimation
4	PCA dan Analisis Faktor	reduksi dimensi dan latent variable model
5	Hierarchical Clustering dan K-Medoids	struktur jarak dan optimisasi cluster
6	K-Means dan Fuzzy C-Means	partitioning objective
7	Gradient Boosting	functional gradient descent
8	UTS Project	workflow analisis
9	Naive Bayes dan K-NN	generative dan instance-based classification
10	Decision Tree	recursive partitioning dan pruning
11	Random Forest	bagging, randomization, OOB risk
12	SVM	margin maximization dan dual optimization
13	SVR	epsilon-insensitive regression
14	Neural Network dan Deep Learning	nonlinear function approximation
15	Validation dan Evaluation	cross-validation, scoring rule, calibration
16	UAS Project	desain studi machine learning

1 Perkenalan Machine Learning

1.1 Tujuan Akademik Bab

Bab ini bertujuan menempatkan machine learning di dalam kerangka statistika modern. Mahasiswa S2 sering kali sudah mengenal regresi, inferensi klasik, dan teori probabilitas, sehingga tantangan pada tahap ini bukan memahami bahwa machine learning “belajar dari data”, tetapi memahami objek statistik apa yang dipelajari, bagaimana risiko didefinisikan, apa bedanya target inferensial dan target prediktif, serta mengapa regularisasi dan validasi menjadi komponen metodologis yang tidak dapat dipisahkan dari pembelajaran mesin.

Secara lebih spesifik, setelah menyelesaikan bab ini mahasiswa diharapkan mampu:

menuliskan masalah machine learning sebagai masalah minimisasi risk,
membedakan population risk, empirical risk, dan regularized empirical risk,
menjelaskan Bayes rule sebagai benchmark teoritis,
membedakan orientasi inferensi, prediksi, dan keputusan statistik,
menjelaskan trade-off bias-variance dan implikasinya pada desain model.

1.2 Machine Learning sebagai Masalah Statistik

Misalkan \((X, Y)\) adalah pasangan peubah acak pada ruang \(\mathcal{X} \times \mathcal{Y}\) dengan distribusi bersama \(P\). Tujuan umum pembelajaran statistik adalah memilih fungsi \(f\) dari suatu kelas fungsi \(\mathcal{F}\) sehingga performa terhadap data baru optimum menurut loss tertentu. Secara formal, jika \(L(Y, f(X))\) adalah loss akibat menggunakan aturan \(f\), maka risk populasi didefinisikan sebagai

\[ R(f) = \mathbb{E}_{(X,Y)\sim P}[L(Y, f(X))]. \]

Karena distribusi \(P\) tidak diketahui, kita bekerja dengan sampel i.i.d.

\[ \mathcal{D}_n = \{(x_i, y_i)\}_{i=1}^{n}. \]

Risk empirisnya adalah

\[ \widehat{R}_n(f) = \frac{1}{n}\sum_{i=1}^{n}L(y_i, f(x_i)). \]

Jika kelas fungsi \(\mathcal{F}\) terlalu kaya, minimisasi langsung atas \(\widehat{R}_n(f)\) dapat menghasilkan overfitting. Karena itu, formulasi yang lebih umum adalah

\[ \widehat{f} = \arg\min_{f \in \mathcal{F}} \left\{ \widehat{R}_n(f) + \lambda \Omega(f) \right\}, \]

dengan \(\Omega(f)\) adalah penalti kompleksitas. Kerangka ini mencakup hampir semua metode dalam buku ini, dari Ridge regression hingga support vector machine dan neural network dengan weight decay.

1.3 Apa yang Sebenarnya Diestimasi?

Dalam tradisi statistika, pertanyaan “apa yang diestimasi?” adalah pertanyaan fundamental. Pada machine learning, pertanyaan ini sering terabaikan karena fokus terlalu cepat bergeser ke algoritma. Padahal, objek statistik yang dipelajari sangat bergantung pada loss function.

1.3.1 Regresi

Jika loss adalah squared error,

\[ L(y, f(x)) = (y-f(x))^2, \]

maka fungsi optimal populasi adalah conditional mean:

\[ f^*(x) = \mathbb{E}[Y \mid X=x]. \]

Jika loss adalah absolute error,

\[ L(y, f(x)) = |y-f(x)|, \]

maka solusi Bayes menjadi conditional median.

1.3.2 Klasifikasi

Jika \(Y \in \{1,\dots,K\}\) dan loss adalah 0-1 loss,

\[ L(y, g(x)) = I(y \neq g(x)), \]

maka aturan optimal adalah Bayes classifier:

\[ g^*(x) = \arg\max_{k} P(Y=k \mid X=x). \]

Jika model mengeluarkan probabilitas, objek yang dipelajari lebih alami adalah posterior probabilities

\[ \eta_k(x) = P(Y=k \mid X=x). \]

1.3.3 Clustering dan Reduksi Dimensi

Pada clustering, tidak ada target \(Y\). Objek yang dipelajari bukan conditional mean atau posterior, tetapi struktur partisi atau representasi laten yang membuat data lebih mudah dijelaskan. Pada PCA, objeknya adalah kombinasi linear yang memaksimalkan variasi. Pada analisis faktor, objeknya adalah latent structure yang menjelaskan kovarians bersama.

1.4 Perspektif Inferensi versus Prediksi

Salah satu sumber kebingungan mahasiswa statistika ketika masuk ke machine learning adalah pencampuran antara tujuan inferensial dan tujuan prediktif. Dalam regresi inferensial klasik, koefisien \(\beta_j\) sering diperlakukan sebagai objek ilmiah utama: apakah signifikan, seberapa besar efeknya, bagaimana interval kepercayaannya. Dalam machine learning prediktif, koefisien atau parameter model dapat menjadi sekunder, sedangkan yang utama adalah performa out-of-sample.

Namun pemisahan ini tidak absolut. Banyak model machine learning tetap dapat diinterpretasikan, dan banyak model statistik klasik dapat dipakai untuk prediksi. Titik perbedaannya bukan semata pada metode, melainkan pada:

target analisis,
cara mengevaluasi model,
jenis keputusan yang hendak diambil.

1.5 Loss Function dan Keputusan Statistik

Loss function adalah jantung machine learning. Ia menerjemahkan tujuan substantif ke dalam ukuran matematis kesalahan. Pilihan loss bukan hal teknis kecil, tetapi keputusan metodologis.

Beberapa contoh loss:

squared loss untuk prediksi mean,
absolute loss untuk prediksi median,
log-loss untuk probabilistic classification,
hinge loss untuk support vector machine,
0-1 loss untuk klasifikasi murni,
custom asymmetric loss untuk situasi di mana underprediction dan overprediction memiliki biaya berbeda.

Dengan demikian, machine learning dapat dibaca sebagai teori estimasi aturan keputusan di bawah loss function tertentu.

1.6 Population Risk, Empirical Risk, dan Generalization Gap

Karena \(R(f)\) tidak dapat dihitung langsung, kita memakai \(\widehat{R}_n(f)\). Tetapi tujuan akhirnya bukan membuat \(\widehat{R}_n(f)\) minimum, melainkan membuat \(R(f)\) kecil. Selisih

\[ R(f) - \widehat{R}_n(f) \]

sering disebut generalization gap. Jika gap besar, model mungkin terlalu mengikuti noise sampel. Di sinilah kompleksitas model menjadi relevan. Kelas fungsi yang sangat fleksibel dapat menghasilkan empirical risk sangat kecil, tetapi dengan generalization gap besar.

Dalam teori statistik pembelajaran, konsep seperti VC dimension, Rademacher complexity, dan uniform convergence dipakai untuk memformalkan hal ini. Dalam buku ini kita tidak masuk terlalu dalam ke teori learning bounds, tetapi intuisi praktisnya harus jelas: semakin fleksibel kelas fungsi, semakin penting regularisasi dan validasi.

1.7 Overfitting, Underfitting, dan Bias-Variance Trade-Off

Jika model terlalu sederhana, ia tidak mampu menangkap struktur data sehingga bias tinggi. Jika model terlalu kompleks, ia sangat sensitif terhadap variasi sampel sehingga variance tinggi. Pada regresi dengan squared loss, decomposition klasik menyatakan

\[ \mathbb{E}\left[(Y-\widehat{f}(X))^2\right] = \sigma^2 + \left(\mathbb{E}[\widehat{f}(X)] - f^*(X)\right)^2 + \mathbb{V}[\widehat{f}(X)]. \]

Interpretasinya:

\(\sigma^2\) adalah irreducible noise,
suku kedua adalah bias squared,
suku ketiga adalah variance estimator.

Regularisasi sering menaikkan bias sedikit tetapi menurunkan variance secara substansial. Itulah sebabnya model regularized sering mengungguli estimator tak bias pada konteks prediksi.

1.8 Bayes Rule sebagai Batas Teoretis

Dalam statistika, Bayes rule adalah benchmark normatif. Untuk 0-1 classification,

\[ g^*(x) = \arg\max_k P(Y=k \mid X=x) \]

adalah aturan yang meminimalkan misclassification probability. Untuk squared loss,

\[ f^*(x) = \mathbb{E}[Y \mid X=x] \]

adalah fungsi terbaik. Maka setiap model praktis dapat dibaca sebagai kompromi antara:

kemampuan menghampiri Bayes rule,
keterbatasan data,
kompleksitas kelas fungsi,
biaya komputasi,
kebutuhan interpretasi.

1.9 No Free Lunch dan Pentingnya Struktur Masalah

Teorema no free lunch menyatakan bahwa tidak ada algoritma tunggal yang unggul untuk semua distribusi data bila dirata-ratakan atas semua kemungkinan masalah. Implikasi praktisnya sangat penting: pemilihan metode harus selalu mempertimbangkan struktur data dan tujuan inferensi/prediksi. Bagi mahasiswa statistika, ini berarti pemahaman data lebih penting daripada mengikuti mode algoritma.

1.10 Tiga Pilar Analisis Machine Learning untuk Statistika

Untuk menjaga kualitas metodologi, buku ini mengusulkan tiga pilar:

1.10.1 Pilar 1: Model

Apa kelas fungsi yang dipakai? Linear, additive, tree-based, kernel, atau compositional?

1.10.2 Pilar 2: Estimasi

Bagaimana parameter ditaksir? Closed-form, gradient descent, coordinate descent, recursive partitioning, EM-like iteration, atau prosedur dual optimization?

1.10.3 Pilar 3: Evaluasi

Bagaimana generalisasi diperiksa? Train-test split, cross-validation, calibration, scoring rule, OOB error, atau time-series backtesting?

Tanpa tiga pilar ini, machine learning mudah berubah menjadi sekadar eksekusi fungsi perangkat lunak.

1.11 Alur Kerja Analisis untuk Mahasiswa S2

Urutan kerja yang sehat adalah:

rumuskan pertanyaan substantif,
identifikasi target statistik,
pilih loss function,
lakukan EDA dan pemeriksaan kualitas data,
tentukan baseline,
pilih model kandidat,
estimasi parameter dan tuning hyperparameter,
evaluasi generalisasi,
interpretasikan hasil dalam konteks ilmiah.

1.12 Contoh R: Train-Test Split dan Baseline Awal

set.seed(123)
idx <- sample(seq_len(nrow(iris)), size = 0.7 * nrow(iris))
train_data <- iris[idx, ]
test_data  <- iris[-idx, ]

dim(train_data)

#> [1] 105   5

dim(test_data)

#> [1] 45  5

prop.table(table(train_data$Species))

#> 
#>     setosa versicolor  virginica 
#>  0.3428571  0.3047619  0.3523810

prop.table(table(test_data$Species))

#> 
#>     setosa versicolor  virginica 
#>  0.3111111  0.4000000  0.2888889

baseline_class <- names(which.max(table(train_data$Species)))
baseline_acc <- mean(test_data$Species == baseline_class)
baseline_class

#> [1] "virginica"

baseline_acc

#> [1] 0.2888889

1.13 Tugas

Tuliskan risk function untuk satu masalah regresi dan satu masalah klasifikasi.
Jelaskan bagaimana pilihan loss memengaruhi objek statistik yang diestimasi.
Berikan contoh situasi di mana model dengan training error sangat kecil justru buruk pada data baru.
Jelaskan mengapa baseline model penting secara metodologis.

1.14 Ringkasan

Machine learning harus dibaca sebagai teori pembelajaran statistik berbasis risk minimization. Fokus utamanya adalah hubungan antara loss function, kelas fungsi, regularisasi, dan performa generalisasi.

2 Analisis Regresi dan Gradient Descent

2.1 Posisi Regresi Linear dalam Machine Learning

Regresi linear adalah fondasi analitik untuk banyak metode machine learning. Di satu sisi, regresi linear adalah model statistik klasik dengan interpretasi parameter yang kuat. Di sisi lain, ia adalah contoh paling bersih dari masalah optimisasi convex. Karena itu, bab ini penting bukan hanya untuk mengulas OLS, tetapi untuk menjembatani analisis parametrik dengan prosedur estimasi numerik.

2.2 Formulasi Model

Misalkan kita memiliki target kontinu \(Y\) dan vektor prediktor \(X=(1,X_1,\ldots,X_p)^\top\). Model linear dinyatakan sebagai

\[ Y_i = X_i^\top \beta + \varepsilon_i, \qquad i=1,\ldots,n, \]

dengan:

\(X_i^\top\) adalah baris ke-\(i\) dari matriks desain,
\(\beta\) adalah vektor parameter,
\(\varepsilon_i\) adalah error acak.

Dalam notasi matriks:

\[ \mathbf{y} = \mathbf{X}\beta + \varepsilon. \]

Jika \(\mathbb{E}[\varepsilon \mid \mathbf{X}] = 0\), maka conditional mean model adalah

\[ \mathbb{E}[\mathbf{y}\mid \mathbf{X}] = \mathbf{X}\beta. \]

Objek statistik yang diestimasi di sini adalah mean bersyarat.

2.3 OLS sebagai M-Estimator

OLS dapat dibaca sebagai M-estimator yang meminimalkan fungsi objektif

\[ Q_n(\beta) = \frac{1}{n}(\mathbf{y}-\mathbf{X}\beta)^\top(\mathbf{y}-\mathbf{X}\beta). \]

Turunan pertama:

\[ \nabla Q_n(\beta) = -\frac{2}{n}\mathbf{X}^\top(\mathbf{y}-\mathbf{X}\beta). \]

Hessian:

\[ \nabla^2 Q_n(\beta) = \frac{2}{n}\mathbf{X}^\top\mathbf{X}. \]

Karena Hessian positif semi-definit, objective bersifat convex. Jika \(\mathbf{X}^\top\mathbf{X}\) full rank, convexity menjadi strict convexity, dan solusi unik.

2.4 Normal Equations dan Solusi Tertutup

Menetapkan gradien sama dengan nol menghasilkan

\[ \mathbf{X}^\top\mathbf{X}\widehat{\beta} = \mathbf{X}^\top\mathbf{y}. \]

Jika matriks invertibel:

\[ \widehat{\beta}_{OLS} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}. \]

Dalam konteks inferensial klasik, di bawah asumsi homoskedastisitas dan eksogenitas, estimator ini tak bias dengan kovarians

\[ \text{Var}(\widehat{\beta}\mid \mathbf{X}) = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}. \]

Dalam konteks prediktif, perhatian beralih pada kestabilan dan error out-of-sample.

2.5 Derivasi Formal Normal Equations

Untuk melihat struktur estimasi secara eksplisit, objective least squares dapat dikembangkan menjadi

\[ Q_n(\beta) = \frac{1}{n} \left( \mathbf{y}^\top\mathbf{y} -2\beta^\top \mathbf{X}^\top \mathbf{y} + \beta^\top \mathbf{X}^\top \mathbf{X}\beta \right). \]

Dengan kalkulus matriks,

\[ \nabla_\beta (\beta^\top A \beta) = (A + A^\top)\beta, \]

dan karena \(\mathbf{X}^\top \mathbf{X}\) simetris, maka

\[ \nabla Q_n(\beta) = -\frac{2}{n}\mathbf{X}^\top \mathbf{y} + \frac{2}{n}\mathbf{X}^\top \mathbf{X}\beta. \]

Menetapkan gradien sama dengan nol menghasilkan

\[ \mathbf{X}^\top \mathbf{X}\widehat{\beta} = \mathbf{X}^\top \mathbf{y}. \]

Pada level magister, derivasi ini penting karena memperlihatkan bahwa normal equations bukan rumus terpisah yang datang dari luar, tetapi konsekuensi langsung dari optimisasi kuadratik convex.

2.6 Asumsi Model dan Implikasinya

Pada level S2, mahasiswa perlu memisahkan asumsi yang diperlukan untuk:

interpretasi mean model,
sifat tak bias,
efisiensi,
inferensi parametrik,
prediksi.

Asumsi \(\mathbb{E}[\varepsilon \mid \mathbf{X}] = 0\) penting untuk unbiasedness. Asumsi homoskedastisitas penting untuk formula varians klasik, tetapi tidak mutlak untuk prediksi. Asumsi normalitas error terutama penting untuk inferensi finite-sample klasik, bukan untuk keberadaan estimator OLS.

2.7 Regresi untuk Inferensi versus Regresi untuk Prediksi

Jika tujuan utama adalah inferensi, kita bertanya:

apakah \(\beta_j\) berbeda dari nol,
seberapa besar efeknya,
bagaimana interval kepercayaan.

Jika tujuan utama adalah prediksi, kita bertanya:

seberapa kecil RMSE pada data baru,
apakah model stabil,
bagaimana trade-off bias-variance.

Perbedaan fokus ini penting karena model yang “baik” untuk inferensi belum tentu optimal untuk prediksi, dan sebaliknya.

2.8 Basis Expansion dan Linearity dalam Parameter

Banyak mahasiswa salah mengira bahwa regresi linear hanya dapat menangkap hubungan garis lurus. Padahal linearity yang dimaksud adalah linear dalam parameter, bukan harus linear dalam prediktor mentah. Misalnya,

\[ Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \varepsilon \]

tetap merupakan model linear dalam parameter \(\beta\). Dengan basis expansion, regresi linear dapat menangkap non-linearitas sederhana sambil tetap mempertahankan struktur optimisasi convex.

2.9 Gradient Descent sebagai Prosedur Estimasi Numerik

Ketika closed-form solution tersedia, mengapa perlu gradient descent? Ada tiga alasan. Pertama, gradient descent adalah fondasi untuk model yang tidak lagi memiliki solusi eksplisit. Kedua, pada data besar, inversi matriks dapat mahal atau tidak stabil. Ketiga, memahami gradient descent berarti memahami cara kerja regularized regression, neural network, dan boosting.

Aturan update gradient descent:

\[ \beta^{(t+1)} = \beta^{(t)} - \eta \nabla Q_n(\beta^{(t)}). \]

Untuk least squares:

\[ \beta^{(t+1)} = \beta^{(t)} + \eta \frac{2}{n}\mathbf{X}^\top(\mathbf{y}-\mathbf{X}\beta^{(t)}). \]

2.10 Konvergensi Gradient Descent

Karena objective least squares convex dan smooth, gradient descent akan konvergen untuk learning rate yang sesuai. Secara teoritis, jika \(\eta\) terlalu besar relatif terhadap eigenvalue terbesar Hessian, iterasi dapat divergen. Secara praktis, scaling prediktor dan pemilihan learning rate konservatif membantu.

2.10.1 Batas Learning Rate dari Spektrum Hessian

Misalkan \(\widehat{\beta}\) adalah solusi OLS dan definisikan error iterasi

\[ e^{(t)} = \beta^{(t)} - \widehat{\beta}. \]

Karena gradien least squares bersifat linear dalam \(\beta\), update gradient descent dapat ditulis sebagai

\[ e^{(t+1)} = \left(I - \eta H\right)e^{(t)}, \qquad H = \frac{2}{n}\mathbf{X}^\top\mathbf{X}. \]

Jika \(\lambda_{\max}(H)\) adalah eigenvalue terbesar Hessian, maka syarat cukup untuk konvergensi adalah

\[ 0 < \eta < \frac{2}{\lambda_{\max}(H)}. \]

Hasil ini menjelaskan dua fakta praktis. Pertama, standardisasi prediktor memperbaiki spektrum Hessian sehingga rentang learning rate yang aman menjadi lebih lebar. Kedua, pada objective kuadratik, perilaku algoritma sepenuhnya ditentukan oleh geometri matriks desain.

2.11 Batch, Stochastic, dan Mini-Batch

2.11.1 Batch Gradient Descent

Memakai seluruh data pada setiap langkah. Stabil, tetapi mahal pada data besar.

2.11.2 Stochastic Gradient Descent

Memakai satu observasi tiap langkah. Lebih noisy, tetapi murah dan bisa lolos dari saddle region pada optimisasi non-convex.

2.11.3 Mini-Batch Gradient Descent

Kompromi antara stabilitas dan efisiensi. Sangat penting pada deep learning.

2.12 Kondisi Numerik dan Standardisasi

Jika kolom \(\mathbf{X}\) memiliki skala sangat berbeda, maka \(\mathbf{X}^\top\mathbf{X}\) dapat poorly conditioned. Secara numerik, ini menyebabkan langkah gradien bergerak sangat lambat di satu arah dan terlalu cepat di arah lain. Standardisasi mengurangi masalah tersebut. Ini salah satu alasan mengapa preprocessing bukan hal kosmetik.

2.13 Regresi Linear sebagai Baseline Wajib

Dalam banyak proyek machine learning, regresi linear harus tetap diperlakukan sebagai baseline yang kuat. Alasannya:

interpretasi jelas,
estimasi sederhana,
menjadi pembanding natural bagi metode lebih kompleks,
sering cukup kompetitif pada data tabular berukuran kecil hingga sedang.

2.14 Simulasi Sederhana untuk Menunjukkan Bias dan Nonlinearitas

Pada praktik S2, penting menunjukkan kapan regresi linear gagal secara sistematis. Salah satu cara adalah membuat data simulasi dari fungsi non-linear dan membandingkan fit linear dengan fit yang lebih fleksibel. Simulasi ini membantu membedakan keterbatasan kelas fungsi dan kualitas estimator.

2.15 Contoh R: Regresi Linear dan Diagnostik Prediktif

fit_lm <- lm(mpg ~ wt + hp + disp, data = mtcars)
summary(fit_lm)

#> 
#> Call:
#> lm(formula = mpg ~ wt + hp + disp, data = mtcars)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -3.891 -1.640 -0.172  1.061  5.861 
#> 
#> Coefficients:
#>              Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept) 37.105505   2.110815  17.579 < 0.0000000000000002 ***
#> wt          -3.800891   1.066191  -3.565              0.00133 ** 
#> hp          -0.031157   0.011436  -2.724              0.01097 *  
#> disp        -0.000937   0.010350  -0.091              0.92851    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2.639 on 28 degrees of freedom
#> Multiple R-squared:  0.8268, Adjusted R-squared:  0.8083 
#> F-statistic: 44.57 on 3 and 28 DF,  p-value: 0.0000000000865

pred_lm <- predict(fit_lm, newdata = mtcars)
res_lm  <- mtcars$mpg - pred_lm

rmse_lm <- sqrt(mean(res_lm^2))
mae_lm  <- mean(abs(res_lm))
r2_lm   <- 1 - sum(res_lm^2) / sum((mtcars$mpg - mean(mtcars$mpg))^2)

c(RMSE = rmse_lm, MAE = mae_lm, R2 = r2_lm)

#>      RMSE       MAE        R2 
#> 2.4684932 1.9070264 0.8268361

2.16 Contoh R: Gradient Descent Manual dan Jejak Konvergensi

dat <- mtcars[, c("mpg", "wt", "hp", "disp")]
y <- dat$mpg
x <- scale(as.matrix(dat[, c("wt", "hp", "disp")]))
X <- cbind(1, x)

beta <- rep(0, ncol(X))
eta <- 0.05
epochs <- 4000
loss_history <- numeric(epochs)

for (t in seq_len(epochs)) {
  error <- X %*% beta - y
  grad  <- (2 / nrow(X)) * t(X) %*% error
  beta  <- beta - eta * grad
  loss_history[t] <- mean(error^2)
}

beta

#>            [,1]
#>      20.0906250
#> wt   -3.7190097
#> hp   -2.1361825
#> disp -0.1161317

tail(loss_history)

#> [1] 6.093459 6.093459 6.093459 6.093459 6.093459 6.093459

plot(loss_history, type = "l", lwd = 2, col = "steelblue",
     xlab = "Iterasi", ylab = "MSE",
     main = "Konvergensi Gradient Descent")

2.17 Contoh R: Closed Form versus Gradient Descent pada Data Hold-Out

set.seed(321)
idx_gd <- sample(seq_len(nrow(mtcars)), size = 0.75 * nrow(mtcars))
train_gd <- mtcars[idx_gd, c("mpg", "wt", "hp", "disp")]
test_gd  <- mtcars[-idx_gd, c("mpg", "wt", "hp", "disp")]

x_train_gd <- scale(train_gd[, c("wt", "hp", "disp")])
y_train_gd <- train_gd$mpg
X_train_gd <- cbind(1, x_train_gd)

beta_cf <- solve(t(X_train_gd) %*% X_train_gd, t(X_train_gd) %*% y_train_gd)

beta_gd <- rep(0, ncol(X_train_gd))
eta_gd <- 0.05

for (iter in 1:5000) {
  err_gd <- X_train_gd %*% beta_gd - y_train_gd
  grad_gd <- (2 / nrow(X_train_gd)) * t(X_train_gd) %*% err_gd
  beta_gd <- beta_gd - eta_gd * grad_gd
}

x_test_gd <- scale(
  test_gd[, c("wt", "hp", "disp")],
  center = attr(x_train_gd, "scaled:center"),
  scale = attr(x_train_gd, "scaled:scale")
)
X_test_gd <- cbind(1, x_test_gd)

pred_cf <- as.numeric(X_test_gd %*% beta_cf)
pred_gd <- as.numeric(X_test_gd %*% beta_gd)

c(
  RMSE_closed_form = sqrt(mean((test_gd$mpg - pred_cf)^2)),
  RMSE_gradient_descent = sqrt(mean((test_gd$mpg - pred_gd)^2)),
  max_abs_coef_diff = max(abs(beta_cf - beta_gd))
)

#>       RMSE_closed_form  RMSE_gradient_descent      max_abs_coef_diff 
#> 1.88486025970556991815 1.88486025970556192455 0.00000000000001421085

2.18 Contoh R: Simulasi Nonlinearitas

set.seed(123)
n <- 150
x <- runif(n, -2, 2)
y <- 1 + 2 * x^2 + rnorm(n, sd = 0.5)

fit_lin <- lm(y ~ x)
fit_quad <- lm(y ~ x + I(x^2))

plot(x, y, pch = 19, col = "gray40")
ord <- order(x)
lines(x[ord], fitted(fit_lin)[ord], col = "red", lwd = 2)
lines(x[ord], fitted(fit_quad)[ord], col = "blue", lwd = 2)
legend("topleft", legend = c("Linear", "Quadratic"),
       col = c("red", "blue"), lwd = 2, bty = "n")

2.19 Tugas

Turunkan normal equations dari objective least squares.
Tunjukkan bahwa Hessian dari least squares adalah positif semi-definit.
Jelaskan mengapa scaling prediktor membantu gradient descent, tetapi tidak mengubah fitted values OLS jika transformasi dikembalikan dengan benar.
Lakukan simulasi data non-linear dan bandingkan model linear dengan model polinomial.

2.20 Ringkasan

Regresi linear adalah fondasi statistik dan komputasional bagi machine learning. Dari bab ini mahasiswa seharusnya memahami perbedaan antara model, estimator, objective function, dan prosedur optimisasi.

3 Ridge, Lasso, dan Elastic Net

3.1 Mengapa Regularisasi Diperlukan

Ketika \(p\) meningkat, ketika prediktor berkorelasi tinggi, atau ketika tujuan utama adalah prediksi out-of-sample, estimator OLS dapat menjadi tidak stabil. Varians estimator membesar dan koefisien mudah berubah drastis akibat fluktuasi sampel. Regularisasi menambahkan penalti kompleksitas untuk menstabilkan estimasi.

Dalam kerangka umum, kita memecahkan

\[ \widehat{\beta} = \arg\min_{\beta_0,\beta} \left\{ \sum_{i=1}^{n}(y_i-\beta_0-x_i^\top\beta)^2 + \lambda P(\beta) \right\}, \]

dengan \(P(\beta)\) adalah penalti.

3.2 Ridge Regression

Pada Ridge,

\[ P(\beta) = \sum_{j=1}^{p}\beta_j^2. \]

Maka objective function menjadi

\[ \widehat{\beta}^{ridge} = \arg\min_{\beta_0, \beta} \left\{ \sum_{i=1}^{n}(y_i-\beta_0-x_i^\top\beta)^2 + \lambda \sum_{j=1}^{p}\beta_j^2 \right\}. \]

Dengan asumsi data telah dipusatkan sehingga intercept terpisah, solusi tertutup adalah

\[ \widehat{\beta}^{ridge} = (\mathbf{X}^\top\mathbf{X} + \lambda I)^{-1}\mathbf{X}^\top \mathbf{y}. \]

Secara aljabar, tambahan \(\lambda I\) memperbaiki conditioning matriks. Secara statistik, ini menukar sedikit bias dengan penurunan variance.

3.2.1 Derivasi Solusi Tertutup Ridge

Dengan mengabaikan intercept yang tidak dipenalti, objective Ridge dapat ditulis sebagai

\[ Q_{ridge}(\beta) = (\mathbf{y}-\mathbf{X}\beta)^\top(\mathbf{y}-\mathbf{X}\beta) + \lambda \beta^\top \beta. \]

Ekspansi aljabar memberi

\[ Q_{ridge}(\beta) = \mathbf{y}^\top\mathbf{y} - 2\beta^\top \mathbf{X}^\top\mathbf{y} + \beta^\top \mathbf{X}^\top\mathbf{X}\beta + \lambda \beta^\top\beta. \]

Turunan pertama terhadap \(\beta\) adalah

\[ \nabla Q_{ridge}(\beta) = -2\mathbf{X}^\top\mathbf{y} +2\mathbf{X}^\top\mathbf{X}\beta +2\lambda \beta. \]

Setelah gradien disamakan dengan nol, diperoleh

\[ (\mathbf{X}^\top\mathbf{X}+\lambda I)\widehat{\beta}^{ridge} = \mathbf{X}^\top\mathbf{y}, \]

yang memberi solusi tertutup

\[ \widehat{\beta}^{ridge} = (\mathbf{X}^\top\mathbf{X}+\lambda I)^{-1}\mathbf{X}^\top\mathbf{y}. \]

Tambahan \(\lambda I\) memastikan matriks yang diinvers lebih stabil bahkan ketika \(\mathbf{X}^\top\mathbf{X}\) mendekati singular.

3.2.2 Representasi Spektral Ridge

Jika \(\mathbf{X}^\top\mathbf{X} = UDU^\top\), maka

\[ \widehat{\beta}^{ridge} = U(D+\lambda I)^{-1}U^\top\mathbf{X}^\top\mathbf{y}. \]

Terlihat bahwa komponen dengan eigenvalue kecil mengalami shrinkage lebih kuat. Ini menjelaskan mengapa Ridge efektif pada multikolinearitas.

3.2.3 Pandangan Bayesian

Jika \(\beta_j \overset{i.i.d.}{\sim} \mathcal{N}(0,\tau^2)\), maka MAP estimator di bawah error Gaussian setara dengan Ridge dengan \(\lambda = \sigma^2/\tau^2\). Interpretasi ini penting karena menunjukkan regularisasi bukan semata trik komputasi, tetapi juga prior belief tentang kecilnya koefisien.

3.3 Lasso

Pada Lasso,

\[ P(\beta) = \sum_{j=1}^{p}|\beta_j|. \]

Sehingga

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

Karena penalti tidak terdiferensialkan di nol, solusi tertutup umum tidak tersedia. Namun sifat sudut pada norma \(L_1\) menyebabkan banyak koefisien menjadi tepat nol.

3.3.1 KKT Conditions

Kondisi KKT untuk Lasso membantu memahami solusi. Untuk setiap \(j\),

\[ 2x_j^\top(\mathbf{y}-\mathbf{X}\widehat{\beta}) = \lambda s_j, \]

dengan

\[ s_j \in \begin{cases} \{+1\}, & \widehat{\beta}_j > 0 \\ [-1,1], & \widehat{\beta}_j = 0 \\ \{-1\}, & \widehat{\beta}_j < 0. \end{cases} \]

Kondisi ini menjelaskan mengapa nol menjadi titik yang “melekat”.

3.3.2 Coordinate Descent dan Soft Thresholding

Untuk menjelaskan bagaimana Lasso diestimasi, pertimbangkan pembaruan satu koefisien pada saat koefisien lain dianggap tetap. Definisikan partial residual

\[ r_j = y - \sum_{\ell \neq j}x_\ell \beta_\ell. \]

Masalah parsial untuk \(\beta_j\) menjadi

\[ \min_{\beta_j} \sum_{i=1}^{n}(r_{ij} - x_{ij}\beta_j)^2 + \lambda |\beta_j|. \]

Jika kolom \(x_j\) telah distandardisasi sehingga \(\sum_i x_{ij}^2 = n\), solusi pembaruan parsial adalah

\[ \widehat{\beta}_j = S\left(\frac{1}{n}\sum_{i=1}^{n}x_{ij}r_{ij}, \frac{\lambda}{2n}\right), \]

dengan operator soft-thresholding

\[ S(z,\gamma) = \text{sign}(z)(|z|-\gamma)_+. \]

Formula ini menjelaskan mengapa glmnet dapat sangat efisien: ia cukup mengulang pembaruan satu koordinat sampai objective stabil.

3.3.3 Soft Thresholding pada Kasus Ortonormal

Jika kolom \(\mathbf{X}\) ortonormal, solusi Lasso menyederhana menjadi

\[ \widehat{\beta}^{lasso}_j = \text{sign}(\widehat{\beta}^{OLS}_j)\left(|\widehat{\beta}^{OLS}_j|-\lambda\right)_+. \]

Ini adalah operator soft-thresholding. Secara konseptual, Lasso tidak hanya mengecilkan, tetapi juga menghapus koefisien kecil.

3.3.4 Pandangan Bayesian

Lasso berkaitan dengan prior Laplace pada koefisien. Berbeda dari prior Gaussian pada Ridge, prior Laplace lebih tajam di nol sehingga mendorong sparsity.

3.4 Elastic Net

Elastic Net memadukan penalti \(L_1\) dan \(L_2\):

\[ \widehat{\beta}^{EN} = \arg\min_{\beta_0,\beta} \left\{ \sum_{i=1}^{n}(y_i-\beta_0-x_i^\top\beta)^2 + \lambda \left[ \alpha \sum_{j=1}^{p}|\beta_j| + \frac{1-\alpha}{2}\sum_{j=1}^{p}\beta_j^2 \right] \right\}. \]

Parameter \(\alpha\) mengontrol derajat sparsity versus shrinkage global. Pada data dengan grup prediktor berkorelasi, Elastic Net sering lebih stabil daripada Lasso murni karena tidak memaksa pemilihan satu variabel secara arbitrer.

3.5 Degrees of Freedom dan Kompleksitas Model

Dalam statistika modern, penalized regression juga dapat dibaca melalui effective degrees of freedom. Pada Ridge, degrees of freedom dapat ditulis sebagai

\[ df(\lambda) = \text{tr}\left\{\mathbf{X}(\mathbf{X}^\top\mathbf{X}+\lambda I)^{-1}\mathbf{X}^\top\right\}. \]

Semakin besar \(\lambda\), effective model complexity menurun. Ide yang sama, meskipun lebih rumit, berlaku juga pada Lasso dan Elastic Net.

3.6 Estimasi Parameter

3.6.1 Ridge

Dapat diselesaikan melalui closed-form, QR/SVD, atau gradient-based optimization.

3.6.2 Lasso dan Elastic Net

Umumnya diestimasi melalui coordinate descent. Pada setiap langkah, satu koefisien diperbarui dengan meminimalkan objective secara parsial, sementara koefisien lain dianggap tetap. Untuk Lasso, pembaruan parsial terkait langsung dengan soft-thresholding.

Coordinate descent sangat efisien karena:

memanfaatkan separabilitas penalti,
cocok untuk jalur solusi berbagai \(\lambda\),
stabil pada dimensi tinggi.

3.7 Pemilihan Parameter Regularisasi

Parameter \(\lambda\) biasanya dipilih melalui \(K\)-fold cross-validation. Dalam glmnet, dua nilai populer adalah:

lambda.min: nilai \(\lambda\) yang meminimalkan error CV,
lambda.1se: nilai \(\lambda\) terbesar yang masih berada dalam satu standard error dari minimum.

Pilihan lambda.1se sering menghasilkan model lebih sederhana dan stabil.

3.8 Hubungan dengan Bias-Variance Trade-Off

Regularisasi memperkenalkan bias:

\[ \mathbb{E}[\widehat{\beta}^{ridge}] \neq \beta \]

dalam pengertian klasik. Namun bias tambahan ini sering dibayar dengan penurunan variance yang cukup besar sehingga error prediksi total menurun. Inilah inti argumen machine learning terhadap estimator tak bias yang tidak stabil.

3.9 Kapan Ridge, Lasso, atau Elastic Net Dipakai

3.9.1 Ridge

Dipakai ketika:

sebagian besar variabel dipercaya relevan,
multikolinearitas tinggi,
tujuan utama adalah prediksi stabil.

3.9.2 Lasso

Dipakai ketika:

diyakini hanya sebagian kecil prediktor penting,
interpretasi model ringkas diinginkan,
seleksi variabel eksplisit dibutuhkan.

3.9.3 Elastic Net

Dipakai ketika:

terdapat grup prediktor berkorelasi,
tetap diinginkan sparsity,
Lasso tampak tidak stabil.

3.10 Contoh R: Jalur Regularisasi dan Perbandingan Koefisien

library(glmnet)

x <- model.matrix(mpg ~ ., data = mtcars)[, -1]
y <- mtcars$mpg

cv_ridge <- cv.glmnet(x, y, alpha = 0)
cv_lasso <- cv.glmnet(x, y, alpha = 1)
cv_enet  <- cv.glmnet(x, y, alpha = 0.5)

coef(cv_ridge, s = "lambda.min")

#> 11 x 1 sparse Matrix of class "dgCMatrix"
#>               lambda.min
#> (Intercept) 21.214906654
#> cyl         -0.364744539
#> disp        -0.005091067
#> hp          -0.011794167
#> drat         1.049747184
#> wt          -1.294217685
#> qsec         0.167934827
#> vs           0.740157846
#> am           1.687219524
#> gear         0.549966315
#> carb        -0.572119610

coef(cv_lasso, s = "lambda.min")

#> 11 x 1 sparse Matrix of class "dgCMatrix"
#>             lambda.min
#> (Intercept) 35.7308062
#> cyl         -0.8801087
#> disp         .        
#> hp          -0.0110762
#> drat         .        
#> wt          -2.6636936
#> qsec         .        
#> vs           .        
#> am           .        
#> gear         .        
#> carb         .

coef(cv_enet, s = "lambda.min")

#> 11 x 1 sparse Matrix of class "dgCMatrix"
#>               lambda.min
#> (Intercept) 31.573809843
#> cyl         -0.661421522
#> disp        -0.002187216
#> hp          -0.013652413
#> drat         0.539328435
#> wt          -2.031340986
#> qsec         .          
#> vs           0.186610837
#> am           0.972770448
#> gear         .          
#> carb        -0.271853921

plot(cv_ridge)

plot(cv_lasso)

plot(cv_enet)

3.11 Contoh R: Membandingkan RMSE In-Sample secara Sederhana

pred_ridge <- as.numeric(predict(cv_ridge, newx = x, s = "lambda.min"))
pred_lasso <- as.numeric(predict(cv_lasso, newx = x, s = "lambda.min"))
pred_enet  <- as.numeric(predict(cv_enet,  newx = x, s = "lambda.min"))

c(
  RMSE_Ridge = sqrt(mean((y - pred_ridge)^2)),
  RMSE_Lasso = sqrt(mean((y - pred_lasso)^2)),
  RMSE_ENet  = sqrt(mean((y - pred_enet)^2))
)

#> RMSE_Ridge RMSE_Lasso  RMSE_ENet 
#>   2.298226   2.540335   2.379778

3.12 Contoh R: Hold-Out Comparison dan Kompleksitas Model

set.seed(456)
idx_pen <- sample(seq_len(nrow(mtcars)), size = 0.75 * nrow(mtcars))
train_pen <- mtcars[idx_pen, ]
test_pen  <- mtcars[-idx_pen, ]

x_train_pen <- model.matrix(mpg ~ ., data = train_pen)[, -1]
y_train_pen <- train_pen$mpg
x_test_pen  <- model.matrix(mpg ~ ., data = test_pen)[, -1]
y_test_pen  <- test_pen$mpg

cv_ridge_hold <- cv.glmnet(x_train_pen, y_train_pen, alpha = 0)
cv_lasso_hold <- cv.glmnet(x_train_pen, y_train_pen, alpha = 1)
cv_enet_hold  <- cv.glmnet(x_train_pen, y_train_pen, alpha = 0.5)

pred_ridge_hold <- as.numeric(predict(cv_ridge_hold, newx = x_test_pen, s = "lambda.1se"))
pred_lasso_hold <- as.numeric(predict(cv_lasso_hold, newx = x_test_pen, s = "lambda.1se"))
pred_enet_hold  <- as.numeric(predict(cv_enet_hold,  newx = x_test_pen, s = "lambda.1se"))

data.frame(
  model = c("Ridge", "Lasso", "Elastic Net"),
  lambda = c(cv_ridge_hold$lambda.1se, cv_lasso_hold$lambda.1se, cv_enet_hold$lambda.1se),
  rmse = c(
    sqrt(mean((y_test_pen - pred_ridge_hold)^2)),
    sqrt(mean((y_test_pen - pred_lasso_hold)^2)),
    sqrt(mean((y_test_pen - pred_enet_hold)^2))
  ),
  nonzero = c(
    sum(abs(as.matrix(coef(cv_ridge_hold, s = "lambda.1se"))[-1, 1]) > 0),
    sum(abs(as.matrix(coef(cv_lasso_hold, s = "lambda.1se"))[-1, 1]) > 0),
    sum(abs(as.matrix(coef(cv_enet_hold, s = "lambda.1se"))[-1, 1]) > 0)
  )
)

#>         model    lambda     rmse nonzero
#> 1       Ridge 16.054248 4.020630      10
#> 2       Lasso  1.845848 4.129358       3
#> 3 Elastic Net  3.064911 4.206565       5

3.13 Tugas

Turunkan solusi Ridge dari objective function-nya.
Jelaskan secara geometris mengapa Lasso menghasilkan sparsity.
Tunjukkan bagaimana KKT conditions menjelaskan kemungkinan koefisien nol pada Lasso.
Bandingkan lambda.min dan lambda.1se pada satu dataset pilihan Anda.

3.14 Ringkasan

Ridge, Lasso, dan Elastic Net adalah contoh utama penalized estimation dalam statistika modern. Ketiganya memperlihatkan bagaimana bias yang dikendalikan dapat menghasilkan penurunan variance dan perbaikan prediksi.

4 Reduksi Dimensi: PCA dan Analisis Faktor

4.1 Mengapa Reduksi Dimensi Penting

Dalam data multivariat, jumlah variabel sering besar dan korelasi antarprediktor tinggi. Hal ini menimbulkan beberapa persoalan:

visualisasi menjadi sulit,
multikolinearitas mengganggu estimasi,
dimensi tinggi meningkatkan risiko overfitting,
interpretasi struktur data menjadi rumit.

Reduksi dimensi bertujuan mengganti himpunan variabel asal dengan sejumlah kecil representasi yang tetap menangkap informasi utama.

4.2 Principal Component Analysis

PCA berangkat dari pertanyaan: kombinasi linear mana dari variabel asal yang menjelaskan variasi terbesar? Misalkan \(X\) telah dipusatkan. Komponen utama pertama adalah

\[ z_1 = a_1^\top X \]

dengan \(a_1\) dipilih untuk memaksimalkan

\[ \text{Var}(a^\top X) = a^\top \Sigma a \]

di bawah kendala \(\|a\|^2 = 1\). Maka masalah optimisasi menjadi

\[ \max_{a:\|a\|=1} a^\top \Sigma a. \]

Dengan metode Lagrange diperoleh

\[ \Sigma a = \lambda a, \]

sehingga \(a_1\) adalah eigenvector dari \(\Sigma\) yang berasosiasi dengan eigenvalue terbesar. Komponen kedua dicari dengan kendala ortogonal terhadap komponen pertama, dan seterusnya.

4.2.1 Derivasi Rayleigh Quotient

Masalah Lagrange untuk komponen pertama adalah

\[ \mathcal{L}(a,\lambda) = a^\top \Sigma a - \lambda(a^\top a - 1). \]

Turunan terhadap \(a\) memberi

\[ 2\Sigma a - 2\lambda a = 0 \quad \Longrightarrow \quad \Sigma a = \lambda a. \]

Nilai objektif pada solusi adalah

\[ a^\top \Sigma a = \lambda, \]

sehingga memaksimalkan varians proyeksi ekuivalen dengan memilih eigenvalue terbesar. Hasil ini penting karena memberi justifikasi formal mengapa PCA identik dengan eigen decomposition.

4.3 PCA sebagai Masalah Aproksimasi

PCA tidak hanya dapat dibaca sebagai maksimisasi varians, tetapi juga sebagai minimisasi reconstruction error. Jika data diproyeksikan ke subruang berdimensi \(k\), PCA menghasilkan subruang yang meminimalkan total squared projection error. Formulasi ini penting karena menunjukkan hubungan PCA dengan reduced-rank approximation dan singular value decomposition.

4.4 Singular Value Decomposition

Jika matriks data terpusat adalah \(\mathbf{X}\), maka

\[ \mathbf{X} = UDV^\top. \]

Kolom-kolom \(V\) adalah principal directions, sedangkan \(\mathbf{X}V\) memberi scores. Dalam praktik komputasi, SVD sering lebih stabil daripada langsung mendekomposisi matriks kovarians.

4.5 Standarisasi dan Pilihan Matriks

Jika variabel memiliki satuan berbeda, PCA pada matriks kovarians akan didominasi variabel dengan varians besar. Karena itu, pada banyak aplikasi sosial dan ekonomi, PCA dilakukan pada matriks korelasi atau data yang telah distandardisasi.

4.6 Interpretasi Loading, Score, dan Proporsi Varian

Loading adalah koefisien kombinasi linear. Score adalah posisi observasi dalam sistem koordinat baru. Eigenvalue menyatakan varians yang dijelaskan komponen. Proporsi kumulatif varian membantu memutuskan berapa komponen yang dipertahankan.

4.7 Analisis Faktor

Berbeda dari PCA, analisis faktor berangkat dari model probabilistik:

\[ X = \mu + \Lambda F + \epsilon, \]

dengan:

\(F\) faktor laten,
\(\Lambda\) matriks loading,
\(\epsilon\) error spesifik.

Asumsi umum:

\[ \mathbb{E}[F]=0,\qquad \text{Cov}(F)=\Phi, \]

\[ \mathbb{E}[\epsilon]=0,\qquad \text{Cov}(\epsilon)=\Psi, \]

dengan \(\Psi\) diagonal.

Maka kovarians total adalah

\[ \Sigma = \Lambda \Phi \Lambda^\top + \Psi. \]

Jika \(\Phi=I\), kita memiliki orthogonal factor model.

4.8 Common Variance dan Unique Variance

Konsep kunci pada analisis faktor adalah pemisahan antara common variance dan unique variance. Communality variabel ke-\(j\) adalah proporsi variansi yang dijelaskan faktor bersama, sedangkan uniqueness adalah bagian variansi yang spesifik pada variabel tersebut.

4.9 Identifikasi Model Faktor

Mahasiswa S2 perlu memperhatikan bahwa model faktor tidak sepenuhnya teridentifikasi tanpa pembatasan tambahan, karena rotasi tertentu pada faktor dapat menghasilkan kovarians yang sama. Itulah sebabnya interpretasi faktor selalu terkait dengan isu rotational indeterminacy.

4.10 Estimasi Parameter pada Analisis Faktor

4.10.1 Principal Factor Method

Metode ini menggunakan pendekatan berbasis reduksi kovarians bersama.

4.10.2 Maximum Likelihood

Jika data diasumsikan normal multivariat, likelihood dapat ditulis dan parameter \(\Lambda\) serta \(\Psi\) diestimasi melalui optimisasi numerik. Uji kecocokan model faktor sering dibangun dari pendekatan likelihood ini.

4.10.3 Likelihood dan Identifikasi Praktis

Jika \(X_i \sim \mathcal{N}_p(\mu, \Sigma)\) dengan

\[ \Sigma = \Lambda \Lambda^\top + \Psi, \]

maka log-likelihood terpusat, hingga konstanta aditif, adalah

\[ \ell(\Lambda,\Psi) = -\frac{n}{2} \left\{ \log |\Sigma| + \text{tr}(S\Sigma^{-1}) \right\}, \]

dengan \(S\) kovarians sampel. Estimasi parameter pada factanal() dilakukan dengan optimisasi numerik terhadap fungsi ini. Secara konseptual, kompleksitas utama bukan hanya mencari maksimum likelihood, tetapi memastikan solusi tetap teridentifikasi dan substantif setelah rotasi.

4.10.4 Factor Scores

Setelah loading diperoleh, skor faktor dapat dihitung dengan metode regression score atau Bartlett score. Perlu diingat bahwa faktor laten tidak teramati, sehingga score yang diperoleh adalah estimator, bukan observasi langsung.

4.11 Rotasi Faktor

Rotasi orthogonal seperti varimax bertujuan membuat struktur loading lebih “sederhana” sehingga interpretasi faktor menjadi substantif. Rotasi oblique mengizinkan faktor saling berkorelasi, yang sering lebih realistis pada aplikasi sosial.

4.12 PCA versus Analisis Faktor

Perbedaan mendasar:

PCA adalah dekomposisi variasi total,
analisis faktor adalah model laten untuk kovarians bersama.

Karena itu, PCA cocok untuk kompresi data dan feature extraction, sedangkan analisis faktor lebih cocok untuk pengukuran konstruk laten.

4.13 Pemilihan Jumlah Komponen atau Faktor

Untuk PCA:

scree plot,
Kaiser rule,
proporsi varian kumulatif,
parallel analysis.

Untuk analisis faktor:

interpretabilitas,
likelihood-based fit,
residual correlation structure,
parallel analysis atau kriteria lain.

4.14 Contoh R: PCA dengan Interpretasi Loading

pca_fit <- prcomp(USArrests, scale. = TRUE)
summary(pca_fit)

#> Importance of components:
#>                           PC1    PC2     PC3     PC4
#> Standard deviation     1.5749 0.9949 0.59713 0.41645
#> Proportion of Variance 0.6201 0.2474 0.08914 0.04336
#> Cumulative Proportion  0.6201 0.8675 0.95664 1.00000

pca_fit$rotation

#>                 PC1        PC2        PC3         PC4
#> Murder   -0.5358995 -0.4181809  0.3412327  0.64922780
#> Assault  -0.5831836 -0.1879856  0.2681484 -0.74340748
#> UrbanPop -0.2781909  0.8728062  0.3780158  0.13387773
#> Rape     -0.5434321  0.1673186 -0.8177779  0.08902432

pca_fit$x[1:10, ]

#>                     PC1         PC2         PC3          PC4
#> Alabama     -0.97566045 -1.12200121  0.43980366  0.154696581
#> Alaska      -1.93053788 -1.06242692 -2.01950027 -0.434175454
#> Arizona     -1.74544285  0.73845954 -0.05423025 -0.826264240
#> Arkansas     0.13999894 -1.10854226 -0.11342217 -0.180973554
#> California  -2.49861285  1.52742672 -0.59254100 -0.338559240
#> Colorado    -1.49934074  0.97762966 -1.08400162  0.001450164
#> Connecticut  1.34499236  1.07798362  0.63679250 -0.117278736
#> Delaware    -0.04722981  0.32208890  0.71141032 -0.873113315
#> Florida     -2.98275967 -0.03883425  0.57103206 -0.095317042
#> Georgia     -1.62280742 -1.26608838  0.33901818  1.065974459

biplot(pca_fit, cex = 0.7)

4.15 Contoh R: Analisis Faktor Satu Faktor

fa_fit <- factanal(
  USArrests,
  factors = 1,
  rotation = "varimax",
  scores = "regression"
)

print(fa_fit, digits = 2, cutoff = 0.30)

#> 
#> Call:
#> factanal(x = USArrests, factors = 1, scores = "regression", rotation = "varimax")
#> 
#> Uniquenesses:
#>   Murder  Assault UrbanPop     Rape 
#>     0.33     0.04     0.93     0.53 
#> 
#> Loadings:
#>          Factor1
#> Murder   0.82   
#> Assault  0.98   
#> UrbanPop        
#> Rape     0.68   
#> 
#>                Factor1
#> SS loadings       2.16
#> Proportion Var    0.54
#> 
#> Test of the hypothesis that 1 factor is sufficient.
#> The chi square statistic is 9.15 on 2 degrees of freedom.
#> The p-value is 0.0103

fa_fit$loadings

#> 
#> Loadings:
#>          Factor1
#> Murder   0.818  
#> Assault  0.979  
#> UrbanPop 0.262  
#> Rape     0.683  
#> 
#>                Factor1
#> SS loadings      2.162
#> Proportion Var   0.540

fa_fit$scores[1:10, ]

#>     Alabama      Alaska     Arizona    Arkansas  California    Colorado 
#>   0.7901515   1.1161087   1.3553508   0.2025243   1.2423296   0.4472463 
#> Connecticut    Delaware     Florida     Georgia 
#>  -0.7724846   0.6409789   1.9416162   0.6412841

4.16 Contoh R: Rekonstruksi Aproksimasi PCA

x_scaled <- scale(USArrests)
pca2 <- prcomp(x_scaled)
scores2 <- pca2$x[, 1:2]
loadings2 <- pca2$rotation[, 1:2]
approx2 <- scores2 %*% t(loadings2)
head(approx2)

#>               Murder   Assault   UrbanPop       Rape
#> Alabama    0.9920554 0.7799093 -0.7078698  0.3424735
#> Alaska     1.4788608 1.3255791 -0.3902348  0.8713524
#> Arizona    0.6265723 0.8790939  1.1300983  1.0720877
#> Arkansas   0.3885458 0.1267449 -1.0064890 -0.2615597
#> California 0.7002647 1.1700159  2.0282388  1.6133934
#> Colorado   0.3946699 0.6906107  1.2703841  0.9783655

4.17 Contoh R: Proporsi Varian dan Communality

prop_var <- pca_fit$sdev^2 / sum(pca_fit$sdev^2)
round(prop_var, 4)

#> [1] 0.6201 0.2474 0.0891 0.0434

communalities <- 1 - fa_fit$uniquenesses
round(communalities, 4)

#>   Murder  Assault UrbanPop     Rape 
#>   0.6685   0.9585   0.0686   0.4664

data.frame(
  variable = names(communalities),
  communality = round(communalities, 4),
  uniqueness = round(fa_fit$uniquenesses, 4)
)

#>          variable communality uniqueness
#> Murder     Murder      0.6685     0.3315
#> Assault   Assault      0.9585     0.0415
#> UrbanPop UrbanPop      0.0686     0.9314
#> Rape         Rape      0.4664     0.5336

4.18 Tugas

Turunkan masalah optimisasi PCA dengan kendala Lagrange.
Jelaskan mengapa PCA dapat dipandang sebagai masalah minimisasi reconstruction error.
Tuliskan perbedaan fundamental PCA dan analisis faktor dari sudut objek statistik.
Terapkan PCA dan analisis faktor pada satu dataset multivariat dan bandingkan interpretasinya.

4.19 Ringkasan

PCA dan analisis faktor sama-sama mereduksi dimensi, tetapi bertumpu pada landasan statistik yang berbeda. PCA bersifat geometrik dan berbasis variasi total, sedangkan analisis faktor bersifat model-based dan berfokus pada struktur laten kovarians.

5 Unsupervised Learning I: Hierarchical Clustering dan K-Medoids

5.1 Tujuan Akademik Bab

Bab ini memperkenalkan clustering dari sudut pandang statistika eksploratif dan optimisasi. Karena tidak ada target teramati, kita tidak berbicara tentang unbiasedness atau prediction error dalam arti supervised learning. Sebaliknya, kita berbicara tentang struktur dissimilarity, definisi partisi, sifat algoritma pengelompokan, dan validitas interpretasi cluster. Untuk mahasiswa S2, poin pentingnya adalah memahami bahwa hasil clustering sangat tergantung pada:

definisi jarak,
skala variabel,
bentuk cluster yang diasumsikan secara implisit,
algoritma optimisasi yang dipakai.

5.2 Clustering sebagai Masalah Statistik

Misalkan data kita adalah \(\{x_1,\ldots,x_n\}\), dengan \(x_i \in \mathbb{R}^p\). Tujuan clustering bukan memprediksi \(Y\), melainkan menemukan partisi

\[ \mathcal{C} = \{C_1,\ldots,C_K\} \]

yang menangkap struktur kedekatan dalam data. Karena struktur ini tidak teramati langsung, clustering lebih dekat ke unsupervised model selection daripada ke inferensi parametris klasik.

Secara umum, banyak algoritma clustering dapat ditulis sebagai minimisasi objective:

\[ \min_{\mathcal{C}} \Phi(\mathcal{C}; \mathbf{D}), \]

dengan \(\mathbf{D}\) matriks jarak atau dissimilarity. Maka pemilihan jarak adalah bagian inti dari model, bukan keputusan teknis kecil.

5.3 Matriks Jarak dan Dissimilarity

Untuk data numerik, jarak Euclidean didefinisikan sebagai

\[ d_E(x_i,x_j) = \left(\sum_{m=1}^{p}(x_{im}-x_{jm})^2\right)^{1/2}. \]

Jarak Manhattan:

\[ d_M(x_i,x_j) = \sum_{m=1}^{p}|x_{im}-x_{jm}|. \]

Jarak Mahalanobis:

\[ d_{Mah}(x_i,x_j) = \left\{ (x_i-x_j)^\top \Sigma^{-1}(x_i-x_j) \right\}^{1/2}. \]

Mahalanobis penting secara statistik karena memperhitungkan korelasi antarprediktor. Untuk data campuran, jarak Gower memberi solusi lebih realistis.

5.4 Standardisasi dan Invariansi Skala

Tanpa standardisasi, variabel dengan skala besar akan mendominasi jarak. Karena itu, untuk banyak aplikasi multivariat, standar awal yang sehat adalah:

\[ z_{ij} = \frac{x_{ij}-\bar{x}_j}{s_j}. \]

Namun, mahasiswa S2 harus menyadari bahwa standardisasi juga mengandung asumsi substantif: semua variabel diberi bobot yang setara dalam penentuan cluster.

5.5 Hierarchical Clustering

Hierarchical clustering membentuk pohon pengelompokan. Pada agglomerative clustering, proses dimulai dari \(n\) singleton cluster lalu menggabungkan cluster secara bertahap sampai semua observasi berada dalam satu cluster.

Output utamanya bukan satu partisi, melainkan struktur bertingkat yang dapat dipotong pada level tertentu. Ini membuat hierarchical clustering lebih eksploratif daripada metode partisi tunggal.

5.6 Linkage Criterion

Jika \(A\) dan \(B\) dua cluster, beberapa definisi jarak antarkluster adalah:

5.6.1 Single Linkage

\[ d(A,B) = \min_{i \in A, j \in B} d(x_i, x_j). \]

Ini cenderung menghasilkan chaining effect, yaitu cluster panjang memanjang karena cukup ada satu pasangan titik dekat.

5.6.2 Complete Linkage

\[ d(A,B) = \max_{i \in A, j \in B} d(x_i, x_j). \]

Metode ini cenderung menghasilkan cluster lebih kompak.

5.6.3 Average Linkage

\[ d(A,B) = \frac{1}{|A||B|}\sum_{i \in A}\sum_{j \in B} d(x_i, x_j). \]

Ini memberi kompromi antara single dan complete linkage.

5.6.4 Ward Linkage

Ward linkage meminimalkan kenaikan within-cluster sum of squares:

\[ \Delta(A,B) = \frac{|A||B|}{|A|+|B|} \|\bar{x}_A - \bar{x}_B\|^2. \]

Karena itu, Ward linkage sering dipandang paling dekat ke prinsip variance minimization.

5.7 Lance-Williams Recurrence

Banyak metode hierarchical clustering dapat ditulis dalam bentuk Lance-Williams recurrence, yang memperbarui jarak cluster baru \((A \cup B)\) terhadap cluster \(C\) sebagai kombinasi dari jarak lama. Hal ini penting secara komputasional dan menunjukkan bahwa linkage method sebenarnya adalah aturan pembaruan jarak yang berbeda-beda.

5.8 Dendrogram dan Interpretasi

Dendrogram merekam sejarah penggabungan cluster. Tinggi cabang dapat ditafsirkan sebagai tingkat dissimilarity saat dua cluster bergabung. Dalam praktik, dendrogram membantu menjawab:

apakah ada kelompok yang sangat terpisah,
apakah struktur cluster stabil,
pada level berapa partisi menjadi substantif.

Namun dendrogram bukan bukti bahwa jumlah cluster tertentu adalah “benar”; ia hanya alat eksplorasi.

5.9 K-Medoids

K-Medoids mencari partisi dan medoid

\[ \{m_1,\ldots,m_K\} \]

dengan objective

\[ \min_{C_1,\ldots,C_K} \sum_{k=1}^{K}\sum_{i \in C_k} d(x_i,m_k). \]

Karena medoid harus merupakan observasi aktual, K-Medoids lebih robust terhadap outlier dibanding K-Means. Dari sudut pandang robust statistics, ini cukup natural karena pusat cluster tidak ditentukan oleh mean.

5.10 PAM: Partitioning Around Medoids

PAM bekerja dalam dua fase besar:

build phase: memilih medoid awal,
swap phase: mencoba pertukaran medoid dengan non-medoid dan menerima pertukaran yang menurunkan objective.

PAM adalah algoritma optimisasi diskrit. Ia tidak memiliki jaminan global optimum, tetapi relatif stabil untuk data berukuran sedang.

5.11 Evaluasi Clustering

Karena tidak ada target, evaluasi clustering biasanya bersifat internal atau eksternal.

5.11.1 Silhouette Width

Untuk observasi \(i\),

\[ s(i) = \frac{b(i)-a(i)}{\max\{a(i),b(i)\}}, \]

dengan:

\(a(i)\): rata-rata jarak observasi \(i\) ke anggota cluster sendiri,
\(b(i)\): jarak rata-rata minimum ke cluster lain.

Nilai \(s(i)\) mendekati 1 menunjukkan alokasi cluster yang baik.

5.11.2 Validitas Substantif

Nilai numerik tinggi tidak selalu berarti cluster ilmiah bermakna. Mahasiswa S2 harus selalu bertanya apakah cluster dapat dijelaskan dalam konteks domain.

5.12 Contoh R: Hierarchical Clustering

scaled_data <- scale(USArrests)
d <- dist(scaled_data, method = "euclidean")
hc_ward <- hclust(d, method = "ward.D2")
plot(hc_ward, main = "Hierarchical Clustering dengan Ward")

cluster_hc <- cutree(hc_ward, k = 3)
table(cluster_hc)

#> cluster_hc
#>  1  2  3 
#> 19 19 12

5.13 Contoh R: K-Medoids

library(cluster)
pam_fit <- pam(scaled_data, k = 3)
pam_fit$medoids

#>                   Murder    Assault   UrbanPop       Rape
#> New Mexico     0.8292944  1.3708088  0.3081225  1.1603196
#> Oklahoma      -0.2727580 -0.2371077  0.1699510 -0.1315342
#> New Hampshire -1.3059321 -1.3650491 -0.6590781 -1.2525644

table(pam_fit$clustering)

#> 
#>  1  2  3 
#> 19 21 10

plot(pam_fit)

5.14 Contoh R: Silhouette Analysis

sil <- silhouette(pam_fit$clustering, dist(scaled_data))
summary(sil)

#> Silhouette of 50 units in 3 clusters from silhouette.default(x = pam_fit$clustering, dist = dist(scaled_data)) :
#>  Cluster sizes and average silhouette widths:
#>        19        21        10 
#> 0.2757843 0.2797126 0.4604416 
#> Individual silhouette widths:
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -0.1723  0.2192  0.3339  0.3144  0.4227  0.5910

plot(sil)

5.15 Tugas

Bandingkan complete linkage dan Ward linkage pada satu dataset multivariat.
Jelaskan bagaimana pilihan jarak dapat mengubah interpretasi cluster.
Terapkan K-Medoids dan bandingkan silhouette-nya dengan hierarchical clustering.
Tuliskan kritik substantif atas hasil cluster yang Anda peroleh.

5.16 Ringkasan

Hierarchical clustering dan K-Medoids menunjukkan dua strategi utama dalam unsupervised learning berbasis jarak: struktur bertingkat dan partisi robust. Inti statistiknya terletak pada definisi jarak, objective function, dan interpretasi hasil.

6 Unsupervised Learning II: K-Means dan Fuzzy C-Means

6.1 Tujuan Akademik Bab

Bab ini membahas partitional clustering dari dua sudut: assignment keras melalui K-Means dan assignment lunak melalui Fuzzy C-Means. Mahasiswa S2 perlu memahami bahwa kedua metode ini dapat dibaca sebagai prosedur optimisasi dengan asumsi bentuk cluster tertentu.

6.2 K-Means sebagai Masalah Optimisasi

K-Means mencari partisi \(\{C_1,\ldots,C_K\}\) dan centroid \(\{\mu_1,\ldots,\mu_K\}\) yang meminimalkan within-cluster sum of squares:

\[ \min_{C_1,\ldots,C_K,\mu_1,\ldots,\mu_K} \sum_{k=1}^{K}\sum_{i \in C_k}\|x_i-\mu_k\|^2. \]

Objective ini tidak convex secara simultan terhadap assignment dan centroid, tetapi conditional convex terhadap masing-masing blok. Itulah dasar Lloyd’s algorithm.

6.3 Lloyd’s Algorithm

Dengan centroid tetap, assignment optimal adalah memberi setiap observasi ke centroid terdekat. Dengan assignment tetap, centroid optimal adalah rata-rata cluster:

\[ \mu_k = \frac{1}{|C_k|}\sum_{i \in C_k}x_i. \]

Karena setiap langkah menurunkan objective, algoritma konvergen ke local optimum. Karena itu, inisialisasi menjadi penting, dan nstart di R bukan sekadar parameter komputasi, melainkan strategi mengurangi risiko local optimum buruk.

6.4 K-Means dan Mixture Gaussian

Secara statistik, K-Means dapat diturunkan sebagai kasus batas dari mixture Gaussian dengan:

kovarians sama untuk semua cluster,
kovarians proporsional terhadap matriks identitas,
probabilitas assignment keras.

Interpretasi ini penting karena menunjukkan asumsi implisit K-Means: cluster cenderung berbentuk bulat dan ukuran varian relatif sebanding.

6.5 K-Means++ dan Inisialisasi

Masalah utama K-Means adalah sensitivitas pada centroid awal. K-Means++ memilih centroid awal dengan probabilitas yang proporsional terhadap jarak kuadrat dari centroid yang sudah dipilih. Secara teoritis, ini meningkatkan kualitas solusi awal dan sering mempercepat konvergensi.

6.6 Between dan Within Variation

Jika total sum of squares dapat ditulis sebagai

\[ TSS = WSS + BSS, \]

maka K-Means secara langsung meminimalkan \(WSS\), sehingga secara tidak langsung memaksimalkan separasi antarcluster. Namun hal ini hanya bermakna ketika bentuk cluster yang dicari sejalan dengan asumsi spherical clusters.

6.7 Fuzzy C-Means

Berbeda dari K-Means, FCM mengizinkan observasi memiliki keanggotaan parsial. Objective function-nya adalah

\[ J_m(U,V) = \sum_{i=1}^{n}\sum_{k=1}^{K}u_{ik}^{m}\|x_i-v_k\|^2, \]

dengan kendala

\[ \sum_{k=1}^{K}u_{ik}=1, \qquad u_{ik}\in[0,1]. \]

Parameter \(m>1\) adalah fuzzifier.

6.8 Update Rules FCM

Dengan metode Lagrange, update pusat cluster:

\[ v_k = \frac{\sum_{i=1}^{n}u_{ik}^{m}x_i}{\sum_{i=1}^{n}u_{ik}^{m}}, \]

dan update keanggotaan:

\[ u_{ik} = \left[ \sum_{\ell=1}^{K} \left( \frac{\|x_i-v_k\|}{\|x_i-v_\ell\|} \right)^{\frac{2}{m-1}} \right]^{-1}. \]

Jika \(m \to 1\), keanggotaan cenderung menjadi keras. Jika \(m\) besar, membership menjadi lebih menyebar.

6.9 Interpretasi Statistik dari FCM

FCM cocok ketika batas cluster tidak tajam. Dalam banyak fenomena sosial, ekonomi, dan biologis, observasi memang bisa memiliki karakter campuran. Karena itu, membership fuzzy sering lebih realistis daripada assignment keras.

6.10 Evaluasi Jumlah Cluster

Selain elbow method, dapat dipakai:

average silhouette,
gap statistic,
Calinski-Harabasz index,
validitas substantif.

Mahasiswa S2 sebaiknya tidak memperlakukan satu indeks sebagai “penentu mutlak”, tetapi sebagai alat bantu keputusan.

6.11 Contoh R: K-Means

scaled_iris <- scale(iris[, -5])
kmeans_fit <- kmeans(scaled_iris, centers = 3, nstart = 25)
kmeans_fit$centers

#>   Sepal.Length Sepal.Width Petal.Length Petal.Width
#> 1  -0.05005221 -0.88042696    0.3465767   0.2805873
#> 2   1.13217737  0.08812645    0.9928284   1.0141287
#> 3  -1.01119138  0.85041372   -1.3006301  -1.2507035

table(kmeans_fit$cluster, iris$Species)

#>    
#>     setosa versicolor virginica
#>   1      0         39        14
#>   2      0         11        36
#>   3     50          0         0

6.12 Contoh R: Fuzzy C-Means

library(e1071)
fcm_fit <- cmeans(scaled_iris, centers = 3, m = 2)
head(fcm_fit$membership)

#>              1           2           3
#> [1,] 0.9915694 0.003134665 0.005295911
#> [2,] 0.8289409 0.051838401 0.119220740
#> [3,] 0.9293791 0.023314905 0.047305979
#> [4,] 0.8695180 0.041294487 0.089187559
#> [5,] 0.9739416 0.010010993 0.016047411
#> [6,] 0.8188019 0.079204296 0.101993772

table(apply(fcm_fit$membership, 1, which.max), iris$Species)

#>    
#>     setosa versicolor virginica
#>   1     50          0         0
#>   2      0         11        37
#>   3      0         39        13

6.13 Contoh R: Elbow Sederhana

wss <- numeric(8)
for (k in 1:8) {
  wss[k] <- kmeans(scaled_iris, centers = k, nstart = 20)$tot.withinss
}
plot(1:8, wss, type = "b", pch = 19,
     xlab = "Jumlah Cluster", ylab = "Total Within SS")

6.14 Tugas

Tunjukkan bahwa mean cluster adalah solusi optimal untuk update centroid K-Means.
Bandingkan K-Means dan FCM pada dataset yang sama.
Interpretasikan observasi yang memiliki membership fuzzy paling ambigu.
Jelaskan kapan K-Means tidak cocok secara substantif.

6.15 Ringkasan

K-Means dan FCM sama-sama mempartisi data, tetapi asumsi dan interpretasinya berbeda. K-Means menghasilkan cluster keras berbasis centroid, sedangkan FCM menghasilkan struktur keanggotaan parsial yang lebih fleksibel.

7 Gradient Boosting untuk Analisis Regresi dan Time Series

Pendekatan Gradient Boosting dalam analisis regresi dan time series pada dasarnya dibangun di atas kemampuan untuk menangkap pola yang kompleks dalam data. Tidak seperti model linear klasik yang mengasumsikan hubungan sederhana antara variabel, Gradient Boosting secara iteratif membangun model yang mampu merepresentasikan struktur nonlinear, interaksi variabel, serta dinamika yang tidak stasioner dalam data time series.

Namun, kekuatan ini sekaligus membawa tantangan: model menjadi semakin sulit diinterpretasikan secara langsung. Ketika hubungan antar variabel tidak lagi linear, kita tidak cukup hanya melihat koefisien—kita perlu memahami bagaimana perubahan input memengaruhi output secara lebih fleksibel dan kontekstual.

Di sinilah visualisasi regresi nonlinear menjadi krusial. Baik dalam konteks regresi maupun time series, visualisasi membantu kita:

Mengidentifikasi pola nonlinear yang ditangkap model
Memahami efek marginal dan interaksi antar variabel
Mengevaluasi apakah model benar-benar mencerminkan struktur data
Menghindari kesalahan interpretasi akibat asumsi linear yang tidak terpenuhi

Dengan kata lain, sebelum melangkah lebih jauh dalam interpretasi model berbasis machine learning seperti Gradient Boosting, kita perlu memastikan bahwa kita memahami bentuk hubungan nonlinear yang mendasari data. Oleh karena itu, pembahasan selanjutnya akan difokuskan pada bagaimana memvisualisasikan regresi nonlinear secara efektif, sebagai fondasi untuk interpretasi model yang lebih robust dan bermakna.

Visualisasi Regresi Nonlinear

set.seed(123)
x <- runif(200)
y <- sin(6*pi*x) + rnorm(200,0,0.2)

plot(x,y, main="Data Nonlinear", pch=16, cex=0.6)

library(gbm)
df <- data.frame(x,y)

fit <- gbm(y ~ x,
  data = df,
  distribution = "gaussian",
  n.trees = 2000,
  shrinkage = 0.01,
  interaction.depth = 2,
  verbose = FALSE
)

pred <- predict(fit, df, n.trees=2000)

ord <- order(x)
lines(x[ord], pred[ord], lwd=2)

Visualisasi Time Series

ts_data <- AirPassengers
plot(ts_data, main="Time Series: AirPassengers")

library(gbm)

ap <- as.numeric(ts_data)
lag1 <- ap[-length(ap)]
y <- ap[-1]

df_ts <- data.frame(y=y, lag1=lag1)

fit_ts <- gbm(y ~ lag1,
  data = df_ts,
  distribution = "gaussian",
  n.trees = 1000,
  shrinkage = 0.01,
  interaction.depth = 2,
  verbose = FALSE
)

pred_ts <- predict(fit_ts, df_ts, n.trees=1000)

plot(y, type='l', main="Actual vs Prediction")
lines(pred_ts, col=2)
legend("topleft", legend=c("Actual","Prediction"), col=c(1,2), lty=1)

\[ \begin{array}{ll} x_i & \text{vektor kovariat ke-}i \\ y_i & \text{respon ke-}i \\ F(x) & \text{fungsi prediksi} \\ L(y,F) & \text{fungsi loss} \\ g_i & \text{gradien loss terhadap } F(x_i) \\ h_i & \text{Hessian (turunan kedua)} \\ r_{im} & \text{pseudo-residual pada iterasi } m \\ h_m(x) & \text{weak learner ke-}m \\ \gamma_{jm} & \text{update pada region } j \\ \nu & \text{learning rate} \end{array} \]

Gradient Boosting untuk Regresi dan Time Series: Derivasi Formal, Newton Boosting, Konvergensi, dan Interpretasi Statistik

7.1 Pendahuluan

Gradient boosting adalah salah satu kerangka pembelajaran statistik yang paling penting karena ia menyatukan tiga gagasan besar ke dalam satu konstruksi yang koheren: empirical risk minimization, ekspansi aditif bertahap, dan optimisasi gradien dalam ruang fungsi. Dalam praktik terapan, metode ini sering direduksi menjadi slogan yang terlalu sederhana, misalnya “membangun banyak pohon kecil” atau “memperbaiki residual berulang kali”. Kedua pernyataan tersebut berguna sebagai intuisi awal, tetapi tidak cukup untuk tujuan akademik tingkat lanjut. Pada level pascasarjana, gradient boosting perlu dipahami sebagai prosedur optimisasi yang eksplisit: kita tidak sedang sekadar menumpuk model, melainkan membangun aproksimasi fungsi prediksi melalui lintasan penurunan risiko empiris.

Secara historis, kontribusi Friedman sangat penting karena ia menunjukkan bahwa boosting dapat dirumuskan sebagai steepest descent dalam ruang fungsi. Perspektif ini menjelaskan mengapa boosting sangat umum: selama suatu fungsi loss dapat didiferensiasikan terhadap prediktor, maka prosedur boosting dapat dirancang. Dengan demikian, kerangka ini melampaui regresi Gaussian dan dapat diperluas ke klasifikasi, regresi robust, kuantil, survival, dan bentuk loss lain. Dalam praktik komputasi modern, generalisasi lebih lanjut menghasilkan apa yang sekarang dikenal sebagai second-order boosting atau Newton boosting, yaitu pembaruan yang tidak hanya memanfaatkan gradien, tetapi juga kelengkungan lokal melalui Hessian. Banyak perangkat lunak cepat modern mengandalkan prinsip ini karena ia mempercepat konvergensi dan menstabilkan pembelajaran.

Bab ini disusun untuk memenuhi kebutuhan yang lebih formal: bukan ringkasan populer, melainkan derivasi matematis yang bertahap dan lengkap. Fokus utama bab adalah lima hal. Pertama, menjelaskan boosting dalam kerangka empirical risk minimization. Kedua, membuktikan hubungan formal antara boosting dan gradient descent dalam ruang fungsi. Ketiga, menurunkan pembaruan Newton per region ketika weak learner berupa pohon regresi, termasuk bentuk berbobot Hessian. Keempat, mendiskusikan konvergensi dan regularisasi. Kelima, memberikan interpretasi statistik agar boosting tidak hanya dipahami sebagai algoritma komputasi, tetapi juga sebagai prosedur estimasi fungsi yang memiliki implikasi inferensial dan prediktif.

Selain regresi tabular, bab ini juga menempatkan gradient boosting dalam kerangka time series. Ini penting karena dalam banyak aplikasi modern, forecasting tidak lagi diperlakukan secara eksklusif dalam keluarga ARIMA atau state space tradisional. Sebaliknya, deret waktu dapat diubah ke bentuk supervised learning melalui fitur lag, musiman, tren, dan kovariat eksternal. Dalam konteks ini, gradient boosting menjadi alat yang sangat fleksibel. Namun, fleksibilitas tersebut datang dengan syarat metodologis: validasi harus menghormati urutan waktu, dan interpretasi tidak boleh mengabaikan sifat dinamis data.

Secara keseluruhan, tujuan bab ini adalah mengangkat boosting dari status “alat prediksi populer” menjadi objek kajian statistik yang utuh. Dengan pemahaman seperti ini, pembaca tidak hanya dapat menjalankan fungsi perangkat lunak, tetapi juga menjelaskan mengapa algoritma itu bekerja, bagaimana ia harus diatur, kapan ia sesuai, dan batas-batas inferensial apa yang harus dijaga.

Ilustrasi Perhitungan Manual Sederhana

Agar proses boosting tidak terasa terlalu abstrak, bagian ini memberikan ilustrasi manual yang sangat sederhana untuk kasus regresi dengan squared error loss dan weak learner berupa decision stump satu split. Tujuannya bukan menghasilkan model terbaik, tetapi memperlihatkan secara konkret bagaimana residual dihitung, bagaimana weak learner dipasang, dan bagaimana model diperbarui langkah demi langkah.

Misalkan tersedia empat observasi berikut.

\[ \begin{array}{c|c|c} i & x_i & y_i \\ \hline 1 & 1 & 2 \\ 2 & 2 & 4 \\ 3 & 3 & 5 \\ 4 & 4 & 8 \\ \end{array} \]

Kita gunakan loss Gaussian: \[ L(y,F)=\frac{1}{2}(y-F)^2. \]

Untuk loss ini, pseudo-residual sama dengan residual biasa: \[ r_{im}=y_i-F_{m-1}(x_i). \]

Kita juga gunakan learning rate \[ \nu = 0.5. \]

Langkah awal: inisialisasi model

Untuk squared error loss, konstanta awal optimal adalah rata-rata respons: \[ F_0(x)=\arg\min_c \sum_{i=1}^{n}\frac{1}{2}(y_i-c)^2 = \bar y. \]

Hitung rata-ratanya: \[ \bar y=\frac{2+4+5+8}{4}=\frac{19}{4}=4.75. \]

Jadi pada iterasi awal, \[ F_0(x)=4.75 \qquad \text{untuk semua } x. \]

Prediksi awal dan residual:

\[ \begin{array}{c|c|c|c} i & y_i & F_0(x_i) & r_{i1}=y_i-F_0(x_i) \\ \hline 1 & 2 & 4.75 & -2.75 \\ 2 & 4 & 4.75 & -0.75 \\ 3 & 5 & 4.75 & 0.25 \\ 4 & 8 & 4.75 & 3.25 \\ \end{array} \]

Secara intuitif, dua observasi pertama diprediksi terlalu tinggi, sedangkan observasi terakhir diprediksi terlalu rendah.

Iterasi pertama: fit weak learner ke residual

Sekarang kita fit satu stump ke residual. Karena \(x\) berurutan, split alami yang masuk akal adalah memisahkan data menjadi dua region: \[ R_{11}=\{x: x \le 2.5\}, \qquad R_{21}=\{x: x > 2.5\}. \]

Dengan split ini:

Region kiri berisi observasi 1 dan 2
Region kanan berisi observasi 3 dan 4

Untuk squared error, nilai optimal per region adalah rata-rata residual pada region tersebut.

Region kiri

Residual di kiri: \[ -2.75,\; -0.75 \]

Rata-ratanya: \[ \gamma_{11}=\frac{-2.75+(-0.75)}{2}=-1.75. \]

Region kanan

Residual di kanan: \[ 0.25,\; 3.25 \]

Rata-ratanya: \[ \gamma_{21}=\frac{0.25+3.25}{2}=1.75. \]

Maka weak learner pertama adalah \[ h_1(x)= \begin{cases} -1.75, & x \le 2.5,\\ 1.75, & x > 2.5. \end{cases} \]

Update model iterasi pertama

Rumus update: \[ F_1(x)=F_0(x)+\nu h_1(x). \]

Karena \(\nu=0.5\), maka: \[ F_1(x)= \begin{cases} 4.75 + 0.5(-1.75)=3.875, & x \le 2.5,\\ 4.75 + 0.5(1.75)=5.625, & x > 2.5. \end{cases} \]

Prediksi baru menjadi:

\[ \begin{array}{c|c|c} i & x_i & F_1(x_i) \\ \hline 1 & 1 & 3.875 \\ 2 & 2 & 3.875 \\ 3 & 3 & 5.625 \\ 4 & 4 & 5.625 \\ \end{array} \]

Residual baru: \[ \begin{array}{c|c|c|c} i & y_i & F_1(x_i) & r_{i2}=y_i-F_1(x_i) \\ \hline 1 & 2 & 3.875 & -1.875 \\ 2 & 4 & 3.875 & 0.125 \\ 3 & 5 & 5.625 & -0.625 \\ 4 & 8 & 5.625 & 2.375 \\ \end{array} \]

Perhatikan bahwa besar residual secara umum sudah mengecil dibanding iterasi awal. Model mulai bergerak ke arah yang benar, tetapi belum selesai.

Iterasi kedua: fit weak learner lagi

Sekarang kita ulangi prosedur yang sama, tetapi menggunakan residual baru: \[ -1.875,\; 0.125,\; -0.625,\; 2.375. \]

Untuk ilustrasi sederhana, kita gunakan split yang sama: \[ R_{12}=\{x: x \le 2.5\}, \qquad R_{22}=\{x: x > 2.5\}. \]

Region kiri

Residual kiri: \[ -1.875,\; 0.125 \]

Rata-ratanya: \[ \gamma_{12}=\frac{-1.875+0.125}{2}=-0.875. \]

Region kanan

Residual kanan: \[ -0.625,\; 2.375 \]

Rata-ratanya: \[ \gamma_{22}=\frac{-0.625+2.375}{2}=0.875. \]

Weak learner kedua: \[ h_2(x)= \begin{cases} -0.875, & x \le 2.5,\\ 0.875, & x > 2.5. \end{cases} \]

Update model iterasi kedua

\[ F_2(x)=F_1(x)+\nu h_2(x). \]

Dengan \(\nu=0.5\): \[ F_2(x)= \begin{cases} 3.875 + 0.5(-0.875)=3.4375, & x \le 2.5,\\ 5.625 + 0.5(0.875)=6.0625, & x > 2.5. \end{cases} \]

Prediksi baru: \[ \begin{array}{c|c|c} i & x_i & F_2(x_i) \\ \hline 1 & 1 & 3.4375 \\ 2 & 2 & 3.4375 \\ 3 & 3 & 6.0625 \\ 4 & 4 & 6.0625 \\ \end{array} \]

Residual baru: \[ \begin{array}{c|c|c|c} i & y_i & F_2(x_i) & r_{i3}=y_i-F_2(x_i) \\ \hline 1 & 2 & 3.4375 & -1.4375 \\ 2 & 4 & 3.4375 & 0.5625 \\ 3 & 5 & 6.0625 & -1.0625 \\ 4 & 8 & 6.0625 & 1.9375 \\ \end{array} \]

Sekali lagi residual berubah dan secara umum mengecil dibanding awal. Jika iterasi dilanjutkan, model akan terus memperbaiki dirinya sedikit demi sedikit.

Bentuk model setelah dua iterasi

Karena boosting adalah model aditif, setelah dua iterasi kita punya: \[ F_2(x)=F_0(x)+\nu h_1(x)+\nu h_2(x). \]

Substitusikan nilai-nilainya: \[ F_2(x)=4.75 + 0.5h_1(x)+0.5h_2(x). \]

Ini adalah contoh sangat konkret bahwa boosting bukan “satu pohon besar”, melainkan penjumlahan dari beberapa weak learners kecil.

Apa yang dipelajari dari ilustrasi manual ini

Ilustrasi kecil ini memberi beberapa pelajaran penting.

Pertama, pada loss Gaussian, boosting memang sama dengan fitting residual secara bertahap. Kedua, learning rate memperkecil kontribusi setiap learner sehingga model tidak langsung melompat terlalu jauh. Ketiga, model akhir adalah penjumlahan aditif dari banyak koreksi lokal. Keempat, bahkan dengan split yang sangat sederhana, boosting tetap bisa memperbaiki prediksi secara progresif.

Dalam praktik nyata, stump atau pohon pada setiap iterasi tidak harus memakai split yang sama. Algoritma akan memilih split yang paling baik untuk residual atau gradien saat itu. Namun ilustrasi manual ini cukup untuk memperlihatkan mekanisme inti algoritma secara transparan.

7.2 Empirical Risk Minimization dan Ruang Fungsi

Misalkan tersedia data pelatihan \[ \mathcal{D}=\{(x_i,y_i)\}_{i=1}^{n}, \] dengan \(x_i \in \mathbb{R}^p\) dan respons \(y_i \in \mathcal{Y}\). Dalam regresi, biasanya \(\mathcal{Y}=\mathbb{R}\), sedangkan pada klasifikasi biner dapat dipakai \(\mathcal{Y}=\{-1,+1\}\) atau \(\{0,1\}\).

Tujuan dasar pembelajaran statistik adalah membangun fungsi prediksi \[ F:\mathbb{R}^p \to \mathbb{R} \] yang meminimalkan risiko populasi \[ R(F)=\mathbb{E}[L(Y,F(X))]. \] Karena distribusi populasi umumnya tidak diketahui, kita menggantinya dengan risiko empiris \[ \widehat{R}(F)=\sum_{i=1}^{n}L(y_i,F(x_i)). \] Masalah estimasi kemudian ditulis sebagai \[ \widehat{F}=\arg\min_{F \in \mathcal{F}} \widehat{R}(F) =\arg\min_{F \in \mathcal{F}} \sum_{i=1}^{n}L(y_i,F(x_i)), \] dengan \(\mathcal{F}\) suatu kelas fungsi kandidat.

Dalam model parametrik klasik, kita akan menulis \(F(x)=x^\top\beta\) lalu mengoptimalkan terhadap \(\beta\). Dalam boosting, objek yang dioptimalkan adalah fungsi itu sendiri. Inilah perbedaan konseptual yang sangat penting. Karena \(\mathcal{F}\) dapat sangat besar atau bahkan tak hingga dimensi, optimisasi tidak lagi dilakukan dengan menghitung gradien terhadap vektor parameter yang kecil, tetapi terhadap nilai fungsi pada titik-titik pelatihan. Oleh sebab itu, boosting hidup di ranah function space optimization.

Ada dua implikasi metodologis dari formulasi ini. Pertama, pilihan fungsi loss menjadi sangat sentral. Fungsi loss tidak lagi sekadar konsekuensi dari asumsi distribusi, tetapi merupakan definisi eksplisit tentang apa yang dianggap sebagai kesalahan. Kedua, kelas fungsi yang dipilih melalui weak learners akan membatasi arah pembaruan yang dapat diambil selama optimisasi. Jadi, boosting selalu merupakan interaksi antara objective function dan function class.

Jika kita menulis masalahnya secara naif sebagai minimisasi langsung di ruang fungsi besar, maka optimisasi menjadi sangat sulit. Friedman menunjukkan bahwa persoalan ini dapat ditangani dengan strategi bertahap: mulai dari fungsi sederhana, lalu tambahkan koreksi sedikit demi sedikit dalam arah yang menurunkan risiko paling cepat. Inilah jembatan menuju gradient descent dalam ruang fungsi.

7.3 Ekspansi Aditif dan Forward Stagewise Modeling

Boosting membangun model secara bertahap sebagai ekspansi aditif \[ F_M(x)=F_0(x)+\sum_{m=1}^{M}\nu \rho_m h_m(x), \] dengan: - \(F_0(x)\) adalah inisialisasi, - \(h_m(x)\) weak learner pada iterasi ke-\(m\), - \(\rho_m\) besar langkah atau koefisien pembaruan, - \(\nu \in (0,1]\) learning rate atau shrinkage parameter.

Bentuk ini sering disederhanakan menjadi \[ F_M(x)=\sum_{m=0}^{M}\rho_m h_m(x), \] tetapi untuk analisis boosting, lebih informatif mempertahankan \(F_0\) dan faktor shrinkage \(\nu\) secara eksplisit.

Pendekatan ini disebut forward stagewise additive modeling karena model dibangun dengan menambahkan satu komponen pada satu waktu. Pada iterasi ke-\(m\), semua komponen sebelumnya dianggap tetap, dan kita mencari learner baru yang paling efektif menurunkan risiko saat ini. Tidak ada optimisasi simultan atas seluruh ekspansi. Prosesnya sepenuhnya bertahap.

Mengapa pendekatan ini masuk akal? Karena optimisasi simultan di atas kelas fungsi aditif yang kaya sering sangat sulit. Dengan membangun model selangkah demi selangkah, kita memperoleh prosedur komputasi yang lebih sederhana, dan secara tidak langsung juga memperoleh regularisasi. Setiap langkah hanya menambah kompleksitas sedikit, sehingga lintasan pembelajaran dapat dikontrol.

Interpretasi statistiknya juga penting. Ekspansi aditif berarti model akhir adalah agregat dari banyak komponen lokal. Jika weak learner berupa pohon dangkal, maka fungsi prediksi akhir dapat dipandang sebagai penjumlahan banyak aturan lokal sederhana. Masing-masing aturan mungkin lemah, tetapi gabungannya dapat sangat fleksibel. Dalam konteks ini, kekuatan model tidak datang dari satu struktur global tunggal, melainkan dari akumulasi koreksi lokal terhadap kesalahan sebelumnya.

7.4 Functional Gradient Descent:

Intuisi: Boosting = gradient descent, tapi di ruang fungsi. Kita tidak update parameter, tapi update fungsi secara bertahap. Hubungan Formal dengan Gradient Descent

Bagian ini adalah jantung teoritis boosting. Kita ingin membuktikan bahwa boosting merupakan analog gradient descent ketika objek yang dioptimalkan adalah fungsi, bukan vektor parameter berhingga dimensi.

7.4.1 Gradient descent biasa

Misalkan kita memiliki fungsi objektif \(Q(\theta)\) dengan parameter \(\theta \in \mathbb{R}^d\). Langkah gradient descent adalah \[ \theta_m = \theta_{m-1} - \nu \nabla Q(\theta_{m-1}). \] Interpretasi langkah ini adalah bergerak ke arah penurunan tercepat (steepest descent) dari fungsi objektif.

7.4.2 Analog dalam ruang fungsi

Sekarang objektif kita adalah \[ \widehat{R}(F)=\sum_{i=1}^{n}L(y_i,F(x_i)). \] Karena fungsi hanya muncul melalui nilai-nilainya pada titik pelatihan, kita dapat memandang \(\widehat{R}\) sebagai fungsi dari vektor \[ \big(F(x_1),F(x_2),\ldots,F(x_n)\big). \] Dengan demikian, gradien empirisnya terhadap nilai fungsi pada sampel adalah \[ \nabla_F \widehat{R}(F)= \left( \frac{\partial L(y_1,F(x_1))}{\partial F(x_1)}, \ldots, \frac{\partial L(y_n,F(x_n))}{\partial F(x_n)} \right)^\top. \]

Definisikan pseudo-residual pada iterasi ke-\(m\) sebagai \[ r_{im} = -\left[ \frac{\partial L(y_i,F(x_i))}{\partial F(x_i)} \right]_{F=F_{m-1}}. \] Vektor residual semu \[ \mathbf{r}_m=(r_{1m},\ldots,r_{nm})^\top \] adalah negatif gradien risiko empiris pada fungsi saat ini. Dengan kata lain, jika kita dapat bergerak bebas di ruang fungsi, maka pembaruan steepest descent ideal adalah \[ F_m(x_i)=F_{m-1}(x_i)+\nu r_{im}, \qquad i=1,\ldots,n. \]

Namun boosting tidak bebas memilih sembarang koreksi. Ia dibatasi oleh kelas weak learners \(\mathcal{H}\). Oleh karena itu, boosting memilih \[ h_m \approx \mathbf{r}_m \] dalam arti terbaik yang dapat direpresentasikan oleh weak learner. Jika weak learner dipasang dengan least squares terhadap pseudo-residual, maka kita menyelesaikan \[ h_m = \arg\min_{h \in \mathcal{H}} \sum_{i=1}^{n}\big(r_{im}-h(x_i)\big)^2. \]

7.4.3 Pernyataan formal

7.4.4 Bukti

Untuk notasi singkat, definisikan \[ g_i(F)=\frac{\partial L(y_i,F(x_i))}{\partial F(x_i)}. \] Negatif gradien adalah \(-g_i(F_{m-1})=r_{im}\). Jika ruang fungsi tak dibatasi, steepest descent satu langkah akan memberi koreksi persis \(r_{im}\) pada tiap titik pelatihan. Karena kita dibatasi oleh \(\mathcal{H}\), maka kita pilih fungsi \(h_m\) yang meminimalkan jarak kuadrat \[ \sum_{i=1}^{n}\left(r_{im}-h(x_i)\right)^2. \] Ini setara dengan memproyeksikan negatif gradien ke span atau subkelas fungsi yang tersedia. Pembaruan \[ F_m=F_{m-1}+\nu h_m \] berarti kita bergerak ke arah proyeksi negatif gradien dengan ukuran langkah \(\nu\). Karena \(h_m\) hanya aproksimasi arah gradien ideal, prosedur ini adalah gradient descent terproyeksi, bukan gradient descent penuh. Selesai.

7.4.5 Makna pembuktian

Pembuktian di atas sangat penting karena menjelaskan bahwa boosting bukan heuristik ad hoc. Ia adalah prosedur optimisasi yang sah, dengan satu nuansa: arah turun tidak diambil secara bebas, melainkan diproyeksikan ke kelas learner yang kita izinkan. Dengan demikian, kualitas learner menentukan kualitas arah pembaruan. Jika learner terlalu miskin, pembelajaran lambat. Jika learner terlalu kaya, pembelajaran dapat menjadi agresif dan rentan overfitting.

7.5 Derivasi Step-by-Step untuk Loss Gaussian

Pertimbangkan loss Gaussian \[ L(y,F)=\frac{1}{2}(y-F)^2. \] Kita turunkan pseudo-residual secara lengkap.

Turunan terhadap prediktor: \[ \frac{\partial L(y,F)}{\partial F} = \frac{\partial}{\partial F}\left[\frac{1}{2}(y-F)^2\right] = \frac{1}{2}\cdot 2(y-F)(-1) = -(y-F). \] Maka, \[ r_{im} = -\left[\frac{\partial L(y_i,F(x_i))}{\partial F(x_i)}\right]_{F=F_{m-1}} = y_i-F_{m-1}(x_i). \]

Jadi, pada loss Gaussian, pseudo-residual sama dengan residual biasa. Inilah alasan mengapa boosting regresi sering tampak seperti prosedur iterative residual fitting. Namun penting dicatat bahwa identitas ini spesifik untuk squared error; pada loss lain, pseudo-residual tidak lagi sama dengan residual klasik.

Jika weak learner adalah pohon regresi, kita menyesuaikan pohon terhadap pasangan data \[ \{(x_i,r_{im})\}_{i=1}^{n}. \] Misalkan pohon menghasilkan partisi \[ R_{1m}, R_{2m}, \ldots, R_{J_m m}. \] Pada masing-masing region, kita mencari nilai konstan \(\gamma_{jm}\) yang meminimalkan \[ \sum_{x_i \in R_{jm}} \frac{1}{2}\left(y_i - \big(F_{m-1}(x_i)+\gamma\big)\right)^2. \] Turunkan terhadap \(\gamma\): \[ \frac{\partial}{\partial \gamma} \sum_{x_i \in R_{jm}} \frac{1}{2} \left(y_i - F_{m-1}(x_i)-\gamma\right)^2 = -\sum_{x_i \in R_{jm}} \left(y_i - F_{m-1}(x_i)-\gamma\right). \] Setarakan dengan nol: \[ \sum_{x_i \in R_{jm}} \left(y_i - F_{m-1}(x_i)-\gamma\right)=0. \] Sehingga \[ |R_{jm}|\gamma = \sum_{x_i \in R_{jm}} \left(y_i - F_{m-1}(x_i)\right) = \sum_{x_i \in R_{jm}} r_{im}, \] dan akhirnya \[ \gamma_{jm} = \frac{1}{|R_{jm}|} \sum_{x_i \in R_{jm}} r_{im}. \]

Dengan learning rate \(\nu\), pembaruan model adalah \[ F_m(x) = F_{m-1}(x) + \nu \sum_{j=1}^{J_m} \gamma_{jm}\mathbf{1}(x \in R_{jm}). \]

Derivasi ini menunjukkan tiga lapis struktur boosting Gaussian:

hitung residual,
aproksimasi residual dengan pohon,
koreksi prediktor per region dengan rata-rata residual lokal.

7.6 Derivasi Step-by-Step untuk Logistic Loss

Sekarang anggap klasifikasi biner dengan pengkodean \(y_i \in \{-1,+1\}\) dan fungsi loss \[ L(y,F)=\log\big(1+\exp(-2yF)\big). \] Loss ini berhubungan dengan model logit aditif.

Turunkan terhadap \(F\): \[ \frac{\partial L(y,F)}{\partial F} = \frac{1}{1+\exp(-2yF)}\cdot \exp(-2yF)\cdot (-2y). \] Karena \[ \frac{\exp(-2yF)}{1+\exp(-2yF)} = \frac{1}{1+\exp(2yF)}, \] maka \[ \frac{\partial L(y,F)}{\partial F} = -\frac{2y}{1+\exp(2yF)}. \] Akibatnya pseudo-residual menjadi \[ r_{im} = -\left[\frac{\partial L(y_i,F(x_i))}{\partial F(x_i)}\right]_{F=F_{m-1}} = \frac{2y_i}{1+\exp(2y_iF_{m-1}(x_i))}. \]

Bentuk ini menjelaskan mengapa residual semu pada klasifikasi tidak lagi sama dengan selisih observasi dan prediksi seperti pada regresi Gaussian. Besarnya residual semu dipengaruhi oleh keyakinan model saat ini. Jika pengamatan sudah diprediksi dengan margin besar dan benar, maka residual semu kecil. Jika salah atau margin lemah, residual semu besar. Dengan kata lain, boosting lebih fokus pada titik-titik yang saat ini masih sulit.

Untuk pembaruan per region, biasanya solusi tertutup tidak sesederhana kasus Gaussian. Karena itu, implementasi sering memakai line search numerik, atau pendekatan second-order berbasis Newton seperti yang akan dibahas pada bagian berikutnya.

7.7 Huber Loss dan Regresi Robust

Huber loss didefinisikan sebagai \[ L_\delta(y,F) = \begin{cases} \frac{1}{2}(y-F)^2, & |y-F|\le \delta,\\[4pt] \delta |y-F| - \frac{1}{2}\delta^2, & |y-F|>\delta. \end{cases} \] Tuliskan residual \(e=y-F\). Maka \[ L_\delta(e) = \begin{cases} \frac{1}{2}e^2, & |e|\le \delta,\\ \delta |e| - \frac{1}{2}\delta^2, & |e|>\delta. \end{cases} \]

Turunan terhadap \(F\) harus memperhatikan bahwa \(e=y-F\), sehingga \(\partial e/\partial F=-1\).

Jika \(|e|\le \delta\), \[ \frac{\partial L_\delta}{\partial F} = \frac{\partial}{\partial F}\left(\frac{1}{2}e^2\right) = e(-1) = -(y-F). \]

Jika \(e>\delta\), maka \(|e|=e\) dan \[ L_\delta=\delta e - \frac{1}{2}\delta^2, \] sehingga \[ \frac{\partial L_\delta}{\partial F} = \delta \frac{\partial e}{\partial F} = -\delta. \]

Jika \(e<-\delta\), maka \(|e|=-e\) dan \[ L_\delta=-\delta e - \frac{1}{2}\delta^2, \] sehingga \[ \frac{\partial L_\delta}{\partial F} = -\delta \frac{\partial e}{\partial F} = \delta. \]

Maka pseudo-residual adalah \[ r_{im} = \begin{cases} y_i-F_{m-1}(x_i), & |y_i-F_{m-1}(x_i)|\le \delta,\\[4pt] \delta, & y_i-F_{m-1}(x_i)>\delta,\\[4pt] -\delta, & y_i-F_{m-1}(x_i)<-\delta. \end{cases} \]

Bentuk ini memperlihatkan sifat robust Huber loss: observasi dengan error sangat besar tidak lagi menghasilkan residual semu yang terus tumbuh secara linear, tetapi dipotong pada \(\pm \delta\). Dengan demikian, outlier tidak mendominasi arah pembelajaran. Dari sudut statistik, ini setara dengan mengurangi sensitivitas estimator terhadap pencilan.

7.8 Representasi Pohon dan Pembaruan Per Region

Dalam boosting berbasis tree, weak learner pada iterasi ke-\(m\) ditulis sebagai \[ h_m(x)=\sum_{j=1}^{J_m} b_{jm}\mathbf{1}(x\in R_{jm}), \] dengan \(R_{jm}\) region terminal dan \(b_{jm}\) prediksi konstan pada region tersebut. Dalam formulasi boosting yang lebih tepat, yang kita update sebenarnya bukan pohon mentah, melainkan langkah optimal per region: \[ F_m(x)=F_{m-1}(x)+\nu \sum_{j=1}^{J_m}\gamma_{jm}\mathbf{1}(x\in R_{jm}). \] Jadi, struktur pohon menentukan partisi ruang, sedangkan \(\gamma_{jm}\) menentukan besar koreksi lokal.

Pada level komputasi, ada dua tahap: 1. bangun pohon menggunakan pseudo-residual atau informasi gradien, 2. hitung pembaruan optimal per region berdasarkan loss asli.

Pembedaan dua tahap ini penting. Banyak penjelasan populer langsung menganggap nilai terminal pohon adalah pembaruan akhir, padahal dalam formulasi umum, setelah region ditetapkan kita masih melakukan optimisasi lokal untuk menentukan \(\gamma_{jm}\).

7.9 Newton Boosting:

Catatan penting: Newton boosting memakai informasi curvature (Hessian), sehingga langkah lebih adaptif dibanding gradient boosting biasa. Motivasi dan Derivasi Umum

Gradient boosting klasik memakai aproksimasi first-order: arah pembaruan ditentukan oleh gradien loss. Newton boosting melangkah lebih jauh dengan memakai informasi second-order, yaitu Hessian lokal. Ide dasarnya sangat mirip dengan metode Newton-Raphson dalam optimisasi parametrik.

Misalkan untuk satu observasi kita definisikan \[ g_i = \left[ \frac{\partial L(y_i,F(x_i))}{\partial F(x_i)} \right]_{F=F_{m-1}}, \qquad h_i = \left[ \frac{\partial^2 L(y_i,F(x_i))}{\partial F(x_i)^2} \right]_{F=F_{m-1}}. \] Untuk pembaruan lokal \(\eta_i\), lakukan ekspansi Taylor orde dua: \[ L(y_i,F_{m-1}(x_i)+\eta_i) \approx L(y_i,F_{m-1}(x_i)) + g_i\eta_i + \frac{1}{2}h_i\eta_i^2. \] Menjumlahkan atas semua \(i\) menghasilkan aproksimasi kuadratik terhadap risiko empiris lokal.

Jika kita bebas memilih \(\eta_i\) per observasi, minimisasi sederhana memberi \[ \eta_i^*=-\frac{g_i}{h_i}, \] selama \(h_i>0\). Ini adalah langkah Newton per titik data.

Namun dalam boosting berbasis pohon, kita tidak memilih pembaruan bebas per observasi. Kita membatasi bentuk pembaruan menjadi konstan per region: \[ \eta(x)=\sum_{j=1}^{J_m} \gamma_{jm}\mathbf{1}(x\in R_{jm}). \] Masukkan bentuk ini ke aproksimasi Taylor. Kontribusi pada region \(R_{jm}\) menjadi \[ \sum_{x_i \in R_{jm}} \left( g_i \gamma_{jm} + \frac{1}{2}h_i\gamma_{jm}^2 \right). \] Turunkan terhadap \(\gamma_{jm}\): \[ \frac{\partial}{\partial \gamma_{jm}} \sum_{x_i \in R_{jm}} \left( g_i \gamma_{jm} + \frac{1}{2}h_i\gamma_{jm}^2 \right) = \sum_{x_i \in R_{jm}} g_i + \gamma_{jm} \sum_{x_i \in R_{jm}} h_i. \] Setarakan dengan nol: \[ \sum_{x_i \in R_{jm}} g_i + \gamma_{jm} \sum_{x_i \in R_{jm}} h_i = 0. \] Sehingga \[ \gamma_{jm} = -\frac{\sum_{x_i \in R_{jm}} g_i} {\sum_{x_i \in R_{jm}} h_i}. \]

Inilah bentuk Newton update per region berbobot Hessian. Ini adalah salah satu formula paling penting dalam boosting modern.

7.9.1 Interpretasi

Bentuk di atas menunjukkan bahwa pembaruan lokal bukan sekadar rata-rata gradien, tetapi rasio antara akumulasi gradien dan akumulasi kelengkungan. Jika Hessian besar, loss sangat melengkung di region itu, sehingga langkah optimal mengecil. Jika Hessian kecil, langkah dapat lebih besar. Dengan demikian, Newton boosting secara adaptif menyesuaikan ukuran langkah lokal terhadap geometri loss.

7.10 Newton Boosting untuk Logistic Loss

Untuk logistic loss \[ L(y,F)=\log(1+\exp(-2yF)), \quad y\in\{-1,+1\}, \] kita sudah memperoleh gradien \[ g_i = -\frac{2y_i}{1+\exp(2y_iF_i)}, \qquad F_i=F_{m-1}(x_i). \] Sekarang hitung Hessian.

Tuliskan \[ g_i=-2y_i(1+\exp(2y_iF_i))^{-1}. \] Turunkan terhadap \(F_i\): \[ h_i = -2y_i \cdot \left[-(1+\exp(2y_iF_i))^{-2}\right] \cdot \exp(2y_iF_i)\cdot 2y_i. \] Karena \(y_i^2=1\), \[ h_i = \frac{4\exp(2y_iF_i)}{(1+\exp(2y_iF_i))^2}. \] Maka Newton update per region adalah \[ \gamma_{jm} = -\frac{\sum_{x_i \in R_{jm}} g_i} {\sum_{x_i \in R_{jm}} h_i} = \frac{\sum_{x_i \in R_{jm}} \frac{2y_i}{1+\exp(2y_iF_i)}} {\sum_{x_i \in R_{jm}} \frac{4\exp(2y_iF_i)}{(1+\exp(2y_iF_i))^2}}. \]

Dalam banyak implementasi praktis, bentuk ini ekuivalen dengan IRLS (iteratively reweighted least squares) lokal: gradien bertindak sebagai pseudo-response, Hessian sebagai bobot. Di sinilah boosting second-order mempunyai hubungan erat dengan metode Newton-Raphson untuk GLM.

7.11 Pembuktian Hubungan Boosting dan Gradient Descent Secara Formal

Bagian ini merumuskan hubungan tersebut sedikit lebih formal dalam bahasa optimisasi.

7.11.1 Definisi arah steepest descent

Pada ruang Euclidean, arah steepest descent dari fungsi terdiferensialkan adalah negatif gradien. Dalam ruang fungsi empiris yang dinilai pada titik sampel, kita dapat memandang setiap fungsi \(f\) sebagai vektor \[ \mathbf{f}=(f(x_1),\ldots,f(x_n))^\top. \] Dengan inner product empiris \[ \langle f,g\rangle_n = \sum_{i=1}^{n} f(x_i)g(x_i), \] arah steepest descent untuk \(\widehat{R}(F)\) pada \(F_{m-1}\) adalah fungsi \(\delta\) yang meminimalkan \[ D\widehat{R}(F_{m-1})[\delta] \quad \text{dengan kendala} \quad \|\delta\|_n = 1, \] di mana \(D\widehat{R}(F)[\delta]\) adalah directional derivative. Dari teori dasar, solusi kendala ini proporsional terhadap negatif gradien empiris.

7.11.2 Proyeksi ke kelas weak learners

Karena boosting membatasi arah ke \(\mathcal{H}\), ia menyelesaikan aproksimasi: \[ h_m = \arg\min_{h\in\mathcal{H}} \left\| h - \big(-\nabla_F \widehat{R}(F_{m-1})\big) \right\|_n^2. \] Jika \(\mathcal{H}\) kaya, maka \(h_m\) mendekati arah steepest descent aktual. Jika \(\mathcal{H}\) sempit, maka pembaruan adalah versi terbatas dari steepest descent.

7.11.3 Konsekuensi

Dengan learning rate \(\nu\), boosting memberi \[ F_m = F_{m-1} + \nu h_m. \] Ini identik dengan projected gradient descent. Maka semua intuisi dasar gradient descent tetap berlaku: - learning rate terlalu besar dapat menyebabkan osilasi atau overfitting, - learning rate kecil memberi lintasan lebih stabil, - kualitas aproksimasi arah turun mempengaruhi efisiensi penurunan loss.

Secara ringkas, boosting adalah gradient descent di ruang fungsi, tetapi dibatasi oleh kelas arah yang dapat direpresentasikan weak learner.

7.12 Regularisasi dan Mengapa Shrinkage Bekerja

Pada pandangan pertama, boosting tampak paradoksal. Ia dapat menghasilkan model sangat kompleks, dengan ratusan atau ribuan pohon, tetapi justru sering memiliki generalisasi yang baik. Kuncinya adalah regularisasi implisit.

Learning rate \(\nu\) memaksa setiap pembaruan menjadi kecil: \[ F_m = F_{m-1} + \nu \Delta_m. \] Jika \(\nu\) kecil, model bergerak perlahan sepanjang lintasan optimisasi. Hal ini memiliki dua akibat: 1. model tidak langsung menyesuaikan noise, 2. jumlah iterasi optimal biasanya meningkat.

Dari sudut pandang jalur regularisasi, boosting dengan learning rate kecil dapat dipahami sebagai prosedur pathwise estimation: kompleksitas model tumbuh sedikit demi sedikit. Ini mirip semangat regularisasi penalti, walaupun mekanismenya berbeda.

Regularisasi lain muncul dari: - kedalaman pohon, - minimum observasi per node, - subsampling (stochastic boosting), - early stopping.

Pada pohon dangkal, setiap pembaruan hanya menangkap struktur lokal yang sederhana. Jika pembaruan kecil tetapi banyak, model bisa fleksibel namun tetap terkontrol. Inilah alasan mengapa kombinasi small trees + small learning rate + many iterations sering sangat efektif.

7.13 Konvergensi: Diskusi Teoretis

Analisis konvergensi boosting penuh dapat sangat teknis karena bergantung pada fungsi loss, kelas learner, dan strategi line search. Namun beberapa gagasan kunci dapat dirangkum secara jelas.

7.13.1 Penurunan monotonic lokal

Jika pada setiap iterasi kita memilih arah \(h_m\) yang memiliki korelasi positif dengan negatif gradien, dan learning rate cukup kecil, maka terdapat \(\nu\) kecil sehingga \[ \widehat{R}(F_m) \le \widehat{R}(F_{m-1}). \] Alasannya berasal dari ekspansi Taylor orde satu: \[ \widehat{R}(F_{m-1}+\nu h_m) = \widehat{R}(F_{m-1}) + \nu D\widehat{R}(F_{m-1})[h_m] + o(\nu). \] Jika \(h_m\) dipilih ke arah negatif gradien, maka directional derivative bernilai negatif. Maka untuk \(\nu\) cukup kecil, risiko turun.

7.13.2 Konvergensi ke titik stasioner

Jika: 1. loss convex dan cukup halus, 2. ukuran langkah memenuhi syarat standar, 3. aproksimasi arah gradien cukup baik, maka lintasan boosting dapat dipandang mendekati titik stasioner dari risiko empiris pada closure kelas fungsi aditif yang dihasilkan weak learners.

Namun, karena boosting bekerja pada kelas fungsi terbatas dan sering dengan line search aproksimatif, konvergensi global ke minimizer absolut tidak selalu dijamin dalam pengertian paling kuat. Yang lebih realistis adalah: boosting menuruni risiko secara sistematis dan, pada kondisi regularitas tertentu, mendekati solusi optimal dalam kelas representasi yang tersedia.

7.13.3 Overfitting dan early stopping

Menariknya, walaupun training loss cenderung terus turun, test error sering mencapai minimum pada iterasi terbatas. Ini berarti konvergensi optimisasi tidak sama dengan generalisasi optimal. Dalam konteks statistik, jumlah iterasi bertindak sebagai parameter regularisasi. Jadi, dari perspektif risiko prediktif, kita sering sengaja berhenti sebelum algoritma “selesai” menurunkan training loss.

7.13.4 Newton boosting dan konvergensi lebih cepat

Karena Newton boosting memakai informasi kelengkungan, ia biasanya mencapai penurunan loss lebih cepat daripada first-order boosting, terutama untuk loss yang cukup halus. Secara lokal, metode Newton memiliki laju konvergensi lebih cepat karena menyesuaikan skala pembaruan dengan Hessian. Namun, jika Hessian sangat kecil atau tidak stabil, regularisasi tambahan sering diperlukan.

7.14 Interpretasi Statistik

Bagian ini penting agar boosting tidak berhenti sebagai algoritma komputasi semata.

7.14.1 Boosting sebagai estimator fungsi

Secara statistik, boosting adalah estimator untuk fungsi regresi atau skor klasifikasi: \[ F^*(x) = \arg\min_f \mathbb{E}[L(Y,f(X))]. \] Pada squared error, target populasi adalah mean kondisional: \[ F^*(x)=\mathbb{E}[Y\mid X=x]. \] Pada logistic loss, targetnya berhubungan dengan log-odds kondisional. Jadi boosting dapat dipahami sebagai prosedur estimasi fungsi target tertentu, bukan hanya “prediktor tanpa arti”.

7.14.2 Bias dan variance

Boosting terutama kuat dalam menurunkan bias karena model tumbuh semakin fleksibel. Namun jika pembelajaran diteruskan terlalu jauh atau learner terlalu kompleks, variance meningkat. Inilah alasan regularisasi penting. Jadi, boosting bukan model yang bebas dari bias-variance trade-off; ia hanya mengelola trade-off itu dengan cara yang berbeda dari model parametrik klasik.

7.14.3 Interpretasi variabel

Ukuran seperti variable importance atau partial dependence bersifat prediktif, bukan kausal. Mereka menjelaskan bagaimana model menggunakan informasi variabel untuk meminimalkan loss, bukan bagaimana variabel tersebut menyebabkan perubahan pada respons. Jika desain studi kausal tidak tersedia, maka interpretasi boosting harus tetap berada pada level asosiasi prediktif.

7.14.4 Ketidakpastian

Boosting standar tidak secara otomatis memberikan interval kepercayaan parameter seperti regresi linear. Hal ini bukan kekurangan fatal, tetapi mencerminkan bahwa objek yang diestimasi adalah fungsi kompleks. Ketidakpastian dapat didekati melalui bootstrap, conformal prediction, atau teknik Bayesian/hibrida, tetapi tidak hadir secara langsung dalam bentuk koefisien klasik.

7.15 Formulasi Time Series sebagai Masalah Supervised Learning

Untuk deret waktu univariat, misalkan tersedia observasi \(\{y_t\}_{t=1}^{T}\). Pendekatan supervised learning menulis \[ y_t = f(y_{t-1},y_{t-2},\ldots,y_{t-p},z_t) + \varepsilon_t, \] dengan \(z_t\) dapat mencakup tren, dummy musiman, hari libur, atau kovariat eksternal. Maka setiap waktu \(t\) menjadi satu observasi dengan fitur \[ x_t = (y_{t-1},\ldots,y_{t-p}, z_t). \] Sekarang boosting dapat digunakan untuk mengestimasi fungsi regresi \[ f(x_t)\approx \mathbb{E}[y_t \mid x_t]. \]

Keuntungan pendekatan ini adalah fleksibilitas. Hubungan nonlinier antara lag dan respons dapat dipelajari tanpa spesifikasi parametrik ketat. Interaksi antara fitur musiman dan lag juga dapat ditangkap otomatis oleh pohon.

Namun pendekatan ini hanya valid jika validasi menjaga urutan waktu. Cross-validation acak standar tidak cocok karena ia membocorkan informasi masa depan ke masa lalu.

7.16 Implementasi Regresi Nonlinier di R

library(gbm)

set.seed(123)

n <- 400
x1 <- runif(n, -2, 2)
x2 <- runif(n, -1, 1)
y  <- 2*sin(pi*x1) + x2^2 - 1.5*x1*x2 + rnorm(n, 0, 0.4)

dat <- data.frame(y = y, x1 = x1, x2 = x2)

fit_gbm <- gbm(
  y ~ x1 + x2,
  data = dat,
  distribution = "gaussian",
  n.trees = 3000,
  interaction.depth = 3,
  shrinkage = 0.01,
  n.minobsinnode = 5,
  bag.fraction = 0.7,
  cv.folds = 5,
  verbose = FALSE
)

best_iter <- gbm.perf(fit_gbm, method = "cv", plot.it = FALSE)
best_iter

#> [1] 2905

Evaluasi training:

pred_gbm <- predict(fit_gbm, newdata = dat, n.trees = best_iter)
rmse_gbm <- sqrt(mean((dat$y - pred_gbm)^2))
rmse_gbm

#> [1] 0.2833577

Perbandingan dengan model linear:

fit_lm <- lm(y ~ x1 + x2 + I(x1*x2), data = dat)
pred_lm <- predict(fit_lm, newdata = dat)
rmse_lm <- sqrt(mean((dat$y - pred_lm)^2))
rmse_lm

#> [1] 1.428002

Variable importance:

summary(fit_gbm, n.trees = best_iter, plotit = FALSE)

#>    var  rel.inf
#> x1  x1 66.63984
#> x2  x2 33.36016

7.17 Implementasi Time Series Serius di R

Contoh berikut memakai deret bulanan AirPassengers, tetapi dengan struktur fitur yang lebih serius daripada contoh lag minimal. Kita tambahkan lag jangka pendek, lag musiman, tren, serta faktor bulan.

ap <- as.numeric(AirPassengers)
tt <- time(AirPassengers)
mo <- cycle(AirPassengers)

lag_construct <- data.frame(
  y     = ap[25:length(ap)],
  lag1  = ap[24:(length(ap)-1)],
  lag2  = ap[23:(length(ap)-2)],
  lag3  = ap[22:(length(ap)-3)],
  lag12 = ap[13:(length(ap)-12)],
  lag24 = ap[1:(length(ap)-24)],
  trend = 25:length(ap),
  month = factor(mo[25:length(ap)])
)

head(lag_construct)

#>     y lag1 lag2 lag3 lag12 lag24 trend month
#> 1 145  140  114  133   115   112    25     1
#> 2 150  145  140  114   126   118    26     2
#> 3 178  150  145  140   141   132    27     3
#> 4 163  178  150  145   135   129    28     4
#> 5 172  163  178  150   125   121    29     5
#> 6 178  172  163  178   149   135    30     6

Pisahkan data secara kronologis:

n_all <- nrow(lag_construct)
n_train <- floor(0.8 * n_all)

train_ts <- lag_construct[1:n_train, ]
test_ts  <- lag_construct[(n_train+1):n_all, ]

Pasang boosting:

fit_ts <- gbm(
  y ~ lag1 + lag2 + lag3 + lag12 + lag24 + trend + month,
  data = train_ts,
  distribution = "gaussian",
  n.trees = 2500,
  interaction.depth = 3,
  shrinkage = 0.01,
  n.minobsinnode = 4,
  bag.fraction = 0.8,
  verbose = FALSE
)

pred_ts <- predict(fit_ts, newdata = test_ts, n.trees = 2500)
rmse_ts <- sqrt(mean((test_ts$y - pred_ts)^2))
mae_ts  <- mean(abs(test_ts$y - pred_ts))

c(RMSE = rmse_ts, MAE = mae_ts)

#>     RMSE      MAE 
#> 54.28827 42.92418

7.18 Rolling Origin Evaluation untuk Time Series

Validasi yang benar perlu mempertahankan arah waktu. Salah satu pendekatan adalah rolling origin atau expanding window backtesting.

errors_rmse <- c()

for (i in 80:(nrow(lag_construct)-1)) {
  tr <- lag_construct[1:i, ]
  te <- lag_construct[i+1, , drop = FALSE]

  fit_roll <- gbm(
    y ~ lag1 + lag2 + lag3 + lag12 + lag24 + trend + month,
    data = tr,
    distribution = "gaussian",
    n.trees = 500,
    interaction.depth = 2,
    shrinkage = 0.01,
    n.minobsinnode = 4,
    bag.fraction = 0.8,
    verbose = FALSE
  )

  pr <- predict(fit_roll, newdata = te, n.trees = 500)
  errors_rmse <- c(errors_rmse, te$y - pr)
}

sqrt(mean(errors_rmse^2))

#> [1] 32.64292

Contoh ini lebih serius daripada single split, karena ia mensimulasikan penggunaan model dari waktu ke waktu. Dalam penelitian forecasting, prosedur seperti ini jauh lebih kredibel.

7.19 Catatan Implementasi Newton Boosting dalam Praktik Modern

Paket gbm klasik pada R terutama mengekspresikan boosting dalam semangat Friedman awal. Namun banyak implementasi modern, seperti yang ada di keluarga XGBoost, memanfaatkan ekspansi orde dua. Inti matematikanya tetap formula \[ \gamma_{jm} = -\frac{\sum_{x_i \in R_{jm}} g_i} {\sum_{x_i \in R_{jm}} h_i + \lambda}, \] dengan \(\lambda\) sebagai penalti ridge lokal untuk stabilisasi numerik. Tambahan ini sangat penting ketika jumlah Hessian kecil atau region sempit.

Dalam konteks ini, pembagian split pohon juga sering dinilai berdasarkan gain second-order: \[ \text{Gain} = \frac{1}{2} \left[ \frac{G_L^2}{H_L+\lambda} + \frac{G_R^2}{H_R+\lambda} - \frac{(G_L+G_R)^2}{H_L+H_R+\lambda} \right] -\gamma, \] dengan \(G\) agregat gradien, \(H\) agregat Hessian, dan \(\gamma\) penalti kompleksitas split. Walaupun bab ini fokus pada derivasi boosting dasar, formula ini menunjukkan bagaimana teori Newton boosting diterjemahkan ke implementasi modern yang efisien.

7.20 Ringkasan Akademik

Bagian-bagian sebelumnya dapat diringkas sebagai berikut.

Boosting meminimalkan risiko empiris dalam ruang fungsi.
Ia membangun model secara aditif dan bertahap.
Pada setiap iterasi, ia bergerak ke arah negatif gradien risiko; jika dibatasi oleh weak learner, ia melakukan proyeksi ke kelas fungsi tersebut.
Untuk loss Gaussian, pseudo-residual sama dengan residual biasa.
Untuk loss non-Gaussian seperti logistic dan Huber, pseudo-residual memiliki bentuk berbeda yang langsung berasal dari gradien loss.
Newton boosting menggunakan gradien dan Hessian sekaligus, menghasilkan pembaruan per region \[ \gamma_{jm}= -\frac{\sum g_i}{\sum h_i}. \]
Konvergensi optimisasi tidak identik dengan generalisasi terbaik; early stopping tetap diperlukan.
Dalam time series, boosting sah digunakan setelah data dibentuk sebagai supervised learning dan dievaluasi dengan prosedur berbasis waktu.

7.21 Penutup

Gradient boosting layak dipandang sebagai salah satu kontribusi besar dalam statistik komputasional modern karena ia memberi jembatan yang sangat elegan antara optimisasi dan pembelajaran prediktif. Dari sisi teori, boosting menunjukkan bahwa ide gradient descent dapat diperluas dari ruang parameter berhingga ke ruang fungsi. Dari sisi praktik, boosting menyediakan cara yang sangat efektif untuk membangun model nonlinier dan interaktif tanpa harus menentukan bentuk fungsi global secara kaku sejak awal.

Namun, kekuatan ini hanya bermanfaat bila dipadukan dengan disiplin metodologis. Learning rate, kompleksitas weak learner, subsampling, dan early stopping bukan aksesori tambahan, melainkan inti regularisasi. Demikian pula, pada time series, validasi yang menghormati arah waktu bukan pilihan opsional, tetapi syarat sah. Dan yang paling penting, interpretasi boosting harus selalu dibedakan antara makna prediktif dan klaim kausal.

Dengan fondasi derivasi formal, Newton update berbobot Hessian, serta diskusi konvergensi dan interpretasi statistik yang telah dibahas dalam bab ini, pembaca diharapkan tidak lagi melihat boosting sekadar sebagai “mesin akurasi tinggi”, tetapi sebagai prosedur estimasi fungsi yang serius, kaya secara teori, dan sangat relevan untuk penelitian modern.

7.22 Latihan Teoretis

Buktikan kembali bahwa pada loss Gaussian, pseudo-residual sama dengan residual biasa.
Turunkan Hessian untuk logistic loss dan tunjukkan bahwa nilainya selalu non-negatif.
Dengan menggunakan ekspansi Taylor orde dua, turunkan kembali formula Newton update per region.
Jelaskan perbedaan konseptual antara gradient boosting first-order dan second-order.
Mengapa early stopping dapat dipandang sebagai regularisasi? Jelaskan dalam bahasa bias-variance.
Pada time series bulanan, mengapa lag 12 dan dummy bulan sering sama-sama informatif tetapi tidak identik?

Ilustrasi Manual Sederhana untuk Logistic Loss

Bagian ini memberi ilustrasi manual untuk boosting pada klasifikasi biner dengan logistic loss. Tujuannya adalah memperlihatkan bahwa pada klasifikasi, residual semu tidak lagi sama dengan residual biasa. Sebaliknya, residual semu berasal langsung dari gradien fungsi loss terhadap prediktor.

Misalkan terdapat empat observasi dengan respons biner yang dikodekan sebagai \[ y_i \in \{-1,+1\}. \]

Gunakan data sederhana berikut: \[ \begin{array}{c|c|c} i & x_i & y_i \\ \hline 1 & 1 & -1 \\ 2 & 2 & -1 \\ 3 & 3 & +1 \\ 4 & 4 & +1 \\ \end{array} \]

Kita gunakan logistic loss: \[ L(y,F)=\log\big(1+\exp(-2yF)\big). \]

Inisialisasi

Untuk ilustrasi paling sederhana, misalkan kita mulai dari \[ F_0(x)=0. \]

Karena \(F_0(x)=0\), probabilitas awal kelas positif adalah \[ p_0(x)=\frac{1}{1+\exp(-2F_0(x))}=\frac{1}{2}. \]

Artinya, sebelum belajar apa pun, model memberi probabilitas 0.5 untuk semua observasi.

Hitung pseudo-residual logistic

Pseudo-residual untuk logistic loss adalah \[ r_{im} = \frac{2y_i}{1+\exp(2y_iF_{m-1}(x_i))}. \]

Karena \(F_0(x_i)=0\), maka \[ r_{i1} = \frac{2y_i}{1+\exp(0)} = \frac{2y_i}{2} = y_i. \]

Jadi pada iterasi pertama: \[ \begin{array}{c|c|c} i & y_i & r_{i1} \\ \hline 1 & -1 & -1 \\ 2 & -1 & -1 \\ 3 & +1 & +1 \\ 4 & +1 & +1 \\ \end{array} \]

Menariknya, pada iterasi pertama residual semu sama dengan label kelas, tetapi ini hanya terjadi karena model awal nol. Pada iterasi berikutnya, residual semu akan tergantung pada margin prediksi saat ini.

Fit weak learner sederhana

Pakai stump dengan split: \[ R_{11}=\{x: x \le 2.5\}, \qquad R_{21}=\{x: x > 2.5\}. \]

Pada region kiri, semua label adalah \(-1\), sedangkan pada region kanan semua label adalah \(+1\). Agar ilustrasi tetap sederhana, kita gunakan nilai terminal mentah: \[ h_1(x)= \begin{cases} -1, & x \le 2.5,\\ +1, & x > 2.5. \end{cases} \]

Gunakan learning rate \[ \nu=0.5. \]

Update model

Maka \[ F_1(x)=F_0(x)+\nu h_1(x) = \begin{cases} -0.5, & x \le 2.5,\\ +0.5, & x > 2.5. \end{cases} \]

Probabilitas kelas positif menjadi: - untuk kiri: \[ p_1(x)=\frac{1}{1+\exp(-2(-0.5))}=\frac{1}{1+e^1}\approx 0.269, \] - untuk kanan: \[ p_1(x)=\frac{1}{1+\exp(-2(0.5))}=\frac{1}{1+e^{-1}}\approx 0.731. \]

Jadi setelah satu iterasi, model sudah mulai memisahkan dua kelas.

Pseudo-residual iterasi kedua

Sekarang hitung residual semu lagi.

Region kiri: \(y=-1\), \(F_1=-0.5\)

\[ r_{i2} = \frac{2(-1)}{1+\exp(2(-1)(-0.5))} = \frac{-2}{1+\exp(1)} \approx \frac{-2}{3.718} \approx -0.538. \]

Region kanan: \(y=+1\), \(F_1=0.5\)

\[ r_{i2} = \frac{2(1)}{1+\exp(2(1)(0.5))} = \frac{2}{1+\exp(1)} \approx 0.538. \]

Maka residual semu iterasi kedua adalah: \[ \begin{array}{c|c|c|c} i & y_i & F_1(x_i) & r_{i2} \\ \hline 1 & -1 & -0.5 & -0.538 \\ 2 & -1 & -0.5 & -0.538 \\ 3 & +1 & +0.5 & +0.538 \\ 4 & +1 & +0.5 & +0.538 \\ \end{array} \]

Perhatikan bahwa residual semu sekarang lebih kecil dalam nilai absolut dibanding iterasi pertama. Ini menunjukkan bahwa ketika model mulai yakin dan benar, besar koreksi berikutnya mengecil. Inilah perilaku khas logistic boosting: observasi yang sudah terklasifikasi dengan baik memberi kontribusi gradien yang makin kecil.

Pelajaran dari ilustrasi logistic

Ada beberapa poin penting.

Pertama, logistic boosting tidak memperbarui model berdasarkan residual biasa, melainkan berdasarkan gradien loss. Kedua, besar residual semu bergantung pada margin prediksi saat ini. Ketiga, jika suatu observasi sudah diprediksi dengan benar dan cukup yakin, maka gradiennya mengecil sehingga algoritma lebih fokus pada observasi yang masih sulit.

Dengan demikian, boosting untuk klasifikasi dapat dipahami sebagai prosedur yang berulang kali menyesuaikan skor klasifikasi agar margin makin baik.

Ilustrasi Manual Sederhana untuk Newton Boosting

Bagian ini memberi ilustrasi manual sederhana untuk Newton boosting. Tujuannya adalah memperlihatkan bagaimana pembaruan orde dua berbeda dari gradient boosting biasa. Pada Newton boosting, kita tidak hanya memakai gradien, tetapi juga memakai Hessian untuk mengatur ukuran langkah.

Agar tetap sederhana, kita gunakan logistic loss yang sama: \[ L(y,F)=\log\big(1+\exp(-2yF)\big). \]

Gradien: \[ g_i = -\frac{2y_i}{1+\exp(2y_iF_i)}, \] dan Hessian: \[ h_i = \frac{4\exp(2y_iF_i)}{(1+\exp(2y_iF_i))^2}. \]

Newton update per region: \[ \gamma_j = -\frac{\sum_{i \in R_j} g_i}{\sum_{i \in R_j} h_i}. \]

Kita gunakan data yang sama seperti sebelumnya: \[ \begin{array}{c|c|c} i & x_i & y_i \\ \hline 1 & 1 & -1 \\ 2 & 2 & -1 \\ 3 & 3 & +1 \\ 4 & 4 & +1 \\ \end{array} \]

dan ambil model awal \[ F_0(x)=0. \]

Hitung gradien dan Hessian pada iterasi awal

Karena \(F_0=0\), maka:

Gradien

Untuk \(y=-1\): \[ g_i = -\frac{2(-1)}{1+\exp(0)} = \frac{2}{2} = 1. \]

Untuk \(y=+1\): \[ g_i = -\frac{2(1)}{1+\exp(0)} = -1. \]

Jadi: \[ \begin{array}{c|c|c} i & y_i & g_i \\ \hline 1 & -1 & 1 \\ 2 & -1 & 1 \\ 3 & +1 & -1 \\ 4 & +1 & -1 \\ \end{array} \]

Hessian

Karena \(F_0=0\), maka untuk semua observasi: \[ h_i = \frac{4\exp(0)}{(1+\exp(0))^2} = \frac{4}{(1+1)^2} = \frac{4}{4} = 1. \]

Jadi semua Hessian awal sama dengan 1.

Bentuk region

Gunakan split yang sama: \[ R_1=\{x: x \le 2.5\}, \qquad R_2=\{x: x > 2.5\}. \]

Region kiri

Di kiri, gradiennya: \[ g_1=1,\quad g_2=1. \]

Hessiannya: \[ h_1=1,\quad h_2=1. \]

Maka Newton update region kiri: \[ \gamma_1 = -\frac{1+1}{1+1} = -\frac{2}{2} = -1. \]

Region kanan

Di kanan, gradiennya: \[ g_3=-1,\quad g_4=-1. \]

Hessiannya: \[ h_3=1,\quad h_4=1. \]

Maka Newton update region kanan: \[ \gamma_2 = -\frac{-1+(-1)}{1+1} = -\frac{-2}{2} = 1. \]

Jadi weak learner Newton menghasilkan update: \[ \eta_1(x)= \begin{cases} -1, & x \le 2.5,\\ +1, & x > 2.5. \end{cases} \]

Update model Newton

Jika learning rate diambil \[ \nu=1, \] maka: \[ F_1(x)=F_0(x)+\eta_1(x) = \begin{cases} -1, & x \le 2.5,\\ +1, & x > 2.5. \end{cases} \]

Bandingkan dengan ilustrasi gradient boosting biasa tadi yang, dengan learning rate 0.5, hanya bergerak ke \(\pm 0.5\). Newton boosting bergerak lebih langsung karena ia memanfaatkan informasi kelengkungan.

Probabilitas setelah Newton update

Untuk kiri: \[ p_1(x)=\frac{1}{1+\exp(-2(-1))}=\frac{1}{1+e^2}\approx 0.119. \]

Untuk kanan: \[ p_1(x)=\frac{1}{1+\exp(-2(1))}=\frac{1}{1+e^{-2}}\approx 0.881. \]

Jadi setelah satu langkah Newton, model menjadi jauh lebih yakin dibanding langkah gradient boosting pertama.

Mengapa Hessian penting?

Formula \[ \gamma_j = -\frac{\sum g_i}{\sum h_i} \] menunjukkan bahwa update tidak hanya melihat arah gradien, tetapi juga skala kelengkungan lokal.

Jika Hessian besar, maka loss sangat melengkung dan update diperkecil.
Jika Hessian kecil, maka update bisa lebih besar.
Jadi Newton boosting menyesuaikan langkah secara adaptif terhadap geometri lokal loss.

Dalam contoh awal ini semua Hessian sama dengan 1, sehingga update Newton tampak mirip dengan rata-rata gradien bertanda negatif. Tetapi pada iterasi-iterasi berikutnya, ketika \(F(x)\) tidak lagi nol, Hessian akan berbeda antarobservasi, dan pembobotan menjadi sangat penting.

Pelajaran dari ilustrasi Newton boosting

Ilustrasi ini memberi beberapa pelajaran inti.

Pertama, Newton boosting adalah generalisasi second-order dari gradient boosting. Kedua, pembaruan per region diperoleh dari agregat gradien dibagi agregat Hessian, sehingga bentuknya seperti rata-rata tertimbang. Ketiga, Newton boosting biasanya bergerak lebih cepat menuju solusi karena ia menggunakan informasi kelengkungan. Keempat, pada implementasi modern seperti XGBoost, bentuk update inilah yang menjadi fondasi utama.

Dengan demikian, jika gradient boosting biasa dapat dilihat sebagai “mengikuti kemiringan lereng”, maka Newton boosting lebih seperti “mengikuti kemiringan sambil juga membaca tingkat kelengkungan tanah”, sehingga langkahnya sering lebih efisien.

7.23 Bukti Alternatif melalui Directional Derivative

Bagian ini memberi pembuktian yang sedikit lebih formal mengenai hubungan boosting dan gradient descent melalui turunan arah (directional derivative). Misalkan \(F\) adalah fungsi saat ini dan \(h\) adalah fungsi kandidat arah pembaruan. Definisikan kurva \[ \phi(\alpha)=\widehat{R}(F+\alpha h) = \sum_{i=1}^{n}L\big(y_i,F(x_i)+\alpha h(x_i)\big). \] Turunan pertama terhadap \(\alpha\) di titik nol adalah \[ \phi'(0) = \sum_{i=1}^{n} \frac{\partial L(y_i,F(x_i))}{\partial F(x_i)} h(x_i). \] Jika kita definisikan gradien empiris \[ g_i(F)=\frac{\partial L(y_i,F(x_i))}{\partial F(x_i)}, \] maka \[ \phi'(0)=\sum_{i=1}^{n} g_i(F)h(x_i) = \langle \nabla_F \widehat{R}(F), h \rangle_n. \] Oleh karena itu, arah yang meminimalkan turunan arah di bawah kendala norm empiris satu adalah \[ h^* = - \frac{\nabla_F \widehat{R}(F)}{\|\nabla_F \widehat{R}(F)\|_n}. \] Jadi, negatif gradien adalah arah steepest descent dalam geometri empiris biasa.

Sekarang, boosting tidak mencari \(h^*\) di seluruh ruang fungsi, tetapi hanya pada kelas \(\mathcal{H}\). Maka kita ingin memilih \[ h_m = \arg\min_{h \in \mathcal{H}, \|h\|_n=1} \langle \nabla_F \widehat{R}(F_{m-1}), h \rangle_n. \] Karena meminimalkan inner product terhadap gradien setara dengan mendekatkan \(h\) ke negatif gradien, kita memperoleh kembali prinsip aproksimasi residual semu. Perspektif directional derivative ini menunjukkan bahwa boosting bukan hanya analog informal dari gradient descent, melainkan realisasi tepat dari konsep steepest descent yang diproyeksikan.

Jika line search dilakukan secara eksplisit, maka setelah arah \(h_m\) dipilih kita menyelesaikan \[ \rho_m = \arg\min_{\rho} \widehat{R}(F_{m-1}+\rho h_m). \] Bila line search tidak dilakukan penuh, maka learning rate tetap \(\nu\) menghasilkan pembaruan \[ F_m = F_{m-1} + \nu h_m. \] Dalam kedua kasus, justifikasi optimisasinya sama: kita memilih arah yang menurunkan turunan arah, lalu mengambil langkah di sepanjang arah itu.

7.24 Newton Update Per Region dalam Bentuk Berbobot

Walaupun formula \[ \gamma_{jm} = -\frac{\sum g_i}{\sum h_i} \] sudah sangat penting, masih berguna untuk menuliskannya dalam bentuk yang lebih “statistik”, yaitu sebagai rata-rata tertimbang. Misalkan pada region \(R_{jm}\) didefinisikan pseudo-response Newton \[ z_i = -\frac{g_i}{h_i}, \] dengan bobot \[ w_i = h_i. \] Maka jumlah kuadratik lokal orde dua dapat ditulis, hingga konstanta yang tidak bergantung pada \(\gamma_{jm}\), sebagai \[ Q(\gamma_{jm}) = \sum_{x_i \in R_{jm}} \left[ g_i \gamma_{jm} + \frac{1}{2} h_i \gamma_{jm}^2 \right]. \] Lengkapi kuadrat: \[ g_i \gamma + \frac{1}{2}h_i \gamma^2 = \frac{1}{2}h_i \left( \gamma^2 + 2\frac{g_i}{h_i}\gamma \right) = \frac{1}{2}h_i \left( \gamma - z_i \right)^2 - \frac{1}{2}h_i z_i^2. \] Maka meminimalkan \(Q(\gamma_{jm})\) ekuivalen dengan meminimalkan \[ \sum_{x_i \in R_{jm}} w_i (\gamma_{jm}-z_i)^2. \] Solusi minimisasi kuadrat berbobot adalah \[ \gamma_{jm} = \frac{\sum_{x_i \in R_{jm}} w_i z_i} {\sum_{x_i \in R_{jm}} w_i}. \] Substitusi \(w_i=h_i\) dan \(z_i=-g_i/h_i\) memberi \[ \gamma_{jm} = \frac{\sum h_i(-g_i/h_i)}{\sum h_i} = -\frac{\sum g_i}{\sum h_i}. \]

Representasi ini sangat penting karena menunjukkan bahwa Newton boosting dapat dipandang sebagai weighted least squares fitting lokal pada pseudo-response \(z_i\) dengan bobot Hessian \(w_i\). Secara konseptual, ini menjelaskan hubungan yang sangat dekat antara Newton boosting, IRLS, dan beberapa algoritma optimisasi GLM. Dari sisi implementasi, inilah sebabnya banyak perangkat lunak modern menghitung agregat gradien dan Hessian per node terminal lalu memperbarui nilai node melalui rasio tersebut.

7.25 Analisis Konvergensi yang Lebih Rinci

Konvergensi boosting dapat dibedakan menjadi tiga level: penurunan risiko empiris, konvergensi algoritmik ke titik stasioner, dan performa generalisasi.

7.25.1 Penurunan risiko empiris

Misalkan loss terdiferensialkan dan Lipschitz-smooth dalam prediktor, yaitu terdapat konstanta \(C>0\) sehingga \[ \widehat{R}(F+\Delta) \le \widehat{R}(F) + D\widehat{R}(F)[\Delta] + \frac{C}{2}\|\Delta\|_n^2. \] Jika pada iterasi ke-\(m\) kita ambil \(\Delta_m = \nu h_m\), maka \[ \widehat{R}(F_m) \le \widehat{R}(F_{m-1}) + \nu D\widehat{R}(F_{m-1})[h_m] + \frac{C\nu^2}{2}\|h_m\|_n^2. \] Jika \(h_m\) berkorelasi negatif dengan gradien, maka suku linear bernilai negatif. Untuk \(\nu\) cukup kecil, suku kuadrat tidak mendominasi, sehingga risiko turun. Ini memberi justifikasi formal untuk fakta bahwa learning rate kecil cenderung menstabilkan boosting.

7.25.2 Kondisi stasioner

Anggap kelas weak learner cukup kaya sehingga jika gradien tidak nol, maka ada \(h \in \mathcal{H}\) yang memiliki inner product negatif dengan gradien. Jika algoritma berhenti pada fungsi \(\widetilde{F}\) sehingga tidak ada lagi arah dalam \(\mathcal{H}\) yang menurunkan risiko, maka \[ \langle \nabla_F \widehat{R}(\widetilde{F}), h \rangle_n \ge 0 \qquad \forall h \in \mathcal{H}. \] Ini adalah bentuk kondisi stasioner relatif terhadap kelas learner. Jika \(\mathcal{H}\) cukup kaya untuk mendekati seluruh arah yang relevan, maka kondisi ini mendekati stasioneritas penuh.

7.25.3 Mengapa overfitting tetap mungkin

Konvergensi optimisasi berarti training loss terus menurun. Akan tetapi, target statistik yang kita inginkan adalah risiko populasi atau error prediksi di luar sampel. Training loss yang terus menurun dapat disertai peningkatan variance. Itulah sebabnya jumlah iterasi optimal untuk generalisasi sering jauh lebih kecil daripada jumlah iterasi yang dibutuhkan untuk menurunkan training loss hingga sangat kecil. Dengan kata lain, optimisasi dan regularisasi adalah dua tujuan yang berbeda. Early stopping menjembatani keduanya.

7.25.4 Newton boosting dan kestabilan numerik

Pada Newton boosting, jika Hessian sangat kecil di suatu region, maka langkah \[ -\frac{\sum g_i}{\sum h_i} \] dapat menjadi sangat besar. Karena itu implementasi praktis sering menambahkan penalti ridge lokal: \[ \gamma_{jm} = -\frac{\sum g_i}{\sum h_i + \lambda}. \] Penalti \(\lambda>0\) bukan hanya trik komputasi, tetapi bentuk regularisasi yang mencegah pembaruan liar ketika kelengkungan lokal lemah. Ini juga membantu ketika region berisi sangat sedikit observasi.

7.26 Interpretasi Statistik yang Lebih Kaya

7.26.1 Apa yang sesungguhnya diestimasi?

Untuk setiap loss, boosting menargetkan fungsi populasi yang berbeda. Pada squared error, boosting menaksir mean kondisional. Pada absolute error, ia menaksir median kondisional. Pada quantile loss, ia menaksir kuantil kondisional tertentu. Pada logistic loss, ia menaksir skor yang terkait dengan probabilitas klasifikasi. Jadi, fungsi loss menentukan parameter fungsional populasi yang menjadi target. Ini memberi makna statistik yang kuat: boosting bukan sekadar algoritma komputasi, tetapi prosedur estimasi untuk functional target tertentu.

7.26.2 Hubungan dengan likelihood

Pada banyak loss yang berasal dari negative log-likelihood, boosting dapat dipahami sebagai minimisasi deviance empiris. Misalnya, logistic loss terkait dengan likelihood Bernoulli. Dalam pengertian ini, boosting mendekati semangat estimasi likelihood, tetapi dengan kelas fungsi jauh lebih fleksibel. Jadi, boosting tidak memutus hubungan dengan statistik klasik; ia memperluasnya.

7.26.3 Mengapa interpretasi koefisien hilang?

Dalam regresi linear, interpretasi koefisien datang dari bentuk model global yang eksplisit. Dalam boosting, fungsi akhir adalah penjumlahan ribuan aturan lokal. Oleh sebab itu, “efek” suatu variabel tidak lagi ringkas dalam satu angka. Kita perlu ringkasan seperti importance, partial dependence, atau ALE plots. Kehilangan interpretasi koefisien ini bukan kecelakaan; itu konsekuensi langsung dari fleksibilitas representasi fungsi.

7.26.4 Ketidakpastian prediksi versus ketidakpastian parameter

Pada model parametrik, perhatian sering tertuju pada interval kepercayaan parameter. Dalam boosting, objek alami justru prediksi atau fungsi yang dipelajari. Karena itu, pendekatan ketidakpastian yang lebih relevan sering berupa interval prediksi, conformal prediction, atau bootstrap fungsional, bukan interval untuk “koefisien” yang memang tidak eksis secara global.

7.27 Time Series: Formulasi Multi-Step dan Masalah Validasi

Dalam forecasting, ada dua strategi utama untuk horizon lebih dari satu langkah.

7.27.1 Strategi recursive

Kita membangun model satu langkah \[ y_{t+1} = f(x_t) + \varepsilon_{t+1}, \] lalu menggunakan prediksi langkah pertama sebagai input untuk meramal langkah kedua, dan seterusnya. Strategi ini hemat model, tetapi error dapat terakumulasi.

7.27.2 Strategi direct

Untuk setiap horizon \(h\), kita bangun model berbeda: \[ y_{t+h} = f_h(x_t) + \varepsilon_{t+h}. \] Pendekatan ini menghindari akumulasi error recursive, tetapi membutuhkan lebih banyak model.

Dalam boosting, kedua strategi bisa dipakai. Yang penting, desain validasi harus meniru cara model akan dipakai. Jika targetnya adalah horizon tiga bulan ke depan, maka evaluasi harus menghitung error tiga langkah, bukan hanya satu langkah.

7.27.3 Konstruksi fitur waktu yang lebih serius

Selain lag mentah, fitur time series yang sering penting meliputi: - tren linear dan nonlinear, - dummy bulan atau kuartal, - Fourier terms untuk musiman halus, - indikator hari libur, - kovariat eksogen, - rolling statistics seperti moving average, - interaksi antara musim dan tren.

Secara statistik, ini berarti boosting pada time series sering lebih tepat dipandang sebagai model dynamic regression with nonlinear function approximation daripada sekadar model autoregresif.

7.28 Implementasi Direct Multi-Step Forecast di R

# contoh direct 3-step ahead forecasting dari AirPassengers

ap <- as.numeric(AirPassengers)
mo <- cycle(AirPassengers)

make_direct_data <- function(series, h = 3, max_lag = 12) {
  idx <- (max_lag + 1):(length(series) - h)
  out <- data.frame(
    y_h   = series[idx + h],
    lag1  = series[idx - 1],
    lag2  = series[idx - 2],
    lag3  = series[idx - 3],
    lag12 = series[idx - 12],
    trend = idx,
    month = factor(mo[idx])
  )
  out
}

direct_df <- make_direct_data(ap, h = 3, max_lag = 12)

n_dir <- nrow(direct_df)
train_dir <- direct_df[1:floor(0.8 * n_dir), ]
test_dir  <- direct_df[(floor(0.8 * n_dir)+1):n_dir, ]

fit_dir <- gbm(
  y_h ~ lag1 + lag2 + lag3 + lag12 + trend + month,
  data = train_dir,
  distribution = "gaussian",
  n.trees = 1500,
  interaction.depth = 3,
  shrinkage = 0.01,
  bag.fraction = 0.8,
  verbose = FALSE
)

pred_dir <- predict(fit_dir, newdata = test_dir, n.trees = 1500)
sqrt(mean((test_dir$y_h - pred_dir)^2))

#> [1] 75.81008

Contoh ini menunjukkan bahwa boosting untuk time series tidak harus berhenti pada peramalan satu langkah sederhana. Dengan membangun data untuk horizon tertentu, kita bisa merancang evaluasi yang lebih relevan terhadap kebutuhan substantif.

7.29 Penutup Tambahan

Dengan tambahan bukti melalui turunan arah, representasi Newton berbobot, serta pembahasan konvergensi dan multi-step forecasting, terlihat bahwa boosting bukan metode yang berdiri di luar tradisi statistik. Sebaliknya, ia merupakan perluasan yang sangat alami dari optimisasi dan estimasi fungsi. Kekuatan utamanya terletak pada kemampuan menyatukan fleksibilitas representasi, objektif loss yang jelas, dan regularisasi bertahap. Itulah sebabnya boosting tetap relevan, baik sebagai materi buku ajar tingkat lanjut maupun sebagai fondasi metodologis untuk penelitian modern.

7.30 Referensi Singkat

Friedman, J. H. Greedy Function Approximation: A Gradient Boosting Machine.
Hastie, Tibshirani, dan Friedman. The Elements of Statistical Learning.
Bühlmann dan Hothorn. Boosting Algorithms: Regularization, Prediction and Model Fitting.

Latihan Terstruktur

Level 1 (Konsep Dasar) 1. Jelaskan mengapa pseudo-residual berbeda antara Gaussian dan logistic loss. 2. Apa peran learning rate dalam boosting?

Level 2 (Matematis) 3. Turunkan ulang formula \(\gamma_{jm}\) untuk Newton boosting. 4. Buktikan bahwa gradient boosting adalah projected gradient descent.

Level 3 (Implementasi) 5. Implementasikan boosting sederhana tanpa package (loop manual di R). 6. Bandingkan performa Gaussian vs Huber loss pada data dengan outlier.

Level 4 (Riset) 7. Diskusikan hubungan boosting dengan penalized regression. 8. Bagaimana mengintegrasikan boosting dengan model spatiotemporal (misal INLA)?

Tugas

Turunkan pseudo-residual untuk squared error loss.
Jelaskan peran learning rate dalam boosting dari sudut regularisasi.
Ubah satu data deret waktu menjadi supervised learning dengan beberapa lag.
Bandingkan boosting dengan model linear pada satu dataset regresi.

Ringkasan

Gradient boosting adalah implementasi functional gradient descent untuk membangun model aditif yang sangat fleksibel. Kekuatan utamanya datang dari kombinasi stagewise fitting, weak learners, dan regularisasi melalui shrinkage.

8 UTS: Project Kelompok 2 Orang

8.1 Tujuan Akademik

UTS pada level S2 harus dilihat sebagai mini-research project. Mahasiswa tidak cukup menunjukkan bahwa mereka bisa menjalankan fungsi paket R, tetapi harus menunjukkan kemampuan:

merumuskan pertanyaan statistik,
menentukan target prediksi atau objek unsupervised dengan jelas,
memilih model dan validasi yang sejalan dengan tujuan,
menulis argumen metodologis yang defensibel.

8.2 Struktur Proposal

Proposal UTS minimal memuat:

latar belakang substantif,
pertanyaan statistik,
sumber data,
definisi variabel target,
kandidat model,
metrik evaluasi,
risiko metodologis seperti imbalance, missing value, leakage.

8.3 Standar Reproduksibilitas

Setiap kelompok sebaiknya:

memakai seed yang dicatat,
memisahkan script cleaning, modeling, dan visualisasi,
menyimpan hasil penting dalam format yang bisa ditinjau ulang,
menulis keputusan praproses secara eksplisit.

8.4 Rubrik Penilaian yang Disarankan

kejelasan masalah statistik,
kualitas EDA dan pembersihan data,
ketepatan formulasi model,
kualitas evaluasi,
ketepatan interpretasi,
kerapian komputasi.

8.5 Ringkasan

UTS adalah latihan membangun workflow machine learning yang ilmiah, reprodusibel, dan argumentatif.

9 Supervised Learning I: Naive Bayes dan K-NN

9.1 Tujuan Akademik Bab

Bab ini mempertemukan dua paradigma yang sangat berbeda: Naive Bayes sebagai model generatif dan K-NN sebagai local non-parametric predictor. Pada tingkat S2, pembelajaran di sini harus menekankan:

objek probabilistik yang diestimasi,
asumsi model,
implikasi dimensionalitas,
trade-off bias-variance.

9.2 Naive Bayes sebagai Model Generatif

Dengan Bayes theorem,

\[ P(Y=c \mid X=x) = \frac{P(Y=c)f(x \mid Y=c)} {\sum_{k=1}^{K}P(Y=k)f(x \mid Y=k)}. \]

Asumsi naive menyatakan

\[ f(x \mid Y=c) = \prod_{j=1}^{p}f(x_j \mid Y=c). \]

Jika fitur kontinu Gaussian:

\[ X_j \mid Y=c \sim \mathcal{N}(\mu_{jc}, \sigma_{jc}^2). \]

Maka parameter yang diestimasi adalah:

prior kelas \(\pi_c = P(Y=c)\),
mean per fitur per kelas \(\mu_{jc}\),
varians per fitur per kelas \(\sigma_{jc}^2\).

Estimator alami:

\[ \widehat{\pi}_c = \frac{n_c}{n}, \qquad \widehat{\mu}_{jc} = \frac{1}{n_c}\sum_{i:y_i=c}x_{ij}, \qquad \widehat{\sigma}_{jc}^2 = \frac{1}{n_c-1}\sum_{i:y_i=c}(x_{ij}-\widehat{\mu}_{jc})^2. \]

9.3 Mengapa Asumsi Naive Bisa Tetap Efektif

Meskipun independensi bersyarat sering tidak realistis, Naive Bayes tetap sering bekerja baik karena:

jumlah parameter relatif sedikit,
variance estimator rendah,
tujuan utama sering klasifikasi, bukan model probabilistik sempurna.

Ini contoh klasik bagaimana model dengan misspecification moderat bisa tetap unggul prediktif.

9.4 K-Nearest Neighbors

K-NN mengestimasi fungsi secara lokal. Untuk klasifikasi:

\[ \widehat{g}(x) = \arg\max_c \sum_{i \in N_k(x)} I(y_i=c). \]

Untuk regresi:

\[ \widehat{f}(x) = \frac{1}{k}\sum_{i \in N_k(x)} y_i. \]

K-NN tidak memiliki parameter global. Maka kompleksitas model tidak terletak pada banyaknya parameter, tetapi pada pilihan \(k\), metrik jarak, dan dimensi ruang fitur.

9.5 Curse of Dimensionality

Pada dimensi tinggi, konsep “tetangga dekat” melemah karena banyak titik menjadi hampir sama jauhnya. Secara statistik, inilah sebab utama K-NN dapat menurun drastis pada data high-dimensional tanpa reduksi dimensi atau feature selection.

9.6 Konsistensi dan Bias-Variance

Secara asimtotik, K-NN dapat konsisten jika:

\[ k \to \infty \quad \text{dan} \quad \frac{k}{n} \to 0. \]

Secara intuitif:

\(k\) harus cukup besar agar variance turun,
tetapi tidak boleh terlalu besar agar neighborhood tetap lokal.

9.7 Contoh R: Naive Bayes

library(e1071)
set.seed(123)
idx <- sample(seq_len(nrow(iris)), size = 0.7 * nrow(iris))
train <- iris[idx, ]
test  <- iris[-idx, ]

nb_fit <- naiveBayes(Species ~ ., data = train)
nb_pred <- predict(nb_fit, newdata = test)
table(nb_pred, test$Species)

#>             
#> nb_pred      setosa versicolor virginica
#>   setosa         14          0         0
#>   versicolor      0         18         0
#>   virginica       0          0        13

mean(nb_pred == test$Species)

#> [1] 1

9.8 Contoh R: K-NN

library(class)
train_x <- scale(train[, -5])
test_x  <- scale(test[, -5],
                 center = attr(train_x, "scaled:center"),
                 scale = attr(train_x, "scaled:scale"))
knn_pred <- knn(train_x, test_x, train$Species, k = 5)
table(knn_pred, test$Species)

#>             
#> knn_pred     setosa versicolor virginica
#>   setosa         14          0         0
#>   versicolor      0         17         0
#>   virginica       0          1        13

mean(knn_pred == test$Species)

#> [1] 0.9777778

9.9 Contoh R: Membandingkan Berbagai Nilai k

k_grid <- c(1, 3, 5, 7, 9, 11)
acc_knn <- numeric(length(k_grid))

for (i in seq_along(k_grid)) {
  pred_i <- knn(train_x, test_x, train$Species, k = k_grid[i])
  acc_knn[i] <- mean(pred_i == test$Species)
}

data.frame(k = k_grid, accuracy = acc_knn)

#>    k  accuracy
#> 1  1 0.9333333
#> 2  3 0.9555556
#> 3  5 0.9777778
#> 4  7 0.9777778
#> 5  9 0.9777778
#> 6 11 0.9777778

9.10 Tugas

Turunkan classifier Naive Bayes dari Bayes theorem.
Jelaskan mengapa K-NN sensitif terhadap skala fitur.
Jelaskan kondisi asimtotik sederhana agar K-NN konsisten.
Bandingkan Naive Bayes dan K-NN pada dataset lain dan tulis trade-off-nya.

9.11 Ringkasan

Naive Bayes dan K-NN mewakili dua pendekatan ekstrem: model probabilistik global yang sederhana versus estimator lokal non-parametrik. Keduanya sangat penting untuk memahami bias-variance dan dampak asumsi model.

10 Supervised Learning II: Decision Tree

10.1 Tujuan Akademik Bab

Decision tree sangat penting karena:

mudah diinterpretasikan,
mampu menangkap non-linearitas dan interaksi,
menjadi fondasi bagi Random Forest dan boosting.

Di level S2, tree harus dibaca sebagai prosedur recursive partitioning dengan objective function dan regularisasi tersendiri.

10.2 Model Piecewise Constant

Untuk regresi, decision tree membentuk partisi

\[ \{R_1,\ldots,R_M\} \]

dan prediksi

\[ \widehat{f}(x) = \sum_{m=1}^{M} c_m I(x \in R_m), \]

dengan

\[ c_m = \arg\min_c \sum_{i:x_i \in R_m}(y_i-c)^2 = \bar{y}_{R_m}. \]

Untuk klasifikasi, leaf memuat kelas mayoritas atau probabilitas empiris kelas.

10.3 Kriteria Split

Pada klasifikasi, ukuran impurity:

\[ H(S) = -\sum_{c=1}^{C}p_c \log p_c, \]

\[ G(S) = 1 - \sum_{c=1}^{C}p_c^2. \]

Pada regresi, split dipilih untuk meminimalkan penjumlahan RSS antar-child nodes:

\[ \sum_{i \in R_1}(y_i-\bar{y}_{R_1})^2 + \sum_{i \in R_2}(y_i-\bar{y}_{R_2})^2. \]

10.4 CART dan Greedy Search

Decision tree dibangun secara greedy:

cari split terbaik pada root,
pecah data,
ulangi pada tiap child node.

Ia tidak menjamin global optimum, tetapi efektif secara komputasi. Kelemahannya adalah instabilitas: split awal sangat menentukan seluruh pohon.

10.5 Cost-Complexity Pruning

Regularisasi tree dilakukan dengan pruning:

\[ R_\alpha(T) = R(T) + \alpha |T|, \]

dengan \(|T|\) jumlah terminal nodes. Ini analog dengan penalized estimation pada model parametrik: kompleksitas pohon dihukum agar generalisasi membaik.

10.6 Penurunan Impurity sebagai Fungsi Objektif Lokal

Misalkan sebuah node \(t\) dipecah menjadi child kiri \(t_L\) dan child kanan \(t_R\). Jika \(I(t)\) adalah impurity pada node \(t\), maka kualitas split dapat ditulis sebagai

\[ \Delta I = I(t) - \frac{N_L}{N_t}I(t_L) - \frac{N_R}{N_t}I(t_R). \]

Tree memilih split yang memaksimalkan \(\Delta I\). Karena pemilihan dilakukan node demi node tanpa melihat konsekuensi global ke seluruh pohon, CART bersifat greedy. Di sisi lain, formulasi ini menunjukkan dengan jelas bahwa tree tetap berbasis objective function yang terukur.

Untuk regresi, jika impurity didefinisikan sebagai mean squared error di dalam node, maka leaf predictor optimal diperoleh dari

\[ c_t = \arg\min_c \sum_{i:x_i \in t}(y_i-c)^2 = \bar{y}_t. \]

Untuk klasifikasi, jika loss-nya 0-1 pada leaf, estimator optimal adalah kelas mayoritas. Jadi tree dapat dilihat sebagai deret masalah optimisasi lokal yang sangat sederhana tetapi disusun secara rekursif.

10.7 Surrogate Splits dan Missing Values

Beberapa implementasi tree dapat menggunakan surrogate splits, yaitu split pengganti saat fitur utama hilang. Ini menunjukkan fleksibilitas tree dalam menangani missingness relatif lebih baik dibanding banyak metode berbasis jarak.

10.8 Bias-Variance pada Tree

Single tree biasanya:

bias moderat,
variance tinggi.

Karena itu ia sangat diuntungkan oleh bagging. Secara konseptual, memahami kelemahan single tree adalah pintu masuk natural ke Random Forest.

10.9 Contoh R: Tree Klasifikasi

library(rpart)
library(rpart.plot)
tree_cls <- rpart(Species ~ ., data = iris, method = "class")
rpart.plot(tree_cls, type = 2, extra = 104)

printcp(tree_cls)

#> 
#> Classification tree:
#> rpart(formula = Species ~ ., data = iris, method = "class")
#> 
#> Variables actually used in tree construction:
#> [1] Petal.Length Petal.Width 
#> 
#> Root node error: 100/150 = 0.66667
#> 
#> n= 150 
#> 
#>     CP nsplit rel error xerror     xstd
#> 1 0.50      0      1.00   1.17 0.050735
#> 2 0.44      1      0.50   0.62 0.060310
#> 3 0.01      2      0.06   0.11 0.031927

10.10 Contoh R: Tree Regresi

tree_reg <- rpart(mpg ~ wt + hp + disp + cyl, data = mtcars, method = "anova")
rpart.plot(tree_reg, type = 2, extra = 101)

10.11 Contoh R: Pruning

best_cp <- tree_cls$cptable[which.min(tree_cls$cptable[, "xerror"]), "CP"]
tree_pruned <- prune(tree_cls, cp = best_cp)
rpart.plot(tree_pruned, type = 2, extra = 104)

10.12 Contoh R: Akurasi Full Tree versus Pruned Tree pada Hold-Out

set.seed(321)
idx_tree <- sample(seq_len(nrow(iris)), size = 0.7 * nrow(iris))
train_tree <- iris[idx_tree, ]
test_tree  <- iris[-idx_tree, ]

tree_full <- rpart(Species ~ ., data = train_tree, method = "class", cp = 0)
cp_star <- tree_full$cptable[which.min(tree_full$cptable[, "xerror"]), "CP"]
tree_small <- prune(tree_full, cp = cp_star)

pred_full <- predict(tree_full, newdata = test_tree, type = "class")
pred_small <- predict(tree_small, newdata = test_tree, type = "class")

table(pred_small, test_tree$Species)

#>             
#> pred_small   setosa versicolor virginica
#>   setosa         17          0         0
#>   versicolor      0         13         1
#>   virginica       0          0        14

c(
  accuracy_full_tree = mean(pred_full == test_tree$Species),
  accuracy_pruned_tree = mean(pred_small == test_tree$Species),
  terminal_nodes_full = sum(tree_full$frame$var == "<leaf>"),
  terminal_nodes_pruned = sum(tree_small$frame$var == "<leaf>")
)

#>    accuracy_full_tree  accuracy_pruned_tree   terminal_nodes_full 
#>             0.9777778             0.9777778             3.0000000 
#> terminal_nodes_pruned 
#>             3.0000000

10.13 Tugas

Bandingkan entropy dan Gini secara konsep.
Jelaskan mengapa recursive partitioning bersifat greedy.
Tunjukkan peran pruning sebagai regularisasi.
Bangun tree klasifikasi dan tree regresi pada dua dataset berbeda.

10.14 Ringkasan

Decision tree adalah model partisi lokal yang sangat interpretatif, tetapi secara statistik tidak stabil. Pruning dan ensemble methods muncul sebagai respon terhadap kelemahan ini.

11 Supervised Learning III: Random Forest

11.1 Intuisi Statistik

Random Forest dibangun dari ide bahwa estimator tidak stabil dapat diperbaiki dengan averaging. Jika satu tree memiliki variance tinggi, maka rata-rata banyak tree yang saling tidak terlalu berkorelasi akan menghasilkan prediktor dengan variance lebih kecil.

11.2 Bagging

Jika \(\widehat{f}^{(b)}\) adalah tree dari bootstrap sample ke-\(b\), maka untuk regresi:

\[ \widehat{f}^{RF}(x) = \frac{1}{B}\sum_{b=1}^{B}\widehat{f}^{(b)}(x). \]

Untuk klasifikasi, prediksi diperoleh melalui mayoritas suara.

11.3 Randomization pada Split

Random Forest tidak hanya melakukan bagging. Pada tiap split, hanya subset acak prediktor berukuran mtry yang dipertimbangkan. Ini mengurangi korelasi antar tree, yang secara langsung meningkatkan manfaat averaging.

11.4 Out-of-Bag Error

Sekitar sepertiga observasi tidak masuk bootstrap sample untuk satu tree. Observasi ini menjadi pseudo test set untuk tree tersebut. OOB error adalah estimasi generalization error internal yang sangat berguna dan efisien.

11.5 Varians Ensemble dan Korelasi antar Tree

Jika setiap tree memiliki varians \(\sigma^2\) dan korelasi pasangan rata-rata \(\rho\), maka secara heuristik varians rata-rata \(B\) tree adalah

\[ \text{Var}\left(\frac{1}{B}\sum_{b=1}^{B}T_b\right) = \rho \sigma^2 + \frac{1-\rho}{B}\sigma^2. \]

Formula ini sangat penting secara pedagogis. Ia menunjukkan bahwa menambah jumlah tree \(B\) memang menurunkan komponen \((1-\rho)\sigma^2/B\), tetapi tidak dapat menghilangkan komponen \(\rho \sigma^2\). Inilah alasan random feature selection penting: ia menurunkan korelasi antar tree, bukan hanya menambah banyaknya tree.

11.6 Variable Importance

Dua ukuran umum:

mean decrease in impurity,
permutation importance.

Permutation importance secara konseptual lebih dekat pada kontribusi prediktif karena mengukur penurunan performa saat satu prediktor diacak.

11.7 Random Forest dan Konsistensi

Secara teori, konsistensi Random Forest cukup rumit dan bergantung pada versi algoritma. Tetapi secara praktis, model ini sangat kuat untuk data tabular karena:

menangani non-linearitas,
menangani interaksi,
relatif tahan terhadap scaling,
memiliki OOB error internal.

11.8 Contoh R

library(randomForest)
rf_fit <- randomForest(
  Species ~ .,
  data = iris,
  ntree = 500,
  mtry = 2,
  importance = TRUE
)
print(rf_fit)

#> 
#> Call:
#>  randomForest(formula = Species ~ ., data = iris, ntree = 500,      mtry = 2, importance = TRUE) 
#>                Type of random forest: classification
#>                      Number of trees: 500
#> No. of variables tried at each split: 2
#> 
#>         OOB estimate of  error rate: 4.67%
#> Confusion matrix:
#>            setosa versicolor virginica class.error
#> setosa         50          0         0        0.00
#> versicolor      0         47         3        0.06
#> virginica       0          4        46        0.08

varImpPlot(rf_fit)

11.9 Contoh R: OOB Error

tail(rf_fit$err.rate)

#>               OOB setosa versicolor virginica
#> [495,] 0.04666667      0       0.06      0.08
#> [496,] 0.04666667      0       0.06      0.08
#> [497,] 0.04666667      0       0.06      0.08
#> [498,] 0.04666667      0       0.06      0.08
#> [499,] 0.04666667      0       0.06      0.08
#> [500,] 0.04666667      0       0.06      0.08

1 - tail(rf_fit$err.rate[, "OOB"], 1)

#> [1] 0.9533333

11.10 Contoh R: Sensitivitas `mtry`

mtry_grid <- 1:4
oob_acc <- numeric(length(mtry_grid))

for (i in seq_along(mtry_grid)) {
  rf_i <- randomForest(
    Species ~ .,
    data = iris,
    ntree = 300,
    mtry = mtry_grid[i]
  )
  oob_acc[i] <- 1 - tail(rf_i$err.rate[, "OOB"], 1)
}

data.frame(mtry = mtry_grid, oob_accuracy = round(oob_acc, 4))

#>   mtry oob_accuracy
#> 1    1       0.9467
#> 2    2       0.9467
#> 3    3       0.9600
#> 4    4       0.9533

11.11 Tugas

Jelaskan mengapa bagging menurunkan variance.
Jelaskan mengapa random feature selection diperlukan selain bootstrap.
Bandingkan akurasi single tree dan Random Forest pada dataset yang sama.
Interpretasikan importance plot secara hati-hati.

11.12 Ringkasan

Random Forest adalah ensemble tree yang menurunkan variance melalui bootstrap dan randomization. Secara statistik, ia adalah jawaban langsung terhadap instabilitas single tree.

12 Supervised Learning IV: Support Vector Machine

12.1 Tujuan Akademik Bab

SVM penting karena memperlihatkan bagaimana masalah klasifikasi dapat diformulasikan sebagai optimisasi convex dengan regularisasi eksplisit. Bagi mahasiswa S2, nilai pedagogiknya sangat tinggi karena menggabungkan:

large margin principle,
hinge loss,
constrained optimization,
duality,
kernel methods.

12.2 Hard Margin SVM

Untuk data yang separable,

\[ \min_{\beta,b} \frac{1}{2}\|\beta\|^2 \]

dengan kendala

\[ y_i(\beta^\top x_i + b) \ge 1,\qquad i=1,\ldots,n. \]

Meminimalkan \(\|\beta\|^2\) ekuivalen dengan memaksimalkan margin geometrik.

12.3 Soft Margin SVM

Untuk data nyata yang tidak separable, diperkenalkan slack variables:

\[ \min_{\beta,b,\xi} \frac{1}{2}\|\beta\|^2 + C\sum_{i=1}^{n}\xi_i \]

dengan

\[ y_i(\beta^\top x_i + b) \ge 1-\xi_i, \qquad \xi_i \ge 0. \]

Parameter \(C\) mengatur trade-off antara margin besar dan penalti misclassification.

12.4 Hinge Loss

Soft margin SVM dapat dibaca sebagai minimisasi regularized empirical risk:

\[ \min_{\beta,b} \frac{1}{n}\sum_{i=1}^{n}\max\{0,1-y_i(\beta^\top x_i+b)\} + \lambda \|\beta\|^2. \]

Ini menempatkan SVM langsung ke dalam kerangka statistika pembelajaran.

12.5 Dual Formulation

Lagrangian menghasilkan dual problem:

\[ \max_{\alpha} \sum_{i=1}^{n}\alpha_i -\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} \alpha_i\alpha_j y_i y_j K(x_i,x_j) \]

dengan kendala

\[ 0 \le \alpha_i \le C, \qquad \sum_{i=1}^{n}\alpha_i y_i = 0. \]

Observasi dengan \(\alpha_i>0\) adalah support vectors. Ini sangat penting secara interpretatif: boundary hanya ditentukan oleh sebagian observasi.

12.6 Kondisi KKT dan Margin Efektif

Kondisi KKT memberi klasifikasi observasi ke dalam tiga kategori. Jika \(\alpha_i = 0\), observasi tidak aktif dalam pembentukan boundary. Jika \(0 < \alpha_i < C\), observasi tepat berada di margin. Jika \(\alpha_i = C\), observasi berada di sisi yang salah atau berada di dalam margin lunak. Secara geometris, inilah inti mengapa SVM sering sparse pada representasi dual.

Untuk hyperplane linear,

\[ \beta^\top x + b = 0, \]

lebar margin adalah

\[ \frac{2}{\|\beta\|}. \]

Karena objective hard margin meminimalkan \(\|\beta\|^2/2\), ia ekuivalen dengan memaksimalkan lebar margin.

12.7 Kernel Trick

Dengan kernel \(K(x_i,x_j)\), kita tidak perlu secara eksplisit membangun feature map \(\phi(x)\). Kernel umum:

linear,
polynomial,
radial basis function,
sigmoid.

RBF kernel:

\[ K(x_i,x_j)=\exp(-\gamma \|x_i-x_j\|^2). \]

12.8 Tuning

Parameter penting:

cost,
gamma untuk RBF,
pilihan kernel.

Karena SVM sensitif terhadap skala fitur, standardisasi merupakan bagian esensial dari estimasi, bukan tambahan opsional.

12.9 Contoh R

library(e1071)
iris_bin <- subset(iris, Species != "virginica")
set.seed(123)
idx <- sample(seq_len(nrow(iris_bin)), size = 0.7 * nrow(iris_bin))
train_bin <- iris_bin[idx, ]
test_bin  <- iris_bin[-idx, ]

svm_linear <- svm(Species ~ ., data = train_bin, kernel = "linear", cost = 1, scale = TRUE)
svm_rbf <- svm(Species ~ ., data = train_bin, kernel = "radial", cost = 1, gamma = 0.5, scale = TRUE)

pred_linear <- predict(svm_linear, newdata = test_bin)
pred_rbf <- predict(svm_rbf, newdata = test_bin)

mean(pred_linear == test_bin$Species)

#> [1] 1

mean(pred_rbf == test_bin$Species)

#> [1] 1

12.10 Contoh R: Grid Tuning dan Jumlah Support Vector

grid_svm <- expand.grid(cost = c(0.1, 1, 10), gamma = c(0.1, 0.5, 1))
grid_svm$accuracy <- NA_real_
grid_svm$n_support <- NA_integer_

for (i in seq_len(nrow(grid_svm))) {
  fit_i <- svm(
    Species ~ .,
    data = train_bin,
    kernel = "radial",
    cost = grid_svm$cost[i],
    gamma = grid_svm$gamma[i],
    scale = TRUE
  )
  pred_i <- predict(fit_i, newdata = test_bin)
  grid_svm$accuracy[i] <- mean(pred_i == test_bin$Species)
  grid_svm$n_support[i] <- sum(fit_i$nSV)
}

grid_svm[order(-grid_svm$accuracy, grid_svm$n_support), ]

#>   cost gamma accuracy n_support
#> 3 10.0   0.1        1         6
#> 2  1.0   0.1        1         9
#> 5  1.0   0.5        1        20
#> 6 10.0   0.5        1        20
#> 8  1.0   1.0        1        27
#> 9 10.0   1.0        1        27
#> 1  0.1   0.1        1        41
#> 4  0.1   0.5        1        47
#> 7  0.1   1.0        1        58

12.11 Tugas

Turunkan soft-margin objective sebagai regularized hinge-loss minimization.
Jelaskan arti support vector secara geometris.
Bandingkan kernel linear dan radial pada satu dataset biner.
Jelaskan mengapa scaling wajib pada SVM.

12.12 Ringkasan

SVM adalah model margin-based dengan fondasi optimisasi yang sangat elegan. Dari sudut statistika, ia adalah regularized classifier dengan loss convex dan representasi dual yang kuat.

13 Supervised Learning V: Support Vector Regression

13.1 Tujuan Akademik Bab

SVR memperluas ide margin maksimum ke regresi. Konsep sentralnya adalah epsilon-insensitive loss, yaitu error kecil diabaikan, sedangkan error besar dikenai penalti.

13.2 Formulasi Primal

SVR linear:

\[ f(x)=\beta^\top x + b. \]

Estimasi diperoleh dari

\[ \min_{\beta,b,\xi,\xi^*} \frac{1}{2}\|\beta\|^2 + C\sum_{i=1}^{n}(\xi_i+\xi_i^*) \]

dengan kendala

\[ y_i - (\beta^\top x_i + b) \le \epsilon + \xi_i, \]

\[ (\beta^\top x_i + b) - y_i \le \epsilon + \xi_i^*, \]

\[ \xi_i,\xi_i^* \ge 0. \]

13.3 Epsilon-Insensitive Loss

Loss-nya dapat ditulis sebagai

\[ L_\epsilon(y,f(x)) = \max\{0, |y-f(x)|-\epsilon\}. \]

Ini membuat estimator tidak terlalu bereaksi terhadap error kecil, yang dapat dipandang sebagai bentuk regularisasi robust.

13.4 Dual dan Kernel

Seperti SVM klasifikasi, SVR memiliki representasi dual sehingga kernel dapat dipakai. Prediktor akhirnya berbentuk

\[ f(x) = \sum_{i=1}^{n}(\alpha_i-\alpha_i^*)K(x_i,x) + b. \]

Hanya sebagian observasi dengan koefisien dual tak nol yang berperan sebagai support vectors.

13.5 Struktur Solusi dan Tabung \(\epsilon\)

SVR dapat dibaca sebagai upaya mencari fungsi paling datar yang tetap menempatkan sebanyak mungkin observasi di dalam tabung \(\epsilon\). Jika \(|y_i - f(x_i)| \le \epsilon\), observasi tersebut tidak menambah loss. Hanya observasi di luar tabung yang membentuk support vectors dan memengaruhi solusi dual. Maka, seperti pada SVM klasifikasi, solusi SVR juga cenderung sparse dalam representasi dual.

13.6 Tuning Parameter

Parameter utama:

cost,
epsilon,
gamma untuk kernel radial.

Pemilihan \(\epsilon\) memengaruhi tingkat toleransi error kecil dan sparsity dari solusi dual.

13.7 Contoh R

library(e1071)
set.seed(123)
idx <- sample(seq_len(nrow(mtcars)), size = 0.75 * nrow(mtcars))
train_reg <- mtcars[idx, ]
test_reg  <- mtcars[-idx, ]

svr_fit <- svm(
  mpg ~ wt + hp + disp,
  data = train_reg,
  type = "eps-regression",
  kernel = "radial",
  cost = 10,
  epsilon = 0.1,
  gamma = 0.1,
  scale = TRUE
)

pred_svr <- predict(svr_fit, newdata = test_reg)
sqrt(mean((test_reg$mpg - pred_svr)^2))

#> [1] 2.607152

13.8 Contoh R: Perbandingan SVR dan Regresi Linear

lm_reg <- lm(mpg ~ wt + hp + disp, data = train_reg)
pred_lm_reg <- predict(lm_reg, newdata = test_reg)

c(
  RMSE_SVR = sqrt(mean((test_reg$mpg - pred_svr)^2)),
  RMSE_LM = sqrt(mean((test_reg$mpg - pred_lm_reg)^2)),
  n_support_vectors = length(svr_fit$index)
)

#>          RMSE_SVR           RMSE_LM n_support_vectors 
#>          2.607152          2.119652         21.000000

13.9 Tugas

Jelaskan hubungan antara SVM klasifikasi dan SVR.
Tuliskan bentuk epsilon-insensitive loss dan interpretasikan.
Bandingkan SVR dengan regresi linear pada satu dataset.
Jelaskan pengaruh perubahan epsilon secara konseptual.

13.10 Ringkasan

SVR menunjukkan bahwa large margin principle tidak terbatas pada klasifikasi. Ia adalah metode regularized regression yang fleksibel dan berbasis optimisasi convex.

14 Neural Network dan Deep Learning

14.1 Tujuan Akademik Bab

Neural network perlu diajarkan secara statistik, bukan mistik. Pada inti terdalamnya, neural network adalah fungsi komposisi non-linear berparameter banyak yang diestimasi melalui minimisasi loss dengan gradien.

14.2 Perceptron dan Unit Dasar

Satu neuron menghitung

\[ z = w_0 + \sum_{j=1}^{p}w_jx_j, \qquad a = g(z), \]

dengan \(g\) fungsi aktivasi. Jika \(g\) identitas, seluruh jaringan berlapis akan kolaps menjadi model linear. Karena itu, non-linearity pada aktivasi adalah syarat esensial.

14.3 Multilayer Perceptron

Untuk satu hidden layer:

\[ h = g(W^{(1)}x + b^{(1)}), \]

\[ f(x) = W^{(2)}h + b^{(2)}. \]

Untuk klasifikasi multikelas, output sering memakai softmax:

\[ \hat{p}_k(x) = \frac{\exp(z_k)}{\sum_{\ell=1}^{K}\exp(z_\ell)}. \]

14.4 Loss Function

14.4.1 Regresi

\[ L = \frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2. \]

14.4.2 Klasifikasi

\[ L = -\frac{1}{n}\sum_{i=1}^{n}\sum_{k=1}^{K}y_{ik}\log \hat{p}_{ik}. \]

Loss ini kemudian ditambah regularisasi seperti weight decay jika diperlukan.

14.5 Backpropagation

Backpropagation adalah prosedur menghitung gradien seluruh bobot melalui chain rule. Jika \(\theta\) adalah semua bobot jaringan, maka update umum:

\[ \theta^{(t+1)} = \theta^{(t)} - \eta \nabla_\theta L(\theta^{(t)}). \]

Secara statistik komputasional, ini adalah generalisasi score-based optimization pada model non-linear berlapis.

14.6 Struktur Gradien pada Output Layer

Untuk klasifikasi multikelas dengan softmax dan cross-entropy, gradien pada output layer memiliki bentuk sangat elegan:

\[ \frac{\partial L}{\partial z_k} = \hat{p}_k - y_k. \]

Hasil ini penting secara pedagogis karena menunjukkan bahwa backpropagation bukan prosedur misterius. Ia hanyalah aplikasi sistematis chain rule pada komposisi fungsi yang panjang. Ketika gradien ini dikalikan dengan representasi hidden layer, kita memperoleh update bobot output; sisanya diteruskan ke layer lebih awal dengan aturan rantai yang sama.

14.7 Regularisasi

Neural network memiliki kapasitas tinggi, sehingga regularisasi penting:

weight decay (\(L_2\) penalty),
dropout,
early stopping,
data augmentation,
batch normalization.

14.8 Universal Approximation Theorem

Teorema ini menyatakan bahwa jaringan satu hidden layer yang cukup besar dapat mengaproksimasi banyak fungsi kontinu. Namun teorema itu tidak memberi informasi tentang:

ukuran sampel yang dibutuhkan,
kemudahan optimisasi,
stabilitas generalisasi,
interpretabilitas.

Karena itu, mahasiswa S2 harus berhati-hati agar tidak menyalahartikan teorema ini sebagai jaminan praktis otomatis.

14.9 Deep Learning

Ketika hidden layer banyak, jaringan memasuki rezim deep learning. Keunggulannya adalah pembelajaran representasi bertingkat. Kelemahannya adalah kebutuhan data, tuning, dan komputasi yang jauh lebih besar.

14.10 Contoh R

library(nnet)
set.seed(123)
idx <- sample(seq_len(nrow(iris)), size = 0.7 * nrow(iris))
train_nn <- iris[idx, ]
test_nn  <- iris[-idx, ]

nn_fit <- nnet(
  Species ~ .,
  data = train_nn,
  size = 5,
  decay = 0.01,
  maxit = 500,
  trace = FALSE
)

pred_nn <- predict(nn_fit, newdata = test_nn, type = "class")
mean(pred_nn == test_nn$Species)

#> [1] 0.9777778

14.11 Tugas

Jelaskan mengapa aktivasi non-linear diperlukan.
Tuliskan cross-entropy loss untuk klasifikasi multikelas.
Jelaskan peran weight decay sebagai regularisasi.
Bandingkan neural network kecil dengan satu model klasik pada satu dataset.

14.12 Ringkasan

Neural network adalah model fungsi komposisi yang sangat fleksibel. Bagi statistika, inti pembahasannya adalah loss function, regularisasi, kapasitas fungsi, dan optimisasi gradien.

15 Validation dan Evaluation

15.1 Tujuan Akademik Bab

Bab ini menegaskan bahwa evaluasi bukan tahap administratif, melainkan bagian inti dari metodologi statistik machine learning. Tanpa validasi yang benar, seluruh hasil prediksi bisa optimistik dan tidak dapat dipertanggungjawabkan.

15.2 Train, Validation, Test

Data umumnya dibagi menjadi:

training set untuk estimasi parameter,
validation set untuk tuning hyperparameter,
test set untuk evaluasi akhir.

Jika test set dipakai berulang-ulang untuk memilih model, maka estimasi performa menjadi bias optimistik.

15.3 K-Fold Cross-Validation

Jika data dibagi ke \(K\) fold:

\[ \widehat{CV} = \frac{1}{K}\sum_{k=1}^{K}\text{Err}_k. \]

CV memberi estimasi generalization error yang lebih stabil daripada satu split. Namun ia juga memiliki variasi acak, sehingga pelaporannya sebaiknya disertai konteks dan, bila mungkin, pengulangan.

15.4 Nested Cross-Validation

Jika hyperparameter dituning di dalam proses evaluasi, nested CV adalah prosedur yang lebih bersih:

outer loop untuk evaluasi,
inner loop untuk tuning.

Pada level S2, mahasiswa sebaiknya mengenal konsep ini sebagai standar metodologis yang lebih kuat.

15.5 Metrik Klasifikasi

Untuk klasifikasi biner:

\[ \text{Accuracy} = \frac{TP+TN}{TP+TN+FP+FN}, \]

\[ \text{Precision} = \frac{TP}{TP+FP}, \]

\[ \text{Recall} = \frac{TP}{TP+FN}, \]

\[ F_1 = 2\frac{\text{Precision}\cdot\text{Recall}}{\text{Precision}+\text{Recall}}. \]

Pada data tidak seimbang, accuracy sering menyesatkan, sehingga recall, precision, F1, ROC-AUC, dan PR-AUC lebih relevan.

15.6 Kalibrasi dan Proper Scoring Rule

Jika model menghasilkan probabilitas, kita ingin probabilitas itu terkalibrasi. Proper scoring rules seperti log-loss dan Brier score penting karena menghukum probabilitas yang terlalu percaya diri tetapi salah.

\[ \text{LogLoss} = -\frac{1}{n}\sum_{i=1}^{n} \left[ y_i \log \hat{p}_i + (1-y_i)\log(1-\hat{p}_i) \right]. \]

15.7 Metrik Regresi

\[ \text{MAE} = \frac{1}{n}\sum_{i=1}^{n}|y_i-\hat{y}_i|, \]

\[ \text{RMSE} = \left(\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat{y}_i)^2\right)^{1/2}, \]

\[ R^2 = 1-\frac{\sum_{i=1}^{n}(y_i-\hat{y}_i)^2} {\sum_{i=1}^{n}(y_i-\bar{y})^2}. \]

RMSE lebih sensitif pada outlier besar; MAE lebih robust. Pilihan metrik harus mengikuti tujuan substantif.

15.8 Time Series Evaluation

Pada time series, splitting acak umumnya salah. Evaluasi harus menjaga urutan waktu, misalnya melalui rolling forecast origin atau expanding window.

15.9 Optimism, Resubstitution, dan Bootstrap

Error training hampir selalu optimistik karena model dievaluasi pada data yang dipakai untuk estimasi. Jika \(\widehat{Err}_{train}\) terlalu kecil dibanding error generalisasi, mahasiswa harus bertanya apakah model sedang menangkap struktur yang stabil atau hanya noise sampel.

Bootstrap dapat dipakai untuk mengukur kestabilan performa atau koefisien. Walaupun tidak menggantikan test set eksternal, bootstrap membantu memahami variabilitas estimator dan metrik evaluasi, terutama ketika ukuran sampel terbatas.

15.10 Contoh R: Cross-Validation Sederhana

set.seed(123)
k <- 5
fold_id <- sample(rep(1:k, length.out = nrow(iris)))
acc <- numeric(k)

for (i in 1:k) {
  train_cv <- iris[fold_id != i, ]
  test_cv  <- iris[fold_id == i, ]

  fit_cv <- glm(Species == "setosa" ~ Sepal.Length + Sepal.Width +
                  Petal.Length + Petal.Width,
                data = train_cv, family = binomial())

  prob_cv <- predict(fit_cv, newdata = test_cv, type = "response")
  pred_cv <- ifelse(prob_cv > 0.5, TRUE, FALSE)
  acc[i] <- mean(pred_cv == (test_cv$Species == "setosa"))
}

mean(acc)

#> [1] 1

15.11 Contoh R: Metrik Klasifikasi Manual

actual <- c(1,1,1,0,0,0,1,0)
pred   <- c(1,1,0,0,0,1,1,0)

TP <- sum(actual == 1 & pred == 1)
TN <- sum(actual == 0 & pred == 0)
FP <- sum(actual == 0 & pred == 1)
FN <- sum(actual == 1 & pred == 0)

accuracy  <- (TP + TN) / (TP + TN + FP + FN)
precision <- TP / (TP + FP)
recall    <- TP / (TP + FN)
f1        <- 2 * precision * recall / (precision + recall)

c(accuracy = accuracy, precision = precision, recall = recall, f1 = f1)

#>  accuracy precision    recall        f1 
#>      0.75      0.75      0.75      0.75

15.12 Contoh R: Brier Score dan Log-Loss

fit_prob <- glm(
  Species == "setosa" ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
  data = iris,
  family = binomial()
)

prob_hat <- predict(fit_prob, type = "response")
y_bin <- as.integer(iris$Species == "setosa")

brier_score <- mean((y_bin - prob_hat)^2)
log_loss <- -mean(y_bin * log(prob_hat) + (1 - y_bin) * log(1 - prob_hat))

c(Brier = brier_score, LogLoss = log_loss)

#>                         Brier                       LogLoss 
#> 0.000000000000000000005012336 0.000000000010980096094116084

15.13 Contoh R: Repeated Cross-Validation Ringan

set.seed(99)
cv_rep <- replicate(20, {
  fold_id_rep <- sample(rep(1:5, length.out = nrow(iris)))
  acc_rep <- numeric(5)

  for (k_rep in 1:5) {
    train_rep <- iris[fold_id_rep != k_rep, ]
    test_rep  <- iris[fold_id_rep == k_rep, ]

    fit_rep <- glm(
      Species == "setosa" ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
      data = train_rep,
      family = binomial()
    )

    prob_rep <- predict(fit_rep, newdata = test_rep, type = "response")
    pred_rep <- prob_rep > 0.5
    acc_rep[k_rep] <- mean(pred_rep == (test_rep$Species == "setosa"))
  }

  mean(acc_rep)
})

c(mean_cv_accuracy = mean(cv_rep), sd_cv_accuracy = sd(cv_rep))

#> mean_cv_accuracy   sd_cv_accuracy 
#>                1                0

15.14 Tugas

Jelaskan mengapa nested CV lebih bersih daripada CV tunggal untuk tuning model.
Bandingkan MAE dan RMSE pada satu model regresi.
Jelaskan situasi di mana recall lebih penting daripada precision.
Rancang skema evaluasi yang tepat untuk satu dataset time series.

15.15 Ringkasan

Validasi dan evaluasi adalah bagian inti dari inferensi prediktif. Pemilihan metrik, skema splitting, dan prosedur tuning harus diselaraskan dengan struktur data dan tujuan keputusan.

16 UAS: Project Kelompok 2 Orang

16.1 Tujuan Akademik

UAS merupakan sintesis seluruh mata kuliah. Pada tahap ini mahasiswa harus menunjukkan bukan hanya pemahaman terhadap algoritma, tetapi juga kematangan metodologis dalam merancang, mengestimasi, mengevaluasi, dan mengomunikasikan sebuah studi machine learning.

16.2 Struktur Studi yang Diharapkan

Laporan akhir sebaiknya memuat:

rumusan masalah,
definisi target statistik,
sumber dan kualitas data,
justifikasi preprocessing,
formulasi model kandidat,
strategi tuning dan validasi,
evaluasi dan interpretasi,
keterbatasan dan saran pengembangan.

16.3 Standar untuk Mahasiswa S2

Pada tingkat magister, penilaian akan lebih tinggi jika mahasiswa mampu:

menghubungkan model ke teori statistik,
menjelaskan objective function,
menilai trade-off antara akurasi dan interpretabilitas,
menunjukkan kesadaran terhadap leakage, imbalance, dan calibration,
menyajikan analisis yang reprodusibel.

16.4 Kriteria Penilaian

problem framing,
kualitas data understanding,
ketepatan model,
kualitas validasi,
kualitas interpretasi,
reproducibility,
kualitas presentasi akademik.

16.5 Ringkasan

UAS adalah latihan menyusun studi machine learning yang utuh, defensibel, dan setara dengan mini-penelitian statistik terapan.

Referensi Inti

Hastie, T., Tibshirani, R., and Friedman, J. The Elements of Statistical Learning.
James, G., Witten, D., Hastie, T., and Tibshirani, R. An Introduction to Statistical Learning.
Bishop, C. M. Pattern Recognition and Machine Learning.
Murphy, K. P. Machine Learning: A Probabilistic Perspective.
Kuhn, M., and Johnson, K. Applied Predictive Modeling.
Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization Paths for Generalized Linear Models via Coordinate Descent.
Breiman, L. (2001). Random Forests.
Vapnik, V. Statistical Learning Theory.

Penutup

Versi buku ini disusun agar machine learning dibaca sebagai bagian dari statistika modern: model, loss function, regularisasi, estimasi parameter, validasi, dan interpretasi. Jika buku ini akan diperdalam lebih lanjut, pengembangan paling tepat adalah memperpanjang tiap bab dengan studi kasus nyata, pembuktian yang lebih formal, dan pembahasan inferensi modern untuk model machine learning.

Ebook Machine Learning untuk S2 Statistika

Formulasi Model, Estimasi Parameter, dan Implementasi dengan R

I Gede Nyoman Mindra Jaya dan Irlandia Ginanjar

13 April 2026

Prakata

Cara Membaca Buku Ini

0.1 Sasaran

0.2 Notasi Umum

0.3 Paket R yang Sering Dipakai

0.4 Peta Buku

1 Perkenalan Machine Learning

1.1 Tujuan Akademik Bab

1.2 Machine Learning sebagai Masalah Statistik

1.3 Apa yang Sebenarnya Diestimasi?

1.3.1 Regresi

1.3.2 Klasifikasi

1.3.3 Clustering dan Reduksi Dimensi

1.4 Perspektif Inferensi versus Prediksi

1.5 Loss Function dan Keputusan Statistik

1.6 Population Risk, Empirical Risk, dan Generalization Gap

1.7 Overfitting, Underfitting, dan Bias-Variance Trade-Off

1.8 Bayes Rule sebagai Batas Teoretis

1.9 No Free Lunch dan Pentingnya Struktur Masalah

1.10 Tiga Pilar Analisis Machine Learning untuk Statistika

1.10.1 Pilar 1: Model

1.10.2 Pilar 2: Estimasi

1.10.3 Pilar 3: Evaluasi

1.11 Alur Kerja Analisis untuk Mahasiswa S2

1.12 Contoh R: Train-Test Split dan Baseline Awal

1.13 Tugas

1.14 Ringkasan

2 Analisis Regresi dan Gradient Descent

2.1 Posisi Regresi Linear dalam Machine Learning

2.2 Formulasi Model

2.3 OLS sebagai M-Estimator

2.4 Normal Equations dan Solusi Tertutup

2.5 Derivasi Formal Normal Equations

2.6 Asumsi Model dan Implikasinya

2.7 Regresi untuk Inferensi versus Regresi untuk Prediksi

2.8 Basis Expansion dan Linearity dalam Parameter

2.9 Gradient Descent sebagai Prosedur Estimasi Numerik

2.10 Konvergensi Gradient Descent

2.10.1 Batas Learning Rate dari Spektrum Hessian

2.11 Batch, Stochastic, dan Mini-Batch

2.11.1 Batch Gradient Descent

2.11.2 Stochastic Gradient Descent

2.11.3 Mini-Batch Gradient Descent

2.12 Kondisi Numerik dan Standardisasi

2.13 Regresi Linear sebagai Baseline Wajib

2.14 Simulasi Sederhana untuk Menunjukkan Bias dan Nonlinearitas

2.15 Contoh R: Regresi Linear dan Diagnostik Prediktif

2.16 Contoh R: Gradient Descent Manual dan Jejak Konvergensi

2.17 Contoh R: Closed Form versus Gradient Descent pada Data Hold-Out

2.18 Contoh R: Simulasi Nonlinearitas

2.19 Tugas

2.20 Ringkasan

3 Ridge, Lasso, dan Elastic Net

3.1 Mengapa Regularisasi Diperlukan

3.2 Ridge Regression

3.2.1 Derivasi Solusi Tertutup Ridge

3.2.2 Representasi Spektral Ridge

3.2.3 Pandangan Bayesian

3.3 Lasso

3.3.1 KKT Conditions

3.3.2 Coordinate Descent dan Soft Thresholding

3.3.3 Soft Thresholding pada Kasus Ortonormal

3.3.4 Pandangan Bayesian

3.4 Elastic Net

3.5 Degrees of Freedom dan Kompleksitas Model

3.6 Estimasi Parameter

3.6.1 Ridge

3.6.2 Lasso dan Elastic Net

3.7 Pemilihan Parameter Regularisasi

3.8 Hubungan dengan Bias-Variance Trade-Off

3.9 Kapan Ridge, Lasso, atau Elastic Net Dipakai

3.9.1 Ridge

3.9.2 Lasso

3.9.3 Elastic Net

3.10 Contoh R: Jalur Regularisasi dan Perbandingan Koefisien

3.11 Contoh R: Membandingkan RMSE In-Sample secara Sederhana