Memakai Linear Regression Di R


Linear regression adalah sebuah teknik klasik di statistika untuk mempelajari hubungan antar-variabel dan memprediksi masa depan. Walaupun tidak seakurat teknik yang lebih modern, kelebihan linear regression adalah mudah dimengerti dan tidak mensyaratkan data harus dalam bentuk tertentu. Sebagai latihan, pada artikel ini saya akan memakai data bengkel motor untuk mencoba mencari tahu pengaruh antara jumlah pekerjaan OR+ (ganti oli) dan CS (servis) pada L/R (penggantian suku cadang). Secara logika, bukankah motor yang rutin dirawat berkala akan semakin awet (sedikit mengganti suku cadang)? Tapi ada juga rumor bahwa semakin sering motor di-servis, maka motor tersebut akan semakin ‘tidak awet’? Apakah rumor ini benar? Saya akan mencoba membuktikannya melalui perintah R!

Langkah pertama yang saya lakukan adalah membaca data dari file CSV dengan menggunakan perintah berikut ini:

> servis <- read.csv("servis.csv")
> head(servis)
  nomorKB   tipeMotor jenisPekerjaan
1  XXXXXX Supra X 125            L/R
2  YYYYYY     Supra X             CS
3  ZZZZZZ       Grand            L/R
4  AAAAAA     Supra X             CS
5  BBBBBB     Megapro            L/R
6  CCCCCC        Beat             CS

Saya terlebih dahulu perlu mengolah data tersebut dengan melakukan agregasi (mirip seperti operasi GROUP BY pada SQL) berdasarkan nomorKB yang sama. Cara yang paling mudah adalah dengan memakai ddply() dari package plyr seperti yang terlihat pada perintah berikut ini:

> require(plyr)
> data <- ddply(servis, c('nomorKB'), summarise,
+   qty.or=length(jenisPekerjaan[jenisPekerjaan=='OR+']),
+   qty.cs=length(jenisPekerjaan[jenisPekerjaan=='CS']),
+   qty.lr=length(jenisPekerjaan[jenisPekerjaan=='L/R'])
+ )
> head(data)
  nomorKB qty.or qty.cs qty.lr
1  XXXXXX      1      0      0
2  YYYYYY      0      1      0
3  ZZZZZZ      0      1      2
4  AAAAAA      0      8      2
5  BBBBBB      1      2      3
6  CCCCCC      0      1      0

Berbeda dengan kebanyakan bahasa, pada R, nama variabel boleh mengandung titik seperti qty.or atau qty.cs.

Untuk melakukan fitting ke persamaan linear, saya dapat menggunakan function lm() seperti berikut ini:

> f <- lm(qty.lr ~ qty.or + qty.cs, data=data)
> summary(f)

Call:
lm(formula = qty.lr ~ qty.or + qty.cs, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5087 -0.6134 -0.4590  0.5410 10.3349 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.45904    0.04320  10.626   <2e-16 ***
qty.or       0.48561    0.02029  23.928   <2e-16 ***
qty.cs       0.15435    0.01836   8.408   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.356 on 1482 degrees of freedom
Multiple R-squared:  0.322, Adjusted R-squared:  0.3211 
F-statistic: 351.9 on 2 and 1482 DF,  p-value: < 2.2e-16

Formula yang saya berikan pada lm() adalah qty.lr ~ qty.or + qty.cs. Ini menunjukkan bahwa qty.or dan qty.cs adalah predictor, sementara itu qty.lr adalah response. Linear regression yang memiliki lebih dari satu predictor seperti ini disebut juga sebagai multiple linear regression.

Karena ini adalah multiple linear regression, nilai yang pertama kali saya perhatikan pada hasil summary() adalah nilai Adjusted R-squared. Pada statistika, R-squared disebut sebagai coefficient of determination. Nilai R-squared yang berkisar mulai dari 0 hingga 1 merupakan sesuatu yang penting karena ia menunjukkan seberapa ‘benar’ model yang saya pakai. Nilai 0.3211 menunjukkan bahwa 32,11% dari data dapat diprediksi oleh persamaan linear yang saya pakai. Berdasarkan jenis data yang saya pakai, angka ini adalah sesuatu yang bisa diterima.

Untuk melihat informasi kesalahan secara visual, saya dapat menggunakan perintah plot() seperti berikut ini:

> plot(f)
Hit <Return> to see next plot: 
Plot residuals versus fitted

Plot residuals versus fitted

Semakin dekat sebuah titik dengan nilai 0 di sumbu Y menunjukkan prediksi model yang semakin baik. Nilai yang berbeda dari prediksi disebut sebagai residual.

Saya dapat melakukan cross-validation untuk memeriksa seberapa tepat model linear yang saya pakai. Sebagai latihan, saya akan memakai metode Leave-One-Out Cross-Validation (LOOCV). LOOCV adalah cross-validation dilakukan dengan mengambil 1 baris data (observation) untuk diuji kebenarannya (data validasi) dimana sisa data lainnya dipakai untuk menghitung nilai koefisien (data training). Proses ini dilakukan berulang kali mulai dari baris pertama hingga baris terakhir sebagai data validasi. Pada R, saya dapat melakukan LOOCV dengan menggunakan cv.glm() di package boot. Karena method ini hanya menerima hasil dari glm(), saya perlu menghitung ulang output dari lm() dalam bentuk glm() (seharusnya identik). Berikut adalah perintah yang saya berikan:

> f <- glm(qty.lr ~ qty.or + qty.cs, data=data)
> coef(f)
(Intercept)      qty.or      qty.cs 
  0.4590442   0.4856098   0.1543520 
> require(boot)
> cv <- cv.glm(data, f)
> cv$delta
[1] 1.863092 1.863082

Nilai yang dikembalikan adalah nilai rata-rata dari mean squared error (MSE) untuk n cross-validation yang telah dilakukan. Nilai MSE yang semakin mendekati nol (0) menunjukkan model yang semakin akurat.

Saya juga dapat menggunakan metode subset selection untuk mencari variabel apa saja yang perlu disertakan pada persamaan guna memperoleh nilai MSE paling kecil. Sebagai contoh, kali ini saya akan membuat dataframe yang berisi qty untuk seluruh jenis pekerjaan yang ada. Agar gampang, saya akan membuat sebuah function baru seperti berikut ini:

buatQtyJenisPekerjaan <- function(.data) {
  t <- xtabs(~jenisPekerjaan, .data)
  kolom <- c('qty.ass1', 'qty.ass2', 'qty.ass3', 'qty.ass4', 
             'qty.cs', 'qty.hr', 'qty.lr', 'qty.or')
  hasil <- data.frame(t[1:8], row.names=kolom)
  return(t(hasil))
}

Pada function di atas, saya memakai xtabs untuk menghitung jumlah per jenis pekerjaan. Ini disebut juga sebagai crosstab. Setelah itu, saya menggunakan t() untuk menerjemahkan baris menjadi kolom. Function ini nantinya bisa saya lewatkan sebagai parameter untuk summarise() seperti pada perintah berikut ini:

> data <- ddply(servis, .(nomorKB), buatQtyJenisPekerjaan)
> head(data)
  nomorKB qty.ass1 qty.ass2 qty.ass3 qty.ass4 qty.cs qty.hr qty.lr qty.or
1  XXXXXX        0        0        0        0      0      0      0      1
2  YYYYYY        0        0        0        0      1      0      0      0
3  ZZZZZZ        0        0        0        0      1      0      2      0
4  AAAAAA        0        0        0        0      8      0      2      0
5  BBBBBB        0        0        0        0      2      1      3      1
6  CCCCCC        0        0        0        0      1      0      0      0

Proses pengolahan data ini terlihat sangat mudah, bukan? Yup! Hal ini karena R adalah sebuah bahasa pemograman sehingga ia dapat diprogram secara fleksibel.

Sekarang, saya akan menerapkan subset regression dengan menggunakan regsubsets() dari package leaps untuk mencari kombinasi variabel yang paling mempengaruhi qty.lr. Untuk itu, saya memberikan perintah seperti berikut ini:

> f <- regsubsets(qty.lr ~ qty.ass1 + qty.ass2 + qty.ass3 + qty.ass4 + qty.cs + qty.hr + qty.or, data)
> summary(f)
Subset selection object
Call: regsubsets.formula(qty.lr ~ qty.ass1 + qty.ass2 + qty.ass3 + 
    qty.ass4 + qty.cs + qty.hr + qty.or, data)
7 Variables  (and intercept)
         Forced in Forced out
qty.ass1     FALSE      FALSE
qty.ass2     FALSE      FALSE
qty.ass3     FALSE      FALSE
qty.ass4     FALSE      FALSE
qty.cs       FALSE      FALSE
qty.hr       FALSE      FALSE
qty.or       FALSE      FALSE
1 subsets of each size up to 7
Selection Algorithm: exhaustive
         qty.ass1 qty.ass2 qty.ass3 qty.ass4 qty.cs qty.hr qty.or
1  ( 1 ) " "      " "      " "      " "      " "    " "    "*"   
2  ( 1 ) " "      " "      " "      " "      "*"    " "    "*"   
3  ( 1 ) " "      " "      " "      " "      "*"    "*"    "*"   
4  ( 1 ) " "      " "      " "      "*"      "*"    "*"    "*"   
5  ( 1 ) " "      "*"      " "      "*"      "*"    "*"    "*"   
6  ( 1 ) " "      "*"      "*"      "*"      "*"    "*"    "*"   
7  ( 1 ) "*"      "*"      "*"      "*"      "*"    "*"    "*"   

Hasil di atas memperlihatkan bahwa bila harus menggunakan 1 predictor, maka predictor yang paling mempengaruhi mempengaruhi qty.lr adalah qty.or. Bila harus menggunakan 2 predictor, maka yang paling berpengaruh adalah qty.cs dan qty.or. Bila harus menggunakan 3 predictor, yang paling berpengaruh adalah qty.cs, qty.hr dan qty.or. Dan seterusnya.

Saya bisa melihat perubahan nilai R-squared bila memakai 1, 2, 3, hingga 7 predictor dengan menggunakan perintah berikut ini:

> s <- summary(f)
> s$rsq
[1] 0.2896593 0.3219981 0.3325110 0.3384162 0.3406387 0.3408268 0.3408714

Terlihat bahwa bila saya menggunakan seluruh 7 variabel yang ada, nilai R-squared meningkat hingga 34%. Akan tetapi, setelah menggunakan 5 variabel, tidak ada lagi peningkatan signifikan yang terjadi. Agar lebih jelas, saya dapat menggambarkan nilai R-squared dalam bentuk visual dengan perintah seperti berikut ini:

> plot(s$adjr2, type='l')
Plot adjusted R-square

Plot adjusted R-square

Dengan perintah yang sama, saya juga dapat menampilkan nilai statistik Cp dan Bayesin information criterion (BIC) yang juga dipakai untuk memilih subset yang optimal.

> plot(s$cp, type='l')
Plot Cp

Plot Cp

> plot(s$bic, type='l')
Plot Bic

Plot Bic

Berbeda dengan RSS, nilai Cp dan BIC yang baik adalah nilai yang semakin kecil.

Bila saya memakai 5 predictor, maka masing-masing koefisien untuk predictor tersebut adalah:

> coef(f,5)
(Intercept)    qty.ass2    qty.ass4      qty.cs      qty.hr      qty.or 
  0.4777711  -0.4552900  -0.5897922   0.1488950   0.7905409   0.4748867 

Hasil di atas memperlihatkan bahwa semakin sering perawatan after sales service (perawatan berkala pada sepeda motor baru) yang dilakukan pada sebuah sepeda motor, maka semakin kecil peluang pergantian sparepart di kemudian hari. Pengaruh antara servis rutin (CS), bongkar mesin (H/R) dan ganti oli (OR) tidak begitu besar terhadap pergantian suku cadang. Walaupun demikian, mereka tetap memicu pergantian suku cadang (karena nilai koefisiennya tidak 0). Mungkin saja ini terjadi karena pada saat melakukan perawatan berkala atau mengganti oli, pemilik kendaraan bermotor memperoleh informasi untuk mengganti suku cadang yang telah usang. Dari ketiga variabel tersebut, yang paling mempengaruhi pergantian suku cadang adalah bongkar mesin/pekerjaan berat (H/R). Ini cukup masuk akal karena H/R biasanya dikerjakan pada sepeda motor lama atau yang bermasalah seperti bekas tabrakan.

Tentu saja hasil analisa ini akan semakin akurat bila saya semakin familiar dengan domain ‘bengkel’ dan semakin banyak hasil observasi yang diperoleh dari berbagai bengkel🙂 Saya juga mungkin nanti perlu membuktikan bahwa jumlah servis rutin (CS) yang kecil memicu jumlah perawatan berat (H/R). Saya juga mungkin bisa menunjukkannya tidak dalam bentuk jumlah (qty) melainkan dalam ongkos. Saya juga bisa menambahkan beberapa variabel lain seperti tahun kendaraan agar hasil analisa semakin akurat.

Sebagai alternatif lain dari metode least squared dan subset selection, saya juga bisa mencoba memakai metode lasso dalam melakukan linear regression. Untuk itu saya dapat menggunakan glmnet() dari package glmnet. Sebagai contoh, saya memberikan perintah berikut ini:

> x <- model.matrix(qty.lr ~ qty.ass1 + qty.ass2 + qty.ass3 + qty.ass4 + qty.cs + qty.hr + qty.or, data)[,-1]
> y <- data$qty.lr
> library(glmnet)
Loading required package: Matrix
Loaded glmnet 1.9-8

> lasso <- glmnet(x, y, alpha=1)

Saya dapat melihat nilai lambda yang dipilih secara otomatis dengan menggunakan perintah seperti berikut ini:

> plot(lasso, 'lambda', label=TRUE)
Plot nilai lambda yang dipakai dan kaitannya pada koefisien predictor

Plot nilai lambda yang dipakai dan kaitannya pada koefisien predictor

Strategi pada teknik lasso adalah memilih nilai lambda yang tepat (sama seperti di teknik ridge selection) . Pada grafis di atas terlihat bahwa semakin besar nilai lambda, satu per satu dari nilai koefisien untuk predictor akan berkurang hingga menjadi 0 (nol). Lalu nilai lambda mana yang harus saya pakai? Saya harus melakukan cross-validating untuk mengetahuinya. Saya akan mulai dengan membagi data menjadi dua: untuk training dan untuk validation/test:

> train <- sample(1:nrow(x), nrow(x)/2)
> test <- (-train)
> y.test = y[test]
> lasso <- glmnet(x[train,], y[train], alpha=1)
> cv <- cv.glmnet(x[train,], y[train], alpha=1)
> cv$lambda.min
[1] 0.01750082
> cv$lambda.1se
[1] 0.3130284
> p <- predict(lasso, s=cv$lambda.min, newx=x[test,])
> mean((p-y.test)^2)
[1] 1.909886
> predict(lasso, type='coefficients', s=cv$lambda.min)
8 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept)  0.50056702
qty.ass1    -0.27511508
qty.ass2    -0.38088034
qty.ass3     .         
qty.ass4    -0.55677860
qty.cs       0.13011162
qty.hr       0.06781857
qty.or       0.51628946

Pada hasil di atas, nilai lambda yang menghasilkan MSE paling kecil adalah 0.01750082 (nilai ini ditampung di cv$lambda.min). Pada nilai lambda tersebut, 6 dari 7 variabel akan disertakan dimana seluruh variabel ass memiliki nilai negatif. Ini tidak berbeda dari hasil temuan sebelumnya yang menunjukkan bahwa semakin sering sebuah motor melakukan perawatan after sales service (ass) atau perawatan garansi, maka semakin sedikit penggantian sparepart yang akan terjadi di kemudian hari.

Berdasarkan hasil cross-validation, hasil ini memiliki nilai MSE yang lebih tinggi dibandingkan dengan cross-validation yang saya lakukan pada model yang memakai least squares (1.909886 > 1.863092). Sebagai patokan, nilai MSE yang semakin kecil menunjukkan tingkat kesalahan yang lebih sedikit.

Perihal Solid Snake
I'm nothing...

3 Responses to Memakai Linear Regression Di R

  1. Ping-balik: Memakai Local Regression Di R | The Solid Snake

  2. Ping-balik: Melakukan Unsupervised Learning Dengan R | The Solid Snake

  3. Imelda Wijaya mengatakan:

    mau tanya kalo mau kerja regresi logistik pake metode lasso di R , bgaimna yah?

Apa komentar Anda?

Please log in using one of these methods to post your comment:

Logo WordPress.com

You are commenting using your WordPress.com account. Logout / Ubah )

Gambar Twitter

You are commenting using your Twitter account. Logout / Ubah )

Foto Facebook

You are commenting using your Facebook account. Logout / Ubah )

Foto Google+

You are commenting using your Google+ account. Logout / Ubah )

Connecting to %s

%d blogger menyukai ini: