knitr::opts_chunk$set(echo = TRUE)

Ilustrasi Kasus

Dataset yang akan digunakan sebagai ilustrasi kasus dalam diskusi kita adalah dataset “MIOPI”. Dataset ini memuat informasi dari 618 siswa mengenai kondisi rabun jauh dan faktor resikonya. Setiap siswa diidentifikasi apakah menderita myopi (rabun jauh) atau tidak.

Selain itu terdapat beberapa data faktor resiko miopi yang dikumpulkan dari setiap siswa yaitu: (1) karakteristik fisik mata, (2) durasi aktivitas luar ruangan, (3) durasi aktivitas dalam ruangan, dan (4) kondisi apakah orangtua menderita miopi.

Berikut ini perintah untuk mengimport data menjadi dataframe di R.

miopi = read.csv("C:/bagusco/miopi.csv")

Banyaknya baris dan kolom dari dataset kita adalah sebagai berikut:

dim(miopi)
## [1] 618  19

Nama-nama kolom pada dataset adalah sebagai berikut:

colnames(miopi)
##  [1] "X"         "id"        "studyyear" "myopic"    "age"       "gender"   
##  [7] "spheq"     "al"        "acd"       "lt"        "vcd"       "sporthr"  
## [13] "readhr"    "comphr"    "studyhr"   "tvhr"      "diopterhr" "mommy"    
## [19] "dadmy"

Berapa persen kejadian miopi pada siswa yang ada dalam data kita? Berikut ini adalah perintah untuk menampilkan frekuensi kejadian miopi dan menyajikannya dalam bentuk piechart.

table(miopi$myopic)
## 
##  No Yes 
## 537  81
prop.table(table(miopi$myopic))
## 
##       No      Yes 
## 0.868932 0.131068
pie(table(miopi$myopic))

Tampak bahwa siswa yang menderita miopi terdapat sekitar 13%.

Melihat Hubungan Variabel Penjelas (Faktor Resiko) dengan Kejadian Miopi

Sebelum lebih jauh menyusun model klasifikasi, kita dapat mendeskripsikan secara sederhana bagaimana hubungan antara faktor resiko (yang akan menjadi variabel penjelas dalam model) dengan kejadian miopi.

Perintah di bawah ini menyajikan ilustrasi untuk melihat keterkaitan antara durasi aktivitas outdoor (seperti olahraga) dengan kejadian miopi. Untuk memudahkan visualisasi, data lama melakukan aktivitas outdoor dikategorikan/dikelompokkan ke dalam 4 (empat) rentang nilai yaitu kurang dari 10 jam per minggu, 10-19 jam, 20-29 jam, dan 30 jam atau lebih.

Berdasarkan gambar yang diperoleh, kita mendapatkan informasi bahwa pada kelompok siswa yang memiliki aktivitas outdoor lebih lama, persentase kejadian miopi-nya cenderung lebih kecil. Dengan kata lain, resiko seorang siswa mengalami miopi cenderung lebih rendah pada siswa yang lebih banyak aktivitas outdoor-nya.

mengalami.miopi = miopi$myopic
sporthr.c = ifelse(miopi$sporthr < 10, 1, ifelse(miopi$sporthr < 20, 2, ifelse(miopi$sporthr < 30, 3, 4)))
kejadian.sporthr = prop.table(table(sporthr.c, mengalami.miopi), margin=1)
plot(kejadian.sporthr, col=c(1,2), main="hubungan SPORT HOUR dengan Kejadian Miopi")

Kita bisa lanjutkan juga untuk melihat keterkaitan antara durasi aktivitas dengan pandangan jarak dekat (seperti nonton TV, bermain video games, membaca buku) dengan kejadian miopi. seperti sebelumnya, kita juga membagi nilai durasi ini ke dalam 4 (empat) rentang nilai.

Berdasarkan gambar yang diperoleh, kita mendapatkan informasi bahwa pada kelompok siswa yang memiliki aktivitas pandangan jarak dekat lebih lama, persentase kejadian miopi-nya cenderung lebih besar. Dengan kata lain, resiko seorang siswa mengalami miopi cenderung lebih tinggi pada siswa yang lebih banyak aktivitas pandangan jarak dekat.

diopter.c = ifelse(miopi$diopter < 15, 1, ifelse(miopi$diopter < 30, 2, ifelse(miopi$diopter < 45, 3, 4)))
kejadian.diopter = prop.table(table(diopter.c, mengalami.miopi), margin=1)
plot(kejadian.diopter, col=c(1,2), main="hubungan DIOPTER HOUR dengan Kejadian Miopi")

Gambar di bawah ini menyajikan hubungan antara ukuran “spheq” dengan resiko terjadinya miopi. Nampak bahwa semakin kecil nilai “spheq” mata seorang siswa, resiko menderita miopi lebih tinggi.

spheq.c = ifelse(miopi$spheq < 0.1, 1, ifelse(miopi$spheq < 0.5, 2, ifelse(miopi$spheq < 1, 3, 4)))
kejadian.spheq = prop.table(table(spheq.c, mengalami.miopi), margin=1)
plot(kejadian.spheq, col=c(1,2), main="hubungan kondisi SPHEQ dengan Kejadian Miopi")

Terakhir, kita dapat melihat apakah resiko seorang siswa menderita miopi juga ada kaitannya dengan apakah ibu-nya menderita miopi. Gambar di bawah ini mengindikasikan bahwa siswa yang memiliki ibu miopi akan memiliki resiko miopi yang lebih besar dibandingkan siswa yang ibunya tidak menderita miopi.

kejadian.mommy = prop.table(table(miopi$mommy, mengalami.miopi), margin=1)
plot(kejadian.mommy, col=c(1,2), main="hubungan Ibu Miopi dengan Kejadian Miopi")

Pemodelan Regresi Logistik

Model klasifikasi pertama yang disusun adalah regresi logistik. Model ini menggunakan kejadian miopi pada siswa sebagai variabel respon.

Variabel penjelas yang digunakan adalah

Fungsi di R yang bisa digunakan adalah glm() dengan opsi family="binomial". Sebelum itu, variabel myopic diubah dulu nilainya menjadi 1 (untuk Yes) dan 0 (untuk No)

miopi$myopic.1 = ifelse(miopi$myopic == "Yes", 1, 0)
model.reglog = glm(myopic.1 ~ spheq + readhr + sporthr + mommy + dadmy, data=miopi, family="binomial")
model.reglog
## 
## Call:  glm(formula = myopic.1 ~ spheq + readhr + sporthr + mommy + dadmy, 
##     family = "binomial", data = miopi)
## 
## Coefficients:
## (Intercept)        spheq       readhr      sporthr     mommyYes     dadmyYes  
##    -0.57011     -3.81152      0.06712     -0.05173      0.77121      0.86870  
## 
## Degrees of Freedom: 617 Total (i.e. Null);  612 Residual
## Null Deviance:       480.1 
## Residual Deviance: 314.5     AIC: 326.5

Berdasarkan output di atas dapat dilihat bahwa:

Pengujian secara statistik apakah variabel-variabel penjelas tersebut signifikan efeknya disajikan pada output di bawah ini.

summary(model.reglog)
## 
## Call:
## glm(formula = myopic.1 ~ spheq + readhr + sporthr + mommy + dadmy, 
##     family = "binomial", data = miopi)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.57011    0.45868  -1.243  0.21389    
## spheq       -3.81152    0.43383  -8.786  < 2e-16 ***
## readhr       0.06712    0.04569   1.469  0.14180    
## sporthr     -0.05173    0.02010  -2.573  0.01007 *  
## mommyYes     0.77121    0.31090   2.481  0.01312 *  
## dadmyYes     0.86870    0.30666   2.833  0.00461 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 480.08  on 617  degrees of freedom
## Residual deviance: 314.55  on 612  degrees of freedom
## AIC: 326.55
## 
## Number of Fisher Scoring iterations: 6

Perintah di bawah ini menyajikan salah satu cara bagaimana mengevaluasi kebaikan model klasifikasi yang telah didapatkan. Caranya adalah dengan melihat tingkat ketepatan prediksinya. Tahapan yang dilakukan kira-kira adalah:

  1. lakukan prediksi terhadap anak-anak yang sebenarnya diketahui apakah dia menderita miopi atau tidak menggunakan model regresi logistik yang telah didapatkan –> hasil prediksi berupa peluang terjadinya miopi

  2. lakukan pengkelasan nilai prediksi peluang/resiko miopi. Jika prediksi peluangnya di atas 0.5 maka disimpulkan diprediksi miopi, jika kurang dari 0.5 disimpulkan tidak miopi

  3. bandingkan kelas prediksi dengan kelas aktual

  4. hitung berapa persen ketepatan prediksinya menggunakan ukuran akurasi, sensitivity, atau specificity.

Berdasarkan output di bawah ini, model regresi logistik kita memiliki akurasi 89.3% yang berarti bahwa 89.3% dugaan model sesuai dengan kondisi aktualnya.

library(caret)
prediksi.reglog.p = predict(model.reglog, newdata=miopi, type='response')
prediksi.reglog = ifelse(prediksi.reglog.p < 0.5, 0, 1)
confusionMatrix(as.factor(prediksi.reglog), as.factor(miopi$myopic.1), positive="1")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 521  50
##          1  16  31
##                                           
##                Accuracy : 0.8932          
##                  95% CI : (0.8661, 0.9164)
##     No Information Rate : 0.8689          
##     P-Value [Acc > NIR] : 0.03915         
##                                           
##                   Kappa : 0.4295          
##                                           
##  Mcnemar's Test P-Value : 4.865e-05       
##                                           
##             Sensitivity : 0.38272         
##             Specificity : 0.97020         
##          Pos Pred Value : 0.65957         
##          Neg Pred Value : 0.91243         
##              Prevalence : 0.13107         
##          Detection Rate : 0.05016         
##    Detection Prevalence : 0.07605         
##       Balanced Accuracy : 0.67646         
##                                           
##        'Positive' Class : 1               
## 

Model Pohon Klasifikasi

library(rpart)
library(rpart.plot)
model.pohon = rpart(myopic.1 ~ spheq + readhr + sporthr + mommy + dadmy, data=miopi, method="class")
model.pohon
## n= 618 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 618 81 0 (0.868932039 0.131067961)  
##    2) spheq>=0.2845 519 32 0 (0.938342967 0.061657033)  
##      4) spheq>=0.7235 311  3 0 (0.990353698 0.009646302) *
##      5) spheq< 0.7235 208 29 0 (0.860576923 0.139423077)  
##       10) sporthr>=3.5 186 20 0 (0.892473118 0.107526882) *
##       11) sporthr< 3.5 22  9 0 (0.590909091 0.409090909)  
##         22) dadmy=No 11  2 0 (0.818181818 0.181818182) *
##         23) dadmy=Yes 11  4 1 (0.363636364 0.636363636) *
##    3) spheq< 0.2845 99 49 0 (0.505050505 0.494949495)  
##      6) spheq>=-0.0515 71 26 0 (0.633802817 0.366197183)  
##       12) dadmy=No 35  8 0 (0.771428571 0.228571429)  
##         24) mommy=No 13  0 0 (1.000000000 0.000000000) *
##         25) mommy=Yes 22  8 0 (0.636363636 0.363636364)  
##           50) spheq< 0.205 14  3 0 (0.785714286 0.214285714) *
##           51) spheq>=0.205 8  3 1 (0.375000000 0.625000000) *
##       13) dadmy=Yes 36 18 0 (0.500000000 0.500000000)  
##         26) sporthr>=5.5 29 12 0 (0.586206897 0.413793103)  
##           52) spheq< 0.161 15  4 0 (0.733333333 0.266666667) *
##           53) spheq>=0.161 14  6 1 (0.428571429 0.571428571) *
##         27) sporthr< 5.5 7  1 1 (0.142857143 0.857142857) *
##      7) spheq< -0.0515 28  5 1 (0.178571429 0.821428571) *
rpart.plot(model.pohon)

library(caret)
prediksi.pohon.p = predict(model.pohon, newdata=miopi, type='prob')[,2]
prediksi.pohon = ifelse(prediksi.pohon.p < 0.5, 0, 1)
confusionMatrix(as.factor(prediksi.pohon), as.factor(miopi$myopic.1), positive="1")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 518  32
##          1  19  49
##                                           
##                Accuracy : 0.9175          
##                  95% CI : (0.8929, 0.9379)
##     No Information Rate : 0.8689          
##     P-Value [Acc > NIR] : 0.0001003       
##                                           
##                   Kappa : 0.6112          
##                                           
##  Mcnemar's Test P-Value : 0.0928919       
##                                           
##             Sensitivity : 0.60494         
##             Specificity : 0.96462         
##          Pos Pred Value : 0.72059         
##          Neg Pred Value : 0.94182         
##              Prevalence : 0.13107         
##          Detection Rate : 0.07929         
##    Detection Prevalence : 0.11003         
##       Balanced Accuracy : 0.78478         
##                                           
##        'Positive' Class : 1               
##