Moving Beyond Linearity
Pendahuluan
Pembelajaran kali ini terkait pemodelan pola hubungan antara dua variabel numerik dimana antar kedua variabel tersebut mempunyai hubungan yang tak linier.
Library ISLR
library(ISLR)Import Data
?Auto## starting httpd help server ... done
attach(Auto)
head(cbind(mpg, horsepower))## mpg horsepower
## [1,] 18 130
## [2,] 15 165
## [3,] 18 150
## [4,] 16 150
## [5,] 17 140
## [6,] 15 198
Eksplorasi Data
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon") Ini merupakan scatter plot dari 392 kendaraan dimana dari masing-masing kendaraan diukur dua peubah/variabel, yaitu kekuatan mesin (horse power pada sumbu X) dan berapa jauh kendaraan tersebut berjalan untuk setiap galon bahan bakarnya (mile per gallon pada sumbu y).
Secara umum,kendaraan dengan horse power yang besar membutuhkan bahan bakar yang banyak yang dicirikan jarak tempuh yang semakin pendek ketika konsumsi bahan bakarnya satu galon. Yang paling kanan, kekuatan mesinnya besar, konsumsi bahan bakar besar karena jarak tempuhnya kecil.
Jika digunakan model linier, diperoleh persaman garis lurus: y=39.936 - 0.158x. Koefisien x negatif karena makin ke kanan makin turun sehingga slope (gradient) nya negatif. Hal ini tidak sesuai dengan apa yang kita inginkan.
Beberapa hal yang kurang menyenangkan dari model linier pada kasus ini, yaitu:
Model linier yang kita gunakan tidak mengikuti pola hubungan pada data. Sebelah kiri curam, sebelah kanan cenderung landai. Kurang pas sehingga prediksinya tidak cukup baik.
Dapat diperoleh nilai penduga peubah respon yang “tidak masuk akal”. Jika garis linier diteruskan, untuk kekuatan mesin tertentu, akan diperoleh respon yang negatif sehingga tidak mungkin terjadi.
Sifat dan perilaku sisaan cenderung tidak sesuai harapan. Jika wilayah dibagi menjadi tiga bagian, bagian yang tengah, eror cenderung negatif karena garis berwarna orange (dugaan) cenderung ada di atas titik2 yang sebenarnya. Sedangkan wilayah yang kanan sisaan cenderung positif karena garis berwarna orange (dugaan) cenderung ada di bawah titik2 yang sebenarnya. Padahal biasanya lebih disukai sisaan.residual dari model ada yang positif dan negatif dengan nilai harapan nol. Di kasus ini, jika disekat,ada wilayah yang nilai harapan positif dan negatif sehingga tidak sesuai harapan.
Pendekatan Pola Hubungan Tidak Linier
Jika kita mempunyai data dengan pola yang cenderung tidak linier, ada beberapa pendekatan yang bisa digunakan seperti regresi polinomial atau fungsi tangga/step function (pada nilai2 tertentu penduganya akan konstan) atau regresi spline. Bisa juga pendekatan lain seperti smoothing atau regresi lokal atau pohon regresi.
Regresi Polinomial
Model dari regresi polinomial dapat dinyatakan dalam bentuk:
d (derajat polynomial) yang digunakan biasanya hingga d=4. d yang terlalu besar akan membuat kurva berbentuk tidak teratur dan tampak “liar”, terutama pada bagian ujungnya.
polinom = lm(mpg ~ poly(horsepower,4))
prediksi = predict( polinom , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="Regresi Polinomial Ordo 4")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2)Contoh d=7 yang menghasilkan kurva berbentuk tidak teratur dan tampak “liar”, terutama pada bagian ujungnya.
polinom = lm(mpg ~ poly(horsepower,7))
prediksi = predict( polinom , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="Regresi Polinomial Ordo 7")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2) Dari gambar di atas terlihat bahwa ada pola yang tidak teratur pada ujung kiri dan ujung kanan kurva dengan d=7 tersebut.
Fungsi Tangga (Piecewise Constant)
Pada nilai tertentu bentuknya mendatar (konstan), kondisi yang lain juga konstan dengan nilai yang berbeda , dst. Cara mendapatkan fungsi tangga:
tentukan K buah titik potong: c1 < c2 < … < ck
definisikan K+1 peubah baru berikut
Jika ada K titik potong ,akan ada K+1 selang nilai. Setiap selang nilai definisikan peubah baru. Selang pertama berindeks nol adalah fungsi indikator x<c1 artinya jika x<c1 nilainya 1, jika tidak nilainya 0. Selang kedua berindeks 1 adalah fungsi indikator c1<x<c2 artinya jika c1<x<c2 nilainya 1, jika tidak nilainya 0, dst. I() adalah fungsi indikator, nilainya 1 jika syaratnya (yang ada dalam tanda kurung) terpenuhi, dan 0 jika tidak. Oleh karena itu, akan diperoleh K+1 peubah dummy (dummy variable) yang nilainya 1 dan 0 saja. Tentu saja peubah dummy tersebut jika kita jumlahkan nilainya akan sama dengan 1.
C0(X) +C1(X)+ … + Ck(X) =1
- Berdasarkan peubah dummy, susun model dari fungsi tangga. Model:
Maka, metode kuadrat terkecil dapat digunakan untuk mencari koefisien dari penduga beta.Berikut visualisasi fungsi tangga.
tangga = lm(mpg ~ cut(horsepower,12))
prediksi = predict(tangga , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="Fungsi Tangga dengan 11 knots")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2) - Untuk memperoleh prediksi nilai peubah respon untuk setiap selang nilai/sekat dari peubah x tadi tidak lain adalah rata-rata pada sekatan dari selang yang mendatar (horizontal) untuk peubah tersebut.
Gambar di atas dibagi menjadi 12 selang nilai.Sekat yang pertama adalah rata-rata dari amatan pada selang pertama, dst. Jika nilai amatan pada selang tersebut umumnya kecil, maka rata2nya juga kecil. Begitu sebaliknya.
- Bj adalah selisih rata-rata antara selang ke-j dengan selang ke-0 (selang yang paling kiri). Bentuk fungsi tangga ditentukan oleh banyaknya titik potong.
#K=5
tangga = lm(mpg ~ cut(horsepower,5))
prediksi = predict(tangga , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="Fungsi Tangga dengan 4 knots")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2) Dari kurva di atas dengan K=5 terlihat bahwa jika titik potongnya sedikit, maka fungsi tangga ini penurunan darisatu sekat ke sekat lain perubahannya banyak dan anak tangganya lebar-lebar.
Sedangkan jika dilihat dari kurva di bawah dengan K=19 nampak bahwa semakin banyak titik potong, maka kurva dari fungsi tangga ini akan semakin mulus (anak tangganya kecil-kecil karena anak tangganya banyak), membentuk pola yang mengikuti data tetapi ada bentuk2 yang menjadi tidak teratur.
#K = 19
tangga = lm(mpg ~ cut(horsepower,19))
prediksi = predict(tangga , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="Fungsi Tangga dengan 18 knots")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2)Model Fungsi Basis (Basis Function Model)
Secara umum, model fungsi basis adalah:
bj(X) merupakan fungsi basis.
Pada pembuatan matriks rancangan, peubah penjelas Xi ditransformasi ke dalam K buah fungsi basis.
Dalam regresi polinomial, jika dipandang sebagai fungsi basis, maka bj(X) = X^j
Sedangkan Dalam fungsi tangga (piecewise constant), bj(X)= Ci(X)=I(cj<=X< cj+1). c adalah nilai-nilai titik potong atau knots nya.
Sehingga, regresi polinomial dan fungsi tangga merupakan bentuk khusus dari model fungsi basis. Jadi, funsi basis sendiri bisa kita ubah dalam bentuk apapun.
Regresi Spline
Regresi spline pada dasarnya merupakan gabungan dari dua pendekatan yang sebelumnya telah dibahas yaitu gabungan regresi polinomial dan fungsi piecewise pada titik potong (knots) c1,c2,…,ck.
Yang dilakukan dalam regresi spline:
- Jika wilayah X tadi dipartisi/sekat, maka dilakukan pemodelan fungsi polinomial di setiap sekatan dan diupayakan terhubung/tergabung secara mulus (smooth).
Regresi Polinomial per-bagian (Piecewise Polynomial): Daerah peubah penjelas dibagi berdasarkan nilai titik potong (knot) dan Di masing-masing bagian/partisi dilakukan pemodelan regresi polinomial
Misal, ada satu knot pada titik X = 120, dan regresi kubik di masing-masing bagian. Sehingga modelnya dibagi menjadi 2, yaitu:
Yang warna biru pada gambar di atas adalah bentuk kurva dari fungsi regresi polinomial yang dihasilkan. Terlihat bahwa pada titik knot 120, modelnya tidak kontinu. Jadi, fungsi yang pertama dan kedua tidak terhubung dengan mulus (diskontinu).Inilah yang akan ditangani dengan menggunakan regresi spline.
- Fungsi basis pada regresi spline berordo M dengan K buah knot adalah K + M -1 buah yang masing-masing adalah
Sebanyak M-1 adalah fungsi polinomial.
K buah fungsi basis yang lain adalah truncated power function. j adalah titik knots
Di atas adalah fungsi basis dari regresi spline.
Untuk Cubic Spline (Pangkat 3) m-1 = 3
m=4
Banyaknya fungsi Basis: K + M - 1 = K + 4 - 1 = K + 3 yaitu:
Berikut grafik cubic splines (M=4) dan satu knot yaitu pada posisi c1=120. Pada 120, fungsinya telah terhubung secara mulus, kontinu, dan dapat diturunkan.
library(splines)
regspline = lm(mpg ~ bs(horsepower, knots=c(120)))
prediksi = predict(regspline , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="cubic splines dengan knots 120")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2)
abline(v=c(120), col="grey50", lty=2) Jika knot nya diperbanyak, cubic splines:
library(splines)
regspline = lm(mpg ~ bs(horsepower, knots=c(80, 120, 160, 200)))
prediksi = predict(regspline , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="cubic splines dengan knots 80, 120, 160, 200")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2)
abline(v=c(80, 120, 160, 200), col="grey50", lty=2)Jika knots nya kita perbanyak menjadi 6 menghasilkan gambar berikut:
library(splines)
regspline = lm(mpg ~ bs(horsepower, knots=c(75, 100, 125, 150,175, 200)))
prediksi = predict(regspline , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="cubic splines dengan knots 75, 100, 125, 150,175, 200")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2)
abline(v=c(75, 100, 125, 150,175, 200), col="grey50", lty=2) Terlihat bahwa jika makin banyak knots yang digunakan, maka jika dilihat ujung-ujung selang, bentuk kurva spline nya cenderung “liar” atau memiliki ragam yang besar.Jadi, ujung kiri dan ujung kanan, bentuk kurva nya nya cenderung tidak stabil.
Kalau kita sebut knot sebagai: c0,c1,c2,…,ck,ck+1, sekatan/partisi yang terluar adalah:
x<c0 (bagian kiri)
dan
x > ck+1 (bagian kanan)
Istilah yang biasa digunakan, yaitu: c0 dan ck+1 (knots yang diluar) disebut sebagai boundary knots sedangkan c1,…,ck disebut ordinary knots.
Regresi Natural Spline
Jika cubic splines menghasilkan kurva yang liar di kiri dan di kanan, maka untuk mengatasinya kita gunakan natural splines.
Natural Cubic Splines (Fungsinya Kubik)
Upaya untuk mengurangi “keliaran” atau keragaman tinggi pada partisi ujung, baik ujung kiri maupun ujung kanan
Yang dilakukan:
- Fungsi basis di kiri dan di kanan adalah fungsi linier.
Pada grafik cubic splines (M=4) dengan empat knot yaitu 80,120,160,200 yang telah digambarkan sebelumnya, terlihat bahwa ujung kiri dan ujung kanan kurva terlihat “liar”. Ketika digambarkan dengan menggunakan natural cubic spline pada kurva di bawah menghasilkan kurva di ujung kiri dan ujung kanan yang cenderung lebih stabil.
library(splines)
regspline = lm(mpg ~ ns(horsepower, knots=c(80, 120, 160, 200)))
prediksi = predict(regspline , data.frame(horsepower), interval="predict" )
plot(horsepower, mpg, pch=19, col="coral",
xlab="Horse Power", ylab="mile per gallon",
main="natural cubic splines dengan knots 80, 120, 160, 200")
ix <- sort(horsepower,index.return=T)$ix
lines(horsepower[ix], prediksi[ix],
main = "Polynomial Regression", col="blue",
type="l", lwd=2)
abline(v=c(80, 120, 160, 200), col="grey50", lty=2)Smoothing Spline
Pada data sebelumnya terlihat bahwa terdapat hubungan yang tidak linier antara x dan y. Bentuk fungsi pola hubungan yang tidak linier ini bisa dinotasikan dengan g(x). Sifat yang diinginkan dari g(x) adalah nilainya dekat dengan y (variabel asal) dan bentuk dari g(x) tidak jumpy (wiggly) atau berkelok-kelok. Jadi, cenderung mulus.
Untuk mendapatkan sifat yang diinginkan dari g(x), bisa dilakukan dengan cara mencari fungsi g yang meminimumkan fungsi berikut:
Bagian pertama dari fungsi ini adalah selisih nilai y dengan g(x) sehingga diupayakan g(x) sedekat mungkin dengan y. Jadi, ini ukuran eror dari model yang kita miliki.
Bagian kedua selain lambda adalah besaran yang kita peroleh untuk mengontrol bentuk “jumpy”. Turunan pertama suatu fungsi g’(x) adalah slope/kemiringan dan turunan kedua suatu fungsi g’’(x) adalah perubahan slope. Jadi, perubahan slope diinginkan kecil agar bentuknya tidak bergerigi/tidak berkelok-kelok.
Besarnya lambda adalah smoothing parameter yang mengontrol mana yang penting antara kedua hal ini.Jika lambda = 0, maka cenderung memprioritaskan bagian pertama. Tapi, jika lambda nya makin besar maka semakin ingin kita membuat fungsinya tidak berkelok-kelok atau semakin besar lambda, maka fungsinya akan semakin mulus bentuknya.
Berikut ilustrasi efek pemilihan nilai parameter pemulusan lambda:
Terlihat bahwa semakin kecil lambda, polanya semakin berkelok-kelok. Semakin besar lambda, polanya semakin mulus.
Nilai lambda dapat ditentukan menggunakan proses cross validation atau validasi silang. Salah satunya adalah LOOCV (leave-one-out cross validation), yaitu menyisihkan satu per satu amatan di setiap iterasi sebagai amatan data validasi (data testing), dan selanjutnya menghitung RSS (residual sum of squares) dari hasil prediksi.
Residual Sum of Square dari CV adalah
yi: nilai aktual dari peubah respon.
g.cap : prediksi. pangkat (-1) berarti untuk mendapatkan g kita membuang satu amatan, yaitu amatan ke-i.
Ilustrasi nya untuk menentukan lambda:
Jika lambda kita ubah-ubah, maka Residual Sum of Square dari CV akan berbeda. Kita mencari pada lambda berapa kita menemukan Residual Sum of Square dari CV yang paling kecil. Ini disebut sebagai lambda optimal. Lambda optimal inilah yang kita gunakan untuk membuat model smoothing splines.
Local Regression
Metode yang paling populer adalah LOESS.
Jika kita mempunyai data, untuk setiap nilai x misal x=100, kita tentukan titik-titik terdekat dari x=100 tersebut. Kita gunakan data yang berwarna orange untuk menghasilkan model regresi linier secara lokal karena melibatkan tidak semua amatan tetapi amatan terdekat dari titik x=100. Jika didapat garis regresi lokal, dapat diperoleh prediksi untuk x=100 berdasarkan regresi lokal.
Tahapan tersebut diulang-ulang untuk x yang lain, misal x=110, x=120 dll. Jika ini dikerjakan, akan diperoleh kurva dari model local regression.
Kurva dan tahapan dari model local regression berikut:
Kurva warna ungu didapat Jika kita ulang nilai x sampai dengan 50 atau 200.
Tentukan dulu data yang terdekat dari x, misal x=100.
Gunakan data pada rentang tersebut saja untuk mendapat model regresi linier dengan menggunakan weighted least square.Jika titik lebih dekat dengan x=100, maka diberikan bobot yang lebih besar. Makin jauh dari x=100, bobotnya makin kecil menuju 0.
Jika telah diperoleh model regresi, maka nilai prediksi/pengepasan pada titik tersebut tinggal dimasukkan pada model regresi linier lokal tadi. Maka, akan diperoleh titik-titik yang berwarna ungu seperti kurva di atas.
Penentuan banyaknya data yang digunakan untuk menghasilkan regresi lokal tadi, s=k/n adalah berlaku sebagai parameter pemulusan. Jika s semakin besar (yang digunakan untuk menghasilkan regresi lokal makin banyak ), maka semakin mulus. Jika s semakin kecil, regresi lokal melibatkan amatan yang semakin sedikit, maka bentuk kurvanya semakin jumpy atau berkelok-kelok karena perubahannya semakin banyak.
Berikut contoh jika kita gunakan lebar 20 dan lebar 50. Yang lebarnya 50, kurvanya semakin mulus.
Referensi
Sartono, B.. (November 1, 2020). Moving Beyond Linearity. Retrieved from https://rpubs.com/bagusco/taklinear
P8 & P9 STA581_Sains Data, Mahasiswa Pascasarjana Statistika dan Sains Data,IPB University, reniamelia@apps.ipb.ac.id↩︎