Prediksi Biaya Asuransi Kesehatan Menggunakan Model Linear Regression

1. Import Library

Berikut merupakan packages yang digunakan pada analisis data insurance

library(dplyr) 
library(GGally) 
library(stringr)
library(lmtest)
library(car) 
library(caret)
library(performance) 
library(tidyr)
library(rmdformats)
library(ggplot2)
library(psych)

2. Read Data

insurance <- read.csv("data input/insurance.csv", stringsAsFactors = T)
rmarkdown::paged_table(insurance)

3. Data Wrangling

glimpse(insurance)

#> Rows: 1,338
#> Columns: 7
#> $ age      <int> 19, 18, 28, 33, 32, 31, 46, 37, 37, 60, 25, 62, 23, 56, 27, 1…
#> $ sex      <fct> female, male, male, male, male, female, female, female, male,…
#> $ bmi      <dbl> 27.900, 33.770, 33.000, 22.705, 28.880, 25.740, 33.440, 27.74…
#> $ children <int> 0, 1, 3, 0, 0, 0, 1, 3, 2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0…
#> $ smoker   <fct> yes, no, no, no, no, no, no, no, no, no, no, yes, no, no, yes…
#> $ region   <fct> southwest, southeast, southeast, northwest, northwest, southe…
#> $ charges  <dbl> 16884.924, 1725.552, 4449.462, 21984.471, 3866.855, 3756.622,…

Description
- age : umur nasabah
- sex : jenis kelamin nasabah
- bmi : index massa tubuh
- children : jumlah tanggungan anak nasabah
- smoker : apakah nasabah perokok atau tidak
- region : region tempat tinggal ansabah
- charges : biaya pengobatan (medical cost) yang ditanggung oleh pihak asuransi terhadap nasabah

Inspeksi data unique

unique(insurance$children)

#> [1] 0 1 3 2 5 4

Menghapus data duplikat

insurance <- insurance[!duplicated(insurance),]

Inspeksi missing value

is.na(insurance) %>% colSums()

#>      age      sex      bmi children   smoker   region  charges 
#>        0        0        0        0        0        0        0

Inspeksi range variabel target

range(insurance$charges)

#> [1]  1121.874 63770.428

4. Business Problem

Kita ingin memprediksi biaya asuransi charges berdasarkan seluruh prediktor yang ada di dataset. * variable target (y): charges * variable prediktor (x): age, sex, bmi, children, smoker, region

5. Exploratory Data Analysis (EDA)

Mari kita cek persebaran data numerik pada dataset

pairs.panels(insurance[-c(2,5,6)])

Plot di atas menunjukkan visualisasi secara garis besar persebaran data seluruh variabel numerik pada dataset insurance. Mari kita bedah satu per satu sebagai berikut.

hist(insurance$age)

range(insurance$age)

#> [1] 18 64

💡 Insight:
* Pihak asuransi memiliki nasabah dengan rentang umur 18 - 64 tahun, dan didominasi oleh nasabah yang berumur < 20 tahun

hist(insurance$bmi)

range(insurance$bmi)

#> [1] 15.96 53.13

💡 Insight:
* Nilai bmi nasabah berada pada rentang 15.96 - 53.13, dimana bmi 25 - 35 adalah yang dominan. Berdasarkan klasifikasi obesitas, rentang bmi 25 - 35 masuk ke dalam kelas overweight hingga obesitas level 1

hist(insurance$charges)

range(insurance$charges)

#> [1]  1121.874 63770.428

💡 Insight:
* Biaya pengobatan yang dikeluarkan oleh pihak asuransi berada pada rentang 1121.874 - 63770.428, dimana nilai yang dominan ada di < 15000an

boxplot(insurance$charges, horizontal = T)

💡 Insight:

Variabel charges memiliki outlier atas (outlier di nilai yang besar)
Persebaran data kita ada di nilai kecil (0-30000an)

Mari kita lihat persebaran data charges nasabah yang memiliki nilai bmi < 30 dan > 30.

# subset data dengan bmi<30
bmi_less30 <- insurance[insurance$bmi<30,]
hist(bmi_less30$charges)

💡 Insight:
* Persebaran data charges pada nasabah dengan bmi < 30 cukup stabil

# subset data dengan bmi < 30
bmi_more30 <- insurance[insurance$bmi>30,]
hist(bmi_more30$charges)

💡 Insight:
* Persebaran data charges pada nasabah dengan bmi > 30 cenderung fluktuatif

# subset data dengan bmi 25-35
bmi_range <- insurance[insurance$bmi>25 & insurance$bmi<35,]
hist(bmi_range$charges)

💡 Insight:
* Persebaran data charges pada nasabah dengan bmi 25-35 cukup stabil namun sedikit fluktuatif

Mari kita lihat persebaran data smoker terhadap biaya charges

plot(insurance$smoker, insurance$charges, horizontal=T)

💡 Insight:
* Nasabah perokok memiliki range charges yang lebih tinggi yaitu 12000-60000an dibandingkan bukan perokok yakni 0-20000an

Mari kita lihat persebaran data age terhadap biaya charges.

plot(insurance$age, insurance$charges)

💡 Insight:
* Dari hasil plot di atas, diperoleh informasi yang menarik bahwa bertambahnya umur nasabah, maka nilai charges nya juga bertambah tinggi.

Mari kita lihat persebaran data sex terhadap biaya charges.

plot(insurance$sex, insurance$charges, horizontal=T)

💡 Insight:
* Dari plot di atas, diperoleh informasi bahwa pada dataset kita terdapat jumlah nasabah pria yang lebih banyak dibandingkan wanita, dimana nasabah pria memiliki charges yang lebih tinggi dibandingkan wanita yaitu 0-40000an terhadap 0-20000an, meskipun ada outlier biaya charges pada wanita lebih tinggi dibandingkan pria. Namun, kita perlu fokus terhadap range boxplot data yang menjadi kumpulan data terbanyak.

Mari kita lihat persebaran data region terhadap biaya charges.

plot(insurance$region, insurance$charges, horizontal=T)

💡 Insight:
* Kita dapat melihat bahwa data terbanyak terdapat pada nasabah yang tinggal di region southeast dengan biaya charges berada pada rentang 0-40000an. Sementara, data paling sedikit terdapat pada nasabah yang tinggal di region southwest dengan biaya charges berada pada rentang 0-20000an.

6. Model Building

Simple Linear Regression

Regresi linear sederhana adalah model dengan satu variabel prediktor. Berikut adalah model dengan menggunakan satu prediktor yaitu age.

model_age <- lm(formula=charges~age, data=insurance)

summary(model_age)

#> 
#> Call:
#> lm(formula = charges ~ age, data = insurance)
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#>  -8064  -6684  -5943   5466  47828 
#> 
#> Coefficients:
#>             Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)  3190.02     938.40   3.399             0.000695 ***
#> age           257.23      22.53  11.419 < 0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 11560 on 1335 degrees of freedom
#> Multiple R-squared:  0.08899,    Adjusted R-squared:  0.08831 
#> F-statistic: 130.4 on 1 and 1335 DF,  p-value: < 0.00000000000000022

ggplot(data = insurance, mapping=aes(x = age, y = charges))+
  geom_point()+
  scale_color_brewer(palette="Accent")+
  geom_smooth(method="lm", se=FALSE)+
  labs(title="Plot Regresi Age vs Biaya Charges")+
  theme_minimal()

💡 Insight:

Dengan Menguji Signifikansi Parameter, variabel charges memiliki nilai p-value 2.2e-16 < 0.05, sehingga disimpulkan bahwa age (umur) memiliki pengaruh yang signifikan terhadap besarnya biaya charges asuransi terhadap nasabah
Adj R-Squared = 0.08831, dengan kata lain hanya sebesar 8,83% dari variabel age dapat menjelaskan model (pengaruh biaya asuransi)

Walaupun secara signifikansi parameter, variabel age berpengaruh signifikan, akan tetapi dari nilai adjR2 nilai yang dihasilkan masih sangat kecil, sehingga perlu penambahan faktor variabel lainnya.

Multiple Linear Regression

Akan dibuat model dengan keseluruhan variabel, dengan variabel target adalah charges dan variabel prediktor antara lain age, sex, bmi, children, smoker, region.

model_all <- lm(formula=charges~.,
                data=insurance)
summary(model_all)

#> 
#> Call:
#> lm(formula = charges ~ ., data = insurance)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -11305.1  -2850.3   -979.9   1395.0  29992.8 
#> 
#> Coefficients:
#>                  Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)     -11936.56     988.23 -12.079 < 0.0000000000000002 ***
#> age                256.76      11.91  21.555 < 0.0000000000000002 ***
#> sexmale           -129.48     333.20  -0.389             0.697630    
#> bmi                339.25      28.61  11.857 < 0.0000000000000002 ***
#> children           474.82     137.90   3.443             0.000593 ***
#> smokeryes        23847.33     413.35  57.693 < 0.0000000000000002 ***
#> regionnorthwest   -349.23     476.82  -0.732             0.464053    
#> regionsoutheast  -1035.27     478.87  -2.162             0.030804 *  
#> regionsouthwest   -960.08     478.11  -2.008             0.044836 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 6064 on 1328 degrees of freedom
#> Multiple R-squared:  0.7507, Adjusted R-squared:  0.7492 
#> F-statistic:   500 on 8 and 1328 DF,  p-value: < 0.00000000000000022

💡 Hasil interpretasi dari summary(model_all):

Variabel yang signifikan mempengaruhi charges adalah age, bmi, children, smoker, region
Ketika variabel age naik 1 satuan, maka charges naik sebesar 256.76 dengan variabel prediktor lainnya tetap
Ketika variabel bmi naik 1 satuan, maka charges naik sebesar 339.25 dengan variabel prediktor lainnya tetap
sexfemale tidak masuk sebagai prediktor, artinya dijadikan basis, koefisien sexmale = -129.48, artinya pria menurunkan charges sebesar 129.48 apabila dibandingan dengan wanita
smokerno dijadikan basis, koefisien smokeryes = 23847.33, artinya perokok menaikkan charges sebesar 23847.33 apabila dibandingkan bukan perokok
regionnortheast dijadikan basis, koefisien regionnorthwest = -349.23, artinya regionnorthwest menurunkan charges sebesar 349.23 apabila dibandingkan regionnortheast
regionnortheast dijadikan basis, koefisien regionsoutheast = -1035.27, artinya regionsoutheast menurunkan charges sebesar 1035.27 apabila dibandingkan regionnortheast
regionnortheast dijadikan basis, koefisien regionsouthwest = -960.08, artinya regionsouthwest menurunkan charges sebesar 960.08 apabila dibandingkan regionnortheast
Adjusted R-squared: 0.7492, artinya model kita berhasil menangkap informasi/varians sebesar 74.92% untuk memprediksi variabel target

7. Evaluasi Model Keseluruhan

Performa model

Kinerja model (seberapa baik model memprediksi variabel target) dapat dihitung dengan menggunakan nilai RMSE (Root Mean Square Error) dengan formula berikut:

\[ RMSE = \sqrt{\frac{1}{n} \sum (\hat prediction - actual)^2} \] RMSE merupakan suatu nilai pengujian performa yang lebih baik dibandingkan MAE (Mean Absolute Error). Hal ini disebabkan nilai RMSE mengkuadratkan selisih antara nilai aktual dan nilai prediksi, artinya prediksi dengan error yang lebih tinggi akan dikenakan sanksi yang besar.

# prediksi
insurance$pred_charges <- predict(model_all, insurance)

# menghitung RMSE
sqrt(mean((insurance$pred_charge - insurance$charges)^2))

#> [1] 6043.85

range(insurance$charges)

#> [1]  1121.874 63770.428

Berdasarkan hasil ini, model_all masih memberikan kesalahan yang cukup tinggi untuk prediksi dan oleh karena itu masih perlu perbaikan. Selanjutnya, kita akan melakukan pengujian asumsi terhadap model_all.

Pengujian Asumsi

Uji asumsi terdiri dari 5 yaitu:
- Linearitas
- Normalitas Residual
- Autokorelasi
- Homoskedastisitas Residual
- Multikolinearitas

Note: Hanya perlu cek asumsi kepada model yang akan digunakan, tidak perlu semua kombinasi model dicek asumsinya. Pada bagian ini, kita akan cek asumsi untuk model_all

Linearitas

Pengecekan asumsi linearitas sudah diperiksa sebelum membangun model menggunakan matriks korelasi pada langkah EDA dan juga scatter plot yang dibuat menggunakan variabel prediktor dan target. Akan tetapi akan dilakukan pemeriksaan cor.test() pada salah satu variabel prediktor dengan target.

Hipotesis

H0 : tidak ada korelasi yang signifikan (p=0)

H1 : terdapat korelasi signifikan (p=!0)

α= 0.05

Statistik Uji: Pearson dengan nilai p-value

cor.test(insurance$bmi, insurance$charges)

#> 
#>  Pearson's product-moment correlation
#> 
#> data:  insurance$bmi and insurance$charges
#> t = 7.3961, df = 1335, p-value = 0.0000000000002468
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  0.1463465 0.2493595
#> sample estimates:
#>       cor 
#> 0.1984008

Kesimpulan: nilai p-value = 0.0000000000002468 < alpha = 0.95, sehingga H0 ditolak. Sehingga terdapat korelasi yang signifikan antara variabel bmi dan charges.

Selain itu, akan dilakukan pengecekan asumsi linearitas dengan membuat plot residual vs fitted untuk model_all.

resfit <- data.frame(residual = model_all$residuals, fitted = model_all$fitted.values)

resfit %>% ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0, colour= "red")) +
  theme_minimal() + 
  labs(title = "Residual vs Fitted Plot") + 
  theme(legend.position = "none")

💡 Insight:
* Dari plot di atas dapat dilihat bahwa residual/error tidak memantul secara acak di sekitar 0 (tidak seragam) sehingga akan menghasilkan rata-rata tidak sama dengan 0. Maka syarat E(ϵ)=0 tidak terpenuhi. Disimpulkan bahwa data kita mungkin tidak linear, sehingga asumsi tidak terpenuhi.

Normalitas Residual

Berikut adalah plot histogram dan normal probability plot (normal Q-Q plot) untuk melihat distribusi dari error.

par(mfrow=c(1,2))
hist(model_all$residuals, breaks = 10)
qqnorm(residuals(model_all), main="Normal Q-Q Plot")
qqline(residuals(model_all))

💡 Insight:
* Dapat dilihat dari plot, bahwa terdapat beberapa residual yang tidak berada pada garis. Hal itu mengindikasikan bahwa mungkin saja distribusi dari eror tidak normal. Untuk lebih meyakinkan, akan dilakukan pengujian normalitas dengan menggunakan shapiro wilk test.

Hipotesis:

H0: error berdistribusi normal

H1: error TIDAK berdistribusi normal

α= 0.05

Statistik Uji: Shapiro Wilk dengan nilai p-value

# shapiro test dari residual
shapiro.test(model_all$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_all$residuals
#> W = 0.89909, p-value < 0.00000000000000022

karena p-value < 0.05, sehingga kita tolak H0 atau asumsi normality tidak terpenuhi.

Autokorelasi

Hipotesis

H0 : tidak ada autokorelasi

H1 : ada auatokorelasi

α= 0.05

Statistik Uji: Durbin Watson dengan nilai p-value

durbinWatsonTest(model_all)

#>  lag Autocorrelation D-W Statistic p-value
#>    1     -0.04585675      2.088973   0.106
#>  Alternative hypothesis: rho != 0

Aturan Keputusan Nilai p-value = 0.102 > 0.05 sehingga H0 gagal ditolak

Kesimpulan: Karena H0 gagal ditolak, sehingga tidak terdapat autokorelasi pada model. Asumsi terpenuhi

Homoskedastisitas

Hipotesis

H0: variansi residual seragam (homoskedastisitas)

H1: variansi residual tidak seragam (heteroskedastisitas)

α= 0.05

Statistik Uji: Breusch Pagan dengan nilai p-value

bptest(model_all)

#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  model_all
#> BP = 121.57, df = 8, p-value < 0.00000000000000022

Aturan Keputusan Nilai p-value = 0.00000000000000022 < 0.05 sehingga H0 ditolak

Kesimpulan Karena H0 ditolak, sehingga residual model yang kita miliki memiliki variansi yang tidak seragam (heteroskedastisitas). Asumsi tidak terpenuhi.

Multikolinearitas

Multikolinearitas adalah kondisi adanya korelasi antar prediktor yang kuat. Hal ini tidak diinginkan karena menandakan prediktor redundan pada model, yang seharusnya dapat dipilih salah satu saja dari variable yang hubungannya amat kuat tersebut. Harapannya tidak terjadi multikolinearitas.

Kondisi yang diharapkan: VIF < 10

vif(model_all)

#>              GVIF Df GVIF^(1/(2*Df))
#> age      1.016794  1        1.008362
#> sex      1.008944  1        1.004462
#> bmi      1.106742  1        1.052018
#> children 1.004017  1        1.002006
#> smoker   1.012100  1        1.006032
#> region   1.099037  3        1.015864

Dari uji VIF, prediktor di model backward lolos uji asumsi multikolinearitas (tidak ada nilai VIF > 10)

Kesimpulan:
* Dari kelima uji asumsi di atas, dapat disimpulkan bahwa model_all tidak memenuhi uji Normalitas, uji linearitas dan uji homoskedastisitas. Sehingga, perlu dilukan improvement pada model.

8. Model Improvement

Dari analisis residual yang dilakukan, disimpulkan bahwa residual dari model_all tidak memenuhi uji asumsi. Dengan demikian akan dilakukan transformasi data untuk variabel prediktor dan respon. Salah satu jenis transformasi yang dapat dilakukan adalah transformasi logaritmik.

Berdasarkan hasil summary model_all diketahui bahwa nilai p-value terbesar terdapat pada variabel sex, namun karena variabel sex merupakan tipe data kategorik sehingga tidak dapat dilakukan logaritmik, kemudian akan dipertimbangkan untuk melakukan transformasi logaritmik pada variabel target charges karena nilai nya yang terlalu besar.

Tuning Model

Transformasi dengan Komponen Logaritmik

trans_model_all <- lm(formula=log(charges)~., data=insurance)
summary(trans_model_all)

#> 
#> Call:
#> lm(formula = log(charges) ~ ., data = insurance)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.07416 -0.19619 -0.04969  0.06655  2.16494 
#> 
#> Coefficients: (1 not defined because of singularities)
#>                   Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)      7.0314797  0.0723842  97.141 < 0.0000000000000002 ***
#> age              0.0345390  0.0008725  39.585 < 0.0000000000000002 ***
#> sexmale         -0.0745638  0.0244054  -3.055              0.00229 ** 
#> bmi              0.0134013  0.0020957   6.395       0.000000000222 ***
#> children         0.1015405  0.0101005  10.053 < 0.0000000000000002 ***
#> smokeryes        1.5537619  0.0302763  51.319 < 0.0000000000000002 ***
#> regionnorthwest -0.0620490  0.0349257  -1.777              0.07586 .  
#> regionsoutheast -0.1573100  0.0350753  -4.485       0.000007923624 ***
#> regionsouthwest -0.1289664  0.0350196  -3.683              0.00024 ***
#> pred_charges            NA         NA      NA                   NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.4442 on 1328 degrees of freedom
#> Multiple R-squared:  0.7676, Adjusted R-squared:  0.7662 
#> F-statistic: 548.4 on 8 and 1328 DF,  p-value: < 0.00000000000000022

Dari trans_model_all didapatkan persamaan regresi sebagai berikut:

ln(charges)=7.0314 + 0.03453(age) - 0.0745(sexmale) + 0.0134(bmi) + 0.1015(children) + 1.5537(smokeryes)- 0.06204(regionnorthwest) - 0.1573(regionsoutheast) - 0.1289(regionsouthwest)

Berdasarkan summary model_all, nilai p-value dari seluruh variable prediktor < 0.05 sehingga dapat disimpulkan bahwa seluruh variabel tersebut memiliki pengaruh terhadap variabel target.

trans_model_all <- lm(formula=log(charges)~., data=insurance)
summary(trans_model_all)

#> 
#> Call:
#> lm(formula = log(charges) ~ ., data = insurance)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.07416 -0.19619 -0.04969  0.06655  2.16494 
#> 
#> Coefficients: (1 not defined because of singularities)
#>                   Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)      7.0314797  0.0723842  97.141 < 0.0000000000000002 ***
#> age              0.0345390  0.0008725  39.585 < 0.0000000000000002 ***
#> sexmale         -0.0745638  0.0244054  -3.055              0.00229 ** 
#> bmi              0.0134013  0.0020957   6.395       0.000000000222 ***
#> children         0.1015405  0.0101005  10.053 < 0.0000000000000002 ***
#> smokeryes        1.5537619  0.0302763  51.319 < 0.0000000000000002 ***
#> regionnorthwest -0.0620490  0.0349257  -1.777              0.07586 .  
#> regionsoutheast -0.1573100  0.0350753  -4.485       0.000007923624 ***
#> regionsouthwest -0.1289664  0.0350196  -3.683              0.00024 ***
#> pred_charges            NA         NA      NA                   NA    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.4442 on 1328 degrees of freedom
#> Multiple R-squared:  0.7676, Adjusted R-squared:  0.7662 
#> F-statistic: 548.4 on 8 and 1328 DF,  p-value: < 0.00000000000000022

Model Backward

Selanjutnya, akan dilakukan backward regression dengan tujuan menemukan model terbaik pada trans_model_all. Dilakukan metode backward regression untuk pemilihan model menggunakan kriteria AIC. Backward regression ialah metode yang akan mengeluarkan variabel bebas yang tidak berpengaruh signifikan terhadap variabel respon secara satu per satu dilihat dari nilai p-value masing-masing variabel bebas, variabel bebas yang pertama dikeluarkan ialah yang memiliki p-value paling besar dan melebihi taraf signifikansi.

set.seed(123)
trans_model_backward <- step(object=trans_model_all,
                             direction="backward",
                             trace=F)

summary(trans_model_backward)

#> 
#> Call:
#> lm(formula = log(charges) ~ age + sex + bmi + children + smoker + 
#>     region, data = insurance)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.07416 -0.19619 -0.04969  0.06655  2.16494 
#> 
#> Coefficients:
#>                   Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)      7.0314797  0.0723842  97.141 < 0.0000000000000002 ***
#> age              0.0345390  0.0008725  39.585 < 0.0000000000000002 ***
#> sexmale         -0.0745638  0.0244054  -3.055              0.00229 ** 
#> bmi              0.0134013  0.0020957   6.395       0.000000000222 ***
#> children         0.1015405  0.0101005  10.053 < 0.0000000000000002 ***
#> smokeryes        1.5537619  0.0302763  51.319 < 0.0000000000000002 ***
#> regionnorthwest -0.0620490  0.0349257  -1.777              0.07586 .  
#> regionsoutheast -0.1573100  0.0350753  -4.485       0.000007923624 ***
#> regionsouthwest -0.1289664  0.0350196  -3.683              0.00024 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.4442 on 1328 degrees of freedom
#> Multiple R-squared:  0.7676, Adjusted R-squared:  0.7662 
#> F-statistic: 548.4 on 8 and 1328 DF,  p-value: < 0.00000000000000022

Berdasarkan hasil model step backward, diperoleh bahwa tidak terjadi pengurangan variabel prediktor.

Model Forward

Untuk dapat melakukan pemodelan step forward, pertama kita perlu melakukan model tanpa prediktor seperti berikut ini.

trans_model_nopredictor <- lm(formula=log(charges)~1, data=insurance)

Setelah model tanpa prediktor dibuat, mari lakukan model step forward

set.seed(123)
trans_model_forward <- step(object=trans_model_nopredictor,
                      direction="forward",
                      scope=list(lower=trans_model_nopredictor, upper=trans_model_all),
                      trace=F)

summary(trans_model_forward)

#> 
#> Call:
#> lm(formula = log(charges) ~ pred_charges + age + children + bmi + 
#>     smoker + sex, data = insurance)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.07838 -0.19878 -0.04986  0.06552  2.15866 
#> 
#> Coefficients:
#>                 Estimate  Std. Error t value             Pr(>|t|)    
#> (Intercept)   8.70846907  0.35654448  24.425 < 0.0000000000000002 ***
#> pred_charges  0.00014064  0.00002902   4.846           0.00000141 ***
#> age          -0.00156135  0.00752835  -0.207             0.835732    
#> children      0.03479435  0.01705488   2.040             0.041534 *  
#> bmi          -0.03442271  0.00956830  -3.598             0.000333 ***
#> smokeryes    -1.80068913  0.69201373  -2.602             0.009368 ** 
#> sexmale      -0.05631231  0.02466837  -2.283             0.022601 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.4439 on 1330 degrees of freedom
#> Multiple R-squared:  0.7676, Adjusted R-squared:  0.7665 
#> F-statistic:   732 on 6 and 1330 DF,  p-value: < 0.00000000000000022

Model Both

set.seed(123)
trans_model_both <- step(object=trans_model_nopredictor,
                      direction="both",
                      scope=list(upper=trans_model_all),
                      trace=F)

summary(trans_model_both)

#> 
#> Call:
#> lm(formula = log(charges) ~ pred_charges + children + bmi + smoker + 
#>     sex, data = insurance)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.07810 -0.20000 -0.05044  0.06386  2.15701 
#> 
#> Coefficients:
#>                  Estimate   Std. Error t value             Pr(>|t|)    
#> (Intercept)   8.635774819  0.065303360 132.241 < 0.0000000000000002 ***
#> pred_charges  0.000134659  0.000003357  40.108 < 0.0000000000000002 ***
#> children      0.037617670  0.010270102   3.663             0.000259 ***
#> bmi          -0.032500687  0.002379846 -13.657 < 0.0000000000000002 ***
#> smokeryes    -1.658250130  0.084789200 -19.557 < 0.0000000000000002 ***
#> sexmale      -0.057058873  0.024395548  -2.339             0.019488 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.4438 on 1331 degrees of freedom
#> Multiple R-squared:  0.7675, Adjusted R-squared:  0.7667 
#> F-statistic:   879 on 5 and 1331 DF,  p-value: < 0.00000000000000022

Performa Model

Untuk semua model kandidat dilakukan perbandingan performa untuk semua model tersebut. Tujuannya adalah mencari model regresi terbaik untuk memodelkan data.

model_comparison <- compare_performance(trans_model_nopredictor, trans_model_all, trans_model_backward, trans_model_forward, trans_model_both)
as.data.frame(model_comparison) %>% 
  select("Name","AIC","BIC","R2","R2_adjusted","RMSE")

#>                      Name      AIC      BIC        R2 R2_adjusted      RMSE
#> 1 trans_model_nopredictor 27903.70 27914.10 0.0000000   0.0000000 0.9183550
#> 2         trans_model_all 25968.45 26020.43 0.7676303   0.7662305 0.4426907
#> 3    trans_model_backward 25968.45 26020.43 0.7676303   0.7662305 0.4426907
#> 4     trans_model_forward 25964.88 26006.46 0.7675559   0.7665073 0.4427615
#> 5        trans_model_both 25962.92 25999.31 0.7675484   0.7666751 0.4427687

Dapat dilihat bahwa model yang terbaik ialah trans_model_both karena memiliki tiga performa kriteria yang lebih unggul dibandingkan performa kriteria pada model yang lain, yaitu nilai Adj R-Square terbesar, nilai AIC terkecil, dan nilai BIC terkecil. Sehingga disimpulkan bahwa model regresi terbaik untuk memodelkan data pada kasus ini ialah trans_model_both. Berikut detailnya:
* RMSE terkecil: trans_model_all, trans_model_backward
* Adj R-Square terbesar: trans_model_both
* AIC terkecil: trans_model_both
* BIC terkecil: trans_model_both

Dari model trans_model_both terdapat pengurangan variabel prediktor sehingga yang tersisa adalah children, bmi, smoker, dan sex. Keempat variabel tersebut kemudian dimasukkan ke dalam step_both2.

step_both2 <- lm(formula=log(charges) ~ children + bmi + smoker + sex, data=insurance)
summary(step_both2)

#> 
#> Call:
#> lm(formula = log(charges) ~ children + bmi + smoker + sex, data = insurance)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.99441 -0.43829  0.05817  0.45328  1.74731 
#> 
#> Coefficients:
#>             Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)  8.10434    0.09499  85.313 < 0.0000000000000002 ***
#> children     0.11790    0.01496   7.879  0.00000000000000683 ***
#> bmi          0.01970    0.00296   6.655  0.00000000004125752 ***
#> smokeryes    1.52011    0.04480  33.933 < 0.0000000000000002 ***
#> sexmale     -0.09735    0.03621  -2.688              0.00727 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.6592 on 1332 degrees of freedom
#> Multiple R-squared:  0.4866, Adjusted R-squared:  0.4851 
#> F-statistic: 315.6 on 4 and 1332 DF,  p-value: < 0.00000000000000022

Pengujian Asumsi

Linearitas

resfit1 <- data.frame(residual = step_both2$residuals, fitted = step_both2$fitted.values)

resfit1 %>% ggplot(aes(fitted, residual)) + 
  geom_point() + 
  geom_hline(aes(yintercept = 0, colour= "red")) +
  theme_minimal() + 
  labs(title = "Residual vs Fitted Plot") + 
  theme(legend.position = "none")

💡 Insight:
* Dari plot di atas dapat dilihat bahwa residual/error cukup memantul secara acak di sekitar 0 (seragam) sehingga akan menghasilkan rata-rata dengan 0. Maka syarat E(ϵ)=0 terpenuhi. Disimpulkan bahwa model kita linear. Asumsi terpenuhi.

Normalitas Residual

Berikut adalah plot histogram dan normal probability plot (normal Q-Q plot) untuk melihat distribusi dari error.

par(mfrow=c(1,2))
hist(step_both2$residuals, breaks = 10)
qqnorm(residuals(step_both2), main="Normal Q-Q Plot")
qqline(residuals(step_both2))

💡 Insight:
* Jika dibandingkan dengan plot normalitas residual pada model_all, distribusi error nya jauh lebih baik, namun masih terdapat beberapa residual yang tidak berada pada garis. Untuk lebih meyakinkan, akan dilakukan pengujian normalitas dengan menggunakan shapiro wilk test.

Hipotesis:

H0: error berdistribusi normal

H1: error TIDAK berdistribusi normal

α= 0.05

Statistik Uji: Shapiro Wilk dengan nilai p-value

#Shapiro wilk test
shapiro.test(step_both2$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  step_both2$residuals
#> W = 0.99525, p-value = 0.0003208

Aturan Keputusan p-value > 0.05; H0 ditolak.

Kesimpulan Maka dari itu, error/residual tidak berdistribusi normal. Asumsi normalitas residual tidak terpenuhi.

Autokorelasi

Hipotesis

H0 : tidak ada autokorelasi

H1 : ada auatokorelasi

α= 0.05

Statistik Uji: Durbin Watson dengan nilai p-value

durbinWatsonTest(step_both2)

#>  lag Autocorrelation D-W Statistic p-value
#>    1    -0.007029496      2.013713    0.77
#>  Alternative hypothesis: rho != 0

Aturan Keputusan Nilai p-value = 0.758 > 0.05 sehingga H0 gagal ditolak

Kesimpulan: Karena H0 gagal ditolak, sehingga tidak terdapat autokorelasi pada model. Asumsi terpenuhi

Homoskedastisitas

Hipotesis

H0: variansi residual seragam (homoskedastisitas)

H1: variansi residual tidak seragam (heteroskedastisitas)

α= 0.05

Statistik Uji: Breusch Pagan dengan nilai p-value

bptest(step_both2)

#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  step_both2
#> BP = 214.66, df = 4, p-value < 0.00000000000000022

Aturan Keputusan Nilai p-value 2.2e-16 < 0.05 sehingga H0 ditolak

Kesimpulan: Karena H0 ditolak, sehingga residual model yang kita miliki memiliki variansi yang tidak seragam (heteroskedastisitas). Asumsi tidak terpenuhi.

Multikolinearitas

vif(step_both2)

#> children      bmi   smoker      sex 
#> 1.000497 1.002300 1.005938 1.008341

Kesimpulan:: Dari output diketahui bahwa semua nilai VIF < 10. Dengan demikian tidak terdapat gejala multikolinearitas pada masing-masing variabel bebas yang terdapat pada step_both2. Asumsi terpenuhi

9. Interpretasi Model Akhir

Dari pengujian asumsi step_both2 diketahui bahwa terdapat beberapa asumsi yang tidak terpenuhi yaitu uji normalitas dan homoskedastisitas. Sehingga, model tidak bisa digunakan di dataset untuk penelitian, namun masih dapat digunakan untuk keperluan bisnis. Selanjutnya akan dilakukan interpretasi model berikut:

ln(charges)= 8.10434 + 0.11790(children) + 0.01970(bmi) + 1.52011(smokeryes) - 0.09735(sex)

Apabila variabel jumlah tanggungan anak bertambah 1 satuan, maka rata-rata banyaknya biaya yang dikeluarkan oleh pihak asuransi per individu naik sebesar e^(0.11790) = 1.12513 kali dengan variabel lain tetap. Apabila variabel bmi bertambah 1 satuan, maka rata-rata banyaknya biaya yang dikeluarkan oleh pihak asuransi per individu naik sebesar e^(0.01970) = 1.01989 kali dengan variabel lain tetap.

10. Kesimpulan

Berdasarkan proses analisis faktor-faktor pada data banyaknya biaya yang dikeluarkan oleh pihak asuransi per individu, didapatkan hasil bahwa faktor yang berpengaruh secara signifikan terhadap biaya charges yaitu jumlah tanggungan anak, nilai bmi, status perokok dan jenis kelamin.