DATA624: Homework 8

library(AppliedPredictiveModeling)
library(caret)
library(mlbench)
library(dplyr)
library(corrplot)

Objective

Do problems 7.2 and 7.5 in Kuhn and Johnson. There are only two but they have many parts. Please submit both a link to your Rpubs and the .rmd file.

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(\pi x_1 x_2) + 20(x_3 - 0.5)^2 + 10x_4 + 5x_5 + N(0, \sigma^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)

## or other methods.

## 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 these data. For example:

set.seed(200)
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.654912  0.4779838  2.958475
##    7  3.529432  0.5118581  2.861742
##    9  3.446330  0.5425096  2.780756
##   11  3.378049  0.5723793  2.719410
##   13  3.332339  0.5953773  2.692863
##   15  3.309235  0.6111389  2.663046
##   17  3.317408  0.6201421  2.678898
##   19  3.311667  0.6333800  2.682098
##   21  3.316340  0.6407537  2.688887
##   23  3.326040  0.6491480  2.705915
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 15.

knnPred <- predict(knnModel, newdata = testData$x)
## The function 'postResample' can be used to get the test set
## perforamnce values
postResample(pred = knnPred, obs = testData$y)

##      RMSE  Rsquared       MAE 
## 3.1750657 0.6785946 2.5443169

Which models appear to give the best performance? Does MARS select the informative predictors (those named X1–X5)?

Models

I developed different nonlinear algorithms to try to fit the data created above. I created KNN, SVM, MARS, and Neural Network models to compare testing performance and identify the best model for the dataset.

KNN

set.seed(200)
knnModel <- train(x = trainingData$x,
                  y = trainingData$y,
                  method = "knn",
                  preProc = c("center", "scale"),
                  tuneLength = 10,
                  trControl = trainControl(method = "cv"))
knnModel

## k-Nearest Neighbors 
## 
## 200 samples
##  10 predictor
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 180, 180, 180, 180, 180, 180, ... 
## Resampling results across tuning parameters:
## 
##   k   RMSE      Rsquared   MAE     
##    5  3.238598  0.5836232  2.705822
##    7  3.117335  0.6295372  2.561052
##    9  3.100423  0.6590940  2.524483
##   11  3.086639  0.6822198  2.506584
##   13  3.094904  0.6902613  2.504433
##   15  3.116059  0.7045172  2.516131
##   17  3.129874  0.7133067  2.529370
##   19  3.151840  0.7183283  2.546422
##   21  3.175787  0.7209301  2.574113
##   23  3.208213  0.7146199  2.611285
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 11.

knnPred <- predict(knnModel, newdata = testData$x)

knnMetrics <- postResample(pred = knnPred, obs = testData$y)

knnMetrics

##      RMSE  Rsquared       MAE 
## 3.1222641 0.6690472 2.4963650

SVM

set.seed(200)
svmModel <- train(x = trainingData$x,
                  y = trainingData$y,
                  method = "svmRadial",
                  preProc = c("center", "scale"),
                  tuneLength = 10,
                  trControl = trainControl(method = "cv"))

svmModel

## Support Vector Machines with Radial Basis Function Kernel 
## 
## 200 samples
##  10 predictor
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 180, 180, 180, 180, 180, 180, ... 
## Resampling results across tuning parameters:
## 
##   C       RMSE      Rsquared   MAE     
##     0.25  2.525164  0.7810576  2.010680
##     0.50  2.270567  0.7944850  1.794902
##     1.00  2.099356  0.8155574  1.659376
##     2.00  2.005858  0.8302852  1.578799
##     4.00  1.934650  0.8435677  1.528373
##     8.00  1.915665  0.8475605  1.528648
##    16.00  1.923914  0.8463074  1.535991
##    32.00  1.923914  0.8463074  1.535991
##    64.00  1.923914  0.8463074  1.535991
##   128.00  1.923914  0.8463074  1.535991
## 
## Tuning parameter 'sigma' was held constant at a value of 0.06299324
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were sigma = 0.06299324 and C = 8.

svmPred <- predict(svmModel, newdata = testData$x)

svmMetrics <- postResample(pred = svmPred, obs = testData$y)

svmMetrics

##      RMSE  Rsquared       MAE 
## 2.0541197 0.8290353 1.5586411

MARS

set.seed(200)
marsGrid <- expand.grid(.degree = 1:2, .nprune = 2:38)

marsModel <- train(x = trainingData$x,
                  y = trainingData$y,
                  method = "earth",
                  preProc = c("center", "scale"),
                  tuneGrid = marsGrid,
                  tuneLength = 10,
                  trControl = trainControl(method = "cv"))

## Loading required package: earth

## Loading required package: Formula

## Loading required package: plotmo

## Loading required package: plotrix

## Loading required package: TeachingDemos

marsModel

## Multivariate Adaptive Regression Spline 
## 
## 200 samples
##  10 predictor
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 180, 180, 180, 180, 180, 180, ... 
## Resampling results across tuning parameters:
## 
##   degree  nprune  RMSE      Rsquared   MAE      
##   1        2      4.188280  0.3042527  3.4606894
##   1        3      3.551182  0.4999832  2.8371162
##   1        4      2.653143  0.7167280  2.1282215
##   1        5      2.405769  0.7562160  1.9481607
##   1        6      2.295006  0.7754603  1.8531995
##   1        7      1.771950  0.8611767  1.3913569
##   1        8      1.647182  0.8774867  1.2995643
##   1        9      1.609816  0.8837307  1.2997049
##   1       10      1.635035  0.8798236  1.3094359
##   1       11      1.571915  0.8896147  1.2607106
##   1       12      1.571561  0.8898750  1.2530769
##   1       13      1.567577  0.8906927  1.2507948
##   1       14      1.571673  0.8909652  1.2455080
##   1       15      1.571673  0.8909652  1.2455080
##   1       16      1.571673  0.8909652  1.2455080
##   1       17      1.571673  0.8909652  1.2455080
##   1       18      1.571673  0.8909652  1.2455080
##   1       19      1.571673  0.8909652  1.2455080
##   1       20      1.571673  0.8909652  1.2455080
##   1       21      1.571673  0.8909652  1.2455080
##   1       22      1.571673  0.8909652  1.2455080
##   1       23      1.571673  0.8909652  1.2455080
##   1       24      1.571673  0.8909652  1.2455080
##   1       25      1.571673  0.8909652  1.2455080
##   1       26      1.571673  0.8909652  1.2455080
##   1       27      1.571673  0.8909652  1.2455080
##   1       28      1.571673  0.8909652  1.2455080
##   1       29      1.571673  0.8909652  1.2455080
##   1       30      1.571673  0.8909652  1.2455080
##   1       31      1.571673  0.8909652  1.2455080
##   1       32      1.571673  0.8909652  1.2455080
##   1       33      1.571673  0.8909652  1.2455080
##   1       34      1.571673  0.8909652  1.2455080
##   1       35      1.571673  0.8909652  1.2455080
##   1       36      1.571673  0.8909652  1.2455080
##   1       37      1.571673  0.8909652  1.2455080
##   1       38      1.571673  0.8909652  1.2455080
##   2        2      4.308195  0.2578397  3.6207010
##   2        3      3.706763  0.4604448  2.9757957
##   2        4      2.735073  0.6839796  2.2293125
##   2        5      2.492961  0.7285910  1.9701005
##   2        6      2.376226  0.7515868  1.8509493
##   2        7      2.056766  0.7880683  1.6268337
##   2        8      1.777202  0.8531962  1.3935838
##   2        9      1.689355  0.8660853  1.3332428
##   2       10      1.644239  0.8756652  1.3024389
##   2       11      1.515215  0.8908390  1.1976900
##   2       12      1.384404  0.9129573  1.1044863
##   2       13      1.315403  0.9237719  1.0512295
##   2       14      1.293663  0.9269840  1.0376277
##   2       15      1.244900  0.9320943  0.9875052
##   2       16      1.238511  0.9327014  0.9875981
##   2       17      1.244898  0.9321291  0.9943780
##   2       18      1.244898  0.9321291  0.9943780
##   2       19      1.244898  0.9321291  0.9943780
##   2       20      1.244898  0.9321291  0.9943780
##   2       21      1.244898  0.9321291  0.9943780
##   2       22      1.244898  0.9321291  0.9943780
##   2       23      1.244898  0.9321291  0.9943780
##   2       24      1.244898  0.9321291  0.9943780
##   2       25      1.244898  0.9321291  0.9943780
##   2       26      1.244898  0.9321291  0.9943780
##   2       27      1.244898  0.9321291  0.9943780
##   2       28      1.244898  0.9321291  0.9943780
##   2       29      1.244898  0.9321291  0.9943780
##   2       30      1.244898  0.9321291  0.9943780
##   2       31      1.244898  0.9321291  0.9943780
##   2       32      1.244898  0.9321291  0.9943780
##   2       33      1.244898  0.9321291  0.9943780
##   2       34      1.244898  0.9321291  0.9943780
##   2       35      1.244898  0.9321291  0.9943780
##   2       36      1.244898  0.9321291  0.9943780
##   2       37      1.244898  0.9321291  0.9943780
##   2       38      1.244898  0.9321291  0.9943780
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nprune = 16 and degree = 2.

marsPred <- predict(marsModel, newdata = testData$x)

marsMetrics <- postResample(pred = marsPred, obs = testData$y)

marsMetrics

##      RMSE  Rsquared       MAE 
## 1.2793868 0.9343367 1.0091132

Neural Network

set.seed(200)
nnetGrid <- expand.grid(.decay = c(0, 0.01, .1),
                        .size = c(1:10),
                        .bag = FALSE)

nnetModel <- train(trainingData$x,
                   trainingData$y,
                   method = "avNNet",
                   tuneGrid = nnetGrid,
                   trControl = trainControl(method = "cv"),
                   preProc = c("center", "scale"),
                   linout = TRUE,
                   trace = FALSE,
                   MaxNWts = 10 * (ncol(trainingData$x) + 1) + 10 + 1,
                   maxit = 500)

## Warning: executing %dopar% sequentially: no parallel backend registered

nnetModel

## Model Averaged Neural Network 
## 
## 200 samples
##  10 predictor
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 180, 180, 180, 180, 180, 180, ... 
## Resampling results across tuning parameters:
## 
##   decay  size  RMSE      Rsquared   MAE     
##   0.00    1    2.409283  0.7621183  1.899438
##   0.00    2    2.422970  0.7596059  1.940174
##   0.00    3    2.050680  0.8164736  1.639116
##   0.00    4    1.945570  0.8359932  1.553962
##   0.00    5    2.459608  0.7843926  1.852817
##   0.00    6    3.619609  0.6518455  2.367746
##   0.00    7    3.732346  0.5465485  2.505487
##   0.00    8    5.341921  0.4925666  3.027228
##   0.00    9    4.448181  0.5628776  2.637284
##   0.00   10    3.736545  0.6189292  2.529286
##   0.01    1    2.380838  0.7641924  1.871203
##   0.01    2    2.456920  0.7487966  1.925584
##   0.01    3    2.152614  0.8037282  1.690702
##   0.01    4    1.926277  0.8453345  1.547266
##   0.01    5    2.129094  0.8103702  1.698914
##   0.01    6    2.140650  0.8117168  1.698805
##   0.01    7    2.414589  0.7646608  1.911319
##   0.01    8    2.366556  0.7741779  1.873354
##   0.01    9    2.368192  0.7641654  1.781225
##   0.01   10    2.336267  0.7855705  1.857180
##   0.10    1    2.392300  0.7614548  1.873846
##   0.10    2    2.437038  0.7557138  1.918843
##   0.10    3    2.136582  0.8043180  1.702665
##   0.10    4    2.009698  0.8245209  1.574401
##   0.10    5    2.015296  0.8345946  1.586740
##   0.10    6    2.038841  0.8283220  1.586032
##   0.10    7    2.129954  0.8133844  1.707177
##   0.10    8    2.148353  0.8099883  1.690601
##   0.10    9    2.254190  0.7942883  1.759231
##   0.10   10    2.359693  0.7719593  1.873586
## 
## Tuning parameter 'bag' was held constant at a value of FALSE
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were size = 4, decay = 0.01 and bag = FALSE.

nnetPred <- predict(nnetModel, newdata = testData$x)

nnetMetrics <- postResample(pred = nnetPred, obs = testData$y)

nnetMetrics

##      RMSE  Rsquared       MAE 
## 2.0603901 0.8320669 1.5289876

Best Model

The best model across all metrics (RMSE, Rsquared, and MAE) is the Multivariate Adaptive Regression Spline model. The neural network and support vector machine models performed slightly worse but were comparable. The worst model for this dataset is the KNN model.

kableExtra::kable(rbind(knnMetrics,
      svmMetrics,
      marsMetrics,
      nnetMetrics))

	RMSE	Rsquared	MAE
knnMetrics	3.122264	0.6690472	2.496365
svmMetrics	2.054120	0.8290353	1.558641
marsMetrics	1.279387	0.9343367	1.009113
nnetMetrics	2.060390	0.8320669	1.528988

MARS Variable Selection

As shown by the variable importance function, X1 - X5 were selected by the MARS algorithm.

varImp(marsModel)

## earth variable importance
## 
##    Overall
## X1  100.00
## X4   75.25
## X2   48.76
## X5   15.54
## X3    0.00

7.5

Exercise 6.3 describes data for a chemical manufacturing process. Use the same data imputation, data splitting, and pre-processing steps as before and train several nonlinear regression models.

Which nonlinear regression model gives the optimal resampling and test set performance?
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?
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?

Data Preparation

data("ChemicalManufacturingProcess")

imputed.knn <- preProcess(ChemicalManufacturingProcess,
           method = "knnImpute",
           k = sqrt(nrow(ChemicalManufacturingProcess))
           )

imputed.data <- predict(imputed.knn, ChemicalManufacturingProcess)

near_zero <- nearZeroVar(imputed.data)

imputed.data <- imputed.data[, -near_zero]

set.seed(115)
train_test_split <- createDataPartition(ChemicalManufacturingProcess$Yield, p = 0.7, list = FALSE)

train.data <- imputed.data[train_test_split,]
test.data <- imputed.data[-train_test_split,]

Part A

KNN

set.seed(200)
knnModel <- train(Yield ~ .,
                   data = train.data,
                  method = "knn",
                  preProc = c("center", "scale"),
                  tuneLength = 10,
                  trControl = trainControl(method = "cv"))

knnPred <- predict(knnModel, newdata = test.data)

knnMetrics <- postResample(pred = knnPred, obs = test.data$Yield)

knnMetrics

##      RMSE  Rsquared       MAE 
## 0.6876584 0.3738299 0.5525322

SVM

set.seed(200)
svmModel <- train(Yield ~ .,
                   data = train.data,
                  method = "svmRadial",
                  preProc = c("center", "scale"),
                  tuneLength = 10,
                  trControl = trainControl(method = "cv"))

svmPred <- predict(svmModel, newdata = test.data)

svmMetrics <- postResample(pred = svmPred, obs = test.data$Yield)

svmMetrics

##      RMSE  Rsquared       MAE 
## 0.5882593 0.5430571 0.4684853

MARS

set.seed(200)
marsGrid <- expand.grid(.degree = 1:2, .nprune = 2:38)

marsModel <- train(Yield ~ .,
                   data = train.data,
                  method = "earth",
                  preProc = c("center", "scale"),
                  tuneGrid = marsGrid,
                  tuneLength = 10,
                  trControl = trainControl(method = "cv"))

marsPred <- predict(marsModel, newdata = test.data)

marsMetrics <- postResample(pred = marsPred, obs = test.data$Yield)

marsMetrics

##      RMSE  Rsquared       MAE 
## 0.6335728 0.5732152 0.5112639

Neural Network

set.seed(200)

nnetModel <- train(Yield ~ .,
                   data = train.data,
                   method = "avNNet",
                   trControl = trainControl(method = "cv"),
                   preProc = c("center", "scale"),
                   linout = TRUE,
                   trace = FALSE,
                   MaxNWts = 10 * (ncol(trainingData$x) + 1) + 10 + 1,
                   maxit = 500)

nnetPred <- predict(nnetModel, newdata = test.data)

nnetMetrics <- postResample(pred = nnetPred, obs = test.data$Yield)

nnetMetrics

##      RMSE  Rsquared       MAE 
## 0.5851825 0.6648197 0.4962526

Model Comparison

The SVM and Neural Network models produce similar results, though the Neural Network model appears to be a slightly superior fit for the data as it is better in terms of both R squared and RMSE and only slightly inferior in MAE.

kableExtra::kable(rbind(knnMetrics,
      svmMetrics,
      marsMetrics,
      nnetMetrics))

	RMSE	Rsquared	MAE
knnMetrics	0.6876584	0.3738299	0.5525322
svmMetrics	0.5882593	0.5430571	0.4684853
marsMetrics	0.6335728	0.5732152	0.5112639
nnetMetrics	0.5851825	0.6648197	0.4962526

Part B

Both the Nonlinear and Linear optimal models appear to weight the importance of each variable similarly.

Nonlinear: Neural Network Variable Importance

There are a combination of ManufacturingProcess and BiologicalMaterial variables that are within the top 10. 4 of the top 6 are BiologicalMaterial though the most important variable is related to the ManufacturingProcess.

varImp(nnetModel)

## loess r-squared variable importance
## 
##   only 20 most important variables shown (out of 56)
## 
##                        Overall
## ManufacturingProcess32  100.00
## BiologicalMaterial06     81.83
## ManufacturingProcess13   80.71
## BiologicalMaterial12     76.99
## BiologicalMaterial03     76.27
## BiologicalMaterial02     71.31
## ManufacturingProcess17   66.10
## ManufacturingProcess36   63.06
## ManufacturingProcess31   56.27
## ManufacturingProcess09   52.94
## BiologicalMaterial11     51.43
## BiologicalMaterial04     48.41
## ManufacturingProcess29   47.09
## ManufacturingProcess33   44.02
## BiologicalMaterial09     43.87
## BiologicalMaterial08     43.36
## ManufacturingProcess06   43.23
## ManufacturingProcess02   39.59
## BiologicalMaterial01     38.21
## ManufacturingProcess15   35.47

Linear: Ridge Regression Variable Importance

The top 5 variables are a combination of biological and process variables. The top 20 contains 11 process variables and 9 biological variables. It appears that there is a fairly equal split of importance between the two variable types.

ridgeGrid <- data.frame(.lambda = seq(0.001, 1.1, length = 20))

ridgeTune <- train(Yield~., train.data,
                   method = "ridge",
                   metric = "RMSE",
                   tuneGrid = ridgeGrid,
                   tuneLength = 20,
                   trControl = trainControl(method = "cv", number = 10),
                   preProc = c("center", "scale"))

varImp(ridgeTune)

## loess r-squared variable importance
## 
##   only 20 most important variables shown (out of 56)
## 
##                        Overall
## ManufacturingProcess32  100.00
## BiologicalMaterial06     81.83
## ManufacturingProcess13   80.71
## BiologicalMaterial12     76.99
## BiologicalMaterial03     76.27
## BiologicalMaterial02     71.31
## ManufacturingProcess17   66.10
## ManufacturingProcess36   63.06
## ManufacturingProcess31   56.27
## ManufacturingProcess09   52.94
## BiologicalMaterial11     51.43
## BiologicalMaterial04     48.41
## ManufacturingProcess29   47.09
## ManufacturingProcess33   44.02
## BiologicalMaterial09     43.87
## BiologicalMaterial08     43.36
## ManufacturingProcess06   43.23
## ManufacturingProcess02   39.59
## BiologicalMaterial01     38.21
## ManufacturingProcess15   35.47

Part C

It appears that the ManufacturingProcess variables have a higher correlation coefficent with Yield compared to the top BiologicalMaterial variables. To increase Yield, a combination of postively correlated variables like ManufacturingProcess32 and BiologicalMaterial02 should be used whereas decreasing ManufacturingProcess36 and ManufacturingProcess13 would prove to be beneficial.

important.vars <- imputed.data %>%
  select("Yield", "ManufacturingProcess32", "BiologicalMaterial06", "ManufacturingProcess13",
         "BiologicalMaterial12", "BiologicalMaterial03", "BiologicalMaterial02", "ManufacturingProcess17", "ManufacturingProcess36",
         "ManufacturingProcess09", "ManufacturingProcess31")

cor.matrix <- cor(important.vars)

corrplot(corr = cor.matrix, tl.col = 'black',  type = 'lower', diag = FALSE)