R Notebook

Use the MLR solution built during class to predict the insurance quote given the following case scenarios:

Case 1: Age=22, Children=3,bmi=24,sex=female,bmi30=0,smoker=no, region=Northwest. Case 1: Age=22, Children=1,bmi=27,sex=male,bmi30=0,smoker=yes, region=Southeast.

exploring insurance data

insurance <- read.csv("insurance.csv", stringsAsFactors = TRUE)
str(insurance)

'data.frame':   1338 obs. of  7 variables:
 $ age     : int  19 18 28 33 32 31 46 37 37 60 ...
 $ sex     : Factor w/ 2 levels "female","male": 1 2 2 2 2 1 1 1 2 1 ...
 $ bmi     : num  27.9 33.8 33 22.7 28.9 25.7 33.4 27.7 29.8 25.8 ...
 $ children: int  0 1 3 0 0 0 1 3 2 0 ...
 $ smoker  : Factor w/ 2 levels "no","yes": 2 1 1 1 1 1 1 1 1 1 ...
 $ region  : Factor w/ 4 levels "northeast","northwest",..: 4 3 3 2 2 3 3 2 1 2 ...
 $ expenses: num  16885 1726 4449 21984 3867 ...

summary(insurance$expenses)

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

histogram of insurance charges

hist(insurance$expenses)

table(insurance$region)

northeast northwest southeast southwest 
      324       325       364       325 

cor(insurance[c("age", "bmi", "children", "expenses")])

               age        bmi   children
age      1.0000000 0.10934101 0.04246900
bmi      0.1093410 1.00000000 0.01264471
children 0.0424690 0.01264471 1.00000000
expenses 0.2990082 0.19857626 0.06799823
           expenses
age      0.29900819
bmi      0.19857626
children 0.06799823
expenses 1.00000000

pairs(insurance[c("age", "bmi", "children", "expenses")])

training model on the data

ins_model <- lm(expenses ~ age + children + bmi + sex + smoker + region,
                data = insurance)
ins_model <- lm(expenses ~ ., data = insurance) 
ins_model

Call:
lm(formula = expenses ~ ., data = insurance)

Coefficients:
    (Intercept)              age  
       -11941.6            256.8  
        sexmale              bmi  
         -131.4            339.3  
       children        smokeryes  
          475.7          23847.5  
regionnorthwest  regionsoutheast  
         -352.8          -1035.6  
regionsouthwest  
         -959.3  

summary(ins_model)

Call:
lm(formula = expenses ~ ., data = insurance)

Residuals:
     Min       1Q   Median       3Q      Max 
-11302.7  -2850.9   -979.6   1383.9  29981.7 

Coefficients:
                Estimate Std. Error t value
(Intercept)     -11941.6      987.8 -12.089
age                256.8       11.9  21.586
sexmale           -131.3      332.9  -0.395
bmi                339.3       28.6  11.864
children           475.7      137.8   3.452
smokeryes        23847.5      413.1  57.723
regionnorthwest   -352.8      476.3  -0.741
regionsoutheast  -1035.6      478.7  -2.163
regionsouthwest   -959.3      477.9  -2.007
                Pr(>|t|)    
(Intercept)      < 2e-16 ***
age              < 2e-16 ***
sexmale         0.693255    
bmi              < 2e-16 ***
children        0.000574 ***
smokeryes        < 2e-16 ***
regionnorthwest 0.458976    
regionsoutheast 0.030685 *  
regionsouthwest 0.044921 *  
---
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.9 on 8 and 1329 DF,  p-value: < 2.2e-16

improving model performance, adding age term, BMI

insurance$age2 <- insurance$age^2

insurance$bmi30 <- ifelse(insurance$bmi >= 30, 1, 0)

creating final model

ins_model2 <- lm(expenses ~ age + age2 + children + bmi + sex +
                   bmi30*smoker + region, data = insurance)

summary(ins_model2)

Call:
lm(formula = expenses ~ age + age2 + children + bmi + sex + bmi30 * 
    smoker + region, data = insurance)

Residuals:
     Min       1Q   Median       3Q      Max 
-17297.1  -1656.0  -1262.7   -727.8  24161.6 

Coefficients:
                  Estimate Std. Error t value
(Intercept)       139.0053  1363.1359   0.102
age               -32.6181    59.8250  -0.545
age2                3.7307     0.7463   4.999
children          678.6017   105.8855   6.409
bmi               119.7715    34.2796   3.494
sexmale          -496.7690   244.3713  -2.033
bmi30            -997.9355   422.9607  -2.359
smokeryes       13404.5952   439.9591  30.468
regionnorthwest  -279.1661   349.2826  -0.799
regionsoutheast  -828.0345   351.6484  -2.355
regionsouthwest -1222.1619   350.5314  -3.487
bmi30:smokeryes 19810.1534   604.6769  32.762
                Pr(>|t|)    
(Intercept)     0.918792    
age             0.585690    
age2            6.54e-07 ***
children        2.03e-10 ***
bmi             0.000492 ***
sexmale         0.042267 *  
bmi30           0.018449 *  
smokeryes        < 2e-16 ***
regionnorthwest 0.424285    
regionsoutheast 0.018682 *  
regionsouthwest 0.000505 ***
bmi30:smokeryes  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4445 on 1326 degrees of freedom
Multiple R-squared:  0.8664,    Adjusted R-squared:  0.8653 
F-statistic: 781.7 on 11 and 1326 DF,  p-value: < 2.2e-16

making predictions with the regression model

insurance$pred <- predict(ins_model2, insurance)
cor(insurance$pred, insurance$expenses)

[1] 0.9307999

plot(insurance$pred, insurance$expenses)
abline(a = 0, b = 1, col = "red", lwd = 3, lty = 2)

Case 1

predict(ins_model2,
        data.frame(age = 22, age2 = 22^2, children = 3,
                   bmi = 24, sex = "female", bmi30 = 0,
                   smoker = "no", region = "northwest"))

       1 
5858.241 

insurance cost

Case 2

predict(ins_model2,
        data.frame(age = 22, age2 = 22^2, children = 1,
                   bmi = 27, sex = "male", bmi30 = 0,
                   smoker = "yes", region = "southeast"))

       1 
17219.31 

insurance cost

tee <- c(1, 1, 1, 2, 2, 3, 4, 5, 5, 6, 6, 7, 7, 7, 7)
at1 <- c(1, 1, 1, 2, 2, 3, 4, 5, 5)
at2 <- c(6, 6, 7, 7, 7, 7)
bt1 <- c(1, 1, 1, 2, 2, 3, 4)
bt2 <- c(5, 5, 6, 6, 7, 7, 7, 7)

sdr_a <- sd(tee) - (length(at1) / length(tee) * sd(at1) + length(at2) / length(tee) * sd(at2))
sdr_b <- sd(tee) - (length(bt1) / length(tee) * sd(bt1) + length(bt2) / length(tee) * sd(bt2))

sdr_a

[1] 1.202815

sdr_b

[1] 1.392751