Introduction: Machine Learning

Tujuan dari machine learning adalah membuat mesin yang dapat belajar sendiri dalam memahami pola data hingga mengestimasi apa yang akan terjadi di masa depan.

Supervised vs Unsupervised Learning

  1. Supervised Learning: memiliki target variabel

    • Regression: target variabelnya numerik

    • Classification: target variabelnya kategorik

  2. Unsupervised Learning: tidak memiliki target variabel

    • Clustering

    • Dimensionality Reduction

Linear Regression

Regresi linier adalah suatu metode dalam machine learning yang digunakan untuk memodelkan hubungan linier antara variabel input (fitur) dan variabel output (target). Tujuan utama regresi linier adalah untuk menemukan garis lurus terbaik (linear) yang dapat memprediksi nilai target berdasarkan nilai fitur-fitur yang diberikan.

y=mx+b

  • y adalah variabel target yang ingin diprediksi,
  • x adalah variabel input (fitur),
  • m adalah koefisien kemiringan (slope) garis regresi,
  • b adalah intercept (konstanta).
copiers <- read.csv("data_input/copiers.csv")
head(copiers,2)

Inspeksi Data

str(copiers)
#> 'data.frame':    62 obs. of  15 variables:
#>  $ Row.ID      : int  336 393 407 516 596 754 1151 1234 1550 1645 ...
#>  $ Order.ID    : chr  "CA-2015-137946" "US-2014-135972" "CA-2017-117457" "CA-2017-127432" ...
#>  $ Order.Date  : chr  "9/1/15" "9/21/14" "12/8/17" "1/22/17" ...
#>  $ Ship.Date   : chr  "9/4/15" "9/23/14" "12/12/17" "1/27/17" ...
#>  $ Ship.Mode   : chr  "Second Class" "Second Class" "Standard Class" "Standard Class" ...
#>  $ Customer.ID : chr  "DB-13615" "JG-15115" "KH-16510" "AD-10180" ...
#>  $ Segment     : chr  "Consumer" "Consumer" "Consumer" "Home Office" ...
#>  $ Product.ID  : chr  "TEC-CO-10001449" "TEC-CO-10002313" "TEC-CO-10004115" "TEC-CO-10003236" ...
#>  $ Category    : chr  "Technology" "Technology" "Technology" "Technology" ...
#>  $ Sub.Category: chr  "Copiers" "Copiers" "Copiers" "Copiers" ...
#>  $ Product.Name: chr  "Hewlett Packard LaserJet 3310 Copier" "Canon PC1080F Personal Copier" "Sharp AL-1530CS Digital Copier" "Canon Image Class D660 Copier" ...
#>  $ Sales       : num  960 1800 1200 3000 1200 ...
#>  $ Quantity    : int  2 3 3 5 3 3 2 2 1 7 ...
#>  $ Discount    : num  0.2 0 0.2 0 0.2 0.2 0 0.4 0.2 0 ...
#>  $ Profit      : num  336 702 435 1380 435 ...

BQ: Misalkan kita ingin memprediksi profit berdasarkan nilai sales. Tentukan:

  • variable target (y): Profit
  • variable prediktor (x): sales

Eksplorasi data:

1. Cek persebaran data

Cek persebaran variable Profit:

boxplot(copiers$Profit)

Insight: - outlier -> ada -> bisa memengaruhi - distribusi -> rendah, berkisar antara 200-800 an

Cek persebaran variable Sales:

boxplot(copiers$Sales)

## nilai korelasi

cor(copiers$Sales, copiers$Profit)
#> [1] 0.9395785

Korelasi antara Profit dan Sales kuat positif.

# visualisasi scatter plot
plot(copiers$Sales, copiers$Profit)

Indikasi bahwa sales adalah prediktor yang baik untuk memprediksi Profit.

Simple Linear Regression

Model tanpa variable prediktor:

lm(y ~ x, data)

# buat model
model_base <- lm(formula = Profit ~ 1, data = copiers) # 1 menandakan tanpa prediktor

# print summary model
summary(model_base)
#> 
#> Call:
#> lm(formula = Profit ~ 1, data = copiers)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -388.8 -286.8 -135.3  162.4 1854.2 
#> 
#> Coefficients:
#>             Estimate Std. Error t value        Pr(>|t|)    
#> (Intercept)   448.81      55.24   8.124 0.0000000000269 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 435 on 61 degrees of freedom

Model dengan variable prediktor:

# buat model
model_ols <- lm(formula = Profit ~ Sales, data = copiers)

# summary model
summary(model_ols)
#> 
#> Call:
#> lm(formula = Profit ~ Sales, data = copiers)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -407.07  -70.08   22.95   76.56  345.05 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept) -114.06251   32.62743  -3.496             0.000895 ***
#> Sales          0.42286    0.01989  21.260 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 150.1 on 60 degrees of freedom
#> Multiple R-squared:  0.8828, Adjusted R-squared:  0.8809 
#> F-statistic:   452 on 1 and 60 DF,  p-value: < 0.00000000000000022

Dihasilkan koefisien intercept dan koefisien untuk setiap prediktor, sehingga formula model regresi:

\[\hat{y}=\beta_0+\beta_1.x_1+...+\beta_n.x_n\].

dimana, \(\beta_0\) adalah intercept, \(\beta_1, ..., \beta_n\) adalah coefficient prediktor, dan \(x_1,...,x_n\) merupakan variable prediktor yang digunakan.

Sehingga formula model yang diperoleh:

\[Profit = -114.0625 + 0.4229*{Sales}\]

Bila garis regresi diplot kan maka menghasilkan sebagai berikut:

plot(copiers$Sales, copiers$Profit)
abline(model_ols, col = "red")

  • intercept: titik dimana garis bersinggungan dengan sumbu y
  • coefficient prediktor: slope/kemiringan garis

Important Points:

  • Pada linear regression, prediktor yang baik adalah yang memiliki korelasi kuat dengan target.
  • Model Linear Regression membuat suatu garis lurus yang menangkap pola data sehingga menghasilkan error terkecill.
  • Formula garis regresi dapat digunakan untuk prediksi.

Ordinary Least Squares (Concept)

Ordinary least squares bekerja dengan mencari suatu garis lurus yang dapat merepresentasikan pola data, dalam kata lain menghasilkan error terkecil. Error yang dimaksud adalah Sum of Squared Error (SSE)

  • Error: selisih dari nilai prediksi (\(\hat{y}\)) dengan nilai actual (\(y\))

\[error = y - \hat{y}\]

  • Sum of Squared Error (SSE):

\[\sum^n_{i=1}(y_i - \hat{y}_i)^2\]

Penjelasan formula:

  1. Error dihitung untuk masing-masing data
  2. Dikuadratkan agar nilai error positif dan error negatif tidak saling menghilangkan
  3. Dijumlahkan untuk mendapatkan suatu nilai yang merepresentasikan error keseluruhan data

Interpretasi Model

Model Machine Learning memiliki beberapa sifat:

  • robustness: performa model yang tinggi, errornya rendah
  • interpretability: bisa tau pengaruh tiap prediktor (meningkatkan/menurunkan target)

Umumnya:

  • model sederhana: interpretability tinggi, robustness rendah
  • model kompleks: interpretability rendah, robustness tinggi

Salah satu keunggulan model linear regression adalah interpretability.

summary(model_ols)
#> 
#> Call:
#> lm(formula = Profit ~ Sales, data = copiers)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -407.07  -70.08   22.95   76.56  345.05 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept) -114.06251   32.62743  -3.496             0.000895 ***
#> Sales          0.42286    0.01989  21.260 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 150.1 on 60 degrees of freedom
#> Multiple R-squared:  0.8828, Adjusted R-squared:  0.8809 
#> F-statistic:   452 on 1 and 60 DF,  p-value: < 0.00000000000000022

Profit = -114.06251 + 0.42286 * Sales

  1. Intercept: titik awal garis regresi terbentuk, menunjukkan nilai target ketika nilai prediktor = 0

    Saat Sales 0, profit = -114.06251 (rugi)

  2. Coefficient/Slope: kenaikan variable target setiap 1 satuan

  • Koefisien positif = korelasi positif, meningkatkan nilai variable target
  • Koefisien negatif = korelasi negatif, menurunkan nilai variable target

Sales menaikkan Profit

  1. Signifikansi prediktor: mengetahui apakah setiap prediktor berpengaruh signifikan terhadap variable targetnya.
  • Sebuah prediktor dikatakan signifikan ketika p-value < 0.05 (alpha)
  • Bisa juga dilihat dari jumlah bintang setiap prediktor

Variable yg signifikan: Sales

  1. R-squared: ukuran kebaikan model. Seberapa baik model dapat menjelaskan target.
  • rentang nilai 0-1, mendekati 1 semakin baik

Prediktor yg kita gunakan pada model bisa menjelaskan sebanyak 88.28% variansi dari target variable, sedangkan sisanya dijelaskan oleh variable lain di luar model.

Latihan 1

Buat model prediksi Profit dengan prediktor Discount

# buat model
model_ols_disc <- lm(Profit ~ Discount, copiers)

# summary
summary(model_ols_disc)
#> 
#> Call:
#> lm(formula = Profit ~ Discount, data = copiers)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -608.26 -248.27  -66.55   87.73 1604.71 
#> 
#> Coefficients:
#>             Estimate Std. Error t value        Pr(>|t|)    
#> (Intercept)   698.25      83.87   8.326 0.0000000000137 ***
#> Discount    -1487.10     400.41  -3.714         0.00045 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 395.5 on 60 degrees of freedom
#> Multiple R-squared:  0.1869, Adjusted R-squared:  0.1734 
#> F-statistic: 13.79 on 1 and 60 DF,  p-value: 0.0004496

Interpretasi:

  1. Intercept: > Saat tidak ada Discount, maka Profit senilai 689.25

  2. Coef Discount: > Korelasi negatif, discount menurunkan Profit. > Ketika Discount naik 1 satuan maka menurunkan Profit sebesar 1487.19

  3. Prediktor Significant: Discount

  4. R-squared: Model bisa menjelaskan sebanyak 18.69% dari target.

Prediksi nilai profit misalkan kita punya data sales sebagai berikut:

sales <- data.frame(Sales = c(300,290,320,450))
sales

Agar bisa menghitung langsung semuanya, gunakan function predict():

predict(object = model_ols, newdata = sales)
#>         1         2         3         4 
#> 12.795125  8.566537 21.252301 76.223945

Leverage vs Influence

Outlier adalah data yang nilainya jauh dari kebanyakan data lainnya, dan umumnya mengganggu pemodelan.

Outlier dimana nilai sumbu x nya jauh dari kebanyakan data disebut high leverage. Outlier berpotensi untuk sangat mempengaruhi koefisien (intercept, slope) dan r-squared dari model regresi (high influence). Untuk membuktikannya kita perlu melihat summary model.

  • leverage mempengaruhi model: high influence
  • leverage tidak mempengaruhi model: low influence

Practical Notes:

  1. Buat model dengan outlier (data awal) -> cek performa model
  2. Kalau performa masih buruk -> cek ada outlier/tidak
  3. Kalau ada outlier, coba dibuang -> buat model tanpa outlier
  4. Cek performa model lagi
    • bila outlier tidak merubah/meningkatkan R-squared, sebaiknya dipertahankan.
    • bila outlier menurunkan R-squared, sebaiknya dibuang.

Sejauh ini, kita sudah membuat model_ols yang menggunakan seluruh observasi. Lakukan exploratory data dengan melihat apakah terdapat outlier pada variable Sales:

# boxplot sales
boxplot(copiers$Sales)

# plot target dan prediktor
plot(copiers$Sales, copiers$Profit)
abline(model_ols, col = "red")

Coba buang data outlier:

# versi base
copiers_new <- copiers[copiers$Sales < 4000,]

# versi dplyr
library(dplyr)
copiers_new <- copiers %>% 
            filter(Sales < 4000)

Buatlah model linear regression dengan data tanpa outlier:

model_no_outlier <- lm(Profit ~ Sales, copiers_new)

Bandingkan garis regresi yang dibentuk dari model_ols dan model_no_outlier, apakah data outlier menghasilkan output yang jauh berbeda?

plot(copiers$Sales, copiers$Profit)
abline(model_ols, col = "red")
abline(model_no_outlier, col = "green")

Insight: Secara visual, garis prediksi model tanpa outlier (garis hijau) tidak jauh berbeda dengan garis regresi dengan outlier (garis merah). Indikasi bahwa outlier tersebut Low influence

Mari pastikan pula perbandingan nilai goodness of fit (R-squared) nya:

summary(model_ols)$r.squared
#> [1] 0.8828077
summary(model_no_outlier)$r.squared
#> [1] 0.8513778

Insight: r-squared lebih bagus pada model_ols(model dengan outlier), maka outlier tidak perlu dibuang.

Multiple Linear Regression

Penggunaan lebih dari 1 prediktor dapat meningkatkan performa model karena lebih banyak informasi yang dapat menjelaskan target.

Pemilihan prediktor (feature selection):

  • berdasarkan bisnis

  • berdasarkan statistik:

    • correlation
    • stepwise regression
  1. Identifikasi Proporsi ship mode
table(copiers$Ship.Mode)
#> 
#>    First Class       Same Day   Second Class Standard Class 
#>             10              5             15             32

Penentuan variable:

  • target: Profit
  • prediktor (berdasarkan business):
    • buang: row.id, order.id, product.id, customer.id, category & sub category
    • buang product.name (yang memiliki jumlah kategori terlalu banyak, atau kemungkinan di lapangan akan ada data baru) -> buang aja. karena model kita jadi terlalu spesifik, hanya bisa memprediksi yang ada pada data awal
    • buang tipe date (full date): jarang digunakan sebagai prediktor: pada data baru pasti tanggalnya berbeda.
  1. Data Wrangling:
# seleksi dan ubah tipe data
copiers <- copiers %>% 
  select(Ship.Mode, Segment, Sales, Quantity, Discount, Profit) %>% 
  mutate(Ship.Mode = as.factor(Ship.Mode),
         Segment = as.factor(Segment))

# cek missing value = tidak ada
anyNA(copiers)
#> [1] FALSE
# cek data
head(copiers)
  1. Eksplorasi data (Exploratory Data Analysis)

Cek korelasi prediktor - target:

# cara manual
cor(copiers$Profit, copiers$Sales)
#> [1] 0.9395785
cor(copiers$Profit, copiers$Quantity)
#> [1] 0.6859931
cor(copiers$Profit, copiers$Discount)
#> [1] -0.4323383
# cara 2: menggunakan fungsi `ggcorr()`
library(GGally)
ggcorr(copiers, label = T)

Insight: prediktor potensial berdasarkan kuatnya nilai korelasi: sales - quantity paling kuat. discount agak lemah, tapi tidak terlalu lemah jadi masih bisa dicoba utk jadi prediktor.

  1. Modelling

Bentuk model menggunakan seluruh variable prediktor:

# buat model
model_ols_multi <- lm(Profit ~ ., copiers)  

# summary model
summary(model_ols_multi)
#> 
#> Call:
#> lm(formula = Profit ~ ., data = copiers)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -285.64  -56.21    6.49   66.51  246.10 
#> 
#> Coefficients:
#>                           Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)               63.45200   55.06367   1.152                0.254    
#> Ship.ModeSame Day          9.13795   57.48532   0.159                0.874    
#> Ship.ModeSecond Class     42.97158   43.92567   0.978                0.332    
#> Ship.ModeStandard Class   13.92121   38.33509   0.363                0.718    
#> SegmentCorporate           9.30122   31.29427   0.297                0.767    
#> SegmentHome Office       -30.70992   39.86033  -0.770                0.444    
#> Sales                      0.42124    0.02522  16.701 < 0.0000000000000002 ***
#> Quantity                 -13.40975   13.54512  -0.990                0.327    
#> Discount                -874.83099  117.02350  -7.476       0.000000000774 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 102.7 on 53 degrees of freedom
#> Multiple R-squared:  0.9515, Adjusted R-squared:  0.9442 
#> F-statistic: 130.1 on 8 and 53 DF,  p-value: < 0.00000000000000022

Formula regresi untuk model di atas:

\[ Profit = 63.452 \\ + 9.13795 * Ship.ModeSame Day \\ + 42.97158 * Ship.ModeSecond Class \\ + 13.92121 * Ship.ModeStandard Class \\ + 9.30122 * SegmentCorporate \\ - 30.70992 * SegmentHome Office \\ + 0.42124 * Sales \\ - 13.40975 * Quantity \\ - 874.83099 * Discount \]

Insight:

  • interpretasi prediktor (numerik):
    • meningkatkan profit: Sales (karena coef bernilai positif)
    • menurunkan profit: Quantity, Discount (karena coef bernilai negatif)
  • interpretasi prediktor (kategorik):
    • setiap nilai dari kolom kategorik akan dijadikan 1 entitias (dummy variable)
    • akan ada 1 nilai yang menjadi basis (di intercept), tidak ditampilkan koefisiennya
    • cara interpretasi:
      • Segmen = Corporate, maka target meningkat
      • Segmen = Home Office maka target menurun
      • Segmen = Consumer, maka target tidak ditambah/dikurang (mengikuti intercept)
  • signifikansi prediktor: Sales & Discount.
    • bila ada variable yg tidak signifikan, kita bisa coba hapus saja dr model. karena mungkin pengaruhnya tidak banyak ke target kita. mungkin hanya memberikan noise/data random yang bisa memperburuk data kita.
    • untuk memastikannya harus dicoba menggunakan pembuatan model baru
  • adj. r-squared: 0.9442
summary(model_ols)$r.squared
#> [1] 0.8828077
summary(model_ols_multi)$adj.r.squared
#> [1] 0.9442254

Note:

  • kalau simple linear regression: multiple r squared

  • kalau multiple linear regression: adj. r. squared

    • semakin banyak prediktor maka R-squared otomatis meningkat
    • adjusted R-squared menghitung R-squared sambil memperhitungkan/melakukan penalty untuk jumlah prediktor yang digunakan. Adj. R-squared hanya akan meningkat bila prediktor memang menghasilkan prediksi yang lebih baik, bukan hanya menambahkan informasi random.

  • Interpretasi model: summary(model)

    • koefisien prediktor: positif/negatif -> meningkatkan/menurunkan target

    • prediktor signifikan: lihat dari p-value/jumlah bintang

    • r-squared:

      • nilai kebaikan model dalam menjelaskan target
      • rentang 0-1: mendekati 1 semakin baik

  • Leverage: outliers pada prediktor

    • Leverage dapat memengaruhi model secara signifikan
    • Leverage tetap dimasukkan ke dalam model apabila berpengaruh baik pada model

  • Multiple linear regression: prediktor > 1

    • r.squared:

      • utk multiple linear regression gunakan adj.r.squared
      • utk simple linear regression gunakan multiple r.squared

  • Prediksi: predict(model, newdata)

Karena terdapat variable tidak signifikan, mari coba gunakan variable yang signifikan saja:

model_ols_multi2 <- lm(Profit ~ Sales + Discount, copiers) 

# summary model
summary(model_ols_multi2)
#> 
#> Call:
#> lm(formula = Profit ~ Sales + Discount, data = copiers)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -277.69  -62.73    6.89   66.79  275.86 
#> 
#> Coefficients:
#>               Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)   66.21587   30.34849   2.182               0.0331 *  
#> Sales          0.40019    0.01363  29.370 < 0.0000000000000002 ***
#> Discount    -894.82460  104.13831  -8.593     0.00000000000549 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 100.9 on 59 degrees of freedom
#> Multiple R-squared:  0.9479, Adjusted R-squared:  0.9462 
#> F-statistic: 537.2 on 2 and 59 DF,  p-value: < 0.00000000000000022

Insight: adj. r-squared adalah 0.9462

Prediksi

Kita akan gunakan model yang telah dibuat:

  • model_ols: 1 prediktor (Sales)
  • model_ols_multi: Semua prediktor
  • model_ols_multi2: Prediktor signifikan (Sales + Discount)
# simpan hasil prediksi ke kolom baru `prediction` di data `copiers_multi`
pred_ols <- predict(model_ols, newdata = copiers)
pred_ols_multi <- predict(model_ols_multi, newdata = copiers)
pred_ols_multi2 <- predict(model_ols_multi2, newdata = copiers)

Evaluasi Model

Kita akan bandingkan dari ketiga model, manakah yang memiliki performa terbaik melalui evaluasi model! Diinginkan error prediksi yang terkecil.

error = data aktual - hasil prediksi

Tiap baris observasi akan menghasilkan error, sehingga butuh 1 nilai untuk mewakili nilai error tersebut, umumnya adalah rata-rata error. Berikut beberapa metrics error yang bisa digunakan untuk regresi:

  1. MAE (Mean Absolute Error)

  2. MSE (Mean Squared Error)

  3. RMSE (Root Mean Square Error)

  4. Pilih error yang dipakai: berdasarkan benefit dan limitasinya

  • mae: kalau harus menjelaskan ke orang yg bukan berlatar belakang statistik dan bisa mengabaikan error outlier yang besar
  • mse: tidak bisa diinterpretasi. biasanya cuma untuk membandingkan model terbaik untuk kita pribadi, lebih sulit untuk dijelaskan ke khalayak.
  • rmse: rmse lebih merepresentasikan error (kalau ada error besar bisa terdeteksi), dan bisa diinterpretasikan juga karena sudah nilai asli. biasanya presentasi untuk khalayak akademisi/analis data.

Error yang dipilih: RMSE -> bisa interpretasikan dan sensitif dengan oulier.

Berdasarkan nilai R-squared:

summary(model_ols)$r.squared
#> [1] 0.8828077
summary(model_ols_multi)$adj.r.squared
#> [1] 0.9442254
summary(model_ols_multi2)$adj.r.squared
#> [1] 0.9461829

Based on error:

library(MLmetrics)
# your code
RMSE(y_pred = pred_ols, y_true = copiers$Profit)
#> [1] 147.6997
RMSE(y_pred = pred_ols_multi, y_true = copiers$Profit)
#> [1] 94.97767
RMSE(y_pred = pred_ols_multi2, y_true = copiers$Profit)
#> [1] 98.43544

Kesimpulan: Dengan mengedepankan error, model yang dipilih adalah model_ols_multi

Bila kita ingin mengetahui apakah RMSE tersebut sudah cukup kecil untuk data kita, bandingkan dengan range variable target:

range(copiers$Profit)
#> [1]   59.998 2302.967

Interpretasi RMSE: hasil prediksi yang dihasilkan model memiliki kemungkinan salah +- 94.997 dari titik prediksinya.

Note: interpretasi MAE pun caranya sama (+- dari hasil prediksi)

Case Study: Inequality Prediction

  1. Read data crime:
crime <- read.csv("data_input/crime.csv") %>% select(-X)
names(crime) <- c("percent_m", "is_south", "mean_education", "police_exp60", "police_exp59", "labour_participation", "m_per1000f", "state_pop", "nonwhites_per1000", "unemploy_m24", "unemploy_m39", "gdp", "inequality", "prob_prison", "time_prison", "crime_rate")
head(crime,2 )

Deskripsi kolom:

  • percent_m: percentage of males aged 14-24
  • is_south: whether it is in a Southern state. 1 for Yes, 0 for No.
  • mean_education: mean years of schooling
  • police_exp60: police expenditure in 1960
  • police_exp59: police expenditure in 1959
  • labour_participation: labour force participation rate
  • m_per1000f: number of males per 1000 females
  • state_pop: state population
  • nonwhites_per1000: number of non-whites resident per 1000 people
  • unemploy_m24: unemployment rate of urban males aged 14-24
  • unemploy_m39: unemployment rate of urban males aged 35-39
  • gdp: gross domestic product per head
  • inequality: income inequality
  • prob_prison: probability of imprisonment
  • time_prison: average time served in prisons
  • crime_rate: crime rate in an unspecified category

Target = inequality

  1. Data Wrangling:
# ubah tipe data
crime <- crime %>% 
  mutate(is_south = as.factor(is_south))

# cek missing value
anyNA(crime)
#> [1] FALSE
  1. EDA - Cek korelasi variable:
# hint: bila tulisan tidak rapih, tambahkan argumen di fungsi:
# `hjust = 1` agar tulisan tidak menyentuh kotak
# `layout.exp = 3` agar panel membesar 3 poin

ggcorr(crime, label = T, hjust = 1, layout.exp = 5)

Insight: Variable dengan korelasi kuat dengan inequality adalah gdp (-0.9) dan mean_education (-0.8)

  1. Buat model regresi linear untuk memprediksi inequality berdasarkan 1 variable paling potensial:
# buat model
model_ineq <- lm(inequality ~ gdp, crime)
# summary model
summary(model_ineq)
#> 
#> Call:
#> lm(formula = inequality ~ gdp, data = crime)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -42.029 -11.754   1.714  12.438  30.006 
#> 
#> Coefficients:
#>              Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept) 386.03058   15.38651   25.09 <0.0000000000000002 ***
#> gdp          -0.36551    0.02881  -12.69 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 18.86 on 45 degrees of freedom
#> Multiple R-squared:  0.7815, Adjusted R-squared:  0.7766 
#> F-statistic: 160.9 on 1 and 45 DF,  p-value: < 0.00000000000000022

Interpretasi model:

  • variable prediktor = menurunkan nilai target
  • signifikansi prediktor: significant
  • r-squared: 0.7815
  1. Buat model linear menggunakan seluruh variable:
# buat model
model_ineq_all <- lm(inequality ~ . , crime)

# summary model
summary(model_ineq_all)
#> 
#> Call:
#> lm(formula = inequality ~ ., data = crime)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -22.2069  -8.3637   0.5959   9.5326  28.0975 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)          330.041730 120.913748   2.730  0.01036 *  
#> percent_m             -0.383670   0.299996  -1.279  0.21042    
#> is_south1             18.794908   9.696399   1.938  0.06173 .  
#> mean_education        -1.067204   0.448779  -2.378  0.02375 *  
#> police_exp60          -0.374454   0.767441  -0.488  0.62904    
#> police_exp59          -0.010789   0.822096  -0.013  0.98961    
#> labour_participation   0.114750   0.099665   1.151  0.25839    
#> m_per1000f             0.079852   0.141404   0.565  0.57634    
#> state_pop              0.132739   0.086233   1.539  0.13388    
#> nonwhites_per1000      0.002937   0.045031   0.065  0.94842    
#> unemploy_m24           0.100185   0.298897   0.335  0.73975    
#> unemploy_m39          -0.303733   0.602692  -0.504  0.61785    
#> gdp                   -0.239312   0.058431  -4.096  0.00028 ***
#> prob_prison           51.196445 167.739663   0.305  0.76224    
#> time_prison            0.064390   0.496290   0.130  0.89761    
#> crime_rate             0.033667   0.010822   3.111  0.00398 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 14.43 on 31 degrees of freedom
#> Multiple R-squared:  0.9118, Adjusted R-squared:  0.8692 
#> F-statistic: 21.38 on 15 and 31 DF,  p-value: 0.000000000002959

Interpretasi:

  • signifikansi prediktor: gdp, crime rate, mean education
  • r-squared: 0.8692
  1. Evaluasi Model

Bandingkan R-squared model_ineq dan model_ineq_all:

summary(model_ineq)$r.squared
#> [1] 0.7814512
summary(model_ineq_all)$adj.r.squared
#> [1] 0.8691871

Bandingkan RMSE model_ineq dan model_ineq_all:

# prediksi
pred_model_ineq <- predict(model_ineq, crime)
pred_model_ineq_all <- predict(model_ineq_all, crime)

# evaluasi
RMSE(y_pred = pred_model_ineq, y_true = crime$inequality)
#> [1] 18.45161
RMSE(y_pred = pred_model_ineq_all, y_true = crime$inequality)
#> [1] 11.71891

Kesimpulan: Berdasarkan R-squared dan error dipilih model model_ineq_all

Step-wise Regression

Pemilihan prediktor pada model, selain didasari business insight bisa didukung secara statistik.

Stepwise regression mencari kombinasi prediktor yang menghasilkan model terbaik berdasarkan nilai AIC (Akaike Information Criterion / Information Loss). AIC menunjukkan banyak informasi yang hilang pada model. Diinginkan AIC yang semakin kecil.

  • backward elimination
  • forward selection
  • both

Backward

Mulai dari menggunakan seluruh prediktor, kemudian dikurangi 1 per 1 variable prediktor sehingga diperoleh AIC terkecil.

  1. Menggunakan seluruh variable

  2. Memilih variable yang ketika dihilangkan, menghasilkan AIC terbaik (rendah)

  3. Menghilangkan variable yang dipilih dari tahap kedua

  4. Mengulangi step 2 dan step 3

  5. Proses berhenti bila saat dilakukan pengurangan variable lagi malah menghasilkan AIC yang lebih tinggi.

# persiapan: buat model seluruh prediktor
model_ineq_all
#> 
#> Call:
#> lm(formula = inequality ~ ., data = crime)
#> 
#> Coefficients:
#>          (Intercept)             percent_m             is_south1  
#>           330.041730             -0.383670             18.794908  
#>       mean_education          police_exp60          police_exp59  
#>            -1.067204             -0.374454             -0.010789  
#> labour_participation            m_per1000f             state_pop  
#>             0.114750              0.079852              0.132739  
#>    nonwhites_per1000          unemploy_m24          unemploy_m39  
#>             0.002937              0.100185             -0.303733  
#>                  gdp           prob_prison           time_prison  
#>            -0.239312             51.196445              0.064390  
#>           crime_rate  
#>             0.033667
# stepwise regression - backward
model_backward <- step(object = model_ineq_all, direction = "backward", trace = F)
# cek summary model
summary(model_backward)
#> 
#> Call:
#> lm(formula = inequality ~ percent_m + is_south + mean_education + 
#>     police_exp60 + labour_participation + state_pop + gdp + crime_rate, 
#>     data = crime)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -23.7645  -7.3648  -0.2304   9.5643  30.5817 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value     Pr(>|t|)    
#> (Intercept)          383.751874  55.466625   6.919 0.0000000317 ***
#> percent_m             -0.322926   0.236604  -1.365     0.180332    
#> is_south1             18.885387   6.700763   2.818     0.007621 ** 
#> mean_education        -0.931265   0.338854  -2.748     0.009114 ** 
#> police_exp60          -0.378395   0.155792  -2.429     0.019987 *  
#> labour_participation   0.137360   0.064315   2.136     0.039206 *  
#> state_pop              0.101107   0.066067   1.530     0.134207    
#> gdp                   -0.250171   0.048611  -5.146 0.0000083974 ***
#> crime_rate             0.032878   0.008102   4.058     0.000238 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 13.27 on 38 degrees of freedom
#> Multiple R-squared:  0.9086, Adjusted R-squared:  0.8894 
#> F-statistic: 47.23 on 8 and 38 DF,  p-value: < 0.00000000000000022

Forward

Mulai dari model tanpa prediktor, kemudian ditambahkan 1 per 1 variable prediktor sehingga diperoleh model AIC terkecil.

  1. Tidak menggunakan variable apapun

  2. Memilih variable yang ketika ditambahkan menghasilkan AIC terbaik (rendah)

  3. Menambahkan variable yang dipilih dari tahap kedua

  4. Mengulangi step 2 dan 3

  5. Proses berhenti bila saat dilakukan penambahan variable lagi malah menghasilkan AIC yang lebih tinggi.

# persiapan: buat model tanpa prediktor
model_ineq_none <- lm(inequality ~ 1, data = crime)

Untuk forward selection, kita perlu mendefinisikan parameter scope untuk menandakan batas atas maksimal kombinasi prediktor.

# stepwise regression - forward
model_forward <- step(object = model_ineq_none, scope = list(upper = model_ineq_all), direction = "forward", trace = F )
summary(model_forward)
#> 
#> Call:
#> lm(formula = inequality ~ gdp + crime_rate + mean_education + 
#>     police_exp59 + is_south + labour_participation + state_pop + 
#>     percent_m, data = crime)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -25.4435  -7.4800  -0.0233   9.0979  30.1273 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value     Pr(>|t|)    
#> (Intercept)          382.103157  55.503952   6.884 0.0000000353 ***
#> gdp                   -0.252009   0.048453  -5.201 0.0000070726 ***
#> crime_rate             0.031544   0.007828   4.029     0.000259 ***
#> mean_education        -0.874630   0.338545  -2.583     0.013752 *  
#> police_exp59          -0.386989   0.161516  -2.396     0.021605 *  
#> is_south1             19.208354   6.714136   2.861     0.006830 ** 
#> labour_participation   0.129018   0.065030   1.984     0.054516 .  
#> state_pop              0.098028   0.065852   1.489     0.144845    
#> percent_m             -0.312760   0.236543  -1.322     0.194003    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 13.29 on 38 degrees of freedom
#> Multiple R-squared:  0.9083, Adjusted R-squared:  0.889 
#> F-statistic: 47.05 on 8 and 38 DF,  p-value: < 0.00000000000000022

Both

Pengurangan dan penambahan prediktor dilakukan bersamaan hingga didapatkan AIC terkecil.

  1. Menggunakan model tanpa prediktor.

  2. Memilih penambahan variable atau pengurangan variable yang menghasilkan AIC terbaik.

  3. Menambahkan atau mengurangi variable yang dipilih dari tahap kedua

  4. Meneruskan step 2 dan 3 hingga nilai AIC terbaik dihasilkan.

# buat stepwise regression - both
model_both <- step(object = model_ineq_none, scope = list(upper = model_ineq_all), direction = "both", trace = F)
summary(model_both)
#> 
#> Call:
#> lm(formula = inequality ~ gdp + crime_rate + mean_education + 
#>     police_exp59 + is_south + labour_participation + state_pop + 
#>     percent_m, data = crime)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -25.4435  -7.4800  -0.0233   9.0979  30.1273 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value     Pr(>|t|)    
#> (Intercept)          382.103157  55.503952   6.884 0.0000000353 ***
#> gdp                   -0.252009   0.048453  -5.201 0.0000070726 ***
#> crime_rate             0.031544   0.007828   4.029     0.000259 ***
#> mean_education        -0.874630   0.338545  -2.583     0.013752 *  
#> police_exp59          -0.386989   0.161516  -2.396     0.021605 *  
#> is_south1             19.208354   6.714136   2.861     0.006830 ** 
#> labour_participation   0.129018   0.065030   1.984     0.054516 .  
#> state_pop              0.098028   0.065852   1.489     0.144845    
#> percent_m             -0.312760   0.236543  -1.322     0.194003    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 13.29 on 38 degrees of freedom
#> Multiple R-squared:  0.9083, Adjusted R-squared:  0.889 
#> F-statistic: 47.05 on 8 and 38 DF,  p-value: < 0.00000000000000022

Nilai adjusted R-squared pada ketiga model dapat dibandingkan:

summary(model_backward)$adj.r.squared
#> [1] 0.8893899
summary(model_forward)$adj.r.squared
#> [1] 0.888989
summary(model_both)$adj.r.squared
#> [1] 0.888989
summary(model_ineq_all)$adj.r.squared
#> [1] 0.8691871

Hasil ketiga step-wise regression kurang lebih sama. Bisa digunakan yang paling praktis untuk kita.

Pertanyaan: Kalau ada prediktor yang secara bisnis penting, namun dibuang pada stepwise regression, apa yang harus dilakukan?

Ditinjau dulu keyakinan terhadap konsiderasi bisnis - bila sangat yakin bahwa prediktor ini penting secara bisnis, sumbernya dari yg paham data, maka prediktor tersebut bisa dimasukkan kembali ke model - bila tidak begitu yakin, bisa dicek dahulu seberapa besar pengaruhnya dengan output yang ingin diprediksi. Buat model dengan & tanpa prediktor tersebut, cek dari nilai error atau R-squared nya.

NOTE: Stepwise regression merupakan greedy algorithm yang fokus pada mencari hasil baik dengan waktu paling cepat, namun belum tentu memberikan hasil paling optimal.

Prediction Interval

Setelah model terbuat, kita bisa gunakan untuk prediksi ke data baru.

Model dari hasil stepwise: model_backward

# prediction
pred_model_step <- predict(model_backward, newdata = crime)

Prediksi ke data baru menggunakan rentang jauh lebih aman dibandingkan hanya memprediksi 1 titik. Kita bisa cari prediksi dengan confidence interval:

# untuk menambahkan batas atas-bawah
pred_model_stepci <- predict(model_backward, newdata = crime, 
                             interval = "prediction", # menambahkan ci utk prediksi 
                             level = 0.95) # default confidence level

Insight: Observasi pertama diprediksi memiliki inequality di rentang 220.03144 - 275.9937.

Setelah lakukan prediksi, bisa lanjut ke evaluasi model (cek error):

Evaluasi Model:

# evaluation
summary(model_ineq)$r.squared # prediktor gdp
#> [1] 0.7814512
summary(model_ineq_all)$adj.r.squared # smua prediktor
#> [1] 0.8691871
summary(model_backward)$adj.r.squared # bberapa prediktor 
#> [1] 0.8893899
# evaluation
RMSE(y_pred = pred_model_ineq, y_true = crime$inequality) # prediktor gdp
#> [1] 18.45161
RMSE(y_pred = pred_model_ineq_all, y_true = crime$inequality) # smua prediktor
#> [1] 11.71891
RMSE(y_pred = pred_model_step, y_true = crime$inequality) # bberapa prediktor 
#> [1] 11.93082

Discussion: model mana yang terbaik diantara model_ineq, model_ineq_all dan model hasil stepwise regression?

Error terkecil itu dimiliki oleh model_ineq_all, namun R-squared terbaik adalah model_backward

Note: kalau kita lihat error antara model_ineq_all dengan model_backward tidak jauh berbeda. Bila seperti itu bisa dipertimbangkan dari segi banyaknya prediktor (efisiensi model) yang kita pakai:

# model all
summary(model_ineq_all)$call # 15 prediktor untuk error 11.71 poin
#> lm(formula = inequality ~ ., data = crime)
# model backward 
summary(model_backward)$call # 8 prediktor untuk error 11.93 poin
#> lm(formula = inequality ~ percent_m + is_south + mean_education + 
#>     police_exp60 + labour_participation + state_pop + gdp + crime_rate, 
#>     data = crime)

Kalau ketersediaan data/kemudahan dalam mengkoleksi data prediktornya mudah, bisa digunakan seluruh prediktor. Namun, kalau ada pertimbangan berupa kesulitan mengambil data, bisa gunakan model hasil stepwise.

Assumption

Untuk memastikan kebaikan model Linear Regression, ada beberapa cek asumsi yang perlu dilakukan:

  1. Linearity
  2. Normality of Residuals
  3. Homoscedasticity of Residuals
  4. No Multicollinearity

Note: Hanya perlu cek asumsi kepada model yang akan digunakan, tidak perlu semua kombinasi model dicek asumsinya.

Mari cek model terpilih untuk ketiga asumsi tersebut!

Linearity

Linearity artinya target variabel dengan prediktornya memiliki hubungan yang linear atau hubungannya bersifat garis lurus. Asumsi ini biasanya dilakukan sebelum membuat model pada tahap EDA.

  1. Berdasarkan nilai korelasi
ggcorr(crime, label = T, hjust = 1, layout.exp = 3)

  1. Uji statistik dengan cor.test()

Linearity hypothesis test:

# misal kita cek untuk 1 prediktor saja
# yaitu inequality dengan gdp
cor.test(crime$inequality, crime$gdp)
#> 
#>  Pearson's product-moment correlation
#> 
#> data:  crime$inequality and crime$gdp
#> t = -12.685, df = 45, p-value < 0.00000000000000022
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  -0.9340509 -0.7998850
#> sample estimates:
#>        cor 
#> -0.8839973

p-value >= alpha, maka gagal tolak H0 p-value < alpha, maka tolak H0

  • H0: korelasi tidak signifikan
  • H1: korelasi signifikan

Diinginkan H1

Kesimpulan: variabel berkorelasi signifikan, uji asumsi kita terpenuhi

Normality of Residuals

Harapannya ketika membuat model linear regression, error/residual yang dihasilkan berdistribusi normal. Artinya error banyak berkumpul disekitar angka 0. Untuk menguji asumsi ini dapat dilakukan:

  1. Visualisasi histogram residual menggunakan fungsi hist()
hist(model_backward$residuals) # error

  1. Uji statistik dengan shapiro.test()
shapiro.test(model_backward$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_backward$residuals
#> W = 0.98628, p-value = 0.8493

Shapiro-Wilk hypothesis test:

Note: dibandingkan p-value dengan alpha (0.05)

p-value >= 0.05 = terima H0 (gagal tolak H0) p-value < 0.05 = terima H1 (tolak H0)

  • H0: error/residual berdistribusi normal
  • H1: error/residual tidak berdistribusi normal

Diinginkan H0

Kesimpulan: error/residual berdistribusi normal, uji asumsi kita terpenuhi

Homoscedasticity

Visualisasi persebaran error, terkait homoscedasticity vs heteroscedasticity:

Highlight:

  • diinginkan error yg tidak berpola/menyebar random/konstan (homoscedasticity)
  • tak diinginkan error yg berpola (heteroscedasticity) -> masih ada pola yang belum tertangkap model
  1. Visualisasi scatterplot: model’s fitted value - model’s error:
# fitted values = hasil prediksi ke data untuk training model
plot(x = model_backward$fitted.values, 
     y = model_backward$residuals)
abline(h = 0, col = "red")

  1. Uji statistik dengan bptest() dari package lmtest
library(lmtest)
bptest(model_backward)
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_backward
#> BP = 8.7653, df = 8, p-value = 0.3625

Breusch-Pagan hypothesis test:

Note: dibandingkan p-value dengan alpha (0.05)

p-value >= 0.05 = terima H0 (gagal tolak H0) p-value < 0.05 = terima H1 (tolak H0)

  • H0: Variansi error menyebar konstan (Homoscedasticity)
  • H1: Variansi error menyebar tidak konstan/membentuk pola (Heteroscedasticity)

Diinginkan H0

Kesimpulan: variansi error menyebar konstan (Homoscedasticity), uji asumsiterpenuhi

Multicollinearity

Multicollinearity adalah kondisi dimana ada korelasi antar prediktor yang amat kuat (hingga 1/-1). Hal ini menandakan ada variable redundan pada model, yang seharusnya dapat dipilih salah satu saja dari variable yang hubungannya amat kuat tersebut.

Harapannya tidak terjadi multicollinearity.

Uji VIF (Variance Inflation Factor) dengan fungsi vif() dari package car:

library(car)
vif(model_backward)
#>            percent_m             is_south       mean_education 
#>             2.310224             2.691396             3.754538 
#>         police_exp60 labour_participation            state_pop 
#>             5.600947             1.765022             1.652978 
#>                  gdp           crime_rate 
#>             5.748421             2.565567

Penentuan:

  • nilai VIF > 10 : ada multicollinearity
  • nilai VIF < 10 : tidak ada multicollinearity

Diinginkan VIF < 10.

Kesimpulan: tidak ada multicollinearity, uji asumsi terpenuhi.

Pertanyaan: bila ada uji asumsi yang tidak terpenuhi apakah model kita masih bisa dianggap baik/bisa digunakan?

Bisa, tapi artinya ada resiko ketika menggunakan model tersebut:

  • Normality tidak terpenuhi: ada kemungkinan error yg besar (jauh dari 0). Besarnya bisa dilihat dari plot distribusi error.
  • Homoscedasticity tidak terpenuhi: ada kemungkinan pola error dari hasil prediksi
  • Multicollinearity tidak terpenuhi: model tidak efisien

Notes:

  • Step-wise Regression: mencoba kombinasi prediktor hingga dihasilkan AIC (information loss) terkecil.
    • tipe stepwise: backward, forward, both
    • hasil ketiganya tidak jauh berbeda, gunakan yang paling praktis
  • Prediction Interval
    • dapat ditambahkan saat predict()
    • agar model tidak hanya prediksi 1 nilai tapi menyediakan batas atas-bawah
  • Asumsi Linear Model
    1. linearity
    • uji signifikansi korelasi antara variabel prediktor dan target
    • uji dengan cor.test(): diinginkan H1
    1. normality of residual:
    • ingin residual berdistribusi normal
    • uji Shapiro-Wilk shapiro.test(): diinginkan H0
    1. homoscedasticity
    • ingin residual menyebar konstan/random/tak berpola
    • uji Breusch-Pagan bptest(): diinginkan H0 4 multicollienarity
    • ingin antar prediktor tidak ada yg berkorelasi amat kuat (-1/1)
    • uji vif(): diinginkan nilai <10