Load required packages

library(kernlab)
library(mlbench)
library(caret)
library(lattice)
library(ggplot2)

7.2

  • After loading the dataset, convert the ‘x’ data from a matrix to a data frame
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)

  • creating a list with a vector ‘y’ and a matrix of predictors ‘x’. Also simulating a large test set to estimate the true error rate with good precision:
## 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:

Tuning Knn Model:

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.
  • Using function ‘postResample’ to get the test set perforamnce values
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.2040595 0.6819919 2.5683461

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

Answer

Now, we are going to build different models whose datasets are scaled and centered for model building.

SVM

Tuning a SVM model with method as svmRadial

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.545335  0.7804647  2.015121
##     0.50  2.319786  0.7965148  1.830009
##     1.00  2.188357  0.8119624  1.726031
##     2.00  2.103655  0.8241314  1.655842
##     4.00  2.066890  0.8294297  1.631062
##     8.00  2.052688  0.8313917  1.623563
##    16.00  2.049883  0.8318288  1.621842
##    32.00  2.049883  0.8318288  1.621842
##    64.00  2.049883  0.8318288  1.621842
##   128.00  2.049883  0.8318288  1.621842
## 
## Tuning parameter 'sigma' was held constant at a value of 0.06802164
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were sigma = 0.06802164 and C = 16.

Function used for resampling is postResample. After a model is tuned, the test set performance is calculated.

svmRadialPred <- predict(svmRadialModel, newdata = testData$x)
postResample(pred = svmRadialPred, obs = testData$y)
##      RMSE  Rsquared       MAE 
## 2.0864652 0.8236735 1.5854649

RMSE for SVM is 2.0864652.

Averaged Neural Networks

Now, let’s tune the model using averaged neural networks over decay and 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)

Displaying model results:

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.595816  0.7312476  2.014857
##   0.00   10    3.032906  0.6630178  2.322010
##   0.00   15    2.666682  0.7202745  2.135429
##   0.00   20    2.634394  0.7273102  2.116201
##   0.01    1    2.572011  0.7356944  1.998411
##   0.01   10    2.719744  0.7111143  2.174586
##   0.01   15    2.434612  0.7634988  1.934567
##   0.01   20    2.346213  0.7802151  1.850714
##   0.10    1    2.580129  0.7336991  2.000463
##   0.10   10    2.533933  0.7489387  2.004893
##   0.10   15    2.312874  0.7875826  1.826406
##   0.10   20    2.289472  0.7922099  1.799901
##   0.50    1    2.620984  0.7251650  2.034070
##   0.50   10    2.392990  0.7727536  1.896902
##   0.50   15    2.247981  0.7981173  1.778612
##   0.50   20    2.257949  0.7973880  1.768102
##   0.90    1    2.649164  0.7195326  2.057454
##   0.90   10    2.337911  0.7806326  1.846869
##   0.90   15    2.247166  0.7980863  1.773883
##   0.90   20    2.248641  0.7989799  1.770375
## 
## 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 = 15, decay = 0.9 and bag = FALSE.

Function used for resampling is postResample. After a model is tuned, the test set performance is calculated.

nnetPred <- predict(nnetModel, newdata = testData$x)
postResample(pred = nnetPred, obs = testData$y)
##     RMSE Rsquared      MAE 
## 1.894755 0.856176 1.441820

RMSE for averaged neural networks is 1.894755.

Mars

Using Mars method to tune our model by earth method:

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"))
## Loading required package: earth
## Loading required package: Formula
## Loading required package: plotmo
## Loading required package: plotrix
## Loading required package: TeachingDemos

Displaying results for model:

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.416233  0.2187486  3.630918
##   1        3      3.716956  0.4495842  2.981710
##   1        4      2.829383  0.6783492  2.285456
##   1        5      2.500824  0.7495068  2.002628
##   1        6      2.377484  0.7721684  1.897799
##   1        7      1.964112  0.8402449  1.540566
##   1        8      1.849850  0.8590078  1.446960
##   1        9      1.760287  0.8734284  1.379551
##   1       10      1.745434  0.8758361  1.358502
##   1       11      1.721573  0.8784527  1.333126
##   1       12      1.741003  0.8757510  1.341457
##   1       13      1.762562  0.8729306  1.355525
##   1       14      1.779852  0.8704986  1.376884
##   1       15      1.796118  0.8682289  1.386179
##   1       16      1.801970  0.8673854  1.392551
##   1       17      1.801970  0.8673854  1.392551
##   1       18      1.801970  0.8673854  1.392551
##   1       19      1.801970  0.8673854  1.392551
##   1       20      1.801970  0.8673854  1.392551
##   2        2      4.421087  0.2136318  3.622047
##   2        3      3.738888  0.4424544  3.004778
##   2        4      2.878704  0.6647839  2.315845
##   2        5      2.556082  0.7359640  2.039618
##   2        6      2.448842  0.7590378  1.941403
##   2        7      2.076809  0.8217323  1.631062
##   2        8      1.919449  0.8477097  1.506739
##   2        9      1.750995  0.8737507  1.383005
##   2       10      1.589003  0.8961730  1.262151
##   2       11      1.503925  0.9077370  1.182204
##   2       12      1.459184  0.9131975  1.139692
##   2       13      1.465062  0.9130346  1.140011
##   2       14      1.440347  0.9155789  1.123394
##   2       15      1.462849  0.9131866  1.141915
##   2       16      1.477622  0.9110707  1.148694
##   2       17      1.479198  0.9110804  1.152683
##   2       18      1.480900  0.9108744  1.151674
##   2       19      1.475717  0.9116473  1.146616
##   2       20      1.475717  0.9116473  1.146616
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nprune = 14 and degree = 2.

Using function postResample to resample our model:

marsPred <- predict(marsModel, newdata = testData$x)
postResample(pred = marsPred, obs = testData$y)
##      RMSE  Rsquared       MAE 
## 1.2779993 0.9338365 1.0147070

RMSE for MARS model is 1.2779993.

As per the RMSE of above three models, MARS model has lowest RMSE of 1.277 and it fits best compared to other two models.

Let’s display informative variables picked by MARS model using varImp since it has best fit and also lowest RMSE.

varImp(marsModel)
## earth variable importance
## 
##    Overall
## X1  100.00
## X4   75.40
## X2   49.00
## X5   15.72
## X3    0.00

the MARS model picked the most informative variables, X1 ~ X5.

7.5

Before beginning with analyses, let’s impute the missing data.

In the ChemicalManufacturingProcess dataset, missing data points are imputed using the bagImpute method.

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)

Splitting the dataset in to test and train

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]

In this question, I will train non-linear models and display their resamlped data together as a list.

Now, lets train KNN model

KNN

Tuning KNN model:

set.seed(1)
knnModel <- train(x = X.train,
                  y = y.train,
                  method = "knn",
                  preProc = c("center", "scale"),
                  tuneLength = 10)

Displaying results for KNN model

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.528693  0.3803651  1.197772
##    7  1.506720  0.4002915  1.177320
##    9  1.499989  0.4086278  1.177109
##   11  1.503919  0.4117534  1.184008
##   13  1.515831  0.4079424  1.200587
##   15  1.520598  0.4105664  1.203504
##   17  1.530349  0.4099642  1.211622
##   19  1.538437  0.4098386  1.222425
##   21  1.545944  0.4111556  1.227777
##   23  1.551160  0.4116032  1.231905
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 9.

averaged neural networks

Tuning averaged neural networks model:

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)

Displaying results for averaged neural networks model:

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.717312  0.28722203  1.386006
##   0.00    5     2.523979  0.29234481  1.976184
##   0.00   10    10.825893  0.07974973  6.674825
##   0.01    1     1.708552  0.36129221  1.359911
##   0.01    5     1.987131  0.31701606  1.504025
##   0.01   10     2.603818  0.25946573  1.918135
##   0.10    1     1.853687  0.35583003  1.418522
##   0.10    5     2.431821  0.25944924  1.614498
##   0.10   10     1.888270  0.33945946  1.449473
## 
## 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.01 and bag = FALSE.

MARS

Tuning MARS model:

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"))

Displaying results for MARS model

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.509191  0.4137720  1.156518
##   1        3      1.359584  0.5227552  1.075452
##   1        4      1.336746  0.5391225  1.063582
##   1        5      1.370408  0.5175297  1.090855
##   1        6      1.387837  0.5095898  1.101868
##   1        7      1.609820  0.4724114  1.149884
##   1        8      1.676364  0.4559593  1.179091
##   1        9      1.772344  0.4237835  1.217472
##   1       10      1.788229  0.4113308  1.235919
##   2        2      1.504918  0.4171505  1.153679
##   2        3      1.418399  0.4712292  1.105962
##   2        4      1.390204  0.5077052  1.090698
##   2        5      1.393000  0.5070131  1.097381
##   2        6      1.406442  0.5112263  1.110931
##   2        7      1.462187  0.4835682  1.141567
##   2        8      1.578470  0.4405677  1.173828
##   2        9      1.594006  0.4380953  1.187403
##   2       10      1.673527  0.4224645  1.215365
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nprune = 4 and degree = 1.

SVM

Tuning SVM model using radial bias method:

set.seed(1)
svmRadialModel <- train(x = X.train,
                        y = y.train,
                        method = "svmRadial",
                        tuneLength=10,
                        preProc = c("center", "scale"))

Displaying results for SVM model

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.546282  0.4487672  1.245100
##     0.50  1.441628  0.4982254  1.164037
##     1.00  1.353530  0.5392437  1.090121
##     2.00  1.302065  0.5617075  1.041127
##     4.00  1.277482  0.5706214  1.012269
##     8.00  1.269350  0.5729619  1.004311
##    16.00  1.268799  0.5732754  1.004297
##    32.00  1.268622  0.5733923  1.004236
##    64.00  1.268622  0.5733923  1.004236
##   128.00  1.268622  0.5733923  1.004236
## 
## Tuning parameter 'sigma' was held constant at a value of 0.0134875
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were sigma = 0.0134875 and C = 32.

(a)

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  1.0201113 1.1122607 1.1734457 1.177109 1.222392 1.481468    0
## NNet 1.0651549 1.2238542 1.3447942 1.359911 1.457554 1.724502    0
## MARS 0.8241305 0.9755207 1.0830901 1.063582 1.136216 1.236849    0
## SVM  0.8589540 0.9380778 0.9957457 1.004236 1.059893 1.246323    0
## 
## RMSE 
##          Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
## KNN  1.276010 1.409447 1.492596 1.499989 1.599320 1.867116    0
## NNet 1.270614 1.568322 1.672416 1.708552 1.841286 2.307736    0
## MARS 1.067359 1.244730 1.345101 1.336746 1.423537 1.573359    0
## SVM  1.047830 1.160957 1.271501 1.268622 1.346235 1.516128    0
## 
## Rsquared 
##            Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
## KNN  0.29995933 0.3559468 0.4052142 0.4086278 0.4516805 0.5196701    0
## NNet 0.06069223 0.2791625 0.3553430 0.3612922 0.4420046 0.5874767    0
## MARS 0.37514800 0.4907027 0.5350952 0.5391225 0.6119040 0.6668477    0
## SVM  0.34092856 0.5296469 0.5568335 0.5733923 0.6390080 0.7077512    0

By seeing the above RMSE, SVM model has lowest RMSE compared to others. Therefore, SVM model is optimal.

Now we will calculate test set performance:

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.1219756 0.5533245 0.9734028
## avNNet    1.2985400 0.5130544 1.0525351
## earth     1.1666837 0.5108776 0.9764355
## svmRadial 0.9862389 0.6739056 0.7958088

Test set performance models RMSE is lowest for SVM model.

Displaying SVM models results:

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

b

For fitted neural network model, there is no direct method to caculate most important predictor, the following list is computed from SVM model.

The varImp function with which we get the important variables of a model, it evaluates variable importance by fitting a loess smoother between the outcome and the predictors giving us R^2 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
## ManufacturingProcess13   82.21
## ManufacturingProcess36   79.27
## BiologicalMaterial06     75.61
## BiologicalMaterial03     71.87
## ManufacturingProcess17   70.62
## BiologicalMaterial12     66.86
## ManufacturingProcess09   62.20
## ManufacturingProcess06   55.39
## BiologicalMaterial02     53.61
## ManufacturingProcess31   46.69
## ManufacturingProcess33   45.77
## BiologicalMaterial11     42.39
## BiologicalMaterial04     39.70
## ManufacturingProcess29   36.70
## ManufacturingProcess11   36.39
## ManufacturingProcess12   35.94
## BiologicalMaterial08     31.86
## BiologicalMaterial09     30.98
## ManufacturingProcess30   30.12

As we see from the above results, out of first 10 predictors, 6 predictors are ManufacturingProcess predictors.

ManufacturingProcess32 is top ranked predictor. Therefore, it seems that ManufacturingProcess are more important predictors compares to BiologicalMaterial.

c

  • Exploring relationship between tp predictors and response for predictors.

  • Extract the support vectors determined by the model and sort the top features determined by loess R^2

vectorIndex <- SVindex(svmModel)
supportVectorX <- X.train[vectorIndex,]
supportVectorY <- y.train[vectorIndex]


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')
}

As we see, above scatter plots show relationships between predictors of the SVM model and the response variable Yield .

  • As we know there are 122 support vectors in the SVM model, I extract those and plotting these vectors with the top 10 variables ranked by the loess R^2 method.

  • We can infer from the plots, that top features have significant relationship with the target variable.

  • In ManufacturingProcess32 plot we can see that, as the variable value increases,Yield also increases