Klasifikasi bertujuan untuk memprediksi kelas (target variable kategorik):

  • binary classification: 2 kelas
  • multiclass classification: > 2 kelas

Logistic Regression

Berapa range peluang?

  • min: 0
  • max: 1

Berapa range hasil prediksi model regresi?

  • min: -inf
  • max: inf

Dibutuhkan suatu jembatan agar regression dapat digunakan untuk memprediksi peluang. Jembatan itu adalah Odds dan Log of Odds.

Odds & Log of Odds

Odds adalah bentuk lain dari peluang, yaitu perbandingan peluang antara peluang kejadian terjadi/peluang kejadian tidak terjadi.

\[\frac{p}{(1-p)}\]

p = peluang suatu kejadian terjadi

Anda hendak berpergian menggunakan pesawat dari Soekarno Hatta Airport. Bila diketahui dari 100 penerbangan di Soekarno Hatta, terdapat 25 pesawat Delay. Berapa odds pesawat Anda On Time?

# peluang
p <- 75/100
# odds
p/(1-p)
#> [1] 3

Interpretasi:

  • kemungkinan pesawat on time 3 kali lebih mungkin

Note: Kalau oddsnya 1 berarti peluangnya? 0.5

Berapa range nilai dari odds?

# odds: p/1-p
# min
0/(1-0)
#> [1] 0
# max
1/(1-1)
#> [1] Inf
  • Probability: 0 1
  • Odds : 0 Inf

Log of odds adalah nilai odds yang dilogaritmikkan:

\(logit(p) = log(\frac{p}{1-p})\)

# log of odds - pesawat on time
log(3)
#> [1] 1.098612

Berapa range nilai log of odds?

# min
log(0/(1-0))
#> [1] -Inf
# max
log(1/(1-1))
#> [1] Inf
  • Probability: 0 1
  • Odds : 0 inf
  • Log of odds: -inf inf

Odds dan log of odds mampu menjembatani antara nilai yang dihasilkan oleh model regresi, ke rentang nilai peluang. Logistic regression menghasilkan Log of Odds. Nilai log of odds dapat dikembalikan ke bentuk odds (untuk diinterpretasikan) dan peluang sehingga dapat digunakan untuk klasifikasi.

# log of odds -> odds
odds <- exp(1.386294)
odds
#> [1] 3.999999

Terdapat cara lain:

  • logit(): peluang -> log of odds
  • inv.logit(): log of odds -> peluang
library(gtools)

# peluang -> log of odds dengan `logit()`
logit(0.8)
#> [1] 1.386294
# log of odds -> peluang dengan `inv.logit()`
inv.logit(1.386294)
#> [1] 0.7999999

Fungsi inv.logit() ini juga biasa disebut sigmoid function.

# sigmoid function
curve(inv.logit(x), from = -10, to=10, 
      xlab = "Log of Odds", 
      ylab = "Peluang")

glm() & Interpretation

Anda adalah seorang analis performa student di universitas. Anda ditugaskan untuk memprediksi status kelulusan siswa dengan honors (cumlaude).

# read & inspect data
library(dplyr)
honors <- read.csv("data_input/sample.csv") %>% 
  select(-femalexmath)

glimpse(honors)
#> Rows: 200
#> Columns: 5
#> $ female <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
#> $ read   <int> 57, 68, 44, 63, 47, 44, 50, 34, 63, 57, 60, 57, 73, 54, 45, 42,…
#> $ write  <int> 52, 59, 33, 44, 52, 52, 59, 46, 57, 55, 46, 65, 60, 63, 57, 49,…
#> $ math   <int> 41, 53, 54, 47, 57, 51, 42, 45, 54, 52, 51, 51, 71, 57, 50, 43,…
#> $ hon    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, …

Deskripsi kolom:

  • female: gender of student (1 for female)
  • read: score in reading test
  • write: score in writing test
  • math: score in math test
  • hon: status of graduating in honors (1 for honors)
# cek missing value
anyNA(honors)
#> [1] FALSE

Cara membuat model logistic regression:

glm(target ~ prediktor, data, family = "binomial")

honors.logit <- glm(formula = hon ~ 1, data = honors, family = "binomial")

summary(honors.logit)
#> 
#> Call:
#> glm(formula = hon ~ 1, family = "binomial", data = honors)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -0.7497  -0.7497  -0.7497  -0.7497   1.6772  
#> 
#> Coefficients:
#>             Estimate Std. Error z value         Pr(>|z|)    
#> (Intercept)  -1.1255     0.1644  -6.845 0.00000000000762 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 222.71  on 199  degrees of freedom
#> Residual deviance: 222.71  on 199  degrees of freedom
#> AIC: 224.71
#> 
#> Number of Fisher Scoring iterations: 4

Intercept: log of odds dari target (student mendapatkan honors (1))

Berikut pembuktiannya:

# peluang
table(honors$hon)
#> 
#>   0   1 
#> 151  49
# log of odds student honors
logit(49/200)
#> [1] -1.12546

Interpretasi: Log of odds tidak dapat diinterpretasikan. Untuk interpretasi, nilai log of odds kita ubah ke odds.

# log of odds -> odds
exp(-1.12546)
#> [1] 0.3245032

berarti kemungkinan calon siswa mendapatkan honors 0.32 kali lebih mungkin dari pada tidak mendapatkan honors

kemungkinan calon siswa lebih kecil mendapatkan honor dari pada tidak mendapatkan honors

1 Prediktor Kategorik

Buat model untuk memprediksi honors berdasarkan gender female:

honors.logit2 <- glm(formula = hon ~ female  , data = honors , family = "binomial")

summary(honors.logit2)
#> 
#> Call:
#> glm(formula = hon ~ female, family = "binomial", data = honors)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -0.8337  -0.8337  -0.6431  -0.6431   1.8317  
#> 
#> Coefficients:
#>             Estimate Std. Error z value     Pr(>|z|)    
#> (Intercept)  -1.4709     0.2690  -5.469 0.0000000453 ***
#> female        0.5928     0.3414   1.736       0.0825 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 222.71  on 199  degrees of freedom
#> Residual deviance: 219.61  on 198  degrees of freedom
#> AIC: 223.61
#> 
#> Number of Fisher Scoring iterations: 4

Female: log of odds ratio dari student female mendapatkan honors dibandingkan student male mendapatkan honors.

# proportion
table(female = honors$female, honors = honors$hon)
#>       honors
#> female  0  1
#>      0 74 17
#>      1 77 32
# peluang
p_female <- 32/(32+77)
p_male <- 17/91

# odds 
o_female <- p_female / (1-p_female)
o_male <- p_male / (1-p_male)

# log of odds
log(o_female/o_male) 
#> [1] 0.5927822
# log dari odds female mendapakan honors dibagi dengan log dari odd male mendapakan honors

log dari odds fema

Intercept: log of odds dari student male yang mendapatkan honors (basis)

log(o_male)
#> [1] -1.470852

Interpretasi:

# odds female dapat honors
exp(0.5927822)
#> [1] 1.809014

kemungkinan perempuan mendapatkan honors 1.8 kali lebih mungkin dari pada pria mendapatkan honors

Buat model untuk memprediksi honors berdasarkan nilai math:

honors.logit3 <- glm(formula = hon ~ math, data = honors, family = "binomial")

summary(honors.logit3)
#> 
#> Call:
#> glm(formula = hon ~ math, family = "binomial", data = honors)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.0332  -0.6785  -0.3506  -0.1565   2.6143  
#> 
#> Coefficients:
#>             Estimate Std. Error z value        Pr(>|z|)    
#> (Intercept) -9.79394    1.48174  -6.610 0.0000000000385 ***
#> math         0.15634    0.02561   6.105 0.0000000010294 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 222.71  on 199  degrees of freedom
#> Residual deviance: 167.07  on 198  degrees of freedom
#> AIC: 171.07
#> 
#> Number of Fisher Scoring iterations: 5

Intercept: -9.79394

Math: 0.15634

Contoh:

hon = -9.79394 + 0.15634 * math

Student A memiliki nilai math 52, student B 53. Hitung masing-masing log of oddsnya, berapa selisihnya?

# hint: substitusi formula model saja
# log of odds
hon52 <-  -9.79394 + 0.15634 * 52
hon53 <- -9.79394 + 0.15634 * 53

hon53-hon52
#> [1] 0.15634

Interpretasi:

# log of odds -> odds

exp(hon53-hon52)
#> [1] 1.169224

dengan menigkatnya 1 point dari nilai matematika maka siswa 1.17 (1.169224) kali lebih mungkin mendapatkan honors

Banyak Prediktor

Buat model untuk memprediksi honors berdasarkan gender female dan nilai math:

honors.logit4 <- glm(hon ~ female + math, data = honors , family= "binomial")


summary(honors.logit4)
#> 
#> Call:
#> glm(formula = hon ~ female + math, family = "binomial", data = honors)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -1.8494  -0.6506  -0.3471  -0.1361   2.5105  
#> 
#> Coefficients:
#>              Estimate Std. Error z value        Pr(>|z|)    
#> (Intercept) -10.80595    1.61654  -6.685 0.0000000000232 ***
#> female        0.96531    0.41599   2.321          0.0203 *  
#> math          0.16422    0.02665   6.161 0.0000000007206 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 222.71  on 199  degrees of freedom
#> Residual deviance: 161.35  on 197  degrees of freedom
#> AIC: 167.35
#> 
#> Number of Fisher Scoring iterations: 5

Interpretasi koefisien:

cari odds dari masing masing predictor:

# female
exp(0.96531 )
#> [1] 2.625601
# math
exp(0.16422)
#> [1] 1.178474

female = 2.625601

kemungkinan female mendapatkan honors 2.63 kali lebih mungkin dari pada pria mendapatkan honors ketika nilai matematikanya sama/konstan

math = 1.178474

setiap pengingkatan 1 point pada matematika mengakibatkan 1.18 kali lebih mungkin mendapatkan honors ketika variable lainnya konstan

jika seorang siswa memiliki nilai 1 point lebih tinggi nilai matematikanya dari siswa lainnya maka dia akan 1.18 kali lebih mungkin mendapatkan honors ketika varible lainnya konstan

Aplikasi:

Final formula: y = -10.80595 + 0.96531 * female1 + 0.16422 * math

case

Joe adalah seorang male yang nilai math-nya 60, berapa peluang dia mendapatkan honors? Apakah dia akan lulus dengan honors?

joe <- -10.80595 + 0.96531 * 0 + 0.16422 *60

joe
#> [1] -0.95275
inv.logit(joe)
#> [1] 0.2783321

probability untuk dia lulus dengan honors adalah 27.8 %

ifelse(inv.logit(joe) < 0.5 , "tidak honors", "honors" )
#> [1] "tidak honors"

Summary 1

  1. Logistic regression menghasilkan log of odds
  • fungsi: glm( formula, data, family = "binomial")
  1. Untuk interpretasi model logistic regression,
  • dilakukan: dirubah ke odds
  • menggunakan fungsi: exp()

Interpretasi koefisien dapat dilakukan dan berbeda untuk masing-masing kondisi:

  • tanpa prediktor : log of odds dari target variable nya
  • 1 prediktor kategorik : odds rasio dari salah satu levels
  • 1 prediktor numerik : odds dari setiap kenaikan 1 point
  • banyak prediktor : sama seperti biasa namun asumsi prediktor yang lain konstan

Bila koefisien variable:

  • positif: meningkatkan kemungkinan ke target variable 1
  • negatif: menurunkan kemungkinan ke target variable 1
  1. Untuk menentukan kelas (klasifikasi) dari hasil logistic regression, nilai log of odds diubah kebentuk probability kemudian ditentukan kelasnya berdasarkan batas tertentu (misal = 0.5).

Perfect Separation

  • Null deviance: deviasi model saat tanpa prediktor (model terburuk).
  • Residual deviance: deviasi model saat menggunakan prediktor.

Umumnya semakin banyak prediktor maka residual deviance akan semakin kecil.

summary(honors.logit)
#> 
#> Call:
#> glm(formula = hon ~ 1, family = "binomial", data = honors)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -0.7497  -0.7497  -0.7497  -0.7497   1.6772  
#> 
#> Coefficients:
#>             Estimate Std. Error z value         Pr(>|z|)    
#> (Intercept)  -1.1255     0.1644  -6.845 0.00000000000762 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 222.71  on 199  degrees of freedom
#> Residual deviance: 222.71  on 199  degrees of freedom
#> AIC: 224.71
#> 
#> Number of Fisher Scoring iterations: 4
# null deviance
honors.logit$null.deviance
#> [1] 222.71
# residual deviance
honors.logit$deviance # wo/ predictor
#> [1] 222.71
honors.logit2$deviance # w/ female
#> [1] 219.6062
honors.logit3$deviance # w/ math
#> [1] 167.0732
honors.logit4$deviance # w/ female + math
#> [1] 161.3454

Mari buat model honors.logit5 untuk memprediksi honors berdasarkan semua prediktor yang ada:

honors.logit5 <- glm(hon ~ . , data = honors, family = "binomial")
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(honors.logit5)
#> 
#> Call:
#> glm(formula = hon ~ ., family = "binomial", data = honors)
#> 
#> Deviance Residuals: 
#>          Min            1Q        Median            3Q           Max  
#> -0.000245237  -0.000000021  -0.000000021  -0.000000021   0.000176979  
#> 
#> Coefficients:
#>                 Estimate   Std. Error z value Pr(>|z|)
#> (Intercept)  -2176.87191 282900.61006  -0.008    0.994
#> female          -4.44248  18078.08748   0.000    1.000
#> read             0.07478    489.40463   0.000    1.000
#> write           36.30917   4837.98255   0.008    0.994
#> math            -0.36898   1362.81267   0.000    1.000
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 222.71004668627  on 199  degrees of freedom
#> Residual deviance:   0.00000013527  on 195  degrees of freedom
#> AIC: 10
#> 
#> Number of Fisher Scoring iterations: 25

NOTE:

  • glm.fit: fitted probabilities numerically 0 or 1 occurred -> warning bahwa bisa dihasilkan probability yang tepat 1 atau 0 (indikasi kondisi perfect separation)
  • glm.fit: algorithm did not converge -> warning bahwa algoritmanya tidak mencapai kondisi stabil hingga iterasi ke-25 (default), dapat terjadi salah satunya karena kondisi perfect separation.

Perfect Separation adalah sebuah kondisi dimana ada 1 variabel yang dapat memisahkan kelas target secara sempurna. Cara mendeteksi:

  • Lihat dari estimate, kalau ada yang paling beda
  • p value, kalau p value terlalu banyak yang tidak signifikan
  • null dan residual deviance.
  • iteration lebih dari 25
# log of odds -> odds

data.frame(odds = exp(honors.logit5$coefficients))

Pada kasus ini, nilai write dapat memisahkan kelas honor dengan sempurna:

table(honors$hon, honors$write)
#>    
#>     31 33 35 36 37 38 39 40 41 42 43 44 45 46 47 49 50 52 53 54 55 57 59 60 61
#>   0  4  4  2  2  3  1  5  3 10  2  1 12  1  9  2 11  2 15  1 17  3 12 25  4  0
#>   1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  4
#>    
#>     62 63 65 67
#>   0  0  0  0  0
#>   1 18  4 16  7
plot(honors$hon, honors$write)

Tidak disarankan menggunakan model dengan perfect separation, karena model amat bias pada salah satu variable dan tidak mempertimbangkan variable lain. Hal ini dapat membuat model tidak akurat (buruk) dalam memprediksi ke data baru.

Apa yang kita lakukan bila bertemu kondisi perfect separation:

  • kalau kasus seperti ini kita terima, tidak usah membuat machine learning, cukup ifelse saja.
  • kalau kasus ini tidak kita terima, maka jangan gunakan variabel ini sebagai prediktor.
  • observasi (data) nya kita tambah
# revisi model
honors.logit6 <- glm(hon ~ female + read + math , data = honors, family = "binomial")
summary(honors.logit6)
#> 
#> Call:
#> glm(formula = hon ~ female + read + math, family = "binomial", 
#>     data = honors)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -1.8305  -0.6327  -0.3300  -0.1258   2.3896  
#> 
#> Coefficients:
#>              Estimate Std. Error z value         Pr(>|z|)    
#> (Intercept) -11.77025    1.71068  -6.880 0.00000000000597 ***
#> female        0.97995    0.42163   2.324           0.0201 *  
#> read          0.05906    0.02655   2.224           0.0261 *  
#> math          0.12296    0.03128   3.931 0.00008442731719 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 222.71  on 199  degrees of freedom
#> Residual deviance: 156.17  on 196  degrees of freedom
#> AIC: 164.17
#> 
#> Number of Fisher Scoring iterations: 5

AIC

AIC = Jumlah informasi yang hilang. Semakin kecil AIC, semakin baik model. Area Under the Curve (AUC) dan Akaike Information Criterion (AIC) adalah dua metrik yang digunakan untuk mengevaluasi model, tetapi keduanya memiliki tujuan yang berbeda dan diterapkan pada konteks yang berbeda. Berikut adalah perbedaan antara AUC dan AIC pada model regresi logistik:

Area Under the Curve (AUC):

Tujuan: AUC digunakan untuk mengevaluasi kinerja model klasifikasi, seperti regresi logistik. AUC mengukur seberapa baik model membedakan antara dua kelas dengan memplot kurva karakteristik operasi penerima (ROC). Interpretasi: Nilai AUC berkisar antara 0 dan 1. Model yang sempurna memiliki AUC setara dengan 1, sementara model yang tidak dapat membedakan antara kelas memiliki AUC sekitar 0.5 (sebagai acak). Akaike Information Criterion (AIC):

Tujuan: AIC digunakan untuk membandingkan model statistik yang berbeda berdasarkan kemungkinan (likelihood) model dan kompleksitasnya. Tujuannya adalah untuk memilih model yang memberikan keseimbangan terbaik antara kemungkinan model dan kompleksitasnya. Interpretasi: Pada AIC, model yang memiliki nilai AIC yang lebih rendah dianggap lebih baik. AIC memperhitungkan baik kecocokan model terhadap data (melalui likelihood) maupun kompleksitas model (melalui jumlah parameter). Dengan kata lain, AUC fokus pada kemampuan model untuk membedakan kelas (klasifikasi), sementara AIC fokus pada keseimbangan antara kecocokan model terhadap data dan kompleksitas model.

honors.logit$aic # wo/ predictor
#> [1] 224.71
honors.logit2$aic # w/ female
#> [1] 223.6062
honors.logit3$aic # w/ math
#> [1] 171.0732
honors.logit4$aic # w/ female + math
#> [1] 167.3454
honors.logit6$aic # w/ all predictor except write
#> [1] 164.1696
honors.logit5$aic # w/ write
#> [1] 10

Important Notes:

Dalam menseleksi model, model yang baik adalah:

  • model dengan nilai AIC rendah
  • model tanpa kondisi Perfect Separation

Assumption

Logistic Regression menganut 3 asumsi:

  • Linearity of Predictor & Log of Odds: cara interpretasi model mengacu pada asumsi ini (contoh: untuk variabel numerik, peningkatan 1 nilai akan meningkatkan log of odds)
  • Multicollinearity: antar prediktor tidak saling berkorelasi kuat (hingga nilai 1 / -1) -> uji vif()
  • Independence of Observations: antar observasi saling independen & tidak berasal dari pengukuran berulang (repeated measurement) -> kita harus ambil data secara random sampling

Asumsi logistic regression menuntut kita untuk memahami data secara mendalam dan memastikan data sudah siap dipakai untuk membuat model. Coba analisis kasus di bawah:

Case 2 Flight prediction

Berikut data penerbangan pesawat dalam flight_sm.csv:

flight <- read.csv("data_input/flight_sm.csv") %>% 
  mutate(DepDel15 = as.factor(DepDel15))
glimpse(flight)
#> Rows: 538,363
#> Columns: 10
#> $ Year        <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013…
#> $ Month       <int> 9, 9, 9, 7, 5, 7, 10, 10, 5, 6, 7, 8, 7, 10, 4, 10, 10, 8,…
#> $ DayofMonth  <int> 16, 23, 7, 15, 16, 28, 6, 9, 12, 30, 21, 4, 17, 20, 15, 2,…
#> $ DayofWeek   <int> 1, 1, 6, 1, 4, 7, 7, 3, 7, 7, 7, 7, 3, 7, 1, 3, 5, 6, 2, 1…
#> $ Carrier     <chr> "DL", "WN", "AS", "OO", "DL", "UA", "WN", "AA", "UA", "OO"…
#> $ CRSDepTime  <int> 1539, 1400, 810, 804, 805, 1251, 2000, 1010, 1522, 1404, 9…
#> $ DepDel15    <fct> 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
#> $ CRSArrTime  <int> 1824, 1425, 1614, 1027, 1117, 1700, 2135, 1240, 2017, 1507…
#> $ OriginState <chr> "FL", "PA", "WA", "IL", "NY", "TX", "CA", "IL", "NV", "CA"…
#> $ DestState   <chr> "NY", "IL", "DC", "OH", "FL", "DC", "NV", "TX", "TX", "CA"…

Dekspripsi kolom:

  • Year, Month, DayofMonth, DayofWeek: self-explanatory
  • Carrier: maskapai
  • CRSDepTime & CRSArrTime: jadwal departure & arrival (hhmm)
  • DepDel15: status delay (1 = delay)
  • OriginState, DestState: lokasi keberangkatan & tujuan

Buat model flight.model untuk memprediksi DepDel15 berdasarkan Month + DayofWeek, kemudian tampilkan summary-nya:

flight.model <- glm(DepDel15 ~ Month + DayofWeek , data = flight, family = "binomial")
  
summary(flight.model)
#> 
#> Call:
#> glm(formula = DepDel15 ~ Month + DayofWeek, family = "binomial", 
#>     data = flight)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -0.7308  -0.6911  -0.6518  -0.6168   1.8786  
#> 
#> Coefficients:
#>              Estimate Std. Error z value             Pr(>|z|)    
#> (Intercept) -0.935641   0.013900 -67.312 < 0.0000000000000002 ***
#> Month       -0.060895   0.001718 -35.454 < 0.0000000000000002 ***
#> DayofWeek   -0.004576   0.001710  -2.676              0.00744 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 541635  on 538362  degrees of freedom
#> Residual deviance: 540365  on 538360  degrees of freedom
#> AIC: 540371
#> 
#> Number of Fisher Scoring iterations: 4
flight_clean <- flight %>% 
  mutate(Month = as.factor(Month),
         DayofWeek = as.factor(DayofWeek))

glimpse(flight_clean)
#> Rows: 538,363
#> Columns: 10
#> $ Year        <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013…
#> $ Month       <fct> 9, 9, 9, 7, 5, 7, 10, 10, 5, 6, 7, 8, 7, 10, 4, 10, 10, 8,…
#> $ DayofMonth  <int> 16, 23, 7, 15, 16, 28, 6, 9, 12, 30, 21, 4, 17, 20, 15, 2,…
#> $ DayofWeek   <fct> 1, 1, 6, 1, 4, 7, 7, 3, 7, 7, 7, 7, 3, 7, 1, 3, 5, 6, 2, 1…
#> $ Carrier     <chr> "DL", "WN", "AS", "OO", "DL", "UA", "WN", "AA", "UA", "OO"…
#> $ CRSDepTime  <int> 1539, 1400, 810, 804, 805, 1251, 2000, 1010, 1522, 1404, 9…
#> $ DepDel15    <fct> 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
#> $ CRSArrTime  <int> 1824, 1425, 1614, 1027, 1117, 1700, 2135, 1240, 2017, 1507…
#> $ OriginState <chr> "FL", "PA", "WA", "IL", "NY", "TX", "CA", "IL", "NV", "CA"…
#> $ DestState   <chr> "NY", "IL", "DC", "OH", "FL", "DC", "NV", "TX", "TX", "CA"…
flight.model_baru <- glm(DepDel15 ~ Month + DayofWeek , data = flight_clean, family = "binomial")
  
summary(flight.model_baru)
#> 
#> Call:
#> glm(formula = DepDel15 ~ Month + DayofWeek, family = "binomial", 
#>     data = flight_clean)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -0.8740  -0.7155  -0.6314  -0.5429   2.0641  
#> 
#> Coefficients:
#>             Estimate Std. Error  z value             Pr(>|z|)    
#> (Intercept) -1.34623    0.01196 -112.608 < 0.0000000000000002 ***
#> Month5      -0.10637    0.01292   -8.235 < 0.0000000000000002 ***
#> Month6       0.35314    0.01222   28.899 < 0.0000000000000002 ***
#> Month7       0.23991    0.01225   19.585 < 0.0000000000000002 ***
#> Month8      -0.01661    0.01270   -1.308                0.191    
#> Month9      -0.36450    0.01382  -26.368 < 0.0000000000000002 ***
#> Month10     -0.32892    0.01351  -24.350 < 0.0000000000000002 ***
#> DayofWeek2  -0.16516    0.01275  -12.954 < 0.0000000000000002 ***
#> DayofWeek3  -0.05737    0.01251   -4.586      0.0000045233803 ***
#> DayofWeek4   0.22770    0.01206   18.879 < 0.0000000000000002 ***
#> DayofWeek5   0.13083    0.01224   10.690 < 0.0000000000000002 ***
#> DayofWeek6  -0.29298    0.01384  -21.172 < 0.0000000000000002 ***
#> DayofWeek7  -0.08550    0.01280   -6.677      0.0000000000243 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 541635  on 538362  degrees of freedom
#> Residual deviance: 534556  on 538350  degrees of freedom
#> AIC: 534582
#> 
#> Number of Fisher Scoring iterations: 4

Classification Workflow

  1. Business Question
  2. Read Data
  3. Data Wrangling
  4. EDA (untuk logistic regression selesaikan pengecekan asumsi disini)
  5. Cross Validation
  6. Data Pre-Processing
  7. Build Model
  8. Predict
  9. Model Evaluation
  10. Model Tuning -> Final Model

Case 3 Credit Risk Analysis

Studi Kasus: Credit Risk Analysis

Buat model untuk memprediksi peluang customer akan gagal bayar pinjaman (loan default), untuk mengindikasikan apakah customer tersebut baik atau tidak untuk diberikan pinjaman.

loans <- read.csv("data_input/loan2017Q4.csv", stringsAsFactors = T)

Data Wrangling

glimpse(loans)
#> Rows: 1,556
#> Columns: 16
#> $ initial_list_status <fct> w, f, w, w, w, w, w, w, w, w, w, w, w, f, w, w, w,…
#> $ purpose             <fct> debt_consolidation, debt_consolidation, debt_conso…
#> $ int_rate            <dbl> 14.08, 9.44, 28.72, 13.59, 15.05, 10.91, 15.05, 10…
#> $ installment         <dbl> 675.99, 480.08, 1010.30, 484.19, 476.33, 130.79, 3…
#> $ annual_inc          <dbl> 156700, 50000, 25000, 175000, 109992, 49000, 65000…
#> $ dti                 <dbl> 19.11, 19.35, 65.58, 12.60, 10.00, 5.12, 22.38, 33…
#> $ verification_status <fct> Source Verified, Not Verified, Verified, Not Verif…
#> $ grade               <fct> C, B, F, C, C, B, C, B, D, D, F, C, C, E, B, C, C,…
#> $ revol_bal           <int> 21936, 5457, 23453, 31740, 2284, 2016, 14330, 2758…
#> $ inq_last_12m        <int> 3, 1, 0, 0, 3, 5, 0, 1, 8, 1, 0, 12, 4, 8, 1, 3, 0…
#> $ delinq_2yrs         <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,…
#> $ home_ownership      <fct> MORTGAGE, RENT, OWN, MORTGAGE, MORTGAGE, MORTGAGE,…
#> $ not_paid            <int> 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1,…
#> $ log_inc             <dbl> 11.962088, 10.819778, 10.126631, 12.072541, 11.608…
#> $ verified            <int> 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ grdCtoA             <int> 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,…

Target: not_paid (paid = 0, not_paid = 1)

Adakah variabel yang tipe datanya belum sesuai?

  • grdCtoA
  • not_paid
  • verified

Adakah variabel yang dapat dibuang?

  • verified_status
  • grade
  • annual_inc
library(dplyr)

loan_clean <- loans %>% 
  select(-c(verification_status, grade, annual_inc)) %>% 
  mutate(grdCtoA = as.factor(grdCtoA),
         not_paid = as.factor(not_paid),
         verified = as.factor(verified))

glimpse(loan_clean)
#> Rows: 1,556
#> Columns: 13
#> $ initial_list_status <fct> w, f, w, w, w, w, w, w, w, w, w, w, w, f, w, w, w,…
#> $ purpose             <fct> debt_consolidation, debt_consolidation, debt_conso…
#> $ int_rate            <dbl> 14.08, 9.44, 28.72, 13.59, 15.05, 10.91, 15.05, 10…
#> $ installment         <dbl> 675.99, 480.08, 1010.30, 484.19, 476.33, 130.79, 3…
#> $ dti                 <dbl> 19.11, 19.35, 65.58, 12.60, 10.00, 5.12, 22.38, 33…
#> $ revol_bal           <int> 21936, 5457, 23453, 31740, 2284, 2016, 14330, 2758…
#> $ inq_last_12m        <int> 3, 1, 0, 0, 3, 5, 0, 1, 8, 1, 0, 12, 4, 8, 1, 3, 0…
#> $ delinq_2yrs         <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,…
#> $ home_ownership      <fct> MORTGAGE, RENT, OWN, MORTGAGE, MORTGAGE, MORTGAGE,…
#> $ not_paid            <fct> 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1,…
#> $ log_inc             <dbl> 11.962088, 10.819778, 10.126631, 12.072541, 11.608…
#> $ verified            <fct> 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ grdCtoA             <fct> 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,…

Cek missing value

anyNA(loan_clean)
#> [1] FALSE

Cek persebaran/pattern data

# explore with summary

summary(loan_clean)
#>  initial_list_status               purpose       int_rate      installment     
#>  f: 323              credit_card       :321   Min.   : 5.32   Min.   :  31.04  
#>  w:1233              debt_consolidation:996   1st Qu.:10.42   1st Qu.: 247.26  
#>                      home_improvement  :161   Median :14.08   Median : 381.23  
#>                      major_purchase    : 51   Mean   :14.71   Mean   : 468.72  
#>                      small_business    : 27   3rd Qu.:18.06   3rd Qu.: 635.86  
#>                                               Max.   :30.94   Max.   :1503.89  
#>       dti           revol_bal       inq_last_12m     delinq_2yrs    
#>  Min.   :  0.00   Min.   :     0   Min.   : 0.000   Min.   :0.0000  
#>  1st Qu.: 11.00   1st Qu.:  4974   1st Qu.: 1.000   1st Qu.:0.0000  
#>  Median : 17.18   Median : 10163   Median : 2.000   Median :0.0000  
#>  Mean   : 18.86   Mean   : 15258   Mean   : 2.359   Mean   :0.3059  
#>  3rd Qu.: 24.53   3rd Qu.: 18094   3rd Qu.: 3.000   3rd Qu.:0.0000  
#>  Max.   :198.56   Max.   :258897   Max.   :18.000   Max.   :8.0000  
#>   home_ownership not_paid    log_inc       verified grdCtoA
#>  MORTGAGE:776    0:778    Min.   : 7.824   0:578    0:961  
#>  OWN     :201    1:778    1st Qu.:10.779   1:978    1:595  
#>  RENT    :579             Median :11.127                   
#>                           Mean   :11.123                   
#>                           3rd Qu.:11.469                   
#>                           Max.   :13.998

Literature: Higher debt-to-income ratio (dti) and amount of credit card debts are both associated with a greater likelihood of loan defaults.

# numeric predictor vs target variable
boxplot(loan_clean$dti , loan_clean$not_paid, ylim = c(0,10))

table(loan_clean$grdCtoA, loan_clean$not_paid)
#>    
#>       0   1
#>   0 422 539
#>   1 356 239

Cek class-imbalance

prop.table(table(loan_clean$not_paid))
#> 
#>   0   1 
#> 0.5 0.5

Proporsi yang seimbang penting agar model dapat mempelajari karakteristik kelas positif maupun negatif secara seimbang, tidak belajar dari satu kelas saja. Hal ini mencegah model dari hanya baik memprediksi 1 kelas saja.

Proporsi yang imbalance umumnya 90/10 atau 95/5.

Kalau datanya imbalance:

  • tambah data
  • downSampling -> buang observasi dari kelas mayoritas, sehingga seimbang
  • upSampling -> duplicate observasi dari kelas minoritas, sehingga seimbang

Cross Validation

  • split data menjadi 2 bagian yaitu data train dan data test.
  • data train akan digunakan untuk training model.
  • data test akan digunakan untuk pengujian performa model. model akan diuji untuk memprediksi data test. hasil prediksi dan data aktual dari data test akan dibandingkan untuk validasi performa model.

Analogi:

  • 100 soal
  • 80 soal saya pakai untuk belajar (data train)
  • 20 soal saya pakai untuk ujian (data test)

tujuan dari cross validation adalah untuk mengetahui seberapa baik model yg sudah kita buat.

# # intuisi set seed: mengunci random number kita

set.seed(417) # pakai set.seed -> random number dikunci, hasil sampling selalu sama
sample(c("Angela", "Anthony", "Rizky", "Kevin"), 2)
#> [1] "Kevin"  "Angela"
RNGkind(sample.kind = "Rounding") # tambahan khusus u/ R 3.6 ke atas 
set.seed(417) # mengunci random number yang dipilih

# index sampling
index <- sample(x = nrow(loan_clean) , size = nrow(loan_clean) * 0.8)

# splitting
loans.train <- loan_clean[index , ]
loans.test <- loan_clean[-index , ]

NOTE: Proporsi 0.8/0.2 tidak mutlak, tergantung kebutuhan kita. Umumnya yang lebih banyak adalah untuk data train.

# re-check class imbalance
prop.table(table(loans.train$not_paid))
#> 
#>         0         1 
#> 0.4959807 0.5040193
prop.table(table(loans.test$not_paid))
#> 
#>         0         1 
#> 0.5160256 0.4839744

proporsi kelas yang balance penting untuk data train karena kita akan melatih model menggunakan data train.

library(rsample)

set.seed(417)

index_1 <- initial_split(data = loan_clean, prop = 0.8 , strata ="not_paid")

loan_train <- training(index_1)
loan_test <- testing(index_1)
  
prop.table(table(loan_train$not_paid))
#> 
#>   0   1 
#> 0.5 0.5
prop.table(table(loan_test$not_paid))
#> 
#>   0   1 
#> 0.5 0.5

Build Model

Buatlah model logistic regression untuk memprediksi status loan (not_paid). Silahkan lakukan feature selection berdasarkan pertimbangan bisnis atau/dan statistik!

pakai data loans.train

model.loans <- glm(not_paid ~ installment + delinq_2yrs +  home_ownership + log_inc + verified + grdCtoA , data = loans.train, family = "binomial" )
summary(model.loans)
#> 
#> Call:
#> glm(formula = not_paid ~ installment + delinq_2yrs + home_ownership + 
#>     log_inc + verified + grdCtoA, family = "binomial", data = loans.train)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.0236  -1.1157   0.6261   1.1203   1.6859  
#> 
#> Coefficients:
#>                      Estimate Std. Error z value    Pr(>|z|)    
#> (Intercept)         3.1346656  1.1955650   2.622    0.008744 ** 
#> installment         0.0011201  0.0002211   5.066 0.000000406 ***
#> delinq_2yrs         0.2564169  0.0850741   3.014    0.002578 ** 
#> home_ownershipOWN   0.3857017  0.1875300   2.057    0.039711 *  
#> home_ownershipRENT  0.1647821  0.1296008   1.271    0.203566    
#> log_inc            -0.3413549  0.1092992  -3.123    0.001789 ** 
#> verified1           0.2274757  0.1250104   1.820    0.068812 .  
#> grdCtoA1           -0.4525553  0.1240481  -3.648    0.000264 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 1724.5  on 1243  degrees of freedom
#> Residual deviance: 1644.5  on 1236  degrees of freedom
#> AIC: 1660.5
#> 
#> Number of Fisher Scoring iterations: 4

Pilih masing-masing 1 untuk prediktor kategorik dan prediktor numerik, kemudian interpretasikan:

# numerik:
exp(-0.3413549)
#> [1] 0.7108066
# kategorik:

exp(-0.4525553)
#> [1] 0.6360009

setiap peningkatan 1 point dari log_inc debitur akan 0.71 kali lebih mungkin untuk not paid kelas 1 (gagal bayar).

untuk debitur yang memiliki grad B dan A 0.64 kali lebih mungkin untuk notpaid kelas 1 (gagal bayar)

Predict

predict(model, newdata, type)

pada type terdapat pilihan:

  • link: menghasilkan log of odds
  • response: menghasilkan peluang

Prediksi log of odds not_paid untuk 6 data teratas:

predict(object = model.loans , 
        newdata = loans.test[1:6,], 
        type = "link")
#>           6           8           9          10          26          38 
#> -0.85787892 -0.04013106  0.09345395  0.11022883  0.55817539  0.57478151

Prediksi probability not_paid untuk 6 data teratas:

predict(object = model.loans , 
        newdata = loans.test[1:6,], 
        type = "response")
#>         6         8         9        10        26        38 
#> 0.2977827 0.4899686 0.5233465 0.5275293 0.6360303 0.6398658

Lakukan prediksi probability not_paid untuk data loans.test dan disimpan pada kolom baru bernama pred.Risk.

loans.test$pred.Risk <- predict(object = model.loans , 
        newdata = loans.test, 
        type = "response")

Klasifikasikan data loans.test berdasarkan pred.Risk dan simpan pada kolom baru bernama pred.Label.

# ifelse(kondisi, benar, salah)

loans.test$pred.Label <- ifelse(loans.test$pred.Risk > 0.5 , yes = 1 , no = 0)

# pastikan kelas target (aktual dan prediksi) bertipe factor

loans.test$pred.Label <- as.factor(loans.test$pred.Label)

head(loans.test)

Note:

Penentuan label yang menjadi angka 1 pada model logistic regression adalah berdasarkan levels.

kelas “0”, “1” -> basis = 0, peluang mendekati 0 -> 0 peluang mendekati 1 -> 1

kelas “honors” “non-honors” -> basis = honors peluang mendekati 0 -> honors peluang mendekati 1 -> non-honors

# lihat hasil prediksi
loans.test %>% 
  select(not_paid, pred.Risk, pred.Label) %>% 
  head(6)

Notes

  1. Seleksi model logistic regression:
  • AIC : aic yang paling kecil
  • Perfect separator -> kolom yang sangat menentukan hasil prediksi, tidak bersinggungan
  1. Asumsi model logistic regression:
  • linearity logit
  • multicol
  • independence of observaton
  1. Workflow klasifikasi:
  • business question
  • read data
  • data wrangling
  • EDA
  • Cross validation
  • Data Preprocessing
  • build model
  • predict
  • evaluasi
  • tuning model

Model Evaluation

Setelah dilakukan prediksi menggunakan model, masih ada saja prediksi yang salah. Pada klasifikasi, kita mengevaluasi model berdasarkan confusion matrix:

  • true positive (TP): prediksi positif; aktual positif
  • true negative (TN): prediksi negatif; aktual negatif
  • false positive (FP): prediksi positif; aktual negatif
  • false negative (FN): prediksi negatif; aktual positif
library(caret)

confusionMatrix(data = loans.test$pred.Label, reference = loans.test$not_paid , positive = "1")
#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction  0  1
#>          0 93 56
#>          1 68 95
#>                                           
#>                Accuracy : 0.6026          
#>                  95% CI : (0.5459, 0.6573)
#>     No Information Rate : 0.516           
#>     P-Value [Acc > NIR] : 0.001291        
#>                                           
#>                   Kappa : 0.2063          
#>                                           
#>  Mcnemar's Test P-Value : 0.323236        
#>                                           
#>             Sensitivity : 0.6291          
#>             Specificity : 0.5776          
#>          Pos Pred Value : 0.5828          
#>          Neg Pred Value : 0.6242          
#>              Prevalence : 0.4840          
#>          Detection Rate : 0.3045          
#>    Detection Prevalence : 0.5224          
#>       Balanced Accuracy : 0.6034          
#>                                           
#>        'Positive' Class : 1               
#>

4 metrics performa model: Accuracy, Sensitivity/Recall, Precision, Specificity

Accuracy

Seberapa banyak yang benar diprediksi dari keseluruhan data (positif maupun negatif).

TP+TN/TOTAL

(95 + 93) / nrow(loans.test)
#> [1] 0.6025641

Digunakan ketika:

  • kelas target sama penting
  • data balance

Ada kondisi ketika accuracy bukanlah metrics terpenting. Umumnya ketika:

  • kita mementingkan salah satu kelas (misal, kelas positif)
  • data kita imbalance

Saat kita mementingkan kelas tertentu (kelas positif), maka kita dapat memilih antara menggunakan metrics Recall / Precision:

Sensitivity/Recall

Seberapa banyak yang benar diprediksi positif, dari yang reality-nya (aktualnya) positif.

TP/(TP+FN)

95 / (95 + 56)
#> [1] 0.6291391

Pos Pred Value/Precision

Seberapa banyak yang benar diprediksi positif, dari yang diprediksi positif.

TP/(TP+FP)

95/ (95 + 68)
#> [1] 0.5828221

Specificity

Seberapa banyak yang tepat diprediksi negatif, dari yang reality-nya negatif. Jarang dipakai karena kita tidak sering fokus pada kelas negatif.

TN/(TN+FP)

93/(93+68)
#> [1] 0.5776398

Metrics yang baik untuk kasus loans adalah? (positif: not paid)

recall

Bila hasil evaluasi (nilai metrics) belum memuaskan, dapat dilakukan Model Tuning:

  1. ganti prediktor
  2. ubah data pre-processingnya: misal upsample/downsample, scaling dll.
  3. ganti modelnya (pakai model yang lebih robust)
  4. ganti treshold prediction (tidak terlalu dianjurkan, karena bisa memaksakan):
  • geser mendekati 0: meningkatkan recall
  • geser mendekati 1: meningkatkan precision
# label baru dari threshold baru
loans.test$pred.Label <- ifelse(loans.test$pred.Risk > 0.4 , yes = 1,no =  0)
loans.test$pred.Label <- as.factor(loans.test$pred.Label)

# confusion matrix
confusionMatrix(data =loans.test$pred.Label  , reference = loans.test$not_paid, positive = "1")
#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction   0   1
#>          0  42  22
#>          1 119 129
#>                                                
#>                Accuracy : 0.5481               
#>                  95% CI : (0.491, 0.6042)      
#>     No Information Rate : 0.516                
#>     P-Value [Acc > NIR] : 0.1409               
#>                                                
#>                   Kappa : 0.1129               
#>                                                
#>  Mcnemar's Test P-Value : 0.0000000000000006234
#>                                                
#>             Sensitivity : 0.8543               
#>             Specificity : 0.2609               
#>          Pos Pred Value : 0.5202               
#>          Neg Pred Value : 0.6562               
#>              Prevalence : 0.4840               
#>          Detection Rate : 0.4135               
#>    Detection Prevalence : 0.7949               
#>       Balanced Accuracy : 0.5576               
#>                                                
#>        'Positive' Class : 1                    
#>

k-NN

k-NN adalah K-nearest neighboor. Metode ini akan mengkasifikasi data baru dengan membandingkan karakteristik data baru (data test) dengan data yang ada (data train). Kedekatan karakteristik tersebut diukur dengan Euclidean Distance hingga didapatkan jarak. Kemudian akan dipilih k tetangga terdekat dari data baru tersebut, kemudian ditentukan kelasnya menggunakan majority voting.

Picking Optimum k

  • jangan terlalu besar: pemilihan kelas hanya berdasarkan kelas yang dominan dan mengabaikan pola kecil yang ternyata penting.
  • jangan terlalu kecil: rentan mengklasifikasikan data baru ke kelas outlier.
  • k optimum adalah akar dari jumlah data kita: sqrt(nrow(data))
  • untuk menghindari seri ketika majority voting:
    • k harus ganjil bila jumlah kelas target genap
    • k harus genap bila jumlah kelas target ganjil
    • k tidak boleh angka kelipatan jumlah kelas target
  • bila hasil majority voting seri, maka kelas akan dipilih secara random.

Karakteristik k-NN

  • tidak membuat model: langsung mengklasifikasi saat itu juga, tidak belajar dari data, setiap ingin mengklasifikasi harus menyediakan data train lagi.
  • tidak ada asumsi
  • dapat memprediksi multiclass
  • baik untuk prediktor numerik (karena mengklasifikasikan berdasarkan jarak), tidak baik untuk prediktor kategorik
  • robust
  • tidak interpretable

Breast Cancer Prediction

Business Question

Kanker payudara adalah kanker yang paling umum menyerang wanita di dunia. Kanker payudara dapat berupa kanker jinak (benign) atau sudah ganas (malignant). Kanker ganas dapat menyebar ke organ-organ tubuh lainnya. Ingin dibuat model prediksi untuk memprediksi apakah kanker masih jinak (benign) atau sudah ganas (malignant).

Read Data

wbcd <- read.csv("data_input/wisc_bc_data.csv")

wbcd_clean <- wbcd[,-1]
head(wbcd_clean)
# inspect data
str(wbcd_clean)
#> 'data.frame':    569 obs. of  31 variables:
#>  $ diagnosis        : chr  "B" "B" "B" "B" ...
#>  $ radius_mean      : num  12.3 10.6 11 11.3 15.2 ...
#>  $ texture_mean     : num  12.4 18.9 16.8 13.4 13.2 ...
#>  $ perimeter_mean   : num  78.8 69.3 70.9 73 97.7 ...
#>  $ area_mean        : num  464 346 373 385 712 ...
#>  $ smoothness_mean  : num  0.1028 0.0969 0.1077 0.1164 0.0796 ...
#>  $ compactness_mean : num  0.0698 0.1147 0.078 0.1136 0.0693 ...
#>  $ concavity_mean   : num  0.0399 0.0639 0.0305 0.0464 0.0339 ...
#>  $ points_mean      : num  0.037 0.0264 0.0248 0.048 0.0266 ...
#>  $ symmetry_mean    : num  0.196 0.192 0.171 0.177 0.172 ...
#>  $ dimension_mean   : num  0.0595 0.0649 0.0634 0.0607 0.0554 ...
#>  $ radius_se        : num  0.236 0.451 0.197 0.338 0.178 ...
#>  $ texture_se       : num  0.666 1.197 1.387 1.343 0.412 ...
#>  $ perimeter_se     : num  1.67 3.43 1.34 1.85 1.34 ...
#>  $ area_se          : num  17.4 27.1 13.5 26.3 17.7 ...
#>  $ smoothness_se    : num  0.00805 0.00747 0.00516 0.01127 0.00501 ...
#>  $ compactness_se   : num  0.0118 0.03581 0.00936 0.03498 0.01485 ...
#>  $ concavity_se     : num  0.0168 0.0335 0.0106 0.0219 0.0155 ...
#>  $ points_se        : num  0.01241 0.01365 0.00748 0.01965 0.00915 ...
#>  $ symmetry_se      : num  0.0192 0.035 0.0172 0.0158 0.0165 ...
#>  $ dimension_se     : num  0.00225 0.00332 0.0022 0.00344 0.00177 ...
#>  $ radius_worst     : num  13.5 11.9 12.4 11.9 16.2 ...
#>  $ texture_worst    : num  15.6 22.9 26.4 15.8 15.7 ...
#>  $ perimeter_worst  : num  87 78.3 79.9 76.5 104.5 ...
#>  $ area_worst       : num  549 425 471 434 819 ...
#>  $ smoothness_worst : num  0.139 0.121 0.137 0.137 0.113 ...
#>  $ compactness_worst: num  0.127 0.252 0.148 0.182 0.174 ...
#>  $ concavity_worst  : num  0.1242 0.1916 0.1067 0.0867 0.1362 ...
#>  $ points_worst     : num  0.0939 0.0793 0.0743 0.0861 0.0818 ...
#>  $ symmetry_worst   : num  0.283 0.294 0.3 0.21 0.249 ...
#>  $ dimension_worst  : num  0.0677 0.0759 0.0788 0.0678 0.0677 ...

Variabel target: diagnosis

B = benign M = Malignant

Data Wrangling

wbcd_clean$diagnosis <- as.factor(wbcd_clean$diagnosis)

Exploratory Data Analysis

  • cek proporsi kelas
prop.table(table(wbcd_clean$diagnosis))
#> 
#>         B         M 
#> 0.6274165 0.3725835

Insight: data masih cukup balance

  • cek range nilai tiap variable prediktor: range harus sama karena knn mengklasifikasikan berdasarkan jarak. kalau ada nilai yang tinggi sendiri dibanding yang lain, maka variable tersebut akan sangat mempengaruhi hasil klasifikasi dan mengabaikan variable yang lain.
# cek range nilai tiap variable
summary(wbcd_clean)
#>  diagnosis  radius_mean      texture_mean   perimeter_mean     area_mean     
#>  B:357     Min.   : 6.981   Min.   : 9.71   Min.   : 43.79   Min.   : 143.5  
#>  M:212     1st Qu.:11.700   1st Qu.:16.17   1st Qu.: 75.17   1st Qu.: 420.3  
#>            Median :13.370   Median :18.84   Median : 86.24   Median : 551.1  
#>            Mean   :14.127   Mean   :19.29   Mean   : 91.97   Mean   : 654.9  
#>            3rd Qu.:15.780   3rd Qu.:21.80   3rd Qu.:104.10   3rd Qu.: 782.7  
#>            Max.   :28.110   Max.   :39.28   Max.   :188.50   Max.   :2501.0  
#>  smoothness_mean   compactness_mean  concavity_mean     points_mean     
#>  Min.   :0.05263   Min.   :0.01938   Min.   :0.00000   Min.   :0.00000  
#>  1st Qu.:0.08637   1st Qu.:0.06492   1st Qu.:0.02956   1st Qu.:0.02031  
#>  Median :0.09587   Median :0.09263   Median :0.06154   Median :0.03350  
#>  Mean   :0.09636   Mean   :0.10434   Mean   :0.08880   Mean   :0.04892  
#>  3rd Qu.:0.10530   3rd Qu.:0.13040   3rd Qu.:0.13070   3rd Qu.:0.07400  
#>  Max.   :0.16340   Max.   :0.34540   Max.   :0.42680   Max.   :0.20120  
#>  symmetry_mean    dimension_mean      radius_se        texture_se    
#>  Min.   :0.1060   Min.   :0.04996   Min.   :0.1115   Min.   :0.3602  
#>  1st Qu.:0.1619   1st Qu.:0.05770   1st Qu.:0.2324   1st Qu.:0.8339  
#>  Median :0.1792   Median :0.06154   Median :0.3242   Median :1.1080  
#>  Mean   :0.1812   Mean   :0.06280   Mean   :0.4052   Mean   :1.2169  
#>  3rd Qu.:0.1957   3rd Qu.:0.06612   3rd Qu.:0.4789   3rd Qu.:1.4740  
#>  Max.   :0.3040   Max.   :0.09744   Max.   :2.8730   Max.   :4.8850  
#>   perimeter_se       area_se        smoothness_se      compactness_se    
#>  Min.   : 0.757   Min.   :  6.802   Min.   :0.001713   Min.   :0.002252  
#>  1st Qu.: 1.606   1st Qu.: 17.850   1st Qu.:0.005169   1st Qu.:0.013080  
#>  Median : 2.287   Median : 24.530   Median :0.006380   Median :0.020450  
#>  Mean   : 2.866   Mean   : 40.337   Mean   :0.007041   Mean   :0.025478  
#>  3rd Qu.: 3.357   3rd Qu.: 45.190   3rd Qu.:0.008146   3rd Qu.:0.032450  
#>  Max.   :21.980   Max.   :542.200   Max.   :0.031130   Max.   :0.135400  
#>   concavity_se       points_se         symmetry_se        dimension_se      
#>  Min.   :0.00000   Min.   :0.000000   Min.   :0.007882   Min.   :0.0008948  
#>  1st Qu.:0.01509   1st Qu.:0.007638   1st Qu.:0.015160   1st Qu.:0.0022480  
#>  Median :0.02589   Median :0.010930   Median :0.018730   Median :0.0031870  
#>  Mean   :0.03189   Mean   :0.011796   Mean   :0.020542   Mean   :0.0037949  
#>  3rd Qu.:0.04205   3rd Qu.:0.014710   3rd Qu.:0.023480   3rd Qu.:0.0045580  
#>  Max.   :0.39600   Max.   :0.052790   Max.   :0.078950   Max.   :0.0298400  
#>   radius_worst   texture_worst   perimeter_worst    area_worst    
#>  Min.   : 7.93   Min.   :12.02   Min.   : 50.41   Min.   : 185.2  
#>  1st Qu.:13.01   1st Qu.:21.08   1st Qu.: 84.11   1st Qu.: 515.3  
#>  Median :14.97   Median :25.41   Median : 97.66   Median : 686.5  
#>  Mean   :16.27   Mean   :25.68   Mean   :107.26   Mean   : 880.6  
#>  3rd Qu.:18.79   3rd Qu.:29.72   3rd Qu.:125.40   3rd Qu.:1084.0  
#>  Max.   :36.04   Max.   :49.54   Max.   :251.20   Max.   :4254.0  
#>  smoothness_worst  compactness_worst concavity_worst   points_worst    
#>  Min.   :0.07117   Min.   :0.02729   Min.   :0.0000   Min.   :0.00000  
#>  1st Qu.:0.11660   1st Qu.:0.14720   1st Qu.:0.1145   1st Qu.:0.06493  
#>  Median :0.13130   Median :0.21190   Median :0.2267   Median :0.09993  
#>  Mean   :0.13237   Mean   :0.25427   Mean   :0.2722   Mean   :0.11461  
#>  3rd Qu.:0.14600   3rd Qu.:0.33910   3rd Qu.:0.3829   3rd Qu.:0.16140  
#>  Max.   :0.22260   Max.   :1.05800   Max.   :1.2520   Max.   :0.29100  
#>  symmetry_worst   dimension_worst  
#>  Min.   :0.1565   Min.   :0.05504  
#>  1st Qu.:0.2504   1st Qu.:0.07146  
#>  Median :0.2822   Median :0.08004  
#>  Mean   :0.2901   Mean   :0.08395  
#>  3rd Qu.:0.3179   3rd Qu.:0.09208  
#>  Max.   :0.6638   Max.   :0.20750

Range tiap variabel berbeda sehingga perlu dilakukan feature rescaling di tahap data pre-processing.

Cross Validation

RNGkind(sample.kind = "Rounding")
set.seed(417)

index <- sample(nrow(wbcd_clean) , nrow(wbcd_clean)*0.8)

wbcd_train <- wbcd_clean[index,]
wbcd_test <- wbcd_clean[-index,]
# recheck class balance
prop.table(table(wbcd_train$diagnosis))
#> 
#>         B         M 
#> 0.6175824 0.3824176

Insight: proporsi masih termasuk balance

Data Pre-processing

Scaling

Scaling bisa menggunakan min-max normalization atau z-score standarization.

  1. Min-Max Normalization

Rumus = x-min(x) / max(x)-min(x)

normalize <- function(x){
  return ( 
    (x - min(x))/(max(x) - min(x)) 
           )
}
# contoh:
normalize(c(1,2,3,4,5)) # memampatkan range nilai menjadi 0-1
#> [1] 0.00 0.25 0.50 0.75 1.00
normalize(c(100,7600,3000,4000,5000))
#> [1] 0.0000000 1.0000000 0.3866667 0.5200000 0.6533333

-> digunakan ketika tau angka pasti min dan max nya. misalnya nilai ujian matematika pasti nilai min-max nya 0 - 100.

  1. Z-score Standarization: dapat menggunakan function scale()

Rumus = x-mean(x) / sd(x)

# contoh:
scale(c(1,2,3,4,5)) # data kita seberapa menyimpang (sd) dari pusatnya (mean)
#>            [,1]
#> [1,] -1.2649111
#> [2,] -0.6324555
#> [3,]  0.0000000
#> [4,]  0.6324555
#> [5,]  1.2649111
#> attr(,"scaled:center")
#> [1] 3
#> attr(,"scaled:scale")
#> [1] 1.581139
mean(c(1,2,3,4,5))
#> [1] 3
sd(c(1,2,3,4,5))
#> [1] 1.581139
(1 - 3)/1.581139
#> [1] -1.264911

-> digunakan ketika tidak diketahui angka min dan max pastinya. misalnya temperature bisa dari kisaran -inf s.d +inf

Untuk k-NN, dipisahkan antara prediktor dan label (target variabelnya).

# prediktor
wbcd_train_x <- wbcd_train[,-1]

wbcd_test_x <- wbcd_test[,-1]

# target
wbcd_train_y <- wbcd_train$diagnosis

wbcd_test_y <- wbcd_test$diagnosis

Data prediktor akan discaling menggunakan z-score standarization. Data test juga harus discaling menggunakan parameter dari data train (karena menganggap data test adalah unseen data).

# scaling data prediktor
wbcd_train_xs <- scale(wbcd_train_x)
wbcd_test_xs <- scale(wbcd_test_x , 
                      center = attr(wbcd_train_xs,"scaled:center") , 
                      scale = attr(wbcd_train_xs,"scaled:scale") )

k-NN tidak membuat model sehingga langsung ke predict.

Predict

# find optimum k
sqrt(nrow(wbcd_train_xs))
#> [1] 21.33073
  • jumlah kelas target: 2
  • k: 19 21 23
library(class) # package untuk fungsi `knn()`

wbcd_pred <- knn(train = wbcd_train_xs , 
                 test = wbcd_test_xs , 
                 cl = wbcd_train_y, 
                 k = 23
                )
# cek hasil prediksi
wbcd_pred
#>   [1] B B B B B B B B B M B B B B B B M B B B B B B B B B B M M B B B B B B B M
#>  [38] B B M M B B B B M M M B B M B B B B M B M B M B M B M B B M B M M B B M B
#>  [75] M B B B M B B M B M M B M B B B B B M B B B M M B M B M B M M B B B M B B
#> [112] M B B
#> Levels: B M

Model evaluation

# confusion matrix
library(caret)

confusionMatrix(data = wbcd_pred, reference = wbcd_test_y, positive = "M")
#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction  B  M
#>          B 76  3
#>          M  0 35
#>                                                
#>                Accuracy : 0.9737               
#>                  95% CI : (0.925, 0.9945)      
#>     No Information Rate : 0.6667               
#>     P-Value [Acc > NIR] : 0.0000000000000002673
#>                                                
#>                   Kappa : 0.9396               
#>                                                
#>  Mcnemar's Test P-Value : 0.2482               
#>                                                
#>             Sensitivity : 0.9211               
#>             Specificity : 1.0000               
#>          Pos Pred Value : 1.0000               
#>          Neg Pred Value : 0.9620               
#>              Prevalence : 0.3333               
#>          Detection Rate : 0.3070               
#>    Detection Prevalence : 0.3070               
#>       Balanced Accuracy : 0.9605               
#>                                                
#>        'Positive' Class : M                    
#> 
pairs(wbcd_clean[,2:6])

k-NN:

k-NN mengkasifikasi berdasarkan jarak - menghitung jarak ke semua datanya - ditentukan berapa k tetangga terdekat (data terdekat) - dilakukan majority voting dari kelas k tetangga terdekatnya