DATA 624 HW 7
Christopher Ayre
April 10, 2022
suppressPackageStartupMessages(library(AppliedPredictiveModeling))
suppressPackageStartupMessages(library(caret))
suppressPackageStartupMessages(library(Amelia))
suppressPackageStartupMessages(library(MASS))
suppressPackageStartupMessages(library(RANN))
suppressPackageStartupMessages(library(pls))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(elasticnet))
suppressPackageStartupMessages(library(glmnet))
suppressPackageStartupMessages(library(psych))
suppressPackageStartupMessages(library(glmnet))Home Work 7
Question 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)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?
nzv <- nearZeroVar(fingerprints)
not_nzv <- fingerprints[, -nzv]
ncol(not_nzv)## [1] 388388 predictors are left for modeling
(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 re-sampled estimate of R2?
set.seed(456)
split <- permeability %>%
createDataPartition(p = 0.8, list = FALSE, times = 1)
# Separate test and training data
#fingerprints train
Xtrain.data <- not_nzv[split, ]
#fingerprints test
xtest.data <- not_nzv[-split, ]
#permeability train
Ytrain.data <- permeability[split, ]
#permeability test
ytest.data <- permeability[-split, ]
# PLS Model
ctrl <- trainControl(method = "cv", number = 10)
plsmod <- train(x = Xtrain.data, y = Ytrain.data, method = "pls", tuneLength = 20, trControl = ctrl, preProc = c("center", "scale"))
plsmod## Partial Least Squares
##
## 133 samples
## 388 predictors
##
## Pre-processing: centered (388), scaled (388)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 119, 121, 119, 120, 121, 121, ...
## Resampling results across tuning parameters:
##
## ncomp RMSE Rsquared MAE
## 1 12.59309 0.3188302 9.598363
## 2 11.36710 0.4595255 7.749247
## 3 11.07388 0.4746904 8.057209
## 4 11.16464 0.4785656 8.382748
## 5 11.13931 0.4855600 8.043630
## 6 11.02535 0.5135675 7.981960
## 7 11.06765 0.5073874 8.127469
## 8 10.98450 0.5229158 8.197463
## 9 11.23015 0.5068022 8.431096
## 10 11.49945 0.4834073 8.591727
## 11 11.83014 0.4725379 8.710300
## 12 11.89082 0.4662821 8.809819
## 13 12.04049 0.4616730 8.933141
## 14 12.42299 0.4413192 9.361306
## 15 12.56859 0.4383558 9.479632
## 16 12.96051 0.4211730 9.635792
## 17 13.27607 0.4045419 9.846898
## 18 13.53374 0.4034482 10.047704
## 19 13.57828 0.4068468 10.111798
## 20 13.62695 0.4053035 10.088053
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was ncomp = 8.plot(plsmod)
The optimal number of principal components included in the PLS model is 8
(D)
Predict the response for the test set. What is the test set estimate of R2?
predictions <- plsmod %>% predict(xtest.data)
cbind(
RMSE = RMSE(predictions, ytest.data),
R_squared = caret::R2(predictions, ytest.data)
)## RMSE R_squared
## [1,] 11.84771 0.5312741plot(predictions, col = "green", main = "Observed (Permeability) vs. Predicted", xlab = "", ylab = "Predictions")
par(new = TRUE)
plot(ytest.data, col = "blue", axes=F, ylab = "", xlab="Observed")
abline(0, 1, col='red')
(E)
Try building other models discussed in this chapter. Do any have better predictive performance?
glm.model <- glmnet(Xtrain.data, Ytrain.data, family="gaussian", alpha=0.5, lambda=0.001)
permeability.glm.predict <- predict(glm.model, xtest.data)
plot(permeability.glm.predict, ytest.data, 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(ytest.data, permeability.glm.predict)^2, 2)))
text(0, 25, paste("RMSE = ", round(sqrt(sum((ytest.data - permeability.glm.predict)^2)), 2)))
(F)
Would you recommend any of your models to replace the permeability laboratory experiment?
I would not recommend any other model to replace the permeability laboratory experiment. The PLS model works best on this data.
Question 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).
na.counts = as.data.frame(((sapply(ChemicalManufacturingProcess,
function(x) sum(is.na(x))))/ nrow(ChemicalManufacturingProcess))*100)
names(na.counts) = "counts"
na.counts = cbind(variables = rownames(na.counts), data.frame(na.counts, row.names = NULL))
na.counts[na.counts$counts > 0,] %>% arrange(counts) %>% mutate(name = factor(variables, levels = variables)) %>%
ggplot(aes(x = name, y = counts)) + geom_segment( aes(xend = name, yend = 0)) +
geom_point(size = 4, color = "steelblue2") + coord_flip() + theme_bw() +
labs(title = "Proportion of Missing Data", x = "Variables", y = "% of Missing data") +
scale_y_continuous(labels = scales::percent_format(scale = 1))
pre.process = preProcess(ChemicalManufacturingProcess[, -c(1)], method = "knnImpute")
chemical = predict(pre.process, ChemicalManufacturingProcess[, -c(1)])(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?
corr = cor(chemical)
corr[corr == 1] = NA
corr[abs(corr) < 0.85] = NA
corr = na.omit(reshape::melt(corr))
head(corr[order(-abs(corr$value)),], 10)## X1 X2 value
## 2090 ManufacturingProcess26 ManufacturingProcess25 0.9975339
## 2146 ManufacturingProcess25 ManufacturingProcess26 0.9975339
## 2148 ManufacturingProcess27 ManufacturingProcess26 0.9960721
## 2204 ManufacturingProcess26 ManufacturingProcess27 0.9960721
## 2091 ManufacturingProcess27 ManufacturingProcess25 0.9934932
## 2203 ManufacturingProcess25 ManufacturingProcess27 0.9934932
## 1685 ManufacturingProcess20 ManufacturingProcess18 0.9917474
## 1797 ManufacturingProcess18 ManufacturingProcess20 0.9917474
## 2095 ManufacturingProcess31 ManufacturingProcess25 0.9706780
## 2431 ManufacturingProcess25 ManufacturingProcess31 0.9706780tooHigh = findCorrelation(cor(chemical), 0.90)
chemical = chemical[, -tooHigh]
(pre.process = preProcess(chemical, method = c("YeoJohnson", "center", "scale")))## Created from 176 samples and 47 variables
##
## Pre-processing:
## - centered (47)
## - ignored (0)
## - scaled (47)
## - Yeo-Johnson transformation (41)
##
## Lambda estimates for Yeo-Johnson transformation:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.8404 0.7085 0.8795 0.9698 1.2484 2.7992set.seed(525)
intrain = createDataPartition(ChemicalManufacturingProcess$Yield, p = 0.8, list = FALSE)
train.p = chemical[intrain, ]
train.r = ChemicalManufacturingProcess$Yield[intrain]
test.p = chemical[-intrain, ]
test.r = ChemicalManufacturingProcess$Yield[-intrain]
# Elastic Net Regression
elastic.model = train(x = train.p, y = train.r, method = "glmnet",
trControl = trainControl("cv", number = 10),
tuneLength = 10, metric = "Rsquared")ggplot(elastic.model) + labs(title = "R-squared Error of Elastic Model") + theme(legend.position = "top")## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 10. Consider
## specifying shapes manually if you must have them.## Warning: Removed 40 rows containing missing values (geom_point).
data.frame(rsquared = elastic.model[["results"]][["Rsquared"]][as.numeric(rownames(elastic.model$bestTune))],
rmse = elastic.model[["results"]][["RMSE"]][as.numeric(rownames(elastic.model$bestTune))])## rsquared rmse
## 1 0.6287452 1.123853(D)
Predict the response for the test set. What is the value of the performance metric and how does this compare with the re-sampled performance metric on the training set?
prediction = predict(elastic.model, test.p)
xyplot(test.r ~ prediction, type = c("p", "g"),
main = "Predicted vs Observed",
xlab = "Predicted",
ylab = "Observed",
panel = function(x, y) {
panel.xyplot(x, y)
panel.abline(lm(y ~ x))
})
postResample(prediction, test.r)## RMSE Rsquared MAE
## 0.9383793 0.7499969 0.8097391rbind(actual = describe(ChemicalManufacturingProcess$Yield)[,-c(1,6,7)],
prediction = describe(prediction)[,-c(1,6,7)])## n mean sd median min max range skew kurtosis se
## actual 176 40.18 1.85 39.97 35.25 46.34 11.09 0.31 -0.11 0.14
## prediction 32 40.23 1.64 40.18 36.28 44.39 8.11 0.07 0.14 0.29(E)
Which predictors are most important in the model you have trained? Do either the biological or process predictors dominate the list?
The top 20 most important variables are ranked below. Manufacturing Processes dominate the list.
varImp(elastic.model)## glmnet variable importance
##
## only 20 most important variables shown (out of 47)
##
## Overall
## ManufacturingProcess32 100.000
## ManufacturingProcess09 53.321
## ManufacturingProcess13 31.462
## ManufacturingProcess17 23.242
## ManufacturingProcess06 21.093
## ManufacturingProcess37 20.358
## ManufacturingProcess34 19.434
## BiologicalMaterial05 16.262
## ManufacturingProcess07 14.693
## ManufacturingProcess39 12.715
## BiologicalMaterial06 12.287
## ManufacturingProcess04 10.545
## ManufacturingProcess15 10.541
## ManufacturingProcess45 8.553
## BiologicalMaterial03 8.374
## ManufacturingProcess28 5.871
## ManufacturingProcess35 5.623
## ManufacturingProcess12 4.488
## ManufacturingProcess43 3.731
## ManufacturingProcess01 2.991ManufacturingProcess32 is the most important variable out of 47, because it has the largest coefficient.
coeffs = coef(elastic.model$finalModel, elastic.model$bestTune$lambda)
(var = data.frame(cbind(variables = coeffs@Dimnames[[1]][coeffs@i+1], coef = coeffs@x)))## variables coef
## 1 (Intercept) 40.1970608100663
## 2 BiologicalMaterial03 0.0704416803906508
## 3 BiologicalMaterial05 0.136787219379557
## 4 BiologicalMaterial06 0.103348569136921
## 5 BiologicalMaterial07 -0.00171902488229934
## 6 ManufacturingProcess01 0.0251590704175465
## 7 ManufacturingProcess04 0.0887033333312616
## 8 ManufacturingProcess06 0.177423181152678
## 9 ManufacturingProcess07 -0.123589777547673
## 10 ManufacturingProcess08 -0.00494927194017869
## 11 ManufacturingProcess09 0.448507043730052
## 12 ManufacturingProcess12 0.0377525554422428
## 13 ManufacturingProcess13 -0.264639735934071
## 14 ManufacturingProcess15 0.0886662868461521
## 15 ManufacturingProcess17 -0.19550032258035
## 16 ManufacturingProcess28 -0.0493841629391022
## 17 ManufacturingProcess32 0.841151058535199
## 18 ManufacturingProcess34 0.163467256666831
## 19 ManufacturingProcess35 -0.0472995610106374
## 20 ManufacturingProcess37 -0.171240962591606
## 21 ManufacturingProcess39 0.106951946343177
## 22 ManufacturingProcess43 0.0313807369041294
## 23 ManufacturingProcess45 0.071945124575606(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?
cor(ChemicalManufacturingProcess$Yield,
ChemicalManufacturingProcess$ManufacturingProcess32)## [1] 0.6083321cor(ChemicalManufacturingProcess$Yield,
ChemicalManufacturingProcess$ManufacturingProcess13)## [1] -0.5036797ManufacturingProcess32 has the highest positive correlaion. Interesting, ManufacturingProcess13 has a negative correlation.
The team should focus on improving predictors ManufacturingProcess32 and lowering predictor ManufacturingProcess13.