Regression Models - Graduate Admission

Intro

About Data

Data wholesale.csv diperoleh dari: https://archive.ics.uci.edu/ml/datasets/wholesale+customers, mencatat penjualan beberapa kategori produk kepada distributor. Kali ini, kita akan melakukan prediksi channel dari tiap observasi data wholesale untuk melihat variabel apa saja yang mempengaruhi sebuah observasi bisa digolongkan ke dalam channel tertentu. Ketika variabel target kelompok-kelompok tertentu, maka algoritma yang akan digunakan dalam prediksi data ini adalah logistic regression dan k-nearest neighbors (KNN).

Steps

Sebagai acuan, berikut ini beberapa hal yang perlu diperhatikan selam proses pemodelan prediksi menggunakan logistic regression dan k-nearest neighbors:

• Variabel target : Channel
• Prediktor : variabel-variabel informatif lain selain target

Step by step:

1. Preparation & Exploratory : select variabel yang diperlukan, cek NA, data outliers, inspect, transformasi tipe data, spliting/cross validation, dll
2. Modeling and Analysing
3. Model Evaluation : Confusion Matrix
4. Conclusion

Preparation & Exploratory

Library

# change scientific number
options(scipen = 99)
library(dplyr)       #wrangling
library(ggplot2)     #visualization
library(MASS)        #stepAIC
library(caret)       #confusionMatrix
library(class)       #knn

Input

Pertama, panggil data wholesale kemudian simpan sebagai sale.

sale <- read.csv("data_c1/wholesale.csv")
head(sale)

str(sale)

## 'data.frame':    440 obs. of  8 variables:
##  $ Channel         : int  2 2 2 1 2 2 2 2 1 2 ...
##  $ Region          : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ Fresh           : int  12669 7057 6353 13265 22615 9413 12126 7579 5963 6006 ...
##  $ Milk            : int  9656 9810 8808 1196 5410 8259 3199 4956 3648 11093 ...
##  $ Grocery         : int  7561 9568 7684 4221 7198 5126 6975 9426 6192 18881 ...
##  $ Frozen          : int  214 1762 2405 6404 3915 666 480 1669 425 1159 ...
##  $ Detergents_Paper: int  2674 3293 3516 507 1777 1795 3140 3321 1716 7425 ...
##  $ Delicassen      : int  1338 1776 7844 1788 5185 1451 545 2566 750 2098 ...

Keterangan variabel pada data:

1. FRESH: annual spending (m.u.) on fresh products (Continuous);
2. MILK: annual spending (m.u.) on milk products (Continuous);
3. GROCERY: annual spending (m.u.)on grocery products (Continuous);
4. FROZEN: annual spending (m.u.)on frozen products (Continuous)
5. DETERGENTS_PAPER: annual spending (m.u.) on detergents and paper products (Continuous)
6. DELICATESSEN: annual spending (m.u.)on and delicatessen products (Continuous);
7. CHANNEL: customer Channel - Horeca (Hotel/Restaurant/Cafe) or Retail channel (Nominal)
8. REGION: customer Region Lisnon, Oporto or Other (Nominal)

Any NA

# anyNA(sale)
colSums(is.na(sale))

##          Channel           Region            Fresh             Milk 
##                0                0                0                0 
##          Grocery           Frozen Detergents_Paper       Delicassen 
##                0                0                0                0

Data Transformation

Kita perlu melakukan transformasi tipe data pada variabel yang sebenarnya menyimpan data kategorik tapi tampak numerik, yaitu channel dan region. Kedua variabel ini dipahami R sebagai integer karena nilainya berupa angka, padahal angka tersebut adalah representasi dari suatu kelas. Misalnya channel terbagi menjadi kelas Horeca dan Retail, sedangkan region terdiri dari Lisbon, Oporto, dan Other Region.

sale <- sale %>% 
   mutate(Channel = factor(Channel, levels = c(1,2), labels = c("Horeca", "Retail")),
          Region = factor(Region, levels = c(1,2,3), labels = c("Lisbon", "Oporto", "Other")))

str(sale)

## 'data.frame':    440 obs. of  8 variables:
##  $ Channel         : Factor w/ 2 levels "Horeca","Retail": 2 2 2 1 2 2 2 2 1 2 ...
##  $ Region          : Factor w/ 3 levels "Lisbon","Oporto",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Fresh           : int  12669 7057 6353 13265 22615 9413 12126 7579 5963 6006 ...
##  $ Milk            : int  9656 9810 8808 1196 5410 8259 3199 4956 3648 11093 ...
##  $ Grocery         : int  7561 9568 7684 4221 7198 5126 6975 9426 6192 18881 ...
##  $ Frozen          : int  214 1762 2405 6404 3915 666 480 1669 425 1159 ...
##  $ Detergents_Paper: int  2674 3293 3516 507 1777 1795 3140 3321 1716 7425 ...
##  $ Delicassen      : int  1338 1776 7844 1788 5185 1451 545 2566 750 2098 ...

Target Proportion

Sebelum melakukan cross validation atau splitting data menjadi data train dan test, kita perlu melihat proporsi dari target variabel channel. Proporsi ini dilakukan supaya hasil modeling dan prediksi dari tiap kelas punya peluang yang sama/tidak timpang.

table(sale$Channel)

## 
## Horeca Retail 
##    298    142

prop.table(table(sale$Channel))

## 
##    Horeca    Retail 
## 0.6772727 0.3227273

Proporsi kedua kelas tidak seimbang, dimana Horeca lebih dominan dibanding Retail. Untuk itu, nantinya akan dilakukan downsampling pada data train agar proporsi kelas dari variabel target lebih seimbang. Proporsi yang seimbang baik bagi model agar bisa menghasilkan prediksi yang lebih akurat.

Cross Validation

Splitting data sale menjadi train dan test dengan bobot 70:30. Bobot data train lebih besar karena data ini akan digunakan dalam proses modeling nanti, tujuannya supaya model kita bisa mempelajari pola data yang beragam.

set.seed(123)
idx <- sample(nrow(sale), nrow(sale)*0.7)
train <- sale[idx,]
test <- sale[-idx,]

prop.table(table(train$Channel))

## 
##    Horeca    Retail 
## 0.6915584 0.3084416

Berhubung proporsi kelas target tidak seimbang, kita akan melakukan downsampling yakni mengurangi data pada kelas dominan secara acak agar menghasilkan proporsi kedua kelas yang imbang. Berikut beberapa parameter yang perlu diperhatikan dalam menggunakan fungsi downSample() - x : variabel prediktor data train - y : variabel target data train - yname : nama kolom variabel target dapa data train

train.ds <- downSample(x = train[, -1], y = train[, 1], yname = "Channel")
prop.table(table(train.ds$Channel))

## 
## Horeca Retail 
##    0.5    0.5

Sekarang kita sudah punya proporsi kelas target yang seimbang. Waktunya membuat model!

Logistic Regression Model

model1.log <- glm(Channel ~ ., train.ds, family = "binomial")
summary(model1.log)

## 
## Call:
## glm(formula = Channel ~ ., family = "binomial", data = train.ds)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.17085  -0.30867  -0.05534   0.09805   2.50080  
## 
## Coefficients:
##                     Estimate  Std. Error z value   Pr(>|z|)    
## (Intercept)      -4.75021311  1.13570370  -4.183 0.00002882 ***
## RegionOporto      1.82583897  1.32096592   1.382     0.1669    
## RegionOther       1.10833975  1.03672078   1.069     0.2850    
## Fresh            -0.00002633  0.00002494  -1.056     0.2911    
## Milk             -0.00013018  0.00011606  -1.122     0.2620    
## Grocery           0.00024823  0.00010195   2.435     0.0149 *  
## Frozen           -0.00003555  0.00005395  -0.659     0.5099    
## Detergents_Paper  0.00106321  0.00023404   4.543 0.00000555 ***
## Delicassen        0.00029287  0.00027052   1.083     0.2790    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 263.396  on 189  degrees of freedom
## Residual deviance:  78.045  on 181  degrees of freedom
## AIC: 96.045
## 
## Number of Fisher Scoring iterations: 8

Significancy codes pada summary(model1.log) menandakan hanya ada sedikit prediktor yang berpengaruh terhadap target, bahkan beberapa di antaranya tidak terlalu signifikan. Nilai AIC pun mencapai 102.2.

Maka dari itu, kita akan coba melakukan modeling dengan fungsi stepAIC() dari library(MASS), dengan sumber datanya dari model yang sudah kita buat.

model2.log <- stepAIC(model1.log, direction = "backward")

## Start:  AIC=96.04
## Channel ~ Region + Fresh + Milk + Grocery + Frozen + Detergents_Paper + 
##     Delicassen
## 
##                    Df Deviance     AIC
## - Region            2   80.154  94.154
## - Frozen            1   78.581  94.581
## - Fresh             1   79.172  95.172
## - Delicassen        1   79.270  95.270
## - Milk              1   79.417  95.417
## <none>                  78.045  96.045
## - Grocery           1   85.467 101.467
## - Detergents_Paper  1  101.092 117.092
## 
## Step:  AIC=94.15
## Channel ~ Fresh + Milk + Grocery + Frozen + Detergents_Paper + 
##     Delicassen
## 
##                    Df Deviance     AIC
## - Frozen            1   80.600  92.600
## - Fresh             1   81.422  93.422
## - Delicassen        1   81.629  93.629
## - Milk              1   82.092  94.092
## <none>                  80.154  94.154
## - Grocery           1   88.593 100.593
## - Detergents_Paper  1  103.766 115.766
## 
## Step:  AIC=92.6
## Channel ~ Fresh + Milk + Grocery + Detergents_Paper + Delicassen
## 
##                    Df Deviance     AIC
## - Delicassen        1   81.825  91.825
## - Fresh             1   82.101  92.101
## <none>                  80.600  92.600
## - Milk              1   83.139  93.139
## - Grocery           1   89.240  99.240
## - Detergents_Paper  1  108.418 118.418
## 
## Step:  AIC=91.82
## Channel ~ Fresh + Milk + Grocery + Detergents_Paper
## 
##                    Df Deviance     AIC
## - Fresh             1   82.670  90.670
## - Milk              1   83.337  91.337
## <none>                  81.825  91.825
## - Grocery           1   89.935  97.935
## - Detergents_Paper  1  108.426 116.426
## 
## Step:  AIC=90.67
## Channel ~ Milk + Grocery + Detergents_Paper
## 
##                    Df Deviance     AIC
## <none>                  82.670  90.670
## - Milk              1   84.941  90.941
## - Grocery           1   91.623  97.623
## - Detergents_Paper  1  110.775 116.775

stepAIC() berupaya mengembalikan model dengan prediktor yang dianggap signifikan berpengaruh terhadap terget, selain itu juga mencari nilai AIC (information loss) terkecil. Semakin kecil nilai AIC dari suatu model berarti informasi yang hilang dari model tersebut makin kecil juga. Dari modeling di atas, kita memperoleh AIC terkecil yaitu 96.94 dengan variabel-variabel prediktor yang dianggap berpengaruh yaitu Region, Milk, dan Detergents_Paper

summary(model2.log)

## 
## Call:
## glm(formula = Channel ~ Milk + Grocery + Detergents_Paper, family = "binomial", 
##     data = train.ds)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.27704  -0.33045  -0.08621   0.12976   2.40103  
## 
## Coefficients:
##                    Estimate Std. Error z value        Pr(>|z|)    
## (Intercept)      -3.7823925  0.5840656  -6.476 0.0000000000942 ***
## Milk             -0.0001213  0.0000793  -1.529         0.12628    
## Grocery           0.0002665  0.0001017   2.621         0.00878 ** 
## Detergents_Paper  0.0010211  0.0002101   4.860 0.0000011750205 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 263.40  on 189  degrees of freedom
## Residual deviance:  82.67  on 186  degrees of freedom
## AIC: 90.67
## 
## Number of Fisher Scoring iterations: 8

Prediction

Model regresi logistik yang sudah dibuat, kita uji untuk melakukan prediksi pada data test. Hasil prediksi akan mengembalikan nilai probability (kemungkinan) dari range 0-1, dari nilai tersebut kita perlu memberi label dengan batas treshold 0.5. Jika di atas 0.5 maka observasi termasuk kelas Retail, jika di bawah 0.5 termasuk kelas Horeca.

test$prob <- predict(model2.log, type = "response", newdata = test)
test$problabel <- factor(ifelse(test$prob > 0.5, "Retail", "Horeca"))

# perbandingan data asli (test) dengan hasil prediksi
test[1:10, c("Channel", "problabel")]

Secara visualisasi, dapat dilihat kecenderungan hasil prediksi:

ggplot(test, aes(x=prob)) +
  geom_density(lwd=0.5) +
  labs(title = "Distribusi Probabilitas") +
  theme_minimal()

Hasil visualisasi di atas menjelaskan nilai probabilitas dari prediksi terhadap channel. Kebanyakan nilainya mendekati 0 maupun 1, dan sedikit yang berada di sekitar nilai 0.5. Nilai probabilitas 0.5 ini nantinya akan menjadi batasan untuk mengklasifikasikan suatu observasi ke kelas tertentu.

Model Evaluation

Secara kasat mata, model kita mampu memprediksi data asli dengan cukup baik dimana pada sepuluh data pertama menghasilkan beberapa klasifikasi observasi secara tepat. Namun evaluasi kasus klasifikasi yang tepat haruslah dilakukan dengan menggunakan confusionMatrix() dari library(caret).

conf.log2 <- confusionMatrix(test$problabel, test$Channel, positive = "Retail")
conf.log2

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Horeca Retail
##     Horeca     79      9
##     Retail      6     38
##                                          
##                Accuracy : 0.8864         
##                  95% CI : (0.8195, 0.935)
##     No Information Rate : 0.6439         
##     P-Value [Acc > NIR] : 0.0000000002248
##                                          
##                   Kappa : 0.7486         
##                                          
##  Mcnemar's Test P-Value : 0.6056         
##                                          
##             Sensitivity : 0.8085         
##             Specificity : 0.9294         
##          Pos Pred Value : 0.8636         
##          Neg Pred Value : 0.8977         
##              Prevalence : 0.3561         
##          Detection Rate : 0.2879         
##    Detection Prevalence : 0.3333         
##       Balanced Accuracy : 0.8690         
##                                          
##        'Positive' Class : Retail         
## 

Hal yang menjadi highlight dari summary model di atas:

• Accuracy
  Seberapa baik model kita memprediksi kelas target secara global (positif maupun negatif). dipakai ketika kelas positif dan negatif sama pentingnya atau ketika proporsi kelas seimbang. Rumus : (TP+TN/TOTAL)
  
• Sensitivity/ Recall
  Ukuran kebaikan model terhadap kelas positif. Recall digunakan saat kita memperhatikan kelas positif dan kita ingin mengurangi resiko adanya salah prediksi untuk setiap kelas positif yang ada. Rumus : (TP/(TP+FN))
  
• Specificity
  Ukuran kebaikan model terhadap kelas negatif, seberapa banyak yang tepat diprediksi negatif, dari yang reality-nya negatif. Metrics ini tidak menjadi fokus karena umumnya kita fokus pada kelas positif. Rumus : (TN/(TN+FP))
  
• Pos Pred Value/Precision
  Seberapa presisi model memprediksi kelas positif. Precision digunakan saat kita ingin model yang kita buat benar-benar tepat menebak kelas positif / negatif dan tidak ingin ada miss-clasify pada kedua kelas tersebut. Rumus : (TP/(TP+FP))
  

Maka diperoleh nilai:

Accuracy <- round((38+79)/nrow(test),2)
Recall <- round((38)/(38+9),2)
Specificity <- round((79)/(79+6),2)
Precision <- round((38)/(38+6),2)

logreg.perf <- cbind.data.frame(Accuracy, Recall, Precision, Specificity)
logreg.perf

Interpretasi dari hasil evaluasi adalah model kita mampu memprediksi kelas target secara keseluruhan sebesar 89%. Dibanding kelas negatif, model mampu memprediksi benar untuk kelas positif (Retail) sebesar 81% dan sebalinya untuk kelas negatif (Horeca) adalah 93%. Kemudian secara ketepatan prediksi kelas positif maupun negatif model kita punya kebaikan sebesar 86%.

Dengan performa model secara keseluruhan di atas 80%, model2.log bisa dikatakan memiliki kemampuan yang cukup baik dalam memprediksi data dan menghasilkan kelas klasifikasi yang tepat.

K-Nearest Neighbors

k <- sqrt(nrow(sale))
k

## [1] 20.97618

# splitting x, y (tanpa variabel kategorik)
train_x <- train %>% 
  select_if(is.integer)
train_y <- train$Channel
test_x <- test %>% 
  select_if(is.integer)

#knn
knn <- knn(train = train_x, test = test_x, cl = train_y, k = k-1)

Model Evaluation

confusionMatrix(as.factor(knn), test$Channel, positive = "Retail")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Horeca Retail
##     Horeca     80     10
##     Retail      5     37
##                                          
##                Accuracy : 0.8864         
##                  95% CI : (0.8195, 0.935)
##     No Information Rate : 0.6439         
##     P-Value [Acc > NIR] : 0.0000000002248
##                                          
##                   Kappa : 0.7462         
##                                          
##  Mcnemar's Test P-Value : 0.3017         
##                                          
##             Sensitivity : 0.7872         
##             Specificity : 0.9412         
##          Pos Pred Value : 0.8810         
##          Neg Pred Value : 0.8889         
##              Prevalence : 0.3561         
##          Detection Rate : 0.2803         
##    Detection Prevalence : 0.3182         
##       Balanced Accuracy : 0.8642         
##                                          
##        'Positive' Class : Retail         
## 

Jika dibandingkan dengan regresi logistik, hasil prediksi kelas positif dan negatif pada target tidak terlalu berbeda. Namun jika dilihat dari accuracy-nya, model regresi logistik mencapai nilai accuracy 92.42%, sedangkan dengan k-nearest neighbors memperoleh nilai accuracy sebesar 87.88%.

Conclusion

Kesimpulan dari permodelan yang dilakukan, model dengan regresi logistik dengan metode stepAIC() memiliki performa terbaik dan nilai AIC terkecil. Namun KNN masih bisa menjadi pilihan yang tidak terlalu buruk karena performanya tidak jauh berbeda dengan regresi logistik.

Jika mengacu pada hasil evaluasi confusionMatrix, saya akan menyarankan untuk berfokus pada nilai specificity. Specificity memang menitikberatkan pada kelas negatif dimana pada kasus ini kelas negatif adalah Horeca yang terdiri dari hotel, restaurant, dan cafe. Distribusi produk yang disalurkan ke tempat sejenis Horeca rasanya akan “berputar” lebih cepat. Artinya produk yang didistribusikan cenderung akan lebih sering digunakan dan di-restock sehingga penjualan dan proses distribusi bisa berjalan dengan cepat juga. Maka ketika nilai True Negative di confusionMatrix tergolong besar justru keadaan itu cukup menguntungkan. Sebab sekali lagi, kelas positif dan negatif dalam klasifikasi kali ini tidak menandakan kebaikan atau keburukan suatu kelas.