InsurancePricePredictor

Library Imports

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(corrplot)
## corrplot 0.92 loaded
library(caTools) 
library(glmnet) 
## Loading required package: Matrix
## Loaded glmnet 4.1-8
library(grid)
library(ggplot2)
library(broom)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Data Loading and Initial Exploration

# Loading the dataset
insurance_data = read.csv('Medical_insurance.csv')

# A glimse of data
head(insurance_data)
##   age    sex    bmi children smoker    region   charges
## 1  19 female 27.900        0    yes southwest 16884.924
## 2  18   male 33.770        1     no southeast  1725.552
## 3  28   male 33.000        3     no southeast  4449.462
## 4  33   male 22.705        0     no northwest 21984.471
## 5  32   male 28.880        0     no northwest  3866.855
## 6  31 female 25.740        0     no southeast  3756.622

Data Dimensions

dim(insurance_data)
## [1] 2772    7

Basic Data Summary

# Checking basic data distribution
summary(insurance_data)
##       age            sex                 bmi           children    
##  Min.   :18.00   Length:2772        Min.   :15.96   Min.   :0.000  
##  1st Qu.:26.00   Class :character   1st Qu.:26.22   1st Qu.:0.000  
##  Median :39.00   Mode  :character   Median :30.45   Median :1.000  
##  Mean   :39.11                      Mean   :30.70   Mean   :1.102  
##  3rd Qu.:51.00                      3rd Qu.:34.77   3rd Qu.:2.000  
##  Max.   :64.00                      Max.   :53.13   Max.   :5.000  
##     smoker             region             charges     
##  Length:2772        Length:2772        Min.   : 1122  
##  Class :character   Class :character   1st Qu.: 4688  
##  Mode  :character   Mode  :character   Median : 9333  
##                                        Mean   :13261  
##                                        3rd Qu.:16578  
##                                        Max.   :63770

Structure

# Printing the data types of all columns
str(insurance_data)
## 'data.frame':    2772 obs. of  7 variables:
##  $ age     : int  19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : chr  "female" "male" "male" "male" ...
##  $ bmi     : num  27.9 33.8 33 22.7 28.9 ...
##  $ children: int  0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr  "yes" "no" "no" "no" ...
##  $ region  : chr  "southwest" "southeast" "southeast" "northwest" ...
##  $ charges : num  16885 1726 4449 21984 3867 ...
# No need of any data type conversion

Preliminary Exploratory Data Analysis

#EDA
par(mfrow = c(2, 2))

Distribution Analysis of Continuous Variables

hist_colors <- c("#99CCFF", "#99FF99", "#FFCC33", "#FF99CC")
hist(insurance_data$age, main = "Age Distribution", xlab = "Age", col = hist_colors[1])

hist(insurance_data$bmi, main = "BMI Distribution", xlab = "BMI", col = hist_colors[2])

hist(insurance_data$charges, main = "Charges Distribution", xlab = "Charges", col = hist_colors[4])

Categorical Variable Frequency Analysis

# Categorical Variables Analysis

table(insurance_data$sex)
## 
## female   male 
##   1366   1406
table(insurance_data$smoker)
## 
##   no  yes 
## 2208  564
table(insurance_data$region)
## 
## northeast northwest southeast southwest 
##       658       664       766       684

Visual Analysis of Categorical Variables

# Visualization of Categorical Variables

barplot_colors <- c("#99CCFF", "#99FF99", "#FFFFCC")

par(mfrow = c(2, 2))

barplot(table(insurance_data$sex), main = "Sex Distribution", xlab = "Sex", ylab = "Frequency", col = barplot_colors[1])
barplot(table(insurance_data$smoker), main = "Smoker Distribution", xlab = "Smoker", ylab = "Frequency", col = barplot_colors[2])
barplot(table(insurance_data$region), main = "Region Distribution", xlab = "Region", ylab = "Frequency", col = barplot_colors[3])
barplot(table(insurance_data$children), main = "Children Distribution", xlab = "children", ylab = "Frequency")

Bivariate Analysis: BMI vs. Charges

par(mfrow = c(3, 2))
# Scatter plot of BMI vs. Charges
plot(insurance_data$bmi, insurance_data$charges, 
     xlab = "BMI", ylab = "Charges", main = "BMI vs. Charges", col = "#0000FF")

# Scatter plot of Age vs. Charges
plot(insurance_data$bmi, insurance_data$charges, 
     xlab = "age", ylab = "Charges", main = "Age vs. Charges", col = "#FF0000")

# Box plot of Charges by Smoker
boxplot(charges ~ smoker, data = insurance_data, 
        xlab = "Smoker", ylab = "Charges", main = "Charges by Smoker", col = c("#99CCFF", "#99FF99"))

# Box plot of Charges by Region
boxplot(charges ~ region, data = insurance_data, 
        xlab = "Region", ylab = "Charges", main = "Charges by Region", col = c("#99CCFF", "#99FF99", "#FFFFCC", "#FF99CC"))

# Box plot of Charges by Sex
boxplot(charges ~ sex, data = insurance_data, 
        xlab = "Sex", ylab = "Charges", main = "Charges by Sex", col = c("#99CCFF", "#99FF99"))

# Box plot of Charges by Number of Children
boxplot(charges ~ children, data = insurance_data, 
        xlab = "Number of Children", ylab = "Charges", main = "Charges by Number of Children", col = "#99CCFF")

Correlation Analysis Among Numeric Variables

# Creating the correlation matrix
correlation_matrix <- round(cor(insurance_data[, c('charges', 'age', 'bmi', 'children')]), 2)

# Visualizing the correlation matrix plot
corrplot(correlation_matrix, method = "color")

Multivariate Scatter Plot Analysis

# Scatterplot to view relationship between 2 numeric variables
multi_scatter <- ggplot(insurance_data, aes(x = bmi, y = charges, col = smoker)) + 
  geom_point() + 
  ggtitle('Charges Vs BMI (along with smoker division)')

# For 2 plots
grid.arrange( multi_scatter, nrow = 1, ncol = 1, top="Multivariate Analysis")

Data Cleaning

# Checking missing values
# Counting null values in each column
colSums(is.na(insurance_data))
##      age      sex      bmi children   smoker   region  charges 
##        0        0        0        0        0        0        0

Data Preparation

# Converting categorical columns to numerical columns
# Label encoding for 'sex' column
insurance_data$sex <- as.numeric(factor(insurance_data$sex, levels = c("female", "male")))

# Label encoding for 'region' column
insurance_data$region <- as.numeric(factor(insurance_data$region, levels = unique(insurance_data$region)))

# Label encoding for 'smoker' column
insurance_data$smoker <- as.numeric(factor(insurance_data$smoker, levels = c("no", "yes")))

Data Standardization

# Standardisation for numerical columns
insurance_data$age <- scale(insurance_data$age)
insurance_data$bmi <- scale(insurance_data$bmi)

Data Splitting for Model Training

# Spliting the data into training and testing sets (80% train, 20% test)
split <- sample.split(insurance_data$charges, SplitRatio = 0.8)
train_data <- subset(insurance_data, split == TRUE)
test_data <- subset(insurance_data, split == FALSE)

Model Selection

# Loading necessary libraries
library(caret)
## Loading required package: lattice
# Defining training control for cross-validation
train_control <- trainControl(method = "cv", number = 10)  # 10-fold cross-validation

# Training the base model using all predictors
lm_base <- train(charges ~ ., data = train_data, method = "lm", trControl = train_control)

# Training a reduced model using selected predictors based on domain knowledge
lm_reduced <- train(charges ~ age + bmi + smoker, data = train_data, method = "lm", trControl = train_control)

# Comparing models using resampling
results <- resamples(list(base = lm_base, reduced = lm_reduced))

# Summary of the results
print(summary(results))
## 
## Call:
## summary.resamples(object = results)
## 
## Models: base, reduced 
## Number of resamples: 10 
## 
## MAE 
##             Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## base    3945.955 4137.483 4161.493 4177.120 4251.962 4353.188    0
## reduced 3798.133 4130.982 4210.757 4226.513 4320.720 4557.378    0
## 
## RMSE 
##             Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## base    5518.170 5890.773 6119.195 6065.842 6266.604 6387.128    0
## reduced 5397.041 5854.595 6057.641 6101.341 6163.082 6863.696    0
## 
## Rsquared 
##              Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
## base    0.7010563 0.7331599 0.7393472 0.7503225 0.7787187 0.8142940    0
## reduced 0.6745369 0.7409773 0.7548370 0.7475536 0.7566856 0.7994139    0
# Box and whisker plot to visualize results
bwplot(results)

* The base model outperforms the reduced model across all three primary metrics (MAE, RMSE, and Rsquared), suggesting it is the better model for predicting charges in this case.

Linear Regression Model Fitting

# Fitting a linear regression model to predict charges
lm_model <- lm(charges ~ ., data = train_data)

# Printing summary of the linear regression model
summary(lm_model)
## 
## Call:
## lm(formula = charges ~ ., data = train_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -11474  -2785  -1013   1375  29691 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -16821.29     563.65 -29.844  < 2e-16 ***
## age           3613.42     117.31  30.804  < 2e-16 ***
## sex            -92.88     233.51  -0.398  0.69084    
## bmi           2008.77     118.89  16.896  < 2e-16 ***
## children       495.55      96.15   5.154 2.73e-07 ***
## smoker       23883.45     289.29  82.559  < 2e-16 ***
## region         371.05     106.63   3.480  0.00051 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6064 on 2716 degrees of freedom
## Multiple R-squared:  0.7506, Adjusted R-squared:  0.7501 
## F-statistic:  1363 on 6 and 2716 DF,  p-value: < 2.2e-16

Analysis of Residual Distributions

# Distribution of residuals
resid_hist <- ggplot(lm_model, aes(x=lm_model$residuals)) + geom_histogram()

# Variance of residuals
resid_fitted <- ggplot(lm_model, aes(x=lm_model$fitted.values, y=lm_model$residuals)) + 
  geom_point() + 
  ggtitle('Residuals vs. Fitted Values')

# For 2 plots
grid.arrange(resid_hist, resid_fitted, nrow=2, ncol=1, top="Regression Assumptions")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Making Predictions on Test Data

# Predicting charges on the test set
predictions <- predict(lm_model, newdata = test_data)

Model Evaluation: RMSE Calculation

# Calculating RMSE
rmse <- sqrt(mean((test_data$charges - predictions)^2))
print(paste("Root Mean Squared Error (RMSE):", round(rmse, 2)))
## [1] "Root Mean Squared Error (RMSE): 6575.4"

Model Evaluation: MAPE Calculation

# Calculating MAPE
mape <- mean(abs((test_data$charges - predictions) / test_data$charges)) * 100
print(paste("Mean Absolute Percentage Error (MAPE):", round(mape, 2)))
## [1] "Mean Absolute Percentage Error (MAPE): 54.55"