Package
Silahkan install jika belum ada
install.packages("glmnet")
install.packages("tidyverse")
install.packages("glmnetUtils")
install.packages("varbvs")
install.packages("leaps")library(tidyverse)## -- Attaching packages ---------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.0
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts ------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(glmnetUtils)
library(leaps)Data
Data yang digunakan pada praktikum kali ini adalah data Leukemia. Berikut adalah informasi singkat mengenai data
Expression levels recorded for 3,571 genes in 72 patients with leukemia (Golub et al, 1999). The binary outcome encodes the disease subtype: acute lymphobastic leukemia (ALL) or acute myeloid leukemia (AML).
data ini bisa diperoleh dalam package varbvs
Import Data ke R
data("leukemia",package = "varbvs")
leukemia = data.frame(y=leukemia$y,leukemia$x)Metode Forward Selection
Metode Forward selection ini dijalankan menggunakan fungsi regsubsets dari package leaps. Argumen nvmax disini adalah untuk menentukan banyaknya peubah/subset yang dimasukan ke model. Karena data leukemia memiliki 72 observasi, maka peubah yang dimasukan ke model tidak bisa lebih dari 71 peubah.
model_forward <- regsubsets(y~.,data=leukemia,
method = "forward",nvmax = 71)## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in =
## force.in, : 3500 linear dependencies foundKemudian, kita akan melihat subset mana yang terbaik. Penentuan ini dilakukan dengan menggunakan nilai BIC. Nilai BIC terkecil yang dianggap sebagai subset terbaik.
best_forward <- data.frame(subset=seq_along(summary(model_forward)$bic),
BIC = summary(model_forward)$bic)Selanjutnya dengan menggunakan plot dibawah ini, kita bisa melihat bahwa nilai yang paling kecil adalah pada subset paling kanan atau dengan kata lain subset ke 70 atau 71.
ggplot(best_forward,aes(subset,BIC))+
geom_line(size=1)+
geom_point(color="red",size=2)+
theme_bw()
Jika diperhatikan pada output dibawah ini, subset 71 nilai BICnya -Inf sehingga kita pilih yang subset 70.
best_forward[70:71,]## subset BIC
## 70 70 -5065.3
## 71 71 -InfSintaks dibawah ini digunakan untuk mengekstract nilai koefisien dari subset 70.
coef_forward = coef(model_forward,70)
length(coef_forward)## [1] 71Jika kita tampilkan, ada sebanyak 70 peubah yang terpilih dengan 1 intercept.
coef_forward ## (Intercept) X183 X218 X235 X287
## 3.472222e-01 -7.404151e-02 -6.090950e-08 7.576329e-03 -2.753443e-03
## X296 X298 X341 X348 X351
## -4.656128e-02 -2.685595e-04 -2.859934e-08 -2.095386e-04 -3.125436e-06
## X418 X446 X478 X598 X604
## 7.108548e-03 7.806248e-05 -1.666704e-04 -5.760329e-06 -4.508718e-04
## X625 X682 X691 X697 X710
## -1.916912e-05 -1.113632e-06 -2.082914e-02 -1.107586e-03 5.790618e-03
## X716 X784 X801 X849 X1053
## 7.966634e-03 2.169628e-02 -5.422339e-02 4.383456e-03 -3.178720e-02
## X1123 X1182 X1192 X1204 X1219
## 3.787290e-02 1.141740e-01 -1.234209e-02 1.736643e-11 -1.698153e-01
## X1267 X1335 X1493 X1551 X1560
## -1.877156e-02 1.451522e-02 2.187258e-04 -2.049213e-03 -5.764137e-10
## X1667 X1713 X1858 X1889 X1909
## -1.533883e-03 5.013101e-04 -1.541831e-05 -4.734743e-06 3.503166e-02
## X1946 X2089 X2102 X2107 X2155
## 7.517935e-02 6.261964e-04 3.353616e-02 1.917570e-02 1.403177e-04
## X2234 X2315 X2363 X2427 X2491
## -1.787946e-04 4.128251e-02 -3.593878e-03 3.409575e-02 2.144905e-02
## X2523 X2558 X2579 X2715 X2751
## -2.823069e-03 -4.467083e-02 6.571662e-05 4.532983e-05 -5.201828e-07
## X2761 X2811 X2861 X2888 X2915
## -1.124179e-03 1.097293e-02 -2.455356e-07 7.308873e-02 -1.769428e-03
## X2954 X3038 X3311 X3478 X3512
## 8.256019e-04 1.414855e-01 8.937411e-09 4.371139e-02 -2.100190e-03
## X3522 X3529 X3530 X3539 X3545
## 1.596075e-02 1.728537e-05 -2.879700e-02 2.584945e-02 -1.363613e-03
## X3557
## -4.631553e-02Metode LASSO
Metode kedua yang kita coba adalah metode LASSO. Metode ini bisa dijalankan menggunakan package glmnet. Package glmnetUtils merupakan package tambahan dari glmnet yang memungkinkan syntax glmnet bisa dinput menggunakan object data.frame.
Fungsi cv.glmnet digunakan untuk menjalankan menjalankan metode LASSO, asalkan nilai argumen alpha=1. Kemudian, dengan fungsi ini juga dimungkinkan memilih parameter penalty lambda terbaik dengan menggunakan cross-validation (validasi-silang) yang didasarkan pada kriteria tertentu, contohnya saja kriteria "deviance" dan "auc".Argumen nfolds digunakan untuk menentukan jumlah folds/lipatan pada cross-validation. Penjelasan cross-validation lebih lengkap bisa dilihat pada video pada link berikut https://www.youtube.com/watch?v=QRlSktU5Cn0&feature=youtu.be.
model_lasso_dev <- cv.glmnet(y~.,data=leukemia,alpha = 1,
type.measure="deviance",
family="binomial",
nfolds=5)
model_lasso_auc <- cv.glmnet(y~.,data=leukemia,alpha = 1,
type.measure="auc",
family="binomial",
nfolds=5)Setelah sintaks diatas dijalankan, selanjutnya dimungkinkan untuk membuat plot nilai lambda dalam bentuk log dengan nilai kriteria yang dipilih.
plot(model_lasso_dev)
plot(model_lasso_auc)
Informasi penting yang bisa dilihat pada plot diatas adalah banyaknya peubah terpilih yang ditulis diatas plot. Kemudian ada garis vertikal putus-putus yang menandakan nilai lambda optimum berdasarkan kriteria yang dipilih. Jika menggunakan AUC maka nilai lambda optimum dipilih berdasarkan nilai AUC terbesar. Sedangkan menggunakan deviance nilai lambda optimum dipilih berdasarkan nilai deviance terkecil. Selain itu, kriteria lambda optimum juga dipilih berdasarkan lembar selang terkecil dari garis pada setiap titik. Sehingga garis vertikal yang satunya menandakan hal tersebut.
Nilai koefisien dari variable terpilih bisa dikeluarkan dengan sintaks dibawah ini. Argumen s pada fungsi coef bisa diisi dengan "lambda.min"(peubah yang terpilih berdasarkan nilai lambda optimal) dan "lambda.1se"(peubah yang terpilih berdasarkan nilai lambda yang memilii selang terkecil).
coef_auc_min <- coef(model_lasso_auc,s = "lambda.min")
broom::tidy(coef_auc_min)## Warning: 'tidy.dgCMatrix' is deprecated.
## See help("Deprecated")## Warning: 'tidy.dgTMatrix' is deprecated.
## See help("Deprecated")## row column value
## 1 (Intercept) 1 -1.44398240
## 2 X158 1 0.03896292
## 3 X456 1 -0.44784427
## 4 X626 1 -0.11057166
## 5 X657 1 -0.02225145
## 6 X672 1 -0.58341501
## 7 X956 1 0.34402770
## 8 X979 1 1.10036406
## 9 X1182 1 0.05155648
## 10 X1219 1 -0.11762869
## 11 X1652 1 0.30744500
## 12 X1946 1 0.28345251
## 13 X2481 1 0.72333138
## 14 X2537 1 0.01597415
## 15 X2888 1 0.11039620
## 16 X3098 1 0.24329050
## 17 X3158 1 0.03852754coef_dev_min <- coef(model_lasso_dev,s = "lambda.min")
broom::tidy(coef_dev_min)## Warning: 'tidy.dgCMatrix' is deprecated.
## See help("Deprecated")## Warning: 'tidy.dgTMatrix' is deprecated.
## See help("Deprecated")## row column value
## 1 (Intercept) 1 -2.105351950
## 2 X158 1 0.052466245
## 3 X219 1 -0.132573105
## 4 X456 1 -0.614156024
## 5 X657 1 -0.161344071
## 6 X672 1 -1.000508499
## 7 X888 1 0.091913710
## 8 X918 1 0.088751241
## 9 X956 1 0.475777266
## 10 X979 1 1.602684717
## 11 X1219 1 -0.266342344
## 12 X1652 1 0.352491511
## 13 X1796 1 0.013808361
## 14 X1946 1 0.566741056
## 15 X2230 1 0.049321276
## 16 X2239 1 -0.014187057
## 17 X2387 1 -0.019321462
## 18 X2481 1 1.034589642
## 19 X2537 1 0.039528516
## 20 X2727 1 -0.030640421
## 21 X2888 1 0.239429232
## 22 X3098 1 0.404944214
## 23 X3158 1 0.051519134
## 24 X3292 1 0.006195577coef_auc_se <- coef(model_lasso_auc,s = "lambda.1se")
broom::tidy(coef_auc_se)## Warning: 'tidy.dgCMatrix' is deprecated.
## See help("Deprecated")## Warning: 'tidy.dgTMatrix' is deprecated.
## See help("Deprecated")## row column value
## 1 (Intercept) 1 -1.44398240
## 2 X158 1 0.03896292
## 3 X456 1 -0.44784427
## 4 X626 1 -0.11057166
## 5 X657 1 -0.02225145
## 6 X672 1 -0.58341501
## 7 X956 1 0.34402770
## 8 X979 1 1.10036406
## 9 X1182 1 0.05155648
## 10 X1219 1 -0.11762869
## 11 X1652 1 0.30744500
## 12 X1946 1 0.28345251
## 13 X2481 1 0.72333138
## 14 X2537 1 0.01597415
## 15 X2888 1 0.11039620
## 16 X3098 1 0.24329050
## 17 X3158 1 0.03852754coef_dev_se <- coef(model_lasso_dev,s = "lambda.1se")
broom::tidy(coef_dev_se)## Warning: 'tidy.dgCMatrix' is deprecated.
## See help("Deprecated")## Warning: 'tidy.dgTMatrix' is deprecated.
## See help("Deprecated")## row column value
## 1 (Intercept) 1 -1.691248517
## 2 X158 1 0.080294631
## 3 X219 1 -0.033533241
## 4 X456 1 -0.520924011
## 5 X626 1 -0.018440407
## 6 X657 1 -0.096450361
## 7 X672 1 -0.757169251
## 8 X888 1 0.012858732
## 9 X918 1 0.003853873
## 10 X956 1 0.388013410
## 11 X979 1 1.296957873
## 12 X1182 1 0.054770890
## 13 X1219 1 -0.181732121
## 14 X1652 1 0.348659359
## 15 X1946 1 0.380613406
## 16 X2481 1 0.850817980
## 17 X2537 1 0.036237227
## 18 X2888 1 0.166043488
## 19 X3098 1 0.308447879
## 20 X3158 1 0.051542961Metode Ridge
Metode terakhir yang kita coba adalah metode Ridge. Metode ini pada dasarnya bukan merupakan metode seleksi, namun metode untuk mengatasi multikolinearitas. Metode ini bisa membantu memperoleh seleksi peubah yang lebih baik jika digabungkan dengan metode LASSO. Penggabungan ini melahirkan metode baru yang bernama Elastic-Net yang akan dibahas pada praktikum berikutnya.
model_ridge_auc <- cv.glmnet(y~.,data=leukemia,alpha = 0,
type.measure="auc",
family="binomial",
nfolds=5)plot(model_ridge_auc)