R Notebook

getwd() launch <- read.csv(“challenger.csv”)

estimate beta manually

b <- cov(launch\(temperature, launch\)distress_ct) / var(launch$temperature) b

estimate alpha manually

a <- mean(launch\(distress_ct) - b * mean(launch\)temperature) a

r <- cov(launch\(temperature, launch\)distress_ct) / (sd(launch\(temperature) * sd(launch\)distress_ct)) r

cor(launch\(temperature, launch\)distress_ct)

computing the slope using correlation

r * (sd(launch\(distress_ct) / sd(launch\)temperature))

confirming the regression line using the lm function (not in text)

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)

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) }

examine the launch data

str(launch)

test regression model with simple linear regression

reg(y = launch$distress_ct, x = launch[2])

use regression model with multiple regression

reg(y = launch$distress_ct, x = launch[2:4])

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

2.240e+00 -3.124e-02 -2.587e-05

flight_num

2.762e-02

summary(model)

Step 2: Exploring and preparing the data —-

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

summarize the charges variable

summary(insurance$expenses)

histogram of insurance charges

hist(insurance$expenses)

table of region

table(insurance$region)

exploring relationships among features: correlation matrix

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

visualing relationships among features: scatterplot matrix

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

Step 3: Training a model on the data —-

ins_model <- lm(expenses ~ age + children + bmi + sex + smoker + region, data = insurance) ins_model <- lm(expenses ~ ., data = insurance) # this is equivalent to above

see the estimated beta coefficients

ins_model

see more detail about the estimated beta coefficients

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

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

making predictions with the regression model

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

plot(insurance\(pred, insurance\)expenses) abline(a = 0, b = 1, col = “red”, lwd = 3, lty = 2)

predict(ins_model2, data.frame(age = 30, age2 = 30^2, children = 2, bmi = 30, sex = “male”, bmi30 = 1, smoker = “no”, region = “northeast”)) predict(ins_model2, data.frame(age = 30, age2 = 30^2, children = 2, bmi = 30, sex = “female”, bmi30 = 1, smoker = “no”, region = “northeast”))

set up the data

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)

compute the SDR

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

compare the SDR for each split

sdr_a

sdr_b

wine <- read.csv(“whitewines.csv”)

examine the wine data

str(wine) # the distribution of quality ratings hist(wine$quality) # summary statistics of the wine data summary(wine)

wine_train <- wine[1:3750, ] wine_test <- wine[3751:4898, ]

regression tree using rpart

library(rpart) m.rpart <- rpart(quality ~ ., data = wine_train) # get basic information about the tree m.rpart # get more detailed information about the tree summary(m.rpart)

install.packages(“rpart.plot”)

use the rpart.plot package to create a visualization

library(rpart.plot) # a basic decision tree diagram rpart.plot(m.rpart, digits = 3) # a few adjustments to the diagram rpart.plot(m.rpart, digits = 4, fallen.leaves = TRUE, type = 3, extra = 101)

generate predictions for the testing dataset

p.rpart <- predict(m.rpart, wine_test)

compare the distribution of predicted values vs. actual values

summary(p.rpart) summary(wine_test\(quality) # compare the correlation cor(p.rpart, wine_test\)quality)

function to calculate the mean absolute error

MAE <- function(actual, predicted) { mean(abs(actual - predicted))
}

MAE(p.rpart, wine_test\(quality) # mean absolute error between actual values and mean value mean(wine_train\)quality) # result = 5.87

MAE(5.87, wine_test$quality)

install.packages(“plyr”) install.packages(“Cubist”) # train a Cubist Model Tree library(Cubist) m.cubist <- cubist(x = wine_train[-12], y = wine_train$quality) # display basic information about the model tree m.cubist

display the tree itself

summary(m.cubist)

generate predictions for the model

p.cubist <- predict(m.cubist, wine_test)

summary statistics about the predictions

summary(p.cubist)

correlation between the predicted and true values

cor(p.cubist, wine_test$quality)

mean absolute error of predicted and true values

(uses a custom function defined above)

MAE(wine_test$quality, p.cubist)