Exercises 7.2 and 7.5 from Kuhn & Johnson, Applied Predictive Modeling, Chapter 7 (Nonlinear Regression Models). Model tuning runs in parallel on 8 cores with doParallel.

Exercise 7.2 — Friedman simulation

The data come from Friedman’s benchmark (1991):

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

There are 10 uniform predictors, but only \(x_1\) to \(x_5\) enter the model; \(x_6\) to \(x_{10}\) are pure noise.

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

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

All models are tuned with 10-fold cross-validation repeated 5 times on the 200 training points, then scored on the 5,000 test points.

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

K-nearest neighbors (the book’s example)

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

## k-Nearest Neighbors 
## 
## 200 samples
##  10 predictor
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 180, 180, 180, 180, 180, 180, ... 
## Resampling results across tuning parameters:
## 
##   k   RMSE      Rsquared   MAE     
##    5  3.199670  0.5910371  2.649336
##    7  3.108444  0.6299684  2.537062
##    9  3.090673  0.6553105  2.506983
##   11  3.070640  0.6798579  2.489940
##   13  3.064958  0.6966689  2.467551
##   15  3.082000  0.7111449  2.491659
##   17  3.102544  0.7142288  2.506028
##   19  3.115994  0.7231515  2.517837
##   21  3.148223  0.7237919  2.560250
##   23  3.177131  0.7236183  2.589908
## 
## 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$x)
postResample(pred = knnPred, obs = testData$y)

##      RMSE  Rsquared       MAE 
## 3.1481557 0.6747755 2.5236041

MARS

set.seed(200)
marsModel <- train(x = trainingData$x, y = trainingData$y, method = "earth",
                   tuneGrid = expand.grid(degree = 1:2, nprune = 2:38),
                   trControl = ctrl)
marsModel$bestTune

##    nprune degree
## 50     14      2

plot(marsModel, main = "MARS: resampled RMSE")

summary(marsModel$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

Support vector machine (radial kernel)

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

##        sigma C
## 6 0.06299324 8

Neural network (model averaging)

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

##    size decay   bag
## 11    4  0.01 FALSE

Comparison

f_models <- list(KNN = knnModel, MARS = marsModel, SVM = svmModel, `Neural net` = nnetModel)

f_compare <- bind_rows(lapply(names(f_models), function(m) {
  fit <- f_models[[m]]
  cv <- getTrainPerf(fit)
  test <- postResample(predict(fit, testData$x), testData$y)
  tibble(model = m,
         tuning = paste(names(fit$bestTune), signif(unlist(fit$bestTune), 3),
                        sep = " = ", collapse = ", "),
         cv_RMSE = cv$TrainRMSE, cv_R2 = cv$TrainRsquared,
         test_RMSE = test[["RMSE"]], test_R2 = test[["Rsquared"]])
})) |>
  arrange(test_RMSE)
knitr::kable(f_compare, digits = 3)

model	tuning	cv_RMSE	cv_R2	test_RMSE	test_R2
MARS	nprune = 14, degree = 2	1.241	0.937	1.172	0.945
SVM	sigma = 0.063, C = 8	1.868	0.857	2.054	0.829
Neural net	size = 4, decay = 0.01, bag = 0	2.018	0.834	2.178	0.815
KNN	k = 13	3.065	0.697	3.148	0.675

bwplot(resamples(f_models), metric = "RMSE", main = "Resampled RMSE")

MARS gives the best performance by a wide margin: test RMSE 1.17 and \(R^2\) 0.945 on the 5,000 test points, against 2.05 for the SVM, 2.18 for the neural network, and 3.15 for KNN. The resampled ranking is the same, so the test set confirms what cross-validation picked. The tuning chose a second-degree MARS model, which makes sense: the \(10\sin(\pi x_1 x_2)\) term is an interaction between \(x_1\) and \(x_2\), and degree-2 hinge products can approximate it directly. The other terms are additive (a quadratic in \(x_3\) and linear in \(x_4\) and \(x_5\)), which hinge functions handle well. KNN does worst because all 10 predictors count equally in its distance, so the 5 noise predictors dilute the neighborhoods.

Does MARS select the informative predictors?

mars_imp <- varImp(marsModel)$importance |>
  rownames_to_column("predictor") |>
  arrange(desc(Overall))
knitr::kable(mars_imp, digits = 1, caption = "MARS variable importance (GCV)")

MARS variable importance (GCV)
predictor	Overall
X1	100.0
X4	75.2
X2	48.7
X5	15.5
X3	0.0

fm <- marsModel$finalModel
mars_used <- colnames(fm$dirs)[colSums(abs(fm$dirs[fm$selected.terms, , drop = FALSE])) > 0]
mars_used

## [1] "X1" "X2" "X3" "X4" "X5"

Yes. The final MARS model uses exactly X1, X2, X3, X4, X5, the 5 informative predictors, and none of \(x_6\) to \(x_{10}\). In the importance table \(x_3\) shows 0 even though it appears in 2 hinge terms, because removing it barely changes the GCV score. Its effect, \(20(x_3 - 0.5)^2\), ranges only from 0 to 5 over \(x_3\) in \([0, 1]\), small next to the other terms. The 2 interaction terms that pair \(x_1\) with \(x_2\) match the true sine interaction.

Exercise 7.5 — Chemical manufacturing, nonlinear models

The data, split, and pre-processing are the same as in Exercise 6.3 (Assignment 7): the same seed and createDataPartition call give the same 80/20 split, missing values are filled with KNN imputation learned on the training set, all predictors are centered and scaled, and the one near-zero-variance predictor is dropped.

data(ChemicalManufacturingProcess)
x_chem <- ChemicalManufacturingProcess |> select(-Yield)
y_chem <- ChemicalManufacturingProcess$Yield

set.seed(624)
in_train_c <- createDataPartition(y_chem, p = 0.8, list = FALSE)
pp <- preProcess(x_chem[in_train_c, ], method = c("knnImpute", "nzv"))
xc_train <- predict(pp, x_chem[in_train_c, ]);  yc_train <- y_chem[in_train_c]
xc_test  <- predict(pp, x_chem[-in_train_c, ]); yc_test  <- y_chem[-in_train_c]
c(train = nrow(xc_train), test = nrow(xc_test), predictors = ncol(xc_train))

##      train       test predictors 
##        144         32         56

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

set.seed(624)
chem_knn <- train(xc_train, yc_train, method = "knn",
                  tuneGrid = data.frame(k = seq(1, 25, 2)), trControl = ctrl_c)

set.seed(624)
chem_mars <- train(xc_train, yc_train, method = "earth",
                   tuneGrid = expand.grid(degree = 1:2, nprune = 2:30),
                   trControl = ctrl_c)

set.seed(624)
chem_svm <- train(xc_train, yc_train, method = "svmRadial",
                  tuneLength = 14, trControl = ctrl_c)

set.seed(624)
chem_nnet <- train(xc_train, yc_train, method = "avNNet",
                   tuneGrid = expand.grid(size = 1:5, decay = c(0.01, 0.1, 0.5, 1), bag = FALSE),
                   linout = TRUE, trace = FALSE, maxit = 500,
                   MaxNWts = 5 * (ncol(xc_train) + 1) + 5 + 1,
                   trControl = ctrl_c)

# The best linear model from Exercise 6.3 (lasso), refit for comparison.
set.seed(624)
chem_lasso <- train(xc_train, yc_train, method = "glmnet",
                    tuneGrid = expand.grid(alpha = c(0.1, 0.5, 1),
                                           lambda = 10^seq(-3, 0, length.out = 20)),
                    trControl = ctrl_c)

(a) Which nonlinear model is best?

chem_models <- list(KNN = chem_knn, MARS = chem_mars, SVM = chem_svm,
                    `Neural net` = chem_nnet, `Lasso (linear, 6.3)` = chem_lasso)

chem_compare <- bind_rows(lapply(names(chem_models), function(m) {
  fit <- chem_models[[m]]
  cv <- getTrainPerf(fit)
  test <- postResample(predict(fit, xc_test), yc_test)
  tibble(model = m,
         tuning = paste(names(fit$bestTune), signif(unlist(fit$bestTune), 3),
                        sep = " = ", collapse = ", "),
         cv_RMSE = cv$TrainRMSE, cv_R2 = cv$TrainRsquared,
         test_RMSE = test[["RMSE"]], test_R2 = test[["Rsquared"]])
}))
knitr::kable(chem_compare, digits = 3)

model	tuning	cv_RMSE	cv_R2	test_RMSE	test_R2
KNN	k = 3	1.285	0.505	1.338	0.650
MARS	nprune = 4, degree = 1	1.218	0.571	1.201	0.713
SVM	sigma = 0.014, C = 16	1.149	0.611	1.149	0.747
Neural net	size = 2, decay = 0.01, bag = 0	1.408	0.483	1.441	0.564
Lasso (linear, 6.3)	alpha = 1, lambda = 0.0785	1.167	0.590	1.171	0.729

bwplot(resamples(chem_models), metric = "RMSE", main = "Resampled RMSE")

nonlinear <- chem_compare |> filter(model != "Lasso (linear, 6.3)")
best_name <- nonlinear$model[which.min(nonlinear$cv_RMSE)]
chem_best <- chem_models[[best_name]]
best_name

## [1] "SVM"

The radial SVM is the best nonlinear model on both counts: resampled RMSE 1.149 (\(R^2\) 0.611) and test RMSE 1.149 (\(R^2\) 0.747). MARS is second; its final model is very small (2 predictors, 3 hinge terms). The neural network does worst, likely because 144 training runs are too few to fit its weights reliably even with weight decay, and KNN suffers from 56 mostly uninformative predictors. The SVM only narrowly beats the lasso from Exercise 6.3 (test RMSE 1.171), so the yield relationship is mostly linear and the nonlinear model adds little.

(b) Most important predictors

importance_table <- function(fit) {
  varImp(fit)$importance |>
    rownames_to_column("predictor") |>
    select(predictor, importance = 2) |>
    arrange(desc(importance)) |>
    mutate(type = ifelse(grepl("^Biological", predictor), "Biological", "Process"))
}

imp_nl  <- importance_table(chem_best)
imp_lin <- importance_table(chem_lasso)

imp_nl |>
  slice_head(n = 10) |>
  ggplot(aes(importance, reorder(predictor, importance), fill = type)) +
  geom_col() +
  labs(title = paste("Top 10 predictors:", best_name), x = "Scaled importance", y = NULL)

table(head(imp_nl$type, 10))

## 
## Biological    Process 
##          4          6

top10 <- tibble(
  rank = 1:10,
  nonlinear = head(imp_nl$predictor, 10),
  linear_lasso = head(imp_lin$predictor, 10)
) |>
  mutate(nonlinear_in_linear_top10 = nonlinear %in% linear_lasso)
knitr::kable(top10, caption = paste("Top 10:", best_name, "vs lasso"))

Top 10: SVM vs lasso
rank	nonlinear	linear_lasso	nonlinear_in_linear_top10
1	ManufacturingProcess32	ManufacturingProcess32	TRUE
2	BiologicalMaterial06	ManufacturingProcess09	TRUE
3	BiologicalMaterial03	ManufacturingProcess17	FALSE
4	ManufacturingProcess36	BiologicalMaterial06	TRUE
5	ManufacturingProcess09	ManufacturingProcess34	TRUE
6	ManufacturingProcess13	ManufacturingProcess37	FALSE
7	ManufacturingProcess31	ManufacturingProcess36	FALSE
8	BiologicalMaterial02	ManufacturingProcess06	FALSE
9	BiologicalMaterial12	ManufacturingProcess04	FALSE
10	ManufacturingProcess06	ManufacturingProcess39	TRUE

unique_nl <- setdiff(top10$nonlinear, top10$linear_lasso)
unique_nl

## [1] "BiologicalMaterial03"   "ManufacturingProcess13" "ManufacturingProcess31"
## [4] "BiologicalMaterial02"   "BiologicalMaterial12"

Caret has no model-specific importance for SVMs, so varImp falls back to a model-free measure: for each predictor, the \(R^2\) of a loess smoother fit of yield on that predictor alone. The lasso ranking is the absolute size of its coefficients. The 2 lists measure different things, so they are compared by rank only.

Process variables still lead, but less heavily. The nonlinear top 10 has 6 process and 4 biological predictors, while the lasso top 10 has 9 process and 1 biological. ManufacturingProcess32 is first in both. Of the nonlinear top 10, 5 also appear in the lasso top 10 (ManufacturingProcess32, 09, 36, 06 and BiologicalMaterial06). The biological predictors rank higher in the nonlinear list because they are strongly correlated with each other and with yield: the lasso keeps one of a correlated group and shrinks the rest to zero, while the per-predictor loess measure gives each of them credit.

(c) Predictors unique to the nonlinear model

xc_train |>
  select(all_of(unique_nl)) |>
  mutate(Yield = yc_train) |>
  pivot_longer(-Yield, names_to = "predictor") |>
  ggplot(aes(value, Yield)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, span = 0.8) +
  facet_wrap(~ predictor, scales = "free_x") +
  labs(title = "Yield vs predictors unique to the nonlinear top 10 (training set, scaled)",
       x = "Scaled value")

knitr::kable(tibble(predictor = unique_nl,
                    cor_with_yield = sapply(unique_nl, function(v) cor(xc_train[[v]], yc_train))),
             digits = 3)

predictor	cor_with_yield
BiologicalMaterial03	0.451
ManufacturingProcess13	-0.446
ManufacturingProcess31	-0.063
BiologicalMaterial02	0.477
BiologicalMaterial12	0.351

The 5 predictors in the nonlinear top 10 but not the lasso top 10 are BiologicalMaterial03, ManufacturingProcess13, ManufacturingProcess31, BiologicalMaterial02, BiologicalMaterial12.

BiologicalMaterial02, 03, and 12 all rise with yield (correlations of about 0.35 to 0.48). BiologicalMaterial02 rises steadily, while 03 and 12 level off and dip slightly at their highest values, a bend a linear model cannot follow. They are also strongly correlated with each other and with BiologicalMaterial06 (pairwise correlations of 0.70 to 0.95). The lasso did not drop them for being unrelated; it dropped them because they carry the same information as BiologicalMaterial06. They describe the quality of the raw material, which is fixed before a run starts, so they help predict yield but cannot be tuned.
ManufacturingProcess13 falls as yield rises (correlation about \(-0.45\)), with the steepest drop at low values, so lower settings look better.
ManufacturingProcess31 has almost no linear correlation with yield (about \(-0.06\)), which is why the lasso ignores it. Almost all of its values sit in a tight cluster near 0, with a single run far below at about \(-11\). That one point bends the loess curve and inflates its importance; within the main cluster, yield drops as the value rises, but over a very narrow range. This predictor’s high rank is mostly an artifact of 1 run.

For improving yield, the takeaway is the same as in 6.3: the controllable levers are the process settings, with ManufacturingProcess32 and 09 the most consistent across both model types. ManufacturingProcess13 is a new candidate from the nonlinear view; the biological predictors are useful for forecasting the yield of a batch, not for changing it.

DATA 624 — Assignment 8 (Applied Predictive Modeling, Chapter 7)

Cailinn

2026-10-27