Data 624 HW 7

library(knitr)
library(rmdformats)

## Global options
options(max.print="85")
opts_chunk$set(cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               comment=NA,
               message=FALSE,
               warning=FALSE)
opts_knit$set(width=35)

library(AppliedPredictiveModeling)
library(caret)

## Loading required package: lattice

## Loading required package: ggplot2

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 4.0-2

library(RANN)
library(tibble)

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:

data("permeability")
str(permeability)

 num [1:165, 1] 12.52 1.12 19.41 1.73 1.68 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:165] "1" "2" "3" "4" ...
  ..$ : chr "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?

fingerprints_df <- as.data.frame(fingerprints)
paste("Total predictors:", ncol(fingerprints_df))

[1] "Total predictors: 1107"

cat(paste("Predictors left for modeling is: ", ncol(fingerprints_df[, -nearZeroVar(fingerprints)])))

Predictors left for modeling is:  388

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

str(permeability)

 num [1:165, 1] 12.52 1.12 19.41 1.73 1.68 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:165] "1" "2" "3" "4" ...
  ..$ : chr "permeability"

perm <- cbind(permeability, fingerprints_df[, -nearZeroVar(fingerprints)])

# Make this reproducible
set.seed(43)

train <- createDataPartition(perm$permeability, times = 1, p = 0.8, list = FALSE)
train_df <- perm[train, ]
test_df <- perm[-train, ]

# Partial Least Squares Model
pls_model <- train(permeability ~ ., data = train_df, method = "pls", tuneLength = 24, 
    trControl = trainControl(method = "repeatedcv", repeats = 8), preProcess = c("center"))

# Plot model RMSE against number of variables to the optimal number of variables
p_title <- paste("Training Set attains its optimal model with the smallest RMSE  at", 
    pls_model$bestTune$ncomp, "variables")
ggplot(pls_model) + xlab("Number of Variables") + ggtitle(p_title) + theme(plot.title = element_text(size = 8))

pls_model$results %>% filter(ncomp == pls_model$bestTune$ncomp) %>% select(ncomp, 
    Rsquared)

  ncomp  Rsquared
1     6 0.5315704

Ans: The corresponding resampled esimate of R² is 0.52 if we pick the number of variables of 3.

(d)

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

# Make predictions
pred <- predict(pls_model, test_df)
# Model performance metrics
postResample(pred = pred, obs = test_df$permeability)

      RMSE   Rsquared        MAE 
14.4140798  0.2448084 10.2056286

Ans: The test set estimate of R² is 0.2413638.

(e)

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

# to find the right lambda using cv.glmnet
x_train <- model.matrix(permeability ~ ., data = train_df)
x_test <- model.matrix(permeability ~ ., data = test_df)
cv.glmnet <- cv.glmnet(x_train, train_df$permeability, alpha = 0)
ridge_model <- glmnet(x_train, train_df$permeability, alpha = 0, lambda = cv.glmnet$lambda.min)
# Make predictions
pred_ridge <- predict(ridge_model, x_test)
# Model performance metrics
postResample(pred = pred_ridge, obs = test_df$permeability)

      RMSE   Rsquared        MAE 
12.9909345  0.2995253  8.8120660

# to find the right lambda using cv.glmnet
cv.lasso_reg <- cv.glmnet(x_train, train_df$permeability, alpha = 1)
lasso_reg_model <- glmnet(x_train, train_df$permeability, alpha = 1, lambda = cv.lasso_reg$lambda.min)
# Make predictions
pred_lassoR <- predict(lasso_reg_model, x_test)
# Model performance metrics
postResample(pred = pred_lassoR, obs = test_df$permeability)

      RMSE   Rsquared        MAE 
13.8538539  0.2857556  9.6811683

# training the elastic net regression model using train
elas_net_model <- train(permeability ~ ., data = train_df, method = "glmnet", trControl = trainControl("repeatedcv", 
    repeats = 8), tuneLength = 4)

# Make predictions
elas_net_pred <- predict(elas_net_model, x_test)
# Model performance metrics
postResample(pred = elas_net_pred, obs = test_df$permeability)

      RMSE   Rsquared        MAE 
12.5521122  0.3566109  8.8283210

Ans: Yes, I do see better predictive performance. Elastic Net Regression does have RSquared of 35.67%, which is the best among all 4 models. What it means the model explains that much of the variance. In terms of RMSE, Elastic Net also has the lowest value among all 4 models. Ridge regression has the lowest MAE among all models and it edged Elastic Net Regression slightly, 8.81 vs 8.83.

(f)

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

Ans: I would not recommend any of my models to replace the permeability laboratory experiment as the highest RSquared is merely ~36%. The models do not explain enough of the variability to make predictions for permeability reliable enough to potentially reduce the need for the assay. So I’d continue to run the experiments.

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)

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

Status Quo of ChemicalManufacturingProcess

visdat::vis_miss(ChemicalManufacturingProcess) + theme(axis.text.x = element_text(color = "grey20", 
    size = 4, angle = 45, hjust = 0, vjust = 0.8, face = "plain"))

Used bagImpute method to impute the missing values.

bagImpute <- preProcess(ChemicalManufacturingProcess[, -c(1)], method = c("bagImpute"))

CMP <- predict(bagImpute, ChemicalManufacturingProcess)
visdat::vis_miss(CMP) + theme(axis.text.x = element_text(color = "grey20", size = 4, 
    angle = 45, hjust = 0, vjust = 0.8, face = "plain"))

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

Ans:

Step 0: Set a seed

# Make this reproducible
set.seed(45)

Step 1: Removing the zero variance variables from the dataset CMP

CMP_2 <- CMP %>% select(-one_of(nearZeroVar(., names = TRUE)))

Step 2: Create train and test split

train <- createDataPartition(CMP_2$Yield, times = 1, p = 0.8, list = FALSE)
CMP2_train_df <- CMP_2[train, ]
CMP2_test_df <- CMP_2[-train, ]

Step 3: Train the model

pls_model <- train(Yield ~ ., data = CMP2_train_df, method = "pls", tuneLength = 24, 
    trControl = trainControl(method = "repeatedcv", repeats = 8), preProcess = c("center", 
        "scale"))

# Plot model RMSE against number of variables to the optimal number of variables
p_title <- paste("Training Set attains its optimal model with the smallest RMSE  at", 
    pls_model$bestTune$ncomp, "variables")
ggplot(pls_model) + xlab("Number of Variables") + ggtitle(p_title) + theme(plot.title = element_text(size = 8))

Step 4: Examining model results

pls_model$results %>% filter(ncomp == pls_model$bestTune$ncomp) %>% select(ncomp, 
    Rsquared)

  ncomp  Rsquared
1     3 0.6237492

The optimal value of the performance metrics is RSquared = 62.48%

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

# Make predictions
pred <- predict(pls_model, CMP2_test_df)
# Model performance metrics
postResample(pred = pred, obs = CMP2_test_df$Yield)

    RMSE Rsquared      MAE 
1.525349 0.382110 1.285310

Ans: The RSquared has dropped off dramatically in the test set, to 38.43%.

(e)

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

Here is a list of top 10 most important predictors:

varImp_2 <- varImp(pls_model)$importance %>% as.data.frame() %>% rownames_to_column() %>% 
    arrange(desc(Overall)) %>% top_n(10)


varImp_2 %>% mutate(rowname = forcats::fct_inorder(rowname)) %>% ggplot() + geom_col(aes(x = rowname, 
    y = Overall)) + coord_flip() + theme_bw()

Ans: Yes, ManufacturingProcess dominated the list.

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

CMP_2 %>% select(as.vector(cbind(c("Yield", varImp_2$rowname)))) %>% cor() %>% corrplot::corrplot(tl.cex = 0.4)

Ans: What stood out the most is there seems to be a strong relationship between the response Yield and the following predictors/variables,

ManufacturingProcess32 (positive corr.)
ManufacturingProcess09 (positive corr.)
ManufacturingProcess13 (inverse correlation)
ManufacturingProcess36 (inverse correlation)
BiologicalMaterial02 (positive corr.)
BiologicalMaterial06 (positive corr.)