DATA624: Homework 7

library(AppliedPredictiveModeling)
library(caret)
library(dplyr)
library(corrplot)

Objective

In Kuhn and Johnson do problems 6.2 and 6.3. There are only two but they consist of many parts. Please submit a link to your Rpubs and submit the .rmd file as well.

Exercise 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:

Part: a

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 65 compounds, while permeability contains permeability response.

Part: 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?

388 predictors are left for modeling.

near_zero <- nearZeroVar(fingerprints)

fingerprints2 <- fingerprints[, -near_zero]

paste0("There are ", length(near_zero), " variables that are sparse, leaving ", ncol(fingerprints2), " predictors for modeling")

## [1] "There are 719 variables that are sparse, leaving 388 predictors for modeling"

Part: 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\)?

The PLS model is optimal with 2 latent variables. The corresponding \(R^2\) estimate is 0.496.

set.seed(15)
train_test_split <- createDataPartition(permeability, p = 0.7, list = FALSE)



train.x <- fingerprints2[train_test_split,]
train.y <- permeability[train_test_split,]

test.x <- fingerprints2[-train_test_split,]
test.y <- permeability[-train_test_split,]


plsTune <- train(train.x, train.y,
                 method = "pls",
                 metric = "Rsquared",
                 tuneLength = 20,
                 trControl = trainControl(method = "cv", number = 10),
                 preProc = c("center", "scale"))

plsTune

## Partial Least Squares 
## 
## 117 samples
## 388 predictors
## 
## Pre-processing: centered (388), scaled (388) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 105, 105, 106, 105, 106, 106, ... 
## Resampling results across tuning parameters:
## 
##   ncomp  RMSE      Rsquared   MAE     
##    1     12.54106  0.3179367  9.621565
##    2     11.14859  0.4962474  8.120444
##    3     11.42797  0.4699825  8.690563
##    4     11.74228  0.4531516  8.856052
##    5     11.83513  0.4455282  8.886757
##    6     11.79029  0.4718753  8.777539
##    7     11.92228  0.4439083  8.985305
##    8     12.01575  0.4421544  9.223336
##    9     12.01172  0.4604221  9.304520
##   10     12.20516  0.4528909  9.374700
##   11     12.14180  0.4619462  9.522388
##   12     12.25611  0.4554170  9.397265
##   13     12.20097  0.4651625  9.477893
##   14     12.49457  0.4470675  9.625209
##   15     12.41119  0.4535839  9.486632
##   16     12.31010  0.4581961  9.407818
##   17     12.22828  0.4717548  9.426583
##   18     12.38673  0.4573004  9.559131
##   19     12.37614  0.4563370  9.559319
##   20     12.35271  0.4529230  9.587720
## 
## Rsquared was used to select the optimal model using the largest value.
## The final value used for the model was ncomp = 2.

ggplot(plsTune) +
  xlab("Number of Predictors")

Part: d

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

The estimated \(R^2\) for the test set is 0.45 which is slightly worse than the training set.

pls.predictions <- predict(plsTune, newdata = test.x)

postResample(pred = pls.predictions, obs = test.y)

##       RMSE   Rsquared        MAE 
## 12.5469842  0.4504197  9.4247413

Part: e

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

None of the models that were discussed in this chapter outperformed the PLS model. The closest model to the PLS model, in terms of predictive performance, was the pcr model with an \(R^2\) of 0.436.

pcrTune <- train(train.x, train.y,
                   method = "pcr",
                   metric = "Rsquared",
                   tuneLength = 10,
                   trControl = trainControl("cv", number = 10),
                   preProc=c('center', 'scale')
)


ridgeGrid <- data.frame(.lambda = seq(0.001, .1, length = 10))

ridgeTune <- train(train.x, train.y,
                   method = "ridge",
                   metric = "Rsquared",
                   tuneGrid = ridgeGrid,
                   tuneLength = 10,
                   trControl = trainControl(method = "cv", number = 10),
                   preProc = c("center", "scale"))



enetGrid <- expand.grid(.lambda = c(0.001, 0.01, .1), 
                        .fraction = seq(.05, 1, length = 10))

enetTune <- train(train.x, train.y,
                  method = "enet",
                  metric = "Rsquared",
                  tuneGrid = enetGrid,
                  tuneLength = 10,
                  trControl = trainControl(method = "cv", number = 10),
                  preProc = c("center", "scale"))

pcr.predictions <- predict(pcrTune, newdata = test.x)
pcr.model <- postResample(pred = pcr.predictions, obs = test.y)

ridge.predictions <- predict(ridgeTune, newdata = test.x)
ridge.model <- postResample(pred = ridge.predictions, obs = test.y)

enet.predictions <- predict(enetTune, newdata = test.x)
enet.model <- postResample(pred = enet.predictions, obs = test.y)

pls.predictions <- predict(plsTune, newdata = test.x)
pls.model <-postResample(pred = pls.predictions, obs = test.y)

rbind(pcr.model, ridge.model, enet.model, pls.model)

##                 RMSE  Rsquared       MAE
## pcr.model   12.81788 0.4359792  9.600574
## ridge.model 15.95178 0.3028220 11.625050
## enet.model  14.20037 0.3647488 10.246233
## pls.model   12.54698 0.4504197  9.424741

Part: f

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

I would not recommend the permeability laboratory experiment to be replaced with any of the models created above. The Rsquared values show that at best, only about 45% of the variability is being explained by the model. Other models only explain about 30% of the variability. This is not enough to confidently replace the experiment.

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

Part: a

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.

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

I chose to impute the missing values using knn imputation.

imputed.knn <- preProcess(ChemicalManufacturingProcess,
           method = "knnImpute",
           k = sqrt(nrow(ChemicalManufacturingProcess))
           )

imputed.data <- predict(imputed.knn, ChemicalManufacturingProcess)

Part: 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?

To preprocess the data, I removed near zero variables. I then created a 70% - 30% train test split. A ridge regression model was created to model the data. The optimal value of the performance metric are provided below. This model achieved the highest R squared and a relatively low RMSE.

near_zero <- nearZeroVar(imputed.data)

imputed.data <- imputed.data[, -near_zero]

set.seed(15)
train_test_split <- createDataPartition(permeability, p = 0.7, list = FALSE)

train.data <- imputed.data[train_test_split,]
test.data <- imputed.data[-train_test_split,]


ridgeGrid <- data.frame(.lambda = seq(0.001, 1.1, length = 20))

ridgeTune <- train(Yield~., train.data,
                   method = "ridge",
                   metric = "Rsquared",
                   tuneGrid = ridgeGrid,
                   tuneLength = 20,
                   trControl = trainControl(method = "cv", number = 10),
                   preProc = c("center", "scale"))

ridgeTune

## Ridge Regression 
## 
## 117 samples
##  56 predictor
## 
## Pre-processing: centered (56), scaled (56) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 105, 106, 105, 105, 105, 105, ... 
## Resampling results across tuning parameters:
## 
##   lambda      RMSE      Rsquared   MAE      
##   0.00100000  3.479676  0.4157273  1.4884010
##   0.05884211  1.314173  0.5333697  0.7854032
##   0.11668421  1.194465  0.5495661  0.7409605
##   0.17452632  1.143323  0.5607429  0.7199977
##   0.23236842  1.119668  0.5683883  0.7120229
##   0.29021053  1.110793  0.5737079  0.7114491
##   0.34805263  1.111056  0.5774683  0.7171152
##   0.40589474  1.117491  0.5801509  0.7252280
##   0.46373684  1.128363  0.5820655  0.7361973
##   0.52157895  1.142572  0.5834179  0.7490055
##   0.57942105  1.159378  0.5843486  0.7636324
##   0.63726316  1.178254  0.5849568  0.7788882
##   0.69510526  1.198810  0.5853141  0.7950401
##   0.75294737  1.220745  0.5854738  0.8125181
##   0.81078947  1.243819  0.5854757  0.8305257
##   0.86863158  1.267837  0.5853508  0.8494645
##   0.92647368  1.292637  0.5851231  0.8686982
##   0.98431579  1.318084  0.5848117  0.8884473
##   1.04215789  1.344063  0.5844319  0.9084056
##   1.10000000  1.370474  0.5839959  0.9282704
## 
## Rsquared was used to select the optimal model using the largest value.
## The final value used for the model was lambda = 0.8107895.

ridgeTune$results[15,]

##       lambda     RMSE  Rsquared       MAE   RMSESD RsquaredSD     MAESD
## 15 0.8107895 1.243819 0.5854757 0.8305257 1.345039  0.2721869 0.5359635

Part: 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?

The test set results are slightly worse than the training set which is to be expected. The model performance seems to be similar on both training and testing, indicating that the model has not been overfit or underfit.

ridge.predictions <- predict(ridgeTune, newdata = test.data)
ridge.model <- postResample(pred = ridge.predictions, obs = test.data$Yield)

ridge.model

##       RMSE   Rsquared        MAE 
## 2.59329505 0.05771302 1.45880691

Part: e

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

Below is a list of the top 20 variables from the model. The top 5 variables are a combination of biological and process variables. The top 20 contains 11 process variables and 9 biological variables. It appears that there is a fairly equal split of importance between the two variable types.

varImp(ridgeTune)

## loess r-squared variable importance
## 
##   only 20 most important variables shown (out of 56)
## 
##                        Overall
## ManufacturingProcess32  100.00
## BiologicalMaterial06     86.35
## ManufacturingProcess13   86.26
## BiologicalMaterial12     82.24
## BiologicalMaterial03     76.27
## ManufacturingProcess17   73.78
## ManufacturingProcess36   72.13
## ManufacturingProcess09   70.38
## ManufacturingProcess06   64.78
## ManufacturingProcess31   64.59
## BiologicalMaterial02     60.68
## ManufacturingProcess11   55.25
## BiologicalMaterial11     54.54
## ManufacturingProcess30   52.98
## BiologicalMaterial04     48.73
## ManufacturingProcess33   48.09
## BiologicalMaterial09     45.89
## BiologicalMaterial08     40.99
## BiologicalMaterial01     37.82
## ManufacturingProcess29   35.88

Part: 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?

Most of the important variables have a moderate correlation with the target variable, Yield. There are a handful of biological material and manufacturing processes that are correlated with eachother. The relationships between the response variable and the other predictor variables can be leveraged into maximizing variables that are positively correlated with Yield and other predictors and reduce variables that have a negative relationship with Yield and other predictors.

important.vars <- imputed.data %>%
  select("Yield", "ManufacturingProcess32", "BiologicalMaterial06", "ManufacturingProcess13",
         "BiologicalMaterial12", "BiologicalMaterial03", "ManufacturingProcess17", "ManufacturingProcess36",
         "ManufacturingProcess09", "ManufacturingProcess06", "ManufacturingProcess31")

cor.matrix <- cor(important.vars)

corrplot(corr = cor.matrix, tl.col = 'black',  type = 'lower', diag = FALSE)