6.2

Developing a model to predict permeability (see Sect. 1.4) could save significant resources for a pharmaceutical company, while at the same time more rapidly identifying molecules that have a sufficient permeability to become a drug:

(a) Start R and use these commands to load the data. The matrix fingerprints contains the 1,107 binary molecular predictors for the 165 compounds, while permeability contains permeability response.

data(permeability)
head(permeability)
##   permeability
## 1       12.520
## 2        1.120
## 3       19.405
## 4        1.730
## 5        1.680
## 6        0.510
set.seed(5)

(b) The fingerprint predictors indicate the presence or absence of substructures of a molecule and are often sparse meaning that relatively few of the molecules contain each substructure. Filter out the predictors that have low frequencies using the nearZeroVar function from the caret package. How many predictors are left for modeling?

filtered_fingerprints <- fingerprints[, -nearZeroVar(fingerprints)]
ncol(filtered_fingerprints)
## [1] 388

We are left with 388 predictors with when filtering out the predictors with near zero variance

(c) Split the data into a training and a test set, pre-process the data, and tune a PLS model. How many latent variables are optimal and what is the corresponding resampled estimate of R2?

train_split_index <- createDataPartition(permeability, p = 0.8, list = FALSE)

trainX <- filtered_fingerprints[train_split_index, ]
testX <- filtered_fingerprints[-train_split_index, ]
trainY <- permeability[train_split_index]
testY <- permeability[-train_split_index]

ctrl <- trainControl(method = "cv", number = 10)
plsTune <- train(trainX, trainY, 
                 method = "pls", 
                 tuneLength = 20, 
                 trControl = ctrl, 
                 preProc = c("center", "scale"))

plot(plsTune)

cat("Optimal latent variables:", plsTune$bestTune$ncomp, " | R squared:", max(plsTune$results$Rsquared), "\n")
## Optimal latent variables: 8  | R squared: 0.488569

Here we see that for the PLS model the RMSE is minimized at the latent variable count of 8. It takes Eight components to get to the lowest RMSE with a 80/20 training test split for this data set and achieves an R squared: 0.488569

(d) Predict the response for the test set. What is the test set estimate of R2?

pred <- predict(plsTune, testX)
R2(pred, testY)
## [1] 0.4755793
plot(pred, testY)
abline(lm(testY ~ pred))

When predicting the response for the test set using the model we get a test set estimate of R squared = 0.4755793 which is very close to the model R squared of 0.488569 that we achieved

(e) Try building other models discussed in this chapter. Do any have better predictive performance?

PCR

pcrFit <- train(trainX, trainY, method = "pcr", tuneLength = 20, trControl = ctrl, preProc = c("center","scale"))
cat("R squared:", max(pcrFit$results$Rsquared), "\n")
## R squared: 0.4849905

Ridge

ridgeFit <- train(trainX, trainY, method = "ridge", tuneLength = 10, trControl = ctrl, preProc = c("center","scale"))
cat("R squared:", max(ridgeFit$results$Rsquared), "\n")
## R squared: 0.5340826

Lasso

lassoFit <- train(trainX, trainY, method = "lasso", tuneLength = 10, trControl = ctrl, preProc = c("center","scale"))
cat("R squared:", max(lassoFit$results$Rsquared), "\n")
## R squared: 0.4323872

(f) Would you recommend any of your models to replace the permeability laboratory experiment?

With this seed and training split Both PCR and Lasso models perform better than the PLS model while ridge performs slightly worse. All three of the models perform similar to the laboratory experiment when comparing R squared values. All these models show that the predictors account for less than half of the variation in the independent variable

6.3.

A chemical manufacturing process for a pharmaceutical product was discussed in Sect. 1.4. In this problem, the objective is to understand the relationship between biological measurements of the raw materials (predictors), 6.5 Computing 139 measurements of the manufacturing process (predictors), and the response of product yield. Biological predictors cannot be changed but can be used to assess the quality of the raw material before processing. On the other hand,manufacturing process predictors can be changed in the manufacturing process. Improving product yield by 1 % will boost revenue by approximately one hundred thousand dollars per batch:

(a) Start R and use these commands to load the data. The matrix processPredictors contains the 57 predictors (12 describing the input biological material and 45 describing the process predictors) for the 176 manufacturing runs. yield contains the percent yield for each run.

PS: the github file shows that the name is ChemicalManufacturingProcess not chemicalManufacturing

data("ChemicalManufacturingProcess")

(b) A small percentage of cells in the predictor set contain missing values. Use an imputation function to fill in these missing values (e.g., see Sect. 3.8).

We can use preprocess to apply K-nearest neighbor imputation

trans <- preProcess(ChemicalManufacturingProcess, method = "knnImpute")
chemicalData <- predict(trans, ChemicalManufacturingProcess)

(c) Split the data into a training and a test set, pre-process the data, and tune a model of your choice from this chapter. What is the optimal value of the performance metric?

train_split_index_2 <- createDataPartition(chemicalData$Yield, p = 0.8, list = FALSE)

train_chem <- chemicalData[train_split_index_2, ]
test_chem <- chemicalData[-train_split_index_2, ]

ctrl_chem <- trainControl(method = "cv", number = 10)

pls_model_chem <- train(
  Yield ~ ., data = train_chem, 
  method = "pls",
  trControl = ctrl_chem,
  tuneLength = 25,
  preProcess = c("center", "scale")
)


cat("Optimal latent variables:", pls_model_chem$bestTune$ncomp, " | R squared:", max(pls_model_chem$results$Rsquared), "\n")
## Optimal latent variables: 4  | R squared: 0.5423709
plot(pls_model_chem)

With Partial least squares we get the highest R squared of 0.566845 at 3 components

(d) Predict the response for the test set. What is the value of the performance metric and how does this compare with the resampled performance metric on the training set?

pred_chem <- predict(pls_model_chem, test_chem)
R2(pred_chem, test_chem$Yield)
## [1] 0.2037987

R squared looks to be hugely changed when testing the model on the test set. We can also try other models

pcr_model_chem <- train(
  Yield ~ ., data = train_chem, method = "pcr",
  trControl = ctrl_chem, tuneLength = 25
)

ridge_model_chem <- train(
  Yield ~ ., data = train_chem, method = "ridge",
  trControl = ctrl_chem, tuneLength = 25,
  preProcess = c("center", "scale")
)

lasso_model_chem <- train(
  Yield ~ ., data = train_chem, method = "lasso",
  trControl = ctrl_chem, tuneLength = 25,
  preProcess = c("center", "scale")
)


train_r2_pcr   <- max(pcr_model_chem$results$Rsquared)
train_r2_ridge <- max(ridge_model_chem$results$Rsquared)
train_r2_lasso <- max(lasso_model_chem$results$Rsquared)

cat("Training (CV) Rsquared values:\n")
## Training (CV) Rsquared values:
cat("  PCR:",   round(train_r2_pcr, 4), "\n")
##   PCR: 0.593
cat("  Ridge:", round(train_r2_ridge, 4), "\n")
##   Ridge: 0.5535
cat("  Lasso:", round(train_r2_lasso, 4), "\n\n")
##   Lasso: 0.6384

Lets see the prediction metrics on the test set

pred_pcr_chem   <- predict(pcr_model_chem, test_chem)
pred_ridge_chem <- predict(ridge_model_chem, test_chem)
pred_lasso_chem <- predict(lasso_model_chem, test_chem)

r2_pcr_chem   <- R2(pred_pcr_chem,   test_chem$Yield)
r2_ridge_chem <- R2(pred_ridge_chem, test_chem$Yield)
r2_lasso_chem <- R2(pred_lasso_chem, test_chem$Yield)

cat("Test Set Rsquared values:\n")
## Test Set Rsquared values:
cat("  PCR:",   round(r2_pcr_chem, 4), "\n")
##   PCR: 0.671
cat("  Ridge:", round(r2_ridge_chem, 4), "\n")
##   Ridge: 0.0243
cat("  Lasso:", round(r2_lasso_chem, 4), "\n")
##   Lasso: 0.3499

(e) Which predictors are most important in the model you have trained? Do either the biological or process predictors dominate the list?

We can print out the predictor Variable Importance Scores for the model to see how each predictor contributes to the model. Here we can see that in the top 10 there seems to be a higher split of Process predictors compared to Biological predictors ( 6 to 4 in favor for Manufacturing ), with the highest contribution coming from Manufacturing process 32

lasso_imp <- varImp(lasso_model_chem)
plot(lasso_imp, top = 20)

(f) Explore the relationships between each of the top predictors and the response. How could this information be helpful in improving yield in future runs of the manufacturing process?

lasso_imp_df <- lasso_imp$importance %>%
  arrange(desc(Overall))

top_vars <- rownames(head(lasso_imp_df,5))
top_vars
## [1] "ManufacturingProcess32" "ManufacturingProcess13" "BiologicalMaterial06"  
## [4] "BiologicalMaterial03"   "ManufacturingProcess09"
for (v in top_vars) {
  p <- ggplot(chemicalData, aes_string(x = v, y = "Yield")) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    theme(plot.title = element_text(size = 12, face = "bold"))
  print(p)
}

Here we can see how the individual predictors from our VIP list affects the yield. We could use these relationship plots to tweak the exposure of these top predictors to the manufacturing process and raw material selection to improve the yield