#Question A: 5858.241 #Question B: 17219.31

#Predicting medical expenses

insurance <- read.csv("insurance.csv", stringsAsFactors = TRUE) #adding insurance csv
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) #sumamry of main variables
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1122    4740    9382   13270   16640   63770
hist(insurance$expenses) #making histogram to visualize expenses

table(insurance$region) #making a table of the region
## 
## northeast northwest southeast southwest 
##       324       325       364       325
cor(insurance[c("age", "bmi", "children", "expenses")]) #a correlation matrix to see how the different variables correlate with each other
##                age        bmi   children   expenses
## age      1.0000000 0.10934101 0.04246900 0.29900819
## bmi      0.1093410 1.00000000 0.01264471 0.19857626
## children 0.0424690 0.01264471 1.00000000 0.06799823
## expenses 0.2990082 0.19857626 0.06799823 1.00000000
pairs(insurance[c("age", "bmi", "children", "expenses")]) #visualizing how the different variabvles relate to each other

ins_model <- lm(expenses ~ age +children+bmi+sex+smoker+region, 
                data = insurance) #making the model
ins_model<- lm(expenses ~ ., data = insurance) #same as above model
ins_model #shows estimated beta coefficants
## 
## Call:
## lm(formula = expenses ~ ., data = insurance)
## 
## Coefficients:
##     (Intercept)              age          sexmale              bmi  
##        -11941.6            256.8           -131.4            339.3  
##        children        smokeryes  regionnorthwest  regionsoutheast  
##           475.7          23847.5           -352.8          -1035.6  
## regionsouthwest  
##          -959.3
summary(ins_model) #more details, many of the variables have p values less than 0.05
## 
## 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 Pr(>|t|)    
## (Intercept)     -11941.6      987.8 -12.089  < 2e-16 ***
## age                256.8       11.9  21.586  < 2e-16 ***
## sexmale           -131.3      332.9  -0.395 0.693255    
## bmi                339.3       28.6  11.864  < 2e-16 ***
## children           475.7      137.8   3.452 0.000574 ***
## smokeryes        23847.5      413.1  57.723  < 2e-16 ***
## regionnorthwest   -352.8      476.3  -0.741 0.458976    
## regionsoutheast  -1035.6      478.7  -2.163 0.030685 *  
## regionsouthwest   -959.3      477.9  -2.007 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
insurance$age2 <- insurance$age^2 #squares the age to improve the results
insurance$bmi30 <- ifelse(insurance$bmi >= 30,1,0) #indicator if greater or equal to 30
ins_model2 <- lm(expenses ~ age+age2+children+bmi+sex+bmi30*smoker+region, data=insurance) #final model
summary(ins_model2) #more info on updated model, still a good amount of p-values less than 0.05
## 
## 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 Pr(>|t|)    
## (Intercept)       139.0053  1363.1359   0.102 0.918792    
## age               -32.6181    59.8250  -0.545 0.585690    
## age2                3.7307     0.7463   4.999 6.54e-07 ***
## children          678.6017   105.8855   6.409 2.03e-10 ***
## bmi               119.7715    34.2796   3.494 0.000492 ***
## sexmale          -496.7690   244.3713  -2.033 0.042267 *  
## bmi30            -997.9355   422.9607  -2.359 0.018449 *  
## smokeryes       13404.5952   439.9591  30.468  < 2e-16 ***
## regionnorthwest  -279.1661   349.2826  -0.799 0.424285    
## regionsoutheast  -828.0345   351.6484  -2.355 0.018682 *  
## regionsouthwest -1222.1619   350.5314  -3.487 0.000505 ***
## bmi30:smokeryes 19810.1534   604.6769  32.762  < 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
insurance$pred <- predict(ins_model2, insurance)
cor(insurance$pred, insurance$expenses) #shows high correlation with model 2 and expenses
## [1] 0.9307999
plot(insurance$pred, insurance$expenses)
abline(a = 0, b = 1, col = "red" , lwd = 3, lty = 2) #plots the correlation with a best fit line

predict(ins_model2, 
        data.frame(age= 30, age2= 30^2, children= 2, bmi = 30, sex = "male", bmi30 = 1, smoker = "no", region = "northeast")) #estimating the insurance expenses for a person who has the above variables
##        1 
## 5973.774
predict(ins_model2, 
        data.frame(age= 30, age2= 30^2, children= 2, bmi = 30, sex = "female", bmi30 = 1, smoker = "no", region = "northeast")) #by just changing the sex to female, the insurance expense increase by around $500
##        1 
## 6470.543
predict(ins_model2, 
        data.frame(age= 30, age2= 30^2, children= 0, bmi = 30, sex = "female", bmi30 = 1, smoker = "no", region = "northeast")) #insurance expenses are much lower with no children for this person
##       1 
## 5113.34

#Question A

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

#Question B

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

#despite having the same age, less children, person B had a much higher insurance expenses. This could be because he had a higher bmi, he is a male, he is from a different region, and he smokes

#regression trees and model trees

tee<-c(1,1,1,2,2,3,4,5,5,6,6,7,7,7,7) #setting up the data
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)) #making SDR
sdr_b <- sd(tee)-(length(bt1)/length(tee)*sd(bt1) +length(bt2)/length(tee) *sd(bt2))
sdr_a #finding sdr of a
## [1] 1.202815
sdr_b #sdr_b had a higher value than a
## [1] 1.392751