Acivity 7

Load Required Libraries

#install.packages("rpart")
#install.packages("rpart.plot")
#install.packages("Cubist")
library(rpart)
library(rpart.plot)
library(Cubist)

Loading required package: lattice

Part 1: Linear Regression

Understanding Regression

getwd()

[1] "/cloud/project"

# Load dataset
launch <- read.csv("challenger2.csv")

# Estimate beta manually
b <- cov(launch$temperature, launch$distress_ct) / var(launch$temperature)
b

[1] -0.03364796

# Estimate alpha manually
a <- mean(launch$distress_ct) - b * mean(launch$temperature)
a

[1] 2.814585

# Calculate the correlation of launch data
r <- cov(launch$temperature, launch$distress_ct) / (sd(launch$temperature) * sd(launch$distress_ct))
r

[1] -0.3359996

# Confirm correlation using cor function
cor(launch$temperature, launch$distress_ct)

[1] -0.3359996

# Compute the slope using correlation
slope <- r * (sd(launch$distress_ct) / sd(launch$temperature))
slope

[1] -0.03364796

# Regression model using lm function
model <- lm(distress_ct ~ temperature, data = launch)
summary(model)

Call:
lm(formula = distress_ct ~ temperature, data = launch)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0649 -0.4929 -0.2573  0.3052  1.7090 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  2.81458    1.24629   2.258   0.0322 *
temperature -0.03365    0.01815  -1.854   0.0747 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7076 on 27 degrees of freedom
Multiple R-squared:  0.1129,    Adjusted R-squared:  0.08004 
F-statistic: 3.436 on 1 and 27 DF,  p-value: 0.07474

# Predicting values and plotting the regression line
plot(launch$temperature, launch$distress_ct, main="Temperature vs. Distress Count")
abline(model, col="red")

Multiple Regression

# Define custom regression function
reg <- function(y, x) {
  x <- as.matrix(x)
  x <- cbind(Intercept = 1, x)
  b <- solve(t(x) %*% x) %*% t(x) %*% y
  colnames(b) <- "estimate"
  print(b)
}

# Examine launch data structure
str(launch)

'data.frame':   29 obs. of  4 variables:
 $ distress_ct         : int  0 1 0 0 0 0 0 0 1 1 ...
 $ temperature         : int  66 70 69 68 67 72 73 70 57 63 ...
 $ field_check_pressure: int  50 50 50 50 50 50 100 100 200 200 ...
 $ flight_num          : int  1 2 3 4 5 6 7 8 9 10 ...

# Test simple linear regression model
reg(y = launch$distress_ct, x = launch[2])

               estimate
Intercept    2.81458456
temperature -0.03364796

# Test multiple regression model
reg(y = launch$distress_ct, x = launch[2:4])

                          estimate
Intercept             2.239817e+00
temperature          -3.124185e-02
field_check_pressure -2.586765e-05
flight_num            2.762455e-02

# Confirm multiple regression using lm function
model <- lm(distress_ct ~ temperature + field_check_pressure + flight_num, data = launch)
summary(model)

Call:
lm(formula = distress_ct ~ temperature + field_check_pressure + 
    flight_num, data = launch)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2744 -0.3335 -0.1657  0.2975  1.5284 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)           2.240e+00  1.267e+00   1.767   0.0894 .
temperature          -3.124e-02  1.787e-02  -1.748   0.0927 .
field_check_pressure -2.587e-05  2.383e-03  -0.011   0.9914  
flight_num            2.762e-02  1.798e-02   1.537   0.1369  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6926 on 25 degrees of freedom
Multiple R-squared:  0.2132,    Adjusted R-squared:  0.1188 
F-statistic: 2.259 on 3 and 25 DF,  p-value: 0.1063

Predicting Medical Expenses

# Load insurance dataset
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 ...

# Summarize expenses variable
summary(insurance$expenses)

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

# Histogram of insurance expenses
hist(insurance$expenses)

# Table of region
table(insurance$region)

northeast northwest southeast southwest 
      324       325       364       325 

# Correlation matrix of selected variables
cor(insurance[c("age", "bmi", "children", "expenses")])

               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

# Scatterplot matrix of selected variables
pairs(insurance[c("age", "bmi", "children", "expenses")])

Training and Evaluating a Model

# Train linear regression model
ins_model <- lm(expenses ~ age + children + bmi + sex + smoker + region, data = insurance)
ins_model <- lm(expenses ~ ., data = insurance) # Equivalent shorthand

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

Improving Model Performance

# Add a higher-order "age" term
insurance$age2 <- insurance$age^2

# Add an indicator for BMI >= 30
insurance$bmi30 <- ifelse(insurance$bmi >= 30, 1, 0)

# Create final improved model
ins_model2 <- lm(expenses ~ age + age2 + children + bmi + sex + smoker + region + bmi30, data = insurance)
summary(ins_model2)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-12493.8  -3360.4    136.9   1310.8  29317.0 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -2948.272   1828.553  -1.612 0.107123    
age               -28.481     80.443  -0.354 0.723357    
age2                3.602      1.004   3.590 0.000343 ***
children          630.455    142.365   4.428 1.03e-05 ***
bmi               154.070     46.073   3.344 0.000849 ***
sexmale          -166.288    328.313  -0.506 0.612595    
smokeryes       23857.060    407.349  58.567  < 2e-16 ***
regionnorthwest  -400.413    469.636  -0.853 0.394033    
regionsoutheast  -888.874    472.837  -1.880 0.060344 .  
regionsouthwest  -947.344    471.206  -2.010 0.044584 *  
bmi30            2725.683    547.814   4.976 7.36e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5977 on 1327 degrees of freedom
Multiple R-squared:  0.7582,    Adjusted R-squared:  0.7564 
F-statistic: 416.2 on 10 and 1327 DF,  p-value: < 2.2e-16