Multiple linear regression is an extension of simple linear regression to predict an outcome based on multiple predictors. It is useful when you want to understand the relationship between a dependent variable and several independent variables. This method is widely used across many fields such as economics, social sciences, biology, and engineering to model and predict outcomes.

The general form of the MLR model is:

\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \dots + \beta_nX_n + \epsilon \]

where:

• \(\beta_0\) is the intercept,

• \(\beta_1, \beta_2, \dots, \beta_n\) are the coefficients of the predictors \(X_1, X_2, \dots, X_n\),

• \(\epsilon\) is the error term, assumed to be normally distributed with mean zero and constant variance.

1. Linearity: The relationship between the dependent and each independent variable must be linear.
2. No Multicollinearity: Independent variables should not be too highly correlated.
3. Independence of Errors: There should be no correlation between the residual (error) terms.
4. Homoscedasticity: The variance of error terms should be constant.
5. Normality of Residuals: The residuals should follow a normal distribution.

knitr::opts_chunk$set(echo = TRUE)


# Load necessary packages
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(readr)
library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

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

The data used in this analysis was obtained from Kaggle, specifically from the following dataset: Insurance Dataset.

data <- read.csv("insurance.csv")
head(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 Preparation

# Prepare the data by converting categorical variables to factors (if they are not already)
data$sex <- as.factor(data$sex)
data$smoker <- as.factor(data$smoker)
data$region <- as.factor(data$region)

# Summary of numerical attributes
summary(data$charges)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1122    4740    9382   13270   16640   63770

summary(data$age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   27.00   39.00   39.21   51.00   64.00

summary(data$bmi)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.96   26.30   30.40   30.66   34.69   53.13

# Box plots for categorical variables against charges
ggplot(data, aes(x = sex, y = charges)) + geom_boxplot() + labs(title = "Charges by Sex")

ggplot(data, aes(x = smoker, y = charges)) + geom_boxplot() + labs(title = "Charges by Smoking Status")

ggplot(data, aes(x = region, y = charges)) + geom_boxplot() + labs(title = "Charges by Region")

# Scatter plots for continuous variables
ggplot(data, aes(x = age, y = charges)) + geom_point() + geom_smooth(method = "lm") + labs(title = "Age vs Charges")

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

ggplot(data, aes(x = bmi, y = charges)) + geom_point() + geom_smooth(method = "lm") + labs(title = "BMI vs Charges")

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

ggplot(data, aes(x = charges)) + geom_histogram(binwidth = 1000, fill = "blue", color = "white") + labs(title = "Distribution of Insurance Charges")

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

## 
## Call:
## lm(formula = charges ~ age + bmi + children + sex + smoker + 
##     region, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11304.9  -2848.1   -982.1   1393.9  29992.8 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -11938.5      987.8 -12.086  < 2e-16 ***
## age                256.9       11.9  21.587  < 2e-16 ***
## bmi                339.2       28.6  11.860  < 2e-16 ***
## children           475.5      137.8   3.451 0.000577 ***
## sexmale           -131.3      332.9  -0.394 0.693348    
## smokeryes        23848.5      413.1  57.723  < 2e-16 ***
## regionnorthwest   -353.0      476.3  -0.741 0.458769    
## regionsoutheast  -1035.0      478.7  -2.162 0.030782 *  
## regionsouthwest   -960.0      477.9  -2.009 0.044765 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6062 on 1329 degrees of freedom
## Multiple R-squared:  0.7509, Adjusted R-squared:  0.7494 
## F-statistic: 500.8 on 8 and 1329 DF,  p-value: < 2.2e-16

coefs <- summary(model)$coefficients
eq <- paste("Charges = ", round(coefs[1, 1], 2), " + ", round(coefs[2, 1], 2), "*(Age) + ", round(coefs[3, 1], 2), "*(BMI) + ", round(coefs[4, 1], 2), "*(Children) + ...", sep="")

plot(model, which = 1)

# Categorizing BMI into recognizable health categories
data$bmi_category <- cut(data$bmi, 
                         breaks = c(0, 18.5, 24.9, 29.9, Inf), 
                         labels = c("Underweight", "Normal", "Overweight", "Obese"),
                         include.lowest = TRUE)

# Create the histogram
p <- plot_ly(data, x = ~bmi, color = ~bmi_category, colors = c("#ADD8E6", "#90EE90", "#FFD700", "#FF6347"),
             type = "histogram",
             histnorm = "percent",
             marker = list(line = list(color = '#000000', width = 2))) %>%
      layout(title = "Distribution of BMI Categories",
             xaxis = list(title = "BMI"),
             yaxis = list(title = "Percentage"),
             barmode = 'overlay',
             hovermode = 'closest')

p

This interactive histogram demonstrates the distribution of BMI across various health-related categories.

Model Diagnostics

par(mfrow=c(2,2))
plot(model)

coef_data <- as.data.frame(summary(model)$coefficients)
coef_data$Variable <- rownames(coef_data)

coef <- ggplot(coef_data, aes(x = reorder(Variable, Estimate), y = Estimate, fill = Estimate > 0)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  labs(x = "Variable", y = "Coefficient Value", title = "Impact of Variables on Insurance Charges") +
  theme_minimal()

coef

The multiple linear regression analysis revealed the following relationship between insurance charges and the explanatory variables:

\[ \text{Charges} = -11938.5 + 256.9 \times \text{Age} + 339.2 \times \text{BMI} + 475.5 \times \text{Children}\] \[- 131.3 \times \text{Sex}_{\text{male}} + 23848.5 \times \text{Smoker}_{\text{yes}} - 353.0 \] \[\times \text{Region}_{\text{nw}} - 1035 \times \text{Region}_{\text{sw}} -960 \times \text{Region}_{\text{sw}} \]

Interpretation

• Age, BMI, and Children positively affect insurance charges.
• Sex (Male) and Regions (Northwest, Southeast, Southwest) are associated with a decrease in charges.
• Smoking is the most significant factor, drastically increasing charges.