library(mlbench)
library(caret)
library(earth)
library(randomForest)
library(gbm)
library(nnet) 
library(e1071)
library(Metrics)
library(ggplot2)

Exercise 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(π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:

Tune several models on these data.

RMSE was used to select the optimal model using the smallest value.

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

set.seed(200)
trainingData <- mlbench.friedman1(200, sd = 1)
trainingData$x <- data.frame(trainingData$x)
testData <- mlbench.friedman1(5000, sd = 1)
testData$x <- data.frame(testData$x)

# cross validation
ctrl <- trainControl(method = "boot", number = 25)

KNN

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

MARS

set.seed(1)
marsModel <- train(x = trainingData$x, y = trainingData$y,
                   method = "earth", tuneLength = 10,
                   trControl = ctrl)

Random forest

set.seed(1)
rfModel <- train(x = trainingData$x, y = trainingData$y,
                 method = "rf", tuneLength = 5, trControl = ctrl)

Gradient boosted machine (GBM)

set.seed(1)
gbmModel <- train(x = trainingData$x, y = trainingData$y,
                  method = "gbm", verbose = FALSE, trControl = ctrl)

Support vector machine (SVM)

set.seed(1)
svmModel <- train(x = trainingData$x, y = trainingData$y,
                  method = "svmRadial", preProc = c("center", "scale"),
                  tuneLength = 10, trControl = ctrl)

Evaluate on the test set

models <- list(knn = knnModel, mars = marsModel, rf = rfModel, gbm = gbmModel, svm = svmModel)
results <- lapply(models, function(model) {
  preds <- predict(model, newdata = testData$x)
  postResample(pred = preds, obs = testData$y)})
results_df <- do.call(rbind, results)
print(round(results_df, 4))

##        RMSE Rsquared    MAE
## knn  3.2041   0.6820 2.5683
## mars 1.8136   0.8677 1.3912
## rf   2.4155   0.7892 1.9092
## gbm  1.8162   0.8676 1.3959
## svm  2.0636   0.8273 1.5678

Assess MARS variable selection

summary(marsModel$finalModel)

## Call: earth(x=data.frame[200,10], y=c(18.46,16.1,17...), keepxy=TRUE, degree=1,
##             nprune=12)
## 
##                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=12)
## 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.9183982

Feature importance

varImp(marsModel)

## earth variable importance
## 
##    Overall
## X1  100.00
## X4   82.08
## X2   62.79
## X5   38.07
## X3   25.80
## X6    0.00

Notes:

Yes, MARS did identify the correct informative predictors, based on both the model summary and variable importance plot.
Selected predictors: X1, X2, X3, X4, X5, and X6.
X1–X5 are informative according to the data-generating process, but X6 us uninformative.

===============================================================================

Exercise 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.

Which nonlinear regression model gives the optimal resampling and test set performance?
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?
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?

library(AppliedPredictiveModeling)
# fixing errors with loading the data
load(system.file("data", "chemicalManufacturingProcess.RData", package = "AppliedPredictiveModeling"))
# str(ChemicalManufacturingProcess)
# ls()


processPredictors <- ChemicalManufacturingProcess[, -1]
yield <- ChemicalManufacturingProcess$Yield


dim(processPredictors)

## [1] 176  57

length(yield)

## [1] 176

sum(is.na(processPredictors))

## [1] 106

preProc <- preProcess(processPredictors, method = c("knnImpute", "center", "scale"))

# Impute and scale
processed_data <- predict(preProc, processPredictors)

# Check
sum(is.na(processed_data))

## [1] 0

Add Yield Back and Prepare Data Split

# yield to the processed predictors
processed_data$Yield <- yield

# Split (75/25)
set.seed(123)
trainIndex <- createDataPartition(processed_data$Yield, p = 0.75, list = FALSE)
trainData <- processed_data[trainIndex, ]
testData  <- processed_data[-trainIndex, ]

# X and y
trainX <- trainData[, -ncol(trainData)]
trainY <- trainData$Yield
testX  <- testData[, -ncol(testData)]
testY  <- testData$Yield

Train Nonlinear Models

ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)

# Random Forest
set.seed(100)
rfModel <- train(trainX, trainY, method = "rf", trControl = ctrl, tuneLength = 5)

# GBM
set.seed(100)
gbmModel <- train(trainX, trainY, method = "gbm", verbose = FALSE, trControl = ctrl, tuneLength = 5)

# MARS
set.seed(100)
marsModel <- train(trainX, trainY, method = "earth", trControl = ctrl, tuneLength = 10)

# SVM (Radial Basis)
set.seed(100)
svmModel <- train(trainX, trainY, method = "svmRadial", trControl = ctrl, tuneLength = 10)

# KNN
set.seed(100)
knnModel <- train(trainX, trainY, method = "knn", trControl = ctrl, tuneLength = 10)

Test Set Performance Comparison

models <- list(rf = rfModel, gbm = gbmModel, mars = marsModel, svm = svmModel, knn = knnModel)

results_test <- lapply(models, function(mod) {
  pred <- predict(mod, newdata = testX)
  df <- data.frame(
    RMSE = rmse(testY, pred),
    R2   = R2(pred, testY),
    MAE  = mae(testY, pred))
  names(df) <- c("RMSE", "R2", "MAE")
  return(df)})

names(results_test) <- names(models)

results_test_df <- do.call(rbind, results_test)
print(round(results_test_df, 4))

##        RMSE     R2    MAE
## rf   1.1869 0.5429 0.8567
## gbm  1.2676 0.4961 0.9857
## mars 1.2116 0.5339 0.9869
## svm  1.1194 0.5899 0.8742
## knn  1.3761 0.3837 1.0971

Notes:

The Support Vector Machine (SVM) is the best nonlinear regression model for this dataset based on all three metrics (RMSE, R2, MAE).

Variable Importance

# SVM
varImp_svm <- varImp(svmModel)
print(varImp_svm)

## loess r-squared variable importance
## 
##   only 20 most important variables shown (out of 57)
## 
##                        Overall
## ManufacturingProcess32  100.00
## BiologicalMaterial06     84.03
## ManufacturingProcess17   81.25
## ManufacturingProcess13   80.20
## ManufacturingProcess36   77.39
## ManufacturingProcess31   77.37
## BiologicalMaterial03     76.02
## ManufacturingProcess09   74.36
## BiologicalMaterial02     67.99
## ManufacturingProcess06   62.05
## BiologicalMaterial12     56.66
## ManufacturingProcess33   54.14
## ManufacturingProcess30   53.56
## ManufacturingProcess11   53.02
## BiologicalMaterial04     49.94
## ManufacturingProcess29   48.85
## ManufacturingProcess02   48.08
## BiologicalMaterial11     45.12
## BiologicalMaterial08     41.55
## BiologicalMaterial01     41.27

# MARS
varImp_mars <- varImp(marsModel)
print(varImp_mars)

## earth variable importance
## 
##                        Overall
## ManufacturingProcess32     100
## ManufacturingProcess09       0

# top 10 important variables
plot(varImp_svm, top = 10, main = "Top 10 Important Variables - SVM")

plot(varImp_mars, top = 10, main = "Top 10 Important Variables - MARS")

Notes:

Out of the top 10 SVM predictors: 7 are manufacturing process variables and 3 are biological material variables.
The best nonlinear model (SVM) heavily favors process variables, though some biological variables (especially BiologicalMaterial06) are also highly influential.
MARS model selected only: ManufacturingProcess32 (important) ManufacturingProcess09 (assigned 0 importance). So MARS failed to identify most informative predictors, which explains its lower predictive performance.

Top Variable Visualizations

topVars <- rownames(
  head(varImp_svm$importance[order(varImp_svm$importance$Overall, decreasing = TRUE), , drop = FALSE], 5))

for (var in topVars) {
  ggplot(data.frame(x = trainX[[var]], y = trainY), aes(x = x, y = y)) +
    geom_point(alpha = 0.5) +
    geom_smooth(method = "loess", se = FALSE, color = "blue") +
    labs(title = paste("Yield vs", var), x = var, y = "Yield") -> p
  print(p)}

Notes:

ManufacturingProcess32: Shows a clear nonlinear S-shaped pattern. Yield is flat (or slightly declining) at low values, increases sharply around 0, and then levels off again. This aligns with its top importance in both SVM and MARS. Suggests a threshold effect—a sweet spot in the mid-range optimizes yield.
BiologicalMaterial06: Displays a mild upward curve—yield increases with increasing values up to ~1.5–2, then may plateau or decline. Indicates that raw material quality (BiologicalMaterial06) is important but non-monotonic in its effect. Linear models might underestimate this effect due to curvature.
ManufacturingProcess17: Strong inverted U-shape, showing peak yield around 0.5. Lower and higher values both reduce yield. This is nonlinear and non-monotonic, justifying the SVM’s flexibility.
ManufacturingProcess13: Similar to Process17, showing an inverse trend: as this variable increases, yield decreases sharply, then flattens. May suggest diminishing returns or harmful effects at higher levels.
ManufacturingProcess36: A step-like pattern with sharp local dips and clustering at discrete values. This might indicate that the process has categorical settings or discrete operational states. Such patterns are difficult to model linearly and benefit from models like SVM with RBF kernels.
These relationships justify the superiority of the SVM model, which captures: Nonlinearity, Thresholds, Plateaus, Non-monotonicity.
ManufacturingProcess variables dominate the top importance rankings and exhibit complex functional forms, which traditional linear models and even MARS may fail to capture fully.

Compare to Linear Model

lmModel <- train(trainX, trainY, method = "lm", trControl = ctrl)
varImp_lm <- varImp(lmModel)

top10_lm  <- rownames(head(varImp_lm$importance[order(varImp_lm$importance$Overall, decreasing = TRUE), , drop = FALSE], 10))
top10_svm <- rownames(head(varImp_svm$importance[order(varImp_svm$importance$Overall, decreasing = TRUE), , drop = FALSE], 10))

setdiff(top10_svm, top10_lm)

## [1] "BiologicalMaterial06"   "ManufacturingProcess17" "ManufacturingProcess13"
## [4] "ManufacturingProcess36" "ManufacturingProcess31" "BiologicalMaterial02"  
## [7] "ManufacturingProcess06"

setdiff(top10_lm, top10_svm)

## [1] "ManufacturingProcess45" "ManufacturingProcess29" "BiologicalMaterial09"  
## [4] "ManufacturingProcess37" "ManufacturingProcess28" "ManufacturingProcess04"
## [7] "BiologicalMaterial12"

Conclusion:

Among the top predictors identified by the SVM model, multiple manufacturing process variables exhibited nonlinear or non-monotonic associations with yield, such as threshold effects (e.g., ManufacturingProcess32), inverse U-shapes (e.g., Process17), and diminishing returns (e.g., Process13). This shows the importance of using nonlinear modeling approaches to uncover such relationships and optimize process settings. MARS identified ManufacturingProcess32 as influential, but it failed to detect several other key predictors likely due to its additive constraints. Biological variables also contributed, notably BiologicalMaterial06, which demonstrated a curved but positive association with yield.

Comparing the top 10 predictors from linear and SVM models revealed only partial overlap, with 7 variables unique to each approach. The SVM model identified several predictors that showed clear nonlinear or non-monotonic relationships with yield in exploratory plots. In contrast, variables like ManufacturingProcess45 and BiologicalMaterial09, ranked highly by the linear model, were deprioritized in the SVM model. This suggests that their linear associations may not reflect true predictive relevance.

624_M11 HW8

Fares A

2025-04-07

Exercise 7.2.

Exercise 7.5