knitr::opts_chunk$set(echo = TRUE)
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%.
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")
Model klasifikasi pertama yang disusun adalah regresi logistik. Model ini menggunakan kejadian miopi pada siswa sebagai variabel respon.
Variabel penjelas yang digunakan adalah
spheq
readhr
sporthr
mommy
dadmy
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:
spheq memiliki koefisien negatif, yang berarti bahwa semakin besar nilai spheq dari mata anak peluang/resiko miopi cenderung mengecil
readhr memiliki koefisien positif, yang berarti bahwa semakin sering seorang anak melakukan aktivitas membaca jarak dekat… maka peluang/resiko miopi cenderung meningkat
sporthr memiliki koefisien negatif, yang berarti bahwa semakin sering seorang anak melakukan aktivitas luar ruangan (outdoor)… maka peluang/resiko miopi cenderung mengecil
Ibu Miopi dan Bapak Miopi memiliki koefisien positif, yang berarti bahwa kalau orangtua menderita miopi… maka peluang/resiko miopi anak cenderung lebih tinggi dari anak yang orangtuanya tidak miopi.
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:
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
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
bandingkan kelas prediksi dengan kelas aktual
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
##
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
##