Assignment 8
Jatin Jain
7.2. Friedman (1991) introduced several benchmark data sets create by sim- ulation. One of these simulations used the following nonlinear equation to create data:
y = 10sin(π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 simula- tion). 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:
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.5683461library(AppliedPredictiveModeling)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 here 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.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.svmRadialPred <- predict(svmRadialModel, newdata = testData$x)
postResample(pred = svmRadialPred, obs = testData$y)## RMSE Rsquared MAE
## 2.0864652 0.8236735 1.5854649Below is an 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)## Warning: executing %dopar% sequentially: no parallel backend registerednnetModel## 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.nnetPred <- predict(nnetModel, newdata = testData$x)
postResample(pred = nnetPred, obs = testData$y)## RMSE Rsquared MAE
## 1.894755 0.856176 1.441820Below 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"))## Loading required package: earth## Loading required package: Formula## Loading required package: plotmo## Loading required package: plotrix## Loading required package: TeachingDemosmarsModel## 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.marsPred <- predict(marsModel, newdata = testData$x)
postResample(pred = marsPred, obs = testData$y)## RMSE Rsquared MAE
## 1.2779993 0.9338365 1.0147070It appears that the MARS model is the best, with the lowest test set RMSE. The variable importance in the MARS model are calculated:
varImp(marsModel)## earth variable importance
##
## Overall
## X1 100.00
## X4 75.40
## X2 49.00
## X5 15.72
## X3 0.00As 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.
data(ChemicalManufacturingProcess)help(ChemicalManufacturingProcess)(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)])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.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:
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.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:
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.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 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.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. 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 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 0Looking at the Mean of the RMSE metric, it appears that 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.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.7958088The 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 = 32
##
## Gaussian Radial Basis kernel function.
## Hyperparameter : sigma = 0.0134874993395927
##
## Number of Support Vectors : 123
##
## Objective Function Value : -98.4469
## Training error : 0.009062b. 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
## 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.126 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 extracted the 128 support vectors as determined by the SVM model, and then plotted those 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.
library(kernlab)##
## Attaching package: 'kernlab'## The following object is masked from 'package:ggplot2':
##
## alphavectorIndex <- 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')
}
