Insurance.Rmd

---
title: "Insurance Healthcare Analysis"
author: "Shardha Ramlal"
date: "2026-05-17"
output:
  word_document: default
  pdf_document: default
editor_options:
  markdown:
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#load packages : https://epirhandbook.com/en/new_pages/cleaning.html
if(!require(pacman)) install.packages("pacman", repos = "http://cran.us.r-project.org")

pacman::p_load(officer, officedown,tidyverse,
               dplyr, RKaggle,skimr, gtsummary, flextable, # tables design
                tidymodels, corrplot, dlookr, tidymodels, workflowsets, finetune,                 doParallel,corrplot )

data <- RKaggle::get_dataset("mirichoi0218/insurance")

# head(data)

# dim(data)
# 1338    7

#explore the dataset

str(data)

data$sex <- as.factor(data$sex)
data$smoker <- as.factor(data$smoker)
data$region <- as.factor(data$region)

summary(data)

# skim(data)

# write dataset to excel
# write.csv(spotify,"spotify.csv")

# we would only be looking at the spotify data here, so will subset.

# drop the missing data from the 
data <- data[complete.cases(data),]

#duplicate cases
data <- data[!duplicated(data), ] # 1duplictae

```

# Table of Contents

```{r, echo=FALSE, warning=FALSE}
block_toc()
```

## Introduction

The research questions:

1.To predict the insurance cost of patients admitted to hospitals

2.To determine the main factors affecting insurance costs

## Methodology

1.  The dataset for analysis was downloaded from Kaggle via URL1:
    <https://www.kaggle.com/datasets/mirichoi0218/insurance>

data is also available via github2
<https://github.com/stedy/Machine-Learning-with-R-datasets/blob/master/insurance.csv>

inspiration -
<https://www.kaggle.com/code/ruslankl/health-care-cost-prediction-w-linear-regression/input3>

2.  Duplicate rows were removed (one present) 4,5

3.  Missing predictor variable (charges) was removed

4.  Exploratory data analysis (EDA) was constructed.

5.  Duplicates in the features were removed

6.  Dataset was split into a 75-25 training-test set combination.
    Charges was log transformed since data was right skewed

7.  Statistical and Machine Learning predictive models were built using
    the tidymodels R framework5. The models run were:

    Linear Regression

    Lasso Regression

    Decision Trees

    Random Forest

    Gradient Boosting

    Neural Network (MLP)

8.  RMSE was used as the comparison metric
    $RMSE = \sqrt{\frac{1}{n} \sum_{i}^{} \left( \hat{y}_{i} - y_{i} \right)^2 }$

# Results

## Data Description and Visualisation

age: age of primary beneficiary

sex: insurance contractor gender, female, male

bmi: Body mass index, providing an understanding of body, weights that
are relatively high or low relative to height, objective index of body
weight (kg / m \^ 2) using the ratio of height to weight, ideally 18.5
to 24.9

children: Number of children covered by health insurance / Number of
dependents

smoker: Smoking

region: the beneficiary's residential area in the US, northeast,
southeast, southwest, northwest.

charges: Individual medical costs billed by health insurance

## Exploratory Data Analysis

```{r, echo=FALSE, warning=FALSE}
# Export outputs and examine html files
# data |>
#    eda_web_report(target = "charges", subtitle = "Insurance Charges",
#                   output_dir = getwd(),
#                   output_file = "EDA_Insurance.html", theme = "blue",
#                   browse=TRUE)


 # data |>
 #   diagnose_web_report(target = "charges", subtitle = "Insurance Charges",
 #                 output_dir = getwd(),
 #                  output_file = "Diagnose_Insurance.html", theme = "blue",
 #                  browse=TRUE)

#scatterplot
# pairs(~., data = data)

```

## Histogram of Insurance Charges

```{r, echo=FALSE, warning=FALSE, fig.width = 7, fig.height = 6}
#drop missing rows
data <- drop_na(data)

#log transform charges
data <- data |> mutate(charges_log = log(charges))

par(mfrow = c(1, 2))

hist(data$charges, prob = TRUE, col = "lightblue", 
     main = "Insurance Charges", xlab = "Charges", cex.main = 0.8)
lines(density(data$charges), col = "red", lwd = 1)

hist(data$charges_log, prob = TRUE, col = "lightblue", 
     main = "Log Transformed Insurance Charges", xlab = "Log(Charges)",cex.main = 0.8)
lines(density(data$charges_log), col = "red", lwd = 1)

par(mfrow = c(1, 1))




```

From the graph, the charges variable has a right skewed distribution -
solution is to log transform, the charges variable

## Summary Statistics of dataset

```{r, echo=FALSE, warning=FALSE}

data |>
  tbl_summary(
    include=c(age, sex, bmi, children, smoker, region, charges, charges_log),
    statistic = list(all_continuous() ~ "{mean} ({sd})" ,
                     all_categorical() ~ "{n} ({p}%)") ,
    digits = list(all_continuous() ~ 2,
                  all_categorical() ~ c(0, 1)), 
      label = list(age ~ "Age",sex ~ "Sex",bmi~ "BMI", smoker~ "Smoker",
    region ~ "Region", charges ~ "Charges" , charges_log ~ "Charges (Log)"),
    missing = "no") |> bold_labels() |> 
  gtsummary::as_flex_table() |>
  flextable::font(fontname = "Times New Roman", part = "all") |>
  flextable::fontsize(size = 12, part = "all") |>
  flextable::autofit()
   


```

### Partition and Split into Training and Test

75 % of the data was used for training, while the remaining was used for the test set. 

```{r partition, echo=FALSE }
# Split dataset into training and test dataset using the tidymodels interface
set.seed(1918)

data_subset <- data |> select(-charges_log)

data_split <- initial_split(data_subset, prop = 0.75)     

# Create data frames for the two sets:
train <- training(data_split)
test  <- testing(data_split)

```

### Dimensionality Reduction

Correlation Plot

```{r, echo=FALSE, warning = FALSE}
library()

tmwr_cols <- colorRampPalette(c("#91CBD765", "#CA225E"))
train |>
  dplyr::select(-c( charges)) |>
   dplyr::select(where(is.numeric)) |>  
    cor() |>
  corrplot(col = tmwr_cols(200), tl.col = "black", method = "ellipse")



train |>
  dplyr::select(-c(charges)) |>
  dplyr::select(where(is.numeric)) |>  
  cor(use = "complete.obs") |>
  as.data.frame() |>
  tibble::rownames_to_column("Variable") |> 
  flextable() |>
  colformat_double(digits = 3) |>
  set_caption("Correlation Matrix") |>
  autofit()


```

## Models

### Build Models using Tidymodels Workflow

The dependent variable (insurance cost) was log transformed for the training dataset, dummy variables were created for categorical variables (smoker, region, sex), numeric predictors were standardised by subtracting the mean and dividing by the standard deviation to ensure the same scale was used for analysis. Predictors with zero variance were removed. 

To ensure speed in computing, models were run in parallel. 

10 fold cross validation was used for tuning parameters in each model. 

```{r, warning=FALSE, echo=FALSE, message=FALSE}


cl <- makePSOCKcluster(parallel::detectCores(logical = FALSE) - 1)
registerDoParallel(cl)


set.seed(1918)

rec <- recipe(charges ~ ., data = train) |>
  step_log(charges, offset = 1) |> # Log transform the outcome
  step_dummy(all_nominal_predictors()) |>
  step_novel(all_nominal_predictors()) |>
  step_zv(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())



#Evaluation metrics set & 10 fold cross validation
metrics   <- metric_set(rmse)

#10 fold cross validation to train tuning parameters 
cv_folds <- vfold_cv(train, v = 10)

lm_model <- linear_reg() |> 
  set_engine("lm") |> 
  set_mode("regression")

lasso_model <- linear_reg(penalty = tune(), mixture = 1) |> 
  set_engine("glmnet") |> 
  set_mode("regression")

dt_model  <-  decision_tree(cost_complexity = tune(), tree_depth = tune(), min_n = tune()) |> 
  set_engine("rpart") |> 
  set_mode("regression")

rf_model  <-  rand_forest(mtry = tune(), trees = 500, min_n = tune()) |> 
  set_engine("ranger") |> 
  set_mode("regression")
 
gbm_model <- boost_tree(trees = tune(), tree_depth = tune(), learn_rate = tune(), loss_reduction = tune())|>
  set_engine("xgboost") |>
  set_mode("regression")

nn_model <-  mlp(hidden_units = tune(), penalty = tune(), epochs = tune()) |>      set_engine("nnet" , MaxNWts = 2000) |> 
  set_mode("regression")

nn_param <- extract_parameter_set_dials(nn_model) |> 
  update(hidden_units = hidden_units(c(1, 27))) #M. Kuhn and Johnson (2013)


all_workflows <- workflow_set(preproc = list(base = rec),
                              models = list(
                                lm = lm_model, 
                                lasso = lasso_model,
                                dt = dt_model, 
                                rf = rf_model,
                                gbm = gbm_model,
                                nn=nn_model))





race_ctrl <-control_race(save_pred = TRUE,parallel_over = "everything",save_workflow = TRUE)

race_results <-all_workflows |>workflow_map("tune_race_anova",
      seed = 1918,
      resamples = cv_folds,
      grid = 5,
      control = race_ctrl,
      metrics= metrics
   )


stopCluster(cl)


race_results |>
  rank_results() |>
  filter(.metric == "rmse") |> 
  select(model, .config, rmse = mean, rank) |> 
  knitr::kable(caption = "Model Ranking by RMSE (Training data)")


#best model is the neural network


```

The best model was seen to be the neural network

# Autoplot of Best Model (Smallest RSME with Training data)

```{r, warning=FALSE, echo=FALSE}

library(scales)

autoplot(race_results,  rank_metric = "rmse", metric = "rmse") +
  theme_minimal() +
  labs(title = "Model Comparison of RMSE across resamples")

 
```

The best model is the neural network. For research question 1, the neural network models appears to fit this dataset the best. 

## Best Results using the test data

```{r, warning=FALSE, echo=FALSE, messages=FALSE}

best_results <- 
   race_results |>
   extract_workflow_set_result("base_nn") |>
   select_best(metric = "rmse") 

knitr::kable(best_results, caption = "Best Neural Network Parameters")

nn_test_results <- 
   race_results |>
   extract_workflow("base_nn") |> 
   finalize_workflow(best_results) |>
   last_fit(split = data_split) #split object - final on training set then test set

test_metrics <- collect_metrics(nn_test_results)
knitr::kable(test_metrics, caption = "Final Test Set Performance")


#top 5 predictors
library(vip)

nn_test_results |>
  extract_fit_parsnip() |>
  vi() |>
  select(Variable, Importance) |>
  knitr::kable(
    digits = 1,
    caption = "Variable Importance for the Neural Network"
  )


```

The RSME yields 42.3 % using the neural network model on the test set. This may indicate, there are predictors missing / unaccounted for that may explain the variation better to predict insurance costs. For example, the dataset did not include diagnosis and clinical information which may have been possible crucial predictors of insurance costs

After extracting the variable importance for each of the predictors, it appeared that the top predictors of insurance cost were smoking, bmi and children. (research question 2)

## Graph showing how the model fits the test data

```{r, , warning=FALSE, echo=FALSE}
nn_test_results |>
   collect_predictions() |> 
   ggplot(aes(x = charges, y = .pred)) + 
   geom_abline(color = "gray50", lty = 2) + 
   geom_point(alpha = 0.5) + 
   coord_obs_pred() +  
   labs(x = "observed", y = "predicted")


```

The predicted and observed insurance costs align mainly on the diagonal line. 

## References

1.  Medical Cost Personal Datasets [Internet]. [cited 2026 May 14].
    Available from:
    <https://www.kaggle.com/datasets/mirichoi0218/insurance>
2.  stedy. GitHub [Internet]. [cited 2026 May 14].
    Machine-Learning-with-R-datasets/insurance.csv at master ·
    stedy/Machine-Learning-with-R-datasets. Available from:
    <https://github.com/stedy/Machine-Learning-with-R-datasets/blob/master/insurance.csv>
3.  Health Care Cost Prediction w/ Linear Regression [Internet]. [cited
    2026 May 14]. Available from:
    <https://kaggle.com/code/ruslankl/health-care-cost-prediction-w-linear-regression>
4.  Identify and Remove Duplicate Data in R - Datanovia [Internet].
    [cited 2026 Apr 25]. Available from:
    <https://www.datanovia.com/en/lessons/identify-and-remove-duplicate-data-in-r/>
5.  How to Remove Duplicate Rows in R DataFrame? - GeeksforGeeks
    [Internet]. [cited 2026 Apr 25]. Available from:
    <https://www.geeksforgeeks.org/r-language/how-to-remove-duplicate-rows-in-r-dataframe/>
6.  R CHARTS \| A collection of charts and graphs made with the R
    programming language [Internet]. 2021 [cited 2026 Jan 29]. Pie chart
    with percentages in ggplot2. Available from:
    <https://r-charts.com/part-whole/pie-chart-percentages-ggplot2/>
7.  Gohel D. Chapter 5 officedown for Word \| officeverse [Internet].
    [cited 2026 Apr 25]. Available from:
    <https://ardata-fr.github.io/officeverse/officedown-for-word.html>
8.  Science B. Make Microsoft Word Reports with R + officedown \|
    R-bloggers [Internet]. 2023 Mar 3 [cited 2026 Apr 25]. Available
    from:
    <https://www.r-bloggers.com/2023/03/make-microsoft-word-reports-with-r-officedown/>
9.  Riederer YX Christophe Dervieux, Emily. R Markdown Cookbook
    [Internet]. [cited 2026 Apr 25]. Available from:
    <https://yihui.org/rmarkdown-cookbook/>
10. 11 Colour scales and legends – ggplot2: Elegant Graphics for Data
    Analysis (3e) [Internet]. [cited 2026 Apr 25]. Available from:
    <https://ggplot2-book.org/scales-colour#brewer-scales>
11. R CHARTS \| A collection of charts and graphs made with the R
    programming language [Internet]. 2021 [cited 2026 Apr 25]. Pie chart
    in ggplot2. Available from:
    <https://r-charts.com/part-whole/pie-chart-ggplot2/>
12. Sequential, diverging and qualitative colour scales from ColorBrewer
    — scale_colour_brewer [Internet]. [cited 2026 Apr 25]. Available
    from: <https://ggplot2.tidyverse.org/reference/scale_brewer.html>
13. RPubs - Reordering a ggplot bar chart axis [Internet]. [cited 2026
    Apr 25]. Available from:
    <https://www.rpubs.com/dvdunne/reorder_ggplot_barchart_axis>
14. Heiss A. Andrew Heiss [Internet]. [cited 2026 Apr 25]. Quick and
    easy ways to deal with long labels in ggplot2. Available from:
    <https://www.andrewheiss.com/blog/2022/06/23/long-labels-ggplot/>
15. DataCamp [Internet]. [cited 2026 Apr 25]. Sorting Data in R.
    Available from: <https://www.datacamp.com/doc/r/sorting>
16. ggplot2: Elegant Graphics for Data Analysis (3e) [Internet]. [cited
    2026 Apr 25]. Available from: <https://ggplot2-book.org/>
17. GeeksforGeeks [Internet]. 15:59:51+00:00 [cited 2026 Apr 25]. Remove
    Rows with Missing Values using R. Available from:
    <https://www.geeksforgeeks.org/r-language/remove-rows-with-missing-values-using-r/>
18. Presentation-Ready Data Summary and Analytic Result Tables
    [Internet]. [cited 2026 Apr 25]. Available from:
    <https://www.danieldsjoberg.com/gtsummary/>
19. Summary table — tbl_summary [Internet]. [cited 2026 Apr 25].
    Available from:
    <https://www.danieldsjoberg.com/gtsummary/reference/tbl_summary.html#ref-usage>
20. Working with model coefficients – tidymodels [Internet]. [cited 2026
    Apr 26]. Available from:
    <https://www.tidymodels.org/learn/models/coefficients/>
21. Fitting and predicting with parsnip [Internet]. [cited 2026 Apr 26].
    Available from:
    <https://parsnip.tidymodels.org/articles/Examples.html#boost_tree-models>
22. RPubs - Train and Evaluate Regression Models using Tidymodels
    [Internet]. [cited 2026 Apr 26]. Available from:
    <https://rpubs.com/eR_ic/regression>
23. Silge MK and J. 15 Screening Many Models \| Tidy Modeling with R
    [Internet]. [cited 2026 Apr 26]. Available from:
    <https://www.tmwr.org/workflow-sets>
24. Sinks LE. Louise E. Sinks - A Tidymodels Tutorial: A Structural
    Approach [Internet]. 2023 [cited 2026 May 9]. Available from:
    <https://lsinks.github.io/posts/2023-04-10-tidymodels/tidymodels_tutorial.html>
25. A Gentle Introduction to tidymodels · R Views [Internet]. [cited
    2026 May 9]. Available from:
    <https://rviews.rstudio.com/2019/06/19/a-gentle-intro-to-tidymodels/>
    
26.	Obtain and format results produced by tuning functions — collect_predictions [Internet]. [cited 2026 May 17]. Available from: <https://tune.tidymodels.org/reference/collect_predictions.html>

27.	beckmw. Variable importance in neural networks | R-bloggers [Internet]. 2013 Aug 12 [cited 2026 May 17]. Available from: <https://www.r-bloggers.com/2013/08/variable-importance-in-neural-networks/>