Belajar Melakukan Penelitian Sederhana Dengan R

Masyarakat awam sering kali menganggap peneliti aneh. Hal ini mungkin karena pola pikir peneliti yang ‘khas’ dalam menyelesaikan masalah. Salah satu contohnya adalah orang awam sering menganggap peneliti suka mencari-cari kesalahan. Mengapa demikian? Banyak masyarakat awan memperoleh kebenaran secara emosional misalnya melalui keyakinan, insting, rayuan, empati, dan sebagainya. Masyarakat awam bisa saja menerima kebenaran tanpa syarat bila hal serupa diakui oleh mayoritas orang disekitarnya.

Sebaliknya, peneliti hanya akan menerima kebenaran setelah terbukti. Membuktikan sesuatu benar, apalagi benar 100% dalam setiap kondisi, adalah sesuatu yang sangat sulit. Akan tetapi, bila peneliti dapat menunjukkan kesalahan dari apa yang dianggap benar, maka apa yang tadinya dianggap benar bisa jadi mungkin adalah salah. Filosofer Karl Popper memiliki ide bahwa sesuatu yang ilmiah (scientific) adalah sesuatu yang bisa disangkal. Penyangkalan (falsifiability) berkaitan dengan pengujian (testability). Sebuah teori ilmiah adalah sesuatu yang bisa diuji (karena ia bisa salah) dan terbukti lolos dari pengujian dengan hasil yang diharapkan. Sebaliknya, sebuah pernyataan yang tidak dapat disangkal akan sulit untuk diuji. Ini hanya menunjukkan bahwa pernyataan tersebut tidak ilmiah (unscientific) tetapi bukan berarti ia salah atau tidak bermanfaat.

Filosofi Karl Popper menjadi landasan dari pengujian null hypothesis. Hipotesa yang hendak dibuktikan kebenarannya disebut sebagai alternative hypothesis (H1). Kebalikannya disebut sebagai null hypothesis (H0). Peneliti biasanya berfokus pada H0. Hasil penelitiannya adalah menolak atau gagal menolak H0! Kenapa peneliti menggunakan kalimat “gagal menolak H0” dan bukan “menerima H0”? Hal ini untuk menghindari kesalahan logika yang disebut sebagai argument from ignorance. Sebagai analogi: bila pengadilan gagal membuktikan kesalahan seorang terdakwa, bukan berarti terdakwa sudah pasti orang benar, bukan? Bisa saja terdakwa sudah menghapus bukti kesalahannya atau pencarian bukti tidak dilakukan di lokasi yang tepat 😀

Pada suatu hari, seorang mahasiswa pernah mengeluhkan harus membaca banyak baris kode program karena ia percaya bahwa kode program yang baik harus singkat. Ia berpendapat bahwa kode program yang jumlah barisnya banyak sudah pasti terlalu kompleks. Apakah benar demikian? Sebelum bisa menguji pendapat mahasiswa tersebut, saya perlu mendefinisikan lebih lanjut apa yang disebut sebagai kompleks. Saya akan menggunakan nilai cyclomatic complexity (https://en.wikipedia.org/wiki/Cyclomatic_complexity) sebagai acuan tingkat kerumitan kode program. Semakin banyak percabangan akibat perintah if, nilai cyclomatic complexity (M) akan semakin meningkat dan kode program akan semakin sulit ditelusuri. Dengan demikian:

H1: semakin besar jumlah baris program (LOC), maka semakin besar nilai cyclomatic complexity (M)

H0: jumlah baris program (LOC) tidak mempengaruhi nilai cyclomatic complexity (M)

Untuk membuktikan mahasiswa tersebut benar, saya perlu menguji dan menolak H0. Saya tidak mungkin melakukan pengujian pada seluruh kode program yang ada di dunia ini. Oleh sebab itu, saya hanya akan mengambil beberapa sampel pengujian, misalnya dari situs yang menyediakan statistik proyek open source seperti http://nemo.sonarqube.org. Saya kemudian mengelompokkan sampel yang ada menjadi 2 jenis: sampel dengan jumlah baris di bawah 10.000 dan sampel dengan jumlah baris program di atas 10.000. Gambar berikut ini memperlihatkan contoh sampel yang terkumpul:

Contoh sampel

Contoh sampel

Untuk menguji H0, saya dapat membandingkan mean dari sampel dengan jumlah baris di bawah 10.000 dan sampel dengan jumlah baris di atas 10.000. Bila nilai keduanya berbeda, saya bisa menolak H0 dengan aman.

Untuk membandingkan dua populasi yang tidak saling berhubungan, saya dapat menggunakan Student’s t-test (http://en.wikipedia.org/wiki/Student’s_t-test), tepatnya pengujian yang disebut sebagai independent sample t-test. Syarat untuk melakukan t-test adalah populasi yang dibandingkan harus mengikuti distribusi normal dan variance untuk kedua kelompok sampel yang diujikan harus sama. Bila variasi kedua kelompok sampel tidak sama, maka metode pengujian yang harus dipakai adalah Welch t-test. Function t.test() pada R secara default mengasumsikan kedua kelompok sampel memiliki variance yang berbeda sehingga akan menerapkan Welch’s t-test bila pengguna tidak memberikan argumen var.equal=TRUE.

Sebagai contoh, saya memberikan perintah seperti berikut ini di R:

> t.test(loc_besar$m, loc_kecil$m)

    Welch Two Sample t-test

data:  loc_besar$m and loc_kecil$m
t = 4.0958, df = 147.153, p-value = 6.92e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4729933 1.3550135
sample estimates:
mean of x mean of y 
 3.233871  2.319868 

Hasil paling penting yang perlu diperhatikan dari nilai di atas adalah p-value atau nilai p. Apa itu nilai p? Peneliti biasanya tidak pernah yakin 100%. Hampir semua peneliti akan skeptis mendengar pernyataan yang terlalu sempurna tanpa celah! Nilai p adalah peluang terjadinya Type I error dimana null hypothesis benar tetapi malah ditolak. Pada hasil di atas, peluang ini adalah 0.0000692. Hal ini berarti dari 2.500.000 kali percobaan dengan sampel berbeda, ada kemungkinan 173 kali hasil pengujian semata-mata terjadi karena kebetulan. Nilai p yang sering kali dijadikan patokan adalah 0.05. Karena nilai p untuk hasil pengujian di bawah 0.05, maka ada peluang besar bahwa mean kedua populasi tersebut memang berbeda.

Hasil pengujian t-test akan akurat bila populasi yang diperbandingkan memiliki distribusi normal. Saya akan mencoba melakukan pengujian dengan menggunakan Mann-Whitney U test yang tidak mengasumsikan distribusi normal, dengan memberikan perintah seperti berikut ini:

> wilcox.test(loc_besar$m, loc_kecil$m)

    Wilcoxon rank sum test with continuity correction

data:  loc_besar$m and loc_kecil$m
W = 12700.5, p-value = 3.544e-07
alternative hypothesis: true location shift is not equal to 0

Nilai p yang dibawah 0.05 kembali menunjukkan bahwa mean populasi kemungkinan besar tidak sama. Berdasarkan hasil ini, saya bisa menolak H0 (jumlah LOC tidak mempengaruhi M). Dengan demikian, bukti yang kuat semakin mengarah ke H1. Saya dapat menambahkan nilai alternative pada perintah sebelumnya untuk menguji arah (one sided test):

> wilcox.test(loc_besar$m, loc_kecil$m, alternative='greater')

    Wilcoxon rank sum test with continuity correction

data:  loc_besar$m and loc_kecil$m
W = 12700.5, p-value = 1.772e-07
alternative hypothesis: true location shift is greater than 0

alternative='greater' akan menguji apakah nilai rata-rata populasi di sebelah kiri (loc_besar$m) lebih besar daripada rata-rata populasi di sebelah kanan (loc_kecil$m). Nilai p yang kecil menunjukkan bahwa ada peluang hasil ini bukan kebetulan.

Salah satu pertanyaan yang mungkin muncul adalah mengapa tidak langsung membandingkan rata-rata (mean) sampel? Hasil akan kurang akurat karena rata-rata sampel tidak dapat dipakai untuk mendeskripsikan rata-rata populasi. Pengujian seperti t-test dan Mann-Whitney U test akan bekerja pada rata-rata (mean) populasi. Dengan demikian, walaupun saya hanya bekerja pada sampel yang jumlahnya terbatas (kode program pada beberapa proyek open-source), hasil yang diperoleh dapat diaplikasikan pada populasi yang saya teliti (yaitu seluruh kode program pada proyek open-source).

Untuk menunjukkan hubungan antar jumlah baris kode program (LOC) dan cyclomatic complexity (M), saya dapat membuat scatterplot untuk seluruh sampel yang ada seperti pada gambar berikut ini:

> require(ggplot2)

> ggplot(d, aes(x=loc, y=m)) + geom_point() + geom_smooth(method=lm)
Scatterplot yang memperlihatkan relasi diantara 2 variabel

Scatterplot yang memperlihatkan relasi diantara 2 variabel

Agar lebih yakin lagi, saya dapat menguji korelasi antar variabel dengan mencari nilai Pearson’s Product-Moment Correlation Coefficient (nilai r). Sebagai contoh, saya memberikan perintah berikut ini di R:

> cor.test(~ m + loc, d)

    Pearson's product-moment correlation

data:  m and loc
t = 3.7321, df = 270, p-value = 0.0002315
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1053281 0.3316868
sample estimates:
      cor 
0.2214889 

Karena nilai p dibawah 0.05, maka hasil pengujian menunjukkan kemungkinan adanya korelasi yang lemah antara jumlah baris dan cyclomatic complexity dengan nilai r berupa 0.2214889. Nilai r yang semakin mendekati 0 menunjukkan semakin tidak ada korelasi, dan nilai r yang semakin mendekati 1 atau -1 menunjukkan korelasi yang semakin erat.

Bila nilai r dipangkat dua, maka saya akan memperoleh nilai coefficient of determination yang berupa 0.04905735 atau 4.9%. Nilai ini menunjukkan bahwa jumlah baris memberikan dampak 4.9% perubahan pada kerumitan. Tentu saja masih ada kemungkinan variabel tersembunyi (yang ikut meningkat seiring peningkatan jumlah baris) yang belum saya temukan yang sesungguhnya memiliki dampak besar.

Setelah percobaan sederhana ini, saya kemudian bisa menyimpulkan bahwa pernyataan si mahasiswa memiliki landasan bukti yang bisa diterima. Sesungguhnya sudah ada beberapa paper yang membahas masalah ini secara lebih ilmiah, misalnya paper Cyclomatic Complexity and Lines of Code: Empirical Evidence of a Stable Linear Relationship yang diterbitkan di Journal of Software Engineering and Applications, Vol. 2 No. 3, 2009 dan paper Empirical Analysis of the Relationship between CC and SLOC in a Large Corpus of Java Methods yang diterbitkan di Software Maintenance and Evolution (ICSME), 2014 IEEE International Conference on.

Iklan

Melakukan Unsupervised Learning Dengan R

Pada artikel sebelumnya, Memakai Local Regression Di R, saya menerapkan sebuah bentuk regression yang tidak dapat digunakan untuk melakukan prediksi. Kali ini, saya juga akan mencoba motode yang tidak dapat digunakan untuk memprediksi yang masuk dalam kategori unsupervised learning. Mengapa unsupervised learning tidak dapat dipakai untuk prediksi? Hal ini karena ia hanya melibatkan predictor (nilai X) dan nilai n (jumlah) untuk nilai predictor yang sama. Lalu apa manfaatnya bila tidak dapat digunakan untuk prediksi? Metode seperti ini dapat dipakai untuk melakukan pengelompokan variabel secara lebih rinci. Sebagai contoh, pada bidang pemasaran (marketing), berbagai metode unsupervised learning sering digunakan untuk mengelompokkan konsumen berdasarkan perilaku belanja mereka.

Saya akan mulai dengan membaca file CSV yang saya peroleh berdasarkan angket undian pada periode tertentu:

> pelanggan <- read.csv("pelanggan.csv")
> tail(pelanggan)
     tipeMotor   ponsel region
192 SUPRAX 125 KARTU AS   KUBU
193    LEGENDA  SIMPATI   KUBU
194 SUPRAX 125  SIMPATI   KUBU
195      GRAND KARTU AS   JAWI
196  SUPRA FIT KARTU AS SERDAM
197  SUPRA FIT KARTU AS SERDAM

Saya akan mulai dengan menerapkan metode coresspondence analysis (CA). Bila seandainya data yang saya miliki adalah data kuantitatif, bukan data kualitatif, maka saya harus menggunakan metode principal component analysis (PCA). Untuk menerapkan CA, saya dapat memakai function CA() dari package FactoMineR. Sebelum mulai, saya perlu terlebih dahulu mengubah data yang saya ke dalam bentuk crosstab. Bila pada artikel Memakai Linear Regression Di R, saya menggunakan cara manual dengan membuat function sendiri, kali ini saya akan menggunakan cara yang lebih mudah dengan memanggil function dcast() dari package reshape2 seperti berikut ini:

> require(reshape2)
Loading required package: reshape2
> data <- dcast(pelanggan, tipeMotor ~ region)
Using region as value column: use value.var to override.
Aggregation function missing: defaulting to length
> data
    tipeMotor AYANI GAMA JAWI JERUJU KOBAR KUBU MERDEKA PURNAMA SERDAM SIANTAN TAPUR
1      ASTREA     0    0    0      0     0    0       0       0      0       1     0
2        BEAT     0    2    2      0     0    2       1       0      2       1     5
3       BLADE     1    0    1      0     0    0       0       2      0       0     2
4         CS1     0    0    0      0     0    0       0       1      0       0     0
5       GRAND     0    2    2      1     0    1       1       0      0       1     1
6    KHARISMA     0    1    2      0     0    1       2       0      0       1     2
7     LEGENDA     0    1    0      0     0    1       0       1      0       0     0
8    MEGA PRO     0    1    0      0     0    1       0       0      0       0     0
9       PRIMA     1    1    1      0     0    2       0       0      0       0     1
10       REVO     0    3    0      0     1    3       1       1      1       0     2
11   REVO 110     0    0    1      0     1    2       2       0      1       1     1
12   REVO FIT     0    0    0      0     0    0       0       0      0       0     3
13     SCOOPY     0    1    0      0     0    0       0       0      0       0     0
14      SONIC     0    0    0      0     0    1       0       0      0       0     0
15       STAR     0    0    0      0     0    0       0       0      0       0     1
16      SUPRA     0    1    0      0     0    2       0       0      1       1     0
17  SUPRA FIT     0    1    2      2     2    4       1       1      4       2     2
18    SUPRA X     3    3    2      1     2    3       2       4      2       2     4
19 SUPRAX 125     3    8    8      2     3   11       5       0      3       7     6
20      TIGER     1    0    1      0     0    0       0       0      0       0     0
21      VARIO     0    1    1      0     0    2       0       1      0       1     0

Setelah memperoleh crosstab, langkah berikutnya adalah mengubah kolom pertama (tipe motor) menjadi nama baris. Untuk itu, saya memberikan perintah berikut ini:

> data <- data.frame(data[,-1], row.names=data[,1])
> data
           AYANI GAMA JAWI JERUJU KOBAR KUBU MERDEKA PURNAMA SERDAM SIANTAN TAPUR
ASTREA         0    0    0      0     0    0       0       0      0       1     0
BEAT           0    2    2      0     0    2       1       0      2       1     5
BLADE          1    0    1      0     0    0       0       2      0       0     2
CS1            0    0    0      0     0    0       0       1      0       0     0
GRAND          0    2    2      1     0    1       1       0      0       1     1
KHARISMA       0    1    2      0     0    1       2       0      0       1     2
LEGENDA        0    1    0      0     0    1       0       1      0       0     0
MEGA PRO       0    1    0      0     0    1       0       0      0       0     0
PRIMA          1    1    1      0     0    2       0       0      0       0     1
REVO           0    3    0      0     1    3       1       1      1       0     2
REVO 110       0    0    1      0     1    2       2       0      1       1     1
REVO FIT       0    0    0      0     0    0       0       0      0       0     3
SCOOPY         0    1    0      0     0    0       0       0      0       0     0
SONIC          0    0    0      0     0    1       0       0      0       0     0
STAR           0    0    0      0     0    0       0       0      0       0     1
SUPRA          0    1    0      0     0    2       0       0      1       1     0
SUPRA FIT      0    1    2      2     2    4       1       1      4       2     2
SUPRA X        3    3    2      1     2    3       2       4      2       2     4
SUPRAX 125     3    8    8      2     3   11       5       0      3       7     6
TIGER          1    0    1      0     0    0       0       0      0       0     0
VARIO          0    1    1      0     0    2       0       1      0       1     0

Sekarang, saya bisa memanggil perintah CA(). Hasilnya adalah sebuah grafis yang terlihat seperti pada gambar berikut ini:

> require(FactoMineR)
Loading required package: FactoMineR
> CA(data)
Hasil dari function CA()

Hasil dari function CA()

Pada grafis di atas, semakin dekat dua buah titik menunjukkan bahwa mereka semakin memiliki hubungan. Sebagai contoh, terlihat bahwa pelanggan yang memakai sepeda motor SupraX 125 memiliki hubungan yang lebih erat dengan wilayah Siantan, Kubu, dan Serdam.

Bagaimana bila saya ingin mencari hubungan untuk lebih dari dua variabel? Saya bisa menggunakan metode multiple correspondence analysis (MCA). Untuk itu, saya bisa memberikan perintah seperti berikut ini:

> pelanggan.mca <- MCA(pelanggan)
Hasil dari function MCA()

Hasil dari function MCA()

Pada grafis di atas, terlihat bahwa variabel ponsel dan region memiliki hubungan yang lebih erat bila dibandingkan dengan hubungan mereka dengan variabel tipeMotor.

Untuk melihat versi grafis yang lebih detail, saya bisa memberikan perintah berikut ini:

> plot(pelanggan.mca)
Plot MCA yang menyertakan nilai individual.

Plot MCA yang menyertakan nilai individual.

Pada grafis di atas, titik biru mewakili setiap baris data. Terlihat bahwa tidak ada pengelompokan data pada kombinasi variabel tertentu.

Operasi lain yang sering dilakukan pada unsupervised learning adalah clustering. Sebagai contoh, pada R, saya dapat melakukan hierarchial clustering dengan menggunakan perintah hclust() seperti berikut ini:

> pelanggan.hclust <- hclust(dist(data))
> plot(pelanggan.hclust)
Dendrogram yang dihasilkan oleh hclust()

Dendrogram yang dihasilkan oleh hclust()

Grafis yang dihasilkan disebut dengan dendrogram. Berbeda dengan K-Means Clustering, saya tidak perlu menentukan jumlah pengelompokan (cluster) bila memakai hierarchial clustering. Sebagai gantinya, saya dapat memilih jumlah cluster berdasarkan dendrogram. Sebagai contoh, pada nilai height=10, saya akan memperoleh 3 kelompok untuk seluruh jenis kendaraan yang ada, yaitu kelompok Supra X 125, kelompok di tengah yang terdiri atas berbagai kendaraan, dan kelompok di kanan (yang terdiri atas Supra X, Suprat Fit, Beat, dan Revo).

Memakai Local Regression Di R

Pada artikel sebelumnya, Memakai Linear Regression Di R, saya mencoba menerapkan linear regression di R. Saya menyebutnya klasik karena teknik tersebut sudah ada sejak lama dan dapat dihitung secara manual. Seiring dengan penggunaan komputer pada bidang statistika, muncul pula berbagai teknik baru yang menawarkan akurasi yang lebih baik. Salah satunya adalah local regression yang memakai algoritma k-Nearest-Neighbor sehingga memperhitungkan titik disekitar saat menghitung posisi sebuah titik. Karena perhitungan ini harus dilakukan pada setiap titik, satu-per-satu sepanjang grafis, local regression boleh dibilang mustahil diproses secara manual tanpa melibatkan komputer.

Sebagai contoh, saya akan memakai data pendapatan dan pengeluaran yang saya baca dengan perintah seperti berikut ini:

> data <- read.csv('finance.csv')
> names(data)
[1] "tgl"    "pendapatan"  "pengeluaran"
> ggplot(data, aes(tgl)) +
+ geom_line(aes(tgl, pendapatan, color='Pendapatan'), size=1) +
+ geom_line(aes(tgl, pengeluaran, color='Pengeluaran'), alpha=0.3) +
+ theme(legend.title=element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank()) +
+ scale_color_manual(values=c('Pengeluaran'='red', 'Pendapatan'='blue'))
Grafis pendapatan per tanggal

Grafis pendapatan per tanggal

Saya akan mencoba mencari persamaan yang paling mendekati untuk memetakan data pendapatan berdasarkan tanggal. Bila memakai simple linear regression, saya akan memperoleh hasil seperti berikut ini:

> grafis <- ggplot(data, aes(tgl,pendapatan)) + theme(legend.title=element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank()) 
> grafis + geom_line() + geom_smooth(method='lm')
Memakai linear regression pada penjualan per tanggal.

Memakai linear regression pada penjualan per tanggal.

Sekarang, bandingkan grafis di atas dengan versi polynomial regression yang saya buat dengan perintah berikut ini:

> grafis + geom_line() + geom_smooth(method='lm', formula=y~poly(x,5))
Memakai polynomial regression pada penjualan per tanggal.

Memakai polynomial regression pada penjualan per tanggal.

Versi yang memakai polynomial regression terlihat lebih baik, bukan? Agar lebih yakin, saya akan berusaha membuktikannya secara numeris. Walaupun polynomial regression tidak lagi linear, ia masih berbasis model linear sehingga teknik pada model linear masih dapat saya gunakan. Oleh sebab itu, saya melakukan perbandingan dengan memberikan perintah berikut ini:

> summary(lm(pendapatan ~ tgl, data=data))

Call:
lm(formula = pendapatan ~ tgl, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-3217705 -1000690   -85050  1097801  3938758 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) 14408908.7  8270127.5   1.742    0.084 .
tgl             -652.6      520.4  -1.254    0.212  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1452000 on 122 degrees of freedom
Multiple R-squared:  0.01273,   Adjusted R-squared:  0.004635 
F-statistic: 1.573 on 1 and 122 DF,  p-value: 0.2122

> summary(lm(pendapatan ~ poly(tgl,5), data=data))

Call:
lm(formula = pendapatan ~ poly(tgl, 5), data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-2950460  -740133  -158111   612332  2881061 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    4038452      99700  40.506  < 2e-16 ***
poly(tgl, 5)1 -1820799    1110207  -1.640    0.104    
poly(tgl, 5)2 -9361247    1110207  -8.432 9.73e-14 ***
poly(tgl, 5)3  -682886    1110207  -0.615    0.540    
poly(tgl, 5)4  1404591    1110207   1.265    0.208    
poly(tgl, 5)5  4652563    1110207   4.191 5.40e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1110000 on 118 degrees of freedom
Multiple R-squared:  0.4416,    Adjusted R-squared:  0.418 
F-statistic: 18.67 on 5 and 118 DF,  p-value: 1.247e-13

Pada hasil di atas, terlihat bahwa versi yang memakai simple linear regression memiliki nilai R-squared hanya sebesar 1,2%. Sementara itu, versi yang memakai polynomial regression memiliki nilai R-squared sebesar 41,8%. Ini adalah perbedaan yang cukup besar!

Tentu saja saya tidak bisa begitu saja memprediksi pendapatan berdasarkan tanggal (terlihat dari p-value). Hal ini ditambah lagi dengan penggunaan persamaan polynomial yang bisa membisa memberikan hasil yang tak terduga. Misalnya, prediksi pendapatan untuk bulan Desember adalah:

> cbind(predict(f,data.frame(tgl=c(seq(as.Date('2014/12/01'), as.Date('2014/12/10'), 'day')))))
       [,1]
1  10171185
2  10307194
3  10444780
4  10583953
5  10724726
6  10867108
7  11011112
8  11156748
9  11304028
10 11452963

Wow, pendapatan meningkat secara drastis! Tapi saya tidak bisa bahagia karena ini kemungkinan besar disebabkan oleh kesalahan model (atau memang efek akhir tahun, ya? 🙂 )

Sekarang, saya akan menggunakan teknik local regression pada data yang sama. Grafis berikut ini memperlihatkan geom_smooth() yang menggunakan method loess:

> grafis + geom_line() + geom_smooth(method='loess')
Memakai local regression dengan nilai alpha = 0.75

Memakai local regression dengan nilai alpha = 0.75

Sebagai informasi, method geom_smooth() dari package ggplot2 secara otomatis akan memakai method loess bila n < 1000. Berdasarkan Wikipedia, LOESS adalah semacam singkatan dari “LOcal regrESSIon”. Grafis yang dihasilkan terlihat lebih lembut bila dibandingkan dengan versi yang memakai polynomial regression. Pada local regression, terdapat nilai alpha yang disebut sebagai smoothing parameter. Grafis di atas memakai dengan nilai default alpha berupa 0.75. Saya dapat menggunakan nilai alpha yang berbeda dengan perintah seperti berikut ini:

> grafis + geom_point() + geom_smooth(method='loess', span=0.25)
Memakai local regression dengan nilai alpha = 0.25

Memakai local regression dengan nilai alpha = 0.25

Saya dapat melakukan fitting yang memakai local regression dengan perintah seperti berikut ini:

> f <- loess(pendapatan ~ as.numeric(tgl), data, span=0.25)

Pada perintah di atas, saya perlu menggunakan as.numeric() untuk menerjemahkan Date menjadi angka agar dapat dihitung. Btw, pada saat memakai geom_smooth() dalam menampilkan grafis, saya tidak perlu melakukan konversi karena hal ini telah dilakukan secara otomatis dimana tanggal diterjemahkan menjadi angka pada sumbu X.

Local regression tidak mengenal nilai seperti R-squared karena ia lebih berguna sebagai model untuk explorasi ketimbang sebagai model prediksi. Hal ini cukup masuk akal karena local regression sangat bergantung pada titik-titik yang sudah disekitarnya. Cara yang paling efektif untuk melihat keakuratan local regression adalah dengan mengamati grafis yang dihasilkan olehnya secara seksama. Walaupun ada function predict() untuk local regression, saya hanya bisa menggunakannya untuk nilai tgl yang sudah ada dan bukan untuk tgl di-luar sampel:

> x.prediksi <- as.numeric(seq(as.Date('2014/12/01'), as.Date('2014/12/10'), 'day'))
> cbind(predict(f, data.frame(tgl=x.prediksi)))
   [,1]
1    NA
2    NA
3    NA
4    NA
5    NA
6    NA
7    NA
8    NA
9    NA
10   NA
> x.prediksi <- as.numeric(seq(as.Date('2013/12/01'), as.Date('2013/12/10'), 'day'))
> cbind(predict(f, data.frame(tgl=x.prediksi)))
      [,1]
1  5617313
2  5627744
3  5637711
4  5647098
5  5655790
6  5663674
7  5670634
8  5676556
9  5681326
10 5684828

Dengan demikian, local regression lebih tepat dipakai untuk mempelajari hubungan antar-variabel dimana hubungan tersebut sangat kompleks dan tidak dapat dijelaskan oleh model yang sudah ada.

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.

Memakai Data Kualitatif Di R

Data kualitatif adalah data yang tidak dalam bentuk angka. Pada R, data kualitatif diwakili oleh tipe data factor. Berbeda dengan tipe data character yang dapat di-isi apa saja, tipe data factor memiliki pilihan atau kategori (misalnya kecil, sedang, dan besar). Programmer akan lebih mengenal factor sebagai enumeration yang identik dengan keyword enum di bahasa pemograman umum. Pada artikel ini, saya akan mencoba melakukan pembelajaran statistik yang melibatkan tipe data kualitatif.

Saya akan mulai dengan membaca data CSV dengan menggunakan perintah berikut ini:

> servis <- read.csv("servis.csv", header=TRUE)
> head(servis)
  nomorKB   tipeMotor jenisPekerjaan
1  XXXXXX Supra X 125            L/R
2  XXXXXX     Supra X             CS
3  XXXXXX       Grand            L/R
4  XXXXXX     Supra X             CS
5  XXXXXX     Megapro            L/R
6  XXXXXX        Beat             CS
> sapply(servis, class)
       nomorKB      tipeMotor jenisPekerjaan 
      "factor"       "factor"       "factor" 

Terlihat bahwa R membaca data yang ada sebagai factor. Untuk melihat kategori apa saja yang diperbolehkan untuk sebuah factor, saya dapat menggunakan perintah levels() seperti pada perintah berikut ini:

> levels(servis$tipeMotor)
 [1] "Absolute Revo 110" "Astrea"            "Beat"              "Blade"             "CS1"              
 [6] "GL Series"         "Grand"             "Karisma"           "Kirana"            "Lain-lain"        
[11] "Legenda"           "Megapro"           "Prima"             "Revo"              "Revo Fit"         
[16] "Scoopy"            "Spacy"             "Supra"             "Supra Fit"         "Supra X"          
[21] "Supra X 125"       "Tiger"             "Vario"             "Win"              
> levels(servis$jenisPekerjaan)
[1] "ASS1" "ASS2" "ASS3" "ASS4" "CS"   "H/R"  "L/R"  "OR+" 

Saya dapat melihat distribusi jumlah data untuk masing-masing factor dengan menggunakan function table() seperti pada perintah berikut ini:

> table(servis$jenisPekerjaan)

ASS1 ASS2 ASS3 ASS4   CS  H/R  L/R  OR+ 
  49   48   46   49 1861   55 1439  968 
> cbind(table(servis$tipeMotor))
                  [,1]
Absolute Revo 110  252
Astrea              20
Beat               378
Blade              121
CS1                 12
GL Series           37
Grand              228
Karisma            185
Kirana               2
Lain-lain          341
Legenda             51
Megapro             50
Prima              120
Revo               175
Revo Fit            52
Scoopy              35
Spacy               53
Supra               57
Supra Fit          325
Supra X            608
Supra X 125       1144
Tiger               55
Vario              192
Win                 22

Agar lebih mudah dipahami, saya juga bisa menggambarkannya dalam bentuk histogram dengan menggunakan gglot() seperti pada perintah berikut ini:

> require(ggplot2)
> ggplot(servis, aes(x=servis$jenisPekerjaan)) + geom_histogram()
Histogram untuk jenis pekerjaan

Histogram untuk jenis pekerjaan

> ggplot(servis, aes(x=servis$tipeMotor)) + geom_histogram()
Histogram untuk tipe motor

Histogram untuk tipe motor

Lalu bagaimana cara menggunakan data kualitatif? Saya tidak bisa dengan mudah memetakannya pada persamaan linear karena mereka bukanlah angka yang dapat dihitung. Salah satu trik untuk memakai data kualitatif di persamaan linear adalah dengan melakukan pengkodean dalam bentuk angka. Misalnya, pada jenisPekerjaan, saya bisa memberi nilai 0 untuk ASS1, 1 untuk ASS2, 2 untuk ASS3, dan seterusnya. Tapi solusi yang lebih baik adalah dengan menggunakan metode logistic regression bila variabel respons terdiri atas 2 kategori dan metode linear discrimant analysis bila variabel respon terdiri atas banyak kategori.

Sebagai contoh, berdasarkan riwayat yang ada, saya ingin memprediksi jenis pekerjaan yang akan dilakukan pada sebuah motor yang hendak di-servis berdasarkan jenis motor tersebut. Karena jenis pekerjaan terdiri atas banyak kategori, maka saya akan memakai metode linear discriminant analysis. Untuk itu, saya dapat memanggil function lda() milik package MASS seperti berikut ini:

> require(MASS)
Loading required package: MASS
> f <- lda(jenisPekerjaan ~ tipeMotor, data=servis)
> summary(f)

Sekarang, saya dapat mencoba melakukan prediksi, misalnya saya ingin mencari tahu jenis pekerjaan apa yang mungkin akan dilakukan bila tipe motor adalah Supra X. Untuk itu, saya memberikan perintah berikut ini:

> predict(f, data.frame(tipeMotor=c('Supra X')))
$class
[1] L/R
Levels: ASS1 ASS2 ASS3 ASS4 CS H/R L/R OR+

$posterior
         ASS1        ASS2        ASS3        ASS4        CS        H/R       L/R       OR+
1 0.001097874 0.001214338 0.001023112 0.001036094 0.3659366 0.01078961 0.3676512 0.2512511

$x
        LD1        LD2        LD3        LD4        LD5      LD6        LD7
1 0.7108102 -0.1152886 -0.7093833 -0.7134284 -0.2333662 -0.20007 0.09106148

Linear discriminant analysis bekerja berdasarkan teori Bayess dan menghasilkan posterior probability. Nilai prior probability diperoleh berdasarkan data training. Pada hasil di atas, terlihat bahwa sepeda motor Supra X yang selanjutnya datang mungkin akan membutuhkan perawatan L/R (penggantian suku cadang) karena nilai posterior probability L/R adalah yang paling besar (0.3659366).

Bagaimana dengan sepeda motor jenis lain? Saya kembali memberikan sebuah perintah seperti berikut ini:

> predict(f, data.frame(tipeMotor=c('Win')))
$class
[1] H/R
Levels: ASS1 ASS2 ASS3 ASS4 CS H/R L/R OR+

$posterior
          ASS1         ASS2         ASS3         ASS4        CS       H/R        L/R        OR+
1 0.0001578041 0.0001753791 0.0001470381 0.0001497294 0.0539218 0.8652549 0.04636545 0.03382785

$x
        LD1       LD2       LD3      LD4       LD5       LD6      LD7
1 0.8467111 0.3420979 -3.655081 6.950073 -1.659424 0.3081269 1.477134

Pada hasil di atas, terlihat bahwa sepeda motor Win yang selanjutnya datang mungkin akan membutuhkan perawatan H/R (bongkar berat) mengingat nilai posterior probability H/R adalah yang paling besar (0.8652549). Ini cukup masuk akal karena sepeda motor Win adalah jenis sepeda motor lama yang tidak diproduksi lagi.

Tentu saja hasil di atas tidak begitu saja bisa saya percaya. Keakuratannya tergantung pada data training yang dipakai karena linear discriminant analysis bekerja dengan baik pada data yang memiliki distribusi normal (Gaussian). Untuk menguji keakuratannya, saya perlu melakukan cross validation (misalnya dengan menambahkan CV=TRUE pada lda() untuk melakukan Leave-one-out Cross Validation/LOOCV). Selain itu, tentu saja prediksi jenis pekerjaan tidak hanya dipengaruhi oleh variabel tipe motor saja melainkan bisa saja ada variabel lain yang belum saya sertakan seperti kondisi motor, tahun beli, dan sebagainya.

Memakai Bahasa Pemograman R Untuk Analisa Statistik

Penerapan statistik dapat dijumpai dengan mudah di berbagai bidang ilmu seperti ekonomi, biologi (biostatistik), geologi, psikologi, ilmu komputer, dan sebagainya. Semakin banyak pula karya ilmiah mahasiswa yang kini dilengkapi dengan presentasi hasil analisa statistik. Salah satu tool yang sering dipakai untuk pembelajaran statistik adalah IBM SPSS (Statistical Package for the Social Science). Walaupun user friendly, tool ini tidak gratis dan tidak open source! Salah satu alternatif yang menarik adalah bahasa pemograman R. R adalah sebuah bahasa pemograman yang dikembangkan khusus untuk hal-hal yang berkaitan dengan statistik. Saya merasa R sangat tepat dipakai oleh mahasiswa ilmu komputer karena mereka sudah terbiasa dengan hal-hal teknis seperti kode program.

Bila IBM SPSS memiliki tampilan GUI, R dipakai dengan memberikan perintah atau kode program. Sebenarnya IBM SPSS juga memiliki kemampuan scripting, tetapi terbatas dan tidak secanggih R. Apa kelebihan scripting dibandingkan GUI? Hal yang sering saya jumpai pada saat ‘belajar kelompok’ IBM SPSS adalah mahasiswa yang satu meminta rekannya mendemonstrasikan langkah demi langkah seluruh menu yang di-klik dan textbox yang di-isi untuk memperoleh sebuah hasil. Sebagai informasi, untuk memperoleh sebuah hasil yang berguna, mahasiswa harus melakukan banyak eksperimen dan membuat berbagai grafis untuk di-analisa (semua langkah ini cukup untuk dijadikan satu soal UAS :). Langkah-langkah ini terkadang sulit di-reproduksi ulang dengan mudah dan ada kemungkinan langkah yang terlupakan oleh mahasiswa. Bila memakai R, para mahasiswa hanya perlu men-‘copy paste’ kode program untuk memperoleh hasil yang sama. Mereka juga bisa menjalankan kembali kode program R yang sama pada data yang berbeda.

R dapat di-download di http://www.r-project.org/. Selain itu, saya juga dapat men-download versi Revolution R Open di http://revolutionanalytics.com/download-r yang dilengkapi dengan tambahan Intel Math Kernel Library (MKL) yang memberikan peningkatan kinerja terutama di Windows dan juga prosesor non-Intel. Selanjutnya, saya butuh sebuah IDE agar lebih nyaman bekerja dengan R. Untuk itu, saya akan memakai R Studio yang dapat di-download secara gratis di http://www.rstudio.com/products/rstudio/download.

Tampilan awal R Studio akan terlihat seperti pada gambar berikut ini:

Tampilan R Studio

Tampilan R Studio

Salah satu kendala awal dalam beralih ke R yang sering jumpai adalah bagaimana cara memasukkan data? IBM SPSS memiliki worksheet mirip Excel dimana pengguna bisa mengisi nilai untuk setiap variabel. R tidak memiliki kemampuan serupa! Walaupun ada beberapa package tambahan untuk membaca data dari file Excel atau database, cara paling mudah untuk ‘mengisi’ data adalah menggunakan file CSV. Hampir semua aplikasi modern sudah bisa men-export data dalam bentuk CSV. Sebagai contoh, saya bisa men-export isi file Excel menjadi file CSV dengan memilih CSV (Command delimited) di bagian Save as type pada saat menyimpan file. Bila ingin menyimpan hasil query MySQL Workbench menjadi file CSV, saya bisa men-klik tombol Export dan memilih CSV di Save as type.

Setelah itu, saya bisa membaca file CSV dengan memberikan perintah seperti berikut ini:

> penjualan <- read.csv("~/Desktop/penjualan_produk1.csv", header=TRUE)

Untuk melihat isi dari variabel penjualan yang telah dibaca di R Studio, saya dapat men-klik tombol tabel di samping kanan nama variabel seperti yang terlihat pada gambar berikut ini:

Melihat Isi DataFrame

Melihat Isi DataFrame

Selain itu, saya juga dapat melihat preview data dengan menggunakan function head() atau tail() seperti:

> head(penjualan_produk1)
     tanggal qty
1 2014-09-02  20
2 2014-09-03  30
3 2014-09-04 100
4 2014-09-05  10
5 2014-09-06  90
6 2014-09-08  40

penjualan adalah sebuah variabel sama seperti variabel pada bahasa pemograman lainnya. Function read.csv() akan menghasilkan nilai dalam tipe data dataframe. Sebagai informasi, karena R adalah bahasa pemograman khusus untuk pengolahan statistik, ia memiliki beberapa tipe data bawaan yang unik seperti vector, matrix, dan dataframe. Sebuah dataframe mirip seperti matrix tetapi masing-masing kolom boleh memiliki tipe data berbeda. Cara paling gampang untuk memahami dataframe adalah dengan menganggapnya sebagai sebuah sheet Excel dimana masing-masing kolom mewakili nilai sebuah variabel (pada statistik) dalam bentuk vector.

Salah satu masalah yang sering timbul pada saat membaca data adalah permasalahan yang berkaitan dengan tipe data. Seperti apa tipe data setiap kolom di penjualan? Saya dapat melihatnya dengan memberikan perintah berikut ini:

> sapply(penjualan_produk1, class)
  tanggal       qty 
 "factor" "integer" 

factor adalah tipe data di R untuk nilai yang terdiri atas pilihan atau kategori, seperti “iya” atau “tidak”. Pada data yang saya baca, nama konsumen dan nama produk sudah tepat bila dianggap sebagai factor karena nama konsumen atau nama produk adalah variabel kualitatif/kategorikal. Tapi hal yang sama tidak berlaku pada kolom tanggal! Kolom tanggal bukan variabel kualitatif melainkan continuous variable yang harus dibaca dalam bentuk Date. Oleh sebab itu, saya perlu membaca file CSV dengan perintah seperti berikut ini:

> penjualan_produk1 <- read.csv('penjualan_produk1.csv', header=TRUE,
+   colClasses=c('Date', 'integer')
+ )
> sapply(penjualan_produk1,class)
  tanggal       qty 
   "Date" "integer" 

Sebagai latihan, saya akan mencoba membandingkan penjualan produk1 dan produk2. Untuk itu, saya kembali membaca sebuah file CSV baru yang berisi data penjualan untuk produk2 dengan memberikan perintah berikut ini:

> penjualan_produk2 <- read.csv('penjualan_produk2.csv', header=TRUE,
+   colClasses=c('Date', 'integer')
+ )
> head(penjualan_produk2)
     tanggal qty
1 2014-09-02  80
2 2014-09-03  90
3 2014-09-04  90
4 2014-09-05  40
5 2014-09-06  30
6 2014-09-08 150

Langkah berikutnya yang perlu saya lakukan adalah membuat sebuah dataframe baru yang mengkombinasikan kolom qty untuk penjualan produk1 dan kolom qty untuk penjualan produk2. Pada istilah statistik, proses ini disebut sebagai data reshaping. Programmer akan lebih mengenalnya dengan istilah sepeti inner join, left join, aggregation, dan sebagainya 🙂 Sebagai contoh, saya memberikan perintah berikut ini:

> require(plyr)
Loading required package: plyr
> p <- join(penjualan_produk1, penjualan_produk2, 'tanggal')
> head(p)
     tanggal qty qty
1 2014-09-02  20  80
2 2014-09-03  30  90
3 2014-09-04 100  90
4 2014-09-05  10  40
5 2014-09-06  90  30
6 2014-09-08  40 150

Saya memakai require() untuk membaca package dan menyertakan function yang ada di dalam package tersebut sehingga dapat dipakai. Salah satu kelebihan R adalah jumlah package-nya yang sangat banyak sekali! Pada repository resmi-nya di http://cran.r-project.org/web/packages disebutkan bahwa saat ini terdapat 6.010 package yang tersedia. Pengguna R dapat men-install package dengan perintah install.packages() atau memilih menu Tools, Install Packages… di R Studio. Bila dirasa belum cukup, siapa saja boleh membuat package baru dan mendistribusikannya di situs pribadi bahkan bisa langsung di-download oleh R dari GitHub.com.

Pada perintah di atas, saya memakai join() untuk menggabungkan kolom qty dari dataframe penjualan_produk1 dan penjualan_produk2 berdasarkan tanggal. Hasil penggabungan kemudian disimpan pada sebuah variabel p. Masalah baru yang muncul adalah kini terdapat dua kolom dengan nama yang sama di p, yaitu qty. Untuk mengubah nama kolom, saya memberikan perintah berikut ini:

> names(p) <- c('tanggal', 'qty_produk1', 'qty_produk2')
> head(p)
     tanggal qty_produk1 qty_produk2
1 2014-09-02          20          80
2 2014-09-03          30          90
3 2014-09-04         100          90
4 2014-09-05          10          40
5 2014-09-06          90          30
6 2014-09-08          40         150

Saya bisa melakukan agregasi untuk memperoleh rata-rata dari qty_produk2 per qty_produk1 dengan menggunakan ddply dari package plyr. Sebagai contoh, perintah berikut ini setara dengan SELECT AVG(qty_produk2) FROM p GROUP BY qty_produk1 di SQL:

> p <- ddply(p, ~ qty_produk1, summarize, qty_produk2 = mean(qty_produk2, na.rm=TRUE))
> head(p)
  qty_produk1 qty_produk2
1           5    57.50000
2          10    52.66667
3          16         NaN
4          20    41.66667
5          23         NaN
6          25    20.00000

Sekarang, saya akan mencoba metode klasik di statistik dengan melakukan fitting data di atas ke sebuah persamaan linear Y = B0 * X + B1 (metode ini disebut linear regression). Karena saya memakai komputer, tentu saja saya tidak perlu menghitung nilai koefisien B0 dan B1 secara manual. Saya dapat memperoleh hasil perhitungan cukup dengan menggunakan perintah lm() seperti yang terlihat pada perintah ini:

> fit1 <- lm(qty_produk1 ~ qty_produk2, data = p)
> summary(fit1)

Call:
lm(formula = qty_produk1 ~ qty_produk2, data = p)

Residuals:
    Min      1Q  Median      3Q     Max 
-46.281 -26.616  -8.119  29.186  57.680 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  64.6446    17.4663   3.701  0.00237 **
qty_produk2  -0.2324     0.3221  -0.722  0.48246   
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 35.12 on 14 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.03585,   Adjusted R-squared:  -0.03301 
F-statistic: 0.5206 on 1 and 14 DF,  p-value: 0.4825

Saya menyimpan hasil dari lm() pada variabel fit1. Setelah itu saya menggunakan summary() untuk menampilkan informasi yang berkaitan dengan hasil fitting tersebut.

Pada perintah lm() di atas, R menghitung nilai coefficients dengan menggunakan metode linear least squares (disebut juga ordinary least squares). Bila persamaan linear adalah Y = B0 * X + B1, maka estimasi nilai B0 adalah -0,2324 dan nilai B1 adalah 64,6446. Ini berarti peningkatan penjualan produk2 akan menyebabkan penurunan penjualan produk1 sebesar 0,2324.

Walaupun demikian, dapat dilihat bahwa nilai B0 memiliki p-value sebesar 0,48246: sebuah nilai yang terlalu besar. Ini menunjukkan bahwa terdapat kemungkinan besar peluang yang mengarah pada null hyphothesis yaitu penjualan produk2 tidak memiliki hubungan dengan penjualan produk1 (hubungan yang terlihat kemungkinan besar hanya kebetulan atau akibat kesalahan). Output dari summary() secara otomatis memberi tanda bintang seperti ***, **, *, dan . pada koefisien yang signifikan secara statistik.

Nilai lainnya yang penting adalah R-squared. Nilai R-squared (disebut juga coefficient of determination) yang semakin mendekati 1 menunjukkan bahwa persamaan linear semakin mendekati data. Pada output summary() di atas, terlihat bahwa persamaan linear hanya mencakup 1% (0.01076) data.

R juga memiliki fasilitas yang mempermudah pengguna dalam menggambarkan grafis. Analisa visual tetap berguna dalam melengkapi analisa numeris. Untuk menggambar grafis, saya bisa menggunakan fungsi bawaan R yaitu plot(). Akan tetapi, pada artikel ini, saya akan memakai function pada package ggplot2 yang perlu di-download secara terpisah. Sebagai contoh, saya dapat menggambarkan nilai p dalam bentuk plot titik dengan menggunakan perintah seperti berikut ini:

> require(ggplot2)
Loading required package: ggplot2
> ggplot(p, aes(x=qty_produk2, y=qty_produk1)) + geom_point()
Hasil plot dengan ggplot()

Hasil plot dengan ggplot()

Saya bisa menambahkan geom_smooth() untuk menampilkan persamaan linear yang dihasilkan oleh lm() dengan menggunakan perintah berikut ini:

> ggplot(p, aes(x=qty_produk2, y=qty_produk1)) + geom_point() + geom_smooth(method='lm')
Hasil plot dengan ggplot()

Hasil plot dengan ggplot()

Sebagai informasi tambahan, R juga memiliki cor() dan cov() untuk menghitung nilai correlation dan covariance secara mudah. Sebagai contoh, saya dapat memberikan perintah seperti berikut ini:

> cov(p$qty_produk1, p$qty_produk2, use='pairwise')
[1] -184.1778
> cor(p$qty_produk1, p$qty_produk2, use='pairwise')
[1] -0.1893496

Nilai covariance yang negatif (-184,1778) sehingga menunjukkan bahwa terdapat hubungan berbanding terbalik (bila jumlah produk1 bertambah maka jumlah produk2 cenderung berkurang). Nilai correlation yang mendekati nol (-0.1893496) menunjukkan variabel yang dibandingkan tidak saling mempengaruhi.