library(AppliedPredictiveModeling)
library(caret)Loading required package: ggplot2Loading required package: latticedata(permeability)DATA_624_Homework_7
Exercise 6.2
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:
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_predictors <- nearZeroVar(fingerprints)
fp_without_nzv <- fingerprints[, -nzv_predictors]
ncol(fp_without_nzv)[1] 388After filtering out the predictors that have low frequencies using the nearZeroVar function, 388 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 resampled estimate of R2?
set.seed(8585)
index <- createDataPartition(permeability, p = 0.8, list = FALSE)
training_x <- fp_without_nzv[index, ]
training_y <- permeability[index]
test_x <- fp_without_nzv[-index, ]
test_y <- permeability[-index]
ctrl <- trainControl(method = "cv", number = 10)set.seed(8585)
pls_model <- train(
x = training_x,
y = training_y,
method = "pls",
tuneLength = 20,
preProcess = c("center", "scale"),
trControl = ctrl
)
pls_modelPartial Least Squares
133 samples
388 predictors
Pre-processing: centered (388), scaled (388)
Resampling: Cross-Validated (10 fold)
Summary of sample sizes: 120, 121, 119, 120, 121, 119, ...
Resampling results across tuning parameters:
ncomp RMSE Rsquared MAE
1 12.71019 0.3569097 9.868381
2 11.50798 0.4639481 8.253118
3 12.07417 0.4328394 8.984111
4 12.36622 0.4260573 9.314963
5 12.59525 0.4227659 9.375107
6 12.49468 0.4356304 9.195829
7 12.25332 0.4326802 9.115110
8 12.07990 0.4383105 9.086909
9 11.86706 0.4507181 8.855191
10 11.98456 0.4521816 8.898481
11 12.27041 0.4421552 9.099205
12 12.50620 0.4295990 9.127911
13 12.72302 0.4195657 9.425178
14 12.82747 0.4154657 9.399555
15 13.17340 0.3976287 9.728923
16 13.52831 0.3805939 9.890032
17 13.93516 0.3655227 10.221197
18 14.33440 0.3547180 10.539780
19 14.53166 0.3477287 10.750755
20 14.82472 0.3420551 11.002697
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 2.The optimal number of latent variables is 2 based on the optimal model having the smallest RMSE value. The corresponding resampled estimate of R2 is 0.4639481.
d.
Predict the response for the test set. What is the test set estimate of R2?
pls_predictions <- predict(pls_model, newdata = test_x)
postResample(pred = pls_predictions, obs = test_y) RMSE Rsquared MAE
12.6187375 0.4094358 8.1671427 The test set estimate of R2 is 0.4094358.
e.
Try building other models discussed in this chapter. Do any have better predictive performance?
set.seed(8585)
pcr_model <- train(
x = training_x,
y = training_y,
method = "pcr",
tuneLength = 20,
preProcess = c("center", "scale"),
trControl = ctrl
)
pcr_modelPrincipal Component Analysis
133 samples
388 predictors
Pre-processing: centered (388), scaled (388)
Resampling: Cross-Validated (10 fold)
Summary of sample sizes: 120, 121, 119, 120, 121, 119, ...
Resampling results across tuning parameters:
ncomp RMSE Rsquared MAE
1 14.67179 0.1337500 11.318228
2 14.67644 0.1344230 11.236392
3 12.82825 0.3374569 9.779824
4 12.01225 0.4310405 9.068880
5 11.89805 0.4394744 8.810441
6 11.95468 0.4354151 8.830675
7 12.06392 0.4255400 8.843129
8 11.66538 0.4545882 8.385389
9 11.47850 0.4604221 8.397701
10 11.52933 0.4565857 8.536617
11 11.75517 0.4312481 8.972941
12 11.70165 0.4387568 8.897261
13 11.77262 0.4427149 8.876303
14 11.83351 0.4403159 8.882031
15 12.13638 0.4465662 9.060966
16 12.35997 0.4398544 9.179170
17 12.53606 0.4199298 9.338353
18 12.49385 0.4204304 9.305853
19 12.61676 0.4147533 9.344162
20 12.98894 0.3998746 9.607979
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 9.This PCR model has a slightly lower RMSE value of 11.47850 and slightly lower Rsquared value of 0.4604221 compared to the PLS model with a RMSE value at 11.50798 and Rsquared value of 0.4639481.
pcr_predictions <- predict(pcr_model, newdata = test_x)
postResample(pred = pcr_predictions, obs = test_y) RMSE Rsquared MAE
12.7798132 0.3915923 8.2976581 set.seed(8585)
ridge_model <- train(
x = training_x,
y = training_y,
method = "ridge",
tuneGrid = data.frame(lambda = seq(0.001, 0.2, length = 15)),
preProcess = c("center", "scale"),
trControl = ctrl)
ridge_modelRidge Regression
133 samples
388 predictors
Pre-processing: centered (388), scaled (388)
Resampling: Cross-Validated (10 fold)
Summary of sample sizes: 120, 121, 119, 120, 121, 119, ...
Resampling results across tuning parameters:
lambda RMSE Rsquared MAE
0.00100000 109389.49232 0.2483295 69157.511421
0.01521429 15.34873 0.3418982 11.301686
0.02942857 14.25789 0.3744960 10.584041
0.04364286 13.65256 0.3952091 10.092600
0.05785714 13.39947 0.4053308 9.861195
0.07207143 13.17027 0.4157110 9.669121
0.08628571 13.02623 0.4232481 9.540821
0.10050000 12.92392 0.4289411 9.442574
0.11471429 12.85584 0.4336074 9.380548
0.12892857 12.81057 0.4376775 9.345393
0.14314286 12.79671 0.4401919 9.316883
0.15735714 12.76484 0.4435462 9.286116
0.17157143 12.76485 0.4451390 9.280212
0.18578571 12.74229 0.4487058 9.267647
0.20000000 12.75276 0.4503415 9.273332
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was lambda = 0.1857857.This ridge regression model had higher RMSE values than the PCR and PLS models with the lowest RMSE value being 12.74229 with a Rsquared value of 0.4487058.
ridge_predictions <- predict(ridge_model, newdata = test_x)
postResample(pred = ridge_predictions, obs = test_y) RMSE Rsquared MAE
11.9117467 0.6036854 8.6675539 resamps <- resamples(list(PLS = pls_model, PCR = pcr_model, Ridge = ridge_model))
summary(resamps)
Call:
summary.resamples(object = resamps)
Models: PLS, PCR, Ridge
Number of resamples: 10
MAE
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
PLS 4.791192 6.390647 8.525104 8.253118 9.358236 12.92699 0
PCR 5.455381 6.925619 8.797594 8.397701 9.130175 12.95542 0
Ridge 5.375295 8.333319 8.892763 9.267647 10.414349 13.76323 0
RMSE
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
PLS 5.668109 7.745063 12.25449 11.50798 14.12683 18.28138 0
PCR 6.853576 8.359330 11.94377 11.47850 13.18725 17.81720 0
Ridge 6.427660 11.152088 12.57308 12.74229 14.76140 19.53033 0
Rsquared
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
PLS 0.01740037 0.2513345 0.4346603 0.4639481 0.7533257 0.8990359 0
PCR 0.06219160 0.2654016 0.4306574 0.4604221 0.7263119 0.8718705 0
Ridge 0.01458501 0.3412260 0.4715230 0.4487058 0.6207161 0.8348688 0f.
Would you recommend any of your models to replace the permeability laboratory experiment?
No, I would not recommend replacing the permeability laboratory experiment entirely with any of my models above. This is because based on the results from the models, they don’t capture enough of the variance in the data. I would recommend using the PCR in conjunction with the permeability laboratory experiment because it had the lowest RMSE value.
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), 6.5 Computing 139 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:
library(AppliedPredictiveModeling)
library(caret)
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.2.0 ✔ readr 2.2.0
✔ forcats 1.0.1 ✔ stringr 1.6.0
✔ lubridate 1.9.5 ✔ tibble 3.3.1
✔ purrr 1.2.1 ✔ tidyr 1.3.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ purrr::lift() masks caret::lift()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(RANN)
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.
dim(ChemicalManufacturingProcess)[1] 176 58b.
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)
sum(is.na(ChemicalManufacturingProcess))[1] 106There are 106 missing values. I will use the knnImpute method to fill in these missing values.
impute_model <- preProcess(ChemicalManufacturingProcess, method = "knnImpute")
chem_imputed <- predict(impute_model, ChemicalManufacturingProcess)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?
set.seed(8585)
index <- createDataPartition(chem_imputed$Yield, p = 0.8, list = FALSE)
training_data <- chem_imputed[index, ]
test_data <- chem_imputed[-index, ]
set.seed(8585)
pls_model <- train(
Yield ~ .,
data = training_data,
method = "pls",
tuneLength = 20,
trControl = trainControl(method = "cv", number = 10)
)
pls_modelPartial Least Squares
144 samples
57 predictor
No pre-processing
Resampling: Cross-Validated (10 fold)
Summary of sample sizes: 129, 129, 130, 129, 129, 130, ...
Resampling results across tuning parameters:
ncomp RMSE Rsquared MAE
1 0.7544749 0.4603976 0.6156528
2 0.6710448 0.5819215 0.5548499
3 0.6447270 0.6190149 0.5299569
4 0.6415626 0.6317688 0.5327480
5 0.6519906 0.6227729 0.5378816
6 0.6539245 0.6150491 0.5409481
7 0.6636368 0.6047773 0.5482924
8 0.6771787 0.5966030 0.5684949
9 0.6865190 0.5943795 0.5716314
10 0.7082968 0.5871221 0.5886683
11 0.7276011 0.5756836 0.6004379
12 0.7346708 0.5704047 0.6010584
13 0.7481320 0.5622787 0.6062104
14 0.7634531 0.5519791 0.6125535
15 0.7707148 0.5538808 0.6131337
16 0.7731536 0.5597650 0.6109230
17 0.7740055 0.5625273 0.6062572
18 0.7824350 0.5627300 0.6102496
19 0.7941523 0.5571840 0.6157749
20 0.8079552 0.5540481 0.6206233
RMSE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 4.The optimal value of the training data is ncomp=4 with an RMSE value of 0.6415626 and Rsquared value of 0.6317688.
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?
pls_predictions <- predict(pls_model, newdata = test_data)
postResample(pred = pls_predictions, obs = test_data$Yield) RMSE Rsquared MAE
0.5425112 0.6650965 0.4630009 The performance metric had a RMSE value of 0.5425112, a Rsquared value of 0.6650965, and a MAE value of 0.4630009. This is better than the training set because it has a lower RMSE value and a higher Rsquared value.
e.
Which predictors are most important in the model you have trained? Do either the biological or process predictors dominate the list?
pls_imp <- varImp(pls_model)
Attaching package: 'pls'The following object is masked from 'package:caret':
R2The following object is masked from 'package:stats':
loadingsprint(pls_imp)pls variable importance
only 20 most important variables shown (out of 57)
Overall
ManufacturingProcess32 100.00
ManufacturingProcess13 87.08
ManufacturingProcess17 86.59
ManufacturingProcess36 83.49
ManufacturingProcess09 75.71
BiologicalMaterial02 61.95
ManufacturingProcess33 60.12
BiologicalMaterial06 58.71
ManufacturingProcess06 58.00
BiologicalMaterial03 54.41
ManufacturingProcess11 54.34
ManufacturingProcess12 53.45
BiologicalMaterial08 51.49
BiologicalMaterial12 49.40
BiologicalMaterial11 48.76
BiologicalMaterial04 47.15
BiologicalMaterial01 46.97
ManufacturingProcess28 44.03
ManufacturingProcess04 42.89
ManufacturingProcess37 40.89plot(pls_imp, top = 20)
The manufacturing process predictors dominate the top of the list. The top 5 most important variables are manufacturing process predictors with ManufacturingProcess32 having a score of 100, while the highest-ranking biological predictor (BiologicalMaterial06) only appears 6th on the list with a significantly lower importance score of 61.95.
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?
library(dplyr)
library(corrplot)corrplot 0.95 loaded
Attaching package: 'corrplot'The following object is masked from 'package:pls':
corrplottop_predictors <- rownames(pls_imp$importance)[order(-pls_imp$importance$Overall)][1:10]
cor_matrix <- chem_imputed |>
select(Yield, all_of(top_predictors)) |>
cor(use = "pairwise.complete.obs")
corrplot(cor_matrix)
This information can be useful in improving yield in future runs of the manufacturing process because it shows the relationship between the predictor and the yield. For example, ManufacturingProcess32 has a positive correlation with Yield so increasing that can improve the yield. On the contrary, ManufacturingProcess36, ManufacturingProcess13, and ManufacturingProcess17 have a negative correlation with Yield so decreasing those can help improve the yield.