From 7bf0aa01c3e92beee3b307e24679b7ee48f008cf Mon Sep 17 00:00:00 2001
From: Paul Adamson <paul.adamson.01@gmail.com>
Date: Wed, 7 Dec 2016 19:55:00 -0500
Subject: [PATCH 1/4] Initial commit of R code for chapter 2.

Initial commit of R code for chapter 2, which includes examples of the
following:

  - model.matrix command for convert categorical features to numerical
  binary features

  - feature extraction with packages dplyr and tidyr

  - feature normalization

  - plotting with packages vcd and ggplot2 (including facets)
---
 .gitignore     |   4 +
 R/Chapter2.Rmd | 232 +++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 236 insertions(+)
 create mode 100644 R/Chapter2.Rmd

diff --git a/.gitignore b/.gitignore
index 348fc88..51e4953 100644
--- a/.gitignore
+++ b/.gitignore
@@ -66,3 +66,7 @@ target/
 
 # Figures
 figures/
+.Rproj.user
+
+# knitr html output
+R/*.html
diff --git a/R/Chapter2.Rmd b/R/Chapter2.Rmd
new file mode 100644
index 0000000..e82226b
--- /dev/null
+++ b/R/Chapter2.Rmd
@@ -0,0 +1,232 @@
+---
+title: "Real-World Machine Learning"
+subtitle: "Chapter 2"
+author: "Paul Adamson"
+date: "December 7, 2016"
+output: html_document
+---
+
+This file contains R code to accompany Chapter 2 of the book 
+["Real-World Machine Learning"](https://www.manning.com/books/real-world-machine-learning),
+by Henrik Brink, Joseph W. Richards, and Mark Fetherolf.  The code was contributed by
+[Paul Adamson](http://github.com/padamson). 
+
+*NOTE: working directory should be set to this file's location.*
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+library(plyr)
+library(dplyr)
+library(tidyr)
+library(ggplot2)
+library(vcd)
+library(sm)
+```
+
+## Listing 2.1 Convert categorical features to numerical binary features
+
+Creating dummy variables in R is extremely easy with the `model.matrix` function.
+In the below code, when the `personData` dataframe is created,
+the `maritalstatus` variable is of type `factor` with levels "single" and
+"married" by default. 
+
+```{r listing2.1}
+personData <- data.frame(person = 1:2, 
+                         name = c("Jane Doe", "John Smith"),
+                         age = c(24, 41),
+                         income = c(81200, 121000),
+                         maritalstatus = c("single","married"))
+personDataNew <- data.frame(personData[,1:4], 
+                            model.matrix(~ maritalstatus - 1, 
+                                         data = personData)) 
+str(personData)
+```
+
+In the call to `model.matrix`, the −1 in the model formula 
+ensures that we create a dummy variable for each of the two marital statuses 
+(technically, it suppresses the creation of an intercept). 
+```{r}
+model.matrix(~ maritalstatus - 1, 
+             data = personData)
+```
+The matrix of dummy variables is then joined to the original dataframe (minus
+the maritalstatus column) with another call to `data.frame`.
+
+## Listing 2.2 Simple feature extraction on Titanic cabins
+
+The packages [`dplyr`](https://cran.rstudio.com/web/packages/dplyr) and 
+[`tidyr`](https://cran.r-project.org/web/packages/tidyr) are excellent 
+for tidying and preprocessing data, 
+including creating new features from existing ones. (Note: `plyr` will
+be used later, but we must load it prior to `dplyr`.)
+
+```{r listing2.2}
+titanic <- read.csv("../data/titanic.csv", 
+                    colClasses = c(
+                      Survived = "factor",
+                      Name = "character",
+                      Ticket = "character",
+                      Cabin = "character"))
+
+titanic$Survived <- revalue(titanic$Survived, c("0"="no", "1"="yes"))
+
+titanicNew <- titanic %>%
+  separate(Cabin, into = "firstCabin", sep = " ", extra = "drop", remove = FALSE) %>%
+  separate(firstCabin, into = c("cabinChar", "cabinNum"), sep = 1) %>%
+  rowwise() %>%
+  mutate(numCabins = length(unlist(strsplit(Cabin, " "))))
+
+str(titanicNew)
+```
+
+In Listing 2.2, `read.csv` is used to read in the comma separated value
+(csv) data file. The `colClasses` argument is used to specify the correct
+class for some features. Then, the `revalue` function changes the levels
+of the `Survived` factor variable so that '0' indicates 'no' and '1'
+indicates 'yes'. The `titanicNew` dataframe is then created by piping
+together `separate` from `tidyr` and `mutate` from `dplyr`. `separate`
+does exactly what its name implies: it separates a single character column
+into multiple columns. `mutate` is used to add a new feature often (as
+in this case) by acting on values of another feature.
+
+## Listing 2.3 Feature normalization
+
+The below code will normalize a feature using the "min-max" method. As
+an example, the `Age` feature of the `titanic` dataframe is normalized
+and a histogram of the new normalized feature is plotted with
+[`ggplot2`](http://docs.ggplot2.org/current/index.html#). 
+
+```{r normalize}
+normalizeFeature <- function(data, fMin=-1.0, fMax=1.0){
+  dMin = min(na.omit(data))
+  dMax = max(na.omit(data))
+  factor = (fMax - fMin) / (dMax - dMin)
+  normalized = fMin + (data - dMin)*factor
+  normalized
+}
+
+titanic$AgeNormalized <- normalizeFeature(titanic$Age)
+ggplot(data=titanic, aes(AgeNormalized)) +
+  geom_histogram()
+```
+
+## Figure 2.12 Mosaic plot for Titanic data: Gender vs. survival
+
+The ["Visualizing Categorical Data" (`vcd`)](https://cran.r-project.org/web/packages/vcd/vcd.pdf) package 
+provides an excellent set of functions for exploring categorical data,
+including mosaic plots.
+
+```{r figure2.12}
+mosaic(
+  ~ Sex + Survived,
+  data = titanic, 
+  main = "Mosaic plot for Titanic data: Gender vs. survival",
+  shade = TRUE,
+  split_vertical = TRUE,
+  labeling_args = list(
+    set_varnames = c(
+      Survived = "Survived?")))
+```
+
+## Figure 2.13 Mosaic plot for Titanic data: Passenger class vs. survival
+
+```{r figure2.13}
+mosaic(
+  ~ Pclass + Survived,
+  data = titanic, 
+  main = "Mosaic plot for Titanic data: Passenger Class vs. survival",
+  shade = TRUE,
+  split_vertical = TRUE,
+  labeling_args = list(
+    set_varnames = c(
+      Pclass = "Passenger class",
+      Survived = "Survived?")))
+```
+
+## Figure 2.14 Box plot for Titanic data: Passenger age vs. survival
+
+The `boxplot` function is provided as part of the standard `graphics` package
+in R. `ggplot2` provides a much nicer version.
+```{r figure2.14a}
+boxplot(Age ~ Survived, 
+        data = titanic,
+        xlab = "Survived?",
+        ylab = "Age\n(years)",
+        las = 1)
+```
+
+```{r figure2.14b}
+ggplot(titanic, aes(Survived, Age)) + 
+  geom_boxplot() +
+  xlab("Survived?") +
+  ylab("Age\n(years)")
+```
+
+## Figure 2.15 Box plots for Titanic data: Passenger fare versus survival
+
+Plots can be combined in rows and columns using the `mfrow` graphical
+parameter set via the `par` function. (Try `help(par)` to learn more.) 
+```{r figure2.15}
+par(mfrow=c(1,2))
+par(mai=c(1,1,.1,.1), las = 1)
+boxplot(Fare ~ Survived, 
+        data = titanic,
+        xlab = "Survived?",
+        ylab = "Fare Amount")
+boxplot(Fare**(1/2) ~ Survived, 
+        data = titanic,
+        xlab = "Survived?",
+        ylab = "sqr\n(fare amount)")
+```
+
+## Figure 2.16 Density plot for MPG data, by region
+
+The `sm.density.compare` function from the `sm` package is useful for
+comparing a set of univariate density estimates.  The first argument is
+a vector of data, and the second argument is a vector of group labels that
+correspond to each value. The `colfill` variable and `legend` function
+are used to place a legend on the plot generated by `sm.density.compare`.
+
+```{r figure2.16}
+par(mfrow=c(1,1))
+auto <- read.csv("../data/auto-mpg.csv",
+                 colClasses = c(origin = "factor"))
+
+auto$origin <- revalue(auto$origin, 
+                       c("1\t"="USA", "2\t"="Europe", "3\t"="Asia"))
+
+sm.density.compare(auto$mpg, auto$origin,
+                   xlab="Miles per gallon",
+                   ylab="Density")
+title(main="Density plot for MPG data, by region")
+
+colfill<-c(2:(2+length(levels(auto$origin)))) 
+legend(x="topleft", 
+       inset=0.05, 
+       text.width=5, 
+       levels(auto$origin), 
+       fill=colfill)
+```
+
+## Figure 2.17 Scatterplots for MPG data
+
+It doesn't get much simpler than the `plot` function in R.
+```{r figure2.17}
+par(mfrow=c(1,2), mai=c(1,1,.1,.1), las = 1)
+plot(auto$weight, auto$mpg,
+     xlab = "Vehicle weight",
+     ylab = "Miles per\ngallon")
+plot(auto$modelyear, auto$mpg,
+     xlab = "Model year",
+     ylab = "Miles per\ngallon")
+
+```
+
+Although not in the book, everyone should become familiar with facets via 
+[`ggplot2`](http://docs.ggplot2.org/current/index.html#).
+
+```{r figure2.17_ggplot2}
+p <- ggplot(auto, aes(mpg, weight)) + 
+  geom_point()
+p + facet_grid(. ~ origin)
+```
\ No newline at end of file

From 5e71a1a015314abe1e3e03a6596b5b3feca1aada Mon Sep 17 00:00:00 2001
From: Paul Adamson <paul.adamson.01@gmail.com>
Date: Wed, 7 Dec 2016 20:25:10 -0500
Subject: [PATCH 2/4] Initial commit of chapter 3 R code

Chapter 3 R code includes

  - knn model of MNIST handwriting data

  - nonlinear random forest model of MPG data
---
 .gitignore     |   4 +-
 R/Chapter3.Rmd | 247 +++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 250 insertions(+), 1 deletion(-)
 create mode 100644 R/Chapter3.Rmd

diff --git a/.gitignore b/.gitignore
index 51e4953..b347ff6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -68,5 +68,7 @@ target/
 figures/
 .Rproj.user
 
-# knitr html output
+# knitr output
 R/*.html
+R/*_cache
+R/*_files
diff --git a/R/Chapter3.Rmd b/R/Chapter3.Rmd
new file mode 100644
index 0000000..f88f884
--- /dev/null
+++ b/R/Chapter3.Rmd
@@ -0,0 +1,247 @@
+---
+title: "Real-World Machine Learning" 
+subtitle: "Chapter 3"
+author: "Paul Adamson"
+date: "December 7, 2016"
+output: html_document
+---
+
+This file contains R code to accompany Chapter 3 of the book 
+["Real-World Machine Learning"](https://www.manning.com/books/real-world-machine-learning),
+by Henrik Brink, Joseph W. Richards, and Mark Fetherolf.  The code was contributed by
+[Paul Adamson](http://github.com/padamson). 
+
+*NOTE: working directory should be set to this file's location.*
+
+```{r setup, include=FALSE}
+library(knitr)
+knitr::opts_chunk$set(echo = TRUE)
+library(plyr)
+library(dplyr)
+library(vcd)
+library(AppliedPredictiveModeling)
+library(caret)
+library(ellipse)
+library(kknn)
+library(gridExtra)
+library(grid)
+library(randomForest)
+set.seed(3456)
+.pardefault <- par()
+```
+
+## Figure 3.4 A subset of the Titanic Passengers dataset
+
+We are going to be interested in predicting survival, so it is useful to specify 
+the `Survived` variable to be of type `factor`. For visualizing the data, 
+it is also useful to use the `revalue` function to specify the `no` and `yes`
+levels for the `factor` variable. The `kable` function is built into the `knitr`
+package.
+
+```{r figure3.4, cache=TRUE}
+titanic <- read.csv("../data/titanic.csv", 
+                    colClasses = c(
+                      Survived = "factor",
+                      Name = "character",
+                      Ticket = "character",
+                      Cabin = "character"))
+titanic$Survived <- revalue(titanic$Survived, c("0"="no", "1"="yes"))
+kable(head(titanic, 6), digits=2)
+```
+
+## Figure 3.5 Mosaic plot for Titanic data: Gender vs. survival
+
+The ["Visualizing Categorical Data" package  (`vcd`)](https://cran.r-project.org/web/packages/vcd/vcd.pdf)
+provides an excellent set of functions for exploring categorical data,
+including mosaic plots.
+
+```{r figure3_5, cache=TRUE, dependson="figure3.4"}
+mosaic(
+  ~ Sex + Survived,
+  data = titanic, 
+  main = "Mosaic plot for Titanic data: Gender vs. survival",
+  shade = TRUE,
+  split_vertical = TRUE,
+  labeling_args = list(
+    set_varnames = c(
+      Survived = "Survived?",
+      Sex = "Gender")))
+```
+
+## Figure 3.6 Processed data
+
+First, we get rid of the variables that we do not want in our model.
+(`Cabin` might actually be useful, but it's not used here.)
+Then we use `is.na` to set missing age values to -1.
+The `mutate` and `select` functions make it easy to take square root of 
+the `Fare` variable and then drop it from the dataset.
+We then drop rows with missing `Embarked` data and remove the unused level 
+`""`. 
+Finally, we convert `factor` variables to dummy variables using the 
+`dummyVars` function in the `caret` package.
+To avoid perfect collinearity (a.k.a. the dummy variable trap), we set
+the `fullRank` parameter to `TRUE`.  `Survived.yes` is then converted back
+to a `factor` variable.
+
+```{r figure3.6, cache=TRUE, dependson="figure3.4"}
+titanicTidy <- subset(titanic, select = -c(PassengerId, Name, Ticket, Cabin))
+
+titanicTidy$Age[is.na(titanicTidy$Age)] <- -1
+
+titanicTidy <- titanicTidy %>%
+  mutate(sqrtFare = sqrt(Fare)) %>%
+  select(-Fare)
+
+titanicTidy <- titanicTidy %>%
+  filter(!(Embarked=="")) %>%
+  droplevels
+
+dummies <- dummyVars(" ~ .", data = titanicTidy, fullRank = TRUE)
+titanicTidyNumeric <- data.frame(predict(dummies, newdata = titanicTidy))
+
+titanicTidyNumeric$Survived.yes <- factor(titanicTidyNumeric$Survived.yes)
+kable(head(titanicTidyNumeric))
+```
+
+## Figure 3.10 Four randomly chosen handwritten digits from the MNIST database
+
+Thanks to [Longhow Lam](https://longhowlam.wordpress.com/2015/11/25/a-little-h2o-deeplearning-experiment-on-the-mnist-data-set/)
+for posting the code used in the `displayMnistSamples` function that display's 
+digits from the MNIST dataset.
+
+```{r figure3.10, cache=TRUE,fig.height=2}
+mnist <- read.csv("../data/mnist_small.csv",
+                  colClasses = c(label = "factor"))
+displayMnistSamples <- function(x) {
+  for(i in x){
+  y = as.matrix(mnist[i, 2:785])
+  dim(y) = c(28, 28)
+  image( y[,nrow(y):1], axes = FALSE, col = gray(0:255 / 255))
+  text( 0.2, 0, mnist[i,1], cex = 3, col = 2, pos = c(3,4))
+  }
+}
+par( mfrow = c(1,4), mai = c(0,0,0,0.1))
+displayMnistSamples(sample(1:length(mnist),4))
+```
+
+## Figure 3.11 Table of predicted probabilities from a k-nearest neighbors classifier, as applied to the MNIST dataset
+
+```{r figure3.11, cache=TRUE, dependson="figure3.10"}
+trainIndex <- createDataPartition(mnist$label, p = .8, 
+                                  list = FALSE, 
+                                  times = 1)
+mnistTrain <- mnist[ trainIndex,]
+mnistTest  <- mnist[-trainIndex,]
+
+mnist.kknn <- kknn(label~., mnistTrain, mnistTest, distance = 1,
+                   kernel = "triangular")
+
+confusionMatrix(fitted(mnist.kknn),mnistTest$label)
+
+rowsInTable <- 1:10
+prob <- as.data.frame(mnist.kknn$prob[rowsInTable,])
+mnistResultsDF <- data.frame(mnistTest$label[rowsInTable],
+                             mnist.kknn$fit[rowsInTable],
+                             prob)
+
+
+kable(mnistResultsDF, digits=2,
+      col.names=c("actual","fit",0:9))
+```
+
+## Figure 3.13 Small subset of the Auto MPG data
+
+```{r figure3.13, cache=TRUE}
+auto <- read.csv("../data/auto-mpg.csv",
+                 colClasses = c(
+                      origin = "factor"))
+
+auto$origin <- revalue(auto$origin, 
+                       c("1\t"="USA", "2\t"="Europe", "3\t"="Asia"))
+
+kable(head(auto,5))
+```
+
+## Figure 3.14 Scatter plots of Vehicle Weight and Model Year versus MPG
+
+```{r figure3.14, cache=TRUE, dependson="figure3.13", warning=FALSE}
+par(.pardefault)
+p1<-ggplot(auto, aes(weight, mpg)) + 
+  geom_point() +
+  labs(y = "Miles per gallon",
+       x = "Vehicle weight")
+p2<-ggplot(auto, aes(modelyear, mpg)) + 
+  geom_point() +
+  theme(axis.title.y=element_blank(),
+        axis.text.y=element_blank(),
+        axis.ticks.y=element_blank()) +
+  labs(x = "Model year")
+grid.arrange(p1,p2,ncol=2, 
+             top=textGrob("Scatterplots for MPG data",
+                          gp=gpar(fontsize=14,font=8)))
+```
+
+## Figure 3.15 The Auto MPG data after expanding the categorical Origin column
+
+Note that the row numbering differs between python and R by 1
+(python starts row numbering at 0, and R starts at 1).
+
+```{r figure3.15, cache=TRUE, dependson="figure3.13"}
+dummies <- dummyVars(" ~ .", data = auto, fullRank = TRUE)
+autoNumeric <- data.frame(predict(dummies, newdata = auto))
+
+kable(tail(autoNumeric,5))
+```
+
+## Figure 3.16 Comparing MPG predictions on a held-out testing set to actual values
+
+```{r figure3.16, cache=TRUE, dependson="figure3.15"}
+trainIndex <- createDataPartition(autoNumeric$mpg, p = .8, 
+                                  list = FALSE, 
+                                  times = 1)
+
+autoTrain <- autoNumeric[ trainIndex,]
+autoTest  <- autoNumeric[-trainIndex,]
+
+lmFit <- train(mpg ~ ., data = autoTrain,
+               method = "lm")
+lmPred <- predict(lmFit, newdata = autoTest)
+
+rowsInTable <- 1:5
+kable(data.frame("Origin.Europe" = autoTest$origin.Europe[rowsInTable],
+                 "Origin.Asia" = autoTest$origin.Asia[rowsInTable],
+                 "MPG" = autoTest$mpg[rowsInTable],
+                 "Predicted MPG" = lmPred[rowsInTable]))
+```
+
+## Figure 3.17 A scatter plot of the actual versus predicted values on the held-out test set. The diagonal line shows the perfect regressor. The closer all of the predictions are to this line, the better the model.
+
+```{r figure3.17, cache=TRUE, dependson="figure3.16"}
+ggplot(autoTest, aes(x=mpg, y=lmPred)) + 
+  geom_point() + 
+  geom_abline(slope = 1, intercept = 0) +
+  labs(x="MPG", y="Predicted MPG")
+```
+
+## Figure 3.18 Table of actual versus predicted MPG values for the nonlinear random forest regression model
+
+```{r figure3.18, cache=TRUE, dependson="figure3.16"}
+rfFit <- train(mpg ~ ., data = autoTrain,
+               method = "rf")
+rfPred <- predict(rfFit, newdata = autoTest)
+
+kable(data.frame("Origin.Europe" = autoTest$origin.Europe[rowsInTable],
+                 "Origin.Asia" = autoTest$origin.Asia[rowsInTable],
+                 "MPG" = autoTest$mpg[rowsInTable],
+                 "Predicted MPG" = rfPred[rowsInTable]))
+
+```
+
+## Figure 3.19 Comparison of MPG data versus predicted values for the nonlinear random forest regression model
+
+```{r figure3.19, cache=TRUE, dependson="figure3.18"}
+ggplot(autoTest, aes(x=mpg, y=rfPred)) + 
+  geom_point() + 
+  geom_abline(slope = 1, intercept = 0) +
+  labs(x="MPG", y="Predicted MPG")
+```

From f2d74a5dc0c64231fc8c4858863801c2b408c785 Mon Sep 17 00:00:00 2001
From: Paul Adamson <paul.adamson.01@gmail.com>
Date: Sun, 25 Dec 2016 21:56:28 -0500
Subject: [PATCH 3/4] Initial commit of Chapter4.Rmd

---
 .gitignore     |   1 +
 R/Chapter4.Rmd | 470 +++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 471 insertions(+)
 create mode 100644 R/Chapter4.Rmd

diff --git a/.gitignore b/.gitignore
index b347ff6..cfe354c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -72,3 +72,4 @@ figures/
 R/*.html
 R/*_cache
 R/*_files
+R/.Rhistory
diff --git a/R/Chapter4.Rmd b/R/Chapter4.Rmd
new file mode 100644
index 0000000..04d2db5
--- /dev/null
+++ b/R/Chapter4.Rmd
@@ -0,0 +1,470 @@
+---
+title: "Real-World Machine Learning"
+subtitle: "Chapter 4"
+author: "Paul Adamson"
+date: "December 25, 2016"
+output: html_document
+---
+
+This notebook contains R code to accompany Chapter 4 of the book 
+["Real-World Machine Learning"](https://www.manning.com/books/real-world-machine-learning),
+by  Henrik Brink, Joseph W. Richards, and Mark Fetherolf.  
+The code was contributed by [Paul Adamson](http://github.com/padamson). 
+
+*NOTE: working directory should be set to this file's location.*
+
+*NOTE: depending on your machine, you may need to adjust the parameter in the call to `makeCluster()` in the code to generate Figure 4.19*
+
+```{r setup, include=FALSE}
+set.seed(1234)
+library(knitr)
+knitr::opts_chunk$set(echo = TRUE)
+library(plyr)
+library(dplyr)
+library(AppliedPredictiveModeling)
+library(caret)
+library(gbm)
+
+library(doParallel)
+library(randomForest)
+
+library(RColorBrewer)
+
+.pardefault <- par()
+```
+
+## Figure 4.9 The first five rows of the Titanic Passengers dataset 
+
+As in Chapter 3, we are going to be interested in predicting survival, so again,
+it is useful to specify 
+the `Survived` variable to be of type `factor`. For visualizing the data, 
+it is also useful to use the `revalue` function to specify the `no` and `yes`
+levels for the `factor` variable. The `kable` function is built into the `knitr`
+package.
+
+```{r figure4.9, cache=TRUE}
+titanic <- read.csv("../data/titanic.csv", 
+                    colClasses = c(
+                      Survived = "factor",
+                      Name = "character",
+                      Ticket = "character",
+                      Cabin = "character"))
+titanic$Survived <- revalue(titanic$Survived, c("0"="no", "1"="yes"))
+kable(head(titanic, 5), digits=2)
+```
+
+## Figure 4.10 Splitting the full dataset into training and testing sets
+
+Here, we follow the same process used for Figure 3.6 to process the data and
+prepare it for our model. First, we get rid of the variables that we do not 
+want in our model.
+(`Cabin` might actually be useful, but it's not used here.)
+Then we use `is.na` to set missing age values to -1.
+The `mutate` and `select` functions make it easy to take square root of 
+the `Fare` variable and then drop it from the dataset.
+We then drop rows with missing `Embarked` data and remove the unused level 
+`""`. 
+Finally, we convert `factor` variables to dummy variables using the 
+`dummyVars` function in the `caret` package.
+To avoid perfect collinearity (a.k.a. the dummy variable trap), we set
+the `fullRank` parameter to `TRUE`.  `Survived.yes` is then converted back
+to a `factor` variable, and its levels are changed back to 'yes' and 'no' again
+(otherwise the `train` function in `caret` will complain later about
+invalid R variable names).
+
+We then make a 80%/20% train/test split 
+using the `Survived` factor
+variable in the `createDataPartition` function to preserve the
+overall class distribution of the data.
+
+```{r figure4.10, cache=TRUE, dependson="figure4.9"}
+titanicTidy <- subset(titanic, select = -c(PassengerId, Name, Ticket, Cabin))
+
+titanicTidy$Age[is.na(titanicTidy$Age)] <- -1
+
+titanicTidy <- titanicTidy %>%
+  mutate(sqrtFare = sqrt(Fare)) %>%
+  select(-Fare)
+
+titanicTidy <- titanicTidy %>%
+  filter(!(Embarked=="")) %>%
+  droplevels
+
+dummies <- dummyVars(" ~ .", data = titanicTidy, fullRank = TRUE)
+titanicTidyNumeric <- data.frame(predict(dummies, newdata = titanicTidy))
+
+titanicTidyNumeric$Survived.yes <- factor(titanicTidyNumeric$Survived.yes)
+titanicTidyNumeric$Survived.yes <- 
+  revalue(titanicTidyNumeric$Survived.yes, c("0"="no", "1"="yes"))
+
+trainIndex <- createDataPartition(titanicTidyNumeric$Survived.yes, p = .8, 
+                                  list = FALSE, 
+                                  times = 1)
+
+titanicTrain <- titanicTidyNumeric[ trainIndex,]
+titanicTest  <- titanicTidyNumeric[-trainIndex,]
+
+kable(head(titanicTidyNumeric, 8), digits=2, caption = "Full dataset")
+kable(head(titanicTrain, 5), digits=2, 
+      caption = "Training set: used only for building the model")
+kable(head(titanicTest, 3), digits=2,
+      caption = "Testing set: used only for evaluating model")
+```
+
+## Figure 4.11 Comparing the testing set predictions with the actual values gives you the accuracy of the model.
+
+```{r figure4.11, eval=TRUE, cache=TRUE, dependson="figure4.10"}
+objControl <- trainControl(method='cv', number=3, 
+                           returnResamp='none', 
+                           summaryFunction = twoClassSummary, 
+                           classProbs = TRUE)
+titanic.gbm <- train(Survived.yes~.,data=titanicTrain, method='gbm',
+                     trControl=objControl,  
+                     metric = "ROC",
+                     preProc = c("center", "scale"))
+
+titanic.gbm.pred <- predict(titanic.gbm,newdata=titanicTest)
+
+figure4.11.df <- data.frame(head(titanicTest$Survived.yes),
+                            head(titanic.gbm.pred))
+kable(figure4.11.df,
+      col.names = c("Test set labels", "Predictions"))
+```
+
+## Figure 4.13 Organizing the class-wise accuracy into a confusion matrix
+```{r figure4.13, eval=TRUE, cache=TRUE, dependson="figure4.11"}
+titanic.gbm.cm <- confusionMatrix(titanic.gbm.pred,titanicTest$Survived.yes)
+titanic.gbm.cm
+```
+## Figure 4.15 A subset of probabilistic predictions from the Titanic test set. After sorting the full table by decreasing survival probability, you can set a threshold and consider all rows above this threshold as survived. Note that the indices are maintained so you know which original row the instance refers to.
+
+```{r figure4.15,eval=TRUE, cache=TRUE, dependson='figure4.11'}
+titanic.gbm.pred.prob <- predict(object=titanic.gbm, 
+                                 titanicTest, type='prob')
+kable(titanic.gbm.pred.prob[15:19,],col.names = c("Died","Survived"))
+titanic.ordered <- titanic.gbm.pred.prob[order(-titanic.gbm.pred.prob$yes),]
+kable(titanic.ordered[titanic.ordered$yes > 0.68 & titanic.ordered$yes < 0.73,])
+```
+
+## Listing 4.3 The ROC curve
+
+```{r listing4.3, eval=TRUE, cache=TRUE}
+# Returns the false-positive and true-positive rates at nPoints thresholds for 
+# the given true and predicted labels
+# trueLabels: 0=FALSE; 1=TRUE
+rocCurve <- function(trueLabels, predictedProbs, nPoints=100, posClass=1){
+  # Allocates the threshold and ROC lists
+  thr <- seq(0,1,length=nPoints)
+  tpr <- numeric(nPoints)
+  fpr <- numeric(nPoints)
+  
+  # Precalculates values for the positive and negative cases, used in the loop
+  pos <- trueLabels == posClass
+  neg <- !pos 
+  nPos <- sum(pos, na.rm=TRUE)
+  nNeg <- sum(neg, na.rm=TRUE)
+  
+  # For each threshold, calculates the rate of true and false positives
+  for (i in 1:length(thr)) {
+    t <- thr[i]
+    meetOrExceedThreshold <- predictedProbs >= t
+    tpr[i] <-  sum((meetOrExceedThreshold & pos), na.rm=TRUE) / nPos
+    fpr[i] <-  sum((meetOrExceedThreshold & neg), na.rm=TRUE) / nNeg
+  }
+  
+  # Create data frame without duplicated fpr's to return
+  duplicatedFPRs <- duplicated(fpr)
+  df <- data.frame(fpr=fpr[!duplicatedFPRs],tpr=tpr[!duplicatedFPRs],thr=thr[!duplicatedFPRs])
+  
+  return(df)
+}
+```
+
+## Figure 4.16 The ROC curve defined by calculating the confusion matrix and ROC metrics at 100 threshold points from 0 to 1. By convention, you plot the false-positive rate on the x-axis and the true-positive rate on the y-axis.
+
+```{r figure4.16, eval=TRUE, cache=TRUE, dependson=c("listing4.3, figure4.15")}
+df<-rocCurve(revalue(titanicTest$Survived.yes, c("no" = 0, "yes" = 1)),
+             titanic.gbm.pred.prob$yes)
+ggplot(df,aes(x=fpr,y=tpr)) +
+  geom_step(direction="vh") +
+  labs(x = "False-positive rate",
+       y = "True-positive rate")
+```
+
+## Listing 4.4 The area under the ROC curve
+
+```{r listing4.4, eval=TRUE, cache=TRUE}
+auc <- function(trueLabels, predictedProbs, nPoints=100, posClass=1){
+  auc <- 0
+  df <- rocCurve(trueLabels = trueLabels,
+                 predictedProbs = predictedProbs,
+                 nPoints = nPoints,
+                 posClass = posClass)
+  
+  for (i in 2:length(df$fpr)) {
+    auc <- auc + 0.5 * (df$fpr[i-1] - df$fpr[i]) * (df$tpr[i-1] + df$tpr[i])
+  }
+  
+  return(auc)
+}
+```
+## Bonus. Relative influence of variables in model.
+
+We can call the `summary` function on our model to get a data frame
+of variables (`var`) and their relative influence (`rel.inf`) in the model.
+Before plotting the data, we reorder by the `rel.inf` variable.
+
+```{r bonus.1, results="hide", eval=TRUE, cache=TRUE, dependson="figure4.11"}
+gbmSummary <- summary(titanic.gbm)
+```
+
+```{r bonus.2, eval=TRUE, cache=TRUE, dependson="bonus.1"}
+gbmSummary <- transform(gbmSummary,
+                        var = reorder(var,rel.inf))
+ggplot(data=gbmSummary, aes(var, rel.inf)) +
+  geom_bar(stat="identity") +
+  coord_flip() + 
+  labs(x="Relative Influence",
+       y="Variable")
+```
+
+## Figure 4.18 Handwritten digits in the MNIST dataset 
+
+Thanks to [Longhow Lam](https://longhowlam.wordpress.com/2015/11/25/a-little-h2o-deeplearning-experiment-on-the-mnist-data-set/)
+for posting the code used in the `displayMnistSamples` function that display's 
+digits from the MNIST dataset.
+
+```{r figure4.18, cache=TRUE,fig.height=5}
+mnist <- read.csv("../data/mnist_small.csv",
+                  colClasses = c(label = "factor"))
+displayMnistSamples <- function(x) {
+  for(i in x){
+  y = as.matrix(mnist[i, 2:785])
+  dim(y) = c(28, 28)
+  image( y[,nrow(y):1], axes = FALSE, col = gray(0:255 / 255))
+  text( 0.2, 0, mnist[i,1], cex = 3, col = 2, pos = c(3,4))
+  }
+}
+par( mfrow = c(4,5), mai = c(0,0,0.1,0.1))
+displayMnistSamples(sample(1:length(mnist),20))
+levels(mnist$label) <- make.names(levels(mnist$label))
+```
+
+## Figure 4.19 The confusion matrix for the 10-class MNIST handwritten digit classification problem
+
+```{r figure4.19, cache=TRUE, dependson="figure4.18", fig.height=4}
+trainIndex <- createDataPartition(mnist$label, p = .8, 
+                                  list = FALSE, 
+                                  times = 1)
+
+mnistTrain <- mnist[ trainIndex,]
+mnistTest  <- mnist[-trainIndex,]
+
+cl = makeCluster(4)
+registerDoParallel(cl)
+
+str(mnistTrain$label)
+myTrainingControl <- trainControl(savePredictions = TRUE, 
+                              classProbs = TRUE, 
+                              verboseIter = FALSE)
+
+mnist.rf <- train(label~., data=mnistTrain, method='rf',
+                  ntree = 50,
+                  trControl = myTrainingControl,
+                  probability=TRUE, allowParallel=TRUE) 
+
+mnist.rf.pred = predict(mnist.rf,newdata=mnistTest)
+
+confusion.matrix <- confusionMatrix(mnist.rf.pred,mnistTest$label)
+
+confusionDF <- data.frame(confusion.matrix$table)
+
+confusionDF$Reference = with(confusionDF, 
+                             factor(Reference, levels = rev(levels(Reference))))
+
+jBuPuFun <- colorRampPalette(brewer.pal(n = 9, "BuPu"))
+paletteSize <- 256
+jBuPuPalette <- jBuPuFun(paletteSize)
+
+confusionPlot <- ggplot(
+  confusionDF, aes(x = Prediction, y = Reference, fill = Freq)) +
+  #theme(axis.text.x = element_text(angle = 0, hjust = 1, vjust = 0.5)) +
+  geom_tile() +
+  labs(x = "Predicted digit", y = "Actual digit") +
+  scale_fill_gradient2(
+    low = jBuPuPalette[1],
+    mid = jBuPuPalette[paletteSize/2],
+    high = jBuPuPalette[paletteSize],
+    midpoint = (max(confusionDF$Freq) + min(confusionDF$Freq)) / 2,
+    name = "") +
+  theme(legend.key.height = unit(2, "cm"))
+confusionPlot
+```
+
+## Figure 4.20 The ROC curves for each class of the MNIST 10-class classifier
+
+```{r figure4.20, cache=TRUE, dependson="figure4.19"}
+rocDF<-NULL
+aucDF<-NULL
+mnist.rf.pred.prob <- predict(object=mnist.rf,
+                              mnistTest, type='prob')
+for (i in 0:9){
+  digitLabel = paste0("X",i)
+  trueLabels <- as.numeric(mnistTest$label==digitLabel)
+  predictedProbs <- mnist.rf.pred.prob[[digitLabel]]
+  rocDF <- rbind(rocDF, 
+                 data.frame(rocCurve(trueLabels = trueLabels, predictedProbs = predictedProbs),digit=i))
+  aucDF <- rbind(aucDF, 
+                 data.frame(auc=auc(trueLabels = trueLabels, predictedProbs = predictedProbs),digit=i))
+}
+rocDF[,'digit'] <- as.factor(rocDF[,'digit'])
+labelVector <- c(paste0("Digit ",0:9,", AUC ",round(aucDF$auc,3)))  
+ggplot(rocDF[rocDF$fpr<0.2,],aes(x=fpr,y=tpr, linetype=digit, colour=digit)) +
+  geom_line() +
+  labs(x = "False-positive rate",
+       y = "True-positive rate") +
+  scale_linetype_discrete(name=NULL,
+                          labels=labelVector, 
+                          guide = guide_legend(keywidth = 3)) +
+  scale_colour_discrete(name=NULL,
+                        labels=labelVector,
+                        guide = guide_legend(keywidth = 3)) +
+  theme(legend.position=c(.8, .40)) 
+```
+
+
+
+## Figure 4.21 A subset of the Auto MPG dataset
+
+```{r figure4.21, cache=TRUE}
+auto <- read.csv("../data/auto-mpg.csv",
+                 colClasses = c(
+                      origin = "factor"))
+
+auto$origin <- revalue(auto$origin, 
+                       c("1\t"="USA", "2\t"="Europe", "3\t"="Asia"))
+
+kable(head(auto,5))
+```
+
+## Figure 4.22 Scatter plot of the predicted MPG versus actual values from the testing set. The diagonal line shows the optimal model.
+
+```{r figure4.22, cache=TRUE, dependson="figure4.21", fig.height=4}
+dummies <- dummyVars(" ~ .", data = auto, fullRank = TRUE)
+autoNumeric <- data.frame(predict(dummies, newdata = auto))
+
+trainIndex <- createDataPartition(autoNumeric$mpg, p = .8, 
+                                  list = FALSE, 
+                                  times = 1)
+
+autoTrain <- autoNumeric[ trainIndex,]
+autoTest  <- autoNumeric[-trainIndex,]
+
+ctrl <- trainControl(method = "repeatedcv", number = 10, savePredictions = TRUE)
+
+auto.glm <- train(mpg ~.,  
+                 data=autoTrain, 
+                 method="glm", 
+                 trControl = ctrl, 
+                 tuneLength = 5)
+
+auto.glm.Pred <- predict(auto.glm, newdata=autoTest)
+
+ggplot(autoTest, aes(x=mpg, y=auto.glm.Pred)) + 
+  geom_point() + 
+  geom_abline(slope = 1, intercept = 0) +
+  labs(x="MPG", y="Predicted MPG")
+```
+
+## Listing 4.5 The root-mean-square error
+
+```{r listing4.5, cache=TRUE}
+rmse <- function(trueValues, predictedValues){
+  return(sqrt(sum((trueValues - predictedValues)**2)/length(trueValues)))
+}
+rmse(autoTest$mpg, auto.glm.Pred)
+```
+## Listing 4.6 The R-squared calculation
+
+```{r listing4.6, cache=TRUE}
+r2 <- function(trueValues, predictedValues) {
+  meanTrueValues <- mean(trueValues)
+  return( 1.0 - (sum((trueValues - predictedValues)**2) / sum((trueValues - meanTrueValues)**2) ))
+}
+r2(autoTest$mpg, auto.glm.Pred)
+```
+
+## Figure 4.24 The residual plot from predictions on the MPG dataset.
+
+```{r figure4.24, cache=TRUE, dependson="figure4.22", fig.height=4}
+autoTest$residuals <- autoTest$mpg - auto.glm.Pred
+
+ggplot(data=autoTest, aes(x=mpg, y=residuals)) +
+  geom_point() + 
+  geom_abline(slope = 0, intercept = 0) +
+  labs(x="MPG", y="Residuals")
+```
+
+## (Substitute for Figure 4.25) Tile plots for Support Vector Machine (SVM) tuning parameter grid search
+
+The `getModelInfo` function in the `caret` package can be used to find
+the available model parameters. By setting `summaryFunction = twoClassSummary`
+in `trainControl`, we request `ROC` to be used to select the optimal model using 
+the largest value.
+
+```{r figure4.25sub1, cache=TRUE, dependson="figure4.10"}
+svmRadialModelInfo <- getModelInfo("svmRadial")
+svmRadialModelInfo$svmRadial$parameters
+
+svmGrid <- expand.grid(sigma = c(0.0001,0.001,0.01,0.1),
+                       C = seq(0.1,2.1,by=0.5)) 
+                         
+fitControl <- trainControl(method = "repeatedcv", 
+                           number = 10, 
+                           savePredictions = TRUE,
+                           summaryFunction=twoClassSummary, 
+                           classProbs=TRUE)
+
+svmFit <- train(Survived.yes ~ ., data = titanicTrain,
+                method = "svmRadial", 
+                trControl = fitControl, 
+                verbose = FALSE, 
+                tuneGrid = svmGrid)
+
+svmFitDF <- data.frame(sigma = as.factor(svmFit$results$sigma), 
+                       C = as.factor(svmFit$results$C),
+                       ROC = svmFit$results$ROC)
+
+ggplot(svmFitDF, aes(x=sigma, y=C)) +
+  geom_tile(aes(fill=ROC))
+
+svmFit$bestTune
+```
+
+```{r figure4.25sub2, cache=TRUE, dependson="figure4.10"}
+svmGrid <- expand.grid(sigma = c(5:1 %o% 10^(-3:-1)),
+                       C = seq(0.1,2.1,by=0.25)) 
+                         
+fitControl <- trainControl(method = "repeatedcv", 
+                           number = 10, 
+                           savePredictions = TRUE,
+                           summaryFunction=twoClassSummary, 
+                           classProbs=TRUE)
+
+svmFit <- train(Survived.yes ~ ., data = titanicTrain,
+                method = "svmRadial", 
+                trControl = fitControl, 
+                verbose = FALSE, 
+                tuneGrid = svmGrid)
+
+svmFitDF <- data.frame(sigma = as.factor(svmFit$results$sigma), 
+                       C = as.factor(svmFit$results$C),
+                       ROC = svmFit$results$ROC)
+
+ggplot(svmFitDF, aes(x=sigma, y=C)) +
+  geom_tile(aes(fill=ROC))
+
+svmFit$bestTune
+
+svmFit$results[row.names(svmFit$bestTune),]
+```

From 7d7e10bae806e12f1ed49dd76a330ffce795fc1b Mon Sep 17 00:00:00 2001
From: Paul Adamson <paul.adamson.01@gmail.com>
Date: Mon, 26 Dec 2016 08:28:48 -0500
Subject: [PATCH 4/4] Update README with info on R files.

---
 README.md | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 38a100d..a5c5e2e 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,7 @@
 # Real-World Machine Learning
 
-This repository includes a set of IPython notebooks with functions and example
+This repository includes a set of IPython and R Markdown notebooks with 
+functions and example
 code from the Real-World Machine Learning book by Henrik Brink, Joseph
 Richards, and Mark Fetherolf published by Manning Books:
 
@@ -10,7 +11,11 @@ The notebooks can be viewed directly in the Github repository:
 
 https://github.com/brinkar/real-world-machine-learning
 
-To run the code yourself, please clone this repo, install IPython (on Python 2.7) and run
+To run the Python code yourself, please clone this repo, install IPython (on Python 2.7) and run
 `ipython notebook` in the repo directory. We highly recommend installing the Anaconda Python distribution, which includes IPython, Scikit-Learn, Pandas, and mostly anything else used throughout the book:
 
 https://www.continuum.io/downloads
+
+To run the R code, please clone this repo, install [RStudio](https://www.rstudio.com/products/rstudio/download/)
+and either `Knit` the entire document (⇧⌘K) or run each code chunk
+interactively.