DATA 624 Homework 7

Author

Jonnathan Zuna

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:

# Loading required packages
library(AppliedPredictiveModeling)
library(caret)
library(pls)
library(elasticnet)
library(RANN)
data(permeability)

dim(fingerprints)
[1]  165 1107
length(permeability)
[1] 165

There are 1,107 predictors for 165 compounds.

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?

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

ncol(fp)
[1] 388

After removing the near zero variance predictors, 388 of the 1,107 remain.

c. Split the data into a training and a test set, preprocess the data, and tune a PLS model. How many latent variables are optimal and what is the corresponding resampled estimate of R squared?

set.seed(624)
idx <- createDataPartition(permeability[, 1], p = 0.8, list = FALSE)

x_train <- fp[idx, ]
x_test <- fp[-idx, ]
y_train <- permeability[idx, 1]
y_test <- permeability[-idx, 1]

133 compounds for training and 32 for testing.

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

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

plot(pls_fit)

pls_fit$results[pls_fit$results$ncomp == pls_fit$bestTune$ncomp, ]
  ncomp     RMSE  Rsquared     MAE   RMSESD RsquaredSD    MAESD
9     9 12.46889 0.4066932 9.76411 1.932523  0.2131999 1.668379

The optimal model uses 9 latent variables, with a resampled R squared of 0.407.

d. Predict the response for the test set. What is the test set estimate of R squared?

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

pls_test
     RMSE  Rsquared       MAE 
10.656898  0.530305  8.388575

The test R squared is 0.53, higher than the resampled 0.407. With only 32 test compounds, the gap is within the resampled R squared SD of 0.21, so the two estimates agree.

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

Ridge regression and an elastic net were tested next.

set.seed(624)
ridge_fit <- train(x_train, y_train,
                   method = "ridge",
                   tuneGrid = data.frame(lambda = c(0.01, 0.05, 0.1, 0.2, 0.5, 1, 2)),
                   preProc = c("center", "scale"),
                   trControl = ctrl)

ridge_fit$bestTune
  lambda
3    0.1
set.seed(624)
enet_fit <- train(x_train, y_train,
                  method = "enet",
                  tuneGrid = expand.grid(lambda = c(0.01, 0.1),
                                         fraction = c(0.01, 0.02, 0.03, 0.05, 0.1, 0.2, 0.35, 0.5, 0.75, 0.95)),
                  preProc = c("center", "scale"),
                  trControl = ctrl)

enet_fit$bestTune
  fraction lambda
6      0.2   0.01

The elastic net selected a fraction of 0.2 and a lambda of 0.01, so the lasso penalty drops many of the 388 predictors.

perf <- function(fit, name) {
  cv <- getTrainPerf(fit)
  te <- postResample(predict(fit, x_test), y_test)
  data.frame(Model = name,
             CV_RMSE = cv$TrainRMSE, CV_R2 = cv$TrainRsquared,
             Test_RMSE = te["RMSE"], Test_R2 = te["Rsquared"])
}

comp <- rbind(perf(pls_fit, "PLS"),
              perf(ridge_fit, "Ridge"),
              perf(enet_fit, "Elastic net"))
rownames(comp) <- NULL

comp
        Model  CV_RMSE     CV_R2 Test_RMSE   Test_R2
1         PLS 12.46889 0.4066932  10.65690 0.5303050
2       Ridge 12.77062 0.3941524  12.03658 0.4529012
3 Elastic net 11.95057 0.4330753  11.26462 0.4655268

The elastic net is best on resampling (RMSE 11.95, R squared 0.433), while PLS is best on the test set (RMSE 10.66, R squared 0.53). Ridge is the weakest on all four measures. The gaps are small, so neither model clearly wins, though resampling over 10 folds is more reliable than 32 test compounds and gives the elastic net a slight edge.

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

No. The models explain only about 40% to 53% of the variation, with an RMSE between 10.7 and 12.8, which is not accurate enough to replace the lab experiment. They could still help screen out compounds likely to have low permeability.

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.

a. Start R and use these commands to load the data:

data(ChemicalManufacturingProcess)

cmp <- ChemicalManufacturingProcess

dim(cmp)
[1] 176  58

176 batches, with Yield and 57 predictors.

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).

x <- cmp[, -1]
yield <- cmp$Yield

sum(is.na(x))
[1] 106
round(100 * mean(is.na(x)), 2)
[1] 1.06
miss <- colSums(is.na(x))

sort(miss[miss > 0], decreasing = TRUE)[1:5]
ManufacturingProcess03 ManufacturingProcess11 ManufacturingProcess10 
                    15                     10                      9 
ManufacturingProcess25 ManufacturingProcess26 
                     5                      5

About 1% of the cells are missing, mostly in ManufacturingProcess03 (15), 11 (10) and 10 (9). They are filled with k nearest neighbor imputation inside train in part c, so the imputation is learned from the training batches only.

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

set.seed(624)
idx <- createDataPartition(yield, p = 0.8, list = FALSE)

x_train <- x[idx, ]
x_test <- x[-idx, ]
y_train <- yield[idx]
y_test <- yield[-idx]
ctrl <- trainControl(method = "cv", number = 10)

set.seed(624)
pls_cmp <- train(x_train, y_train,
                 method = "pls",
                 tuneLength = 15,
                 preProc = c("knnImpute", "nzv", "center", "scale"),
                 trControl = ctrl)

plot(pls_cmp)

pls_cmp$results[pls_cmp$results$ncomp == pls_cmp$bestTune$ncomp, ]
  ncomp     RMSE  Rsquared      MAE    RMSESD RsquaredSD    MAESD
3     3 1.368345 0.5231989 1.059317 0.4447395  0.1702433 0.189138

PLS was used again, with imputation, near zero variance filtering, centering and scaling. The optimal model uses 3 latent variables, with a resampled RMSE of 1.368 and R squared of 0.523.

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?

cmp_test <- postResample(predict(pls_cmp, x_test), y_test)

cmp_test
     RMSE  Rsquared       MAE 
1.2015788 0.6756127 0.9535893

The test RMSE is 1.202 and R squared is 0.676, both better than the resampled values. The RMSE gap of 0.166 is smaller than the resampled RMSE SD of 0.445, so the difference is normal variation.

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

imp <- varImp(pls_cmp)

imp
pls variable importance

  only 20 most important variables shown (out of 56)

                       Overall
ManufacturingProcess32  100.00
ManufacturingProcess09   82.48
ManufacturingProcess36   81.17
ManufacturingProcess13   75.96
ManufacturingProcess17   69.82
ManufacturingProcess06   62.10
BiologicalMaterial06     57.95
BiologicalMaterial03     57.52
BiologicalMaterial02     57.42
ManufacturingProcess11   56.52
BiologicalMaterial08     55.26
ManufacturingProcess34   55.17
ManufacturingProcess33   52.69
BiologicalMaterial12     50.78
BiologicalMaterial01     49.98
BiologicalMaterial11     48.97
BiologicalMaterial04     46.99
ManufacturingProcess12   46.31
ManufacturingProcess28   44.56
ManufacturingProcess04   43.76
top <- rownames(imp$importance)[order(-imp$importance$Overall)][1:10]

table(substr(top, 1, 10))

Biological Manufactur 
         3          7

7 of the top 10 are process predictors and 3 are biological, so the process predictors dominate. The top three are ManufacturingProcess32, 09 and 36.

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?

round(cor(cmp[, top], cmp$Yield, use = "complete.obs"), 3)
                         [,1]
ManufacturingProcess32  0.581
ManufacturingProcess09  0.534
ManufacturingProcess36 -0.506
ManufacturingProcess13 -0.526
ManufacturingProcess17 -0.463
ManufacturingProcess06  0.394
BiologicalMaterial06    0.458
BiologicalMaterial03    0.462
BiologicalMaterial02    0.462
ManufacturingProcess11  0.342
featurePlot(x = cmp[, top[1:4]], y = cmp$Yield,
            plot = "scatter", type = c("p", "smooth"),
            layout = c(2, 2))

ManufacturingProcess32 (0.58) and 09 (0.53) have the strongest positive correlations with yield, and ManufacturingProcess13 (-0.53) and 36 (-0.51) the strongest negative ones. Since process predictors can be changed, raising 32 and 09 and lowering 13 and 36 are worth testing in future runs. BiologicalMaterial02 and 03 (0.46) cannot be changed but could be used to screen raw material lots. These are correlations, so any change should be tested before assuming it improves yield.