This week:

• Wednesday (May 4): Lab 9 due in GitHub (by 5pm); Homework 9 due in D2L (by 11pm)
• Friday (May 6) by 11pm:
  • Take-home final exam (optional) released at 8:00am in D2L
  • Post comments on at least two of your classmate’s RShiny app D2L discussion posts
  • Complete D2L team evaluation survey

Finals week:

• Take-home final exam (optional) due in Gradescope by 11:00pm Monday May 9
• In-class final exam Tuesday May 10 at 12:00-1:50pm in our usual classroom

Your final exam grade will be the higher of the following scores:

1. In-class final exam score.
2. Weighted average of your take-home and in-class final exam scores, weighted as 40% in-class, 60% take-home.

Here are a few questions to consider:

• What does statistical learning mean to you?
• Is statistical learning different from statistics as a whole?
• What about terms like: data science, data mining, data analytics, machine learning, predictive analytics, how are these different from statistics and statistical learning?

Statistical learning refers to a set of tools for modeling and understanding complex datasets. It is a recently developed area in statistics and blends with parallel developments in computer science and, in particular, machine learning. The field encompasses many methods such as the lasso and sparse regression, classification and regression trees, and boosting and support vector machines.

Courtesy of An Introduction to Statistical Learning: with Applications in R, by Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani. Note: a free e-version of this textbook can be obtained through the MSU Library.

• Features - what statisticians typically call “variables” or “predictors”
• Training vs testing data - Use training data to fit the model; use testing data to evaluate the model (e.g., goodness-of-fit, predictive)
• Supervised learning - Outcome variable values (“labels”) are known for your training data, e.g., linear regression, logistic regression, classification
• Unsupervised learning - No clear outcome variable or outcome variable values unknown for training data, e.g., clustering, dimension reduction

“The basic goal of predictive modeling is to find a function that accurately describes how different measured explanatory variables can be combined to make a prediction about a response variable.”

• Baumer, Kaplan and Horton (Section 10.1)

Recall the Seattle housing data set, how would you:

• Build a model to predict housing prices in King County
• Determine if your model was good or useful?

A loss function is a principled way to compare a set of predictive models.

Squared Error: \[ (Price_{pred} - Price_{actual}) ^ 2\]

Zero-One Loss (binary setting): \[ f(x)= \begin{cases} 1,& \text{if } y_{pred} \neq y_{actual}\\ 0, & y_{pred} = y_{actual} \end{cases} \]

Suppose we fit a model using all of the Seattle housing data — can that model be used to predict prices for homes in that data set?

We cannot assess the predictive performance by fitting a model to data and then evaluating the model using the same data.

There are two common options to give valid estimates of model performance:

• Test / Training approach. Generally 70% of the data is used to fit the model and the other 30% is held out for prediction.

• Cross-Validation. Cross validation breaks your data into k groups, or folds. Then a model is fit on the data on the k-1 groups and then used to make predictions on data in the held out k\(^{th}\) group. This process continues until all groups have been held out once.

set.seed(11142017)
num.houses <- nrow(Seattle)
Seattle$zipcode <- as.factor(Seattle$zipcode)
test.ids <- base::sample(1:num.houses, size=round(num.houses*.3))
test.set <- Seattle[test.ids,]
train.set <- Seattle[(1:num.houses)[!(1:num.houses) %in% 
  test.ids],]
dim(Seattle)

## [1] 869  14

dim(test.set)

## [1] 261  14

dim(train.set)

## [1] 608  14

lm.1 <- lm(price ~ bedrooms + bathrooms + sqft_living + 
                   zipcode + waterfront, 
           data=train.set)
summary(lm.1)

mad.lm1 <- mean(abs(test.set$price - predict(lm.1, test.set)))

The mean absolute deviation in housing price predictions using the linear model is $\(1.80149\times 10^{5}\)

Now include squared terms for square foot of living space, too.

train.set$sqft_living2 <- train.set$sqft_living^2
test.set$sqft_living2 <- test.set$sqft_living^2

lm.2 <- lm(price ~ bedrooms + bathrooms + sqft_living + sqft_living2 + 
                   zipcode + waterfront, 
           data=train.set)
summary(lm.2)

mad.lm2 <- mean(abs(test.set$price - predict(lm.2,test.set)))

Including this squared term lowers our predictive error from $\(1.80149\times 10^{5}\) in the first case to $\(1.38583\times 10^{5}\).

rmse.tree1 <- sqrt(mean((test.set$price - predict(tree1,test.set))^2)) 
mad.tree1 <- mean(abs(test.set$price - predict(tree1,test.set)))
mad.tree1

## [1] 205934.8

mad.lm1

## [1] 180149.3

mad.lm2

## [1] 138583

The predictive error for this tree, $\(2.05935\times 10^{5}\) is similar to the first linear model $\(1.80149\times 10^{5}\) and not quite as good as our second linear model $\(1.38583\times 10^{5}\).

Ensemble methods combine a large set of predictive models into a single framework. One example is a random forest, which combines a large number of trees.

While these methods are very effective in a predictive setting, it is often difficult to directly assess the impact of particular variables in the model.

One specific kind of ensemble method is known as a random forest, which combines several decision trees.

rf1 <- randomForest(price ~ ., data = Seattle)

mad.rf <- mean(abs(test.set$price - predict(rf1,test.set)))

The prediction error for the random forest is substantially better than the other models we have identified $\(5.1452\times 10^{4}\).

bikes <- read.csv('http://www.math.montana.edu/ahoegh/teaching/stat408/datasets/Bike.csv')
set.seed(11142017)
num.obs <- nrow(bikes)
test.ids <- base::sample(1:num.obs, size=round(num.obs*.3))
test.bikes <- bikes[test.ids,]
train.bikes <- bikes[(1:num.obs)[!(1:num.obs) %in% 
  test.ids],]
dim(bikes)

## [1] 10886    12

dim(test.bikes)

## [1] 3266   12

dim(train.bikes)

## [1] 7620   12

lm.bikes <- lm(count ~ holiday + atemp,
               data=train.bikes)
lm.mad <- mean(abs(test.bikes$count - predict(lm.bikes,test.bikes)))

Create another predictive model and compare the results to the MAD of the linear model above (\(127\)). However, don’t use both casual and registered in your model as those two will sum to the total count.

rf.bikes <- randomForest(count ~ holiday + atemp +
     humidity + season + workingday + weather,
                data=train.bikes)
tree.mad <- mean(abs(test.bikes$count - predict(rf.bikes,test.bikes)))

The random forest has a prediction error of \(108\).

kable(cbind(pred.points,
    round(predict(tree1,as.data.frame(pred.points))[,2],3)),
    col.names=c('x','y','Prob[Val = 1]'))

tree1 <- rpart(labels ~., data=supervised, method = 'class')
plot(tree1, margin = 1)
text(tree1, fancy = TRUE)

See if you can improve the classification error from the model below.

glm.titanic <- glm(Survived ~ Age, data=train.titanic, family='binomial')
Class.Error <- mean(test.titanic$Survived != round(predict(glm.titanic, test.titanic, type='response')))

The logistic regression model only using age is wrong \(44\)% of the time.

km <- kmeans(combined, 3)
plot(combined,type='n',axes=F, xlab='',ylab='')
box()
points(combined,pch=as.character(km$cluster),
       col=c(rep('dodgerblue',25),
             rep('forestgreen',25),
             rep('firebrick4',50)))
draw.circle(.31,-0.1,.335, border='dodgerblue')
draw.circle(.79,.65,.3, border='firebrick4')
draw.circle(.14,1.05,.3, border='forestgreen')

hc <- hclust(dist(combined))
plot(hc, hang=-1)
plot(combined,type='n',axes=F, xlab='',ylab='')
box()
points(combined,pch=as.character(cutree(hc,4)),
       col=c(rep('dodgerblue',25),
             rep('forestgreen',25),
             rep('firebrick4',50)))

Given these plots that we have seen, how do we choose the number of clusters?

wss <- rep(0,15)
for (i in 1:15) {
  wss[i] <- sum(kmeans(combined,centers=i)$withinss)
}
plot(1:15, wss, type="b", xlab="Number of Clusters", 
     ylab="Within groups sum of squares")

Use the dataset create below for the following questions.

zoo.data <- read.csv('http://www.math.montana.edu/ahoegh/teaching/stat408/datasets/ZooClean.csv')
rownames(zoo.data) <- zoo.data[,1]
zoo.data <- zoo.data[,-1]

• Use multidimensional scaling to visualize the data in two dimensions. What are two animals that are very similar and two that are very different?

• Create a hierachical clustering object for this dataset. Why are a leopard and raccoon clustered together for any cluster size?

• Now add colors corresponding to four different clusters to your MDS plot. Interpret what each of the four clusters correspond to.