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:

library(AppliedPredictiveModeling)
## Warning: package 'AppliedPredictiveModeling' was built under R version 4.1.3
data(permeability)

The matrix fingerprints contains the 1,107 binary molecular predictors for the 165 compounds, while permeability contains permeability response. (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?

fdf <- as.data.frame(fingerprints)
print(paste('Total predictors:', ncol(fdf)))
## [1] "Total predictors: 1107"
 total <- nearZeroVar(fingerprints)
x <- fingerprints[,-total]

length(total) %>% paste('  - number of predictors removed. ', dim(x)[2], ' -number of predictors left .' ,sep='') %>% print()
## [1] "719  - number of predictors removed. 388 -number of predictors left ."
  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 resamp11led estimate of R2?

Using createDataPartition function to split data.

set.seed(666)

index <- createDataPartition(permeability, p=0.8, list=FALSE)
trainx <- x[index, ]
trainy <- permeability[index, ]
testx <- x[-index, ]
testy <- permeability[-index, ]
set.seed(999)
fit <- train(x=trainx,
                y=trainy, 
                method='pls',
                metric='Rsquared',
                tuneLength=20,
                trControl=trainControl(method='cv'),
                preProcess=c('center', 'scale')
                )

results <- fit$results
fit
## Partial Least Squares 
## 
## 133 samples
## 388 predictors
## 
## Pre-processing: centered (388), scaled (388) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 119, 119, 121, 120, 121, 120, ... 
## Resampling results across tuning parameters:
## 
##   ncomp  RMSE      Rsquared   MAE      
##    1     12.53693  0.3609049   9.615214
##    2     11.15229  0.4541524   7.986096
##    3     11.03294  0.4991362   8.334416
##    4     11.02691  0.5046469   8.621354
##    5     10.40376  0.5543281   8.042310
##    6     10.16324  0.5737500   7.932905
##    7     10.13563  0.5759777   8.042927
##    8     10.35977  0.5598915   8.246622
##    9     10.64398  0.5494430   8.410117
##   10     10.88973  0.5384905   8.609805
##   11     11.05270  0.5291800   8.655664
##   12     11.39426  0.5105068   8.972148
##   13     11.88588  0.4870079   9.365672
##   14     12.26186  0.4664294   9.554427
##   15     12.43766  0.4568896   9.706843
##   16     12.87842  0.4354995  10.043766
##   17     13.25981  0.4150657  10.271586
##   18     13.40412  0.4096239  10.314111
##   19     13.81138  0.3938779  10.597023
##   20     14.15448  0.3817126  10.967574
## 
## Rsquared was used to select the optimal model using the largest value.
## The final value used for the model was ncomp = 7.
plot(fit)

ncomp to be 7, with the maximum R^2 being 0.5759777.

  1. Predict the response for the test set. What is the test set estimate of R2?
p <- predict(fit, newdata=testx)
postResample(pred=p, obs=testy)
##       RMSE   Rsquared        MAE 
## 12.3106458  0.4827401  9.4376024

R^2 is 0.4827401.

  1. Try building other models discussed in this chapter. Do any have better predictive performance?
rgrid <- data.frame(.lambda = seq(0, .1, length = 15))
rFit <- train(trainx, trainy, method = "ridge",
                     tuneGrid = rgrid, trControl = trainControl(method='cv'),
                     preProc = c("center", "scale"))

rFit
## Ridge Regression 
## 
## 133 samples
## 388 predictors
## 
## Pre-processing: centered (388), scaled (388) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 119, 121, 119, 121, 120, 120, ... 
## Resampling results across tuning parameters:
## 
##   lambda       RMSE      Rsquared   MAE      
##   0.000000000  15.05851  0.3326972  11.058068
##   0.007142857  33.50735  0.3401164  21.821264
##   0.014285714  67.68943  0.3602361  40.904036
##   0.021428571  12.97709  0.4127463   9.785239
##   0.028571429  12.84830  0.4251292   9.687384
##   0.035714286  12.48182  0.4545915   9.494972
##   0.042857143  12.30221  0.4708278   9.377671
##   0.050000000  12.18945  0.4823071   9.303019
##   0.057142857  12.07435  0.4940735   9.230316
##   0.064285714  11.97154  0.5035379   9.146136
##   0.071428571  11.89493  0.5115853   9.090353
##   0.078571429  11.84263  0.5183763   9.052239
##   0.085714286  11.81776  0.5225097   9.027003
##   0.092857143  11.74631  0.5296352   8.966917
##   0.100000000  11.70609  0.5345796   8.935717
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 0.1.
pred <- predict(rFit, newdata = testx)
postResample(pred = pred, obs = testy)
##       RMSE   Rsquared        MAE 
## 12.4873335  0.4889747  9.3216417

R^2 is 0.4889747.

  1. Would you recommend any of your models to replace the permeability laboratory experiment?
qqnorm(p, main = "PLS")
qqline(p)

qqnorm(pred, main = "Ridge-Regression")
qqline(pred)

Although the R2 for both models are quite low, the examonation of the residuals should be performed and to see how well is each model performing. From the grapghs above it is visible that the PLS model is id doing a better job fitting the data as opposed to the ridge-regression model, which are a bit of above and under, and over predictors. But it is something that can be improved upon. I do not think would replace the permeability laboratory experiment at this point, but potentially there is a possibility.

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:

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

(impute <- preProcess(ChemicalManufacturingProcess[,-c(1)], method=c('bagImpute')))
## Created from 152 samples and 57 variables
## 
## Pre-processing:
##   - bagged tree imputation (57)
##   - ignored (0)
cp <- predict(impute, ChemicalManufacturingProcess[,-c(1)])
  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(999)
inx <- createDataPartition(ChemicalManufacturingProcess$Yield, p=0.8, list=FALSE)
trainx <- cp[inx, ]
trainy <- ChemicalManufacturingProcess$Yield[inx]
testx <- cp[-inx, ]
testy <- ChemicalManufacturingProcess$Yield[-inx]
set.seed(666)
efit <- train(x=trainx,
                 y=trainy,
                 method='enet',
                 metric='RMSE',
                 tuneGrid=expand.grid(.fraction = seq(0, 1, by=0.1), 
                                      .lambda = seq(0, 1, by=0.1)),
                 trControl=trainControl(method='cv'),
                 preProcess=c('center','scale')
                  )
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
plot(efit)

The elastic net model is tuned using 10-fold cross validation with parameters lambda ranging from 0 to 1, and fraction ranging from 0 to 1. The metric used to decide is the RMSE.

  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?
epred <- predict(efit, newdata=testx)
(pred <- postResample(pred=epred, obs=testy))
##      RMSE  Rsquared       MAE 
## 1.0510904 0.7606085 0.8552145

The RMSE is 1.0516956 and comparing it with the resampled performance metric on the training set indicates that it is lower, which suggests that test date does better than the training set.

  1. Which predictors are most important in the model you have trained? Do either the biological or process predictors dominate the list?
library(elasticnet)
## Loading required package: lars
## Loaded lars 1.2
(coff <- predict.enet(efit$finalModel, s=efit$bestTune[1, "fraction"], type="coef", mode="fraction")$coefficients)
##   BiologicalMaterial01   BiologicalMaterial02   BiologicalMaterial03 
##             0.00000000             0.09484693             0.01498162 
##   BiologicalMaterial04   BiologicalMaterial05   BiologicalMaterial06 
##             0.00000000             0.00000000             0.12349005 
##   BiologicalMaterial07   BiologicalMaterial08   BiologicalMaterial09 
##             0.00000000             0.00000000             0.00000000 
##   BiologicalMaterial10   BiologicalMaterial11   BiologicalMaterial12 
##             0.00000000             0.00000000             0.00000000 
## ManufacturingProcess01 ManufacturingProcess02 ManufacturingProcess03 
##             0.00000000             0.00000000             0.00000000 
## ManufacturingProcess04 ManufacturingProcess05 ManufacturingProcess06 
##             0.00000000             0.00000000             0.05914029 
## ManufacturingProcess07 ManufacturingProcess08 ManufacturingProcess09 
##             0.00000000             0.00000000             0.27856224 
## ManufacturingProcess10 ManufacturingProcess11 ManufacturingProcess12 
##             0.00000000             0.02368239             0.00000000 
## ManufacturingProcess13 ManufacturingProcess14 ManufacturingProcess15 
##            -0.27150743             0.00000000             0.00000000 
## ManufacturingProcess16 ManufacturingProcess17 ManufacturingProcess18 
##             0.00000000            -0.13519037             0.00000000 
## ManufacturingProcess19 ManufacturingProcess20 ManufacturingProcess21 
##             0.00000000             0.00000000             0.00000000 
## ManufacturingProcess22 ManufacturingProcess23 ManufacturingProcess24 
##             0.00000000             0.00000000             0.00000000 
## ManufacturingProcess25 ManufacturingProcess26 ManufacturingProcess27 
##             0.00000000             0.00000000             0.00000000 
## ManufacturingProcess28 ManufacturingProcess29 ManufacturingProcess30 
##             0.00000000             0.00000000             0.00000000 
## ManufacturingProcess31 ManufacturingProcess32 ManufacturingProcess33 
##             0.00000000             0.58201107             0.00000000 
## ManufacturingProcess34 ManufacturingProcess35 ManufacturingProcess36 
##             0.00000000             0.00000000            -0.29453743 
## ManufacturingProcess37 ManufacturingProcess38 ManufacturingProcess39 
##            -0.03059759             0.00000000             0.00000000 
## ManufacturingProcess40 ManufacturingProcess41 ManufacturingProcess42 
##             0.00000000             0.00000000             0.00000000 
## ManufacturingProcess43 ManufacturingProcess44 ManufacturingProcess45 
##             0.00000000             0.00000000             0.00000000
coff_sort <- abs(coff)
coff_sort <- coff_sort[coff_sort>0]

coff_sort
##   BiologicalMaterial02   BiologicalMaterial03   BiologicalMaterial06 
##             0.09484693             0.01498162             0.12349005 
## ManufacturingProcess06 ManufacturingProcess09 ManufacturingProcess11 
##             0.05914029             0.27856224             0.02368239 
## ManufacturingProcess13 ManufacturingProcess17 ManufacturingProcess32 
##             0.27150743             0.13519037             0.58201107 
## ManufacturingProcess36 ManufacturingProcess37 
##             0.29453743             0.03059759

The highest coefficients are all from the ManufacturingProcess predictors. It appears that ManufacturingProcess are more important than the Biological Materials.

  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 highest scoring Biological Material was at 0.082433754. Overall, as we know the biological materials cannot be changed during the refinement process. But looking into what ingredients/materials are more vital would assist to make sure a higher yield as the company can focus on obtaining higher quality ingredients of those materials. Also, being aware of the most important manufacturing process steps helps the company to find out where they can start tuning the procedure/model.