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

Exercises 6.2 and 6.3 from Kuhn & Johnson, Applied Predictive Modeling, Chapter 6 (Linear Regression and Its Cousins).

---

Exercise 6.2 — Permeability

(a) Load the data

data(permeability)
dim(fingerprints)

## [1]  165 1107

summary(permeability[, 1])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.06    1.55    4.91   12.24   15.47   55.60

fingerprints has 165 compounds and 1107 binary substructure predictors; permeability is the response.

(b) Remove near-zero-variance predictors

nzv <- nearZeroVar(fingerprints)
fp <- fingerprints[, -nzv]
c(removed = length(nzv), remaining = ncol(fp))

##   removed remaining 
##       719       388

nearZeroVar removes 719 sparse fingerprints, leaving 388 predictors for modeling.

(c) Split, pre-process, and tune a PLS model

set.seed(624)
perm <- permeability[, 1]
in_train <- createDataPartition(perm, p = 0.8, list = FALSE)
x_train <- fp[in_train, ];  y_train <- perm[in_train]
x_test  <- fp[-in_train, ]; y_test  <- perm[-in_train]

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

set.seed(624)
pls_fit <- train(
  x_train, y_train, method = "pls",
  preProcess = c("center", "scale"),
  tuneLength = 20, trControl = ctrl
)
pls_fit$bestTune

##   ncomp
## 9     9

best_pls <- pls_fit$results[pls_fit$results$ncomp == pls_fit$bestTune$ncomp, ]
best_pls[, c("ncomp", "RMSE", "Rsquared", "MAE")]

##   ncomp     RMSE  Rsquared      MAE
## 9     9 11.82483 0.4737577 9.141689

plot(pls_fit, metric = "Rsquared", main = "PLS: resampled R-squared by components")

Resampling (10-fold CV, repeated 5 times) picks 9 latent variables, with a resampled \(R^2\) of 0.474 (RMSE 11.82). Past that point extra components start to fit noise and the resampled \(R^2\) levels off and drops.

(d) Predict the test set

pls_pred <- predict(pls_fit, x_test)
postResample(pls_pred, y_test)

##      RMSE  Rsquared       MAE 
## 10.656898  0.530305  8.388575

The test-set \(R^2\) is 0.53 (RMSE 10.66), close to the resampled estimate. The test set has only 32 compounds, so it is a noisy check and the resampled number is the more reliable one.

(e) Other models

set.seed(624)
ridge_fit <- train(
  x_train, y_train, method = "ridge",
  preProcess = c("center", "scale"),
  tuneGrid = data.frame(lambda = seq(0, 0.2, length.out = 11)),
  trControl = ctrl
)

set.seed(624)
enet_fit <- train(
  x_train, y_train, method = "glmnet",
  preProcess = c("center", "scale"),
  tuneGrid = expand.grid(alpha = c(0.1, 0.5, 1),
                         lambda = 10^seq(-3, 0.5, length.out = 20)),
  trControl = ctrl
)

set.seed(624)
pcr_fit <- train(
  x_train, y_train, method = "pcr",
  preProcess = c("center", "scale"),
  tuneLength = 30, trControl = ctrl
)

perm_models <- list(PLS = pls_fit, PCR = pcr_fit, Ridge = ridge_fit,
                    `Elastic net` = enet_fit)

perm_compare <- bind_rows(lapply(names(perm_models), function(m) {
  fit <- perm_models[[m]]
  cv <- getTrainPerf(fit)
  test <- postResample(predict(fit, x_test), y_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(perm_compare, digits = 3)

model	tuning	cv_RMSE	cv_R2	test_RMSE	test_R2
PLS	ncomp = 9	11.825	0.474	10.657	0.530
PCR	ncomp = 28	12.544	0.425	9.931	0.564
Ridge	lambda = 0.12	11.869	0.483	11.855	0.468
Elastic net	alpha = 0.1, lambda = 1.35	11.360	0.483	11.156	0.475

All four models land in the same range: resampled \(R^2\) between about 0.42 and 0.48 and RMSE between about 11.4 and 12.5. Elastic net has the lowest resampled RMSE and keeps 202 of the 388 fingerprints. PCR does best on the test set, but it needs 28 components and has the worst resampled RMSE, so its test result is most likely luck on 32 compounds. Because model choice is based on resampling, elastic net is the best of these, with PLS a close second.

(f) Would you replace the permeability lab experiment?

best_perm <- perm_models[[perm_compare$model[which.min(perm_compare$cv_RMSE)]]]
tibble(observed = y_test, predicted = predict(best_perm, x_test)) |>
  ggplot(aes(observed, predicted)) +
  geom_point() +
  geom_abline(linetype = "dashed") +
  labs(title = "Test set: predicted vs observed permeability")

No. Even the best model has an RMSE of about 11 on a response whose median is 4.91, and it explains only about half of the variance. Its predictions are too rough to replace the lab measurement for individual compounds. It could still work as a cheap screen: on the test set every compound it predicted below 10 really was below 10, so it could rank candidates and send only the promising ones to the lab assay.

---

Exercise 6.3 — Chemical manufacturing process

(a) Load the data

data(ChemicalManufacturingProcess)
chem <- ChemicalManufacturingProcess
dim(chem)

## [1] 176  58

table(sub("[0-9]+$", "", names(chem)[-1]))

## 
##   BiologicalMaterial ManufacturingProcess 
##                   12                   45

(b) Impute missing values

miss <- colSums(is.na(chem))
c(columns_with_NA = sum(miss > 0), total_NA = sum(miss))

## columns_with_NA        total_NA 
##              28             106

head(sort(miss[miss > 0], decreasing = TRUE), 5)

## ManufacturingProcess03 ManufacturingProcess11 ManufacturingProcess10 
##                     15                     10                      9 
## ManufacturingProcess25 ManufacturingProcess26 
##                      5                      5

Missing values are filled with k-nearest-neighbors imputation (knnImpute, k = 5) inside preProcess. The imputation is learned from the training set only and then applied to the test set, so no test information leaks into the model.

x_chem <- chem |> select(-Yield)
y_chem <- chem$Yield

set.seed(624)
in_train_c <- createDataPartition(y_chem, p = 0.8, list = FALSE)
xc_train_raw <- x_chem[in_train_c, ];  yc_train <- y_chem[in_train_c]
xc_test_raw  <- x_chem[-in_train_c, ]; yc_test  <- y_chem[-in_train_c]

pp <- preProcess(xc_train_raw, method = c("knnImpute", "nzv"))
xc_train <- predict(pp, xc_train_raw)
xc_test  <- predict(pp, xc_test_raw)
c(NA_after_train = sum(is.na(xc_train)), NA_after_test = sum(is.na(xc_test)),
  predictors = ncol(xc_train))

## NA_after_train  NA_after_test     predictors 
##              0              0             56

pp$method$remove

## [1] "BiologicalMaterial07"

knnImpute also centers and scales every predictor, and the nzv step drops the near-zero-variance column listed above.

(c) Split, pre-process, and tune a model

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

set.seed(624)
chem_pls <- train(xc_train, yc_train, method = "pls",
                  tuneLength = 20, trControl = ctrl_c)

set.seed(624)
chem_enet <- 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)

chem_cv <- bind_rows(
  PLS = getTrainPerf(chem_pls),
  `Elastic net` = getTrainPerf(chem_enet),
  .id = "model"
)
knitr::kable(chem_cv[, c("model", "TrainRMSE", "TrainRsquared", "TrainMAE")], digits = 3)

model	TrainRMSE	TrainRsquared	TrainMAE
PLS	1.187	0.589	0.983
Elastic net	1.167	0.590	0.951

chem_pls$bestTune

##   ncomp
## 4     4

chem_enet$bestTune

##    alpha   lambda
## 53     1 0.078476

plot(chem_pls, main = "PLS: resampled RMSE by components")

Both models are tuned with 10-fold CV repeated 5 times. PLS picks 4 components. Elastic net picks alpha = 1 (the lasso) with lambda = 0.0785. Elastic net has the best resampled performance: RMSE 1.167, \(R^2\) 0.59, slightly ahead of PLS.

(d) Predict the test set

chem_test <- bind_rows(
  PLS = postResample(predict(chem_pls, xc_test), yc_test),
  `Elastic net` = postResample(predict(chem_enet, xc_test), yc_test),
  .id = "model"
)
knitr::kable(chem_test, digits = 3)

model	RMSE	Rsquared	MAE
PLS	1.596	0.500	1.077
Elastic net	1.171	0.729	0.894

On the test set the elastic net reaches RMSE 1.171 and \(R^2\) 0.729. That RMSE matches its resampled RMSE, and the \(R^2\) is higher, so the model generalizes well. PLS does worse on the test set (RMSE 1.596) than in resampling, which is more evidence for choosing the elastic net. With only 32 test runs, these numbers are noisy.

(e) Most important predictors

chem_best <- list(PLS = chem_pls, `Elastic net` = chem_enet)[[
  chem_cv$model[which.min(chem_cv$TrainRMSE)]]]
imp <- varImp(chem_best)$importance |>
  rownames_to_column("predictor") |>
  select(predictor, importance = Overall) |>
  arrange(desc(importance)) |>
  mutate(type = ifelse(grepl("^Biological", predictor), "Biological", "Process"))

top_n_imp <- 15
imp |>
  slice_head(n = top_n_imp) |>
  ggplot(aes(importance, reorder(predictor, importance), fill = type)) +
  geom_col() +
  labs(title = "Top 15 predictors by importance", x = "Scaled importance", y = NULL)

table(head(imp$type, top_n_imp))

## 
## Biological    Process 
##          2         13

Importance for the lasso model is the absolute size of each coefficient on the scaled predictors. The lasso keeps 20 of 56 predictors. Process variables dominate the list: 13 of the top 15 are manufacturing process variables, and the top 2 (ManufacturingProcess32 and 09) are well ahead of everything else. Only BiologicalMaterial06 (4th) and BiologicalMaterial05 (13th) make the top 15.

(f) Relationships between top predictors and yield

top_vars <- head(imp$predictor, 6)
cor_tbl <- sapply(top_vars, function(v) cor(xc_train[[v]], yc_train))
knitr::kable(tibble(predictor = top_vars, cor_with_yield = cor_tbl), digits = 3)

predictor	cor_with_yield
ManufacturingProcess32	0.595
ManufacturingProcess09	0.498
ManufacturingProcess17	-0.349
BiologicalMaterial06	0.491
ManufacturingProcess34	0.245
ManufacturingProcess37	-0.140

xc_train |>
  select(all_of(top_vars)) |>
  mutate(Yield = yc_train) |>
  pivot_longer(-Yield, names_to = "predictor") |>
  mutate(predictor = factor(predictor, levels = top_vars)) |>
  ggplot(aes(value, Yield)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ predictor, scales = "free_x") +
  labs(title = "Yield vs top predictors (training set, scaled)", x = "Scaled value")

ManufacturingProcess32 and 09 have the strongest positive relationships with yield (correlations of about 0.6 and 0.5), and ManufacturingProcess17 and 37 are negatively related. BiologicalMaterial06 is also positively related, but raw material quality is fixed before the run starts, so it mostly helps predict yield rather than control it.

Process variables can be changed, which makes them the useful finding. Running future batches with higher values of ManufacturingProcess32, 09 and 34, and lower values of 17 and 37, within their observed ranges, should raise yield. These are observational correlations, so the next step is a designed experiment on these few settings to confirm the effects are causal before changing the process permanently.