DATA 624 — Assignment 9 (Applied Predictive Modeling, Chapter 8)
Cailinn
2026-11-04
Exercises 8.1, 8.2, 8.3 and 8.7 from Kuhn & Johnson, Applied
Predictive Modeling, Chapter 8 (Regression Trees and Rule-Based
Models). Model tuning runs in parallel on 8 cores with
doParallel.
Exercise 8.1 — Friedman data and variable importance
set.seed(200)
simulated <- mlbench.friedman1(200, sd = 1)
simulated <- cbind(simulated$x, simulated$y)
simulated <- as.data.frame(simulated)
colnames(simulated)[ncol(simulated)] <- "y"\(V1\) to \(V5\) enter the true model, and \(V6\) to \(V10\) are noise.
(a) Random forest importance
set.seed(200)
model1 <- randomForest(y ~ ., data = simulated, importance = TRUE, ntree = 1000)
rfImp1 <- varImp(model1, scale = FALSE)
rfImp1 |> rownames_to_column("predictor") |> arrange(desc(Overall)) |>
knitr::kable(digits = 3, caption = "Random forest permutation importance")| predictor | Overall |
|---|---|
| V1 | 8.679 |
| V4 | 7.782 |
| V2 | 6.489 |
| V5 | 2.393 |
| V3 | 0.771 |
| V6 | 0.015 |
| V7 | 0.007 |
| V9 | -0.027 |
| V10 | -0.039 |
| V8 | -0.109 |
No. The 5 informative predictors have the largest importance (V1, V4 and V2 lead, then V5 and V3), while V6 to V10 all have importance close to 0 (between -0.11 and 0.01, some negative). Permuting a noise predictor does not hurt the forest’s predictions, so it did not rely on them. V3 has a low score as well; its true effect, \(20(x_3 - 0.5)^2\), is symmetric around 0.5 and small next to the other terms.
(c) Conditional inference forest (cforest)
cf_imp <- function(data, conditional) {
set.seed(200)
fit <- party::cforest(y ~ ., data = data,
controls = party::cforest_unbiased(ntree = 1000, mtry = 3))
party::varimp(fit, conditional = conditional)
}
base_cols <- c(paste0("V", 1:10), "y")
cf <- bind_rows(
tibble(set = "original", predictor = names(cf_imp(simulated[base_cols], FALSE)),
traditional = cf_imp(simulated[base_cols], FALSE),
conditional = cf_imp(simulated[base_cols], TRUE)),
tibble(set = "+ duplicates", predictor = names(cf_imp(simulated, FALSE)),
traditional = cf_imp(simulated, FALSE),
conditional = cf_imp(simulated, TRUE))
)
cf |>
pivot_wider(names_from = set, values_from = c(traditional, conditional)) |>
knitr::kable(digits = 3, caption = "cforest importance, traditional and conditional")| predictor | traditional_original | traditional_+ duplicates | conditional_original | conditional_+ duplicates |
|---|---|---|---|---|
| V1 | 7.209 | 4.502 | 4.195 | 1.643 |
| V2 | 5.848 | 4.722 | 4.437 | 3.563 |
| V3 | 0.142 | 0.095 | 0.121 | 0.067 |
| V4 | 6.823 | 5.672 | 5.063 | 4.190 |
| V5 | 2.028 | 1.578 | 1.227 | 1.194 |
| V6 | -0.014 | -0.005 | 0.001 | 0.028 |
| V7 | 0.023 | -0.064 | 0.047 | -0.024 |
| V8 | -0.072 | -0.023 | -0.048 | -0.014 |
| V9 | -0.001 | -0.064 | 0.001 | 0.018 |
| V10 | -0.020 | -0.098 | -0.033 | -0.042 |
| duplicate1 | NA | 2.433 | NA | 0.954 |
| duplicate2 | NA | 2.295 | NA | 0.696 |
The traditional cforest importance
shows the same pattern as the random forest: V6 to V10 are near 0, and
adding the 2 copies lowers V1 from about 7.2 to 4.5 while the copies
take about 2.3 to 2.4 each.
The conditional importance (Strobl et al.) permutes each predictor within groups of similar values of its correlated predictors, so it measures what a predictor adds beyond them. It is smaller for V1 even without the copies, and it reacts more strongly to them: V1 drops from about 4.2 to 1.6, and the copies get only about 0.7 to 1.0. In other words, once V1 has near-identical copies, none of the 3 adds much that the others do not already carry. The noise predictors stay near 0 under both measures, so the pattern is the same; conditional importance just reflects the redundancy more honestly.
(d) Boosted trees and Cubist
fit_imp <- function(data) {
set.seed(200)
gb <- gbm(y ~ ., data = data, distribution = "gaussian", n.trees = 1000,
interaction.depth = 3, shrinkage = 0.01, bag.fraction = 0.5)
gb_imp <- summary(gb, plotit = FALSE) |> as_tibble() |> select(predictor = var, gbm = rel.inf)
cu <- cubist(x = data |> select(-y), y = data$y, committees = 100)
cu_imp <- varImp(cu) |> rownames_to_column("predictor") |> rename(cubist = Overall)
full_join(gb_imp, cu_imp, by = "predictor")
}
bind_rows(
fit_imp(simulated[base_cols]) |> mutate(set = "original"),
fit_imp(simulated[c(base_cols, "duplicate1")]) |> mutate(set = "+ duplicate1"),
fit_imp(simulated) |> mutate(set = "+ duplicates")
) |>
filter(predictor %in% c("V1", "V2", "V3", "V4", "V5", "duplicate1", "duplicate2")) |>
pivot_wider(names_from = set, values_from = c(gbm, cubist)) |>
knitr::kable(digits = 1, caption = "Boosted trees (relative influence, %) and Cubist (usage, %)")| predictor | gbm_original | gbm_+ duplicate1 | gbm_+ duplicates | cubist_original | cubist_+ duplicate1 | cubist_+ duplicates |
|---|---|---|---|---|---|---|
| V4 | 28.4 | 28.3 | 27.8 | 48.0 | 48.0 | 49.0 |
| V1 | 26.6 | 20.3 | 18.4 | 71.5 | 71.0 | 48.0 |
| V2 | 22.0 | 21.7 | 21.5 | 58.5 | 59.0 | 55.0 |
| V5 | 11.7 | 11.7 | 12.0 | 33.0 | 32.5 | 40.5 |
| V3 | 7.7 | 7.4 | 7.2 | 47.0 | 46.5 | 45.0 |
| duplicate1 | NA | 7.2 | 6.6 | NA | 0.0 | 7.0 |
| duplicate2 | NA | NA | 3.1 | NA | NA | 19.5 |
Boosted trees show the same pattern, more mildly. V1’s relative influence falls from about 27% to 20% with 1 copy and 18% with 2, and the copies take about 7% and 3%. Boosting grows each tree on the residuals of the previous ones, so once V1 has been used the copy has less left to explain, and V1 keeps more of its share than in the random forest.
Cubist behaves differently. Its importance is the percentage of rules and linear models that use each predictor. With 1 copy, Cubist ignores it completely (0%) and V1 stays at about 71%. With 2 copies, V1 drops to about 48% and the copies are used in 7% and 20% of the rules. All 3 methods keep the noise predictors near 0, apart from Cubist using V6 in about 12% of rules once the copies are added.
Exercise 8.2 — Tree bias toward high-granularity predictors
Simulation: the response is pure noise, and 3 predictors that are all unrelated to it differ only in how many distinct values they take. An unbiased tree would choose each one for the first split about equally often. A second scenario gives the 2-level predictor a real but weak effect.
root_var <- function(effect) {
n <- 200
d <- tibble(
many = runif(n), # about 200 distinct values
some = sample(1:10, n, replace = TRUE) / 10, # 10 distinct values
two = sample(0:1, n, replace = TRUE), # 2 distinct values
y = effect * two + rnorm(n)
)
fit <- rpart(y ~ ., data = d, control = rpart.control(maxdepth = 1, cp = 0, minsplit = 2))
as.character(fit$frame$var[1])
}
set.seed(624)
sims <- bind_rows(
tibble(scenario = "y unrelated to all 3", root = replicate(1000, root_var(0))),
tibble(scenario = "y depends weakly on 'two'", root = replicate(1000, root_var(0.3)))
)
sims |>
count(scenario, root) |>
group_by(scenario) |>
mutate(share = n / sum(n)) |>
knitr::kable(digits = 3, caption = "Predictor chosen for the first split, 1,000 simulations each")| scenario | root | n | share |
|---|---|---|---|
| y depends weakly on ‘two’ | many | 443 | 0.443 |
| y depends weakly on ‘two’ | some | 128 | 0.128 |
| y depends weakly on ‘two’ | two | 429 | 0.429 |
| y unrelated to all 3 | many | 761 | 0.761 |
| y unrelated to all 3 | some | 200 | 0.200 |
| y unrelated to all 3 | two | 39 | 0.039 |
sims |>
count(scenario, root) |>
ggplot(aes(root, n, fill = root)) +
geom_col(show.legend = FALSE) +
facet_wrap(~ scenario) +
labs(title = "First-split variable chosen by rpart", x = NULL, y = "Simulations (of 1,000)")
When none of the predictors is related to the
response, an unbiased tree would pick each of the 3 for the
first split about a third of the time. Instead, rpart picks
the predictor with about 200 distinct values in 76% of simulations, the
10-value predictor in 20%, and the 2-value predictor in only 4%. A
predictor with more distinct values offers more candidate split points,
so by chance alone one of them reduces the error more.
When the 2-value predictor has a real but weak
effect, it still wins only 43% of the time, about the same as
the pure-noise continuous predictor (44%). That is the selection bias
described in Section 8.1: trees favor high-granularity predictors, so
tree-based importance can rank a noisy continuous predictor above a
truly useful categorical one. Conditional inference trees
(ctree, cforest) were designed to remove this
bias.
Exercise 8.3 — Boosting: bagging fraction and learning rate
Figure 8.24 compares boosted trees on the solubility data with bagging fraction and learning rate both at 0.1 (left) and both at 0.9 (right). The importance plots are reproduced below, along with the resampled RMSE of each setting.
data(solubility)
gbm_sol <- function(rate, frac, depth = 3) {
set.seed(624)
gbm(y ~ ., data = cbind(solTrainXtrans, y = solTrainY), distribution = "gaussian",
n.trees = 100, interaction.depth = depth, shrinkage = rate, bag.fraction = frac,
cv.folds = 5)
}
g01 <- gbm_sol(0.1, 0.1)
g09 <- gbm_sol(0.9, 0.9)
g01_d10 <- gbm_sol(0.1, 0.1, depth = 10)
top25 <- function(fit, label) {
summary(fit, plotit = FALSE) |> as_tibble() |> slice_head(n = 25) |>
mutate(rank = row_number(), setting = label)
}
imp_sol <- bind_rows(top25(g01, "rate 0.1, fraction 0.1"),
top25(g09, "rate 0.9, fraction 0.9"))
tibble(setting = c("rate 0.1, fraction 0.1", "rate 0.9, fraction 0.9"),
cv_RMSE = sqrt(c(min(g01$cv.error), min(g09$cv.error))),
test_RMSE = c(RMSE(predict(g01, solTestXtrans, n.trees = gbm.perf(g01, plot.it = FALSE)), solTestY),
RMSE(predict(g09, solTestXtrans, n.trees = gbm.perf(g09, plot.it = FALSE)), solTestY)),
top_predictor_share = c(imp_sol$rel.inf[imp_sol$setting == "rate 0.1, fraction 0.1"][1],
imp_sol$rel.inf[imp_sol$setting == "rate 0.9, fraction 0.9"][1])) |>
knitr::kable(digits = 3)| setting | cv_RMSE | test_RMSE | top_predictor_share |
|---|---|---|---|
| rate 0.1, fraction 0.1 | 0.738 | 0.680 | 27.277 |
| rate 0.9, fraction 0.9 | 0.777 | 0.719 | 39.075 |
ggplot(imp_sol, aes(rank, rel.inf, colour = setting)) +
geom_line() + geom_point(size = 1.5) +
labs(title = "Relative influence of the top 25 predictors", x = "Rank", y = "Relative influence (%)")
(a) The learning rate controls how much of each new tree is added to the model. At 0.9, the first few trees fit most of the signal, and they split on the strongest predictors (molecular weight, number of carbons, surface area). Later trees only fit what is left over, so they rarely need other predictors. The bagging fraction of 0.9 also means that every tree sees almost the same data, so successive trees make similar choices. Together, importance concentrates on a few predictors: here the top predictor takes 39.1% of the influence and the top 5 take 82.4%. With rate and fraction 0.1, each tree adds a small step fit to a different small subsample, so many predictors get a chance to contribute (top predictor 27.3%, top 5 61.1%).
(b) The model on the left (0.1 and 0.1) should predict other samples better. A small learning rate is the usual regularization for boosting (Friedman’s shrinkage), and subsampling adds randomness that reduces variance. A large learning rate tends to overfit. Here the cross-validated RMSE agrees: 0.738 for the left model vs 0.777 for the right. The test set agrees as well (see the table above). With only 100 trees the gap is modest.
(c) Increasing the interaction depth makes each tree larger, so each tree uses more predictors in its splits. Importance is then spread across more predictors, which flattens the slope of the importance plot. Here the effect is modest: with depth 10 instead of 3 (rate and fraction 0.1), the top predictor’s share falls from 27.3% to 24.9%, and the 25th predictor’s rises from 0.54% to 0.57%.
Exercise 8.7 — Chemical manufacturing, tree-based models
Same data, split and pre-processing as Exercises 6.3 and 7.5 (Assignments 7 and 8): KNN imputation and near-zero-variance filtering learned on the training set, the same 80/20 split, and 10-fold cross-validation repeated 5 times.
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]
ctrl_c <- trainControl(method = "repeatedcv", number = 10, repeats = 5)set.seed(624)
chem_rpart <- train(xc_train, yc_train, method = "rpart2",
tuneGrid = data.frame(maxdepth = 1:10), trControl = ctrl_c)
set.seed(624)
chem_bag <- train(xc_train, yc_train, method = "treebag", nbagg = 50, trControl = ctrl_c)
set.seed(624)
chem_rf <- train(xc_train, yc_train, method = "rf", ntree = 500, importance = TRUE,
tuneGrid = data.frame(mtry = c(5, 10, 15, 20, 30, 40, 56)), trControl = ctrl_c)
set.seed(624)
chem_gbm <- train(xc_train, yc_train, method = "gbm", verbose = FALSE,
tuneGrid = expand.grid(n.trees = c(100, 300, 500, 1000),
interaction.depth = c(1, 3, 5),
shrinkage = c(0.01, 0.1),
n.minobsinnode = 10),
trControl = ctrl_c)
set.seed(624)
chem_cubist <- train(xc_train, yc_train, method = "cubist",
tuneGrid = expand.grid(committees = c(1, 10, 50, 100),
neighbors = c(0, 1, 5, 9)),
trControl = ctrl_c)
# The best linear (6.3) and non-linear (7.5) models, 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)
set.seed(624)
chem_svm <- train(xc_train, yc_train, method = "svmRadial", tuneLength = 14, trControl = ctrl_c)(a) Which tree-based model is best?
chem_models <- list(`Single tree` = chem_rpart, `Bagged trees` = chem_bag,
`Random forest` = chem_rf, `Boosted trees` = chem_gbm,
Cubist = chem_cubist,
`Lasso (6.3)` = chem_lasso, `SVM (7.5)` = chem_svm)
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)
tune <- fit$bestTune[, !names(fit$bestTune) %in% c("parameter", "n.minobsinnode"), drop = FALSE]
tibble(model = m,
tuning = if (ncol(tune) == 0) "none" else
paste(names(tune), signif(unlist(tune), 3), sep = " = ", collapse = ", ") |>
sub("interaction.depth", "depth", x = _) |> sub("shrinkage", "rate", x = _),
cv_RMSE = cv$TrainRMSE, cv_R2 = cv$TrainRsquared,
test_RMSE = test[["RMSE"]], test_R2 = test[["Rsquared"]])
}))
knitr::kable(chem_compare, digits = 3,
col.names = c("Model", "Tuning", "CV RMSE", "CV R2", "Test RMSE", "Test R2"))| Model | Tuning | CV RMSE | CV R2 | Test RMSE | Test R2 |
|---|---|---|---|---|---|
| Single tree | maxdepth = 2 | 1.450 | 0.381 | 1.429 | 0.569 |
| Bagged trees | none | 1.219 | 0.578 | 1.238 | 0.720 |
| Random forest | mtry = 15 | 1.137 | 0.654 | 1.256 | 0.728 |
| Boosted trees | n.trees = 1000, depth = 5, rate = 0.1 | 1.089 | 0.658 | 1.151 | 0.737 |
| Cubist | committees = 100, neighbors = 1 | 0.932 | 0.734 | 0.882 | 0.834 |
| Lasso (6.3) | alpha = 1, lambda = 0.0785 | 1.167 | 0.590 | 1.171 | 0.729 |
| SVM (7.5) | sigma = 0.014, C = 16 | 1.149 | 0.611 | 1.149 | 0.747 |
bwplot(resamples(chem_models), metric = "RMSE", main = "Resampled RMSE")
trees <- chem_compare |> filter(!grepl("6.3|7.5", model))
best_tree_name <- trees$model[which.min(trees$cv_RMSE)]
chem_best <- chem_models[[best_tree_name]]
best_tree_name## [1] "Cubist"Cubist gives the best performance of the tree-based models, both in resampling (RMSE 0.932, \(R^2\) 0.734) and on the test set (RMSE 0.882, \(R^2\) 0.834), with committees = 100, neighbors = 1. Boosted trees are second and the random forest third. The single tree is the weakest model in the table, which is the usual pattern: averaging or boosting many trees reduces the high variance of one tree.
Cubist is also clearly better than the best linear model from 6.3 (lasso, test RMSE 1.171) and the best non-linear model from 7.5 (SVM, 1.149). Cubist fits a linear model inside each rule, so it can be linear within a region and still change slope between regions, which suits these data.
(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_tree <- importance_table(chem_best)
imp_lin <- importance_table(chem_lasso)
imp_nl <- importance_table(chem_svm)
imp_tree |>
slice_head(n = 10) |>
ggplot(aes(importance, reorder(predictor, importance), fill = type)) +
geom_col() +
labs(title = paste("Top 10 predictors:", best_tree_name), x = "Scaled importance", y = NULL)
top10 <- tibble(rank = 1:10,
tree = head(imp_tree$predictor, 10),
linear_lasso = head(imp_lin$predictor, 10),
nonlinear_svm = head(imp_nl$predictor, 10))
knitr::kable(top10, caption = "Top 10 predictors by model")| rank | tree | linear_lasso | nonlinear_svm |
|---|---|---|---|
| 1 | ManufacturingProcess32 | ManufacturingProcess32 | ManufacturingProcess32 |
| 2 | ManufacturingProcess17 | ManufacturingProcess09 | BiologicalMaterial06 |
| 3 | ManufacturingProcess09 | ManufacturingProcess17 | BiologicalMaterial03 |
| 4 | ManufacturingProcess33 | BiologicalMaterial06 | ManufacturingProcess36 |
| 5 | ManufacturingProcess01 | ManufacturingProcess34 | ManufacturingProcess09 |
| 6 | ManufacturingProcess37 | ManufacturingProcess37 | ManufacturingProcess13 |
| 7 | ManufacturingProcess45 | ManufacturingProcess36 | ManufacturingProcess31 |
| 8 | ManufacturingProcess15 | ManufacturingProcess06 | BiologicalMaterial02 |
| 9 | ManufacturingProcess21 | ManufacturingProcess04 | BiologicalMaterial12 |
| 10 | BiologicalMaterial02 | ManufacturingProcess39 | ManufacturingProcess06 |
tibble(model = c(best_tree_name, "Lasso", "SVM"),
process_in_top10 = c(sum(head(imp_tree$type, 10) == "Process"),
sum(head(imp_lin$type, 10) == "Process"),
sum(head(imp_nl$type, 10) == "Process")),
shared_with_tree = c(10, length(intersect(top10$tree, top10$linear_lasso)),
length(intersect(top10$tree, top10$nonlinear_svm)))) |>
knitr::kable()| model | process_in_top10 | shared_with_tree |
|---|---|---|
| Cubist | 9 | 10 |
| Lasso | 9 | 4 |
| SVM | 6 | 3 |
Cubist’s importance is the percentage of rule conditions and linear models that use each predictor. Process variables dominate: 9 of the Cubist top 10 are process variables, the same as the lasso and more than the SVM (counts: 9 for the lasso, 6 for the SVM). ManufacturingProcess32 is first in all 3 models, and it is used far more often than anything else in Cubist (twice the usage of the next predictors). The Cubist top 10 shares 4 predictors with the lasso top 10 (ManufacturingProcess32, 17, 09 and 37) and 3 with the SVM top 10 (ManufacturingProcess32, 09 and BiologicalMaterial02). Several process variables appear only in the Cubist list (ManufacturingProcess33, 01, 45, 15 and 21), so the rule-based model finds uses for process settings that the other models did not rank highly. The many tied values (for example ManufacturingProcess17 and 09) are a feature of usage-based importance and should not be over-read.
(c) The optimal single tree
chem_rpart$bestTune## maxdepth
## 2 2plot(as.party(chem_rpart$finalModel), gp = gpar(fontsize = 8),
main = "Optimal single regression tree; boxes show yield in each terminal node")
Resampling selects a tree of depth 2 with 4 terminal nodes. It first splits on ManufacturingProcess32, the most important predictor in every model so far: runs with higher (scaled) values have higher yield. Within each branch, a biological material decides the rest: BiologicalMaterial12 among the low-ManufacturingProcess32 runs and BiologicalMaterial03 among the high ones. Yield rises from left to right across the terminal nodes, from about 38.3 in the worst node to about 42.6 in the best.
Yes, the tree adds useful knowledge, even though it predicts poorly. It shows an interaction the other models only imply: ManufacturingProcess32 is the main lever, and the quality of the starting material matters most once that setting is known. The node boxplots show the overlap between groups and the spread within each group, which explains why one tree cannot match the ensembles. A practical reading is that the highest yields come from a high ManufacturingProcess32 setting combined with good BiologicalMaterial03, so ManufacturingProcess32 is the setting to target and raw-material quality decides how much it can deliver.