diff --git a/8_PREDMACHLEARN/.directory b/8_PREDMACHLEARN/.directory
new file mode 100644
index 0000000..f88ed35
--- /dev/null
+++ b/8_PREDMACHLEARN/.directory
@@ -0,0 +1,4 @@
+[Dolphin]
+Timestamp=2017,2,16,14,12,56
+Version=3
+ViewMode=1
diff --git a/8_PREDMACHLEARN/Practical Machine Learning Course Notes.Rmd b/8_PREDMACHLEARN/Practical Machine Learning Course Notes.Rmd
index 96fae3d..ace370a 100644
--- a/8_PREDMACHLEARN/Practical Machine Learning Course Notes.Rmd
+++ b/8_PREDMACHLEARN/Practical Machine Learning Course Notes.Rmd
@@ -13,6 +13,10 @@ header-includes: \usepackage{graphicx} \usepackage{mathtools}
---
+```{r setup, include=FALSE, error=TRUE, message=FALSE}
+knitr::opts_chunk$set(echo = TRUE, comment = NA, error = TRUE)
+```
+
$\pagebreak$
@@ -427,7 +431,8 @@ folds$test[[1]]
```{r}
# returns the arguments of the default train function
-args(train.default)
+# args(train.default) <- this dosn't work anymore in caret
+args(train)
```
* `train` function has a large set of parameters, below are the default options
@@ -774,8 +779,11 @@ testing <- spam[-inTrain,]
preProc <- preProcess(log10(training[,-58]+1),method="pca",pcaComp=2)
# calculate PCs for training data
trainPC <- predict(preProc,log10(training[,-58]+1))
+# add variable `type` to trainPC
+type <- training$type
+trainPC <- data.frame(trainPC, type) # join variable `type`
# run model on outcome and principle components
-modelFit <- train(training$type ~ .,method="glm",data=trainPC)
+modelFit <- train(type ~ .,method="glm",data=trainPC)
# calculate PCs for test data
testPC <- predict(preProc,log10(testing[,-58]+1))
# compare results
@@ -789,7 +797,7 @@ confusionMatrix(testing$type,predict(modelFit,testPC))
```{r message = FALSE, warning = FALSE}
# construct model
-modelFit <- train(training$type ~ .,method="glm",preProcess="pca",data=training)
+modelFit <- train(type ~ .,method="glm",preProcess="pca",data=training)
# print results of model
confusionMatrix(testing$type,predict(modelFit,testing))
```
@@ -1321,6 +1329,7 @@ pred.lda
- ***example: `caret` package***
```{r message = F, warning = F}
+# package needed: klaR
# using the same data from iris, run naive Bayes on training data
nb <- train(Species ~ ., data=training,method="nb")
# predict test outcomes using naive Bayes model
diff --git a/8_PREDMACHLEARN/Practical_Machine_Learning_Course_Notes.html b/8_PREDMACHLEARN/Practical_Machine_Learning_Course_Notes.html
index 1c9b791..826225d 100644
--- a/8_PREDMACHLEARN/Practical_Machine_Learning_Course_Notes.html
+++ b/8_PREDMACHLEARN/Practical_Machine_Learning_Course_Notes.html
@@ -8,38 +8,58 @@
+
Practical Machine Learning Course Notes
-
+
-
-
-
-
+
+
+
+
+
+
+
@@ -63,17 +110,43 @@
color: inherit;
background-color: rgba(0, 0, 0, 0.04);
}
-img {
- max-width:100%;
- height: auto;
+img {
+ max-width:100%;
+ height: auto;
+}
+.tabbed-pane {
+ padding-top: 12px;
+}
+button.code-folding-btn:focus {
+ outline: none;
}
+
+
+
+
+
+
+
+
+
+
-
-
Practical Machine Learning Course Notes
+
+
+
+
+
+
+
Practical Machine Learning Course Notes
Xing Su
+
@@ -155,28 +228,29 @@
Xing Su
+
\(\pagebreak\)
Prediction
-
process for prediction = population \(\rightarrow\) probability and sampling to pick set of data \(\rightarrow\) split into training and test set \(\rightarrow\) build prediction function \(\rightarrow\) predict for new data \(\rightarrow\) evaluate
+
process for prediction = population \(\rightarrow\) probability and sampling to pick set of data \(\rightarrow\) split into training and test set \(\rightarrow\) build prediction function \(\rightarrow\) predict for new data \(\rightarrow\) evaluate
Note: choosing the right dataset and knowing what the specific question is are paramount to the success of the prediction algorithm (GoogleFlu failed to predict accurately when people’s search habits changed)
-
+
-
components of predictor = question \(\rightarrow\) input data \(\rightarrow\) features (extracting variables/characteristics)\(\rightarrow\) algorithm \(\rightarrow\) parameters (estimate)\(\rightarrow\) evaluation
+
components of predictor = question \(\rightarrow\) input data \(\rightarrow\) features (extracting variables/characteristics)\(\rightarrow\) algorithm \(\rightarrow\) parameters (estimate)\(\rightarrow\) evaluation
relative order of importance = question (concrete/specific) > data (relevant) > features (properly extract) > algorithms
data selection
-
Note: “garbage in = garbage out” \(\rightarrow\) having the correct/relevant data will decide whether the model is successful
+
Note: “garbage in = garbage out” \(\rightarrow\) having the correct/relevant data will decide whether the model is successful
data for what you are trying to predict is most helpful
-
more data \(\rightarrow\) better models (usually)
+
more data \(\rightarrow\) better models (usually)
feature selection
-
good features \(\rightarrow\) lead to data compression, retain relevant information, created based on expert domain knowledge
-
common mistakes \(\rightarrow\) automated feature selection (can yield good results but likely to behave inconsistently with slightly different data), not understanding/dealing with skewed data/outliers, throwing away information unnecessarily
+
good features \(\rightarrow\) lead to data compression, retain relevant information, created based on expert domain knowledge
+
common mistakes \(\rightarrow\) automated feature selection (can yield good results but likely to behave inconsistently with slightly different data), not understanding/dealing with skewed data/outliers, throwing away information unnecessarily
algorithm selection
@@ -218,16 +292,16 @@
In Sample vs Out of Sample Errors
example
-
# load data
+
# load datalibrary(kernlab); data(spam); set.seed(333)
# picking a small subset (10 values) from spam data set
smallSpam <-spam[sample(dim(spam)[1],size=10),]
# label spam = 2 and ham = 1
spamLabel <-(smallSpam$type=="spam")*1 +1# plot the capitalAve values for the dataset with colors differentiated by spam/ham (2 vs 1)
-plot(smallSpam$capitalAve,col=spamLabel)
-
-
# first rule (over-fitting to capture all variation)
+plot(smallSpam$capitalAve,col=spamLabel)
+
+
# first rule (over-fitting to capture all variation)
rule1 <-function(x){
prediction <-rep(NA,length(x))
prediction[x >2.7] <- "spam"
@@ -237,12 +311,12 @@
In Sample vs Out of Sample Errors
return(prediction)
}
# tabulate results of prediction algorithm 1 (in sample error -> no error in this case)
-table(rule1(smallSpam$capitalAve),smallSpam$type)
-
##
-## nonspam spam
-## nonspam 5 0
-## spam 0 5
-
# second rule (simple, setting a threshold)
+table(rule1(smallSpam$capitalAve),smallSpam$type)
+
+ nonspam spam
+ nonspam 5 0
+ spam 0 5
+
# second rule (simple, setting a threshold)
rule2 <-function(x){
prediction <-rep(NA,length(x))
prediction[x >2.8] <- "spam"
@@ -250,32 +324,32 @@
In Sample vs Out of Sample Errors
return(prediction)
}
# tabulate results of prediction algorithm 2(in sample error -> 10% in this case)
-table(rule2(smallSpam$capitalAve),smallSpam$type)
-
##
-## nonspam spam
-## nonspam 5 1
-## spam 0 4
-
# tabulate out of sample error for algorithm 1
-table(rule1(spam$capitalAve),spam$type)
createDataPartition(y=data$var, times=1, p=0.75, list=FALSE)\(\rightarrow\) creates data partitions using given variable
+
createDataPartition(y=data$var, times=1, p=0.75, list=FALSE)\(\rightarrow\) creates data partitions using given variable
y=data$var = specifies what outcome/variable to split the data on
times=1 = specifies number of partitions to create (number of data splitting performed)
@@ -593,7 +667,7 @@
Data Slicing
example
-
# load packages and data
+
# load packages and datalibrary(caret)
# create training set indexes with 75% of data
inTrain <-createDataPartition(y=spam$type,p=0.75, list=FALSE)
@@ -602,18 +676,18 @@
Data Slicing
# subset spam data (the rest) to test
testing <-spam[-inTrain,]
# dimension of original and training dataset
-rbind("original dataset" =dim(spam),"training set" =dim(training))
-
## [,1] [,2]
-## original dataset 4601 58
-## training set 3451 58
[,1] [,2]
+original dataset 4601 58
+training set 3451 58
-
createFolds(y=data$var, k=10, list=TRUE, returnTrain=TRUE) = slices the data in to \(k\) folds for cross validation and returns \(k\) lists of indices
+
createFolds(y=data$var, k=10, list=TRUE, returnTrain=TRUE) = slices the data in to \(k\) folds for cross validation and returns \(k\) lists of indices
y=data$var = specifies what outcome/variable to split the data on
each training set has approximately has \(\frac{k-1}{k} \%\) of the data (in this case 90%)
-
each training set has approximately has \(\frac{1}{k} \%\) of the data (in this case 10%)
+
each training set has approximately has \(\frac{k-1}{k} \%\) of the data (in this case 90%)
+
each training set has approximately has \(\frac{1}{k} \%\) of the data (in this case 10%)
list=TRUE = returns k list of indices that corresponds to the cross-validated sets
@@ -631,39 +705,39 @@
Data Slicing
example
-
# create 10 folds for cross validation and return the training set indices
+
# create 10 folds for cross validation and return the training set indices
folds <-createFolds(y=spam$type,k=10,list=TRUE,returnTrain=TRUE)
# structure of the training set indices
-str(folds)
# return the test set indices instead# note: returnTrain = FALSE is unnecessary as it is the default behavior
folds.test <-createFolds(y=spam$type,k=10,list=TRUE,returnTrain=FALSE)
-str(folds.test)
# return first 10 elements of the first training set
+folds[[1]][1:10]
+
[1] 1 2 3 4 6 7 8 9 10 12
createResample(y=data$var, times=10, list=TRUE) = create 10 resamplings from the given data with replacement
@@ -674,21 +748,21 @@
Data Slicing
times=10 = number of samples to create
-
# create 10 resamples
+
# create 10 resamples
resamples <-createResample(y=spam$type,times=10,list=TRUE)
# structure of the resamples (note some samples are repeated)
-str(resamples)
table(cutVariable, data$var2) = tabulates the cut factor variable vs another variable in the dataset (ie; builds a contingency table using cross-classifying factors)
prop.table(table, margin=1) = converts a table to a proportion table
@@ -867,22 +941,22 @@
some predictors may have strange distributions (i.e. skewed) and may need to be transformed to be more useful for prediction algorithm
-
particularly true for model based algorithms \(\rightarrow\) naive Bayes, linear discriminate analysis, linear regression
+
particularly true for model based algorithms \(\rightarrow\) naive Bayes, linear discriminate analysis, linear regression
centering = subtracting the observations of a particular variable by its mean
scaling = dividing the observations of a particular variable by its standard deviation
-
normalizing = centering and scaling the variable \(\rightarrow\) effectively converting each observation to the number of standard deviations away from the mean
+
normalizing = centering and scaling the variable \(\rightarrow\) effectively converting each observation to the number of standard deviations away from the mean
the distribution of the normalized variable will have a mean of 0 and standard deviation of 1
Note: normalizing data can help remove bias and high variability, but may not be applicable in all cases
# set up BoxCox transforms
preObj <-preProcess(training[,-58],method=c("BoxCox"))
# perform preprocessing on training data
trainCapAveS <-predict(preObj,training[,-58])$capitalAve
# plot histogram and QQ Plot# Note: the transformation definitely helped to# normalize the data but it does not produce perfect result
-par(mfrow=c(1,2)); hist(trainCapAveS); qqnorm(trainCapAveS)
# Make some values NA
training$capAve <-training$capitalAve
selectNA <-rbinom(dim(training)[1],size=1,prob=0.05)==1
training$capAve[selectNA] <-NA
@@ -977,10 +1051,10 @@
generally, more knowledge and understanding you have of the system/data, the easier it will be to extract the summarizing features
-
when in doubt, more features is always safer \(\rightarrow\) lose less information and the features can be filtered during model construction
+
when in doubt, more features is always safer \(\rightarrow\) lose less information and the features can be filtered during model construction
this process can be automated (i.e. PCA) but generally have to be very careful, as one very useful feature in the training data set may not have as much effect on the test data set
Note: science is the key here, Google “feature extraction for [data type]” for more guidance
@@ -1021,31 +1095,31 @@
Covariate Creation/Feature Extraction
Creating Dummy Variables
-
convert factor variables to indicator/dummy variable \(\rightarrow\) qualitative become quantitative
+
convert factor variables to indicator/dummy variable \(\rightarrow\) qualitative become quantitative
dummyVars(outcome~var, data=training) = creates a dummy variable object that can be used through predict function to create dummy variables
predict(dummyObj, newdata=training) = creates appropriate columns to represent the factor variable with appropriate 0s and 1s
-
2 factor variable \(\rightarrow\) two columns which have 0 or 1 depending on the outcome
-
3 factor variable \(\rightarrow\) three columns which have 0, 0, and 1 representing the outcome
+
2 factor variable \(\rightarrow\) two columns which have 0 or 1 depending on the outcome
+
3 factor variable \(\rightarrow\) three columns which have 0, 0, and 1 representing the outcome
Note: only one of the columns can have values of 1 for each observation
-
# setting up data
+
# setting up data
inTrain <-createDataPartition(y=Wage$wage,p=0.7, list=FALSE)
training <-Wage[inTrain,]; testing <-Wage[-inTrain,]
# create a dummy variable object
dummies <-dummyVars(wage ~jobclass,data=training)
# create the dummy variable columns
-head(predict(dummies,newdata=training))
# plot all age vs wage dataplot(training$age,training$wage,pch=19,cex=0.5)
# plot the fitted polynomial function
-points(training$age,predict(lm1,newdata=training),col="red",pch=19,cex=0.5)
-
-
# predict on test values
-head(predict(bsBasis,age=testing$age))
Note: once registered, you should see in your task manager/activity monitor that 4 “R Session” appear when you run your code
-
\(\pagebreak\)
+
\(\pagebreak\)
@@ -1152,7 +1226,7 @@
prcomp Function
pr<-prcomp(data) = performs PCA on all variables and returns a prcomp object that contains information about standard deviations and rotations
-
pr$rotations = returns eigenvectors for the linear combinations of all variables (coefficients that variables are multiplied by to come up with the principal components) \(\rightarrow\) how the principal components are created
+
pr$rotations = returns eigenvectors for the linear combinations of all variables (coefficients that variables are multiplied by to come up with the principal components) \(\rightarrow\) how the principal components are created
often times, it is useful to take the log transformation of the variables and adding 1 before performing PCA
helps to reduce skewness or strange distribution in data
@@ -1165,23 +1239,23 @@
prcomp Function
-
# load spam data
+
# load spam datadata(spam)
# perform PCA on dataset
prComp <-prcomp(log10(spam[,-58]+1))
# print out the eigenvector/rotations first 5 rows and PCs
-head(prComp$rotation[, 1:5], 5)
# create new variable that marks spam as 2 and nospam as 1
typeColor <-((spam$type=="spam")*1 +1)
# plot the first two principal components
-plot(prComp$x[,1],prComp$x[,2],col=typeColor,xlab="PC1",ylab="PC2")
Note: the same PCA must be performed on the test set
-
# create train and test sets
+
# create train and test sets
inTrain <-createDataPartition(y=spam$type,p=0.75, list=FALSE)
training <-spam[inTrain,]
testing <-spam[-inTrain,]
@@ -1209,117 +1283,120 @@
caret Package
preProc <-preProcess(log10(training[,-58]+1),method="pca",pcaComp=2)
# calculate PCs for training data
trainPC <-predict(preProc,log10(training[,-58]+1))
+# add variable `type` to trainPC
+type <-training$type
+trainPC <-data.frame(trainPC, type) # join variable `type`# run model on outcome and principle components
-modelFit <-train(training$type ~.,method="glm",data=trainPC)
+modelFit <-train(type ~.,method="glm",data=trainPC)
# calculate PCs for test data
testPC <-predict(preProc,log10(testing[,-58]+1))
# compare results
-confusionMatrix(testing$type,predict(modelFit,testPC))
alternatively, PCA can be directly performed with the train method
train(outcome ~ ., method="glm", preProcess="pca", data=training) = performs PCA first on the training set and then runs the specified model
-
effectively the same procedures as above (preProcess\(\rightarrow\)predict)
+
effectively the same procedures as above (preProcess\(\rightarrow\)predict)
Note: in both cases, the PCs were able to achieve 90+% accuracy
-
# construct model
-modelFit <-train(training$type ~.,method="glm",preProcess="pca",data=training)
+
# construct model
+modelFit <-train(type ~.,method="glm",preProcess="pca",data=training)
# print results of model
-confusionMatrix(testing$type,predict(modelFit,testing))
where \(\hat \beta_0\) is the estimated intercept (when all variables are 0)
-
\(\hat \beta_1, \hat \beta_2, \ldots, \hat \beta_p\) are the estimated coefficients
-
\(X_{1i}, X_{2i}, \ldots, X_{pi}\) are the variables/covariates
-
\(\hat Y_i\) is the predicted outcome
+
where \(\hat \beta_0\) is the estimated intercept (when all variables are 0)
+
\(\hat \beta_1, \hat \beta_2, \ldots, \hat \beta_p\) are the estimated coefficients
+
\(X_{1i}, X_{2i}, \ldots, X_{pi}\) are the variables/covariates
+
\(\hat Y_i\) is the predicted outcome
R Commands and Examples
-
lm<-lm(y ~ x, data=train) = runs a linear model of outcome y on predictor x \(\rightarrow\) univariate regression
+
lm<-lm(y ~ x, data=train) = runs a linear model of outcome y on predictor x \(\rightarrow\) univariate regression
-
summary(lm) = returns summary of the linear regression model, which will include coefficients, standard errors, \(t\) statistics, and p values
+
summary(lm) = returns summary of the linear regression model, which will include coefficients, standard errors, \(t\) statistics, and p values
lm(y ~ x1+x2+x3, data=train) = run linear model of outcome y on predictors x1, x2, and x3
lm(y ~ ., data=train = run linear model of outcome y on all predictors
-
predict(lm, newdata=df) = use the constructed linear model to predict outcomes (\(\hat Y_i\)) for the new values
+
predict(lm, newdata=df) = use the constructed linear model to predict outcomes (\(\hat Y_i\)) for the new values
newdata data frame must have the same variables (factors must have the same levels) as the training data
newdata=test = predict outcomes for the test set based on linear regression model from the training
@@ -1327,10 +1404,10 @@
R Commands and Examples
RSME can be calculated to measure the accuracy of the linear model
-
Note: \(RSME_{test}\), which estimates the out-of-sample error, is almost always GREATER than \(RSME_{train}\)
+
Note: \(RSME_{test}\), which estimates the out-of-sample error, is almost always GREATER than \(RSME_{train}\)
-
# load data
+
# load datadata(faithful)
# create train and test sets
inTrain <-createDataPartition(y=faithful$waiting, p=0.5, list=FALSE)
@@ -1338,31 +1415,31 @@
R Commands and Examples
# build linear model
lm1 <-lm(eruptions ~waiting,data=trainFaith)
# print summary of linear model
-summary(lm1)
-
##
-## Call:
-## lm(formula = eruptions ~ waiting, data = trainFaith)
-##
-## Residuals:
-## Min 1Q Median 3Q Max
-## -1.24867 -0.36292 0.00002 0.35768 1.19858
-##
-## Coefficients:
-## Estimate Std. Error t value Pr(>|t|)
-## (Intercept) -2.165648 0.227486 -9.52 <2e-16 ***
-## waiting 0.079396 0.003146 25.24 <2e-16 ***
-## ---
-## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## Residual standard error: 0.5013 on 135 degrees of freedom
-## Multiple R-squared: 0.8251, Adjusted R-squared: 0.8238
-## F-statistic: 636.9 on 1 and 135 DF, p-value: < 2.2e-16
-
# predict eruptions for new waiting time
+summary(lm1)
+
+Call:
+lm(formula = eruptions ~ waiting, data = trainFaith)
+
+Residuals:
+ Min 1Q Median 3Q Max
+-1.30246 -0.40746 0.03955 0.40465 1.19221
+
+Coefficients:
+ Estimate Std. Error t value Pr(>|t|)
+(Intercept) -1.858966 0.237636 -7.823 1.33e-12 ***
+waiting 0.075444 0.003287 22.952 < 2e-16 ***
+---
+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+
+Residual standard error: 0.5272 on 135 degrees of freedom
+Multiple R-squared: 0.796, Adjusted R-squared: 0.7945
+F-statistic: 526.8 on 1 and 135 DF, p-value: < 2.2e-16
+
# predict eruptions for new waiting time
newdata <-data.frame(waiting=80)
-predict(lm1,newdata)
-
## 1
-## 4.186
-
# create 1 x 2 panel plot
+predict(lm1,newdata)
+
1
+4.176566
+
# create 1 x 2 panel plotpar(mfrow=c(1,2))
# plot train data with the regression lineplot(trainFaith$waiting,trainFaith$eruptions,pch=19,col="blue",xlab="Waiting",
@@ -1371,30 +1448,30 @@
R Commands and Examples
# plot test data with the regression lineplot(testFaith$waiting,testFaith$eruptions,pch=19,col="blue",xlab="Waiting",
ylab="Duration", main ="Test")
-lines(testFaith$waiting,predict(lm1,newdata=testFaith),lwd=3)
-
-
# Calculate RMSE on training and test sets
+lines(testFaith$waiting,predict(lm1,newdata=testFaith),lwd=3)
+
+
# Calculate RMSE on training and test setsc(trainRMSE =sqrt(sum((lm1$fitted-trainFaith$eruptions)^2)),
- testRMSE =sqrt(sum((predict(lm1,newdata=testFaith)-testFaith$eruptions)^2)))
pi<-predict(lm, newdata=test, interval="prediction") = returns 3 columns for fit (predicted value, same as before), lwr (lower bound of prediction interval), and upr (upper bound of prediction interval)
matlines(x, pi, type="l") = plots three lines, one for the linear fit and two for upper/lower prediction interval bounds
-
# calculate prediction interval
+
# calculate prediction interval
pred1 <-predict(lm1,newdata=testFaith,interval="prediction")
# plot data points (eruptions, waiting)plot(testFaith$waiting,testFaith$eruptions,pch=19,col="blue")
# plot fit line and prediction interval
-matlines(testFaith$waiting,pred1,type="l",,col=c(1,2,2),lty =c(1,1,1), lwd=3)
lm <- train(y ~ x, method="lm", data=train) = run linear model on the training data \(\rightarrow\) identical to lm function
+
lm <- train(y ~ x, method="lm", data=train) = run linear model on the training data \(\rightarrow\) identical to lm function
-
summary(lm$finalModel) = returns summary of the linear regression model, which will include coefficients, standard errors, \(t\) statistics, and p values \(\rightarrow\) identical to summary(lm) for a lm object
+
summary(lm$finalModel) = returns summary of the linear regression model, which will include coefficients, standard errors, \(t\) statistics, and p values \(\rightarrow\) identical to summary(lm) for a lm object
train(y ~ ., method="lm", data=train) = run linear model on all predictors in training data
multiple predictors (dummy/indicator variables) are created for factor variables
@@ -1409,7 +1486,7 @@
R Commands and Examples
-
# create train and test sets
+
# create train and test sets
inTrain <-createDataPartition(y=Wage$wage,p=0.7, list=FALSE)
training <-Wage[inTrain,]; testing <-Wage[-inTrain,]
# fit linear model for age jobclass and education
@@ -1419,14 +1496,14 @@
R Commands and Examples
# set up 2 x 2 panel plotpar(mfrow =c(2, 2))
# construct diagnostic plots for model
-plot(finMod,pch=19,cex=0.5,col="#00000010")
-
+plot(finMod,pch=19,cex=0.5,col="#00000010")
+
plotting residuals by fitted values and coloring with a variable not used in the model helps spot a trend in that variable.
-
# plot fitted values by residuals
-qplot(finMod$fitted, finMod$residuals, color=race, data=training)
-
+
# plot fitted values by residuals
+qplot(finMod$fitted, finMod$residuals, color=race, data=training)
+
plotting residuals by index (ie; row numbers) can be helpful in showing missing variables
@@ -1437,24 +1514,24 @@
R Commands and Examples
-
# plot residual by index
-plot(finMod$residuals,pch=19,cex=0.5)
-
+
# plot residual by index
+plot(finMod$residuals,pch=19,cex=0.5)
+
here the residuals increase linearly with the index, and the highest residuals are concentrated in the higher indexes, so there must be a missing variable
-
\(\pagebreak\)
+
\(\pagebreak\)
Prediction with Trees
-
prediction with trees = iteratively split variables into groups (effectively constructing decision trees) \(\rightarrow\) produces nonlinear model
+
prediction with trees = iteratively split variables into groups (effectively constructing decision trees) \(\rightarrow\) produces nonlinear model
-
the classification tree uses interactions between variables \(\rightarrow\) the ultimate groups/leafs may depend on many variables
+
the classification tree uses interactions between variables \(\rightarrow\) the ultimate groups/leafs may depend on many variables
the result (tree) is easy to interpret, and generally performs better predictions than regression models when the relationships are non-linear
-
transformations less important \(\rightarrow\) monotone transformations (order unchanged, such as \(\log\)) will produce same splits
+
transformations less important \(\rightarrow\) monotone transformations (order unchanged, such as \(\log\)) will produce same splits
trees can be used for regression problems as well and use RMSE as measure of impurity
however, without proper cross-validation, the model can be over-fitted (especially with large number of variables) and results may be variable from one run to the next
\(\hat p_mk\) is the probability of the objects in group \(m\) to take on the classification \(k\)
-
\(N_m\) is the size of the group
-
Misclassification Error\[ 1 - \hat{p}_{m~k(m)}\] where \(k(m)\) is the most common classification/group
+
\(\hat p_mk\) is the probability of the objects in group \(m\) to take on the classification \(k\)
+
\(N_m\) is the size of the group
+
Misclassification Error\[ 1 - \hat{p}_{m~k(m)}\] where \(k(m)\) is the most common classification/group
0 = perfect purity
0.5 = no purity
-
Note: it is not 1 here because when \(\hat{p}_{m~k(m)} < 0.5\) or there’s not predominant classification for the objects, it means the group should be further subdivided until there’s a majority
+
Note: it is not 1 here because when \(\hat{p}_{m~k(m)} < 0.5\) or there’s not predominant classification for the objects, it means the group should be further subdivided until there’s a majority
Information Gai\[ -\sum_{k=1}^K \hat{p}_{mk} \log_2\hat{p}_{mk} \]
+
Information Gai\[ -\sum_{k=1}^K \hat{p}_{mk} \log_2\hat{p}_{mk} \]
0 = perfect purity
1 = no purity
example
-
# set margin and seed
+
# set margin and seedpar(mar=c(1,1,1,1), mfrow =c(1, 2)); set.seed(1234);
# simulate data
x =rep(1:4,each=4); y =rep(1:4,4)
# plot first scenarioplot(x,y,xaxt="n",yaxt="n",cex=3,col=c(rep("blue",15),rep("red",1)),pch=19)
# plot second scenario
-plot(x,y,xaxt="n",yaxt="n",cex=3,col=c(rep("blue",8),rep("red",8)),pch=19)
# load datalibrary(ElemStatLearn); data(ozone,package="ElemStatLearn")
# reorder rows based on ozone variable
ozone <-ozone[order(ozone$ozone),]
@@ -1622,8 +1699,8 @@
Bagging
# plot each prediction model
for(i in 1:10){lines(1:155,ll[i,],col="grey",lwd=2)}
# plot the average in red
-lines(1:155,apply(ll,2,mean),col="red",lwd=2)
-
+lines(1:155,apply(ll,2,mean),col="red",lwd=2)
+
Bagging Algorithms
@@ -1648,7 +1725,7 @@
Bagging Algorithms
example
-
# load relevant package and data
+
# load relevant package and datalibrary(party); data(ozone,package="ElemStatLearn")
# reorder rows based on ozone variable
ozone <-ozone[order(ozone$ozone),]
@@ -1667,9 +1744,9 @@
Bagging Algorithms
# plot the first fitpoints(ozone$ozone,predict(treebag$fits[[1]]$fit,predictors),pch=19,col="red")
# plot the aggregated predictions
-points(ozone$ozone,predict(treebag,predictors),pch=19,col="blue")
# apply random forest
modFit <-train(Species~.,data=training,method="rf",prox=TRUE)
# return the second tree (first 6 rows)
-head(getTree(modFit$finalModel,k=2))
# compute cluster centers
irisP <-classCenter(training[,c(3,4)], training$Species, modFit$finalModel$prox)
# convert irisP to data frame and add Species column
irisP <-as.data.frame(irisP); irisP$Species <-rownames(irisP)
# plot data points
p <-qplot(Petal.Width, Petal.Length, col=Species,data=training)
# add the cluster centers
-p +geom_point(aes(x=Petal.Width,y=Petal.Length,col=Species),size=5,shape=4,data=irisP)
-
-
# predict outcome for test data set using the random forest model
+p +geom_point(aes(x=Petal.Width,y=Petal.Length,col=Species),size=5,shape=4,data=irisP)
+
+
# predict outcome for test data set using the random forest model
pred <-predict(modFit,testing)
# logic value for whether or not the rf algorithm predicted correctly
testing$predRight <-pred==testing$Species
# tabulate results
-table(pred,testing$Species)
# plot data points with the incorrect classification highlighted
-qplot(Petal.Width,Petal.Length,colour=predRight,data=testing,main="newdata Predictions")
# plot data points with the incorrect classification highlighted
+qplot(Petal.Width,Petal.Length,colour=predRight,data=testing,main="newdata Predictions")
+
+
\(\pagebreak\)
@@ -1774,28 +1851,28 @@
Boosting
boosting = one of the most widely used and accurate prediction models, along with random forest
boosting can be done with any set of classifiers, and a well-known approach is gradient boosting
process: take a group of weak predictors \(\rightarrow\) weight them and add them up \(\rightarrow\) result in a stronger predictor
+
process: take a group of weak predictors \(\rightarrow\) weight them and add them up \(\rightarrow\) result in a stronger predictor
-
start with a set of classifiers \(h_1, \ldots, h_k\)
+
start with a set of classifiers \(h_1, \ldots, h_k\)
examples: all possible trees, all possible regression models, all possible cutoffs (divide data into different parts)
-
calculate a weighted sum of classifiers as the prediction value \[f(x) = \sum_i \alpha_i h_i(x)\] where \(\alpha_i\) = coefficient/weight and \(h_i(x)\) = value of classifier
+
calculate a weighted sum of classifiers as the prediction value \[f(x) = \sum_i \alpha_i h_i(x)\] where \(\alpha_i\) = coefficient/weight and \(h_i(x)\) = value of classifier
goal = minimize error (on training set)
-
select one \(h\) at each step (iterative)
+
select one \(h\) at each step (iterative)
calculate weights based on errors
-
up-weight missed classifications and select next \(h\)
+
up-weight missed classifications and select next \(h\)
example
-
+
we start with space with blue + and red - and the goal is to classify all the object correctly
only straight lines will be used for classification
-
+
from the above, we can see that a group of weak predictors (lines in this case), can be combined and weighed to become a much stronger predictor
@@ -1816,7 +1893,7 @@
R Commands and Examples
predict function can be used to apply the model to test data, similar to the rest of the algorithms in caret package
example
-
# load data
+
# load datadata(Wage)
# remove log wage variable (we are trying to predict wage)
Wage <-subset(Wage,select=-c(logwage))
@@ -1826,37 +1903,35 @@
R Commands and Examples
# run the gbm model
modFit <-train(wage ~., method="gbm",data=training,verbose=FALSE)
# print model summary
-print(modFit)
-
## Stochastic Gradient Boosting
-##
-## 2102 samples
-## 10 predictor
-##
-## No pre-processing
-## Resampling: Bootstrapped (25 reps)
-##
-## Summary of sample sizes: 2102, 2102, 2102, 2102, 2102, 2102, ...
-##
-## Resampling results across tuning parameters:
-##
-## interaction.depth n.trees RMSE Rsquared RMSE SD Rsquared SD
-## 1 50 35.64972 0.3317422 1.495165 0.02312183
-## 1 100 34.95593 0.3429594 1.503223 0.02319651
-## 1 150 34.84473 0.3451634 1.496016 0.02324176
-## 2 50 34.91119 0.3462101 1.490004 0.02425481
-## 2 100 34.74433 0.3487227 1.480423 0.02278780
-## 2 150 34.74823 0.3487136 1.472941 0.02314265
-## 3 50 34.83480 0.3467828 1.493846 0.02292988
-## 3 100 34.85342 0.3449018 1.482881 0.02373242
-## 3 150 34.99413 0.3401694 1.544378 0.02498133
-##
-## Tuning parameter 'shrinkage' was held constant at a value of 0.1
-##
-## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
-## RMSE was used to select the optimal model using the smallest value.
-## The final values used for the model were n.trees = 100,
-## interaction.depth = 2, shrinkage = 0.1 and n.minobsinnode = 10.
-
\(\pagebreak\)
+print(modFit)
+
Stochastic Gradient Boosting
+
+2102 samples
+ 10 predictors
+
+No pre-processing
+Resampling: Bootstrapped (25 reps)
+Summary of sample sizes: 2102, 2102, 2102, 2102, 2102, 2102, ...
+Resampling results across tuning parameters:
+
+ interaction.depth n.trees RMSE Rsquared
+ 1 50 35.24716 0.3061446
+ 1 100 34.74137 0.3165154
+ 1 150 34.70108 0.3179425
+ 2 50 34.67703 0.3199690
+ 2 100 34.58430 0.3226554
+ 2 150 34.67104 0.3202411
+ 3 50 34.59842 0.3219575
+ 3 100 34.73009 0.3183701
+ 3 150 34.97108 0.3107605
+
+Tuning parameter 'shrinkage' was held constant at a value of 0.1
+
+Tuning parameter 'n.minobsinnode' was held constant at a value of 10
+RMSE was used to select the optimal model using the smallest value.
+The final values used for the model were n.trees = 100,
+ interaction.depth = 2, shrinkage = 0.1 and n.minobsinnode = 10.
+
\(\pagebreak\)
@@ -1869,34 +1944,34 @@
Model Based Prediction
can be reasonably accurate on real problems
this approach does make additional assumptions about the data, which can lead to model failure/reduced accuracy if they are too far off
-
goal = build parameter-based model (based on probabilities) for conditional distribution \(P(Y = k~|~X = x)\), or the probability of the outcome \(Y\) is equal to a particular value \(k\) given a specific set of predictor variables \(x\)
+
goal = build parameter-based model (based on probabilities) for conditional distribution \(P(Y = k~|~X = x)\), or the probability of the outcome \(Y\) is equal to a particular value \(k\) given a specific set of predictor variables \(x\)
-
Note: \(X\) is the data for the model (observations for all predictor variables), which is also known as the design matrix
+
Note: \(X\) is the data for the model (observations for all predictor variables), which is also known as the design matrix
typical approach/process
-
start with the quantity \(P(Y = k~|~X = x)\)
-
apply Bayes’ Theorem such that \[ P(Y = k ~|~ X=x) = \frac{P(X=x~|~Y=k)P(Y=k)}{\sum_{\ell=1}^K P(X=x ~|~Y = \ell) P(Y=\ell)}\] where the denominator is simply the sum of probabilities for the predictor variables are the set specified in \(x\) for all outcomes of \(Y\)
-
assume the term \(P(X=x~|~Y=k)\) in the numerator follows a parameter-based probability distribution, or \(f_k(x)\)
+
start with the quantity \(P(Y = k~|~X = x)\)
+
apply Bayes’ Theorem such that \[ P(Y = k ~|~ X=x) = \frac{P(X=x~|~Y=k)P(Y=k)}{\sum_{\ell=1}^K P(X=x ~|~Y = \ell) P(Y=\ell)}\] where the denominator is simply the sum of probabilities for the predictor variables are the set specified in \(x\) for all outcomes of \(Y\)
+
assume the term \(P(X=x~|~Y=k)\) in the numerator follows a parameter-based probability distribution, or \(f_k(x)\)
-
common choice = Gaussian distribution\[f_k(x) = \frac{1}{\sigma_k \sqrt{2 \pi}}e^{-\frac{(x-\mu_k)^2}{2\sigma_k^2}}\]
+
common choice = Gaussian distribution\[f_k(x) = \frac{1}{\sigma_k \sqrt{2 \pi}}e^{-\frac{(x-\mu_k)^2}{2\sigma_k^2}}\]
-
assume the probability for the outcome \(Y\) to take on value of \(k\), or \(P(Y=k)\), is determined from the data to be some known quantity \(\pi_k\)
+
assume the probability for the outcome \(Y\) to take on value of \(k\), or \(P(Y=k)\), is determined from the data to be some known quantity \(\pi_k\)
so the quantity \(P(Y = k~|~X = x)\) can be rewritten as \[P(Y = k ~|~ X=x) = \frac{f_k(x) \pi_k}{\sum_{\ell = 1}^K f_{\ell}(x) \pi_{\ell}}\]
-
estimate the parameters (\(\mu_k\), \(\sigma_k^2\)) for the function \(f_k(x)\) from the data
-
calculate \(P(Y = k~|~X = x)\) using the parameters
-
the outcome \(Y\) is where the value of \(P(Y = k ~|~ X = x)\) is the highest
+
so the quantity \(P(Y = k~|~X = x)\) can be rewritten as \[P(Y = k ~|~ X=x) = \frac{f_k(x) \pi_k}{\sum_{\ell = 1}^K f_{\ell}(x) \pi_{\ell}}\]
+
estimate the parameters (\(\mu_k\), \(\sigma_k^2\)) for the function \(f_k(x)\) from the data
+
calculate \(P(Y = k~|~X = x)\) using the parameters
+
the outcome \(Y\) is where the value of \(P(Y = k ~|~ X = x)\) is the highest
prediction models that leverage this approach
-
linear discriminant analysis = assumes \(f_k(x)\) is multivariate Gaussian distribution with same covariance for each predictor variables
+
linear discriminant analysis = assumes \(f_k(x)\) is multivariate Gaussian distribution with same covariance for each predictor variables
effectively drawing lines through “covariate space”
-
quadratic discriminant analysis = assumes \(f_k(x)\) is multivariate Gaussian distribution with different covariance for predictor variables
+
quadratic discriminant analysis = assumes \(f_k(x)\) is multivariate Gaussian distribution with different covariance for predictor variables
effectively drawing curves through “covariate space”
@@ -1910,35 +1985,36 @@
Model Based Prediction
Linear Discriminant Analysis
-
to compare the probability for outcome \(Y = k\) versus probability for outcome \(Y = k\), we can look at the ratio of \[\frac{P(Y=k~|~X=x)}{P(Y=j~|~X=x)}\]
-
take the log of the ratio and apply Bayes’ Theorem, we get \[\log \frac{P(Y = k~|~X=x)}{P(Y = j~|~X=x)} = \log \frac{f_k(x)}{f_j(x)} + \log \frac{\pi_k}{\pi_j}\] which is effectively the log ratio of probability density functions plus the log ratio of prior probabilities
+
to compare the probability for outcome \(Y = k\) versus probability for outcome \(Y = k\), we can look at the ratio of \[\frac{P(Y=k~|~X=x)}{P(Y=j~|~X=x)}\]
+
take the log of the ratio and apply Bayes’ Theorem, we get \[\log \frac{P(Y = k~|~X=x)}{P(Y = j~|~X=x)} = \log \frac{f_k(x)}{f_j(x)} + \log \frac{\pi_k}{\pi_j}\] which is effectively the log ratio of probability density functions plus the log ratio of prior probabilities
-
Note: \(\log\) = monotone transformation, which means taking the \(\log\) of a quantity does not affect implication of the ratio since th \(\log(ratio)\) is directly correlated with ratio
+
Note: \(\log\) = monotone transformation, which means taking the \(\log\) of a quantity does not affect implication of the ratio since th \(\log(ratio)\) is directly correlated with ratio
-
if we substitute \(f_k(x)\) and \(f_l(x)\) with Gaussian probability density functions \[f(x) = \frac{1}{\sigma \sqrt{2 \pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\] so the ratio can be simplified to \[\log \frac{P(Y = k~|~X=x)}{P(Y = j~|~X=x)} = \log \frac{\pi_k}{\pi_j} - \frac{1}{2}(\mu_k + \mu_j)^T \Sigma^{-1}(\mu_k + \mu_j) + x^T \Sigma^{-1} (\mu_k - \mu_j)\] where \(\Sigma^{-1}\) = covariance matrix for the predictor variables, \(x^T\) = set of predictor variables, and \(\mu_k\) / \(\mu_j\) = mean of \(k\), \(j\) respectively
-
as annotated above, the log-ratio is effectively an equation of a line for a set of predictor variables \(x\)\(\rightarrow\) the first two terms are constants and the last term is in the form of \(X \beta\)
+
if we substitute \(f_k(x)\) and \(f_l(x)\) with Gaussian probability density functions \[f(x) = \frac{1}{\sigma \sqrt{2 \pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}\] so the ratio can be simplified to \[\log \frac{P(Y = k~|~X=x)}{P(Y = j~|~X=x)} = \log \frac{\pi_k}{\pi_j} - \frac{1}{2}(\mu_k + \mu_j)^T \Sigma^{-1}(\mu_k + \mu_j) + x^T \Sigma^{-1} (\mu_k - \mu_j)\] where \(\Sigma^{-1}\) = covariance matrix for the predictor variables, \(x^T\) = set of predictor variables, and \(\mu_k\) / \(\mu_j\) = mean of \(k\), \(j\) respectively
+
as annotated above, the log-ratio is effectively an equation of a line for a set of predictor variables \(x\)\(\rightarrow\) the first two terms are constants and the last term is in the form of \(X \beta\)
Note: the lines are also known as decision boundaries
-
therefore, we can classify values based on which side of the line the value is located (\(k\) vs \(j\))
-
discriminant functions are used to determine value of \(k\), the functions are in the form of \[\delta_k(x) = x^T \Sigma^{-1} \mu_k - \frac{1}{2}\mu_k \Sigma^{-1}\mu_k + log(\mu_k)\]
+
therefore, we can classify values based on which side of the line the value is located (\(k\) vs \(j\))
+
discriminant functions are used to determine value of \(k\), the functions are in the form of \[\delta_k(x) = x^T \Sigma^{-1} \mu_k - \frac{1}{2}\mu_k \Sigma^{-1}\mu_k + log(\mu_k)\]
-
plugging in the set of predictor variables, \(x^T\), into the discriminant function, we can find the value of \(k\) that maximizes the function \(\delta_k(x)\)
+
plugging in the set of predictor variables, \(x^T\), into the discriminant function, we can find the value of \(k\) that maximizes the function \(\delta_k(x)\)
the terms of the discriminant function can be estimated using maximum likelihood
-
the predicted value for the outcome is therefore \(\hat{Y}(x) = argmax_k \delta_k(x)\)
+
the predicted value for the outcome is therefore \(\hat{Y}(x) = argmax_k \delta_k(x)\)
example
-
classify a group of values into 3 groups using 2 variables (\(x\), \(y\) coordinates)
+
classify a group of values into 3 groups using 2 variables (\(x\), \(y\) coordinates)
3 lines are draw to split the data into 3 Gaussian distributions
-
each line splits the data into two groups \(\rightarrow\) 1 vs 2, 2 vs 3, 1 vs 3
-
each side of the line represents a region where the probability of one group (1, 2, or 3) is the highest
+
each line splits the data into two groups \(\rightarrow\) 1 vs 2, 2 vs 3, 1 vs 3
+
each side of the line represents a region where the probability of one group (1, 2, or 3) is the highest
+
the result is represented in the the following graph
-
+
R Commands
@@ -1947,7 +2023,7 @@
Linear Discriminant Analysis
example: caret package
-
# load data
+
# load datadata(iris)
# create training and test sets
inTrain <-createDataPartition(y=iris$Species,p=0.7, list=FALSE)
@@ -1958,27 +2034,27 @@
Linear Discriminant Analysis
# predict test outcomes using LDA model
pred.lda <-predict(lda,testing)
# print results
-pred.lda
for predictors \(X_1,\ldots,X_m\), we want to model \(P(Y = k ~|~ X_1,\ldots,X_m)\)
-
by applying Bayes’ Theorem, we get \[P(Y = k ~|~ X_1,\ldots,X_m) = \frac{\pi_k P(X_1,\ldots,X_m~|~ Y=k)}{\sum_{\ell = 1}^K P(X_1,\ldots,X_m ~|~ Y=k) \pi_{\ell}}\]
-
since the denominator is just a sum (constant), we can rewrite the quantity as \[P(Y = k ~|~ X_1,\ldots,X_m) \propto \pi_k P(X_1,\ldots,X_m~|~ Y=k)\] or the probability is proportional to the numerator
+
for predictors \(X_1,\ldots,X_m\), we want to model \(P(Y = k ~|~ X_1,\ldots,X_m)\)
+
by applying Bayes’ Theorem, we get \[P(Y = k ~|~ X_1,\ldots,X_m) = \frac{\pi_k P(X_1,\ldots,X_m~|~ Y=k)}{\sum_{\ell = 1}^K P(X_1,\ldots,X_m ~|~ Y=k) \pi_{\ell}}\]
+
since the denominator is just a sum (constant), we can rewrite the quantity as \[P(Y = k ~|~ X_1,\ldots,X_m) \propto \pi_k P(X_1,\ldots,X_m~|~ Y=k)\] or the probability is proportional to the numerator
Note: maximizing the numerator is the same as maximizing the ratio
-
\(\pi_k P(X_1,\ldots,X_m~|~ Y=k)\) can be rewritten as \[\begin{aligned}
+
\(\pi_k P(X_1,\ldots,X_m~|~ Y=k)\) can be rewritten as \[\begin{aligned}
\pi_k P(X_1,\ldots,X_m~|~ Y=k) & = \pi_k P(X_1 ~|~ Y = k)P(X_2,\ldots,X_m ~|~ X_1,Y=k) \\
& = \pi_k P(X_1 ~|~ Y = k) P(X_2 ~|~ X_1, Y=k) P(X_3,\ldots,X_m ~|~ X_1,X_2, Y=k) \\
& = \pi_k P(X_1 ~|~ Y = k) P(X_2 ~|~ X_1, Y=k)\ldots P(X_m~|~X_1\ldots,X_{m-1},Y=k) \\
@@ -1986,7 +2062,7 @@
Naive Bayes
this is effectively indicating that each of the predictors may be dependent on other predictors
-
however, if we make the assumption that all predictor variables are independent to each other, the quantity can be simplified to \[ \pi_k P(X_1,\ldots,X_m~|~ Y=k) \approx \pi_k P(X_1 ~|~ Y = k) P(X_2 ~|~ Y = k)\ldots P(X_m ~|~,Y=k)\] which is effectively the product of the prior probability for \(k\) and the probability of variables \(X_1,\ldots,X_m\) given that \(Y = k\)
+
however, if we make the assumption that all predictor variables are independent to each other, the quantity can be simplified to \[ \pi_k P(X_1,\ldots,X_m~|~ Y=k) \approx \pi_k P(X_1 ~|~ Y = k) P(X_2 ~|~ Y = k)\ldots P(X_m ~|~,Y=k)\] which is effectively the product of the prior probability for \(k\) and the probability of variables \(X_1,\ldots,X_m\) given that \(Y = k\)
Note: the assumption is naive in that it is unlikely the predictors are completely independent of each other, but this model still produces useful results particularly with large number of binary/categorical variables
@@ -2000,21 +2076,22 @@
Naive Bayes
example: caret package
-
# using the same data from iris, run naive Bayes on training data
+
# package needed: klaR
+# using the same data from iris, run naive Bayes on training data
nb <-train(Species ~., data=training,method="nb")
# predict test outcomes using naive Bayes model
pred.nb <-predict(nb,testing)
# print results
-pred.nb
# create logical variable that returns TRUE for when predictions from the two models match
equalPredictions <-(pred.lda==pred.nb)
# plot the comparison
-qplot(Petal.Width,Sepal.Width,colour=equalPredictions,data=testing)
as we can see from above, only one data point, which is located inbetween the two classes is predicted differently by the two models
-
\(\pagebreak\)
+
\(\pagebreak\)
@@ -2049,7 +2126,7 @@
Model Selection
the error for the prediction model on test set decreases first and then increases as number of predictors used approaches the total number of predictors available
-
+
this is expected since as more predictors used, the model is more likely to overfit the training data
goal in selecting models = avoid overfitting on training data and minimize error on test data
@@ -2077,7 +2154,7 @@
Example: Training vs Test Error for Combination of Predictors
to demonstrate the behavior of training and test errors, the prostate dataset from Elements of Statistical Learning is used
all combinations of predictors are used to produce prediction models, and Residual Squared Error (RSS) is calculated for all models on both the training and test sets
-
# load data and set seed
+
# load data and set seeddata(prostate); set.seed(1)
# define outcome y and predictors x
covnames <-names(prostate[-(9:10)])
@@ -2133,8 +2210,8 @@
Example: Training vs Test Error for Combination of Predictors
# plot line through the minimum test RSS data points in bluelines((1:p), minrss, col="red", lwd=1.7)
# add legend
-legend("topright", c("Train", "Test"), col=c("blue", "red"), pch=1)
from the above, we can clearly that test RSS error approaches the minimum at around 3 predictors and increases slightly as more predictors are used
@@ -2161,16 +2238,16 @@
Split Samples
Decompose Expected Prediction Error
-
the outcome \(Y_i\) can be modeled by \[Y_i = f(X_i) + \epsilon_i\] where \(\epsilon_i\) = error term
-
the expected prediction error is defined as \[EPE(\lambda) = E\left[\left(Y - \hat{f}_{\lambda}(X)\right)^2\right]\] where \(\lambda\) = specific set of tuning parameters
-
estimates from the model constructed with training data can be denoted as \(\hat{f}_{\lambda}(x^*)\) where \(X = x^*\) is the new data point that we would like to predict at
-
the expected prediction error is as follows \[\begin{aligned}
+
the outcome \(Y_i\) can be modeled by \[Y_i = f(X_i) + \epsilon_i\] where \(\epsilon_i\) = error term
+
the expected prediction error is defined as \[EPE(\lambda) = E\left[\left(Y - \hat{f}_{\lambda}(X)\right)^2\right]\] where \(\lambda\) = specific set of tuning parameters
+
estimates from the model constructed with training data can be denoted as \(\hat{f}_{\lambda}(x^*)\) where \(X = x^*\) is the new data point that we would like to predict at
goal of prediction model = minimize overall expected prediction error
-
irreducible error = noise inherent to the data collection process \(\rightarrow\) cannot be reduced
+
irreducible error = noise inherent to the data collection process \(\rightarrow\) cannot be reduced
bias/variance = can be traded in order to find optimal model (least error)
@@ -2186,9 +2263,9 @@
Hard Thresholding
hard thresholding can help estimate the coefficients/model by taking subsets of predictors and building models
process
-
model the outcome as \[Y_i = f(X_i) + \epsilon_i\] where \(\epsilon_i\) = error term
-
assume the prediction estimate has a linear form \[\hat{f}_{\lambda}(x) = x'\beta\] where only \(\lambda\) coefficients for the set of predictors \(x\) are nonzero
-
after setting the value of \(\lambda\), compute models using all combinations of \(\lambda\) variables to find which variables’ coefficients should be set to be zero
+
model the outcome as \[Y_i = f(X_i) + \epsilon_i\] where \(\epsilon_i\) = error term
+
assume the prediction estimate has a linear form \[\hat{f}_{\lambda}(x) = x'\beta\] where only \(\lambda\) coefficients for the set of predictors \(x\) are nonzero
+
after setting the value of \(\lambda\), compute models using all combinations of \(\lambda\) variables to find which variables’ coefficients should be set to be zero
problem
@@ -2199,43 +2276,43 @@
Hard Thresholding
as we can see from the results below, some of the coefficients have values of NA
-
# load prostate data
+
# load prostate datadata(prostate)
# create subset of observations with 10 variables
small =prostate[1:5,]
# print linear regression
-lm(lpsa ~.,data =small)
-
##
-## Call:
-## lm(formula = lpsa ~ ., data = small)
-##
-## Coefficients:
-## (Intercept) lcavol lweight age lbph
-## 9.60615 0.13901 -0.79142 0.09516 NA
-## svi lcp gleason pgg45 trainTRUE
-## NA NA -2.08710 NA NA
+lm(lpsa ~.,data =small)
+
+Call:
+lm(formula = lpsa ~ ., data = small)
+
+Coefficients:
+(Intercept) lcavol lweight age lbph
+ 9.60615 0.13901 -0.79142 0.09516 NA
+ svi lcp gleason pgg45 trainTRUE
+ NA NA -2.08710 NA NA
regularized regression = fit a regression model and adjust for the large coefficients in attempt to help with bias/variance trade-off or model selection
-
when running regressions unconstrained (without specifying any criteria for coefficients), the model may be susceptible to high variance (coefficients explode \(\rightarrow\) very large values) if there are variables that are highly correlated
+
when running regressions unconstrained (without specifying any criteria for coefficients), the model may be susceptible to high variance (coefficients explode \(\rightarrow\) very large values) if there are variables that are highly correlated
controlling/regularizing coefficients may slightly increase bias (lose a bit of prediction capability) but will reduce variance and improve the prediction error
however, this approach may be very demanding computationally and generally does not perform as well as random forest/boosting
-
Penalized Residual Sum of Squares (PRSS) is calculated by adding a penalty term to the prediction squared error \[PRSS(\beta) = \sum_{j=1}^n (Y_j - \sum_{i=1}^m \beta_{1i} X_{ij})^2 + P(\lambda; \beta)\]
+
Penalized Residual Sum of Squares (PRSS) is calculated by adding a penalty term to the prediction squared error \[PRSS(\beta) = \sum_{j=1}^n (Y_j - \sum_{i=1}^m \beta_{1i} X_{ij})^2 + P(\lambda; \beta)\]
penalty shrinks coefficients if their values become too large
penalty term is generally used to reduce complexity and variance for the model, while respecting the structure of the data/relationship
example: co-linear variables
-
given a linear model, \[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon\] where \(X_1\) and \(X_2\) are nearly perfectly correlated (co-linear)
-
the model can then be approximated by this model by omitting \(X_2\), so the model becomes \[Y = \beta_0 + (\beta_1 + \beta_2)X_1 + \epsilon\]
-
with the above model, we can get a good estimate of \(Y\)
+
given a linear model, \[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon\] where \(X_1\) and \(X_2\) are nearly perfectly correlated (co-linear)
+
the model can then be approximated by this model by omitting \(X_2\), so the model becomes \[Y = \beta_0 + (\beta_1 + \beta_2)X_1 + \epsilon\]
+
with the above model, we can get a good estimate of \(Y\)
k=4 = number of variables that should be retained
@@ -2292,26 +2369,26 @@
Regularized Regression - Ridge Regression
this means that length(predictors)-k variables will be eliminated
-
[caret package] predict(model,test) = use the model to predict on test set \(\rightarrow\) similar to all other caret algorithms
+
caret packagepredict(model,test) = use the model to predict on test set \(\rightarrow\) similar to all other caret algorithms
-
example: ridge coefficient paths vs \(\lambda\)
+
example: ridge coefficient paths vs \(\lambda\)
-
using the same prostate dataset, we will run ridge regressions with different values of \(\lambda\) and find the optimum \(\lambda\) value that minimizes test RSS
+
using the same prostate dataset, we will run ridge regressions with different values of \(\lambda\) and find the optimum \(\lambda\) value that minimizes test RSS
-
+
Regularized Regression - LASSO Regression
LASSO (least absolute shrinkage and selection operator) was introduced by Tibshirani (Journal of the Royal Statistical Society 1996)
similar to ridge, with slightly different penalty term
-
the penalized residual sum of squares (PRSS) takes the form of \[PRSS(\beta_j) = \sum_{i=1}^N \left(y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j \right)^2 + \lambda \sum_{j=1}^p |\beta_j|\]
-
this is equivalent to solving the equation \[PRSS(\beta_j) = \sum_{i=1}^N \left(y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j \right)^2\] subject to constraint \(\sum_{j=1}^p |\beta_j| \leq s\) where \(s\) is inversely proportional to \(\lambda\)
-
\(\lambda\) = tuning parameter
+
the penalized residual sum of squares (PRSS) takes the form of \[PRSS(\beta_j) = \sum_{i=1}^N \left(y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j \right)^2 + \lambda \sum_{j=1}^p |\beta_j|\]
+
this is equivalent to solving the equation \[PRSS(\beta_j) = \sum_{i=1}^N \left(y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j \right)^2\] subject to constraint \(\sum_{j=1}^p |\beta_j| \leq s\) where \(s\) is inversely proportional to \(\lambda\)
+
\(\lambda\) = tuning parameter
controls size of coefficients or the amount of regularization
-
large values of \(\lambda\) will set some coefficient equal to zero
+
large values of \(\lambda\) will set some coefficient equal to zero
Note: LASSO effectively performs model selection (choose subset of predictors) while shrinking other coefficients, where as Ridge only shrinks the coefficients
@@ -2323,7 +2400,7 @@
Regularized Regression - LASSO Regression
Note: the predictors are centered and scaled first before the regression is run
as.matrix(x) = the predictors must be in matrix/dataframe format
trace=TRUE = prints progress of the lasso regression
-
lasso$lambda = return the \(\lambda\)s used for each step of the lasso regression
+
lasso$lambda = return the \(\lambda\)s used for each step of the lasso regression
plot(lasso) = prints plot that shows the progression of the coefficients as they are set to zero one by one
predit.lars(lasso, test) = use the lasso model to predict on test data
@@ -2338,14 +2415,14 @@
Regularized Regression - LASSO Regression
[enet package] lasso<-enet(predictors, outcome, lambda = 0) = perform elastic net regression on given predictors and outcome
-
lambda=0 = default value for \(\lambda\)
+
lambda=0 = default value for \(\lambda\)
Note: lasso regression is a special case of elastic net regression, and forcing lambda=0 tells the function to fit a lasso regression
plot(lasso) = prints plot that shows the progression of the coefficients as they are set to zero one by one
predict.ent(lasso, test)= use the lasso model to predict on test data
-
[caret package] train(outcome ~ predictors, data=training, method="lasso") = perform lasso regression with given outcome and predictors
+
caret packagetrain(outcome ~ predictors, data=training, method="lasso") = perform lasso regression with given outcome and predictors
Note: outcome and predictors must be in the same dataframe
preProcess=c("center", "scale") = centers and scales the predictors before the model is built
@@ -2353,7 +2430,7 @@
Regularized Regression - LASSO Regression
Note: this is generally a good idea for building lasso regressions
-
[caret package] train(outcome~predictors,data=train,method="relaxo",lambda=5,phi=0.3) = perform relaxed lasso regression on given predictors and outcome
+
caret packagetrain(outcome~predictors,data=train,method="relaxo",lambda=5,phi=0.3) = perform relaxed lasso regression on given predictors and outcome
# set up data
inBuild <-createDataPartition(y=Wage$wage,p=0.7, list=FALSE)
validation <-Wage[-inBuild,]; buildData <-Wage[inBuild,]
inTrain <-createDataPartition(y=buildData$wage,p=0.7, list=FALSE)
@@ -2483,11 +2560,11 @@
Example - Model Ensembling
# validation data set RMSE Errorsvalidation =c(sqrt(sum((glm.pred.val-validation$wage)^2)),
sqrt(sum((rf.pred.val-validation$wage)^2)),
- sqrt(sum((comb.pred.val-validation$wage)^2))))
simple moving averages = prediction will be made for a time point by averaging together values from a number of prior periods \[ Y_{t}=\frac{1}{2*k+1}\sum_{j=-k}^k {y_{t+j}}\]
-
exponential smoothing/exponential moving average = weight time points that are closer to point of prediction than those that are further away \[\hat{y}_{t+1} = \alpha y_t + (1-\alpha)\hat{y}_{t-1}\]
+
simple moving averages = prediction will be made for a time point by averaging together values from a number of prior periods \[ Y_{t}=\frac{1}{2*k+1}\sum_{j=-k}^k {y_{t+j}}\]
+
exponential smoothing/exponential moving average = weight time points that are closer to point of prediction than those that are further away \[\hat{y}_{t+1} = \alpha y_t + (1-\alpha)\hat{y}_{t-1}\]
Note: many different methods of exponential smoothing are available, more information can be found here
@@ -2568,7 +2645,7 @@
R Commands and Examples
ts(googOpen, frequency=12) = convert data to a time series with frequency observations per time unit
-
frequency=12 = number of observations per unit time (12 in this case because there are 12 months in each year \(\rightarrow\) converts data into yearly time series)
+
frequency=12 = number of observations per unit time (12 in this case because there are 12 months in each year \(\rightarrow\) converts data into yearly time series)
decompose(ts) = decomposes time series into trend, seasonal, and irregular components by using moving averages
@@ -2607,25 +2684,25 @@
R Commands and Examples
quandl package is also used for finance-related predictions
example: decomposed time series
-
# load quantmod package
+
# load quantmod packagelibrary(quantmod);
# specify to and from dates
from.dat <-as.Date("01/01/00", format="%m/%d/%y")
to.dat <-as.Date("3/2/15", format="%m/%d/%y")
# get data for AAPL from Google Finance for the specified dates
-getSymbols("AAPL", src="google", from = from.dat, to = to.dat)
-
## [1] "AAPL"
-
# convert the retrieved daily data to monthly data
+getSymbols("AAPL", src="google", from = from.dat, to = to.dat)
+
[1] "AAPL"
+
# convert the retrieved daily data to monthly data
mAAPL <-to.monthly(AAPL)
# extract the closing price and convert it to yearly time series (12 observations per year)
ts <-ts(Cl(mAAPL), frequency =12)
# plot the decomposed parts of the time series
-plot(decompose(ts),xlab="Years")
-
+plot(decompose(ts),xlab="Years")
+
example: forecast
-
# load forecast library
+
# load forecast librarylibrary(forecast)
# find the number of rows (years)
rows <-ceiling(length(ts)/12)
@@ -2636,24 +2713,24 @@
R Commands and Examples
# plot the training setplot(ts.train)
# add the moving average in red
-lines(ma(ts.train,order=3),col="red")
-
-
# compute the exponential smoothing average
+lines(ma(ts.train,order=3),col="red")
+
+
# compute the exponential smoothing average
ets <-ets(ts.train,model="MMM")
# construct a forecasting model using the exponential smoothing function
fcast <-forecast(ets)
# plot forecast and add actual data in red
-plot(fcast); lines(ts.test,col="red")
-
-
# print the accuracy of the forecast model
-accuracy(fcast,ts.test)
-
## ME RMSE MAE MPE MAPE MASE
-## Training set 0.1188298 2.825883 1.646959 -0.61217 10.68901 0.1924329
-## Test set -7.8132889 16.736910 15.079222 -13.64900 20.31005 1.7618772
-## ACF1 Theil's U
-## Training set 0.09773823 NA
-## Test set 0.84664431 3.360515
-
\(\pagebreak\)
+plot(fcast); lines(ts.test,col="red")
+
+
# print the accuracy of the forecast model
+accuracy(fcast,ts.test)
+
ME RMSE MAE MPE MAPE MASE
+Training set 0.3738734 2.576067 1.438703 0.2377946 10.38942 0.1759433
+Test set 0.7292502 18.507031 15.347659 -3.6663982 19.01227 1.8769111
+ ACF1 Theil's U
+Training set 0.1026281 NA
+Test set 0.8661382 3.200942
+
\(\pagebreak\)
@@ -2700,7 +2777,7 @@
R Commands and Examples
cl_predict function in clue package provides similar functionality
-
# load iris data
+
# load iris datadata(iris)
# create training and test sets
inTrain <-createDataPartition(y=iris$Species,p=0.7, list=FALSE)
@@ -2715,8 +2792,8 @@
