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:

  1. Start R and use these commands to load the data:
library(AppliedPredictiveModeling)
data(permeability)

The matrix fingerprints contains the 1,107 binary molecular predictors for the 165 compounds, while permeability contains permeability response.

  1. 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?
  It would appear that the data was narrowed down from 1107 to 388 columns.
library(caret)
library(dplyr)

modeling <- fingerprints[,-nearZeroVar(fingerprints)]


as.data.frame(modeling)
model_df <- as.data.frame(modeling) %>% bind_cols(permeability)
  1. 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 re sampled estimate of R2?
set.seed(123)
train_index <- createDataPartition(model_df$permeability , p=.8, list=F)

train <-  model_df[ train_index,] 
       
  
test <- model_df[-train_index,]
pls_model <- train(
  permeability ~ ., data = train, method = "pls",
  center = TRUE,
  trControl = trainControl("cv", number = 10),
  tuneLength = 25
)
plot(pls_model)

  The the best tune occurred with 7 components for our model, with an R-squared of 0.5373659
pls_model$results %>%  inner_join(  pls_model$bestTune)
  1. Predict the response for the test set. What is the test set estimate of R2?
  We see blow that our R-squared value is calculated to be 0.349713.
  postResample(pred = predict(pls_model,test), obs =test$permeability)
      RMSE   Rsquared        MAE 
11.9762238  0.3497131  8.4615399
  1. Try building other models discussed in this chapter. Do any have better predictive performance?
ridgeGrid <- data.frame(.lambda = seq(0, .1, length = 15))

set.seed(100)

 ridgeRegFit <- train(
   permeability ~ ., data = train,
   method = "ridge",
   center = TRUE,
   tuneGrid = ridgeGrid,
   trControl = trainControl("cv", number = 10))


 enetGrid <- expand.grid(.lambda = c(0, 0.01, .1),
.fraction = seq(.05, 1, length = 20)) 
 
enetTune <- train(
  
  permeability ~ ., data = train,
   method = "enet",
  center = TRUE,
   tuneGrid = enetGrid,
   trControl = trainControl("cv", number = 10)
 )
pcr_model <- train(
  permeability ~ ., data = train, method = "pcr",
  center = TRUE,
  trControl = trainControl("cv", number = 10),
  tuneLength = 25
)

glmnet_model <-  train(
  permeability~ ., data = train,
  center = TRUE,
  method = "glmnet",
  trControl = trainControl("cv", number = 10)
)
  1. Would you recommend any of your models to replace the permeability laboratory experiment?
It would appear that many of these models model that data very well. It would require more testing to find the best but, but likely any one of these choices would work very well.
postResample(pred = predict(ridgeRegFit,test), obs =test$permeability)
      RMSE   Rsquared        MAE 
12.8430320  0.4004427  9.2116952 
postResample(pred = predict(enetTune,test), obs =test$permeability)
      RMSE   Rsquared        MAE 
11.5871592  0.3519218  7.2484765 
postResample(pred = predict(pcr_model,test), obs =test$permeability)
      RMSE   Rsquared        MAE 
11.7605501  0.3336618  7.5132196 
postResample(pred = predict(glmnet_model,test), obs =test$permeability)
      RMSE   Rsquared        MAE 
11.0800848  0.3581228  7.3179572 
postResample(pred = predict(pls_model,test), obs =test$permeability)
      RMSE   Rsquared        MAE 
11.9762238  0.3497131  8.4615399

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), 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:

  1. Start R and use these commands to load the data:
library(AppliedPredictiveModeling)
data("ChemicalManufacturingProcess")

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.

ChemicalManufacturingProcess
  1. 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).
library(RANN)
impute <- preProcess(ChemicalManufacturingProcess, "knnImpute")

 chem_data <- predict(impute, ChemicalManufacturingProcess)
  1. 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?
set.seed(123)

chem_data <- chem_data %>% select(!nearZeroVar(.))
train_index_chen <- createDataPartition(chem_data$Yield , p=.8, list=F)

train_chem <-  chem_data[ train_index_chen,] 
       
  
test_chem <- chem_data[-train_index_chen,]
pls_model_chem <- train(
  Yield ~ ., data = train_chem, method = "pls",
  trControl = trainControl("cv", number = 10),
  tuneLength = 25
)
  The optimal value is shown at ncomp 3 with a R-squared of 0.6290839.
pls_model_chem$results  %>% inner_join(pls_model_chem$bestTune)
plot(pls_model_chem)

  1. 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?
  We see a lower R-squared value of 0.4838632 on the test data. This indicates that our model is good at working with data beyond what is was trained on and that we likely do not need to worry about over fitting.
  postResample(pred = predict(pls_model_chem,test_chem), obs =test_chem$Yield)
     RMSE  Rsquared       MAE 
0.7313840 0.4838632 0.6282700
  1. Which predictors are most important in the model you have trained? Do either the biological or process predictors dominate the list?
  We can see below that the process predictors dominate the list. Of the 10 shown below, only the bottom three are biological.
plot(varImp(pls_model_chem), top = 10)

  1. 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?
  The correlation plot below allows us to see that there are many processes that can potentially lead to less yield. It will be important to carefully study those variables and see how they are affecting the end result. It may that combinations of those processes are the cause and not the individual process alone as well. The processes that directly lead to more yield are also very important, but the relationships those processes have with the other processes is just as important as well!
library(DataExplorer)


  chem_data %>% select(Yield , ends_with(c("32","13","17","09","36","33","11","08")  )) %>%
  plot_correlation()