Capstone: Regression Model
Husain A. Hidayat
Januay 19, 2021r
Introduction
Penerimaan mahasiswa baru Universitas sangat kompetitif, terutama untuk kampus unggulan. Terkadang rasio antara pendaftar dan kuota yang diterima terlampau besar. Melalui model regresi linear ini, kita akan menganalisa faktor apa saja yang berpengaruh terhadap penerimaan mahasiswa. Selain itu juga model ini akan memprediksi peluang diterimanya berdasarkan profil mereka. Dataset yang digunakan pada model ini adalah data penerimaan pascasarjana di India yang dapat diunduh di site Kaggle.
Data Preparation
About the Data
Berikut penjelasan singkat dari masing-masing variabel di dataset ini:
- GRE Scores (290 to 340)
- TOEFL Scores (92 to 120)
- University Rating (1 to 5)
- Statement of Purpose (1 to 5)
- Letter of Recommendation Strength (1 to 5)
- Undergraduate CGPA (6.8 to 9.92)
- Research Experience (0 or 1)
- Chance of Admit (0.34 to 0.97)
Karena kita akan memprediksi peluang diterima (Chance of Admit) berdasarkan beberapa faktor, maka variabel Chance of Admit sebagai Target Variabel, sedangkan variabel lainnya sebagai variabel prediktor
Import Library and Data
library(dplyr)
library(tidyr)
library(readr)
library(GGally)
library(car)
library(scales)
library(lmtest)
library(ggplot2)
library(plotly)
library(MLmetrics)Load Data
grad_adm <- read_csv("data_input/Admission_Predict_Ver1.1.csv")##
## -- Column specification --------------------------------------------------------
## cols(
## `Serial No.` = col_double(),
## `GRE Score` = col_double(),
## `TOEFL Score` = col_double(),
## `University Rating` = col_double(),
## SOP = col_double(),
## LOR = col_double(),
## CGPA = col_double(),
## Research = col_double(),
## `Chance of Admit` = col_double()
## )Saya rename nama kolom yang mengandung spasi (whitespace) dan hilangkan kolom “Serial No.” dari dataset
#hilangkan kolom "Serial No. dari dataset"
grad_adm_rename <- grad_adm %>%
select(-"Serial No.")
names(grad_adm_rename) <- c('GRE_Score','TOEFL_Score',
'UnivRank_Score','SOP','LOR',
'CGPA','Research','Chance_of_Admit')
colnames(grad_adm_rename)## [1] "GRE_Score" "TOEFL_Score" "UnivRank_Score" "SOP"
## [5] "LOR" "CGPA" "Research" "Chance_of_Admit"10 Data urutan atas
head(grad_adm_rename,10)10 Data urutan bawah
tail(grad_adm_rename,10)EDA
Cek struktur dataset
glimpse(grad_adm_rename)## Rows: 500
## Columns: 8
## $ GRE_Score <dbl> 337, 324, 316, 322, 314, 330, 321, 308, 302, 323, 3...
## $ TOEFL_Score <dbl> 118, 107, 104, 110, 103, 115, 109, 101, 102, 108, 1...
## $ UnivRank_Score <dbl> 4, 4, 3, 3, 2, 5, 3, 2, 1, 3, 3, 4, 4, 3, 3, 3, 3, ...
## $ SOP <dbl> 4.5, 4.0, 3.0, 3.5, 2.0, 4.5, 3.0, 3.0, 2.0, 3.5, 3...
## $ LOR <dbl> 4.5, 4.5, 3.5, 2.5, 3.0, 3.0, 4.0, 4.0, 1.5, 3.0, 4...
## $ CGPA <dbl> 9.65, 8.87, 8.00, 8.67, 8.21, 9.34, 8.20, 7.90, 8.0...
## $ Research <dbl> 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, ...
## $ Chance_of_Admit <dbl> 0.92, 0.76, 0.72, 0.80, 0.65, 0.90, 0.75, 0.68, 0.5...dim(grad_adm_rename)## [1] 500 8Berdasarkan rangkuman data diatas, semua data berjenis numerikal. Terdiri dari 500 baris dan 8 kolom.
cek Missing Value
colSums(is.na(grad_adm_rename))## GRE_Score TOEFL_Score UnivRank_Score SOP LOR
## 0 0 0 0 0
## CGPA Research Chance_of_Admit
## 0 0 0Tidak ada missing value
Cek distribusi data
summary(grad_adm_rename)## GRE_Score TOEFL_Score UnivRank_Score SOP
## Min. :290.0 Min. : 92.0 Min. :1.000 Min. :1.000
## 1st Qu.:308.0 1st Qu.:103.0 1st Qu.:2.000 1st Qu.:2.500
## Median :317.0 Median :107.0 Median :3.000 Median :3.500
## Mean :316.5 Mean :107.2 Mean :3.114 Mean :3.374
## 3rd Qu.:325.0 3rd Qu.:112.0 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :340.0 Max. :120.0 Max. :5.000 Max. :5.000
## LOR CGPA Research Chance_of_Admit
## Min. :1.000 Min. :6.800 Min. :0.00 Min. :0.3400
## 1st Qu.:3.000 1st Qu.:8.127 1st Qu.:0.00 1st Qu.:0.6300
## Median :3.500 Median :8.560 Median :1.00 Median :0.7200
## Mean :3.484 Mean :8.576 Mean :0.56 Mean :0.7217
## 3rd Qu.:4.000 3rd Qu.:9.040 3rd Qu.:1.00 3rd Qu.:0.8200
## Max. :5.000 Max. :9.920 Max. :1.00 Max. :0.9700karena Research nilainya hanya 0 dan 1 maka saya ubah menjadi factor
#grad_adm_rename$Research <- as.factor(grad_adm_rename$Research)
grad_adm_rename$Research <- as.factor(ifelse(grad_adm_rename$Research == 0, "No", "Yes"))Cek distribusi data
boxplot(grad_adm_rename)
Korelasi antar variabel
pairs(grad_adm_rename)
Cek Korelasi antar variabel
ggcorr(grad_adm_rename, label = T)## Warning in ggcorr(grad_adm_rename, label = T): data in column(s) 'Research' are
## not numeric and were ignored
Dari representasi Correlation Map diatas, antar variabel numerik menghasilkan korelasi yang cukup kuat (diatas 0.5)
Modelling
Train-Test Split (Holdout)
Dari 500 data observasi saya bagi 80% (400 data) sebagai data-train dan 20% (100 data) sebagai data-test. Data-train digunakan untuk melakukan pemodelan dan data-test dianggap sebagai unseen data yang digunakan untuk menguji seberapa baik model yang dibuat.
set.seed(123)
row_data <- nrow(grad_adm_rename)
index <- sample(row_data, row_data*0.8)
data_train <- grad_adm_rename[ index, ]
data_test <- grad_adm_rename[ -index, ]Linear Regression Model Full-Predictor
membuat model full prediktor
model1 <- lm(Chance_of_Admit ~ .,data_train)
summary(model1)##
## Call:
## lm(formula = Chance_of_Admit ~ ., data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.270696 -0.022021 0.009249 0.033021 0.154912
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3036332 0.1181235 -11.036 < 2e-16 ***
## GRE_Score 0.0020428 0.0005575 3.664 0.000282 ***
## TOEFL_Score 0.0025839 0.0009604 2.690 0.007442 **
## UnivRank_Score 0.0053167 0.0042205 1.260 0.208512
## SOP 0.0023630 0.0051321 0.460 0.645459
## LOR 0.0152702 0.0047203 3.235 0.001320 **
## CGPA 0.1180706 0.0103240 11.437 < 2e-16 ***
## ResearchYes 0.0236419 0.0073194 3.230 0.001342 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05917 on 392 degrees of freedom
## Multiple R-squared: 0.8186, Adjusted R-squared: 0.8153
## F-statistic: 252.6 on 7 and 392 DF, p-value: < 2.2e-16Interpretasi: - Dengan menggunakan semua variabel, model ini memberikan pengaruh terhadap Chance_of_Admit sebesar 81.53%. sisanya dijelaskan oleh faktor lain. - Terdapat dua variabel (yakni UnivRank_Score dan SOP) yang secara statistik kurang berpengaruh signifikan terhadap variabel target
\[ y = -1.3036332 + 0.0020428 GRE\_Score + 0.0025839 TOEFL\_Score + 0.0053167 UnivRank\_Score + 0.0023630 SOP + 0.0152702 LOR + 0.1180706 CGPA + 0.0236419 Research`Yes` \]
Linear Regression Model Step-Wise
Selanjutnya akan dicoba menentukan variabel prediktor secara otomatis menggunakan step-wise regression dengan metode ‘backward’ elimination.
Metode step-wise regression ini akan menghasilkan formula optimum berdasarkan nilai AIC yang terendah, dimana semakin rendah nilai Akaike Information Criterion (AIC) tersebut, maka model yang dihasilkan semakin baik.
model2 <- step(model1, direction = "backward")## Start: AIC=-2253.96
## Chance_of_Admit ~ GRE_Score + TOEFL_Score + UnivRank_Score +
## SOP + LOR + CGPA + Research
##
## Df Sum of Sq RSS AIC
## - SOP 1 0.00074 1.3731 -2255.8
## - UnivRank_Score 1 0.00556 1.3780 -2254.3
## <none> 1.3724 -2254.0
## - TOEFL_Score 1 0.02534 1.3977 -2248.6
## - Research 1 0.03653 1.4089 -2245.5
## - LOR 1 0.03664 1.4090 -2245.4
## - GRE_Score 1 0.04701 1.4194 -2242.5
## - CGPA 1 0.45791 1.8303 -2140.8
##
## Step: AIC=-2255.75
## Chance_of_Admit ~ GRE_Score + TOEFL_Score + UnivRank_Score +
## LOR + CGPA + Research
##
## Df Sum of Sq RSS AIC
## <none> 1.3731 -2255.8
## - UnivRank_Score 1 0.00859 1.3817 -2255.2
## - TOEFL_Score 1 0.02655 1.3997 -2250.1
## - Research 1 0.03719 1.4103 -2247.1
## - LOR 1 0.04466 1.4178 -2244.9
## - GRE_Score 1 0.04635 1.4195 -2244.5
## - CGPA 1 0.48387 1.8570 -2137.0summary(model2)##
## Call:
## lm(formula = Chance_of_Admit ~ GRE_Score + TOEFL_Score + UnivRank_Score +
## LOR + CGPA + Research, data = data_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.269804 -0.021551 0.009993 0.032356 0.155307
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3066313 0.1178256 -11.090 < 2e-16 ***
## GRE_Score 0.0020210 0.0005549 3.642 0.000307 ***
## TOEFL_Score 0.0026303 0.0009542 2.757 0.006113 **
## UnivRank_Score 0.0060792 0.0038781 1.568 0.117786
## LOR 0.0159677 0.0044661 3.575 0.000393 ***
## CGPA 0.1190048 0.0101126 11.768 < 2e-16 ***
## ResearchYes 0.0238215 0.0073017 3.262 0.001201 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05911 on 393 degrees of freedom
## Multiple R-squared: 0.8185, Adjusted R-squared: 0.8157
## F-statistic: 295.3 on 6 and 393 DF, p-value: < 2.2e-16Dari hasil step-wise backward, nilai AIC optimalnya -2258.25, dimana dari tujuh variabel prediktor di data train, terdapat satu variabel yang tidak dipakai yakni “Statement of Purpose” (SOP)
\[ y = -1.3066313 + 0.002021 GRE\_Score + 0.0026303 TOEFL\_Score + 0.0060792 UnivRank\_Score + 0.0159677 LOR + 0.1190048 CGPA + 0.0238215 Research`Yes` \]
Evaluation
Performance
Performa kedua model tadi kita coba ukur menggunakan MAE, MAPE, RMSE
residual_train_model1 <- predict(model1, data_train)
residual_train_model2 <- predict(model2, data_train)
residual_test_model1 <- predict(model1, data_test)
residual_test_model2 <- predict(model2, data_test)Cek perbandingan evaluasi model1 pada data-train dengan metrik MAE, MAPE,dan RMSE
MAE_Trainmod1 <- MAE(residual_train_model1, data_train$Chance_of_Admit)
MAPE_Trainmod1 <- MAPE(residual_train_model1, data_train$Chance_of_Admit)
RMSE_Trainmod1 <- RMSE(residual_train_model1, data_train$Chance_of_Admit)
#evaluasi `model2` pada data-train dengan metrik MAE, MAPE,dan RMSE
MAE_Trainmod2 <- MAE(residual_train_model2, data_train$Chance_of_Admit)
MAPE_Trainmod2 <- MAPE(residual_train_model2, data_train$Chance_of_Admit)
RMSE_Trainmod2 <- RMSE(residual_train_model2, data_train$Chance_of_Admit)
#evaluasi `model1` pada data-test dengan metrik MAE, MAPE,dan RMSE
MAE_Testmod1 <- MAE(residual_test_model1, data_test$Chance_of_Admit)
MAPE_Testmod1 <- MAPE(residual_test_model1, data_test$Chance_of_Admit)
RMSE_Testmod1 <- RMSE(residual_test_model1, data_test$Chance_of_Admit)
#evaluasi `model2` pada data-test dengan metrik MAE, MAPE,dan RMSE
MAE_Testmod2 <- MAE(residual_test_model2, data_test$Chance_of_Admit)
MAPE_Testmod2 <- MAPE(residual_test_model2, data_test$Chance_of_Admit)
RMSE_Testmod2 <- RMSE(residual_test_model2, data_test$Chance_of_Admit)
df_metric <- data.frame(model = c("model1","model1", "model2", "model2"),
set = c("data_train", "data_test",
"data_train", "data_test"),
MAE_score = c(MAE_Trainmod1, MAE_Testmod1,
MAE_Trainmod2, MAE_Testmod2),
MAPE_score = c(MAPE_Trainmod1, MAPE_Testmod1,
MAPE_Trainmod2, MAPE_Testmod2),
RMSE_score = c(RMSE_Trainmod1, RMSE_Testmod1,
RMSE_Trainmod2, RMSE_Testmod2))
df_metricdf_metric %>% pivot_longer(
cols = c("MAE_score","MAPE_score","RMSE_score"),
names_to = "method",
values_to = "score_in_percentage"
) %>% mutate(score_in_percentage = round(score_in_percentage,5)*100) %>%
ggplot(aes(x= model,
y= score_in_percentage)) +
geom_col(aes(x= model,
y= score_in_percentage,
fill = method), position = "dodge") +
geom_label(aes(label = score_in_percentage, group = method),
position = position_dodge(width = 1), size = 3) +
facet_wrap(~set, ncol = 2, labeller = as_labeller(c("data_train"="Data Train",
"data_test"="Data Test"))) +
labs(
title = "Model Performance Comparation",
subtitle = "using MAE, MAPE, RMSE",
x = "dataset", y = "score (in percent)"
)
#ggplotly(plot_metrics)Dari hasil pengujian performance kedua model, skor uji antara data train dengan data test tidak berbeda jauh sehingga bisa disimpulkan bahwa model tidak mengalami overfit maupun underfit. JIka dibandingkan skor uji di datatest antara kedua model tersebut, performance model2lebih baik dibandingkan model1
Model Assumption
Sebagai model parametrik, least square harus memenuhi empat asumsi: linearity, normality, Homocesdasticity, dan Multicolinearity.
Sebagai model parametrik, terdapat beberapa asumsi yang perlu dipenuhi agar interpretasi model yang diperoleh tidak bersifat bias. Beberapa asumsi yang harus dipenushi model antara lain:
- Linearity : terdapat hubungan linear antara variabel prediktor dengan target variabel
- Normality : error atau residual berdistribusi normal
- Constant Variance of Error (Homocesdasticity) : residual atau error bersifat konstan atau tidak membentuk pola tertentu.
- No-Multikolinearity : tidak ada multikolinearitas antar variabel independen
Kita keempat uji asumsi terhadap model2
1. Linearity
df_res_vs_act <- data.frame(residual = model2$residuals,
fitted = model2$fitted.values)
df_res_vs_act %>% ggplot(aes(fitted, residual)) +
geom_point() + geom_smooth() +
geom_hline(aes(yintercept = 0)) +
theme(panel.grid = element_blank(), panel.background = element_blank())## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
dari residual plot, tidak terlihat pola yang signifikan sehingga bisa disimpulkan bahwa model ini linear.
2. Normality
Residual atau error dari data perlu berdistribusi normal. Ketika kita membuat sebuah model regresi linear, persamaan yang akan kita dapat adalah sebagai berikut:
\[ y = \alpha + \beta\times x + \epsilon \\ \epsilon \sim Normal(0, \sigma) \]
dimana (\(\epsilon\)) adalah error atau residual yang tidak ditangkap model.
Model berasumsi bahwa error ini memiliki distribusi normal yang berpusat di angka 0. Untuk memeriksa distribusi dari nilai residual apakah berdistribusi normal (membentuk kurva lonceng) atau tidak, kita bisa menggunakan histogram.
hist(model2$residuals)
plot(density(model2$residuals))
Terlihat dari plot histogram, distribusi residual banyak disekitar angka nol.
Alternatif lainnya adalah melakukan uji statistik menggunakan uji Shapiro-Wilk
\(H_0\) = Data berdistribusi normal \(H_1\) = Data tidak berdistribusi normal
shapiro.test(model2$residuals)##
## Shapiro-Wilk normality test
##
## data: model2$residuals
## W = 0.92047, p-value = 1.041e-13Hasil uji Shapiro-Wilk, p-value < 0.05 sehingga kesimpulan secara statistik residual model tidak berdistribusi normal (tolak \(H_0\))
Menggunakan QQ-plot, jika data berdistribusi normal maka semua titik ada di dalam area biru dan mengikuti pola garis biru
qqPlot(model2$residuals)
## [1] 194 243Dari hasil QQ-Plot, terdapat beberapa titik diluat pola garis biru
3. Homocesdasticity
Homocesdasticity menunjukkan bahwa residual atau error bersifat konstan atau tidak membentuk pola tertentu. Jika error membentuk pola tertentu seperti garis linear atau mengerucut, maka kita sebut dengan Heterocesdasticity dan akan berpengaruh pada nilai standard error pada estimate/koefisien prediktor yan bias (terlalu sempit atau terlalu lebar).
Untuk pengujian secara statistik kali ini menggunakan uji Bresuch-Pagan \[ H_0 = error\ bersifat\ konstan\ (Homocesdasticity)\\ H_1 = error\ bersifat\ tidak\ konstan\ (Heteroscesdasticity) \]
library(lmtest)
bptest(model2)##
## studentized Breusch-Pagan test
##
## data: model2
## BP = 17.033, df = 6, p-value = 0.009163Hasil uji Bresuch-Pagan, p-value < 0.05 sehingga kesimpulan secara statistik tolak \(H_0\) atau error tidak bersifat constant (heterocesdasticity)
Homocesdasticity bisa dicek secara visual dengan melihat apakah ada pola antara hasil prediksi (fitted values) dengan nilai residualnya.
data.frame(prediction = model2$fitted.values,
residual = model2$residuals) %>%
ggplot(aes(prediction, residual)) +
geom_hline(yintercept = 0) +
geom_point( color = "skyblue4") +
theme_minimal()
Pada plot berikut terlihat bahwa terdapat pola tertentu sehingga kita bisa menyimpulkan bahwa model sudah memiliki error yang konstan (heterocesdasticity)
4. Multicolinearity
Multikolinearitas terjadi ketika antar variabel prediktor yang digunakan pada model memiliki hubungan yang kuat. Ada atau tidak multikolinearitas dapat dilihat dari nilai VIF(Variance Inflation Factor). VIF merupakan ukuran yang menjelaskan seberapa besar variansi koefisien yang meningkat karena multikolinieritas
Ketika nilai VIF lebih dari 10 artinya terjadi multikolinearitas. Jika hal tersebut terjadi, bisa dipilih salah satu variabel yang dibuang dari model yang memiliki VIF > 10.
Asumsi:
- Antara target variabel dan prediktor ada korelasi
- Antara prediktor dengan prediktor lain tidak ada korelasi
library(car)
vif(model2)## GRE_Score TOEFL_Score UnivRank_Score LOR CGPA
## 4.261929 3.703423 2.305448 1.861680 4.103733
## Research
## 1.498096Nilai VIF dari masing-masing variabel prediktor kurang dari 10, sehingga tidak ditemukan Multicollinearity antar variabel (antar variabel prediktor saling independen).
Conclusion
Berdasarkan modelling dan hasil analisis regresi linear yang telah dilakukan dapat disimpulkan sebagai berikut 1. Model regresi linear yang dihasilkan dari proses stepwise (model2) memiliki performa yang lebih baik dibandingnkan model regresi linear menggunakan full prediktor (model1).
Beberapa point uji asumsi tidak terpenuhi (residu tidak berdistribusi normal, hererosscesdacity) sehinggsa interpretasi dari model bisa bersifat bias, namun model secara umum model regresi linear ini masih bisa digunakan untuk prediksi.
Saran perbaikan model agar bisa memenuhi asumsi regresi linear bisa dicoba transformasi data, baik di target variabel maupun variabel prediktor. Selain itu untuk permasalahan ini bisa mencoba pendekatan model lain seperti logistic regression, svm, ann