set.seed(31415)

Exercise 6.2

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

6.2(a)

Start R and use these commands to load the data:

library(AppliedPredictiveModeling) data(permeability)

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

# Code provided by the textbook
library(AppliedPredictiveModeling)
data(permeability)

6.2(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?

library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
print(length(nearZeroVar(fingerprints)))
## [1] 719
print(dim(fingerprints))
## [1]  165 1107

There are 719 predictors that are filtered out. There are 1107 fingerprints predictors, so hopefully when the ones with low frequencies are filtered out, there will be \(1107-719=388\) left for modeling.

fingerprints_filtered <- fingerprints[,-nearZeroVar(fingerprints)]
print(dim(fingerprints_filtered))
## [1] 165 388

There are 388 remaining predictors, as desired.

6.2(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 \(R^2\)?

# Code modified from page 81 of the textbook
training_rows <- createDataPartition(permeability, p = 0.8, list = FALSE)
train_predictors <- fingerprints_filtered[training_rows,]
train_permeability <- permeability[training_rows]
test_predictors <- fingerprints_filtered[-training_rows,]
test_permeability <- permeability[-training_rows]

It is not necessary to pre-process the data further since the predictors are binary variables.

# Model Creation
plsModel <- train(train_predictors, train_permeability, method = 'pls')
print(plsModel$results[plsModel$results$RMSE == min(plsModel$results$RMSE),])
##   ncomp     RMSE  Rsquared     MAE  RMSESD RsquaredSD    MAESD
## 3     3 11.89551 0.4502546 8.59405 1.50281  0.1365067 1.008602

3 components are optimal, resulting in an RMSE of approximately 11.9. This results in an \(R^2 \approx 0.4503\).

6.2(d)

Predict the response for the test set. What is the test set estimate of \(R^2\)?

# Predict the response for the test set
fingerprints_predictions <- predict(plsModel, test_predictors, ncomp = 3)
fingerprints_results <- data.frame(obs = test_permeability, pred = fingerprints_predictions)
defaultSummary(fingerprints_results)
##      RMSE  Rsquared       MAE 
## 12.450558  0.308335  8.768495

Applying the model to the test set results in \(R^2 \approx 0.3083\).

6.2(e)

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

Simple Linear Model:

# Create the model from the training data
lm_train <- as.data.frame(train_predictors)
lm_train$permeability <- train_permeability
lm_fingerprints <- lm(permeability ~ ., data = lm_train)

# Predict the response for the test set
lm_predictions <- predict(lm_fingerprints, as.data.frame(test_predictors))

# Report the accuracy of the predictions
lm_Values <- data.frame(obs = test_permeability, pred = lm_predictions)
defaultSummary(lm_Values)
##        RMSE    Rsquared         MAE 
## 28.42816313  0.01378655 17.33276372

That’s not very good.

Ridge Regression Model (Penalized):

# Create the model from the training data
library(MASS)
library(lmridge)
ridgeModel <- lm.ridge(permeability ~ ., data = lm_train)

# Create a function to use the ridge model to predict a value from a dataframe of predictors
ridge_predict <- function(model, df) {
  prediction <- c()
  for (i in 1:nrow(df)) {
    value_predict <- model$Inter
    for (j in 1:ncol(df)) {
      value_predict <- value_predict + (model$coef[j] * df[i,j])
    }
    prediction[i] <- value_predict
  }
  return(prediction)
}

# Use the function to predict the response for the test set
ridge_predictions <- ridge_predict(ridgeModel, as.data.frame(test_predictors))

# Report the accuracy of the predictions
ridge_Values <- data.frame(obs = test_permeability, pred = ridge_predictions)
defaultSummary(ridge_Values)
##         RMSE     Rsquared          MAE 
## 3.958257e+14 2.275400e-01 3.690929e+14

The R-squared isn’t as bad as the simple linear model, but the RMSE is many orders of magnitude higher. Neither of these models is as good as the original.

6.2(f)

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

None of my models have a high enough value of \(R^2\) to be valuable to the pharmaceutical company. The laboratory experiment should continue, regardless of resource expenditure.

Exercise 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 re- lationship 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 pro- cess. Improving product yield by 1% will boost revenue by approximately one hundred thousand dollars per batch:

6.3(a)

Start R and use these commands to load the data: > library(AppliedPredictiveModeling) > data(chemicalManufacturing)

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.

# Code provided by the textbook
library(AppliedPredictiveModeling)
data(ChemicalManufacturingProcess)

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

# Check for missing values
sum(is.na(ChemicalManufacturingProcess))
## [1] 106

There are 106 missing values in the data.

# Impute missing values using knn
library(impute)
chemicals_impute <- impute.knn(as.matrix(ChemicalManufacturingProcess), rng.seed = 31415)
chemicals_impute <- as.data.frame(chemicals_impute$data)
sum(is.na(chemicals_impute))
## [1] 0

There are no longer any missing values.

6.3(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?

training_rows_chemicals <- createDataPartition(chemicals_impute$Yield, p = 0.8, list = FALSE)
train_predictors_chemicals <- chemicals_impute[training_rows_chemicals,-1]
train_yield <- chemicals_impute[training_rows_chemicals, 1]
test_predictors_chemicals <- chemicals_impute[-training_rows_chemicals,-1]
test_yield <- chemicals_impute[-training_rows_chemicals, 1]

plsModel_2 <- train(train_predictors_chemicals,
                    train_yield,
                    method = 'pls',
                    preProc = c('center', 'scale'))

print(plsModel_2$results[plsModel_2$results$RMSE == min(plsModel_2$results$RMSE),])
##   ncomp     RMSE  Rsquared      MAE    RMSESD RsquaredSD     MAESD
## 3     3 1.410656 0.5410581 1.038902 0.3603746  0.1087564 0.1208303

Like with the fingerprints data the optimal model uses 3 components, this time resulting in an RMSE of approximately 1.4107 with \(R^2 \approx 0.5411\).

6.3(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?

# Predict the response for the test set
yield_predictions <- predict(plsModel_2, test_predictors_chemicals, ncomp = 3)
yield_values <- data.frame(obs = test_yield, pred = yield_predictions)
defaultSummary(yield_values)
##      RMSE  Rsquared       MAE 
## 1.8794555 0.2573087 1.2458304

Predicting the response for the test set results in an RMSE of approximately 1.879 and \(R^2 \approx 0.2573\).These values are worse than those for the training set, suggesting that the model may have overfit for the training data.

6.3(e)

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

varImp(plsModel_2, scale = FALSE)
## Warning: package 'pls' was built under R version 4.3.3
## 
## Attaching package: 'pls'
## The following object is masked from 'package:caret':
## 
##     R2
## The following object is masked from 'package:stats':
## 
##     loadings
## pls variable importance
## 
##   only 20 most important variables shown (out of 57)
## 
##                        Overall
## ManufacturingProcess32 0.11349
## ManufacturingProcess36 0.10170
## ManufacturingProcess13 0.09859
## ManufacturingProcess09 0.09418
## ManufacturingProcess17 0.09218
## ManufacturingProcess06 0.07632
## BiologicalMaterial02   0.07370
## ManufacturingProcess33 0.07049
## BiologicalMaterial06   0.06952
## BiologicalMaterial08   0.06945
## ManufacturingProcess11 0.06720
## BiologicalMaterial03   0.06702
## ManufacturingProcess30 0.06606
## BiologicalMaterial04   0.06534
## BiologicalMaterial12   0.06492
## ManufacturingProcess29 0.06410
## ManufacturingProcess31 0.06392
## BiologicalMaterial11   0.06275
## ManufacturingProcess12 0.06165
## BiologicalMaterial01   0.05866

Of the top 20 most important variables, 12 are Manufacturing Process predictors and 8 are Biological Material predictors. However, this does not tell the whole story, as all of the top 6 predictors (and 7 of the top 8) are Manufacturing Process predictors. From that, one could argue that the Manufacturing Process predictors dominate the list. On the other hand, there are only 12 Biological Material predictors and 8 of those 12 are in the top 20, so from that one could argue that the Biological Material predictors are over-represented among the top predictors.

6.3(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?

ggplot(chemicals_impute, aes(x = ManufacturingProcess32, y = Yield)) + geom_point()

ggplot(chemicals_impute, aes(x = ManufacturingProcess36, y = Yield)) + geom_point()

ggplot(chemicals_impute, aes(x = ManufacturingProcess13, y = Yield)) + geom_point()

ggplot(chemicals_impute, aes(x = ManufacturingProcess09, y = Yield)) + geom_point()

ggplot(chemicals_impute, aes(x = ManufacturingProcess17, y = Yield)) + geom_point()

ggplot(chemicals_impute, aes(x = ManufacturingProcess06, y = Yield)) + geom_point()

Of the 6 manufacturing processes at the top of the most important variables list, 3 appear to have a negative impact on yield (36, 13, and 17). This information could allow the company manufacturing the product to identify ways of improving these processes, which would increase the yield.