Introduction

The goal of this R markdown is to predict insurance costs using medical cost personal datasets provided by https://github.com/stedy/Machine-Learning-with-R-datasets/blob/master/insurance.csv.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(AnalysisLin)

Data Preprocessing

data<-read.csv("/Users/zhiweilin/Downloads/insurance.csv", header =T, na.string=c("","NA"))
desc_stat(data) # desc_stat() is a function from my personal EDA package AnalysisLin. The function is used to generate some key statistics for the dataset. 
## Descriptive Statistics Results:
## =========================================================================
##               age     sex       bmi      children smoker    region   
## count         1338    1338      1338     1338     1338      1338     
## unique        47      2         548      6        2         4        
## duplicate     1       1         1        1        1         1        
## null          0       0         0        0        0         0        
## null_rate     0       0         0        0        0         0        
## type          integer character numeric  integer  character character
## min           18      <NA>      15.96    0        <NA>      <NA>     
## p25           27      <NA>      26.29625 0        <NA>      <NA>     
## mean          39.207  <NA>      30.6634  1.0949   <NA>      <NA>     
## median        39      <NA>      30.4     1        <NA>      <NA>     
## p75           51      <NA>      34.69375 2        <NA>      <NA>     
## max           64      <NA>      53.13    5        <NA>      <NA>     
## sd            14.05   <NA>      6.0982   1.2055   <NA>      <NA>     
## kurtosis      1.7551  <NA>      2.945    3.1972   <NA>      <NA>     
## skewness      0.0556  <NA>      0.2837   0.9373   <NA>      <NA>     
## jarque_pvalue 0       <NA>      0        0        <NA>      <NA>     
##               charges     
## count         1338        
## unique        1337        
## duplicate     1           
## null          0           
## null_rate     0           
## type          numeric     
## min           1121.8739   
## p25           4740.28715  
## mean          13270.4223  
## median        9382.033    
## p75           16639.912515
## max           63770.42801 
## sd            12110.0112  
## kurtosis      4.5958      
## skewness      1.5142      
## jarque_pvalue 0

Some data Insights:

data<-distinct(data) # Remove duplicate rows based on all columns
data <- mutate_at(data, vars(sex,smoker,region), as.factor)

The variables ‘sex’, ‘smoker’, and ‘region’ are converted from character variables to factor variables.

Data Visualization

numeric_plot(data,prob=T,dens=T) # numeric_plot() is a function from my personal EDA package AnalysisLin. The function is used to generate histogram/line plot for all numerical variables. 

## list()

Some takeaways:

categoric_plot(data,n_col = 2,bar_legend = F) # categorical_plot() is a function from my personal EDA package AnalysisLin. The function is used to generate pie/bar plot for all categrical(level) variables.

Some takeaways:

Correlation Analysis

numeric_columns <- sapply(data, is.numeric)
numeric_data <- data[, numeric_columns]
corr_matrix(numeric_data,corrplot=T) #corr_matrix() is a function from my personal EDA package, AnalysisLin. The function is used to generated a table of correlation between variables and a plot of correlations.
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: corrplot
## corrplot 0.92 loaded

##        row   column        cor            p
## 4      age  charges 0.29830821 0.000000e+00
## 5      bmi  charges 0.19840083 2.469136e-13
## 1      age      bmi 0.10934361 6.164372e-05
## 6 children  charges 0.06738935 1.371703e-02
## 2      age children 0.04153621 1.290128e-01
## 3      bmi children 0.01275466 6.412463e-01

For numerical variables only, age is most importnat variable in determining insurance permium, followed by bmi.

corr_cluster(numeric_data)

This is correlation clustering that tells the relationships between four numerical variables. As can be seen, variables age is closest to variable charges.

Model Selection and Training

Random Forest Model

library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:ggplot2':
## 
##     margin
set.seed(123)
training.samples <- data$charges %>%
  createDataPartition(p = 0.8, list = FALSE)
train.data  <- data[training.samples, ]
test.data <- data[-training.samples, ]

split dataset into train and test data, with 80% train data and 20% test data.

set.seed(123)
rf_model <- train(
  charges ~., data = train.data, method = "rf",
  trControl = trainControl("cv", number = 10),
  importance = TRUE
  )

10-fold cross-validation is used

rf_model$bestTune #best set of tuning parameter.
##   mtry
## 2    5

Extreme Gradient Boosting Model

library(xgboost)
## 
## Attaching package: 'xgboost'
## The following object is masked from 'package:dplyr':
## 
##     slice
xgb_model <- train(
  charges ~ ., data = train.data, method = "xgbTree",
  trControl = trainControl(method = "cv", number = 10),
  verbosity = 0
)
xgb_model$bestTune #best set of tuning parameter.
##    nrounds max_depth eta gamma colsample_bytree min_child_weight subsample
## 34      50         2 0.3     0              0.8                1         1

Model Evaluation

Random Forest Result

varImp(rf_model)
## rf variable importance
## 
##                  Overall
## smokeryes       100.0000
## bmi              51.1292
## age              50.7188
## children          4.9639
## regionsoutheast   2.0844
## regionsouthwest   0.3861
## regionnorthwest   0.3150
## sexmale           0.0000

In random forest model, smokeryes is the most significant variable in determining insurance charges, followed by BMI and age. Other factors have minor to no impact on the charges

varImpPlot(rf_model$finalModel, type = 1)

rf_predictions <- rf_model %>% predict(test.data)
head(rf_predictions)
##         5         9        14        19        25        28 
##  4877.633  7903.622 12120.191 13150.115  6576.011 12628.450
RMSE(rf_predictions, test.data$charges)
## [1] 4806.051

Extreme Gradient Boosting Result

varImp(xgb_model)
## xgbTree variable importance
## 
##                   Overall
## smokeryes       100.00000
## bmi              22.97129
## age              14.01767
## children          1.11484
## regionsouthwest   0.14746
## regionsoutheast   0.10335
## sexmale           0.02424
## regionnorthwest   0.00000

In XGBoost model, smokeryes is the most significant variable in determining insurance charges, followed by BMI and age. Other factors have minor to no impact on the charges

xgb_predictions <- xgb_model %>% predict(test.data)
head(xgb_predictions)
## [1]  5323.405  7166.240 12638.623 10981.429  7134.224 12472.825
RMSE(xgb_predictions, test.data$charges)
## [1] 4787.083

Conclusion

Both the Random Forest and XGBoost models yield the same result, such that “smoker,” “bmi,” and “age” are three crucial variables that determines insurance premiums. However, the XGBoost model outperforms the Random Forest model with a lower Root Mean Squared Error (RMSE). In general, if your dataset is small and you value simplicity, the Random Forest model is a solid choice. Conversely, for larger datasets where predictive accuracy is your goal, the XGBoost model is a better fit.