Exercise 6.2: Predicting permeability

(a) Load and examine the permeability data

Question: What are the predictors and the response?

data("permeability", package = "AppliedPredictiveModeling")

x_perm <- as.data.frame(fingerprints)
names(x_perm) <- paste0("Fingerprint", seq_len(ncol(x_perm)))
y_perm <- as.numeric(permeability)

dim(x_perm)
## [1]  165 1107
length(y_perm)
## [1] 165
head(x_perm[, 1:6])
##   Fingerprint1 Fingerprint2 Fingerprint3 Fingerprint4 Fingerprint5 Fingerprint6
## 1            0            0            0            0            0            1
## 2            0            0            0            0            0            0
## 3            0            0            0            0            0            1
## 4            0            0            0            0            0            0
## 5            0            0            0            0            0            0
## 6            0            0            0            0            0            0

Answer: The dataset contains 165 compounds and 1107 binary fingerprint predictors. A 1 indicates that a particular substructure is present; a 0 indicates it is absent. The response is the measured permeability of each compound.

(b) Remove near-zero-variance predictors

Question: How many predictors remain after filtering out the rare fingerprints?

rare <- nearZeroVar(x_perm)

data.frame(
  Before = ncol(x_perm),
  Removed = length(rare),
  Remaining = ncol(x_perm) - length(rare)
)
##   Before Removed Remaining
## 1   1107     719       388

Answer: Using the default settings in nearZeroVar(), I removed 719 predictors, leaving 388. This is the count for the full dataset. In the modeling steps, I have train() repeat the filter using only the fitting data so that the validation and test observations do not influence preprocessing.

(c) Split the data and tune a PLS model

Question: How many latent variables are optimal, and what is the corresponding cross-validation R²?

First, I split the observations into 80% training and 20% test. I save the training folds so the other models will be compared on the same splits.

set.seed(6242)
train_rows_perm <- sort(sample.int(nrow(x_perm), floor(0.8 * nrow(x_perm))))

x_train_perm <- x_perm[train_rows_perm, ]
x_test_perm <- x_perm[-train_rows_perm, ]
y_train_perm <- y_perm[train_rows_perm]
y_test_perm <- y_perm[-train_rows_perm]

set.seed(16242)
folds_perm <- createFolds(y_train_perm, k = 10, returnTrain = TRUE)
cv_perm <- trainControl(method = "cv", number = 10, index = folds_perm)

c(Training = nrow(x_train_perm), Test = nrow(x_test_perm))
## Training     Test 
##      132       33

Next, I use caret::train() to fit partial least squares (PLS). The model tries one through 20 components. Preprocessing removes sparse fingerprints, centers each retained predictor, and scales it. I choose the number of components with the lowest cross-validation RMSE.

prep_perm <- c("nzv", "center", "scale")

set.seed(6242)
pls_fit <- train(
  x = x_train_perm, y = y_train_perm,
  method = "pls",
  preProcess = prep_perm,
  tuneGrid = data.frame(ncomp = 1:20),
  trControl = cv_perm,
  metric = "RMSE"
)

pls_best <- pls_fit$results[
  pls_fit$results$ncomp == pls_fit$bestTune$ncomp,
  c("ncomp", "RMSE", "Rsquared")
]
knitr::kable(pls_best, digits = 3, row.names = FALSE)
ncomp RMSE Rsquared
2 10.63 0.57 
plot(pls_fit, main = "Permeability: tuning PLS")

Answer: The best model uses 2 latent variables. Its mean cross-validation R² is 0.57, with an RMSE of 10.63. The R² reported by caret is based on squared correlation between the observed and predicted values.

(d) Test the PLS model

Question: What is the test-set R²?

pls_pred <- predict(pls_fit, newdata = x_test_perm)
pls_test <- postResample(pred = pls_pred, obs = y_test_perm)
round(pls_test, 3)
##     RMSE Rsquared      MAE 
##   13.907    0.350    9.470
limits <- range(y_test_perm, pls_pred)
plot(y_test_perm, pls_pred, pch = 19, col = "steelblue",
     xlim = limits, ylim = limits,
     xlab = "Observed permeability", ylab = "Predicted permeability",
     main = "PLS predictions on the test set")
abline(0, 1, lty = 2)

Answer: The test-set R² is 0.35, and test RMSE is 13.907. The dashed line shows where the points would fall if the predictions exactly matched the observed values.

(e) Compare PLS with other models from the chapter

Question: Do any other models predict permeability better?

I compare PLS with principal component regression (PCR) and ridge regression using the same cross-validation folds. PCR chooses components from the predictors without looking at permeability. Ridge keeps the predictors but shrinks their coefficients. In glmnet, alpha = 0 specifies ridge regression; lambda controls the amount of shrinkage.

set.seed(6242)
pcr_fit <- train(
  x = x_train_perm, y = y_train_perm,
  method = "pcr",
  preProcess = prep_perm,
  tuneGrid = data.frame(ncomp = 1:20),
  trControl = cv_perm,
  metric = "RMSE"
)
pcr_fit$bestTune
##    ncomp
## 10    10
ridge_grid <- expand.grid(alpha = 0, lambda = 10^seq(-3, 2, length.out = 26))

set.seed(6242)
ridge_perm <- train(
  x = x_train_perm, y = y_train_perm,
  method = "glmnet",
  preProcess = prep_perm,
  tuneGrid = ridge_grid,
  trControl = cv_perm,
  metric = "RMSE",
  standardize = FALSE,
  lambda = sort(ridge_grid$lambda, decreasing = TRUE)
)
ridge_perm$bestTune
##    alpha lambda
## 21     0     10
pcr_pred <- predict(pcr_fit, newdata = x_test_perm)
ridge_pred <- predict(ridge_perm, newdata = x_test_perm)

pcr_test <- postResample(pred = pcr_pred, obs = y_test_perm)
ridge_test <- postResample(pred = ridge_pred, obs = y_test_perm)

comparison <- data.frame(
  Model = c("PLS", "PCR", "Ridge"),
  CV_RMSE = c(min(pls_fit$results$RMSE), min(pcr_fit$results$RMSE),
              min(ridge_perm$results$RMSE)),
  Test_RMSE = c(pls_test[["RMSE"]], pcr_test[["RMSE"]], ridge_test[["RMSE"]]),
  Test_R2 = c(pls_test[["Rsquared"]], pcr_test[["Rsquared"]],
              ridge_test[["Rsquared"]])
)
knitr::kable(comparison, digits = 3, row.names = FALSE)
Model CV_RMSE Test_RMSE Test_R2
PLS 10.630 13.907 0.350
PCR 10.815 13.949 0.350
Ridge 10.150 13.541 0.381

Answer: Ridge had the lowest cross-validation RMSE. On this test split, Ridge had the lowest RMSE and Ridge had the highest R². These comparisons are based on one split, so they do not guarantee that the same ranking will hold on new data.

(f) Decide whether to replace the laboratory experiment

Question: Would I recommend replacing the laboratory experiment with any of these models?

Answer: No. I would use the model to help decide which compounds to test first. The dataset is small, and a model could make a wrong decision about a compound near the permeability cutoff. Before replacing any lab measurements, I would check the model on new compounds and see how often it incorrectly classifies compounds near that cutoff.

Exercise 6.3: Predicting manufacturing yield

(a) Load and examine the manufacturing data

Question: What are the predictors and the response?

data("ChemicalManufacturingProcess", package = "AppliedPredictiveModeling")

manufacturing <- ChemicalManufacturingProcess
x_mfg <- manufacturing[, names(manufacturing) != "Yield"]
y_mfg <- manufacturing$Yield

data.frame(
  Runs = nrow(x_mfg),
  Biological = sum(grepl("^BiologicalMaterial", names(x_mfg))),
  Manufacturing = sum(grepl("^ManufacturingProcess", names(x_mfg))),
  Total_predictors = ncol(x_mfg)
)
##   Runs Biological Manufacturing Total_predictors
## 1  176         12            45               57

Answer: The response is Yield, measured as a percent. The dataset has 176 runs and 57 predictors, including 12 biological measurements and 45 manufacturing process measurements.

(b) Impute missing predictor values

Question: How can I fill in the missing predictor values?

I first count the missing values. Then I split the data before imputing so the test observations do not affect the medians. I use preProcess() here to show what median imputation does. In part (c), train() repeats the imputation inside each cross-validation fold.

set.seed(6243)
train_rows_mfg <- sort(sample.int(nrow(x_mfg), floor(0.8 * nrow(x_mfg))))

x_train_mfg <- x_mfg[train_rows_mfg, ]
x_test_mfg <- x_mfg[-train_rows_mfg, ]
y_train_mfg <- y_mfg[train_rows_mfg]
y_test_mfg <- y_mfg[-train_rows_mfg]

imputer <- preProcess(x_train_mfg, method = "medianImpute")
train_imputed <- predict(imputer, newdata = x_train_mfg)
test_imputed <- predict(imputer, newdata = x_test_mfg)

data.frame(
  Set = c("Training", "Test"),
  Missing_before = c(sum(is.na(x_train_mfg)), sum(is.na(x_test_mfg))),
  Missing_after = c(sum(is.na(train_imputed)), sum(is.na(test_imputed)))
)
##        Set Missing_before Missing_after
## 1 Training             81             0
## 2     Test             25             0

Answer: There were 106 missing predictor values. I replaced each missing value with the median of its predictor, calculated from the training data. There are no missing values in the imputed predictor sets. I did not impute the response.

(c) Tune a regression model

Question: What is the best cross-validation performance for the model I chose?

I use ridge regression because some predictors may be related to each other, and ridge shrinks coefficients without selecting a subset. I remove any predictors that are constant in the fitting data because they provide no information for that fit. I evaluate each candidate lambda using 10-fold cross-validation and select the one with the lowest RMSE.

set.seed(16243)
cv_mfg <- trainControl(method = "cv", number = 10)

ridge_mfg <- train(
  x = x_train_mfg, y = y_train_mfg,
  method = "glmnet",
  preProcess = c("zv", "medianImpute", "center", "scale"),
  tuneGrid = ridge_grid,
  trControl = cv_mfg,
  metric = "RMSE",
  standardize = FALSE,
  lambda = sort(ridge_grid$lambda, decreasing = TRUE)
)

mfg_best <- ridge_mfg$results[
  ridge_mfg$results$lambda == ridge_mfg$bestTune$lambda,
  c("lambda", "RMSE", "Rsquared")
]
knitr::kable(mfg_best, digits = 4, row.names = FALSE)
lambda RMSE Rsquared
3.9811 1.4349 0.4823 
plot(ridge_mfg, main = "Manufacturing yield: tuning ridge")

Answer: The selected lambda is 3.981. The lowest mean cross-validation RMSE is 1.435 percentage points. The corresponding cross-validation R² is 0.482.

(d) Test the manufacturing model

Question: What is the test RMSE, and how does it compare with cross-validation RMSE?

mfg_pred <- predict(ridge_mfg, newdata = x_test_mfg)
mfg_test <- postResample(pred = mfg_pred, obs = y_test_mfg)
rmse_difference <- mfg_test[["RMSE"]] - mfg_best$RMSE

knitr::kable(data.frame(
  Set = c("Cross-validation", "Test"),
  RMSE = c(mfg_best$RMSE, mfg_test[["RMSE"]]),
  R2 = c(mfg_best$Rsquared, mfg_test[["Rsquared"]])
), digits = 3)
Set RMSE R2
Cross-validation 1.435 0.482
Test 1.352 0.655

Answer: The test RMSE is 1.352 percentage points, compared with 1.435 in cross-validation. The test RMSE was 0.083 percentage points lower before rounding the two RMSE values. The test set has 36 runs, so I would not read too much into the difference from this one split.

(e) Identify the most important predictors

Question: Which predictors are most important? Do biological or manufacturing process predictors dominate?

I use varImp() on the ridge model. For this model, predictors are ranked by the absolute sizes of their coefficients at the selected lambda. Centering and scaling the predictors before fitting makes the coefficient sizes more comparable.

importance <- varImp(ridge_mfg, scale = FALSE)$importance
importance <- importance[order(importance$Overall, decreasing = TRUE), ,
                         drop = FALSE]

top_names <- head(rownames(importance), 6)
top_predictors <- data.frame(
  Predictor = top_names,
  Type = ifelse(grepl("^BiologicalMaterial", top_names),
                "Biological", "Manufacturing process"),
  Importance = importance[top_names, "Overall"]
)
process_count <- sum(top_predictors$Type == "Manufacturing process")
knitr::kable(top_predictors, digits = 3, row.names = FALSE)
Predictor Type Importance
ManufacturingProcess32 Manufacturing process 0.184
ManufacturingProcess09 Manufacturing process 0.137
ManufacturingProcess13 Manufacturing process 0.137
ManufacturingProcess36 Manufacturing process 0.134
ManufacturingProcess17 Manufacturing process 0.128
ManufacturingProcess33 Manufacturing process 0.093

Answer: The highest-ranked predictor is ManufacturingProcess32. Of the top six, 6 are manufacturing process predictors and 0 are biological predictors. Manufacturing process predictors dominate this list. This ranking measures importance to the fitted model; it does not show which variables cause yield to change.

(f) Explore the top predictors and relate them to future runs

Question: How does yield relate to each top predictor, and how could the relationships help improve yield?

I plot the observed training values for each top predictor against yield. I leave out missing pairs rather than plotting imputed values as if they were actual measurements. The orange line shows the trend in the observed points.

correlations <- numeric(length(top_names))
old_par <- par(mfrow = c(2, 3), mar = c(4, 4, 2.8, 0.7))

for (i in seq_along(top_names)) {
  variable <- top_names[i]
  observed <- !is.na(x_train_mfg[[variable]])
  x <- x_train_mfg[[variable]][observed]
  y <- y_train_mfg[observed]

  plot(x, y, pch = 19, cex = 0.6, col = "steelblue",
       xlab = "Predictor value", ylab = "Yield (%)",
       main = variable, cex.main = 0.8)
  lines(lowess(x, y), col = "darkorange3", lwd = 2)
  correlations[i] <- cor(x, y)
}

par(old_par)

knitr::kable(data.frame(
  Predictor = top_names,
  Correlation_with_yield = correlations,
  Direction = ifelse(correlations > 0, "Positive",
                     ifelse(correlations < 0, "Negative", "Zero"))
), digits = 3, row.names = FALSE)
Predictor Correlation_with_yield Direction
ManufacturingProcess32 0.617 Positive
ManufacturingProcess09 0.483 Positive
ManufacturingProcess13 -0.477 Negative
ManufacturingProcess36 -0.547 Negative
ManufacturingProcess17 -0.371 Negative
ManufacturingProcess33 0.454 Positive

Answer: The table shows the direction for each predictor. The strongest observed relationship by absolute correlation was for ManufacturingProcess32 (correlation 0.617). A positive relationship means higher predictor values tended to occur with higher yield; a negative one means they tended to occur with lower yield.

I would use the process variables to decide which settings to investigate in controlled trial runs. For biological variables, I would use any useful relationships to screen or monitor incoming material, since these measurements cannot be directly changed during production. The patterns are associations, so I would not change operating settings based only on these plots.

Reference

Kuhn, M., & Johnson, K. (2013). Applied Predictive Modeling. Springer. Chapter 6, exercises 6.2 and 6.3.