Both exercises are from Chapter 6 of Kuhn and Johnson, Applied Predictive Modeling. Exercise 6.2 predicts the permeability of drug compounds from molecular fingerprints, and Exercise 6.3 predicts the yield of a pharmaceutical manufacturing process. The models used come from the chapter: ordinary linear regression, partial least squares (PLS), ridge regression and the elastic net (which becomes the lasso when its ridge penalty is zero).

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) Load the data

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

data(permeability)
dim(fingerprints)
## [1]  165 1107
length(permeability)
## [1] 165

There are 165 compounds. fingerprints holds the 1,107 binary predictors and permeability holds the response. Permeability is the quantity a lab experiment measures, so a model that predicts it well could stand in for some of those experiments.

hist(permeability, breaks = 20, col = "lightblue",
     main = "Permeability", xlab = "Permeability")

summary(permeability)
##   permeability  
##  Min.   : 0.06  
##  1st Qu.: 1.55  
##  Median : 4.91  
##  Mean   :12.24  
##  3rd Qu.:15.47  
##  Max.   :55.60

Most compounds have low permeability and a few have much higher values (the mean is well above the median), so the response is right-skewed. That matters later, because a few large values can dominate the squared-error measures.

(b) Remove low-frequency predictors

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?

# nearZeroVar() returns the column numbers of predictors that are almost constant
nzv_cols <- nearZeroVar(fingerprints)
length(nzv_cols)
## [1] 719
fp_filtered <- fingerprints[, -nzv_cols]
ncol(fp_filtered)
## [1] 388

Of the 1,107 fingerprints, 719 are dropped and 388 remain. A substructure that appears in only a handful of the 165 molecules gives a model almost nothing to learn from, and in a resample it can easily be absent from the whole training portion.

(c) Split the data and tune a PLS model

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?

# createDataPartition() samples within quantile groups of the response, so the
# training and test sets have similar permeability distributions
set.seed(624)
train_rows <- createDataPartition(permeability, p = 0.80, list = FALSE)

perm_train_x <- fp_filtered[train_rows, ]
perm_test_x  <- fp_filtered[-train_rows, ]
perm_train_y <- permeability[train_rows]
perm_test_y  <- permeability[-train_rows]

c(train = length(perm_train_y), test = length(perm_test_y))
## train  test 
##   133    32

The fingerprints are already 0/1, so the only pre-processing is centering and scaling, which PLS needs when predictors are on different scales (Section 6.3). train() applies it inside every resample. Resampling is 10-fold cross-validation repeated 5 times; with only 133 training compounds, a single 10-fold run would be noisy.

perm_ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 5)

set.seed(100)
perm_pls <- train(perm_train_x, perm_train_y,
                  method = "pls",
                  tuneLength = 20,          # try 1 to 20 components
                  trControl = perm_ctrl,
                  preProc = c("center", "scale"))
perm_pls
## Partial Least Squares 
## 
## 133 samples
## 388 predictors
## 
## Pre-processing: centered (388), scaled (388) 
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 120, 120, 119, 121, 120, 119, ... 
## Resampling results across tuning parameters:
## 
##   ncomp  RMSE      Rsquared   MAE      
##    1     13.54431  0.2999465  10.357891
##    2     12.45699  0.4353791   9.034741
##    3     12.55645  0.4408147   9.625860
##    4     12.54530  0.4485399   9.500076
##    5     12.37792  0.4704753   9.163230
##    6     12.15460  0.4880841   8.979011
##    7     12.03877  0.4844514   9.069455
##    8     11.96407  0.4770859   9.183494
##    9     11.97064  0.4823157   9.174287
##   10     12.22528  0.4726105   9.278042
##   11     12.18000  0.4824664   9.256237
##   12     12.09841  0.4894173   9.184731
##   13     12.37982  0.4756196   9.287083
##   14     12.67083  0.4565530   9.444517
##   15     12.93843  0.4406317   9.575829
##   16     13.17834  0.4295027   9.770078
##   17     13.53809  0.4131543   9.985511
##   18     13.68082  0.4041703  10.102511
##   19     13.62904  0.4065538  10.030040
##   20     13.57007  0.4065453  10.012104
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was ncomp = 8.
plot(perm_pls, main = "PLS: resampled RMSE by number of components")

perm_pls$bestTune
##   ncomp
## 8     8
perm_pls$results[perm_pls$bestTune$ncomp, ]
##   ncomp     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
## 8     8 11.96407 0.4770859 9.183494 2.663436  0.1899359 2.079471

The resampled RMSE is lowest with 8 latent variables, and the corresponding resampled R^2 is about 0.48. The curve drops from one to two components, stays roughly level up to about 12, and then rises as more components start fitting noise, so a modest number of directions carries the usable signal.

(d) Test set performance for PLS

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

perm_pls_pred <- predict(perm_pls, perm_test_x)
postResample(pred = perm_pls_pred, obs = perm_test_y)
##       RMSE   Rsquared        MAE 
## 10.1771878  0.5563284  8.0478163
# scale of the test responses, to judge the RMSE against
c(median = median(perm_test_y), sd = sd(perm_test_y))
##   median       sd 
##  5.23750 15.21124

postResample() returns the test set RMSE, R^2 and MAE. The test R^2 of about 0.56 is a little higher than the resampled R^2. The test set has only 32 compounds, so one number from it is noisy, and the resampled estimate is the steadier of the two.

(e) Other models from the chapter

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

Ordinary linear regression is the natural first try. With more predictors than training compounds, the least squares solution is not unique (Section 6.2), which lm() shows directly:

perm_lm_data <- data.frame(perm_train_x, permeability = perm_train_y)
perm_lm <- lm(permeability ~ ., data = perm_lm_data)
c(coefficients = length(coef(perm_lm)), not_estimable = sum(is.na(coef(perm_lm))))
##  coefficients not_estimable 
##           389           293

The number of coefficients that lm() could not estimate is shown above. R fits the largest identifiable model by dropping predictors, so predictions still exist but the coefficients cannot be interpreted and the fit is unstable (Section 6.2). Ridge regression and the elastic net handle this case by penalizing the coefficients. Both are tuned with the same resampling settings and the same seed, so the resamples line up.

# lambda = 0 would be ordinary least squares, which fails here, so the grid starts above zero
ridge_grid <- data.frame(lambda = seq(0.05, 1, length = 20))

set.seed(100)
perm_ridge <- train(perm_train_x, perm_train_y,
                    method = "ridge",
                    tuneGrid = ridge_grid,
                    trControl = perm_ctrl,
                    preProc = c("center", "scale"))
perm_ridge$bestTune
##   lambda
## 3   0.15
# lambda is the ridge penalty; fraction is the position along the lasso path (1 = no lasso penalty, the ridge penalty still applies)
enet_grid <- expand.grid(lambda = c(0.01, 0.1, 1),
                         fraction = seq(0.05, 1, length = 20))

set.seed(100)
perm_enet <- train(perm_train_x, perm_train_y,
                   method = "enet",
                   tuneGrid = enet_grid,
                   trControl = perm_ctrl,
                   preProc = c("center", "scale"))
perm_enet$bestTune
##    fraction lambda
## 29     0.45    0.1
perm_models <- list(PLS = perm_pls, Ridge = perm_ridge, Enet = perm_enet)

# resampled performance at each model's best tuning values: the row of $results with the lowest RMSE
perm_resampled <- t(sapply(perm_models, function(m) {
  best <- m$results[which.min(m$results$RMSE), ]
  c(RMSE = best$RMSE, Rsquared = best$Rsquared)
}))

# test set performance
perm_test <- t(sapply(perm_models, function(m)
  postResample(pred = predict(m, perm_test_x), obs = perm_test_y)))

perm_summary <- data.frame(
  resampled_RMSE = perm_resampled[, "RMSE"],
  resampled_R2   = perm_resampled[, "Rsquared"],
  test_RMSE      = perm_test[, "RMSE"],
  test_R2        = perm_test[, "Rsquared"])
round(perm_summary, 3)
##       resampled_RMSE resampled_R2 test_RMSE test_R2
## PLS           11.964        0.477    10.177   0.556
## Ridge         12.312        0.494    11.664   0.485
## Enet          11.667        0.502    10.928   0.503
perm_resamp <- resamples(perm_models)
bwplot(perm_resamp, metric = "RMSE")

# the models share resamples (same seed), so the differences are paired
summary(diff(perm_resamp))
## 
## Call:
## summary.diff.resamples(object = diff(perm_resamp))
## 
## p-value adjustment: bonferroni 
## Upper diagonal: estimates of the difference
## Lower diagonal: p-value for H0: difference = 0
## 
## MAE 
##       PLS       Ridge     Enet 
## PLS             0.111     0.583
## Ridge 1.0000000           0.472
## Enet  0.0009463 0.0062288      
## 
## RMSE 
##       PLS      Ridge    Enet   
## PLS            -0.3484   0.2968
## Ridge 0.168282           0.6452
## Enet  0.226825 0.006649        
## 
## Rsquared 
##       PLS    Ridge     Enet     
## PLS          -0.016505 -0.024990
## Ridge 0.2169           -0.008485
## Enet  0.1272 1.0000

The three penalized and projection models land close together. The elastic net has the lowest resampled RMSE and the highest resampled R^2, PLS has the best test set numbers, and ridge is last on the test set. The models share the same resamples, so diff() compares them pair by pair, resample by resample; the box plot shows the spread of the resampled RMSE for each model, and the table of p-values is printed after it (upper diagonal: differences, lower diagonal: adjusted p-values). For RMSE the elastic net is distinguishable from ridge (small adjusted p-value), while there is no evidence of a difference between the elastic net and PLS or between ridge and PLS in this comparison. Apart from that one comparison, the ranking in the table above is a ranking of point estimates only, and because the repeated folds overlap these comparisons are exploratory. Ordinary linear regression was not carried forward: with more predictors than compounds its coefficients are not unique, so it cannot be interpreted, and the penalized and projection models are the chapter’s answer to that.

(f) Would any model replace the permeability experiment?

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

I would not recommend replacing the permeability experiment with any of these models. The best test R^2 is about 0.56 (postResample() reports this as the squared correlation between predicted and observed values), and the test RMSE of about 10 is only modestly below the standard deviation of the test responses (about 15, printed in part (d)), so the typical error is still large compared with how much permeability varies across compounds. The models do carry some signal, so their predictions might help prioritize which candidate molecules go to the lab first. With errors this large, a low-scoring compound should not be discarded automatically, and this analysis does not validate any ranking or threshold. Using the predictions in place of the measurement, where a wrong value could send a poor drug candidate forward or discard a good one, would need a much more accurate model than 133 training compounds produced, and no acceptable error tolerance is given for this exercise, so this is a judgment of caution and not a fixed R^2 cutoff.

Exercise 6.3

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) Load the data

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.

data(ChemicalManufacturingProcess)
dim(ChemicalManufacturingProcess)
## [1] 176  58

The exercise text refers to a matrix processPredictors and a vector yield. In the installed package the same data come as one data frame, ChemicalManufacturingProcess, whose first column Yield is the response and whose other 57 columns are the 12 biological and 45 process predictors for the 176 runs. Splitting it back into the two pieces:

process_predictors <- ChemicalManufacturingProcess[, -1]
yield <- ChemicalManufacturingProcess$Yield
summary(yield)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   35.25   38.75   39.97   40.18   41.48   46.34

(b) Impute the missing values

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

missing_per_predictor <- colSums(is.na(process_predictors))
missing_per_predictor[missing_per_predictor > 0]
## ManufacturingProcess01 ManufacturingProcess02 ManufacturingProcess03 
##                      1                      3                     15 
## ManufacturingProcess04 ManufacturingProcess05 ManufacturingProcess06 
##                      1                      1                      2 
## ManufacturingProcess07 ManufacturingProcess08 ManufacturingProcess10 
##                      1                      1                      9 
## ManufacturingProcess11 ManufacturingProcess12 ManufacturingProcess14 
##                     10                      1                      1 
## ManufacturingProcess22 ManufacturingProcess23 ManufacturingProcess24 
##                      1                      1                      1 
## ManufacturingProcess25 ManufacturingProcess26 ManufacturingProcess27 
##                      5                      5                      5 
## ManufacturingProcess28 ManufacturingProcess29 ManufacturingProcess30 
##                      5                      5                      5 
## ManufacturingProcess31 ManufacturingProcess33 ManufacturingProcess34 
##                      5                      5                      5 
## ManufacturingProcess35 ManufacturingProcess36 ManufacturingProcess40 
##                      5                      5                      1 
## ManufacturingProcess41 
##                      1
c(missing_cells = sum(is.na(process_predictors)),
  percent_of_cells = round(100 * mean(is.na(process_predictors)), 2),
  runs_with_a_missing_value = sum(rowSums(is.na(process_predictors)) > 0))
##             missing_cells          percent_of_cells runs_with_a_missing_value 
##                    106.00                      1.06                     24.00

The missing cells are concentrated in some of the process predictors (ManufacturingProcess03 is missing the most), as Section 3.4 describes, and none of the biological predictors are missing.

An imputation model is built from the data, so it has to be built from the training data only (Section 3.4). The split therefore comes before the imputation here, and the test rows are imputed with the model learned from the training rows.

set.seed(624)
cmp_rows <- createDataPartition(yield, p = 0.80, list = FALSE)

cmp_train_x <- process_predictors[cmp_rows, ]
cmp_test_x  <- process_predictors[-cmp_rows, ]
cmp_train_y <- yield[cmp_rows]
cmp_test_y  <- yield[-cmp_rows]

c(train = length(cmp_train_y), test = length(cmp_test_y))
## train  test 
##   144    32

Before imputing, any near-zero-variance predictor is dropped (Section 3.5):

cmp_nzv <- nearZeroVar(cmp_train_x, names = TRUE)
cmp_nzv
## [1] "BiologicalMaterial07"
cmp_train_x <- cmp_train_x[, !(names(cmp_train_x) %in% cmp_nzv)]
cmp_test_x  <- cmp_test_x[, !(names(cmp_test_x) %in% cmp_nzv)]

preProcess() with method = "knnImpute" replaces each missing value with the average of that predictor over the nearest training runs (Sections 3.4 and 3.8). It standardizes the predictors to compute the distances, so the imputed data come back centered and scaled.

cmp_impute <- preProcess(cmp_train_x, method = "knnImpute")

cmp_train_imp <- predict(cmp_impute, cmp_train_x)
cmp_test_imp  <- predict(cmp_impute, cmp_test_x)

c(train_missing = sum(is.na(cmp_train_imp)), test_missing = sum(is.na(cmp_test_imp)))
## train_missing  test_missing 
##             0             0

No missing values remain in either set. This explicit imputation answers the exercise and shows that the training-based model fills the test rows too. For tuning, however, the imputation is repeated inside each resample (next part), because imputing once on all 144 training rows would let the held-out fold influence the values its own predictors receive.

(c) Tune a model

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?

Three of the chapter’s models are tuned with the same resampling, and the one with the lowest resampled RMSE is used for the rest of the exercise. RMSE is the performance metric train() minimizes for a numeric response. The models are tuned on the predictors that still contain missing values, with preProc = "knnImpute" so that the imputation (and the centering and scaling it includes) is fitted inside each resample, as the book describes for train() (Section 3.8). predict() on the test set then applies the training imputation automatically.

cmp_ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 5)

set.seed(100)
cmp_pls <- train(cmp_train_x, cmp_train_y, method = "pls",
                 tuneLength = 20, trControl = cmp_ctrl,
                 preProc = "knnImpute")

set.seed(100)
cmp_ridge <- train(cmp_train_x, cmp_train_y, method = "ridge",
                   tuneGrid = data.frame(lambda = seq(0.01, 0.5, length = 20)),
                   trControl = cmp_ctrl, preProc = "knnImpute")

set.seed(100)
cmp_enet <- train(cmp_train_x, cmp_train_y, method = "enet",
                  tuneGrid = expand.grid(lambda = c(0.01, 0.1, 1),
                                         fraction = seq(0.05, 1, length = 20)),
                  trControl = cmp_ctrl, preProc = "knnImpute")

cmp_models <- list(PLS = cmp_pls, Ridge = cmp_ridge, Enet = cmp_enet)
lapply(cmp_models, function(m) m$bestTune)
## $PLS
##   ncomp
## 4     4
## 
## $Ridge
##       lambda
## 13 0.3194737
## 
## $Enet
##    fraction lambda
## 26      0.3    0.1
cmp_resampled <- t(sapply(cmp_models, function(m) {
  best <- m$results[which.min(m$results$RMSE), ]
  c(RMSE = best$RMSE, Rsquared = best$Rsquared, MAE = best$MAE)
}))
round(cmp_resampled, 3)
##        RMSE Rsquared   MAE
## PLS   1.347    0.532 1.042
## Ridge 1.524    0.529 1.103
## Enet  1.196    0.576 0.969
cmp_resamp <- resamples(cmp_models)
bwplot(cmp_resamp, metric = "RMSE")

summary(diff(cmp_resamp))
## 
## Call:
## summary.diff.resamples(object = diff(cmp_resamp))
## 
## p-value adjustment: bonferroni 
## Upper diagonal: estimates of the difference
## Lower diagonal: p-value for H0: difference = 0
## 
## MAE 
##       PLS      Ridge    Enet    
## PLS            -0.06116  0.07252
## Ridge 0.009262           0.13368
## Enet  0.064502 0.009122         
## 
## RMSE 
##       PLS      Ridge    Enet   
## PLS            -0.1773   0.1509
## Ridge 0.006838           0.3282
## Enet  0.180302 0.024873        
## 
## Rsquared 
##       PLS    Ridge    Enet    
## PLS           0.00343 -0.04358
## Ridge 1.0000          -0.04701
## Enet  0.1511 0.1601

The elastic net has the lowest resampled RMSE, so it is the model used for the rest of the exercise. Its tuning values are in the bestTune output above (ridge penalty 0.1 and path fraction 0.3). The diff() table shows the elastic net is distinguishable from ridge (small adjusted p-value for RMSE), with no evidence of a difference from PLS (large p-value), so the choice between the elastic net and PLS rests on the point estimates alone.

# resampled performance of the chosen model at its best tuning values
cmp_enet$results[which.min(cmp_enet$results$RMSE), ]
##    lambda fraction     RMSE  Rsquared      MAE    RMSESD RsquaredSD     MAESD
## 26    0.1      0.3 1.196128 0.5758582 0.969269 0.3216711  0.1808719 0.2087112

The optimal resampled RMSE of the elastic net is 1.196, with the resample-to-resample spread in RMSESD.

(d) Test set performance

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?

cmp_test_perf <- t(sapply(cmp_models, function(m)
  postResample(pred = predict(m, cmp_test_x), obs = cmp_test_y)))
round(cmp_test_perf, 3)
##        RMSE Rsquared   MAE
## PLS   1.596    0.500 1.077
## Ridge 2.304    0.374 1.249
## Enet  1.220    0.724 0.936

The elastic net’s test set RMSE is 1.22, close to its resampled RMSE of 1.196, so the resampling estimate predicted the test performance well here. It is also the best of the three models on the test set. PLS has a higher test RMSE than resampled RMSE, and ridge is worst on both. The test set has only 32 runs, so the agreement for the elastic net is encouraging but one small sample.

(e) Most important predictors

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

The chosen model is the elastic net, so its importance comes from its own coefficients. The predictors were centered and scaled inside train(), so the coefficients are on a common scale and their absolute sizes can be compared, as the chapter does when it extracts coefficients with predict(..., type = "coefficients", mode = "fraction") (Section 6.5). The coefficients are taken at the selected path fraction of the final model.

cmp_coef <- predict(cmp_enet$finalModel,
                    s = cmp_enet$bestTune$fraction,
                    mode = "fraction",
                    type = "coefficients")$coefficients

# the lasso part of the penalty sets many coefficients exactly to zero
c(nonzero = sum(cmp_coef != 0), zero = sum(cmp_coef == 0))
## nonzero    zero 
##      13      43
cmp_rank <- sort(abs(cmp_coef[cmp_coef != 0]), decreasing = TRUE)
importance_table <- data.frame(
  Predictor = names(cmp_rank),
  Coefficient = round(cmp_coef[names(cmp_rank)], 3),
  Absolute_coefficient = round(as.numeric(cmp_rank), 3),
  row.names = NULL)
importance_table
##                 Predictor Coefficient Absolute_coefficient
## 1  ManufacturingProcess32       0.651                0.651
## 2  ManufacturingProcess09       0.415                0.415
## 3  ManufacturingProcess36      -0.161                0.161
## 4  ManufacturingProcess17      -0.134                0.134
## 5  ManufacturingProcess34       0.133                0.133
## 6    BiologicalMaterial06       0.113                0.113
## 7  ManufacturingProcess11       0.107                0.107
## 8  ManufacturingProcess06       0.100                0.100
## 9  ManufacturingProcess13      -0.095                0.095
## 10   BiologicalMaterial03       0.074                0.074
## 11 ManufacturingProcess37      -0.067                0.067
## 12 ManufacturingProcess44       0.062                0.062
## 13 ManufacturingProcess39       0.028                0.028
top_names <- names(cmp_rank)
table(ifelse(grepl("^Biological", top_names), "Biological", "Process"))
## 
## Biological    Process 
##          2         11
dotchart(rev(importance_table$Absolute_coefficient),
         labels = rev(importance_table$Predictor),
         pch = 19, xlab = "Absolute coefficient (standardized predictors)",
         main = "Elastic net: predictors with nonzero coefficients")

The lasso part of the penalty has set 43 of the 56 coefficients to exactly zero, so the model uses 13 predictors: 11 process predictors and 2 biological ones (BiologicalMaterial06 and 03). Process predictors therefore dominate the retained set and take five of the six largest coefficients, led by ManufacturingProcess32 and 09, which are well ahead of the rest. This is a statement about this fitted model, not about biological and process variables in general, and the starting pool also holds more process predictors (45) than biological ones (11 after filtering). The ranking is by coefficient size on standardized predictors, so it shows how strongly the model leans on each predictor given the others. It is not a causal effect and not a variable-importance-in-projection score. Process predictors also matter more in practice because they can be adjusted in the plant, while the biological ones describe the raw material and can only be used to judge a batch before processing starts.

(f) Relationships between the top predictors and yield

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?

# original (un-imputed, un-scaled) values, so the axes are in the process units
top6 <- head(top_names, 6)
par(mfrow = c(2, 3))
for (v in top6) {
  plot(ChemicalManufacturingProcess[[v]], ChemicalManufacturingProcess$Yield,
       xlab = v, ylab = "Yield", main = v, pch = 19, col = "steelblue")
  abline(lm(Yield ~ ChemicalManufacturingProcess[[v]], data = ChemicalManufacturingProcess),
         col = "firebrick", lwd = 2)
}

par(mfrow = c(1, 1))
round(sapply(top6, function(v)
  cor(ChemicalManufacturingProcess[[v]], ChemicalManufacturingProcess$Yield,
      use = "complete.obs")), 3)
## ManufacturingProcess32 ManufacturingProcess09 ManufacturingProcess36 
##                  0.608                  0.503                 -0.528 
## ManufacturingProcess17 ManufacturingProcess34   BiologicalMaterial06 
##                 -0.426                  0.175                  0.478
round(cor(ChemicalManufacturingProcess[, top6], use = "pairwise.complete.obs"), 2)
##                        ManufacturingProcess32 ManufacturingProcess09
## ManufacturingProcess32                   1.00                   0.04
## ManufacturingProcess09                   0.04                   1.00
## ManufacturingProcess36                  -0.79                  -0.05
## ManufacturingProcess17                   0.02                  -0.72
## ManufacturingProcess34                   0.07                   0.08
## BiologicalMaterial06                     0.60                   0.23
##                        ManufacturingProcess36 ManufacturingProcess17
## ManufacturingProcess32                  -0.79                   0.02
## ManufacturingProcess09                  -0.05                  -0.72
## ManufacturingProcess36                   1.00                  -0.01
## ManufacturingProcess17                  -0.01                   1.00
## ManufacturingProcess34                  -0.03                  -0.10
## BiologicalMaterial06                    -0.53                   0.01
##                        ManufacturingProcess34 BiologicalMaterial06
## ManufacturingProcess32                   0.07                 0.60
## ManufacturingProcess09                   0.08                 0.23
## ManufacturingProcess36                  -0.03                -0.53
## ManufacturingProcess17                  -0.10                 0.01
## ManufacturingProcess34                   1.00                -0.06
## BiologicalMaterial06                    -0.06                 1.00

These plots and correlations use the full data set after the test evaluation. They are descriptive exploration of the supplied runs and were not fed back into tuning. Runs with a missing value for a plotted predictor are left out of that plot, so each plot uses the available observations from the 176 runs (171 for ManufacturingProcess36 and 34).

Among the top six, ManufacturingProcess32, 09, 34 and BiologicalMaterial06 have positive correlations with yield in the table above, and ManufacturingProcess36 and 17 have negative ones. The scatter plots show loose trends, not tight lines, so each is a tendency across the available runs and does not guarantee that moving one setting changes yield. ManufacturingProcess34 shows why the correlation and the model coefficient answer different questions: its correlation with yield is only 0.175, yet it has the fifth-largest coefficient, because the model holds the other predictors fixed. It also takes only four distinct values, so its plot is a set of vertical stripes. ManufacturingProcess36 shows the same coarse stripes with six values.

The correlation matrix shows that some of the top predictors move together. ManufacturingProcess32 and 36 are strongly negatively correlated (-0.79), 09 and 17 are too (-0.72), and BiologicalMaterial06 is correlated with both 32 (0.60) and 36 (-0.53). The individual slopes therefore cannot be read as the separate effect of each setting, and raising 32 may be partly the same signal as lowering 36. Correlation and regularization can affect how the model distributes weight among predictors. In this fitted model ManufacturingProcess32 has a larger absolute coefficient than 36 (0.651 against 0.161), so the two should not be interpreted as independent contributions to yield.

For improving yield, the process predictors are the ones that can be acted on. The useful next step is to test the plant at higher settings of ManufacturingProcess32 and 09 and lower settings of 17 and 36, within ranges the process can run safely, in a designed experiment, since these data only show association. The biological predictors cannot be changed, but BiologicalMaterial06 can be used to flag raw material batches that are likely to give low yield, so that the process settings can be adjusted for them. With a yield gain of 1% worth about one hundred thousand dollars per batch, even a small confirmed improvement from one of these settings would be worth the experiment.