library(kernlab)
library(mlbench)
library(caret)
library(lattice)
library(ggplot2)
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:
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.
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.
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