title: “Assignment Part 2 - Conor Flannery”

output: pdf_document

date: “04-11-2024”

## Loading various packages required for my analysis. 
library(tidyverse)

## Warning: package 'ggplot2' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── 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(GGally) 

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

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

library(dplyr)
library(MASS)

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select

library(tidyr)
library(skimr)
library(summarytools)

## Warning: package 'summarytools' was built under R version 4.3.3

## 
## Attaching package: 'summarytools'
## 
## The following object is masked from 'package:tibble':
## 
##     view

library(broom)

## Warning: package 'broom' was built under R version 4.3.3

library(car)

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

## Loading required package: carData

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

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

library(parameters)

## Warning: package 'parameters' was built under R version 4.3.3

library(plotly)

## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:MASS':
## 
##     select
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

library(ggthemes)

Medical_costs <- read_csv("C:/Users/CFlannery/Downloads/Medical_cost.csv")

## Rows: 1338 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): sex, smoker, region
## dbl (4): age, bmi, children, charges
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Medical_costs <- distinct(Medical_costs)
head(Medical_costs)

## # A tibble: 6 × 7
##     age sex      bmi children smoker region    charges
##   <dbl> <chr>  <dbl>    <dbl> <chr>  <chr>       <dbl>
## 1    19 female  27.9        0 yes    southwest  16885.
## 2    18 male    33.8        1 no     southeast   1726.
## 3    28 male    33          3 no     southeast   4449.
## 4    33 male    22.7        0 no     northwest  21984.
## 5    32 male    28.9        0 no     northwest   3867.
## 6    31 female  25.7        0 no     southeast   3757.

# Displays the first 6 rows of the dataframe Medical_costs. This is useful for getting an initial look at the structure and content of the data, such as column names and a quick glance at some entries.
summary(Medical_costs) # Provides a statistical summary of the numerical variables in Medical_costs. For numeric columns, it will show min, max, mean, median, and quartiles. For categorical variables, it shows the frequency of different levels. This is useful for identifying missing values and understanding the data distribution.

##       age            sex                 bmi           children    
##  Min.   :18.00   Length:1337        Min.   :15.96   Min.   :0.000  
##  1st Qu.:27.00   Class :character   1st Qu.:26.29   1st Qu.:0.000  
##  Median :39.00   Mode  :character   Median :30.40   Median :1.000  
##  Mean   :39.22                      Mean   :30.66   Mean   :1.096  
##  3rd Qu.:51.00                      3rd Qu.:34.70   3rd Qu.:2.000  
##  Max.   :64.00                      Max.   :53.13   Max.   :5.000  
##     smoker             region             charges     
##  Length:1337        Length:1337        Min.   : 1122  
##  Class :character   Class :character   1st Qu.: 4746  
##  Mode  :character   Mode  :character   Median : 9386  
##                                        Mean   :13279  
##                                        3rd Qu.:16658  
##                                        Max.   :63770

dim(Medical_costs) # This command returns the dimensions of the dataframe Medical_costs in the format (number of rows, number of columns). This gives you a quick understanding of the size of the dataset.

## [1] 1337    7

str(Medical_costs) #This displays the structure of the dataframe. It shows each variable’s type (numeric, character, factor, etc.), and a preview of the first few entries. It's useful to ensure data types are correct.

## tibble [1,337 × 7] (S3: tbl_df/tbl/data.frame)
##  $ age     : num [1:1337] 19 18 28 33 32 31 46 37 37 60 ...
##  $ sex     : chr [1:1337] "female" "male" "male" "male" ...
##  $ bmi     : num [1:1337] 27.9 33.8 33 22.7 28.9 ...
##  $ children: num [1:1337] 0 1 3 0 0 0 1 3 2 0 ...
##  $ smoker  : chr [1:1337] "yes" "no" "no" "no" ...
##  $ region  : chr [1:1337] "southwest" "southeast" "southeast" "northwest" ...
##  $ charges : num [1:1337] 16885 1726 4449 21984 3867 ...

skim(Medical_costs) # provides a more detailed summary than summary(). It gives insights into missing values, data types, and key statistics for each column. This is helpful for a quick data exploration, especially for large datasets.

Data summary

Name	Medical_costs
Number of rows	1337
Number of columns	7
_______________________
Column type frequency:
character	3
numeric	4
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
sex	0	1	4	6	0	2	0
smoker	0	1	2	3	0	2	0
region	0	1	9	9	0	4	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
age	0	1	39.22	14.04	18.00	27.00	39.00	51.00	64.00	▇▅▅▆▆
bmi	0	1	30.66	6.10	15.96	26.29	30.40	34.70	53.13	▂▇▇▂▁
children	0	1	1.10	1.21	0.00	0.00	1.00	2.00	5.00	▇▂▂▁▁
charges	0	1	13279.12	12110.36	1121.87	4746.34	9386.16	16657.72	63770.43	▇▂▁▁▁

View(Medical_costs) # allows you to browse through the data in a more visual, spreadsheet-like format, which is convenient for manual inspection.

Manipulating the data to allow for easier interogation of the dataset.

sum(is.na(Medical_costs)) # Checking for missing values, which is none.

## [1] 0

colSums(is.na(Medical_costs)) # checking for missing values again but in each column.

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

Medical_costs[!complete.cases(Medical_costs), ]

## # A tibble: 0 × 7
## # ℹ 7 variables: age <dbl>, sex <chr>, bmi <dbl>, children <dbl>, smoker <chr>,
## #   region <chr>, charges <dbl>

Medical_costs <- Medical_costs %>%
  mutate(obs = case_when(
    bmi < 18.5 ~ "underweight",
    bmi >= 18.5 & bmi <= 24.999 ~ "normal",
    bmi >= 25 & bmi <= 29.999 ~ "overweight",
    bmi >= 30 ~ "obese"
  ))


Medical_cost <- Medical_costs %>%
  mutate(smoking_status = if_else(smoker == "yes", "1", "0"))

head(Medical_cost)

## # A tibble: 6 × 9
##     age sex      bmi children smoker region    charges obs        smoking_status
##   <dbl> <chr>  <dbl>    <dbl> <chr>  <chr>       <dbl> <chr>      <chr>         
## 1    19 female  27.9        0 yes    southwest  16885. overweight 1             
## 2    18 male    33.8        1 no     southeast   1726. obese      0             
## 3    28 male    33          3 no     southeast   4449. obese      0             
## 4    33 male    22.7        0 no     northwest  21984. normal     0             
## 5    32 male    28.9        0 no     northwest   3867. overweight 0             
## 6    31 female  25.7        0 no     southeast   3757. overweight 0

Boxplots allowing for an understanding of the dataset.

boxplot(charges ~ region, 
        data = Medical_costs, 
        col = c("skyblue","coral","pink","honeydew"),
        ylab = "Charges",
        main = "Price by Region",
        xlab = "Region")

# small quartile range for southwest and northwest customers

boxplot(charges ~ smoker, 
        data = Medical_costs, 
        col = c("green","purple"),
        ylab = "price",
        main = "Prices by smoking",
        xlab = "Smoker")# clearly can see higher charges for smokers.

boxplot(charges ~ children, 
        data = Medical_costs, 
        col = c("cyan","cornflowerblue","aliceblue","darksalmon","orchid","forestgreen"),
        ylab = "price",
        main = "Prices by no. children",
        xlab = "No. of children") # Can observe higher range and value of charges for people with 2 to 3 children.

ggplot(Medical_costs,aes(y=charges,x=obs,fill=obs))+
  geom_boxplot()+
  ggtitle("Boxplot of the charges by the category of BMI")+
  ylab("charges")+
  xlab("category of BMI")+
  theme_solarized()+
  scale_fill_brewer(palette="Set1")+
  theme(legend.position = "none")+
  theme(plot.title = element_text(hjust = 0.5)) #clearly can see wider range and higher charges for obese customers

fig <- ggplot(Medical_costs, aes(x = smoker, y = charges, fill = smoker)) +
  geom_boxplot() + 
  ylab("Cost") +
  xlab("Smoker") +
  theme_bw() +
  theme(legend.position = "none")
ggplotly(fig) # Smokers are incurring higher costs than non-smokers

Medical_costs %>%
  ggplot(aes(x = charges, fill = factor(smoker))) +
  geom_histogram(position = "identity", alpha = 0.7, bins = 30) +
  labs(x = "Charges", y = "Frequency") +
  theme_minimal() #distribution of bike rentals based on these weather groups.Rentals are more frequent on clear days (red), with peaks that are higher and more spread out across various rental counts, suggesting that clear weather is more conducive or appealing for bike rentals.Misty or cloudy days (blue) show a moderate frequency of rentals, with the distribution spread broadly across different rental counts, indicating a reasonable level of bike usage despite less ideal weather.Light snow or rain (green) exhibits the lowest frequency, mostly concentrated at the lower end of the rental counts, which implies significantly reduced bike usage during these conditions.

# The following section is to understand the mean of Total, casual and registered rentals per various variables

Medical_costs %>% 
  group_by(sex) %>%
  count(smoker) %>%
  mutate(prop = n/sum(n), pct = prop*100) %>%
  print # allows me to understand % of smokers by sex

## # A tibble: 4 × 5
## # Groups:   sex [2]
##   sex    smoker     n  prop   pct
##   <chr>  <chr>  <int> <dbl> <dbl>
## 1 female no       547 0.826  82.6
## 2 female yes      115 0.174  17.4
## 3 male   no       516 0.764  76.4
## 4 male   yes      159 0.236  23.6

Medical_costs %>%
  group_by(smoker) %>%
  summarise(mean_charges = mean(charges),
            mean_bmi = mean(bmi, na.rm = TRUE),
            sd_charges = sd(charges, na.rm = TRUE),
            sd_bmi = sd(bmi, na.rm = TRUE)) %>%
  print # summary statistics of the Medical_costs data broken down by smoking status

## # A tibble: 2 × 5
##   smoker mean_charges mean_bmi sd_charges sd_bmi
##   <chr>         <dbl>    <dbl>      <dbl>  <dbl>
## 1 no            8441.     30.7      5993.   6.05
## 2 yes          32050.     30.7     11542.   6.32

Medical_costs %>%
  group_by(sex) %>%
  summarise(mean_charges = mean(charges),
            mean_bmi = mean(bmi, na.rm = TRUE),
            sd_charges = sd(charges, na.rm = TRUE),
            sd_bmi = sd(bmi, na.rm = TRUE)) %>%
  print # summary statistics of the Medical_costs data broken down by sex

## # A tibble: 2 × 5
##   sex    mean_charges mean_bmi sd_charges sd_bmi
##   <chr>         <dbl>    <dbl>      <dbl>  <dbl>
## 1 female       12570.     30.4     11129.   6.05
## 2 male         13975.     30.9     12972.   6.14

Medical_costs %>%
  group_by(children) %>%
  summarise(mean_charges = mean(charges),
            sd_charges = sd(charges, na.rm = TRUE)) %>%
  print # summary statistics of the Medical_costs data broken down by number of children

## # A tibble: 6 × 3
##   children mean_charges sd_charges
##      <dbl>        <dbl>      <dbl>
## 1        0       12385.     12025.
## 2        1       12731.     11824.
## 3        2       15074.     12891.
## 4        3       15355.     12331.
## 5        4       13851.      9139.
## 6        5        8786.      3808.

Medical_costs %>%
  group_by(region) %>%
  summarise(mean_charges = mean(charges),
            mean_bmi = mean(bmi, na.rm = TRUE),
            sd_charges = sd(charges, na.rm = TRUE),
            sd_bmi = sd(bmi, na.rm = TRUE)) %>%
  print # summary statistics of the Medical_costs data broken down by region

## # A tibble: 4 × 5
##   region    mean_charges mean_bmi sd_charges sd_bmi
##   <chr>            <dbl>    <dbl>      <dbl>  <dbl>
## 1 northeast       13406.     29.2     11256.   5.94
## 2 northwest       12451.     29.2     11073.   5.14
## 3 southeast       14735.     33.4     13971.   6.48
## 4 southwest       12347.     30.6     11557.   5.69

Medical_costs %>%
  group_by(obs) %>%
  summarise(mean_charges = mean(charges),
            mean_bmi = mean(bmi, na.rm = TRUE),
            sd_charges = sd(charges, na.rm = TRUE),
            sd_bmi = sd(bmi, na.rm = TRUE)) %>%
  print

## # A tibble: 4 × 5
##   obs         mean_charges mean_bmi sd_charges sd_bmi
##   <chr>              <dbl>    <dbl>      <dbl>  <dbl>
## 1 normal            10409.     22.6      7506.  1.64 
## 2 obese             15572.     35.3     14553.  4.07 
## 3 overweight        10988.     27.6      8039.  1.43 
## 4 underweight        8852.     17.6      7735.  0.608

Medical_costs %>%
  ggplot(aes(x = charges, fill = smoker)) +
  geom_histogram(position = "identity", alpha = 0.7, bins = 30) +
  labs(x = "Cost", y = "Frequency") +
  theme_minimal() #We can see smokers have higher Medical_costs costs than non-smokers.

Medical_costs %>%
  ggplot(aes(x = charges, fill = factor(region))) +
  geom_histogram(position = "identity", alpha = 0.7, bins = 30) +
  labs(x = "Cost", y = "Frequency") +
  theme_minimal() # region is not showing much variance in the distribution of costs

Medical_costs %>% ggplot( 
  aes(x = bmi, 
      y = charges,
      colour = sex,
      shape = smoker,
      Size = region)) +
  geom_point() # again can see smokers even by gender incur higher costs.

plot(charges ~ age, data = Medical_costs,
     main = "Distribution of charge by age",
     xlab = "age",
     ylab = "charges",
     col = "#6138CE")

ggplot(Medical_costs,aes(x=age, y = charges)) +
  geom_point() +
  geom_smooth(method = 'lm', formula = y~x)# some indication that as you get older you incur higher costs.

ggplot(data = Medical_costs, aes_string(x = "age", y = "charges", col = "smoker" )) +
  geom_jitter() + 
  scale_color_manual(values=c("#FCC063" , "#236C19" )) +
  geom_smooth(method = 'lm') +
  ggtitle(glue::glue("Charges vs age"))  # more indications that you will incur higher charges as a smoke and older age.

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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

ggplot(data = Medical_costs, aes_string(x = 'bmi', y = 'charges', col = 'smoker')) + 
  geom_jitter() + 
  scale_color_manual(values=c("#FCC063" , "#236C19" )) +
  geom_smooth(method = 'lm') +
  ggtitle(glue::glue("Charges vs bmi")) # can see a relationship between a higher BMI and smokers

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

Fig_9 <- Medical_costs %>%  
  ggplot(aes(y = charges, x = sex,
             fill = sex)) +
  geom_boxplot() + 
  ylab("MonthlyCharges") +
  theme_bw()
ggplotly(Fig_9) # boxplot visual representation of charges by sex and slightly higher costs for males. We can also see outliers in charges for both Females & Males.

fig_7 <- Medical_costs %>% # facet grid to understand relationship between charges, age, smoking status and sex.
  ggplot(aes(age,charges))+geom_point(color="steelblue")+
  facet_grid(sex ~ smoker)
ggplotly(fig_7) 

Medical_costs %>% 
  ggplot(aes(y = charges, x = smoker,fill = smoker)) +
  geom_boxplot() +
  xlab ("") +
  theme_bw() # boxplot of charges and smoker where we can see some outliers for non-smokers.

Medical_costs %>% 
  ggplot(aes(y = charges, x = sex,fill = sex)) +
  geom_boxplot() +
  xlab ("") +
  theme_bw() # boxplot of charges and sex where both genders have a wide range of medical charges, males tend to have slightly higher median charges compared to females, and both groups exhibit several outliers indicating exceptionally high charges for some individuals and other factors influencing the charges suchas smoking and/or BMI.

Medical_costs %>% 
  ggplot(aes(age,charges))+geom_point(color="blue")+
  facet_grid(sex ~ obs) # facet grid to understand relationship between charges, bmi, age and sex.

Medical_costs %>% 
  group_by(children,sex, smoker) %>% 
  summarise(
    avg_cost = mean(charges),
    avg_bmi = mean(bmi)) %>%
  print # summary of mean charges and bmi by number of children, sex and smoking status

## `summarise()` has grouped output by 'children', 'sex'. You can override using
## the `.groups` argument.

## # A tibble: 22 × 5
## # Groups:   children, sex [12]
##    children sex    smoker avg_cost avg_bmi
##       <dbl> <chr>  <chr>     <dbl>   <dbl>
##  1        0 female no        7688.    30.4
##  2        0 female yes      30685.    30.3
##  3        0 male   no        7557.    30.7
##  4        0 male   yes      31902.    30.8
##  5        1 female no        8759.    30.2
##  6        1 female yes      30264.    29.2
##  7        1 male   no        7837.    30.9
##  8        1 male   yes      32905.    32.0
##  9        2 female no        9922.    31.0
## 10        2 female yes      31664.    29.2
## # ℹ 12 more rows

Fig_3 <- Medical_costs %>% 
  ggplot(aes(bmi,charges))+geom_point(color="blue")+
  facet_grid(sex ~ smoker) # facet grid to understand relationship between monthly charges, tenure, payment method and churn
ggplotly(Fig_3)

ggplot(Medical_costs, aes(x = charges)) +
  geom_histogram(binwidth = 1000, fill = "lightcoral", color = "black") +
  labs(title = "Distribution of Charges",
       x = "Charges",
       y = "Frequency") #distributin of charges

ggplot(Medical_costs, aes(x = bmi)) +
  geom_histogram(binwidth = 1, fill = "lightgreen", color = "black") +
  labs(title = "Distribution of BMI",
       x = "BMI",
       y = "Frequency") #distribution of BMI

Fig12 <- Medical_costs %>% ggplot(aes_string(x = 'obs', y = 'charges', col = 'smoker')) + 
  geom_jitter() + 
  scale_color_manual(values=c("#FCC063" , "#236C19" )) +
  geom_smooth(method = 'lm') +
  ggtitle(glue::glue("Charges vs bmi"))  
ggplotly(Fig12) # clearly see higher charges for obese smokers

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

Remove outliers from non-smoker charges

Medical_costs_outliers <- Medical_costs %>%
  group_by(smoker) %>%
  mutate(
    Q1 = quantile(charges, 0.25),
    Q3 = quantile(charges, 0.75),
    IQR = Q3 - Q1,
    Lower_Bound = Q1 - 1.5 * IQR,
    Upper_Bound = Q3 + 1.5 * IQR
  ) %>%
  filter(charges < Lower_Bound | charges > Upper_Bound) %>%
  select(smoker, charges, Lower_Bound, Upper_Bound) %>%
  distinct()

# Print identified outliers
View(Medical_costs_outliers)

out_nosmoker <- boxplot.stats(Medical_costs$charges[Medical_costs$smoker=='no'])$out

#Merge two lists of outliers and get indices:
out_ind <- which(
  ((Medical_costs$charges %in% c(out_nosmoker)) & (Medical_costs$smoker == 'no'))
)

Medical_filt <- Medical_costs[-c(out_ind),]

Medical_filt %>%
  group_by(smoker) %>%
  summarise(mean_charges = mean(charges),
            mean_bmi = mean(bmi, na.rm = TRUE),
            sd_charges = sd(charges, na.rm = TRUE),
            sd_bmi = sd(bmi, na.rm = TRUE)) %>%
  print

## # A tibble: 2 × 5
##   smoker mean_charges mean_bmi sd_charges sd_bmi
##   <chr>         <dbl>    <dbl>      <dbl>  <dbl>
## 1 no            7584.     30.6      4477.   6.09
## 2 yes          32050.     30.7     11542.   6.32

Medical_filt %>% 
  ggplot(aes(y = charges, x = smoker,fill = smoker)) +
  geom_boxplot() +
  xlab ("") +
  theme_bw() # boxplot of charges and smoker with outliers removed for non-smokers.

# Using GGpairs to understand the relationships between variables

# Key Insights:
# Charges and Age: There is a positive correlation suggesting that older individuals tend to incur higher medical charges.
# Charges and Smoker: The box plots show a distinct difference in charges between smokers and non-smokers, with smokers generally facing higher medical charges, indicating a strong effect of smoking on medical costs.
# BMI and Sex: The box plots comparing bmi across sex show slight differences in distribution which might be statistically significant.

ggpairs(Medical_costs, columns = c("charges", "age", "sex", "bmi", "smoker","region"))

ggpairs(Medical_filt, columns = c("charges", "age", "sex", "bmi", "smoker","region"))

The following code is a number of statistical tests using Hypothesis tests, correlation, and regression modelling to understand the strength of the relationships and assess the impact of each independent variable as a predictor of Medical_costs charges.

charges_smoker <- Medical_filt$charges[Medical_filt$smoker == "yes"]
charges_non_smoker <- Medical_filt$charges[Medical_filt$smoker == "no"]
#two-sample t-test to see if there is significant difference in charges between smokers and non-smokers 
t.test(charges_smoker, charges_non_smoker) 

## 
##  Welch Two Sample t-test
## 
## data:  charges_smoker and charges_non_smoker
## t = 34.399, df = 295.45, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  23066.16 25865.64
## sample estimates:
## mean of x mean of y 
##  32050.23   7584.33

#  The p-value is less than 2.2e-16, which is essentially zero. This indicates that the observed difference in means between the two groups (smokers and non-smokers) is highly statistically significant. It's so low that the probability of observing this large of a difference by chance is virtually nonexistent under the null hypothesis (which states that there is no difference in means) The 95% confidence interval for the difference in mean charges between smokers and non-smokers ranges from 23,066.16 to 25,865.64. This interval does not include zero, further supporting that there is a significant difference between the two groups. It suggests that, on average, smokers incur between €23,066.16 and €25,865.64 more in medical charges than non-smokers.


charges_female <- Medical_filt$charges[Medical_filt$sex == "female"]
charges_male <- Medical_filt$charges[Medical_filt$sex == "male"]
t.test(charges_female, charges_male) 

## 
##  Welch Two Sample t-test
## 
## data:  charges_female and charges_male
## t = -2.3908, df = 1262, p-value = 0.01696
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2897.0617  -285.4859
## sample estimates:
## mean of x mean of y 
##  11970.82  13562.10

#two-sample t-test to see if there is significant difference in charges between females and males which indicate a significant difference in the mean charges between female and male patients, with a p-value of 0.01696. This suggests that the average medical charges for males are higher than for females, with the estimated difference in means ranging from approximately $285 to $2897 higher for males

Pearson correlation tests

# correlation test to understand the relationship between bmi and charges
cor(Medical_costs$bmi, Medical_costs$charges, method = "pearson") 

## [1] 0.1984008

# correlation test to understand the relationship between age and charges
cor(Medical_costs$age, Medical_costs$charges, method = "pearson")

## [1] 0.2983082

# correlation test to understand the relationship between number of children and charges
cor(Medical_costs$children, Medical_costs$charges, method = "pearson")

## [1] 0.06738935

Regression model, initially looking at Total rentals which is a combination of Casual and Registered rentals.

# Multiple linear regression model to understand the relationship between all independent variables on Medical_costs costs
lm_model <- lm(charges ~ age + bmi + children + smoker + region + sex, data = Medical_costs)
summary(lm_model)

## 
## Call:
## lm(formula = charges ~ age + bmi + children + smoker + region + 
##     sex, data = Medical_costs)
## 
## 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  < 2e-16 ***
## age                256.76      11.91  21.555  < 2e-16 ***
## bmi                339.25      28.61  11.857  < 2e-16 ***
## children           474.82     137.90   3.443 0.000593 ***
## smokeryes        23847.33     413.35  57.693  < 2e-16 ***
## 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 *  
## sexmale           -129.48     333.20  -0.389 0.697630    
## ---
## 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: < 2.2e-16

lm_filt <- lm(charges ~ age + bmi + children + smoker + region + sex, data = Medical_filt)
summary(lm_filt)

## 
## Call:
## lm(formula = charges ~ age + bmi + children + smoker + region + 
##     sex, data = Medical_filt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11628.4  -2247.5   -310.2   1834.3  29957.4 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -12291.19     855.96 -14.360  < 2e-16 ***
## age                239.18      10.39  23.022  < 2e-16 ***
## bmi                346.67      24.79  13.984  < 2e-16 ***
## children           358.55     119.52   3.000  0.00275 ** 
## smokeryes        24606.53     354.55  69.402  < 2e-16 ***
## regionnorthwest   -292.75     415.89  -0.704  0.48162    
## regionsoutheast   -917.28     417.49  -2.197  0.02819 *  
## regionsouthwest   -595.06     415.20  -1.433  0.15205    
## sexmale            -50.96     289.55  -0.176  0.86033    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5177 on 1282 degrees of freedom
## Multiple R-squared:  0.8152, Adjusted R-squared:  0.814 
## F-statistic: 706.7 on 8 and 1282 DF,  p-value: < 2.2e-16

# Multiple linear regression model to understand the relationship between all independent variables on Medical_costs costs but now replacing BMI with BMI category
lm_model_obs <- lm(charges ~ age + obs + children + smoker + region + sex, data = Medical_costs)
summary(lm_model_obs)

## 
## Call:
## lm(formula = charges ~ age + obs + children + smoker + region + 
##     sex, data = Medical_costs)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -12342  -3590   -156   1569  27590 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4402.66     688.90  -6.391 2.28e-10 ***
## age               258.06      11.83  21.821  < 2e-16 ***
## obsobese         4747.68     469.14  10.120  < 2e-16 ***
## obsoverweight     837.64     506.89   1.653 0.098667 .  
## obsunderweight  -1211.79    1409.24  -0.860 0.390007    
## children          485.51     136.95   3.545 0.000406 ***
## smokeryes       23856.44     410.57  58.106  < 2e-16 ***
## regionnorthwest  -435.95     473.86  -0.920 0.357740    
## regionsoutheast  -658.15     468.39  -1.405 0.160221    
## regionsouthwest  -943.30     475.47  -1.984 0.047472 *  
## sexmale          -153.58     330.92  -0.464 0.642644    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6021 on 1326 degrees of freedom
## Multiple R-squared:  0.7546, Adjusted R-squared:  0.7528 
## F-statistic: 407.8 on 10 and 1326 DF,  p-value: < 2.2e-16

# Multiple linear regression model to understand the relationship between age and Medical_costs costs
lm_model_age <- lm(charges ~ age, data = Medical_costs)
summary(lm_model_age)

## 
## Call:
## lm(formula = charges ~ age, data = Medical_costs)
## 
## 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  < 2e-16 ***
## ---
## 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: < 2.2e-16

# Multiple linear regression model to understand the relationship between number of children and Medical_costs costs
lm_model_children <- lm(charges ~ children, data = Medical_costs)
summary(lm_model_children)

## 
## Call:
## lm(formula = charges ~ children, data = Medical_costs)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -11587  -8760  -4080   3471  51233 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  12537.4      446.8  28.061   <2e-16 ***
## children       676.9      274.3   2.468   0.0137 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12090 on 1335 degrees of freedom
## Multiple R-squared:  0.004541,   Adjusted R-squared:  0.003796 
## F-statistic:  6.09 on 1 and 1335 DF,  p-value: 0.01372

# Multiple linear regression model to understand the relationship between smoking status and Medical_costs costs
lm_model_smoker <- lm(charges ~ smoker, data = Medical_costs)
summary(lm_model_smoker)

## 
## Call:
## lm(formula = charges ~ smoker, data = Medical_costs)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -19221  -5048   -923   3702  31720 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8440.7      229.1   36.84   <2e-16 ***
## smokeryes    23609.6      506.2   46.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7471 on 1335 degrees of freedom
## Multiple R-squared:  0.6197, Adjusted R-squared:  0.6195 
## F-statistic:  2176 on 1 and 1335 DF,  p-value: < 2.2e-16

# Multiple linear regression model to understand the relationship between smoking status and Medical_costs costs
lm_model_smoker_withbmi <- lm(charges ~ smoker*bmi, age, data = Medical_filt)
summary(lm_model_smoker_withbmi)

## 
## Call:
## lm(formula = charges ~ smoker * bmi, data = Medical_filt, subset = age)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -6869  -3488  -2233   3811  12091 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     2324.52     807.12   2.880  0.00404 ** 
## smokeryes     -22977.81    1849.34 -12.425  < 2e-16 ***
## bmi              166.62      25.67   6.492  1.2e-10 ***
## smokeryes:bmi   1542.15      55.57  27.753  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4618 on 1287 degrees of freedom
## Multiple R-squared:  0.9034, Adjusted R-squared:  0.9032 
## F-statistic:  4014 on 3 and 1287 DF,  p-value: < 2.2e-16

lm_model_smoker_bmi <- lm(charges ~ smoker, bmi, age, data = Medical_filt)
summary(lm_model_smoker_bmi)

## 
## Call:
## lm(formula = charges ~ smoker, data = Medical_filt, subset = bmi, 
##     weights = age)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -124962  -33577  -23645   29032   78003 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10112.4      169.7   59.59   <2e-16 ***
## smokeryes    31492.5      344.5   91.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32210 on 1289 degrees of freedom
## Multiple R-squared:  0.8663, Adjusted R-squared:  0.8662 
## F-statistic:  8355 on 1 and 1289 DF,  p-value: < 2.2e-16

#visually represent the distribution of residuals of the smokers regression model
hist(lm_model_smoker$residuals) 

#visually represent the distribution of residuals of the collective regression model
hist(lm_model$residuals)

plot(lm_model_smoker)

Medical_costs$smoker <- as.factor(Medical_costs$smoker) 
model_anova <- aov(charges ~ smoker, data = Medical_costs)
summary(model_anova)

##               Df    Sum Sq   Mean Sq F value Pr(>F)    
## smoker         1 1.214e+11 1.214e+11    2176 <2e-16 ***
## Residuals   1335 7.451e+10 5.581e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# The ANOVA results further reinforce the significant impact of smoking on medical charges. 
# F-value: The F-value of 2178 is very high, which strongly indicates that the mean charges between smokers and non-smokers are statistically significantly different. The higher the F-value, the more certain we can be that any observed differences in means across groups are not due to randomness.
# p-value: The p-value (< 2e-16) suggests the result is highly significant, meaning the likelihood that the observed difference in charges between smokers and non-smokers is due to chance is extremely low.

# ANCOVA model
model_ancova <- aov(charges ~ smoker + age + bmi, data = Medical_costs)
summary(model_ancova)

##               Df    Sum Sq   Mean Sq F value Pr(>F)    
## smoker         1 1.214e+11 1.214e+11  3269.3 <2e-16 ***
## age            1 1.988e+10 1.988e+10   535.3 <2e-16 ***
## bmi            1 5.114e+09 5.114e+09   137.7 <2e-16 ***
## Residuals   1333 4.951e+10 3.714e+07                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#The ANCOVA results suggest very strong evidence for significant effects of smoking status, age, and BMI on medical charges, given the extremely low p-values (<2e-16) for each predictor.
#Smoker: The smoker variable is the most significant predictor, with an F-value of 3269.3. This indicates that there is a very strong relationship between whether a patient smokes and their medical charges, with smokers likely incurring significantly higher charges compared to non-smokers.
# Age: Age also significantly affects medical charges, with an F-value of 535.3. This suggests that as age increases, so do the medical charges, reflecting perhaps the increased medical needs and complications that often accompany aging.
#BMI: BMI is another significant predictor, with an F-value of 137.7. Higher BMI values, which indicate higher body fat, are associated with higher medical charges, potentially due to related health issues such as cardiovascular diseases, diabetes, and joint problems.

#Residuals: The residuals' line shows the variation in medical charges that is not explained by the model. Here, the mean squared error (Mean Sq) is approximately 37,140,000, which reflects the average squared difference between the observed medical charges and the charges predicted by the model.

#Overall, the model indicates that both smoking status, age, and BMI are key factors influencing medical charges, with smoking status being the most potent predictor among the three.

AIC(lm_model)

## [1] 27096.22

AIC(lm_model_smoker)

## [1] 27646.96

AIC(lm_model_smoker_bmi)

## [1] 25948.13

AIC(lm_model_smoker_withbmi)

## [1] 25455.88