DATA 624 Homework 8

Exercise 7.2
Exercise 7.5

knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(cli.unicode = FALSE)
library(AppliedPredictiveModeling)
library(mlbench)
library(caret)
library(earth)
library(kernlab)
library(tidyverse)

Exercise 7.2

Friedman (1991) introduced several benchmark data sets created by simulation. One of these used the nonlinear equation

\[y = 10\sin(\pi x_1 x_2) + 20(x_3 - 0.5)^2 + 10x_4 + 5x_5 + N(0, \sigma^2)\]

where the x values are uniformly distributed on [0, 1]. Five further non-informative variables (X6 to X10) are also created.

Simulating the data

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

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

dim(trainingData$x)

## [1] 200  10

dim(testData$x)

## [1] 5000   10

featurePlot(trainingData$x, trainingData$y)

The feature plot shows what the equation implies. X4 has a clear positive linear relationship with the response and X5 a weaker one, X1 and X2 show looser patterns because they enter through the sine term, and X3 shows the U shape produced by its squared term. X6 to X10 show nothing at all, which is correct since they are noise.

Tuning several models

I tune the four nonlinear methods from this chapter, each with 10-fold cross validation.

ctrl <- trainControl(method = "cv", number = 10)

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

set.seed(624)
mars_model <- train(trainingData$x, trainingData$y, method = "earth",
                    tuneGrid = expand.grid(degree = 1:2, nprune = 2:20), trControl = ctrl)

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

set.seed(624)
nnet_model <- train(trainingData$x, trainingData$y, method = "nnet",
                    preProc = c("center", "scale"),
                    tuneGrid = expand.grid(size = 1:8, decay = c(0, 0.01, 0.1)),
                    trControl = ctrl, linout = TRUE, trace = FALSE, maxit = 500,
                    MaxNWts = 10 * (ncol(trainingData$x) + 1) + 10 + 1)

models <- list(KNN = knn_model, MARS = mars_model, SVM = svm_model, NNet = nnet_model)

map_dfr(names(models), function(nm) {
  m <- models[[nm]]
  test_perf <- postResample(predict(m, testData$x), testData$y)
  tibble(model = nm,
         tuning = paste(names(m$bestTune), round(unlist(m$bestTune), 4), collapse = ", "),
         cv_RMSE = min(m$results$RMSE),
         cv_Rsquared = max(m$results$Rsquared, na.rm = TRUE),
         test_RMSE = test_perf[["RMSE"]],
         test_Rsquared = test_perf[["Rsquared"]])
}) |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  arrange(test_RMSE)

## # A tibble: 4 x 6
##   model tuning              cv_RMSE cv_Rsquared test_RMSE test_Rsquared
##   <chr> <chr>                 <dbl>       <dbl>     <dbl>         <dbl>
## 1 MARS  nprune 14, degree 2    1.28       0.939      1.18         0.945
## 2 SVM   sigma 0.0561, C 4      1.85       0.872      2.08         0.827
## 3 NNet  size 4, decay 0.01     2.21       0.813      2.14         0.822
## 4 KNN   k 9                    3.03       0.738      3.13         0.657

plot(mars_model, main = "MARS tuning")

plot(svm_model, main = "SVM tuning")

Which models appear to give the best performance?

MARS is clearly the best, and the margin is large rather than marginal:

MARS: test RMSE 1.175, test R2 0.945 (degree 2, 14 terms)
SVM: test RMSE 2.079, test R2 0.827
Neural network: test RMSE 2.136, test R2 0.822
KNN: test RMSE 3.131, test R2 0.657

The cross-validated and test set rankings agree exactly, which gives confidence in the ordering. The 5,000-sample test set also means the test estimates are precise, unlike the small test sets in the previous homework.

MARS wins because its structure matches how the data was generated. The true equation is additive in separate pieces, with one interaction between X1 and X2, and MARS builds exactly that: piecewise linear hinge functions plus selected products of them. Its test RMSE of 1.175 is close to the noise floor of the simulation, since the errors were generated with a standard deviation of 1, so there is almost no systematic error left.

The SVM and the neural network do reasonably well but have to approximate the surface smoothly across all ten predictors, including the five that carry no signal. KNN is the weakest, which is expected: it measures distance using all ten predictors, so the five noise variables dilute every neighbourhood it forms. This is the curse of dimensionality in miniature.

Does MARS select the informative predictors (those named X1 to X5)?

varImp(mars_model)

## earth variable importance
## 
##    Overall
## X1  100.00
## X4   75.24
## X2   48.74
## X5   15.53
## X3    0.00

summary(mars_model$finalModel)$coefficients

##                                        y
## (Intercept)                    21.905319
## h(X4-0.734892)                 10.297964
## h(0.734892-X4)                 -9.863526
## h(X1-0.621722)                  9.234027
## h(0.621722-X1)                -15.726181
## h(X2-0.601063)                 10.448545
## h(0.601063-X2)                -18.253527
## h(0.850094-X5)                 -5.324175
## h(0.447442-X3)                  9.700589
## h(0.649253-X1)*h(0.601063-X2)  27.120721
## h(X1-0.218266)*h(X2-0.601063) -29.100250
## h(0.218266-X1)*h(X2-0.601063) -55.358726
## h(X1-0.621722)*h(X2-0.295997) -26.833129
## h(X3-0.606015)                 12.674694

Yes, perfectly. The variable importance list contains only X1 to X5, and the hinge terms in the final model involve only those five predictors. X6 through X10 never appear, so MARS discarded every non-informative variable through its pruning step.

The ordering is also sensible. X1, X4 and X2 lead, followed by X5, with X3 last. X4 and X5 enter the true equation as plain linear terms with coefficients of 10 and 5, so X4 should outrank X5, and it does. X1 and X2 enter through the sine interaction, and the model captures this with product terms such as h(0.218266-X1)*h(X2-0.601063), which is MARS discovering the interaction without being told it exists. X3 enters only as a squared term, which the model approximates with two hinges, h(0.447442-X3) and h(X3-0.606015), forming a V shape that stands in for the parabola.

This feature selection is the main practical advantage of MARS here, and it explains the performance gap: MARS is the only one of the four models that effectively ignores the noise predictors.

Exercise 7.5

Exercise 6.3 describes data for a chemical manufacturing process. Using the same imputation, data splitting and pre-processing steps as before, I train several nonlinear regression models.

Same preparation as Exercise 6.3

data(ChemicalManufacturingProcess)
chem <- ChemicalManufacturingProcess

predictors <- chem[, -1]
yield <- chem$Yield

predictors <- predictors[, -nearZeroVar(predictors)]

impute_model <- preProcess(predictors, method = "knnImpute")
predictors_imputed <- predict(impute_model, predictors)

set.seed(624)
chem_rows <- createDataPartition(yield, p = 0.8, list = FALSE)

X_train <- predictors_imputed[chem_rows, ]
X_test  <- predictors_imputed[-chem_rows, ]
y_train <- yield[chem_rows]
y_test  <- yield[-chem_rows]

c(train = nrow(X_train), test = nrow(X_test))

## train  test 
##   144    32

This repeats the earlier steps exactly: drop BiologicalMaterial07 as near-zero variance, impute the 106 missing cells with K-nearest neighbours, and use the same 80/20 split with the same seed, so the results are directly comparable with the PLS model from Exercise 6.3.

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

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

set.seed(624)
chem_knn <- train(X_train, y_train, method = "knn",
                  preProc = c("center", "scale"), tuneLength = 15, trControl = chem_ctrl)

set.seed(624)
chem_mars <- train(X_train, y_train, method = "earth",
                   tuneGrid = expand.grid(degree = 1:2, nprune = 2:15), trControl = chem_ctrl)

set.seed(624)
chem_svm <- train(X_train, y_train, method = "svmRadial",
                  preProc = c("center", "scale"), tuneLength = 14, trControl = chem_ctrl)

set.seed(624)
chem_nnet <- train(X_train, y_train, method = "nnet",
                   preProc = c("center", "scale"),
                   tuneGrid = expand.grid(size = 1:6, decay = c(0, 0.01, 0.1, 0.5)),
                   trControl = chem_ctrl, linout = TRUE, trace = FALSE, maxit = 500,
                   MaxNWts = 6 * (ncol(X_train) + 1) + 6 + 1)

set.seed(624)
chem_pls <- train(X_train, y_train, method = "pls",
                  tuneLength = 20, trControl = chem_ctrl, preProc = c("center", "scale"))

The PLS model from Exercise 6.3 is included as the linear benchmark.

chem_models <- list(KNN = chem_knn, MARS = chem_mars, SVM = chem_svm,
                    NNet = chem_nnet, `PLS (linear)` = chem_pls)

chem_results <- map_dfr(names(chem_models), function(nm) {
  m <- chem_models[[nm]]
  test_perf <- postResample(predict(m, X_test), y_test)
  tibble(model = nm,
         tuning = paste(names(m$bestTune), round(unlist(m$bestTune), 4), collapse = ", "),
         cv_RMSE = min(m$results$RMSE),
         cv_Rsquared = max(m$results$Rsquared, na.rm = TRUE),
         test_RMSE = test_perf[["RMSE"]],
         test_Rsquared = test_perf[["Rsquared"]])
})

chem_results |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  arrange(cv_RMSE)

## # A tibble: 5 x 6
##   model        tuning             cv_RMSE cv_Rsquared test_RMSE test_Rsquared
##   <chr>        <chr>                <dbl>       <dbl>     <dbl>         <dbl>
## 1 SVM          sigma 0.014, C 8      1.15       0.606      1.14         0.747
## 2 MARS         nprune 4, degree 2    1.21       0.577      1.20         0.713
## 3 KNN          k 5                   1.32       0.487      1.41         0.618
## 4 PLS (linear) ncomp 3               1.36       0.523      1.20         0.678
## 5 NNet         size 1, decay 0.01    1.59       0.402      1.59         0.494

resamples(list(KNN = chem_knn, MARS = chem_mars, SVM = chem_svm,
               NNet = chem_nnet, PLS = chem_pls)) |>
  bwplot(metric = "RMSE", main = "Resampled RMSE across models")

The radial basis function SVM is the best on both measures, which is a reassuring agreement:

SVM (sigma 0.014, C = 8): resampled RMSE 1.154, test RMSE 1.145 with R2 of 0.747
MARS (degree 2, 4 terms): resampled RMSE 1.207, test RMSE 1.201
KNN (k = 5): resampled RMSE 1.325, test RMSE 1.411
Neural network (size 1, decay 0.01): resampled RMSE 1.587, test RMSE 1.586

The SVM also beats the linear PLS benchmark, which had a resampled RMSE of 1.363 and a test RMSE of 1.198. The improvement is meaningful but modest: roughly 15% on the resampled RMSE. The box plot shows the SVM and MARS distributions overlapping, so the gap between those two is not large, while KNN and especially the neural network are clearly behind. The neural network is handicapped here by the sample size, since 144 training runs is very little for a model of that flexibility, and the tuning settles on a single hidden unit.

That the nonlinear models improve on PLS at all suggests there is some curvature in the relationship between the process measurements and yield, which the next part explores.

b. Which predictors are most important in the optimal nonlinear regression model?

svm_imp <- varImp(chem_svm)
plot(svm_imp, top = 10, main = "Top 10 predictors, SVM model")

svm_top10 <- svm_imp$importance |>
  rownames_to_column("predictor") |>
  arrange(desc(Overall)) |>
  slice(1:10) |>
  mutate(type = if_else(str_detect(predictor, "Biological"), "Biological", "Process"),
         Overall = round(Overall, 1))

svm_top10

##                 predictor Overall       type
## 1  ManufacturingProcess32   100.0    Process
## 2    BiologicalMaterial06    89.3 Biological
## 3    BiologicalMaterial03    77.5 Biological
## 4  ManufacturingProcess36    74.6    Process
## 5  ManufacturingProcess09    73.9    Process
## 6  ManufacturingProcess13    73.2    Process
## 7  ManufacturingProcess31    67.9    Process
## 8    BiologicalMaterial02    66.9 Biological
## 9    BiologicalMaterial12    64.9 Biological
## 10 ManufacturingProcess06    59.1    Process

svm_top10 |> count(type)

##         type n
## 1 Biological 4
## 2    Process 6

An important caveat about this list: caret cannot compute a variable importance measure from the internals of a radial basis SVM, so it falls back on a model-free filter, scoring each predictor by how well a loess curve of the response against that predictor fits. The ranking therefore reflects each predictor’s univariate relationship with yield rather than the SVM’s own use of it. It is still informative, but it is not strictly the SVM’s internal view.

Do biological or process predictors dominate? Process predictors still lead, but less decisively than before: 6 of the top 10 are process measurements and 4 are biological. ManufacturingProcess32 is again first, but BiologicalMaterial06 and BiologicalMaterial03 now rank second and third.

Comparing with the top ten from the optimal linear model

pls_top10 <- varImp(chem_pls)$importance |>
  rownames_to_column("predictor") |>
  arrange(desc(Overall)) |>
  slice(1:10) |>
  mutate(type = if_else(str_detect(predictor, "Biological"), "Biological", "Process"),
         Overall = round(Overall, 1))

pls_top10

##                 predictor Overall       type
## 1  ManufacturingProcess32   100.0    Process
## 2  ManufacturingProcess09    82.8    Process
## 3  ManufacturingProcess36    82.3    Process
## 4  ManufacturingProcess13    76.0    Process
## 5  ManufacturingProcess17    70.3    Process
## 6  ManufacturingProcess06    62.5    Process
## 7    BiologicalMaterial06    57.7 Biological
## 8    BiologicalMaterial03    57.3 Biological
## 9    BiologicalMaterial02    57.1 Biological
## 10 ManufacturingProcess34    56.3    Process

pls_top10 |> count(type)

##         type n
## 1 Biological 3
## 2    Process 7

tibble(predictor = union(svm_top10$predictor, pls_top10$predictor)) |>
  mutate(in_SVM = predictor %in% svm_top10$predictor,
         in_PLS = predictor %in% pls_top10$predictor,
         status = case_when(in_SVM & in_PLS ~ "both",
                            in_SVM ~ "SVM only",
                            TRUE ~ "PLS only")) |>
  arrange(status, predictor)

## # A tibble: 12 x 4
##    predictor              in_SVM in_PLS status  
##    <chr>                  <lgl>  <lgl>  <chr>   
##  1 ManufacturingProcess17 FALSE  TRUE   PLS only
##  2 ManufacturingProcess34 FALSE  TRUE   PLS only
##  3 BiologicalMaterial12   TRUE   FALSE  SVM only
##  4 ManufacturingProcess31 TRUE   FALSE  SVM only
##  5 BiologicalMaterial02   TRUE   TRUE   both    
##  6 BiologicalMaterial03   TRUE   TRUE   both    
##  7 BiologicalMaterial06   TRUE   TRUE   both    
##  8 ManufacturingProcess06 TRUE   TRUE   both    
##  9 ManufacturingProcess09 TRUE   TRUE   both    
## 10 ManufacturingProcess13 TRUE   TRUE   both    
## 11 ManufacturingProcess32 TRUE   TRUE   both    
## 12 ManufacturingProcess36 TRUE   TRUE   both

The two lists agree on 8 of 10 predictors, including the top one, ManufacturingProcess32, which both models rank first. The differences are:

Unique to the SVM: ManufacturingProcess31 and BiologicalMaterial12
Unique to PLS: ManufacturingProcess17 and ManufacturingProcess34

The other difference is in emphasis. PLS puts 7 process predictors in its top 10 against the SVM’s 6, and PLS ranks the biological measurements lower (BiologicalMaterial06 is 7th for PLS but 2nd for the SVM). The broad conclusion from Exercise 6.3 survives: the controllable process settings matter most, with ManufacturingProcess32 at the head of both lists, and the raw material measurements form a secondary but real group.

c. Explore the relationships between the top predictors and the response for the predictors unique to the optimal nonlinear model

The two predictors that appear only in the SVM list are ManufacturingProcess31 and BiologicalMaterial12.

unique_preds <- c("ManufacturingProcess31", "BiologicalMaterial12")

tibble(predictor = unique_preds,
       correlation = round(as.numeric(cor(predictors_imputed[, unique_preds], yield)), 3))

## # A tibble: 2 x 2
##   predictor              correlation
##   <chr>                        <dbl>
## 1 ManufacturingProcess31      -0.071
## 2 BiologicalMaterial12         0.367

curvature <- map_dfr(unique_preds, function(v) {
  d <- tibble(x = predictors_imputed[, v], y = yield)
  linear <- lm(y ~ x, d)
  quadratic <- lm(y ~ poly(x, 2), d)
  tibble(predictor = v,
         linear_R2 = round(summary(linear)$r.squared, 3),
         quadratic_R2 = round(summary(quadratic)$r.squared, 3),
         p_value_for_curvature = signif(anova(linear, quadratic)$`Pr(>F)`[2], 3))
})

curvature

## # A tibble: 2 x 4
##   predictor              linear_R2 quadratic_R2 p_value_for_curvature
##   <chr>                      <dbl>        <dbl>                 <dbl>
## 1 ManufacturingProcess31     0.005        0.145           0.000000318
## 2 BiologicalMaterial12       0.135        0.18            0.00233

predictors_imputed[, unique_preds] |>
  as_tibble() |>
  mutate(Yield = yield) |>
  pivot_longer(-Yield, names_to = "predictor", values_to = "value") |>
  ggplot(aes(value, Yield)) +
  geom_point(alpha = 0.6, colour = "steelblue") +
  geom_smooth(method = "lm", se = FALSE, colour = "firebrick", linetype = "dashed") +
  geom_smooth(method = "loess", se = FALSE, colour = "darkgreen") +
  facet_wrap(~ predictor, scales = "free_x") +
  labs(title = "Predictors unique to the SVM model (dashed red = linear, green = loess)",
       x = "Value (centred and scaled)")

These plots explain exactly why a nonlinear model finds these two predictors and a linear one does not.

ManufacturingProcess31 is the striking case. Its correlation with yield is -0.071, essentially zero, and a straight line through it explains just 0.5% of the variation. A linear model has no reason to use it at all, which is why it does not appear in the PLS list. But adding a quadratic term raises the explained variation to 14.5%, and the improvement is highly significant (p = 0.0000003). The loess curve shows why: yield is low at the lowest values of the predictor, rises across the middle of its range, and falls away again at the top. A straight line through that hump has almost no slope, so the linear model sees nothing while a nonlinear model sees a genuine effect.

BiologicalMaterial12 shows a milder version of the same thing. Its linear correlation of 0.367 is already respectable, but the quadratic fit improves the explained variation from 13.5% to 18.0% (p = 0.0023), as the relationship flattens off at higher values.

What this tells us about the process. ManufacturingProcess31 appears to have an optimal operating range rather than a “more is better” relationship. That is a familiar pattern for physical settings such as a temperature, pressure or flow rate, where both too little and too much hurt the result. The practical implication is different from the linear advice in Exercise 6.3: rather than pushing this setting up or down, the goal would be to hold it near the middle of its range, where the curve peaks.

The flattening in BiologicalMaterial12 carries a different lesson. Raw material quality helps yield, but only up to a point, so paying a premium for material above that level would probably not be rewarded.

Both findings come with the same caution as before: this is observational data from past production runs, so these are associations rather than demonstrated causes, and a designed experiment that deliberately varies ManufacturingProcess31 across its range would be the proper way to confirm the optimum before changing production practice.