1. Introduction

Kali ini akan membahas mengenai regression model dari data sebuah perusahaan asuransi untuk melakukan prediksi premi asuransi yang harus dibayarkan oleh seorang pelanggan yang akan melakukan pemilihan asuransi.

1.1 Import Library

Diperlukan langkah import library untuk keperluan pengolahan data atau visualisasi data

#import dplyr untuk data cleansing
library(dplyr) 
#import ggplot2 untuk keperluan visualisasi
library(ggplot2)
#import GGally untuk visualisasi (korelasi)
library(GGally) 
#import MLmetrics untuk Metrics Error
library(MLmetrics) 
library(performance) # Compare performance
library(lmtest) # Uji Asumsi
library(car) # Uji Asumsi

1.2 Input Data

Data diinput dan pastikan data sudah berada pada folder yang sama dengan Rproj.

train_ins <- read.csv("train_insurance.csv")
train_ins

Data yang diinput terdiri dari beberapa column, yaitu: - age: usia dari customer - sex: gender dari customer berupa male atau female - bmi: body mass index customer - children: jumlah anak dari customer - smoker: merupakan seorang perokok atau bukan - region: wilayah tempat tinggal
- charges: biaya medical dalam dolar

2. Exploratory Data Analysis

Exploratory data analysis perlu dilakukan sebelum melanjutkan ke tahap selanjutnya yaitu tahapan permodelan

2.1 Check Data Type

Cek tipe data setiap column/variabel untuk memastikan kebeneran tipe data dengan menggunakan str()

str(train_ins)
#> 'data.frame':    1070 obs. of  7 variables:
#>  $ age     : int  19 62 46 18 18 39 41 62 20 24 ...
#>  $ sex     : chr  "female" "female" "female" "female" ...
#>  $ bmi     : num  35.1 38.1 28.9 33.9 34.4 ...
#>  $ children: int  0 2 2 0 0 5 3 0 0 0 ...
#>  $ smoker  : chr  "no" "no" "no" "no" ...
#>  $ region  : chr  "northwest" "northeast" "southwest" "southeast" ...
#>  $ charges : num  2135 15230 8823 11483 1137 ...

Berdasarkan hasil, masih terdapat tipe data yang kurang tepat sehingga perlu diubah

2.2 Change Data Type

Tipe data yang masih salah yaitu sex, smoker, dan region. Ketiga tersebut tadinya mmiliki tipe data chr, akan tetapi ketiga variabel tersebut diubah menjadi tipe data factor karena terdapat banyak jawaban pengulangan dan sangat spesifik

train_ins <- train_ins %>% 
  mutate(sex = as.factor(sex),
         smoker = as.factor(smoker),
         region = as.factor(region))

str(train_ins)
#> 'data.frame':    1070 obs. of  7 variables:
#>  $ age     : int  19 62 46 18 18 39 41 62 20 24 ...
#>  $ sex     : Factor w/ 2 levels "female","male": 1 1 1 1 2 1 1 2 2 2 ...
#>  $ bmi     : num  35.1 38.1 28.9 33.9 34.4 ...
#>  $ children: int  0 2 2 0 0 5 3 0 0 0 ...
#>  $ smoker  : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 2 2 ...
#>  $ region  : Factor w/ 4 levels "northeast","northwest",..: 2 1 4 3 3 2 4 3 4 1 ...
#>  $ charges : num  2135 15230 8823 11483 1137 ...

Data type sudah berhasil dibenarkan tipe datanya, sehingga dapat dilanjutkan langkah berikutnya

2.3 Check Missing Value

Missing value pada setiap column di data set perlu untuk dilakukan

colSums(is.na(train_ins))
#>      age      sex      bmi children   smoker   region  charges 
#>        0        0        0        0        0        0        0

Tidak terdapat missing value pada setiap column yang ditunjukan oleh angka 0 pada setiap column

2.4 Check Correlation

Melakukan pengecekan korelasi seluruh variabel dengan chargs

ggcorr(train_ins, label = T)

Berdasarkan output korelasi diatas Seluruh variabel memiliki korelasi yang cukup lemah dengan charges

#3.Making a Model

3.1 Model With All Predictor

Pembuatan model yang dilakukan dengan semua prediktor

model_all <- lm(formula = charges ~., data = train_ins)
summary(model_all)
#> 
#> Call:
#> lm(formula = charges ~ ., data = train_ins)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -11297  -2846  -1005   1542  29791 
#> 
#> Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)     -12031.63    1132.08 -10.628  < 2e-16 ***
#> age                243.43      13.65  17.834  < 2e-16 ***
#> sexmale           -216.88     378.52  -0.573 0.566792    
#> bmi                355.06      32.62  10.884  < 2e-16 ***
#> children           559.84     156.28   3.582 0.000356 ***
#> smokeryes        24172.29     459.97  52.552  < 2e-16 ***
#> regionnorthwest   -556.52     539.58  -1.031 0.302594    
#> regionsoutheast   -856.23     541.73  -1.581 0.114281    
#> regionsouthwest   -967.82     541.96  -1.786 0.074421 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 6145 on 1061 degrees of freedom
#> Multiple R-squared:  0.7533, Adjusted R-squared:  0.7514 
#> F-statistic: 404.9 on 8 and 1061 DF,  p-value: < 2.2e-16

Tahapan pembuatan model tlah brhasil dilakukan, maka tahapn selanjutnya dapat dilakukan

3.2 Stepwise Regression

Stepwise regression digunakan dengan metode direction backward untuk mengetahui variabel apa saja yang signifikan berpengaruh terhadap charges

model_backward <- step(model_all, direction = "backward", trace = 0)
summary(model_backward)
#> 
#> Call:
#> lm(formula = charges ~ age + bmi + children + smoker, data = train_ins)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -11858.4  -2937.3   -936.3   1585.6  29425.9 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -12351.07    1080.68 -11.429  < 2e-16 ***
#> age            244.00      13.64  17.893  < 2e-16 ***
#> bmi            341.76      31.10  10.990  < 2e-16 ***
#> children       555.40     156.16   3.557 0.000392 ***
#> smokeryes    24162.30     458.29  52.723  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 6145 on 1065 degrees of freedom
#> Multiple R-squared:  0.7523, Adjusted R-squared:  0.7514 
#> F-statistic: 808.8 on 4 and 1065 DF,  p-value: < 2.2e-16

3.3 Check Perform Two Models

  • Berdasrkan nilai RMSE Periksa performa kedua model di atas untuk meilhat model yang memiliki performa paling baik
compare_performance(model_all, model_backward)

Berdasarkan RMSE, model dengan seluruh prediktor memiliki performa yang lebih baik daripada model backward

  • Berdasarkan nilai MAPE
# Berdasarkan MAPE
MAPE(model_all$fitted.values, train_ins$charges)
#> [1] 0.4256517
MAPE(model_backward$fitted.values, train_ins$charges)
#> [1] 0.4291682

Berdasarkan MAPE, model dengan seluruh prediktor memiliki nilai MAPE yang sedikit lebih rendah jika dibandingkan dengan nilai MAPE dari model backward

4. Interpretation The Bst Model

Melakukan interpretasi untuk model yang memiliki performa paling baik

summary(model_all)
#> 
#> Call:
#> lm(formula = charges ~ ., data = train_ins)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -11297  -2846  -1005   1542  29791 
#> 
#> Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)     -12031.63    1132.08 -10.628  < 2e-16 ***
#> age                243.43      13.65  17.834  < 2e-16 ***
#> sexmale           -216.88     378.52  -0.573 0.566792    
#> bmi                355.06      32.62  10.884  < 2e-16 ***
#> children           559.84     156.28   3.582 0.000356 ***
#> smokeryes        24172.29     459.97  52.552  < 2e-16 ***
#> regionnorthwest   -556.52     539.58  -1.031 0.302594    
#> regionsoutheast   -856.23     541.73  -1.581 0.114281    
#> regionsouthwest   -967.82     541.96  -1.786 0.074421 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 6145 on 1061 degrees of freedom
#> Multiple R-squared:  0.7533, Adjusted R-squared:  0.7514 
#> F-statistic: 404.9 on 8 and 1061 DF,  p-value: < 2.2e-16

Berdasarkan hasil di atas, maka dapat disimpulkan bahwa: > Jika customer tersebut adalah pria maka chargesnya akan berkurang sebesar 216.88 dibandingkan customer wanita

Jika customer anak-anak maka chargesnya bertambah besar 559.84 dibandingka customer bukan anak-anak

Jika customer perokok maka charges bertambah sebesar 24172.29 dibandingkan customer bukan perokok

Jika customer tinggal di northwest maka charges berkurang sebesar 556.52 dibandingkan merekayang tinggal di northeast

Jika customer tinggal di southeast maka charges berkurang sebesar 856.23 dibandingkan merekayang tinggal di northeast

Jika customer tinggal di southwest maka charges berkurang sebesar 967.82 dibandingkan merekayang tinggal di northeast

5. Test Assumption

Mlakukan uji asumsi untuk model yang dibuat dengan salah satu model yang memiliki performance baik

5.1 Linearity

Uji asumsi untuk linearity dengan mnggunakan cor.test()

Hipotesis yang digunakan: H0 : Tidak Linear H1 : Linear

cor.test(train_ins$age, train_ins$charges)
#> 
#>  Pearson's product-moment correlation
#> 
#> data:  train_ins$age and train_ins$charges
#> t = 9.1648, df = 1068, p-value < 2.2e-16
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.2135462 0.3246964
#> sample estimates:
#>       cor 
#> 0.2700206
cor.test(train_ins$bmi, train_ins$charges)
#> 
#>  Pearson's product-moment correlation
#> 
#> data:  train_ins$bmi and train_ins$charges
#> t = 6.7499, df = 1068, p-value = 2.423e-11
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.1440899 0.2590630
#> sample estimates:
#>       cor 
#> 0.2022733
cor.test(train_ins$children, train_ins$charges)
#> 
#>  Pearson's product-moment correlation
#> 
#> data:  train_ins$children and train_ins$charges
#> t = 1.9825, df = 1068, p-value = 0.04768
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.0006231342 0.1200454838
#> sample estimates:
#>        cor 
#> 0.06055099

5.2 Normality of Residual

Uji Statistik dengan menggunakan shapiro.test()

H0: residual berdistribusi normal

H1: residual tidak berdistribusi normal

shapiro.test(model_all$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_all$residuals
#> W = 0.9048, p-value < 2.2e-16

5.3 Homoscedasticity

H0: model homoscedasticity

H1: model heteroscedasticity

plot(train_ins$charges, model_all$residuals)

bptest(model_all)
#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_all
#> BP = 95.768, df = 8, p-value < 2.2e-16

5.4 No Multicollinearity

Uji No Multicollinearity dengan menggunakan vif(). Memeriksa nilai VIF yang dimiliki variabel

vif(model_all)
#>              GVIF Df GVIF^(1/(2*Df))
#> age      1.016119  1        1.008027
#> sex      1.013515  1        1.006735
#> bmi      1.113290  1        1.055126
#> children 1.002896  1        1.001447
#> smoker   1.008439  1        1.004210
#> region   1.099864  3        1.015991

Berdasarkan uji asumsi di atas normality of residual dan Homoscedasticity merupakan asumsi yang dilanggar

6. Prediction Data Test

Langkah selanjutnya adalah memprediksi data test asuransi yang sudah ada menggunakan model terbaik yang telah dibuat, berdasarkan model linear regression yang sudah dibuat,

6.1 Import data test_insurance.csv

Data diinput dan pastikan data sudah berada pada folder yang sama dengan Rproj

test_insurance <- read.csv("test_insurance.csv", stringsAsFactors = T)
test_insurance

6.2 Data Wrangling

Lakukan data wrangling seperti yang dilakukan pada data train

#Melakukan pengecekan tipe data
str(test_insurance)
#> 'data.frame':    268 obs. of  7 variables:
#>  $ age     : int  56 27 60 30 63 19 41 18 36 48 ...
#>  $ sex     : Factor w/ 2 levels "female","male": 1 2 1 1 1 1 2 1 2 2 ...
#>  $ bmi     : num  39.8 42.1 36 32.4 23.1 ...
#>  $ children: int  0 0 0 1 0 5 1 2 1 1 ...
#>  $ smoker  : Factor w/ 2 levels "no","yes": 1 2 1 1 1 1 1 1 2 2 ...
#>  $ region  : Factor w/ 4 levels "northeast","northwest",..: 3 3 1 4 1 4 3 1 3 4 ...
#>  $ charges : num  11091 39612 13229 4150 14452 ...

Terdapat beberapa variable yang masih kurang tepat tipe datanya, yaitu variabel sex, smoker, dan region

# Data wrangling untuk data test (lakukan sama seperti data train)
test_insurance <- test_insurance %>% 
  mutate(sex = as.factor(sex),
         smoker = as.factor(smoker),
         region = as.factor(region))

Check Missing Value Missing value pada setiap column di data set perlu untuk dilakukan

colSums(is.na(test_insurance))
#>      age      sex      bmi children   smoker   region  charges 
#>        0        0        0        0        0        0        0

Tidak terdapat missing value pada setiap column yang ditunjukan oleh angka 0 pada setiap column

6.3 Prediction Data Test

Melakukan prediksi data test menggunakan model yang telah dibuat (model_all)

prediksi <- predict(model_all, test_insurance)

6.4 Checking Error

Melakukan pemeriksaan error (MAPE) dari hasil prediksi data test yang dilakukan sebelumnya

MAPE(prediksi, test_insurance$charges) * 100
#> [1] 45.92472

7. Conclusion

Berdsarkan hasil yang dilakukan, model all memiliki performa yang cukup buruk dalam melakukan prediksi dengan nilai MAPE mencapai 45.92472. Oleh karena itu, perlu dilakukan tahapan-tahapan selanjutnya agar mendapatkan hasil yang lebih baik.

8. Updating Prediction

##8.1 Transformation Data Karena model yang telah dibuat masih belum baik untuk memprediksi nilai charges dan ada asumsi yang dilanggar, sehingga perlu dilakukan beberapa transformasi data pada variabel charges.

predictor_df <- train_ins %>% 
  select(-charges)

# log transformation
log_df <- cbind(predictor_df, log_charges = log(train_ins$charges))
head(log_df)
# sqrt transformation
sqrt_df <- cbind(predictor_df, sqrt_charges = sqrt(train_ins$charges))
head(sqrt_df)

##8.2 Making Some Models Langkah selanjutnya dengan melakukan pemodelan untuk kedua data tersebut

model_log <- lm(log_charges~., data = log_df)
model_sqrt <- lm(sqrt_charges~., data = sqrt_df)

##8.3 Comparing The Models Setelah memperoleh model, performa keduanya dibandingkan dan lihat transformasi yang menghasilkan performa paling baik

compare_performance(model_sqrt, model_log)

Berdasrkan hasil di atas, transformasi log memiliki hasil yang lebih baik daripada sqrt

##8.4 Prediction Data Test Selanjutnya, lakukan prediksi pada data test

log_test <- test_insurance %>% 
  mutate(log_charges = log(charges))

lihat kembaliperformanya meningkat dan bandingkan dengan data yang tidak ditransformasi

prediction_log <- predict(model_log, log_test)
MAPE(prediction_log, log_test$log_charges) * 100
#> [1] 3.175772

Hasil prediksi yang didapatkan menunjukkan hasil yang sudah meningkat dengan nilai 3.175772 dan juga sudah jauh lebih baik dari segi nilai MAPE