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=10sin(\pi { x }_{ 1 }{ x }_{ 2 })+20{ ({ x }_{ 3 }-0.5) }^{ 2 }+10{ x }_{ 4 }+5{ x }_{ 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: 

library(mlbench)
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: 

library(caret)
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.565620  0.4887976  2.886629
##    7  3.422420  0.5300524  2.752964
##    9  3.368072  0.5536927  2.715310
##   11  3.323010  0.5779056  2.669375
##   13  3.275835  0.6030846  2.628663
##   15  3.261864  0.6163510  2.621192
##   17  3.261973  0.6267032  2.616956
##   19  3.286299  0.6281075  2.640585
##   21  3.280950  0.6390386  2.643807
##   23  3.292397  0.6440392  2.656080
## 
## 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)?

All of the models built have the data centered and scaled. The resampling method used is bootstrapping. After a model is tuned, the test set performance is calculated.

Below, a SVM model with radial basis kernel is tuned over the cost penalty.

svmRadialModel <- train(x = trainingData$x,
                        y = trainingData$y,
                        method = "svmRadial",
                        tuneLength=10,
                        preProc = c("center", "scale"))
svmRadialModel
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 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:
## 
##   C       RMSE      Rsquared   MAE     
##     0.25  2.600536  0.7740370  2.070216
##     0.50  2.402691  0.7877362  1.887487
##     1.00  2.277967  0.8018302  1.776410
##     2.00  2.193637  0.8141532  1.712248
##     4.00  2.168900  0.8168553  1.691505
##     8.00  2.170608  0.8164106  1.696634
##    16.00  2.171541  0.8162344  1.696791
##    32.00  2.171541  0.8162344  1.696791
##    64.00  2.171541  0.8162344  1.696791
##   128.00  2.171541  0.8162344  1.696791
## 
## Tuning parameter 'sigma' was held constant at a value of 0.06651357
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were sigma = 0.06651357 and C = 4.
svmRadialPred <- predict(svmRadialModel, newdata = testData$x)
postResample(pred = svmRadialPred, obs = testData$y)
##      RMSE  Rsquared       MAE 
## 2.0682337 0.8274206 1.5638719

Below is a averaged neural networks model, tuned over the decay and the number of hidden units.

nnetGrid <- expand.grid(.decay=c(0, 0.01, 0.1, 0.5, 0.9),
                        .size=c(1, 10, 15, 20),
                        .bag=FALSE)

nnetModel <- train(x = trainingData$x,
                   y = trainingData$y,
                   method = "avNNet",
                   tuneGrid = nnetGrid,
                   preProc = c("center", "scale"),
                   trace=FALSE,
                   linout=TRUE,
                   maxit=500)

nnetModel
## Model Averaged Neural Network 
## 
## 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:
## 
##   decay  size  RMSE      Rsquared   MAE     
##   0.00    1    2.602828  0.7322740  2.039564
##   0.00   10    2.998200  0.6653139  2.309161
##   0.00   15    2.723540  0.7068081  2.158401
##   0.00   20    2.686683  0.7171313  2.120541
##   0.01    1    2.596094  0.7326926  2.024904
##   0.01   10    2.715628  0.7121922  2.150282
##   0.01   15    2.478809  0.7589825  1.963637
##   0.01   20    2.369387  0.7780794  1.842924
##   0.10    1    2.601462  0.7325583  2.022374
##   0.10   10    2.608817  0.7339390  2.055593
##   0.10   15    2.384686  0.7738628  1.875732
##   0.10   20    2.303036  0.7893723  1.805320
##   0.50    1    2.630215  0.7260825  2.048674
##   0.50   10    2.420944  0.7660550  1.889225
##   0.50   15    2.290734  0.7910127  1.783877
##   0.50   20    2.290799  0.7903527  1.777449
##   0.90    1    2.653411  0.7211391  2.071753
##   0.90   10    2.341611  0.7804838  1.815667
##   0.90   15    2.302331  0.7885456  1.793927
##   0.90   20    2.282793  0.7917056  1.772184
## 
## 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 = 20, decay = 0.9 and bag
##  = FALSE.
nnetPred <- predict(nnetModel, newdata = testData$x)
postResample(pred = nnetPred, obs = testData$y)
##      RMSE  Rsquared       MAE 
## 1.9858742 0.8409197 1.5164941

Below is a MARS model, tuned over the number of degree and number of terms to remove.

marsGrid <- expand.grid(.degree=1:2,
                        .nprune=2:20)

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

marsModel
## Multivariate Adaptive Regression Spline 
## 
## 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:
## 
##   degree  nprune  RMSE      Rsquared   MAE     
##   1        2      4.393073  0.2108510  3.635119
##   1        3      3.696621  0.4392012  2.992905
##   1        4      2.768210  0.6850359  2.226188
##   1        5      2.486185  0.7508486  1.975898
##   1        6      2.352806  0.7764305  1.860228
##   1        7      1.951355  0.8432353  1.534278
##   1        8      1.833062  0.8614528  1.449376
##   1        9      1.779539  0.8696732  1.403572
##   1       10      1.773260  0.8706526  1.386477
##   1       11      1.789521  0.8683454  1.403531
##   1       12      1.794487  0.8674361  1.405076
##   1       13      1.796928  0.8665304  1.415398
##   1       14      1.813537  0.8640154  1.417664
##   1       15      1.820343  0.8632130  1.421937
##   1       16      1.832905  0.8617081  1.431074
##   1       17      1.842618  0.8602311  1.437411
##   1       18      1.842618  0.8602311  1.437411
##   1       19      1.842618  0.8602311  1.437411
##   1       20      1.842618  0.8602311  1.437411
##   2        2      4.393073  0.2108510  3.635119
##   2        3      3.683862  0.4441170  2.984109
##   2        4      2.792686  0.6798224  2.241597
##   2        5      2.524880  0.7421610  2.004248
##   2        6      2.368873  0.7713308  1.861675
##   2        7      2.006290  0.8323134  1.584652
##   2        8      1.905067  0.8486125  1.496059
##   2        9      1.677959  0.8822861  1.322076
##   2       10      1.549804  0.9002228  1.228746
##   2       11      1.451818  0.9125002  1.158580
##   2       12      1.414729  0.9174028  1.122700
##   2       13      1.393863  0.9194045  1.107304
##   2       14      1.378216  0.9210424  1.088726
##   2       15      1.376000  0.9218160  1.083421
##   2       16      1.389363  0.9204868  1.092314
##   2       17      1.407142  0.9186494  1.103848
##   2       18      1.415173  0.9177465  1.108418
##   2       19      1.418349  0.9174046  1.110582
##   2       20      1.420031  0.9173037  1.111949
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nprune = 15 and degree = 2.
marsPred <- predict(marsModel, newdata = testData$x)
postResample(pred = marsPred, obs = testData$y)
##      RMSE  Rsquared       MAE 
## 1.1908806 0.9428866 0.9496858

It appears the MARS model is the best, with the lowest test set RMSE. Below the variable importances in the MARS model are calculated:

varImp(marsModel)
## earth variable importance
## 
##     Overall
## X1   100.00
## X4    85.07
## X2    69.06
## X5    48.95
## X3    39.50
## X8     0.00
## X7     0.00
## X6     0.00
## X10    0.00
## X9     0.00

As you can see, the MARS model picked the most informative variables, X1 ~ X5.

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.

The missing values in the ChemicalManufacturingProcess data are imputed using the bagImpute method. The train test set are splitted, with 20% of the data assigned to the test set.

library(AppliedPredictiveModeling)
data(ChemicalManufacturingProcess)

# Impute missing values using `bagImpulte`
(cmpImpute <- preProcess(ChemicalManufacturingProcess[,-c(1)], method=c('bagImpute')))
## Created from 152 samples and 57 variables
## 
## Pre-processing:
##   - bagged tree imputation (57)
##   - ignored (0)
cmp <- predict(cmpImpute, ChemicalManufacturingProcess[,-c(1)])

# Train/test plitting data, 20% testing
set.seed(1)
trainRow <- createDataPartition(ChemicalManufacturingProcess$Yield, p=0.8, list=FALSE)
X.train <- cmp[trainRow, ]
y.train <- ChemicalManufacturingProcess$Yield[trainRow]
X.test <- cmp[-trainRow, ]
y.test <- ChemicalManufacturingProcess$Yield[-trainRow]

Below, 4 nonlinear regression models are trained: KNN, averaged neural networks, MARS, SVM with radial basis kernel. The data are centered and scaled before training. The bootstrapped resampling method is used with 25 repetition.

KNN:

set.seed(1)
knnModel <- train(x = X.train,
                  y = y.train,
                  method = "knn",
                  preProc = c("center", "scale"),
                  tuneLength = 10)
knnModel
## k-Nearest Neighbors 
## 
## 144 samples
##  57 predictor
## 
## Pre-processing: centered (57), scaled (57) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 144, 144, 144, 144, 144, 144, ... 
## Resampling results across tuning parameters:
## 
##   k   RMSE      Rsquared   MAE     
##    5  1.490410  0.3751489  1.171802
##    7  1.452857  0.4035753  1.155866
##    9  1.442266  0.4161293  1.157868
##   11  1.433195  0.4293049  1.153951
##   13  1.444311  0.4193332  1.168129
##   15  1.439420  0.4295644  1.169162
##   17  1.445506  0.4291705  1.177224
##   19  1.452724  0.4282260  1.180844
##   21  1.455932  0.4293823  1.181112
##   23  1.467001  0.4233281  1.188719
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 11.

Averaged neural networks:

nnetGrid <- expand.grid(.decay=c(0, 0.01, 0.1),
                        .size=c(1, 5, 10),
                        .bag=FALSE)
set.seed(1)
nnetModel <- train(x = X.train,
                   y = y.train,
                   method = "avNNet",
                   tuneGrid = nnetGrid,
                   preProc = c("center", "scale"),
                   trace=FALSE,
                   linout=TRUE,
                   maxit=500)

nnetModel
## Model Averaged Neural Network 
## 
## 144 samples
##  57 predictor
## 
## Pre-processing: centered (57), scaled (57) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 144, 144, 144, 144, 144, 144, ... 
## Resampling results across tuning parameters:
## 
##   decay  size  RMSE       Rsquared    MAE     
##   0.00    1     1.642247  0.29873409  1.346165
##   0.00    5     2.777550  0.25160092  2.150528
##   0.00   10    11.602290  0.06743605  6.930266
##   0.01    1     1.696003  0.34540922  1.331958
##   0.01    5     2.144809  0.23060732  1.591667
##   0.01   10     2.547722  0.25807295  1.910411
##   0.10    1     1.875253  0.32218474  1.397362
##   0.10    5     2.352555  0.22997441  1.621383
##   0.10   10     1.956479  0.29067790  1.481998
## 
## 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 = 1, decay = 0 and bag
##  = FALSE.

MARS:

marsGrid <- expand.grid(.degree=1:2,
                        .nprune=2:10)
set.seed(1)
marsModel <- train(x = X.train,
                   y = y.train,
                   method = "earth",
                   tuneGrid = marsGrid,
                   preProc = c("center", "scale"))

marsModel
## Multivariate Adaptive Regression Spline 
## 
## 144 samples
##  57 predictor
## 
## Pre-processing: centered (57), scaled (57) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 144, 144, 144, 144, 144, 144, ... 
## Resampling results across tuning parameters:
## 
##   degree  nprune  RMSE      Rsquared   MAE     
##   1        2      1.400030  0.4459540  1.092897
##   1        3      1.349694  0.4834531  1.057376
##   1        4      1.351728  0.4943238  1.054085
##   1        5      1.414877  0.4583082  1.103210
##   1        6      1.399430  0.4713266  1.081606
##   1        7      1.649951  0.4250764  1.162053
##   1        8      2.251717  0.4191628  1.248845
##   1        9      2.594464  0.3819324  1.327671
##   1       10      2.310011  0.3871793  1.287845
##   2        2      1.408844  0.4415388  1.105878
##   2        3      1.341126  0.4927003  1.044363
##   2        4      1.394660  0.4752776  1.079283
##   2        5      1.397888  0.4815162  1.078715
##   2        6      1.462966  0.4524576  1.112856
##   2        7      1.530919  0.4343312  1.143171
##   2        8      1.540733  0.4325093  1.144439
##   2        9      1.583671  0.4267020  1.164235
##   2       10      3.464534  0.3846292  1.449434
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nprune = 3 and degree = 2.

SVM with radial basis:

set.seed(1)
svmRadialModel <- train(x = X.train,
                        y = y.train,
                        method = "svmRadial",
                        tuneLength=10,
                        preProc = c("center", "scale"))
svmRadialModel
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 144 samples
##  57 predictor
## 
## Pre-processing: centered (57), scaled (57) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 144, 144, 144, 144, 144, 144, ... 
## Resampling results across tuning parameters:
## 
##   C       RMSE      Rsquared   MAE     
##     0.25  1.461355  0.4604820  1.181075
##     0.50  1.379240  0.4875414  1.107229
##     1.00  1.328743  0.5064688  1.057170
##     2.00  1.303212  0.5189233  1.034093
##     4.00  1.296799  0.5220798  1.025622
##     8.00  1.298143  0.5195631  1.028264
##    16.00  1.298396  0.5192017  1.028455
##    32.00  1.298396  0.5192017  1.028455
##    64.00  1.298396  0.5192017  1.028455
##   128.00  1.298396  0.5192017  1.028455
## 
## Tuning parameter 'sigma' was held constant at a value of 0.01501522
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were sigma = 0.01501522 and C = 4.

(a) Which nonlinear regression model gives the optimal resampling and test set performance?

The resampling performance of all the models are calculated below:

resamp <- resamples(list(KNN=knnModel, NNet=nnetModel, MARS=marsModel, SVM=svmRadialModel))
summary(resamp)
## 
## Call:
## summary.resamples(object = resamp)
## 
## Models: KNN, NNet, MARS, SVM 
## Number of resamples: 25 
## 
## MAE 
##           Min.   1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## KNN  0.9392392 1.0746345 1.145797 1.153951 1.232707 1.353296    0
## NNet 1.2216962 1.2559464 1.343618 1.346165 1.425278 1.510521    0
## MARS 0.7546217 0.9792683 1.037590 1.044363 1.112475 1.251543    0
## SVM  0.8595576 0.9575932 1.030041 1.025622 1.119000 1.164606    0
## 
## RMSE 
##           Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## KNN  1.2162112 1.375834 1.419268 1.433195 1.510581 1.680709    0
## NNet 1.4693753 1.568498 1.639158 1.642247 1.730818 1.886740    0
## MARS 0.9390693 1.257275 1.322779 1.341126 1.428352 1.670686    0
## SVM  1.0890359 1.194449 1.327781 1.296799 1.389997 1.506162    0
## 
## Rsquared 
##           Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
## KNN  0.3068508 0.3633298 0.4222442 0.4293049 0.4800884 0.5711319    0
## NNet 0.1340076 0.2453872 0.3038915 0.2987341 0.3498206 0.5200922    0
## MARS 0.2753641 0.3905506 0.5118570 0.4927003 0.5677710 0.7242827    0
## SVM  0.3875883 0.4853334 0.5062911 0.5220798 0.5731451 0.6440086    0

Looking at the Mean of the RMSE metric, it appears the SVM model is optimal.

The test set performance is calculated below:

testPerf <- function(models, testData, testTarget) {
  method <- c()
  res <- data.frame()
  for(model in models){
    method <- c(method, model$method)
    pred <- predict(model, newdata=testData)
    res <- rbind(res, t(postResample(pred=pred, obs=testTarget)))
  }
  row.names(res) <- method
  return(res)
}

models <- list(knnModel, nnetModel, marsModel, svmRadialModel)

performance <- testPerf(models, X.test, y.test)
performance
##               RMSE  Rsquared       MAE
## knn       1.272998 0.5841402 1.0299432
## avNNet    1.409065 0.5734579 1.0933501
## earth     1.226947 0.5926094 1.0064580
## svmRadial 0.955912 0.7397831 0.7254831

The test set performance also suggests that the SVM model is the best, with the least RMSE and MAE, and highest R^2.

(svmModel <- svmRadialModel$finalModel)
## Support Vector Machine object of class "ksvm" 
## 
## SV type: eps-svr  (regression) 
##  parameter : epsilon = 0.1  cost C = 4 
## 
## Gaussian Radial Basis kernel function. 
##  Hyperparameter : sigma =  0.0150152233802607 
## 
## Number of Support Vectors : 129 
## 
## Objective Function Value : -86.019 
## Training error : 0.031149

(b) 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?

The varImp function does not have a model-specific method to calculate feature importance of a SVM radial basis function model. When calling varImp on a SVM RBF model, it uses filterVarImp to evaluate the variable importance by fitting a loess smoother between the outcome and the predictors, and using the R^2 statistic as the metric for variable importance.

topFeatures <- varImp(svmRadialModel)
topFeatures
## loess r-squared variable importance
## 
##   only 20 most important variables shown (out of 57)
## 
##                        Overall
## ManufacturingProcess32  100.00
## BiologicalMaterial06     84.75
## ManufacturingProcess36   73.66
## BiologicalMaterial03     70.69
## ManufacturingProcess13   68.69
## BiologicalMaterial02     64.16
## BiologicalMaterial12     60.11
## ManufacturingProcess17   57.19
## ManufacturingProcess09   55.42
## ManufacturingProcess33   55.39
## BiologicalMaterial04     51.35
## BiologicalMaterial11     46.74
## ManufacturingProcess29   46.47
## ManufacturingProcess06   46.38
## BiologicalMaterial01     41.19
## BiologicalMaterial08     38.75
## ManufacturingProcess31   37.22
## BiologicalMaterial09     26.01
## ManufacturingProcess11   25.25
## ManufacturingProcess15   24.31

6 out of the top 10 ranked predictors are ManufacturingProcess predictors. The top ranking predictor is ManufacturingProcess32. It appears that the ManufacturingProcess predictors are more important.

For the SVM model with radial basis function, there is no direct way of finding the variable importance. Only SVM with linear kernel can be evaluated with regards to variable importance.

(c) 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?

Below, I extract the 128 support vectors as determined by the SVM model, and then plot these vectors with the top 10 variables ranked by the loess R^2 method. These plots show that these top features do have noticeable relationship with the target variable. For example, you can see that as variable ManufacturingProcess32 increases, the Yield increases. However, under the context of SVM RBF, the plots do not reveal intuition about their relationship with yield. The SVM model with radial basis function transforms the original data into a higher dimensional space of transformed features. It is difficult to visualize how each feature contribute to the model in the higher dimensional space.

# Extract the support vectors determined by the model
vectorIndex <- SVindex(svmModel)
supportVectorX <- X.train[vectorIndex,]
supportVectorY <- y.train[vectorIndex]

# Sort the top features determined by loess R^2
topFeatures <- topFeatures$importance
topFeatures$predictor <- row.names(topFeatures)
topFeatures <- topFeatures[order(topFeatures$Overall, decreasing = T), ]
topFeatures <- row.names(topFeatures)


for (i in 1:10){
  plot(x=supportVectorX[, topFeatures[i]], y=supportVectorY, 
       xlab=topFeatures[i], ylab='Yield')
}