Chapter 6: Regression Methods

Linear Regression: Challenger2 Dataset

launch <- read.csv("challenger2.csv")
head(launch)

 #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

cor(launch$temperature, launch$distress_ct)

[1] -0.3359996

# computing the slope using correlation
r * (sd(launch$distress_ct) / sd(launch$temperature))

[1] -0.03364796

model <- lm(distress_ct ~ temperature, data = launch)
model

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

Coefficients:
(Intercept)  temperature  
    2.81458     -0.03365  

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

# creating a simple multiple 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)
}

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 regression model with simple linear regression
reg(y = launch$distress_ct, x = launch[2])

               estimate
Intercept    2.81458456
temperature -0.03364796

# use regression model with multiple regression
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

# confirming the multiple regression result using the lm function (not in text)
model <- lm(distress_ct ~ temperature + field_check_pressure + flight_num, data = launch)
model

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

Coefficients:
         (Intercept)           temperature  field_check_pressure            flight_num  
           2.240e+00            -3.124e-02            -2.587e-05             2.762e-02  

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

Attempting to improve th model was not successful, as the R-squared decreased, therefore reducing the model’s explanatory power.

---

Predicting Medical Expenses: Using “insurance.csv”

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

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 

hist(insurance$expenses, main="Expenses: Histogram", xlab = "Expenses",border="white", col="red")

table(insurance$region)

northeast northwest southeast southwest 
      324       325       364       325 

# exploring relationships among features: correlation matrix
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

# visualing relationships among features: scatterplot matrix
pairs(insurance[c("age", "bmi", "children", "expenses")], col="yellow", main="Feature Relationship Visuals")

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

Evaluating Model Performance

summary(ins_model)

Call:
lm(formula = expenses ~ age + bmi + children + sex + smoker + 
    region, 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 ***
bmi                339.3       28.6  11.864  < 2e-16 ***
children           475.7      137.8   3.452 0.000574 ***
sexmale           -131.3      332.9  -0.395 0.693255    
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

# Create age squared
insurance$age2 <- insurance$age^2

# Create bmi where it is equal to 30
insurance$bmi30 <- ifelse(insurance$bmi >= 30, 1, 0)

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

Call:
lm(formula = expenses ~ age + age2 + bmi + sex + region + children + 
    bmi30 * smoker, 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 ***
bmi               119.7715    34.2796   3.494 0.000492 ***
sexmale          -496.7690   244.3713  -2.033 0.042267 *  
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 ***
children          678.6017   105.8855   6.409 2.03e-10 ***
bmi30            -997.9355   422.9607  -2.359 0.018449 *  
smokeryes       13404.5952   439.9591  30.468  < 2e-16 ***
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

ins_model without age feature

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-17401.6  -1642.4  -1267.9   -747.2  24244.5 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -408.6563   921.2668  -0.444 0.657418    
age2                3.3282     0.1088  30.580  < 2e-16 ***
bmi               119.0291    34.2435   3.476 0.000525 ***
sexmale          -495.6944   244.2987  -2.029 0.042652 *  
regionnorthwest  -277.7215   349.1801  -0.795 0.426550    
regionsoutheast  -826.6339   351.5459  -2.351 0.018847 *  
regionsouthwest -1222.1659   350.4386  -3.488 0.000503 ***
children          661.2122   100.9411   6.550 8.18e-11 ***
bmi30            -983.4079   422.0088  -2.330 0.019939 *  
smokeryes       13404.3214   439.8423  30.475  < 2e-16 ***
bmi30:smokeryes 19809.4575   604.5154  32.769  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4444 on 1327 degrees of freedom
Multiple R-squared:  0.8664,    Adjusted R-squared:  0.8654 
F-statistic: 860.3 on 10 and 1327 DF,  p-value: < 2.2e-16

The R-squared stayed the same for both models, however, the Adjusted R-squared for the model without age feature improved slightly

# Making prediction with in_model2
insurance$pred <- predict(ins_model2, insurance)
cor(insurance$pred, insurance$expenses)

[1] 0.9307838

plot(insurance$pred, insurance$expenses, col="blue", xlab = "Predicted Expenses",
     ylab = "Actual Expenses", main = "Actual vs. Predicted Expenses")
abline(a = 0, b = 1, col = "red", lwd = 2, lty = 2)

Case A: Age= 22, Children= 3, bmi= 24, sex= female, bmi30= 0, smoker= no, region= Northwest.

predict(ins_model2,
        data.frame(age2 = 22, children = 3,
                   bmi = 30, sex = "male", bmi30 = 0,
                   smoker = "no", region = "northwest"))

       1 
4445.657 

The prediction for case A demonstrated that the level of expense is 1 unit of measure, when all features are 0, and when computing the values for all features, expenses increases around 4445.657, on average.

Case B: Age= 22, Children= 1, bmi= 27, sex= male, bmi30= 0, smoker= yes, region= Southeast.

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

Similarly, the prediction for case B demonstrated that the level of expense is 1 unit of measure, when all features are 0, and when computing the values for all features, expenses increases around 15621.55, on average.