Data 624 - Homework 7

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:

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.

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?
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 R2?
Predict the response for the test set. What is the test set estimate of R2?
Try building other models discussed in this chapter. Do any have better predictive performance?
Would you recommend any of your models to replace the permeability laboratory experiment?

a.

library(AppliedPredictiveModeling)
data(permeability)

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

b.

Filter out these predictors:

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

Seems like there are 719 predictors that we can drop because they have near zero variance. This leaves us with 388 predictors.

c.

Lets look at pairwise correlation amongst the predictors:

correlations <- cor(filtered_fingerprints)
plot(correlations)

Lets remove them with threshold .90, noticed when \(R^2\) decreases and RMSE increases
Find highly correlated preditors and remove:

highCorr <- findCorrelation(correlations, cutoff = .9)
filtered_fingerprints <- filtered_fingerprints[, -highCorr]
correlations <- cor(filtered_fingerprints)
corrplot::corrplot(correlations, order = "hclust")

After cleaning the correlation matrix, it looks better now
Split the data into a training and test set:

train_fingerprints <- filtered_fingerprints[1:124, ]
test_fingerprints <- filtered_fingerprints[1:124, ]
train_permeability <- permeability[1:124, ]
test_permeability <- permeability[1:124, ]

Train a PLS Model:

model.pls <- train(train_fingerprints, train_permeability, method = "pls", tuneLength = 10, trControl = trainControl(method = "cv"))

Plot PLS Model to see at which number of components we minimize error:

plot(model.pls, main="RMSE Error of PLS Model vs Number of Components")

model.pls$results[model.pls$results$ncomp == 5, 'Rsquared']

## [1] 0.4746335

We see from the above plot that the number of components is 5 when we minimize error, with an \(R^2\) = 0.4746335

d.

Predict the permeability based on test set of fingerprints and print out predictions + outcomes:

predict_permeability = predict(model.pls, test_fingerprints)
plot(predict_permeability, test_permeability, main="Observed versus Predicted Permeability from PLS Model (n=5)", xlab="Predicted Permeability", ylab="Observed Permeability")
abline(0, 1, col='red')
text(0, 30, paste("R^2 = ", round(cor(test_permeability, predict_permeability)^2, 2)))
text(0, 27, paste("RMSE = ", round(sqrt(sum((test_permeability - predict_permeability)^2)), 2)))

The graph above shows a plot of observed values versus predicted values for permeability.
The red line shown is the y = x line: the closer the values are to the red line, the better it is.
\(R^2\) = 0.57

e.

Using some ideas from this link, we’ll use glmnet to build a model penalizing both the L1-norm and L2-norm.

glm.model <- glmnet(train_fingerprints, train_permeability, family="gaussian", alpha=0.5, lambda=0.001)
permeability.glm.predict <- predict(glm.model, test_fingerprints)
plot(permeability.glm.predict, test_permeability, main="Observed versus Predicted Permeability from GLM Model", xlab="Predicted Permeability", ylab="Observed Permeability")
abline(0, 1, col='red')
text(0, 30, paste("R^2 = ", round(cor(test_permeability, permeability.glm.predict)^2, 2)))
text(0, 25, paste("RMSE = ", round(sqrt(sum((test_permeability - permeability.glm.predict)^2)), 2)))

f.

The GLM Net model since it has a high \(R^2\).
But we need to make sure to tune this to its best and confirm the data preprocessing steps taken make sense.

6.3. A chemical manufacturing process for a pharmaceutical product was discussed in Sect. 1.4. In this problem, the objective is the 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.

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.

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).
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?
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?
Which predictors are most important in the model you have trained? Do either the biological or process predictors dominate the list?
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?

a.

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.

b.

missmap(ChemicalManufacturingProcess)

It appears that the missing values are quite minor and seem to be missing at random.

set.seed(1234)
registerDoParallel(cores=3)
chem_imput <- missForest(ChemicalManufacturingProcess, ntree = 100, parallelize = "forests")$ximp

##   missForest iteration 1 in progress...done!
##   missForest iteration 2 in progress...done!
##   missForest iteration 3 in progress...done!
##   missForest iteration 4 in progress...done!

ds <- list(org_ds = chem_imput)

c.

Lets explore the distribution of the variables:

multi.hist(ds$org_ds[1:8])

multi.hist(ds$org_ds[9:16])

multi.hist(ds$org_ds[17:24])

multi.hist(ds$org_ds[25:32])

The distributions of several of the variables are noticeably skewed, so we will obtain the box cox transformation lambda

variables.to.transform <- c("BiologicalMaterial04", "ManufacturingProcess01",
                            "ManufacturingProcess02","ManufacturingProcess03",
                            "ManufacturingProcess06", "ManufacturingProcess07",
                            "ManufacturingProcess08","ManufacturingProcess12",
                            "ManufacturingProcess16","ManufacturingProcess18") 

results <- NULL
for (i in 1:length(variables.to.transform)){
  results[[i]] <- list(variable = variables.to.transform[[i]],
                       lambda = BoxCox.lambda(ds$org_ds[,variables.to.transform[[i]]]))}
results

## [[1]]
## [[1]]$variable
## [1] "BiologicalMaterial04"
## 
## [[1]]$lambda
## [1] -0.9999242
## 
## 
## [[2]]
## [[2]]$variable
## [1] "ManufacturingProcess01"
## 
## [[2]]$lambda
## [1] 1.999959
## 
## 
## [[3]]
## [[3]]$variable
## [1] "ManufacturingProcess02"
## 
## [[3]]$lambda
## [1] 1.999959
## 
## 
## [[4]]
## [[4]]$variable
## [1] "ManufacturingProcess03"
## 
## [[4]]$lambda
## [1] -0.9999242
## 
## 
## [[5]]
## [[5]]$variable
## [1] "ManufacturingProcess06"
## 
## [[5]]$lambda
## [1] -0.9999242
## 
## 
## [[6]]
## [[6]]$variable
## [1] "ManufacturingProcess07"
## 
## [[6]]$lambda
## [1] -0.9999242
## 
## 
## [[7]]
## [[7]]$variable
## [1] "ManufacturingProcess08"
## 
## [[7]]$lambda
## [1] 1.999924
## 
## 
## [[8]]
## [[8]]$variable
## [1] "ManufacturingProcess12"
## 
## [[8]]$lambda
## [1] 4.102259e-05
## 
## 
## [[9]]
## [[9]]$variable
## [1] "ManufacturingProcess16"
## 
## [[9]]$lambda
## [1] 1.999959
## 
## 
## [[10]]
## [[10]]$variable
## [1] "ManufacturingProcess18"
## 
## [[10]]$lambda
## [1] 1.999959

We have our lambda transformations and we can apply them to our data set then set our test and train splices.

trans.data.set <- ds$org_ds %>% 
                  mutate(BiologicalMaterial04.boxcoxtrans = BoxCox(BiologicalMaterial04, -0.99992424816297),
                         ManufacturingProcess01.boxcoxtrans = BoxCox(ManufacturingProcess01, 1.99995900720725),
                         ManufacturingProcess02.boxcoxtrans = BoxCox(ManufacturingProcess02, 1.99995900720725),
                         ManufacturingProcess03.boxcoxtrans = BoxCox(ManufacturingProcess03, -0.99992424816297),
                         ManufacturingProcess06.boxcoxtrans = BoxCox(ManufacturingProcess06, -0.99992424816297),
                         ManufacturingProcess07.boxcoxtrans = BoxCox(ManufacturingProcess07,-0.99992424816297),
                         ManufacturingProcess08.boxcoxtrans = BoxCox(ManufacturingProcess08, 1.99992424816297),
                         ManufacturingProcess12.boxcoxtrans = BoxCox(ManufacturingProcess12, 0.0000410225926326142),
                         ManufacturingProcess16.boxcoxtrans = BoxCox(ManufacturingProcess16, 1.99995900720725),
                         ManufacturingProcess18.boxcoxtrans = BoxCox(ManufacturingProcess18, 1.99995900720725))

Splitting the test and training data set:

set.seed(123)
ds <- c(ds, list(trans.data.set = trans.data.set))  %>% 
            c(., list(train.test.split = floor(0.75 * nrow(chem_imput)))) %>% 
            c(., train.index = list(sample(seq_len(nrow(chem_imput)), size = .$train.test.split))) %>%
            c(., train = list(as.data.frame(trans.data.set[.$train.index,])), test = list(as.data.frame(trans.data.set[-.$train.index,])))

Now we can proceed with the model of choice which is the Partial Least Squares. Note that the evaluation below relied heavily on this rpubs post- https://rpubs.com/omicsdata/pls.

chem_plsFit <- plsr(Yield ~ ., data =  ds$train, validation = "CV")
validationplot(chem_plsFit, val.type="RMSEP")

We see a steep increase at 10 components followed by many mountains and valleys. Our intention is to select the most optimal number of components to minimize the RMSE.

pls.RMSEP <- RMSEP(chem_plsFit, estimate="CV") 
plot(pls.RMSEP, main="RMSEP PLS Solubility", xlab="components")
min_comp <- which.min(pls.RMSEP$val)
points(min_comp, min(pls.RMSEP$val), pch=1, col="red", cex=1.5)

Our optimal number of components is 4.

d.

Our results from our training set is:

chem_plsPredict_train <-  predict(chem_plsFit, data = ds$train, ncomp = min_comp)
pls.eval <- data.frame(obs = ds$train[,1], pred = chem_plsPredict_train[,1,1])
defaultSummary(pls.eval)

##      RMSE  Rsquared       MAE 
## 1.7105077 0.1358886 1.3526818

Results from the test set:

chem_plsPredict <- predict(chem_plsFit, data = ds$test, ncomp = min_comp)
pls.eval <- data.frame(obs = ds$test[,1], pred = chem_plsPredict[,1,1])
defaultSummary(pls.eval)

##       RMSE   Rsquared        MAE 
## 1.91966600 0.00446297 1.59148120

While our RMSE is within range of both the training and test data set the Rsquared is significantly less for the test set suggesting a model that does not explain much of the variation.
This may not be the most optimal model.

e.

dt <- setDT(varImp(chem_plsFit), keep.rownames = TRUE)[]
options(scipen=99)
data.frame(dt[order(-Overall)])

##                                    rn           Overall
## 1    BiologicalMaterial04.boxcoxtrans  1.31112243939398
## 2              ManufacturingProcess29  0.77341701670702
## 3              ManufacturingProcess45  0.72543078751800
## 4  ManufacturingProcess03.boxcoxtrans  0.41837888147600
## 5                BiologicalMaterial08  0.35690900111193
## 6                BiologicalMaterial12  0.34046191064591
## 7              ManufacturingProcess40  0.32954182222043
## 8              ManufacturingProcess36  0.31434324127400
## 9              ManufacturingProcess30  0.29523141240910
## 10             ManufacturingProcess41  0.26167185391625
## 11             ManufacturingProcess11  0.26116974648025
## 12             ManufacturingProcess21  0.24975733031797
## 13               BiologicalMaterial01  0.23742911190879
## 14             ManufacturingProcess43  0.21228437013414
## 15             ManufacturingProcess17  0.18610069292116
## 16               BiologicalMaterial04  0.17652195449247
## 17               BiologicalMaterial07  0.14858631410337
## 18             ManufacturingProcess02  0.14094019383270
## 19             ManufacturingProcess03  0.12583526118383
## 20             ManufacturingProcess10  0.12152373003526
## 21               BiologicalMaterial06  0.11718234452644
## 22             ManufacturingProcess39  0.11043096613233
## 23             ManufacturingProcess37  0.11010899668671
## 24             ManufacturingProcess13  0.09909003285226
## 25             ManufacturingProcess09  0.09647785649708
## 26             ManufacturingProcess34  0.08982642389392
## 27             ManufacturingProcess38  0.08833845706754
## 28             ManufacturingProcess01  0.08398218556732
## 29             ManufacturingProcess33  0.08218889899008
## 30               BiologicalMaterial05  0.06380298847161
## 31             ManufacturingProcess31  0.05061694573618
## 32             ManufacturingProcess32  0.04650901471839
## 33               BiologicalMaterial02  0.03563216276840
## 34               BiologicalMaterial10  0.03218525642229
## 35             ManufacturingProcess23  0.03028469255186
## 36             ManufacturingProcess42  0.02917992823264
## 37             ManufacturingProcess22  0.02449152377729
## 38               BiologicalMaterial11  0.02393291327013
## 39               BiologicalMaterial03  0.01754929077114
## 40             ManufacturingProcess28  0.01243580803101
## 41             ManufacturingProcess06  0.01100159118131
## 42             ManufacturingProcess27  0.01011411117224
## 43 ManufacturingProcess02.boxcoxtrans  0.00904868059662
## 44 ManufacturingProcess01.boxcoxtrans  0.00820789665616
## 45             ManufacturingProcess35  0.00780177810024
## 46             ManufacturingProcess25  0.00765268202301
## 47             ManufacturingProcess18  0.00758657771746
## 48             ManufacturingProcess24  0.00606977998651
## 49             ManufacturingProcess19  0.00597429840038
## 50             ManufacturingProcess26  0.00523155371168
## 51             ManufacturingProcess04  0.00460861463640
## 52             ManufacturingProcess14  0.00390491278230
## 53             ManufacturingProcess20  0.00351396194673
## 54             ManufacturingProcess16  0.00194251903057
## 55             ManufacturingProcess05  0.00189788227258
## 56 ManufacturingProcess06.boxcoxtrans  0.00016216059574
## 57             ManufacturingProcess15  0.00009526443582
## 58 ManufacturingProcess18.boxcoxtrans  0.00000191740596
## 59 ManufacturingProcess16.boxcoxtrans  0.00000087852302
## 60             ManufacturingProcess08 -0.00000008228568
## 61 ManufacturingProcess07.boxcoxtrans -0.00000360323736
## 62 ManufacturingProcess08.boxcoxtrans -0.00001459629022
## 63 ManufacturingProcess12.boxcoxtrans -0.00002205129960
## 64             ManufacturingProcess12 -0.00004543694794
## 65             ManufacturingProcess44 -0.00728754501795
## 66               BiologicalMaterial09 -0.11186450241894
## 67             ManufacturingProcess07 -0.11347906787007

The Manufacturing Processes dominate the list for variable importance with a few biological materials having importance.

f.

The top manufacturing process is ManufacturingProcess40

ds$org_ds %>% 
  select(Yield, ManufacturingProcess40) %>% 
  cor() %>% 
  .[2]

## [1] -0.04735445

ManufacturingProcess40 represents a process that can decrease the Yield due to negative correlation ManufacturingProcess40 and Yield
Not sure if ManufacturingProcess40 decreases Yield but the negative relationship is there
The top biological process is BiologicalMaterial12

ds$org_ds %>% 
  select(Yield, BiologicalMaterial12) %>% 
  cor() %>% 
  .[2]

## [1] 0.3674976

We notice the opposite for the Biological indicator, there’s a positive correlation to Yield.
Not sure if Yield improvement is because of BiologicalMaterial12 but a positive relationship is noticed

Data 624 - Homework 7

Ohannes (Hovig) Ohannessian

4/5/2019

a.

b.

c.

d.

e.

f.

a.

b.

c.

d.

e.

f.