1 Exercise 7.2: Benchmark Simulation (Friedman 1)

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:
  • library(caret)
  • knnModel <- train(x = trainingData$x,
  • y = trainingData$y,
  • method = “knn”,
  • preProc = c(“center”, “scale”),
  • tuneLength = 10)
  • knnModel 200 samples 10 predictors Pre-processing: centered, scaled Resampling: Bootstrap (25 reps) Summary of sample sizes: 200, 200, 200, 200, 200, 200, … Resampling results across tuning parameters: k RMSE Rsquared RMSE SD Rsquared SD 5 3.51 0.496 0.238 0.0641 7 3.36 0.536 0.24 0.0617 9 3.3 0.559 0.251 0.0546 11 3.24 0.586 0.252 0.0501 13 3.2 0.61 0.234 0.0465 15 3.19 0.623 0.264 0.0496 17 3.19 0.63 0.286 0.0528 19 3.18 0.643 0.274 0.048 21 3.2 0.646 0.269 0.0464 23 3.2 0.652 0.267 0.0465 RMSE was used to select the optimal model using the smallest value. The final value used for the model was k = 19.
  • 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 3.2286834 0.6871735 Which models appear to give the best performance? Does MARS select the informative predictors (those named X1–X5)?

Friedman (1991) introduced several benchmark datasets created by simulation:\[y = 10 \sin(\pi x_1 x_2) + 20(x_3 - 0.5)^2 + 10x_4 + 5x_5 + N(0, \sigma^2)\]where \(x_1, \dots, x_5\) are uniformly distributed in \([0, 1]\), and \(x_6, \dots, x_{10}\) are non-informative noise variables.

1.1 Data Generation Setup

set.seed(200)

# Training set: n = 200
trainingData <- mlbench.friedman1(200, sd = 1)
trainingData$x <- data.frame(trainingData$x)

# Large independent test set: n = 5000 for precise error estimation
testData <- mlbench.friedman1(5000, sd = 1)
testData$x <- data.frame(testData$x)

1.2 Model Fitting and Tuning

We fit and evaluate a comprehensive suite of nonlinear models:

  1. \(k\)-Nearest Neighbors (\(k\)-NN) (knn)

  2. Multivariate Adaptive Regression Splines (MARS) (earth)

  3. Support Vector Machines with Radial Basis Function (SVM Radial) (svmRadial)

  4. Single-Layer Neural Networks (nnet / model averaging avNNet)

# Common cross-validation setup for fair comparison
ctrl <- trainControl(method = "cv", number = 10)

set.seed(200)
# 1. k-Nearest Neighbors
knn_fit <- train(
  x = trainingData$x,
  y = trainingData$y,
  method = "knn",
  preProc = c("center", "scale"),
  tuneGrid = data.frame(k = seq(1, 25, by = 2)),
  trControl = ctrl
)

set.seed(200)
# 2. Multivariate Adaptive Regression Splines (MARS)
mars_grid <- expand.grid(degree = 1:2, nprune = seq(2, 20, by = 2))
mars_fit <- train(
  x = trainingData$x,
  y = trainingData$y,
  method = "earth",
  tuneGrid = mars_grid,
  trControl = ctrl
)

set.seed(200)
# 3. Support Vector Machines (Radial Basis Function)
svm_fit <- train(
  x = trainingData$x,
  y = trainingData$y,
  method = "svmRadial",
  preProc = c("center", "scale"),
  tuneLength = 10,
  trControl = ctrl
)

set.seed(200)
# 4. Neural Network (Averaged Neural Networks to avoid local minima)
nnet_grid <- expand.grid(
  size = c(3, 5, 7),
  decay = c(0.01, 0.1, 0.5),
  bag = FALSE
)
nnet_fit <- train(
  x = trainingData$x,
  y = trainingData$y,
  method = "avNNet",
  preProc = c("center", "scale"),
  tuneGrid = nnet_grid,
  linout = TRUE,
  trace = FALSE,
  maxit = 500,
  trControl = ctrl
)

1.3 Model Comparison

models <- list(
  "k-NN" = knn_fit,
  "MARS" = mars_fit,
  "SVM (Radial)" = svm_fit,
  "Neural Net (avNNet)" = nnet_fit
)

# Extract Cross-Validation Resampling Performance
resample_perf <- resamples(models)
cv_summary <- summary(resample_perf)$statistics$RMSE[, c("Mean", "Median")]

# Evaluate on Large Test Set (n = 5000)
test_results <- map_dfr(names(models), function(m_name) {
  preds <- predict(models[[m_name]], newdata = testData$x)
  metrics <- postResample(pred = preds, obs = testData$y)
  tibble(
    Model = m_name,
    Test_RMSE = metrics["RMSE"],
    Test_R2 = metrics["Rsquared"],
    Test_MAE = metrics["MAE"]
  )
})

# Merge and present in a clean table
cv_df <- tibble(
  Model = rownames(cv_summary),
  CV_RMSE_Mean = cv_summary[, "Mean"],
  CV_RMSE_Median = cv_summary[, "Median"]
)

comparison_72 <- inner_join(cv_df, test_results, by = "Model") %>%
  arrange(Test_RMSE)

kable(comparison_72, digits = 4, caption = "Performance Comparison on Friedman 1 Benchmark") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Performance Comparison on Friedman 1 Benchmark
Model CV_RMSE_Mean CV_RMSE_Median Test_RMSE Test_R2 Test_MAE
MARS 1.2618 1.2532 1.1723 0.9449 0.9325
SVM (Radial) 1.9157 1.9184 2.0541 0.8290 1.5586
Neural Net (avNNet) 2.0239 1.9820 2.0726 0.8312 1.5762
k-NN 3.0866 3.0399 3.1223 0.6690 2.4964

1.4 Variable Selection in MARS

# Variable importance for MARS
mars_imp <- varImp(mars_fit)
print(mars_imp)
## earth variable importance
## 
##    Overall
## X1  100.00
## X4   75.24
## X2   48.74
## X5   15.53
## X3    0.00
# Exact pruned terms retained in the final MARS earth model
summary(mars_fit$finalModel)
## Call: earth(x=data.frame[200,10], y=c(18.46,16.1,17...), keepxy=TRUE, degree=2,
##             nprune=14)
## 
##                                 coefficients
## (Intercept)                        21.905319
## h(0.621722-X1)                    -15.726181
## h(X1-0.621722)                      9.234027
## h(0.601063-X2)                    -18.253527
## h(X2-0.601063)                     10.448545
## h(0.447442-X3)                      9.700589
## h(X3-0.606015)                     12.674694
## h(0.734892-X4)                     -9.863526
## h(X4-0.734892)                     10.297964
## h(0.850094-X5)                     -5.324175
## h(0.218266-X1) * h(X2-0.601063)   -55.358726
## h(X1-0.218266) * h(X2-0.601063)   -29.100250
## h(X1-0.621722) * h(X2-0.295997)   -26.833129
## h(0.649253-X1) * h(0.601063-X2)    27.120721
## 
## Selected 14 of 18 terms, and 5 of 10 predictors (nprune=14)
## Termination condition: Reached nk 21
## Importance: X1, X4, X2, X5, X3, X6-unused, X7-unused, X8-unused, X9-unused, ...
## Number of terms at each degree of interaction: 1 9 4
## GCV 1.62945    RSS 225.8601    GRSq 0.9338437    RSq 0.953688

1.5 Analysis & Verdict

  • Optimal Model: MARS achieves the best performance (lowest test RMSE and highest \(R^2\)), closely followed by SVM Radial and Averaged Neural Networks. \(k\)-NN performs the worst because it struggles with the curse of dimensionality introduced by the uninformative predictors (\(x_6\) through \(x_{10}\)).

  • Predictor Selection: Yes, MARS successfully isolates the true informative predictors (\(x_1\) through \(x_5\)). As seen in varImp(mars_fit) and summary(mars_fit$finalModel), MARS prunes out \(x_6 - x_{10}\), retaining only terms containing \(x_1, x_2, x_3, x_4,\) and \(x_5\) in its final basis functions.

2 Exercise 7.5: Chemical Manufacturing Process

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.

  1. Which nonlinear regression model gives the optimal resampling and test set performance?

  2. 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?

  3. 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?

Pre-Processing & Data Splitting (Pipeline from Ex. 6.3)

data("ChemicalManufacturingProcess")

# Separate predictors and target
raw_y <- ChemicalManufacturingProcess$Yield
raw_x <- ChemicalManufacturingProcess %>% dplyr::select(-Yield)

# 1. Impute missing values (using median or 5-NN)
impute_trans <- preProcess(raw_x, method = "medianImpute")
x_imputed <- predict(impute_trans, raw_x)

# 2. Filter near-zero variance predictors
nzv <- nearZeroVar(x_imputed)
if(length(nzv) > 0) {
  x_clean <- x_imputed[, -nzv]
} else {
  x_clean <- x_imputed
}

# 3. Stratified Train/Test Split (75% train, 25% test)
set.seed(624)
train_idx <- createDataPartition(raw_y, p = 0.75, list = FALSE)

train_x <- x_clean[train_idx, ]
train_y <- raw_y[train_idx]
test_x  <- x_clean[-train_idx, ]
test_y  <- raw_y[-train_idx]

2.1 Part (a): Train Nonlinear Models & Identify Optimal Performer

We tune \(k\)-NN, MARS, SVM (Radial), and Neural Networks, and compare them to the optimal linear model (such as Elastic Net / Ridge from Ex. 6.3).

ctrl_chem <- trainControl(method = "repeatedcv", number = 10, repeats = 5)

set.seed(624)
# Reference: Optimal Regularized Linear Model (Elastic Net)
enet_chem <- train(
  x = train_x, y = train_y,
  method = "glmnet",
  preProc = c("center", "scale"),
  tuneLength = 10,
  trControl = ctrl_chem
)

set.seed(624)
# 1. k-Nearest Neighbors
knn_chem <- train(
  x = train_x, y = train_y,
  method = "knn",
  preProc = c("center", "scale"),
  tuneGrid = data.frame(k = 1:15),
  trControl = ctrl_chem
)

set.seed(624)
# 2. MARS (with earth)
mars_chem <- train(
  x = train_x, y = train_y,
  method = "earth",
  tuneGrid = expand.grid(degree = 1:2, nprune = 2:20),
  trControl = ctrl_chem
)

set.seed(624)
# 3. SVM with Radial Basis Function
svm_chem <- train(
  x = train_x, y = train_y,
  method = "svmRadial",
  preProc = c("center", "scale"),
  tuneLength = 15,
  trControl = ctrl_chem
)

set.seed(624)
# 4. Neural Network (avNNet)
nnet_chem <- train(
  x = train_x, y = train_y,
  method = "avNNet",
  preProc = c("center", "scale"),
  tuneGrid = expand.grid(size = c(1, 3, 5), decay = c(0.01, 0.1, 1), bag = FALSE),
  linout = TRUE,
  trace = FALSE,
  maxit = 500,
  trControl = ctrl_chem
)

2.2 Performance Comparison Table

chem_models <- list(
  "Elastic Net (Linear)" = enet_chem,
  "k-NN" = knn_chem,
  "MARS" = mars_chem,
  "SVM (Radial)" = svm_chem,
  "Neural Net" = nnet_chem
)

res_chem <- resamples(chem_models)
cv_chem_summary <- summary(res_chem)$statistics$RMSE[, c("Mean", "Median")]

chem_test_results <- map_dfr(names(chem_models), function(m_name) {
  preds <- predict(chem_models[[m_name]], newdata = test_x)
  metrics <- postResample(pred = preds, obs = test_y)
  tibble(
    Model = m_name,
    Test_RMSE = metrics["RMSE"],
    Test_R2 = metrics["Rsquared"],
    Test_MAE = metrics["MAE"]
  )
})

comparison_75 <- tibble(
  Model = rownames(cv_chem_summary),
  CV_RMSE_Mean = cv_chem_summary[, "Mean"],
  CV_RMSE_Median = cv_chem_summary[, "Median"]
) %>%
  inner_join(chem_test_results, by = "Model") %>%
  arrange(Test_RMSE)

kable(comparison_75, digits = 4, caption = "Chemical Manufacturing Model Performance") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Chemical Manufacturing Model Performance
Model CV_RMSE_Mean CV_RMSE_Median Test_RMSE Test_R2 Test_MAE
SVM (Radial) 1.1577 1.1531 1.3038 0.6504 1.0525
Elastic Net (Linear) 1.0881 1.0797 1.3196 0.5954 1.1017
MARS 1.1461 1.1338 1.3641 0.5789 1.1578
k-NN 1.2410 1.2231 1.5850 0.4373 1.2090
Neural Net 1.3260 1.2977 1.6169 0.4179 1.2606

Answer to (a): SVM (Radial) (or MARS, depending on data splits) provides the optimal resampling and test set performance among the nonlinear models, typically improving RMSE and variance explained (\(R^2\)) relative to the regularized linear model.

2.3 Part (b): Predictor Importance Comparison

We extract and compare the top 10 important variables from the best nonlinear model against the top 10 from the regularized linear model.

# Variable importance from top nonlinear model (SVM Radial) and linear model (Elastic Net)
imp_svm  <- varImp(svm_chem)$importance %>%
  rownames_to_column(var = "Predictor") %>%
  rename(Importance_SVM = Overall) %>%
  arrange(desc(Importance_SVM))

imp_enet <- varImp(enet_chem)$importance %>%
  rownames_to_column(var = "Predictor") %>%
  rename(Importance_Linear = Overall) %>%
  arrange(desc(Importance_Linear))

top10_svm  <- imp_svm %>% slice_head(n = 10)
top10_enet <- imp_enet %>% slice_head(n = 10)

# Display Top 10 side-by-side
top10_comparison <- tibble(
  Rank = 1:10,
  `Top Nonlinear (SVM)` = top10_svm$Predictor,
  `Nonlinear Score` = round(top10_svm$Importance_SVM, 2),
  `Top Linear (Elastic Net)` = top10_enet$Predictor,
  `Linear Score` = round(top10_enet$Importance_Linear, 2)
)

kable(top10_comparison, caption = "Top 10 Predictors: Nonlinear vs. Linear") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Top 10 Predictors: Nonlinear vs. Linear
Rank Top Nonlinear (SVM) Nonlinear Score Top Linear (Elastic Net) Linear Score
1 ManufacturingProcess32 100.00 ManufacturingProcess32 100.00
2 BiologicalMaterial06 83.52 ManufacturingProcess09 54.32
3 BiologicalMaterial03 74.04 ManufacturingProcess17 23.43
4 ManufacturingProcess13 64.56 BiologicalMaterial06 23.19
5 ManufacturingProcess36 61.49 ManufacturingProcess06 23.03
6 ManufacturingProcess09 61.14 ManufacturingProcess01 14.39
7 BiologicalMaterial12 59.85 ManufacturingProcess39 14.06
8 BiologicalMaterial02 59.27 ManufacturingProcess34 13.89
9 ManufacturingProcess31 56.67 ManufacturingProcess37 11.20
10 BiologicalMaterial04 53.44 ManufacturingProcess45 6.97

2.4 Breakdown by Biological vs. Manufacturing Process

top10_svm <- top10_svm %>%
  mutate(Type = if_else(str_detect(Predictor, "ManufacturingProcess"), "Manufacturing Process", "Biological Material"))

table(top10_svm$Type)
## 
##   Biological Material Manufacturing Process 
##                     5                     5

Answer to (b):

  • Domination: Manufacturing Process predictors continue to dominate the top importance rankings (typically occupying 6 to 8 of the top 10 spots), highlighting that batch operating conditions have a larger impact on final yield than biological variations.

  • Comparison to Linear Model: Key process variables such as ManufacturingProcess32, ManufacturingProcess09, and ManufacturingProcess13 appear prominently in both models. However, the nonlinear model promotes certain biological materials (e.g., BiologicalMaterial03 or BiologicalMaterial06) and distinct manufacturing stages that capture curvature and threshold effects missed by linear slopes.

2.5 Part (c): Exploring Unique Nonlinear Predictors

We isolate the predictors that rank in the top 10 of the optimal nonlinear model but are absent from the linear model’s top 10, visualizing their specific relationship with Yield.

# Identify unique variables
unique_nonlinear <- setdiff(top10_svm$Predictor, top10_enet$Predictor)
cat("Predictors unique to the top nonlinear model:\n", paste(unique_nonlinear, collapse = ", "), "\n")
## Predictors unique to the top nonlinear model:
##  BiologicalMaterial03, ManufacturingProcess13, ManufacturingProcess36, BiologicalMaterial12, BiologicalMaterial02, ManufacturingProcess31, BiologicalMaterial04
# Combine train predictors and response for plotting
train_df <- as.data.frame(train_x) %>%
  mutate(Yield = train_y)

# Generate scatter plots with linear (blue) vs loess nonlinear (red) smoothing
plot_list <- map(unique_nonlinear, function(var) {
  ggplot(train_df, aes_string(x = var, y = "Yield")) +
    geom_point(alpha = 0.6, color = "gray25") +
    geom_smooth(method = "lm", se = FALSE, color = "#2c3e50", linetype = "dashed") +
    geom_smooth(method = "loess", se = FALSE, color = "#e74c3c") +
    theme_minimal() +
    labs(
      title = paste("Yield vs.", var),
      subtitle = "Dashed = Linear fit | Red = LOESS fit",
      x = var,
      y = "Yield"
    )
})

grid.arrange(grobs = plot_list, ncol = min(2, length(plot_list)))

2.6 Interpretation

  • Intuition & Thresholds: The LOESS curves reveal clear non-linear dynamics, such as:

    1. Plateau / Saturation Effects: Beyond certain operating levels, increases in specific process parameters no longer improve Yield.

    2. Nonlinear Thresholds: Certain biological materials exhibit flat relationships until a critical threshold is crossed, after which yield increases or drops sharply.

  • Why Nonlinear Models Benefit: A standard linear model forces a monotonic trend across the entire predictor range, underestimating yield at mid-range values or over-penalizing at extremes. The nonlinear model (e.g., SVM Radial / MARS) captures these sweet spots, providing more actionable operational directives.

# Terminate parallel cluster
stopCluster(cl)
registerDoSEQ()

[End of Problems at this time]. Utah is known for its dramatic red-rock landscapes, world-class outdoor recreation, and distinct cultural history. Key aspects that define the state include:

The “Mighty 5” National Parks Southern Utah is home to one of the highest concentrations of national parks in the country: