Exercise 6.2: Permeability

(a) Load the data

data(permeability)
dim(fingerprints)          # compounds x binary predictors
## [1]  165 1107
length(permeability)       # response, one value per compound
## [1] 165
perm <- as.numeric(permeability)

(b) Filter low-frequency predictors

nzv_idx <- nearZeroVar(fingerprints)
fp      <- fingerprints[, -nzv_idx, drop = FALSE]
n_removed <- length(nzv_idx)
n_left    <- ncol(fp)

Answer (b). The default filter removes 719 predictors and leaves 388 across all compounds. This count describes the full dataset. For the models, I apply the filter separately to the fitting rows in each fold.

(c) Split, pre-process, tune PLS

I split the original fingerprints into about 80% training and 20% test rows. All three models use the same 10 cross-validation folds and select tuning settings by the lowest mean RMSE.

set.seed(624)
train_idx <- createDataPartition(perm, p = 0.8, list = FALSE)

x_train <- fingerprints[train_idx, , drop = FALSE]
y_train <- perm[train_idx]
x_test  <- fingerprints[-train_idx, , drop = FALSE]
y_test  <- perm[-train_idx]

folds <- createFolds(y_train, k = 10, returnTrain = TRUE)
ctrl  <- trainControl(method = "cv", index = folds)
fp_preprocess <- c("nzv", "center", "scale")

PLS combines fingerprints into latent components using their relationship with permeability. I test 1 through 20 components. By putting preProcess inside train(), each validation fold receives the fingerprint filter, means, and standard deviations learned from its fitting rows.

set.seed(624)
pls_fit <- train(x_train, y_train, method = "pls",
                 tuneLength = 20, trControl = ctrl,
                 preProcess = fp_preprocess, metric = "RMSE")

plot(pls_fit)

pls_best <- pls_fit$results[pls_fit$results$ncomp == pls_fit$bestTune$ncomp, ]
pls_best
##    ncomp     RMSE  Rsquared      MAE  RMSESD RsquaredSD    MAESD
## 11    11 11.14376 0.5423326 8.571912 3.30798   0.220232 2.592576

Answer (c). PLS selects 11 latent components by mean cross-validation (CV) RMSE. That setting gives a CV RMSE of 11.144 and a mean CV \(R^2\) of 0.542. The final PLS model retains 407 fingerprints when the filter is fitted on all 133 training compounds.

(d) Test-set performance

predict() applies the final model’s saved preprocessing to the raw test fingerprints. The \(R^2\) reported by postResample() is the squared correlation between observed and predicted values.

pls_test <- postResample(predict(pls_fit, x_test), y_test)
pls_test
##       RMSE   Rsquared        MAE 
## 11.7321881  0.4541521  9.2002902

Answer (d). On 32 test compounds, PLS has a test \(R^2\) of 0.454 and a test RMSE of 11.732 permeability units.

(e) Other models from the chapter

Ridge regression and elastic net shrink the fingerprint coefficients. I tune both on the same folds as PLS and evaluate the selected settings on the same test compounds.

For a faster knit, install doParallel once in the RStudio Console with install.packages("doParallel"). The chunk below uses up to four local workers to run the cross-validation fits concurrently.

fit_penalized_models <- function() {
  cores <- parallel::detectCores(logical = FALSE)
  backend <- foreach::getDoParName()
  if (requireNamespace("doParallel", quietly = TRUE) &&
      (is.null(backend) || backend == "doSEQ") && !is.na(cores) && cores > 1L) {
    workers <- tryCatch(
      parallel::makePSOCKcluster(min(4L, cores), setup_timeout = 15),
      error = function(e) NULL
    )
    if (!is.null(workers)) {
      on.exit({
        foreach::registerDoSEQ()
        parallel::stopCluster(workers)
      })
      parallel::clusterCall(workers, .libPaths, .libPaths())
      doParallel::registerDoParallel(workers)
    }
  }

  set.seed(624)
  ridge <- train(x_train, y_train, method = "ridge",
                 tuneGrid = data.frame(lambda = seq(0.001, 0.2, length = 15)),
                 trControl = ctrl, preProcess = fp_preprocess, metric = "RMSE")

  set.seed(624)
  enet <- train(x_train, y_train, method = "enet",
                tuneGrid = expand.grid(lambda   = c(0.001, 0.01, 0.1),
                                       fraction = seq(0.05, 1, length = 20)),
                trControl = ctrl, preProcess = fp_preprocess, metric = "RMSE")

  list(Ridge = ridge, `Elastic net` = enet)
}

penalized <- fit_penalized_models()
res <- resamples(c(list(PLS = pls_fit), penalized))
cv_stats <- summary(res, metric = c("RMSE", "Rsquared"))$statistics
test_tab <- rbind(
  PLS   = pls_test,
  Ridge = postResample(predict(penalized$Ridge, x_test), y_test),
  `Elastic net` = postResample(predict(penalized[["Elastic net"]], x_test), y_test)
)

cv_rmse <- cv_stats$RMSE[rownames(test_tab), "Mean"]
cv_r2 <- cv_stats$Rsquared[rownames(test_tab), "Mean"]
cv_winner <- names(which.min(cv_rmse))
test_winner <- rownames(test_tab)[which.min(test_tab[, "RMSE"])]
cv_winner_range <- range(res$values[[paste0(cv_winner, "~RMSE")]], na.rm = TRUE)

perm_scores <- data.frame(
  Model = rownames(test_tab),
  CV_RMSE = cv_rmse,
  CV_R2 = cv_r2,
  Test_RMSE = test_tab[, "RMSE"],
  Test_R2 = test_tab[, "Rsquared"]
)
knitr::kable(perm_scores, digits = 3, row.names = FALSE,
             col.names = c("Model", "CV RMSE", "CV R²", "Test RMSE", "Test R²"))
Model CV RMSE CV R² Test RMSE Test R²
PLS 11.144 0.542 11.732 0.454
Ridge 11.487 0.554 11.560 0.490
Elastic net 10.813 0.572 10.509 0.528

Finding (e). Elastic net has the lowest mean CV RMSE (10.813). Elastic net has the lowest test RMSE (10.509). The CV winner’s RMSE varies from 5.03 to 14.93 across folds, so I treat the model differences as estimates from this split.

(f) Replace the laboratory experiment?

summary(y_train)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.060   1.550   4.495  12.181  15.470  55.600

Finding for Exercise 6.2. I would keep the laboratory assay and use the CV-selected model to prioritize compounds for testing. Its test RMSE is 10.51 permeability units, compared with a training interquartile range of 13.92 units. That error matters when a prediction is near the screening cutoff. Replacing the assay would require an acceptable error threshold and independent validation on new compounds.

Exercise 6.3: Chemical manufacturing process

(a) Load the data

data(ChemicalManufacturingProcess)
x_all <- ChemicalManufacturingProcess[, -1]   # 57 predictors
y_all <- ChemicalManufacturingProcess$Yield

dim(x_all)
## [1] 176  57
table(grepl("^Biological", names(x_all)))     # 12 biological, 45 process
## 
## FALSE  TRUE 
##    45    12

(b) Impute missing values

miss <- colSums(is.na(x_all))
miss[miss > 0]
## ManufacturingProcess01 ManufacturingProcess02 ManufacturingProcess03 
##                      1                      3                     15 
## ManufacturingProcess04 ManufacturingProcess05 ManufacturingProcess06 
##                      1                      1                      2 
## ManufacturingProcess07 ManufacturingProcess08 ManufacturingProcess10 
##                      1                      1                      9 
## ManufacturingProcess11 ManufacturingProcess12 ManufacturingProcess14 
##                     10                      1                      1 
## ManufacturingProcess22 ManufacturingProcess23 ManufacturingProcess24 
##                      1                      1                      1 
## ManufacturingProcess25 ManufacturingProcess26 ManufacturingProcess27 
##                      5                      5                      5 
## ManufacturingProcess28 ManufacturingProcess29 ManufacturingProcess30 
##                      5                      5                      5 
## ManufacturingProcess31 ManufacturingProcess33 ManufacturingProcess34 
##                      5                      5                      5 
## ManufacturingProcess35 ManufacturingProcess36 ManufacturingProcess40 
##                      5                      5                      1 
## ManufacturingProcess41 
##                      1
missing_fraction <- mean(is.na(x_all))
missing_fraction
## [1] 0.01056619

Answer (b). The predictors contain 106 missing cells, about 1.06% of all predictor cells. I use k-nearest-neighbor imputation with \(k = 5\). knnImpute centers and scales the predictors for distance calculations, then uses nearby fitting rows to fill each missing measurement.

chem_preprocess <- c("nzv", "knnImpute")

Both models fit their preprocessing inside cross-validation. This keeps the validation rows out of the scaling and imputation for their fold.

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

I tune two models from the chapter (PLS and elastic net) with identical folds and keep the one with the lower resampled RMSE.

set.seed(624)
tr <- createDataPartition(y_all, p = 0.8, list = FALSE)
x_train <- x_all[tr, , drop = FALSE]
y_train <- y_all[tr]
x_test  <- x_all[-tr, , drop = FALSE]
y_test  <- y_all[-tr]

folds <- createFolds(y_train, k = 10, returnTrain = TRUE)
ctrl  <- trainControl(method = "cv", index = folds)

set.seed(624)
pls_chem <- train(x_train, y_train, method = "pls",
                  tuneLength = 20, trControl = ctrl,
                  preProcess = chem_preprocess, metric = "RMSE")

set.seed(624)
enet_chem <- train(x_train, y_train, method = "enet",
                   tuneGrid = expand.grid(lambda   = c(0.001, 0.01, 0.1),
                                          fraction = seq(0.05, 1, length = 20)),
                   trControl = ctrl, preProcess = chem_preprocess, metric = "RMSE")
chem_models <- list(PLS = pls_chem, `Elastic net` = enet_chem)
chem_cv <- do.call(rbind, lapply(chem_models, function(model) {
  getTrainPerf(model)[c("TrainRMSE", "TrainRsquared")]
}))
best_name <- rownames(chem_cv)[which.min(chem_cv$TrainRMSE)]
best_fit  <- chem_models[[best_name]]
best_res  <- chem_cv[best_name, ]

knitr::kable(chem_cv, digits = 3, col.names = c("CV RMSE", "CV R²"))
CV RMSE CV R²
PLS 1.395 0.550
Elastic net 1.201 0.579 
best_fit$bestTune
##   fraction lambda
## 2      0.1  0.001

Answer (c). With 144 training runs and 32 held-out runs, Elastic net has the lowest mean CV RMSE, 1.201 percentage points of yield. Its mean CV \(R^2\) is 0.579.

(d) Test-set performance

test_processed <- predict(best_fit$preProcess, x_test)
c(Missing_before = sum(is.na(x_test)), Missing_after = sum(is.na(test_processed)))
## Missing_before  Missing_after 
##              7              0
test_perf <- postResample(predict(best_fit, x_test), y_test)
rmse_gap <- test_perf[["RMSE"]] - best_res$TrainRMSE
test_perf
##      RMSE  Rsquared       MAE 
## 1.2808703 0.7213274 0.9988901

Finding (d). The test RMSE is 1.281 percentage points. That is 0.08 points higher than the CV RMSE. Test \(R^2\) is 0.721, compared with 0.579 in CV. RMSE measures errors in the magnitude of the predictions, while \(R^2\) measures how closely predictions and observations move together. I interpret the gap in the context of one held-out split.

Some runs share a starting-material batch. The evaluation splits by run, so performance on entirely new batches would require batch IDs and a split that holds out whole batches. Dataset documentation

(e) Most important predictors

I rank predictors using the selected model. For elastic net, the score is the absolute fitted coefficient at the selected tuning values. The preprocessing puts the predictors on comparable scales. For PLS, the score comes from the model-specific varImp() method.

if (best_name == "Elastic net") {
  coefficients <- elasticnet::predict.enet(
    best_fit$finalModel,
    s = best_fit$bestTune$fraction,
    mode = "fraction",
    type = "coefficients"
  )$coefficients
  nonzero <- coefficients[coefficients != 0]
  scores <- abs(nonzero)
} else {
  pls_importance <- varImp(best_fit)$importance
  scores <- setNames(pls_importance$Overall, rownames(pls_importance))
}

scores <- sort(scores, decreasing = TRUE)
vi <- data.frame(Predictor = names(scores), Score = unname(scores))
vi$Type <- ifelse(startsWith(vi$Predictor, "Biological"), "Biological", "Process")

top_vi <- head(vi, 10)
knitr::kable(top_vi, digits = 3, row.names = FALSE)
Predictor Score Type
ManufacturingProcess32 0.760 Process
ManufacturingProcess09 0.520 Process
ManufacturingProcess34 0.098 Process
ManufacturingProcess17 0.089 Process
ManufacturingProcess06 0.085 Process
BiologicalMaterial06 0.055 Biological
ManufacturingProcess44 0.050 Process
ManufacturingProcess36 0.043 Process
ManufacturingProcess37 0.040 Process
ManufacturingProcess11 0.033 Process

Finding (e). ManufacturingProcess32 ranks first. Of the top 10 predictors, 9 come from the process and 1 from the biological starting material. The full predictor pool contains 45 process measurements and 12 biological measurements, so there were more process measurements available to rank.

(f) Relationships between top predictors and yield

I plot the original measurements for the six highest-ranked predictors. The points are observed training runs. Each correlation uses the runs with that predictor measured, and the straight trend line summarizes the direction across its observed range.

top_names <- head(vi$Predictor, 6)
top_raw <- x_train[, top_names, drop = FALSE]

featurePlot(x = top_raw, y = y_train,
            plot = "scatter", type = c("p", "r"),
            labels = c("Recorded measurement", "Yield (%)"),
            layout = c(3, 2), as.table = TRUE, between = list(x = 1, y = 1))
Original training measurements and yield for the six highest-ranked predictors.

Original training measurements and yield for the six highest-ranked predictors. 

correlations <- drop(cor(top_raw, y_train, use = "pairwise.complete.obs"))
relationships <- data.frame(
  Predictor = top_names,
  Observed_runs = colSums(!is.na(top_raw)),
  Pearson_r = correlations,
  Direction = ifelse(correlations > 0, "Positive",
                     ifelse(correlations < 0, "Negative", "Zero"))
)
knitr::kable(relationships, digits = 3, row.names = FALSE,
             col.names = c("Predictor", "Observed runs", "Pearson r", "Direction"))
Predictor Observed runs Pearson r Direction
ManufacturingProcess32 144 0.595 Positive
ManufacturingProcess09 144 0.498 Positive
ManufacturingProcess34 139 0.241 Positive
ManufacturingProcess17 144 -0.349 Negative
ManufacturingProcess06 142 0.393 Positive
BiologicalMaterial06 144 0.491 Positive

Finding for Exercise 6.3. On this split, ManufacturingProcess32 and ManufacturingProcess09 have the largest fitted coefficients in absolute value and show positive relationships with yield. ManufacturingProcess17 trends downward, and ManufacturingProcess34 has a narrower range of observed values, with most runs near 2.5.

I would first check what ManufacturingProcess32 and ManufacturingProcess09 measure and whether either is adjustable. A controlled trial within the observed operating range, accounting for starting-material batch, could test whether a change raises yield. If a measurement is observed during production, it could serve as a monitoring signal. Lower BiologicalMaterial06 readings tend to accompany lower yield, so this measurement could flag incoming material for closer inspection after validation on new batches. Dataset documentation