Homework 8
Seung Min Song
2024-04-07
7.2 and 7.5 in Kuhn and Johnson
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(π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 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
KNN
set.seed(123)
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.548433 0.4919564 2.888283
## 7 3.425531 0.5255725 2.778090
## 9 3.346026 0.5523023 2.704791
## 11 3.252313 0.5875603 2.620492
## 13 3.232552 0.6000482 2.601113
## 15 3.205067 0.6203296 2.586704
## 17 3.172791 0.6408339 2.566738
## 19 3.183306 0.6494300 2.587220
## 21 3.190873 0.6556293 2.596793
## 23 3.202234 0.6597746 2.604279
##
## 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.5683461SVM
set.seed(123)
svmModel <- train(x = trainingData$x,
y = trainingData$y,
method = "svmLinear",
preProc = c("center", "scale"),
tuneLength = 10)
trControl = trainControl(method = "cv", number = 10)
print(svmModel)## Support Vector Machines with Linear 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:
##
## RMSE Rsquared MAE
## 2.512719 0.7487788 2.009062
##
## Tuning parameter 'C' was held constant at a value of 1svmPred <- predict(svmModel, newdata = testData$x)
svmPerformance <- postResample(pred = svmPred, obs = testData$y)
print(svmPerformance)## RMSE Rsquared MAE
## 2.7633860 0.6973384 2.0970616MARS
set.seed(123)
trControl = trainControl(method = "cv", number = 10)
marsModel <- train(x = trainingData$x, y = trainingData$y,
method = "earth",
trControl = trainControl(method = "cv", number = 10))
print(marsModel)## Multivariate Adaptive Regression Spline
##
## 200 samples
## 10 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 180, 180, 180, 180, 180, 180, ...
## Resampling results across tuning parameters:
##
## nprune RMSE Rsquared MAE
## 2 4.311247 0.2748122 3.603533
## 8 1.710582 0.8857656 1.323419
## 15 1.692444 0.8801290 1.317804
##
## Tuning parameter 'degree' was held constant at a value of 1
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nprune = 15 and degree = 1.summary(marsModel$finalModel)## Call: earth(x=data.frame[200,10], y=c(18.46,16.1,17...), keepxy=TRUE, degree=1,
## nprune=15)
##
## coefficients
## (Intercept) 18.451984
## h(0.621722-X1) -11.074396
## h(0.601063-X2) -10.744225
## h(X3-0.281766) 20.607853
## h(0.447442-X3) 17.880232
## h(X3-0.447442) -23.282007
## h(X3-0.636458) 15.150350
## h(0.734892-X4) -10.027487
## h(X4-0.734892) 9.092045
## h(0.850094-X5) -4.723407
## h(X5-0.850094) 10.832932
## h(X6-0.361791) -1.956821
##
## Selected 12 of 18 terms, and 6 of 10 predictors (nprune=15)
## Termination condition: Reached nk 21
## Importance: X1, X4, X2, X5, X3, X6, X7-unused, X8-unused, X9-unused, ...
## Number of terms at each degree of interaction: 1 11 (additive model)
## GCV 2.540556 RSS 397.9654 GRSq 0.8968524 RSq 0.9183982marsPred <- predict(marsModel, newdata = testData$x)
marsPerformance <- postResample(pred = marsPred, obs = testData$y)
print(marsPerformance)## RMSE Rsquared MAE
## 1.8136467 0.8677298 1.3911836Random Forest
numVars <- ncol(trainingData$x)
mtryValue <- floor(numVars / 3)
tuneGrid = expand.grid(.mtry = mtryValue)
set.seed(123)
rfModel <- train(x = trainingData$x,
y = trainingData$y,
method = "rf",
preProc = c("center", "scale"),
tuneLength = 10,
ntree = 100,
trControl = trainControl(method = "cv"),
tuneGrid = tuneGrid)
print(rfModel)## Random Forest
##
## 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:
##
## RMSE Rsquared MAE
## 2.673419 0.7877664 2.235552
##
## Tuning parameter 'mtry' was held constant at a value of 3rfPred <- predict(rfModel, newdata = testData$x)
rfPerformance <- postResample(pred = rfPred, obs = testData$y)
print(rfPerformance)## RMSE Rsquared MAE
## 2.650427 0.803300 2.120929Which models appear to give the best performance?
As a result of the comparison, the MARS model showed the best performance in all aspects of RMSE, R-squared, and MAE. It has the lowest RMSE and MAE and the highest R-squared, making it the best predictor of complex nonlinear relationships in the data.
The next best performing model is Random Forest. In particular, R-squared is high, explaining the data’s volatility relatively well. However, the RMSE and MAE are higher than MARS, so the prediction accuracy is somewhat lower.
Does MARS select the informative predictors (those named X1–X5)?
The MARS model selected highly informative predictor variables, X1 to X5. As you can see in the output of summary(marsModel$finalModel), the model selected a total of 6 predictors (including X6), including X1, X2, X3, X4, and X5. Here, you can see that one or more terms for each variable from X1 to X5 are included in the model. This means that MARS considered these variables important in constructing the model.
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.
library(AppliedPredictiveModeling)
data(ChemicalManufacturingProcess)
ChemicalManufacturingProcess <- data.frame(ChemicalManufacturingProcess)
dim(ChemicalManufacturingProcess)## [1] 176 58(a)
Which nonlinear regression model gives the optimal resampling and test set performance?
\(KNN\)
library(mlbench)
set.seed(123)
splitIndex <- createDataPartition(ChemicalManufacturingProcess$Yield, p = 0.7, list = FALSE)
trainingData <- ChemicalManufacturingProcess[splitIndex, ]
testData <- ChemicalManufacturingProcess[-splitIndex, ]
set.seed(123)
trainingData <- na.omit(trainingData)
testData <- na.omit(testData)
set.seed(123)
knnModel <- train(Yield ~ ., data = trainingData,
method = "knn",
preProc = c("center", "scale"),
tuneLength = 10)
print(knnModel)## k-Nearest Neighbors
##
## 111 samples
## 57 predictor
##
## Pre-processing: centered (57), scaled (57)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 111, 111, 111, 111, 111, 111, ...
## Resampling results across tuning parameters:
##
## k RMSE Rsquared MAE
## 5 1.632973 0.3506610 1.305092
## 7 1.634430 0.3393930 1.320684
## 9 1.623002 0.3477667 1.323579
## 11 1.602604 0.3658009 1.302110
## 13 1.595492 0.3753677 1.294506
## 15 1.602211 0.3766015 1.293872
## 17 1.605805 0.3806470 1.289912
## 19 1.608939 0.3827905 1.290204
## 21 1.612810 0.3905769 1.296412
## 23 1.618294 0.3911850 1.301567
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 13.knnPred <- predict(knnModel, newdata = testData)
knnPerformance <- postResample(pred = knnPred, obs = testData$Yield)
print(knnPerformance)## RMSE Rsquared MAE
## 1.1373193 0.5518414 0.9593058\(SVM\)
set.seed(123)
trainingData <- na.omit(trainingData)
testData <- na.omit(testData)
svmModel <- train(Yield ~ ., data = trainingData,
method = "svmLinear",
preProc = c("center", "scale"),
tuneLength = 10)
trControl = trainControl(method = "cv", number = 10)
print(svmModel)## Support Vector Machines with Linear Kernel
##
## 111 samples
## 57 predictor
##
## Pre-processing: centered (57), scaled (57)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 111, 111, 111, 111, 111, 111, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 6.098155 0.2457446 2.413358
##
## Tuning parameter 'C' was held constant at a value of 1svmPred <- predict(svmModel, newdata = testData)
svmPerformance <- postResample(pred = svmPred, obs = testData$Yield)
print(svmPerformance)## RMSE Rsquared MAE
## 1.5345023 0.4250666 1.2058334\(MARS\)
set.seed(123)
trainingData <- na.omit(trainingData)
testData <- na.omit(testData)
trControl = trainControl(method = "cv", number = 10)
marsModel <- train(Yield ~ ., data = trainingData,
method = "earth",
trControl = trainControl(method = "cv", number = 10))
print(marsModel)## Multivariate Adaptive Regression Spline
##
## 111 samples
## 57 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 100, 99, 99, 101, 99, 100, ...
## Resampling results across tuning parameters:
##
## nprune RMSE Rsquared MAE
## 2 1.530251 0.4451783 1.244308
## 9 1.386700 0.5639891 1.113800
## 17 1.435022 0.5626567 1.150932
##
## Tuning parameter 'degree' was held constant at a value of 1
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nprune = 9 and degree = 1.marsPred <- predict(marsModel, newdata = testData)
marsPerformance <- postResample(pred = marsPred, obs = testData$Yield)
print(marsPerformance)## RMSE Rsquared MAE
## 1.2960534 0.5335829 1.0363911\(Random Forest\)
set.seed(123)
trainingData <- na.omit(trainingData)
testData <- na.omit(testData)
numVars <- ncol(trainingData)
mtryValue <- floor(numVars / 3)
tuneGrid = expand.grid(.mtry = mtryValue)
rfModel <- train(Yield ~ ., data = trainingData,
method = "rf",
preProc = c("center", "scale"),
tuneLength = 10,
ntree = 100,
trControl = trainControl(method = "cv"),
tuneGrid = tuneGrid)
print(rfModel)## Random Forest
##
## 111 samples
## 57 predictor
##
## Pre-processing: centered (57), scaled (57)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 100, 99, 99, 101, 99, 100, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 1.165449 0.6862737 0.907214
##
## Tuning parameter 'mtry' was held constant at a value of 19rfPred <- predict(rfModel, newdata = testData)
rfPerformance <- postResample(pred = rfPred, obs = testData$Yield)
print(rfPerformance)## RMSE Rsquared MAE
## 0.9266758 0.7378862 0.7270000When comparing the performance of each model, the Random Forest model showed the best results in all performance metrics (RMSE, R-squared, MAE) on the test set. In particular, it has the lowest RMSE and the highest R-squared value, indicating that it best explains the variability of the data and has high prediction accuracy.
(b)
Which predictors are most important in the optimal nonlinear regression model?
The most important predictors in the model are ManufacturingProcess13, ManufacturingProcess32, ManufacturingProcess09, BiologicalMaterial03, and ManufacturingProcess31.
Do either the biological or process variables dominate the list?
The list shows a slightly higher number of ManufacturingProcess variables (12) compared to biological variables (8) among the top 20 most important predictors identified by the random forest model.
How do the top ten important predictors compare to the top ten predictors from the optimal linear model?
\(Commonly important variables:\) Both models selected ManufacturingProcess32, ManufacturingProcess09, ManufacturingProcess17, and BiologicalMaterial06 as the top 10 important variables. This indicates that these variables play a particularly important role in the model.
\(Features of the random forest model:\) The random forest model selected ManufacturingProcess13, BiologicalMaterial03, ManufacturingProcess31, BiologicalMaterial11, BiologicalMaterial12, and BiologicalMaterial02 as important variables. This may reflect that random forests are particularly useful for capturing nonlinearities and complex interactions in data.
\(Features unique to the Elastic Net model:\) The Elastic Net model selected ManufacturingProcess06, ManufacturingProcess37, BiologicalMaterial05, ManufacturingProcess45, ManufacturingProcess04, and ManufacturingProcess34 as important variables. As a linear model, Elastic Net likely selects important variables based on linear relationships between the variables.
importance <- varImp(rfModel, scale = FALSE)
importance## rf variable importance
##
## only 20 most important variables shown (out of 57)
##
## Overall
## ManufacturingProcess13 71.516
## ManufacturingProcess32 55.618
## ManufacturingProcess09 27.878
## BiologicalMaterial03 21.560
## ManufacturingProcess31 17.775
## ManufacturingProcess17 16.778
## BiologicalMaterial06 16.553
## BiologicalMaterial11 10.207
## BiologicalMaterial12 9.515
## BiologicalMaterial02 9.116
## BiologicalMaterial05 7.506
## ManufacturingProcess36 7.309
## ManufacturingProcess15 7.210
## ManufacturingProcess06 6.760
## ManufacturingProcess20 6.752
## ManufacturingProcess28 6.271
## BiologicalMaterial08 6.165
## BiologicalMaterial01 5.780
## ManufacturingProcess39 5.656
## ManufacturingProcess27 4.883top15 <- importance$importance %>%
as.data.frame() %>%
mutate(Variable = row.names(.)) %>%
slice_max(order_by = Overall, n = 15) %>%
mutate(Scaled = (Overall - min(Overall)) / (max(Overall) - min(Overall)) * 100) %>%
arrange(desc(Scaled))
ggplot(top15, aes(x = reorder(Variable, Scaled), y = Scaled)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(x = "Variable", y = "Scaled Importance", title = "Top 10 Variable Importance") +
theme_minimal()
(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?
Although ManufacturingProcess32 and ManufacturingProcess09 have the highest correlation with Yield, many other ManufacturingProcesses have negative correlations. On the other hand, almost all Biological Materials have a positive correlation with Yield.
Ultimately, variable importance in random forests reflects a combination of the influence of the variable itself, its interactions with other variables, and its non-linear relationship with the outcome. Unlike correlations, which are useful for initial exploration, random forests capture this whole picture, so some inconsistencies may arise when comparing the two.
selected_variables <- c("Yield","ManufacturingProcess13", "ManufacturingProcess32", "ManufacturingProcess09",
"BiologicalMaterial03", "ManufacturingProcess31", "ManufacturingProcess17",
"BiologicalMaterial06", "BiologicalMaterial11", "BiologicalMaterial12",
"BiologicalMaterial02", "BiologicalMaterial05", "ManufacturingProcess36",
"ManufacturingProcess15", "ManufacturingProcess06", "ManufacturingProcess20",
"ManufacturingProcess28", "BiologicalMaterial08", "BiologicalMaterial01",
"ManufacturingProcess39", "ManufacturingProcess27")
selected_data <- ChemicalManufacturingProcess %>%
select(all_of(selected_variables)) %>%
na.omit()
cor_matrix <- cor(selected_data, use = "complete.obs")
yield_cor <- cor_matrix["Yield", ]
yield_cor_sorted <- sort(yield_cor, decreasing = TRUE)
cor_matrix_sorted <- cor_matrix[names(yield_cor_sorted), names(yield_cor_sorted)]
corrplot(cor_matrix_sorted, method = "color", type = "upper",
order = "original", diag = FALSE,
tl.col = "darkblue", tl.srt = 90)