Libraries
library(mlbench)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.8
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.1.2 ✓ forcats 0.5.1## Warning: package 'tidyr' was built under R version 4.1.2## Warning: package 'readr' was built under R version 4.1.2## Warning: package 'dplyr' was built under R version 4.1.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftlibrary(earth)## Warning: package 'earth' was built under R version 4.1.2## Loading required package: Formula## Loading required package: plotmo## Warning: package 'plotmo' was built under R version 4.1.2## Loading required package: plotrix## Loading required package: TeachingDemosQuestion 7.2
Friedman (1991) introduced several benchmark data sets create by simulation. One of these simulations used the following nonlinear equation to create data: y = 10 sin(πx1x2) + 20(x3 − 0.5)2 + 10x4 + 5x5 + N(0, σ2) where the x values are random variables uniformly distributed between [0, 1] (there are also 5 other non-informative variables also created in the simulation). The package mlbench contains a function called mlbench.friedman1 that simulates these data
set.seed(200)
trainingData <- mlbench.friedman1(200, sd = 1)
##We convert the 'x' data from a matrix to a data frame
##One reason is that this will give the columns names.
trainingData$x <- data.frame(trainingData$x)
##Look at the data using
featurePlot(trainingData$x, trainingData$y)
#This creates a list with a vector 'y' and a matrix
##of predictors 'x'. Also simulate a large test set to
##estimate the true error rate with good precision:
testData <- mlbench.friedman1(5000, sd = 1)
testData$x <- data.frame(testData$x)Tune several models on this data:
# Model 1
knnModel <- train(
x = trainingData$x, y = trainingData$y, method = "knn", preProc = c('center','scale'), tuneLength=10)knnModel## k-Nearest Neighbors
##
## 200 samples
## 10 predictor
##
## Pre-processing: centered (10), scaled (10)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 200, 200, 200, 200, 200, 200, ...
## Resampling results across tuning parameters:
##
## k RMSE Rsquared MAE
## 5 3.466085 0.5121775 2.816838
## 7 3.349428 0.5452823 2.727410
## 9 3.264276 0.5785990 2.660026
## 11 3.214216 0.6024244 2.603767
## 13 3.196510 0.6176570 2.591935
## 15 3.184173 0.6305506 2.577482
## 17 3.183130 0.6425367 2.567787
## 19 3.198752 0.6483184 2.592683
## 21 3.188993 0.6611428 2.588787
## 23 3.200458 0.6638353 2.604529
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 17.knnPred <- predict(knnModel, newdata=testData$x)
## The function postResample can be used to get the test set performance values
postResample(pred=knnPred, obs=testData$y)## RMSE Rsquared MAE
## 3.2040595 0.6819919 2.5683461Tune Several Models on this data
Neural Net
The first model we will try is a Neural Network with averaging
nnetAvg <- avNNet(trainingData$x, trainingData$y,
size = 5,
decay = 0.01,
## Specify how many models to average
repeats = 5,
linout = TRUE,
## Reduce the amount of printed output
trace = FALSE,
## Expand the number of iterations to find
## parameter estimates..
maxit = 500)## Warning: executing %dopar% sequentially: no parallel backend registerednnetPred <- predict(nnetAvg, newdata=testData$x)
## The function postResample can be used to get the test set performance values
postResample(pred=nnetPred, obs=testData$y)## RMSE Rsquared MAE
## 1.5949271 0.8988612 1.2097831Using a neural network model, we can see some immediate gains in RSquared as well as a significant reductin in RSME.
MARS (Multivariate Adaptive Regression Spline)
marsFit <- earth(trainingData$x, trainingData$y)
summary(marsFit)## Call: earth(x=trainingData$x, y=trainingData$y)
##
## coefficients
## (Intercept) 18.451984
## h(0.621722-X1) -11.074396
## h(0.601063-X2) -10.744225
## h(X3-0.281766) 20.607853
## h(0.447442-X3) 17.880232
## h(X3-0.447442) -23.282007
## h(X3-0.636458) 15.150350
## h(0.734892-X4) -10.027487
## h(X4-0.734892) 9.092045
## h(0.850094-X5) -4.723407
## h(X5-0.850094) 10.832932
## h(X6-0.361791) -1.956821
##
## Selected 12 of 18 terms, and 6 of 10 predictors
## Termination condition: Reached nk 21
## Importance: X1, X4, X2, X5, X3, X6, X7-unused, X8-unused, X9-unused, ...
## Number of terms at each degree of interaction: 1 11 (additive model)
## GCV 2.540556 RSS 397.9654 GRSq 0.8968524 RSq 0.9183982marsPred <- predict(marsFit, newdata=testData$x)
## The function postResample can be used to get the test set performance values
postResample(pred=marsPred, obs=testData$y)## RMSE Rsquared MAE
## 1.8136467 0.8677298 1.3911836MARS performs significantly better than KNN, but fails to out-perfom the neural networks attempted earlier. We can also see from looking at the varImp() method below, that the MARS model did in fact select the most informative predictors:
varImp(marsFit)Question 7.5
library(AppliedPredictiveModeling)
data("ChemicalManufacturingProcess")# impute missing valuues
knn_impute_preproc <- preProcess(ChemicalManufacturingProcess, method = c("knnImpute"))
imputed <- predict(knn_impute_preproc, ChemicalManufacturingProcess)# Train test split
set.seed(42)
index <- createDataPartition(imputed$Yield, p=0.7, list=FALSE)
train <- imputed[index, ]
test <- imputed[-index, ]# Scaling
train_scale <- preProcess(train, method = c("center", "scale"))
test_scale <- preProcess(test, method = c("center", "scale"))## Warning in preProcess.default(test, method = c("center", "scale")): These
## variables have zero variances: BiologicalMaterial07train_prep <- predict(train_scale, train)
test_prep <- predict(test_scale, test)Knn Models
# Knn Model
knnModel <- train(
x = select(train_prep, -c('Yield')), y = train_prep$Yield, method = "knn", preProc = c('center','scale'), tuneLength=10)## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07
## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut =
## 10, : These variables have zero variances: BiologicalMaterial07knnModel## k-Nearest Neighbors
##
## 124 samples
## 57 predictor
##
## Pre-processing: centered (57), scaled (57)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 124, 124, 124, 124, 124, 124, ...
## Resampling results across tuning parameters:
##
## k RMSE Rsquared MAE
## 5 0.8192141 0.3402783 0.6496097
## 7 0.8122402 0.3544480 0.6526401
## 9 0.8083955 0.3653803 0.6532791
## 11 0.8086448 0.3710152 0.6532782
## 13 0.8011149 0.3918936 0.6466320
## 15 0.8046287 0.3939598 0.6508065
## 17 0.8068083 0.3959505 0.6525159
## 19 0.8087262 0.4012354 0.6526588
## 21 0.8120203 0.4042662 0.6576347
## 23 0.8155379 0.4072117 0.6612124
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 13.knnPred <- predict(knnModel, newdata = select(test_prep, -c('Yield')))
postResample(pred = knnPred, obs = test$Yield)## RMSE Rsquared MAE
## 0.6660585 0.5112143 0.5451521Neural Networks
nnetAvg <- avNNet(x = select(train_prep, -c('Yield')), y = train_prep$Yield,
size = 5,
decay = 0.01,
## Specify how many models to average
repeats = 5,
linout = TRUE,
## Reduce the amount of printed output
trace = FALSE,
## Expand the number of iterations to find
## parameter estimates..
maxit = 500)nnetPred <- predict(nnetAvg, newdata = select(test_prep, -c('Yield')))
postResample(pred = nnetPred, obs = test$Yield)## RMSE Rsquared MAE
## 0.6785891 0.5293388 0.5735806Interestingly, the neural network actually perfomed worse than the KNN model. Now let’s check out MARS.
MARS
marsFit <- earth(select(train_prep, -c('Yield')), train_prep$Yield)
marsPred <- predict(marsFit, newdata = select(test_prep, -c('Yield')))
postResample(pred = marsPred, obs = test_prep$Yield)## RMSE Rsquared MAE
## 0.6664640 0.5546325 0.5689533While only marginally, MARS performs the best out of the three models tried.
Questions
- Which nonlinear regression model gives the optimal resampling and test set performance?
Of the three non-linear algorithms that I tried, MARS provided the best results from resampling and test set.
- Which predictors are most important in the optimal nonlinear regression model? Do either the biological or process variables dominate the list? How do the top ten important predictors compare to the top ten predictors from the optimal linear model?)
# showing most important predictors in the MARS model
varImp(marsFit)We see that the most important predictors in the MARS model are also present as some of the top predictors in our non-linear model. However, not all are present. For one, we can see that in the non-linear model, none of the BiologicalMaterial flags appear as top predcitors.
Below are the results from last chapter:
ManufacturingProcess32 100.00000
ManufacturingProcess13 94.93985
BiologicalMaterial06 90.30530
ManufacturingProcess36 82.82351
ManufacturingProcess17 80.76723
ManufacturingProcess09 78.44090
BiologicalMaterial03 77.17445
BiologicalMaterial02 76.80806
BiologicalMaterial12 73.91389
ManufacturingProcess06 71.21710
- Explore the relationships between the top predictors and the response for the predictors that are unique to the optimal nonlinear regression model. Do these plots reveal intuition about the biological or process predictors and their relationship with yield?
ggplot(train_prep, aes(ManufacturingProcess32, Yield)) +
geom_point()
For ManufacturingProcess32 there appears to be a close-to-linear relationship with Yield. However, towards the end of the process there is a drop-off in returned yield.
ggplot(train_prep, aes(ManufacturingProcess09, Yield)) +
geom_point()
MP09 quickly creates yield
ggplot(train_prep, aes(ManufacturingProcess39, Yield)) +
geom_point()
Is this process boolean? It seems that Yield is either negative, or
positive depending on the status of the process here.
ggplot(train_prep, aes(ManufacturingProcess01, Yield)) +
geom_point()
This features is similar to the previous, but less pronounced.