Analisis Biaya Asuransi Kesehatan menggunakan Model Linear Regression

1. Input Library

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(GGally)

## Warning: package 'GGally' was built under R version 4.4.1

## Loading required package: ggplot2

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(stringr)
library(lmtest)

## Warning: package 'lmtest' was built under R version 4.4.1

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(tidyr)
library(rmdformats)

## Warning: package 'rmdformats' was built under R version 4.4.1

library(ggplot2)
library(psych)

## Warning: package 'psych' was built under R version 4.4.1

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(car)

## Warning: package 'car' was built under R version 4.4.1

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.4.1

## 
## Attaching package: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

library(caret)

## Warning: package 'caret' was built under R version 4.4.1

## Loading required package: lattice

library(performance)

2. Read Data

insurance<- read.csv("insurance data rv.csv", stringsAsFactors = TRUE)
rmarkdown::paged_table(insurance)

Description - age : umur nasabah - sex : jenis kelamin nasabah - bmi : index massa tubuh - children : jumlah tanggungan anak nasabah - smoker : apakah nasabah perokok atau tidak - Claim Amount : - Past Consultations : - Num of steps: - Hospital Expenditure : - Number of past Hospitalizations : - Anual Salary - region : region tempat tinggal ansabah - charges : biaya pengobatan (medical cost) yang ditanggung oleh pihak asuransi terhadap nasabah

3. Data Wraling

Cek tipe data

glimpse(insurance)

## Rows: 10,008
## Columns: 7
## $ age      <int> 45, 64, 19, 36, 19, 34, 47, 63, 19, 20, 62, 51, 62, 50, 59, 2…
## $ sex      <fct> male, male, female, male, female, male, female, female, femal…
## $ bmi      <dbl> 28.700, 34.500, 32.110, 28.880, 24.605, 34.210, 36.000, 26.22…
## $ children <int> 2, 0, 0, 3, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 3, 0, 0, 4, 3…
## $ smoker   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, yes, …
## $ region   <fct> southwest, southwest, northwest, northeast, northwest, southe…
## $ charges  <dbl> 8027.968, 13822.803, 2130.676, 6748.591, 2709.244, 3935.180, …

Inspeksi data unique

unique(insurance$children)

## [1]  2  0  3  1  4  5 NA

Menghapus data duplikat

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

Inspeksi missing value

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

##      age      sex      bmi children   smoker   region  charges 
##        9        0        3        5        0        0        0

Menghapus baris yang mengandung NA

# Menentukan kolom mana yang tidak memiliki nilai NA
complete_columns <- colSums(is.na(insurance)) == 0

# Subset dataframe hanya dengan kolom-kolom yang tidak memiliki NA
insurance_clean <- subset(insurance, select = complete_columns)

# Menampilkan struktur dataframe baru
str(insurance_clean)

## 'data.frame':    1333 obs. of  4 variables:
##  $ sex    : Factor w/ 2 levels "female","male": 2 2 1 2 1 2 1 1 1 1 ...
##  $ smoker : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ region : Factor w/ 4 levels "northeast","northwest",..: 4 4 2 1 2 3 4 2 4 3 ...
##  $ charges: num  8028 13823 2131 6749 2709 ...

Inspeksi range variable target

range(insurance$charges)

## [1]  1121.874 63770.428

4. Summary Data

summary(insurance)

##       age            sex           bmi           children     smoker    
##  Min.   :18.00   female:660   Min.   :15.96   Min.   :0.000   no :1061  
##  1st Qu.:27.00   male  :673   1st Qu.:26.24   1st Qu.:0.000   yes: 272  
##  Median :39.00                Median :30.33   Median :1.000             
##  Mean   :39.28                Mean   :30.65   Mean   :1.092             
##  3rd Qu.:51.00                3rd Qu.:34.60   3rd Qu.:2.000             
##  Max.   :64.00                Max.   :53.13   Max.   :5.000             
##  NA's   :9                    NA's   :3       NA's   :5                 
##        region       charges     
##  northeast:323   Min.   : 1122  
##  northwest:323   1st Qu.: 4738  
##  southeast:363   Median : 9361  
##  southwest:324   Mean   :13209  
##                  3rd Qu.:16456  
##                  Max.   :63770  
##

5. 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

6. Exploratory Data Analysis (EDA)

Cek persebaran data numerik pada dataset

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

Age Distribution Analysis

hist(insurance$age, col = "lightblue", main = "Histogram of Age")

range(insurance$age, na.rm = TRUE)

## [1] 18 64

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

Bmi Distribution Analysis

hist(insurance$bmi, col = "lightgreen", main = "Histogram of BMI")

range(insurance$bmi, na.rm = TRUE)

## [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

Charges Distribution Analysis

hist(insurance$charges, col = "lightcoral", main = "Histogram of 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 = TRUE, col = "#ffdab9", main = "Boxplot of Insurance Charges")

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

Melihat persebaran data charges nasabah yang memiliki nilai bmi

# subset data dengan bmi<30
bmi_less30 <- insurance[insurance$bmi<30,]
hist(bmi_less30$charges, col = "#BDB76B", main = "Histogram of Charges for BMI < 30")

💡 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, col = "#20B2AA", main = "Histogram of Charges for BMI > 30")

💡 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, col = "#BC8F8F", main = "Histogram of Charges for BMI 25-35")

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

Melihat persebaran data smoker terhadap biaya (charges)

plot(insurance$smoker, insurance$charges, 
     col = ifelse(insurance$smoker == "yes", "#F08080", "#B0C4DE"), 
     pch = 19, xlab = "Smoker Status", ylab = "Charges", 
     main = "Scatter Plot of Charges vs Smoker Status",
     xlim = c(0.8, 2.2), 
     ylim = c(0, 80000))

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

Melihat persebaran data age terhadap biaya (charges)

insurance$smoker <- factor(insurance$smoker, levels = c("no", "yes"))

# Define colors based on smoker status
colors <- ifelse(insurance$smoker == "yes", "#F08080", "#B0C4DE")

# Plotting with color
plot(insurance$age, insurance$charges, col = colors,
     pch = 19, xlab = "Age", ylab = "Charges",
     main = "Scatter Plot of Charges vs Age",
     ylim = c(0, 80000))

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

Melihat 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.

Melihat 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.

7. 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 
##  -8066  -6691  -5865   5252  47886 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3489.21     944.32   3.695 0.000229 ***
## age           249.27      22.64  11.009  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11550 on 1322 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.08398,    Adjusted R-squared:  0.08329 
## F-statistic: 121.2 on 1 and 1322 DF,  p-value: < 2.2e-16

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()

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

💡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.08329, dengan kata lain hanya sebesar 8,32% 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 
## -11227.0  -2854.1   -973.6   1405.1  30078.7 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -11923.07     994.96 -11.983  < 2e-16 ***
## age                252.21      12.04  20.953  < 2e-16 ***
## sexmale            -83.57     335.94  -0.249 0.803586    
## bmi                343.06      28.89  11.873  < 2e-16 ***
## children           485.75     139.07   3.493 0.000494 ***
## smokeryes        23788.23     415.46  57.258  < 2e-16 ***
## regionnorthwest   -319.25     480.43  -0.665 0.506476    
## regionsoutheast  -1032.15     482.06  -2.141 0.032448 *  
## regionsouthwest   -938.87     481.61  -1.949 0.051456 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6066 on 1307 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.7493, Adjusted R-squared:  0.7478 
## F-statistic: 488.3 on 8 and 1307 DF,  p-value: < 2.2e-16

💡Insight: 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 252.21 dengan variabel prediktor lainnya tetap
Ketika variabel bmi naik 1 satuan, maka charges naik sebesar 343.06 dengan variabel prediktor lainnya tetap
Sexfemale tidak masuk sebagai prediktor, artinya dijadikan basis, koefisien sexmale = -83.57, artinya pria menurunkan charges sebesar 83.57 apabila dibandingan dengan wanita
Smokerno dijadikan basis, koefisien smokeryes = 23788.23, artinya perokok menaikkan charges sebesar 23788.23 apabila dibandingkan bukan perokok
Regionnortheast dijadikan basis, koefisien regionnorthwest = -319.25, artinya regionnorthwest menurunkan charges sebesar 319.25 apabila dibandingkan regionnortheast
Regionnortheast dijadikan basis, koefisien regionsoutheast = -1032.15, artinya regionsoutheast menurunkan charges sebesar 1032.15 apabila dibandingkan regionnortheast
Regionsouthwest dijadikan basis, koefisien regionsouthwest = -938.87 , artinya regionsouthwest menurunkan charges sebesar 938.87 apabila dibandingkan regionnortheast
Adjusted R-squared: 0.7478, artinya model kita berhasil menangkap informasi/varians sebesar 74.78% untuk memprediksi variabel target

8. 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 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 dengan mengabaikan NA
rmse <- sqrt(mean((insurance$pred_charges - insurance$charges)^2, na.rm = TRUE))
print(rmse)

## [1] 6045.715

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:

A. 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

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

## 
##  Pearson's product-moment correlation
## 
## data:  insurance$bmi and insurance$charges
## t = 7.289, df = 1328, p-value = 5.342e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1438988 0.2472785
## sample estimates:
##       cor 
## 0.1961336

Kesimpulan: Nilai p-value = 5.342e-13 < 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.

B. 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(model_all$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_all$residuals
## W = 0.90004, p-value < 2.2e-16

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

C. Autokorelasi

Hipotesis:

H0 : tidak ada autokorelasi

H1 : ada autokorelasi

α : 0.05

Statistik Uji: Durbin Watson dengan nilai P-value

durbinWatsonTest(model_all)

##  lag Autocorrelation D-W Statistic p-value
##    1      0.05395633      1.891228   0.044
##  Alternative hypothesis: rho != 0

💡 Insight: * Nilai autocorrelation pada lag 1 adalah 0.054, yang sangat mendekati 0. Ini menunjukkan bahwa tidak ada autokorelasi yang signifikan pada lag 1. * Nilai D-W 1.891 mendekati 2, yang biasanya menunjukkan bahwa tidak ada autokorelasi residual yang signifikan. Namun, karena nilai ini sedikit kurang dari 2, ada kemungkinan kecil adanya autokorelasi positif. * P-Value 0.058 lebih besar dari 0.05 tetapi sangat dekat dengan level signifikansi 0.05, menunjukkan bahwa meskipun tidak ada cukup bukti untuk menolak hipotesis nol pada level signifikansi 5%, ada indikasi bahwa mungkin ada autokorelasi residual.

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

D. 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 = 120.4, df = 8, p-value < 2.2e-16

💡 Insight: * P-Value < 2.2e-16 menunjukkan bahwa P-Value sangat kecil, jauh di bawah level signifikansi 0.05, maka tolak H0.

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

E. 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.

vif(model_all)

##              GVIF Df GVIF^(1/(2*Df))
## age      1.018435  1        1.009176
## sex      1.008837  1        1.004409
## bmi      1.110440  1        1.053774
## children 1.003980  1        1.001988
## smoker   1.012021  1        1.005992
## region   1.101955  3        1.016313

💡 Insight: * 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.

9. 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.

Turning 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.10000 -0.19667 -0.05002  0.06782  2.14864 
## 
## Coefficients: (1 not defined because of singularities)
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      7.034456   0.072822  96.598  < 2e-16 ***
## age              0.034107   0.000881  38.714  < 2e-16 ***
## sexmale         -0.067589   0.024588  -2.749  0.00606 ** 
## bmi              0.013798   0.002115   6.524 9.73e-11 ***
## children         0.099926   0.010179   9.817  < 2e-16 ***
## smokeryes        1.550733   0.030408  50.998  < 2e-16 ***
## regionnorthwest -0.057166   0.035163  -1.626  0.10425    
## regionsoutheast -0.154701   0.035282  -4.385 1.25e-05 ***
## regionsouthwest -0.123386   0.035249  -3.500  0.00048 ***
## pred_charges           NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.444 on 1307 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.7649, Adjusted R-squared:  0.7634 
## F-statistic: 531.5 on 8 and 1307 DF,  p-value: < 2.2e-16

Didapatkan persamaan regresi sebagai berikut: charges = 7.034456 + 0.034107 age − 0.067589 sexmale + 0.013798 bmi + 0.099926 children + 1.550733 smokeryes − 0.057166 regionnorthwest − 0.154701 regionsoutheast − 0.123386 regionsouthwest

A. 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.10000 -0.19667 -0.05002  0.06782  2.14864 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      7.034456   0.072822  96.598  < 2e-16 ***
## age              0.034107   0.000881  38.714  < 2e-16 ***
## sexmale         -0.067589   0.024588  -2.749  0.00606 ** 
## bmi              0.013798   0.002115   6.524 9.73e-11 ***
## children         0.099926   0.010179   9.817  < 2e-16 ***
## smokeryes        1.550733   0.030408  50.998  < 2e-16 ***
## regionnorthwest -0.057166   0.035163  -1.626  0.10425    
## regionsoutheast -0.154701   0.035282  -4.385 1.25e-05 ***
## regionsouthwest -0.123386   0.035249  -3.500  0.00048 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.444 on 1307 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.7649, Adjusted R-squared:  0.7634 
## F-statistic: 531.5 on 8 and 1307 DF,  p-value: < 2.2e-16

B. Model Forward

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

insurance_clean <- na.omit(insurance)
trans_model_nopredictor <- lm(log(charges) ~ 1, data = insurance_clean)

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 = FALSE)

summary(trans_model_forward)

## 
## Call:
## lm(formula = log(charges) ~ pred_charges + age + children + bmi + 
##     smoker + sex, data = insurance_clean)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.10416 -0.19791 -0.04842  0.06799  2.14218 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.681e+00  3.585e-01  24.215  < 2e-16 ***
## pred_charges  1.382e-04  2.929e-05   4.720 2.62e-06 ***
## age          -7.436e-04  7.464e-03  -0.100 0.920662    
## children      3.278e-02  1.743e-02   1.881 0.060208 .  
## bmi          -3.375e-02  9.753e-03  -3.461 0.000556 ***
## smokeryes    -1.738e+00  6.966e-01  -2.495 0.012705 *  
## sexmale      -5.596e-02  2.469e-02  -2.266 0.023598 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4437 on 1309 degrees of freedom
## Multiple R-squared:  0.7648, Adjusted R-squared:  0.7637 
## F-statistic: 709.4 on 6 and 1309 DF,  p-value: < 2.2e-16

C. Model Both

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

summary(trans_model_both)

## 
## Call:
## lm(formula = log(charges) ~ pred_charges + children + bmi + smoker + 
##     sex, data = insurance_clean)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.10403 -0.19716 -0.04798  0.06668  2.14136 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.646e+00  6.596e-02 131.090   <2e-16 ***
## pred_charges  1.353e-04  3.450e-06  39.229   <2e-16 ***
## children      3.418e-02  1.036e-02   3.298   0.0010 ** 
## bmi          -3.281e-02  2.420e-03 -13.558   <2e-16 ***
## smokeryes    -1.669e+00  8.648e-02 -19.304   <2e-16 ***
## sexmale      -5.620e-02  2.457e-02  -2.287   0.0223 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4436 on 1310 degrees of freedom
## Multiple R-squared:  0.7648, Adjusted R-squared:  0.7639 
## F-statistic: 851.9 on 5 and 1310 DF,  p-value: < 2.2e-16

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 27471.61 27481.97 0.0000000   0.0000000 0.9125434
## 2         trans_model_all 25582.52 25634.34 0.7648744   0.7634353 0.4424900
## 3    trans_model_backward 25582.52 25634.34 0.7648744   0.7634353 0.4424900
## 4     trans_model_forward 25578.96 25620.42 0.7647954   0.7637173 0.4425644
## 5        trans_model_both 25576.97 25613.25 0.7647936   0.7638959 0.4425660

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.

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.99975 -0.43355  0.05411  0.45511  1.74647 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.09298    0.09508  85.115  < 2e-16 ***
## children     0.12034    0.01501   8.016 2.38e-15 ***
## bmi          0.01988    0.00296   6.717 2.76e-11 ***
## smokeryes    1.51655    0.04483  33.830  < 2e-16 ***
## sexmale     -0.09038    0.03625  -2.493   0.0128 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6572 on 1320 degrees of freedom
##   (8 observations deleted due to missingness)
## Multiple R-squared:  0.4881, Adjusted R-squared:  0.4865 
## F-statistic: 314.6 on 4 and 1320 DF,  p-value: < 2.2e-16

Pengujian Asumsi

A. 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.

B. 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.test(step_both2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  step_both2$residuals
## W = 0.99522, p-value = 0.0003274

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

C. Autokorelasi

Hipotesis: H0 : tidak ada autokorelasi H1 : ada autokorelasi α : 0.05

Statistik Uji: Durbin Watson dengan nilai P-value

durbinWatsonTest(step_both2)

##  lag Autocorrelation D-W Statistic p-value
##    1     -0.03494602      2.069749   0.196
##  Alternative hypothesis: rho != 0

💡 Insight: * Aturan Keputusan Nilai p-value = 0.196 > 0.05 sehingga H0 gagal ditolak

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

D. 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 = 215.39, df = 4, p-value < 2.2e-16

💡 Insight: * P-Value < 2.2e-16 menunjukkan bahwa P-Value sangat kecil, jauh di bawah level signifikansi 0.05, maka tolak H0.

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

E. Multikolinearitas

vif(step_both2)

## children      bmi   smoker      sex 
## 1.000463 1.001912 1.005604 1.007732

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

10. 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.09298 + 0.12034 children + 0.01988 bmi + 1.51655 smokeryes - 0.09038 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.12034) = 1.1273 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.01988) = 1.0201 kali dengan variabel lain tetap.

11. 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.