Rpubs URL: https://rpubs.com/umaisabdullah/Data624_Homework7

Objective

Complete and submit Exercises 6.2 and 6.3 from Applied Predictive Modeling by Kuhn and Johnson. We have to submit both our Rpubs link as well as attach the .pdf file with our code.

Overview

The earlier homeworks were about looking at data and cleaning it up. This homework is the first one where I actually build models and measure how well they predict.

Both exercises use linear models, but they are not ordinary least squares. The reason is that both datasets have a problem that plain regression cannot handle. In Exercise 6.2 there are 1,107 predictors and only 165 samples, so there are far more predictors than rows and ordinary regression has no unique answer at all. In Exercise 6.3 the predictors are heavily correlated with each other, which makes ordinary regression unstable.

The methods in this chapter fix those problems in two different ways. Partial least squares builds a small number of new predictors out of the originals, choosing combinations that relate to the response. Ridge and the elastic net keep all the predictors but shrink their coefficients toward zero, which trades a little bias for a lot less variance.

The other theme running through this homework is honest measurement. Every model gets tuned by cross validation on the training set and then tested once on data it has never seen. The gap between those two numbers is what tells me whether a model will really work.

Pre-Requisite

Both datasets come from the AppliedPredictiveModeling package, so nothing needs to be downloaded into the working directory. The packages used are:

  1. AppliedPredictiveModeling contains the permeability and chemical manufacturing datasets.
  2. caret provides nearZeroVar(), the data splitting, the cross validation and the model tuning.
  3. pls provides the partial least squares method.
  4. elasticnet provides the ridge and elastic net methods.
  5. ggplot2, tidyr and dplyr are used for the plots and data reshaping.

Exercise 6.2

Question. 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:

library(AppliedPredictiveModeling)
data(permeability)

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

My thought process. The shape of this data decides everything. There are 1,107 predictors and only 165 compounds, so there are nearly seven times more predictors than samples. Ordinary least squares cannot even be fitted here, because with more predictors than rows there is no single best answer. It would fit the training data perfectly and predict nothing.

That is why the chapter’s methods matter. PLS compresses the predictors down into a handful of components, and the penalized methods shrink the coefficients. Both are ways of handling far too many predictors.

(a) Load the data

data(permeability)

dim(fingerprints)
## [1]  165 1107
dim(permeability)
## [1] 165   1
summary(as.numeric(permeability))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.06    1.55    4.91   12.24   15.47   55.60
data.frame(permeability = as.numeric(permeability)) |>
  ggplot(aes(x = permeability)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  labs(title = "Distribution of the permeability response",
       x = "Permeability", y = "Count")

The response is strongly right skewed. The median is about 4.9 but the mean is about 12.2, and the maximum reaches 55.6. When the mean sits well above the median like that, it means a few large values are pulling the average up. Most compounds have low permeability and a few have much higher values.

This matters for how I read the results later. With a skewed response, the error measures get dominated by how well the model does on those few large values.

(b) Filter low frequency predictors

My thought process. The fingerprints are binary, so each one says whether a molecule has a particular substructure. If a substructure appears in only two or three of the 165 compounds, that column is almost entirely zeros and carries almost no usable information. That is exactly what nearZeroVar() finds.

nzv_cols <- nearZeroVar(fingerprints)

length(nzv_cols)
## [1] 719
fingerprints_filtered <- fingerprints[, -nzv_cols]

ncol(fingerprints_filtered)
## [1] 388

How many predictors are left? 388. The filter removes 719 of the original 1,107 predictors, leaving 388 for modeling.

That is a large reduction, and it is worth noting that the problem is not solved by it. Even at 388 predictors there are still more than twice as many predictors as the 165 compounds, so I still cannot use ordinary regression. The filter makes the problem smaller but not easy.

(c) Split, pre-process and tune a PLS model

My thought process. I split before doing anything else so the test set stays genuinely untouched. Centring and scaling go inside the train() call rather than being applied to the whole dataset first, which matters. If I scaled everything up front, information from the test set would leak into the training process through the means and standard deviations. Putting it inside train() means the scaling is recomputed within each resampling fold.

I use 10 fold cross validation to choose the number of components, since with only 125 training rows a single validation split would be too noisy to trust.

set.seed(624)

perm_y <- as.numeric(permeability)

train_idx <- createDataPartition(perm_y, p = 0.75, list = FALSE)

X_train <- fingerprints_filtered[train_idx, ]
X_test  <- fingerprints_filtered[-train_idx, ]
y_train <- perm_y[train_idx]
y_test  <- perm_y[-train_idx]

c(train_rows = nrow(X_train), test_rows = nrow(X_test))
## train_rows  test_rows 
##        125         40
ctrl <- trainControl(method = "cv", number = 10)

set.seed(624)
pls_model <- train(
  x = X_train, y = y_train,
  method = "pls",
  tuneLength = 20,
  trControl = ctrl,
  preProcess = c("center", "scale")
)

pls_model$bestTune
##   ncomp
## 8     8
# The resampled performance at the chosen number of components.
pls_model$results |>
  filter(ncomp == pls_model$bestTune$ncomp) |>
  select(ncomp, RMSE, Rsquared, MAE) |>
  mutate(across(where(is.numeric), \(x) round(x, 4)))
##   ncomp    RMSE Rsquared    MAE
## 1     8 10.9996   0.5294 8.7679
ggplot(pls_model$results, aes(x = ncomp, y = RMSE)) +
  geom_line(colour = "steelblue") +
  geom_point(colour = "steelblue") +
  geom_vline(xintercept = pls_model$bestTune$ncomp,
             linetype = "dashed", colour = "darkred") +
  labs(title = "Cross validated RMSE against number of PLS components",
       subtitle = "Dashed line marks the chosen value",
       x = "Number of components", y = "RMSE (cross validated)")

How many latent variables are optimal? Eight components.

What is the resampled R2? About 0.53, with a cross validated RMSE near 11.

The profile plot shows why eight is chosen. The error falls steeply for the first few components, flattens out, and then starts to drift upward as more components are added. That upward drift at the right hand side is overfitting beginning to show.

The compression here is striking. The model takes 388 predictors and reduces them to 8 components, which is a huge reduction, and those 8 components explain about half the variation in permeability.

(d) Predict the test set

pls_pred <- predict(pls_model, X_test)

pls_test_perf <- postResample(pred = pls_pred, obs = y_test)

round(pls_test_perf, 4)
##     RMSE Rsquared      MAE 
##  12.4751   0.4089   9.6282

What is the test set R2? About 0.41.

That is lower than the resampled estimate of 0.53, and the drop is the honest and expected result. Cross validation was used to choose the number of components, so it is slightly optimistic by construction. The test set was never used for any decision, so its number is the trustworthy one.

A gap of this size is also not surprising with a test set of only 40 compounds. With so few rows, a handful of hard to predict molecules can move the score noticeably.

(e) Other models

My thought process. PLS handles the wide data by building components. The penalized methods take a different approach, keeping all the predictors but shrinking their coefficients. Comparing the two is the point of this part, so I tune both with the same folds and the same pre-processing to make the comparison fair.

set.seed(624)
ridge_model <- train(
  x = X_train, y = y_train,
  method = "ridge",
  tuneGrid = data.frame(lambda = seq(0.01, 0.5, length = 15)),
  trControl = ctrl,
  preProcess = c("center", "scale")
)

ridge_model$bestTune
##   lambda
## 5   0.15
set.seed(624)
enet_model <- train(
  x = X_train, y = y_train,
  method = "enet",
  tuneGrid = expand.grid(lambda   = c(0.001, 0.01, 0.1),
                         fraction = seq(0.05, 1, length = 10)),
  trControl = ctrl,
  preProcess = c("center", "scale")
)

enet_model$bestTune
##     fraction lambda
## 12 0.1555556   0.01

Note that I start the ridge lambda grid at 0.01 rather than 0. A lambda of exactly zero is ordinary least squares, which cannot be fitted when there are more predictors than rows, so including it causes some folds to fail.

get_cv <- function(model, label) {
  best <- model$results[which.min(model$results$RMSE), ]
  data.frame(Model = label,
             CV_RMSE = round(best$RMSE, 3),
             CV_Rsquared = round(best$Rsquared, 4))
}

get_test <- function(model, label) {
  p <- postResample(predict(model, X_test), y_test)
  data.frame(Model = label,
             Test_RMSE = round(p[["RMSE"]], 3),
             Test_Rsquared = round(p[["Rsquared"]], 4))
}

comparison <- bind_rows(
  get_cv(pls_model,   "PLS"),
  get_cv(ridge_model, "Ridge"),
  get_cv(enet_model,  "Elastic net")
) |>
  left_join(
    bind_rows(
      get_test(pls_model,   "PLS"),
      get_test(ridge_model, "Ridge"),
      get_test(enet_model,  "Elastic net")
    ),
    by = "Model"
  )

comparison
##         Model CV_RMSE CV_Rsquared Test_RMSE Test_Rsquared
## 1         PLS  11.000      0.5294    12.475        0.4089
## 2       Ridge  11.603      0.5338    12.666        0.4021
## 3 Elastic net   9.883      0.5993    12.871        0.3315

Do any have better predictive performance? No, and the way the elastic net fails is the most useful lesson in this exercise.

On the cross validated training results the elastic net looks clearly best, with an RMSE of about 9.88 and an R2 of about 0.60, against roughly 11.0 and 0.53 for PLS. On the test set the order reverses. The elastic net is the worst of the three, with an R2 of only 0.33 and the highest RMSE (12.87). PLS is best, at an R2 of 0.41, with ridge just behind it at 0.40.

The drop from cross validation to test is what separates them. PLS falls from 0.53 to 0.41 and ridge from 0.53 to 0.40, similar and modest drops. The elastic net falls from 0.60 to 0.33, roughly twice as far. Its cross validation score promised the most and delivered the least. I think two things explain that.

  • It was tuned over more combinations. The elastic net grid has 30 candidate settings, against 20 for PLS and 15 for ridge. Each cross validated score is itself a noisy estimate, and picking the best of many noisy estimates rewards luck as well as skill: the more settings you try, the more optimistic the winner’s score becomes.
  • The chosen model is sparse. With a fraction of about 0.16, the elastic net keeps only a small subset of the fingerprints. With 125 training compounds, exactly which subset gets picked is unstable, so the model latches onto substructures that happen to matter in the training folds but do not carry over to new compounds. PLS and ridge spread their weight across all 388 predictors, which makes them more stable.

My choice is PLS. It has the best test performance, and its cross validated estimate was the most honest guide to how it would actually do. Ridge is essentially tied with it. The gap between PLS and ridge on the test set is far smaller than the uncertainty that comes with only 40 test compounds, so I would not claim PLS is genuinely better than ridge, only that it is no worse and is the simpler of the two to describe (8 components). The elastic net’s gap is a different matter: it is large, and it shows up as a fall from best to worst.

What all three agree on is the ceiling. Every method explains somewhere between a third and two fifths of the variation on new compounds, and none breaks past that. When three different approaches stop at about the same place, the limit is in the data rather than in the choice of model.

(f) Would I recommend replacing the laboratory experiment?

No, not as a replacement. The test set R2 is about 0.41, meaning the model leaves close to 60 percent of the variation in permeability unexplained. The test RMSE of roughly 12.5 also needs to be read against the response itself, which has a median of only 4.9. So the typical error is larger than the typical value being predicted. A lab deciding whether a compound is promising could not rely on that.

But I would recommend it as a screening filter, which is a different job. The value here is not in replacing the experiment but in ordering the queue. If a company has thousands of candidate compounds and can only test a limited number in the lab, a model that captures even 40 percent of the variation can push the more promising compounds to the front. Some will be wrong, but the lab work still gets done on a better ordered list than random choice.

This is the practical distinction I would put to the company. Replacing the assay requires accuracy the model does not have. Prioritising which assays to run first only requires the model to be better than guessing, and it clears that bar.

Two things would have to be checked before using it even that way. The compounds here are all from one project, so the model may not carry over to chemically different molecules. And the skewed response means the model is weakest on exactly the high permeability compounds that are most worth finding.

What I did not try. The obvious next step for that last problem would be to model log permeability rather than permeability itself. A log transform pulls in the long right tail, so the few very permeable compounds no longer dominate the squared errors, and the model is fitted on a scale where errors are proportional rather than absolute. I kept the original scale so that the R2 and RMSE answer the question as asked, in the units the lab measures, but I would expect a log model to rank compounds more reliably, which is exactly what a screening filter needs. Principal components regression is the other model from this chapter I did not fit; it compresses the predictors like PLS but ignores the response when building components, so I would not expect it to beat PLS here.

Exercise 6.3

Question. 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, 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 percent will boost revenue by approximately one hundred thousand dollars per batch:

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.

My thought process. This exercise has a different shape to 6.2 and a different goal. Here there are 176 runs and only 57 predictors, so I have more rows than columns and ordinary regression could technically be fitted.

The real point is in parts e and f, and the question sets it up in the sentence about biological versus process predictors. Biological predictors describe the raw material and cannot be changed. Process predictors are settings that can be adjusted. So if the important predictors turn out to be process ones, the company can act on them. If they are biological, the company can only screen incoming material. That distinction is worth real money here, since the question says one percent of yield is worth about a hundred thousand dollars a batch.

(a) Load the data

data(ChemicalManufacturingProcess)

cmp <- ChemicalManufacturingProcess

dim(cmp)
## [1] 176  58
# The first column is the response.
names(cmp)[1]
## [1] "Yield"
predictors_raw <- cmp[, -1]
yield          <- cmp[, 1]

c(predictors    = ncol(predictors_raw),
  biological    = sum(grepl("^Biological",    names(predictors_raw))),
  manufacturing = sum(grepl("^Manufacturing", names(predictors_raw))))
##    predictors    biological manufacturing 
##            57            12            45

This matches the question: 57 predictors made up of 12 biological and 45 manufacturing process measurements, across 176 runs.

(b) Impute the missing values

My thought process. Before imputing I want to know how much is missing and where, because that decides whether imputation is reasonable at all. This is the same check I did on the Soybean data in Homework 4, where the missingness turned out to be structural and imputation would have been the wrong answer. Here I expect something different, but it is worth confirming rather than assuming.

c(missing_cells   = sum(is.na(predictors_raw)),
  percent_missing = round(100 * mean(is.na(predictors_raw)), 2),
  rows_affected   = sum(!complete.cases(predictors_raw)),
  cols_affected   = sum(colSums(is.na(predictors_raw)) > 0))
##   missing_cells percent_missing   rows_affected   cols_affected 
##          106.00            1.06           24.00           28.00
colSums(is.na(predictors_raw)) |>
  (\(x) x[x > 0])() |>
  sort(decreasing = TRUE) |>
  head(10)
## ManufacturingProcess03 ManufacturingProcess11 ManufacturingProcess10 
##                     15                     10                      9 
## ManufacturingProcess25 ManufacturingProcess26 ManufacturingProcess27 
##                      5                      5                      5 
## ManufacturingProcess28 ManufacturingProcess29 ManufacturingProcess30 
##                      5                      5                      5 
## ManufacturingProcess31 
##                      5

Only about 1 percent of the values are missing, spread across 28 predictors and 24 rows. The worst single predictor is missing 15 values out of 176.

This is a very different situation from the Soybean data. There the gaps were concentrated in particular classes and were clearly structural. Here the missing values are scattered thinly across many predictors and rows, with no single column badly affected. That is the pattern where imputation makes sense, because there are plenty of complete neighbouring rows to estimate from.

Split first, then impute. In Exercise 6.2 I made the point that any pre-processing fitted on the whole dataset leaks information from the test set into training. Imputation is pre-processing too, and knnImpute also centres and scales the data, so it computes means, standard deviations and nearest neighbours. If I fitted it on all 176 runs, the test runs would help fill the training gaps and set the scaling. So I make the train/test split here, fit the imputation on the training runs only, and then apply that fitted imputation to both sets. The test runs are filled using training neighbours and scaled with training statistics, exactly as a genuinely new run would be.

set.seed(624)
chem_idx <- createDataPartition(yield, p = 0.75, list = FALSE)

c(train_rows = length(chem_idx), test_rows = nrow(predictors_raw) - length(chem_idx))
## train_rows  test_rows 
##        132         44
# Fit the imputation (and its centring and scaling) on the training runs only.
impute_model <- preProcess(predictors_raw[chem_idx, ], method = "knnImpute")

# Apply it to every run. Test runs are filled from training neighbours only.
predictors_imputed <- predict(impute_model, predictors_raw)

sum(is.na(predictors_imputed))
## [1] 0

I use k nearest neighbours imputation, which fills each gap using the most similar runs rather than just putting the column average everywhere. Since the manufacturing predictors are related to each other, similar runs are genuinely informative about what the missing value would have been.

# Near-zero variance is also judged on the training runs only.
nzv_chem <- nearZeroVar(predictors_imputed[chem_idx, ])

length(nzv_chem)
## [1] 1
if (length(nzv_chem) > 0) {
  predictors_imputed <- predictors_imputed[, -nzv_chem]
}

ncol(predictors_imputed)
## [1] 56

1 predictor is nearly constant in the training runs and is dropped, leaving 56.

(c) Split, pre-process and tune a model

My thought process. I fit two models rather than one. The question says to tune a model of my choice, but PLS and the elastic net handle correlated predictors in different ways, and fitting both lets me check whether the conclusion depends on the choice. Doing that turned out to matter here, as the next part shows.

chem_X_train <- predictors_imputed[chem_idx, ]
chem_X_test  <- predictors_imputed[-chem_idx, ]
chem_y_train <- yield[chem_idx]
chem_y_test  <- yield[-chem_idx]

c(train_rows = nrow(chem_X_train), test_rows = nrow(chem_X_test))
## train_rows  test_rows 
##        132         44
set.seed(624)
chem_pls <- train(
  x = chem_X_train, y = chem_y_train,
  method = "pls",
  tuneLength = 15,
  trControl = ctrl,
  preProcess = c("center", "scale")
)

chem_pls$bestTune
##   ncomp
## 1     1
chem_pls$results |>
  select(ncomp, RMSE, Rsquared) |>
  head(8) |>
  mutate(across(where(is.numeric), \(x) round(x, 4)))
##   ncomp   RMSE Rsquared
## 1     1 1.3101   0.4722
## 2     2 1.3793   0.5214
## 3     3 1.3507   0.5562
## 4     4 1.5909   0.5518
## 5     5 1.8760   0.5541
## 6     6 1.9388   0.5582
## 7     7 2.0085   0.5575
## 8     8 2.3051   0.5533
pls_res    <- chem_pls$results
pls_best_n <- chem_pls$bestTune$ncomp
pls_r2_n   <- pls_res$ncomp[which.max(pls_res$Rsquared)]
pls_cv     <- pls_res[pls_res$ncomp == pls_best_n, ]
pls_rmse_8 <- max(pls_res$RMSE[pls_res$ncomp <= 8])

Something worth stopping on. Cross validation chooses 1 component by RMSE. R2, though, is highest at 6 components (0.56 against 0.47 at the chosen value), so the two measures disagree about how many components to use. The more telling pattern is in the RMSE column. It is 1.31 at the chosen value but climbs to 2.31 within the first eight components, while R2 barely moves. RMSE squares each error, so a jump like that while R2 holds steady usually means a few runs in some folds are getting wildly wrong predictions, typically runs whose predictor values are unusual, while most runs are predicted about as well as before. That is a sign the PLS fit is unstable, so I do not accept it at face value, which is why fitting a second model matters.

set.seed(624)
chem_enet <- train(
  x = chem_X_train, y = chem_y_train,
  method = "enet",
  tuneGrid = expand.grid(lambda   = c(0.001, 0.01, 0.1),
                         fraction = seq(0.05, 1, length = 10)),
  trControl = ctrl,
  preProcess = c("center", "scale")
)

chem_enet$bestTune
##     fraction lambda
## 12 0.1555556   0.01
enet_cv <- chem_enet$results[which.min(chem_enet$results$RMSE), ]

bind_rows(
  pls_cv[, c("RMSE", "Rsquared")]  |> mutate(Model = "PLS"),
  enet_cv[, c("RMSE", "Rsquared")] |> mutate(Model = "Elastic net")
) |>
  select(Model, CV_RMSE = RMSE, CV_Rsquared = Rsquared) |>
  mutate(across(where(is.numeric), \(x) round(x, 4)))
##         Model CV_RMSE CV_Rsquared
## 1         PLS  1.3101      0.4722
## 2 Elastic net  1.1110      0.6091
enet_better_cv <- enet_cv$RMSE < pls_cv$RMSE

What is the optimal value of the performance metric? I use RMSE, since it is in the same units as yield. The elastic net is the better model. Its cross validated RMSE is about 1.11 with an R2 of about 0.61, against 1.31 and 0.47 for PLS. So my answer to part (c) is the elastic net, with an optimal resampled RMSE of about 1.11. Yield is measured in percent, so an RMSE near 1.1 is an average error of a little over one percentage point.

(d) Predict the test set

pls_test  <- postResample(predict(chem_pls,  chem_X_test), chem_y_test)
enet_test <- postResample(predict(chem_enet, chem_X_test), chem_y_test)

chem_results <- bind_rows(
  data.frame(Model = "PLS",         t(round(pls_test, 4))),
  data.frame(Model = "Elastic net", t(round(enet_test, 4)))
)

chem_results
##         Model   RMSE Rsquared    MAE
## 1         PLS 1.8245   0.3282 1.3827
## 2 Elastic net 1.3460   0.5862 1.1302

How does the test set compare with the resampled estimate? The elastic net’s cross validated RMSE of about 1.11 becomes 1.35 on the test set, and its R2 goes from about 0.61 to 0.59. PLS goes from an RMSE of 1.31 to 1.82, and from an R2 of 0.47 to 0.33. The elastic net holds up far better. Some drop from the resampled figure is normal, because that figure was used for tuning, but PLS deteriorates much more than the elastic net does, which confirms the instability the tuning table warned about.

The lesson I take is that a single resampled number is not enough on its own. Fitting a second model, and reading the whole tuning table rather than only the chosen row, is what let me judge how far to trust each result before touching the test set.

data.frame(
  observed  = chem_y_test,
  predicted = predict(chem_enet, chem_X_test)
) |>
  ggplot(aes(x = observed, y = predicted)) +
  geom_point(alpha = 0.7, colour = "steelblue") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "darkred") +
  labs(title = "Elastic net: predicted against observed yield on the test set",
       subtitle = "Dashed line is perfect prediction",
       x = "Observed yield (percent)", y = "Predicted yield (percent)")

The points follow the diagonal reasonably closely, which is what a working model should look like. The spread is wider at the extremes, so the model is least reliable for the highest and lowest yielding runs.

(e) Which predictors are most important?

How I measure importance. The obvious tool is caret’s varImp(), but it has to be used with care here.

# Does caret have a model-specific importance method for the elastic net?
enet_no_varimp <- is.null(chem_enet$modelInfo$varImp)
enet_no_varimp
## [1] TRUE

It does not. For a model without its own importance method, varImp() silently falls back to a model-free filter: it fits a loess curve between each predictor and yield, one predictor at a time, and ranks them by that fit. That ranking ignores the elastic net entirely, so it cannot tell me what the model relies on. Instead I rank the predictors by the size of their elastic net coefficients. Every predictor was centred and scaled before fitting, so the coefficient sizes are directly comparable, and the sign of each coefficient gives the direction of its effect with the other predictors held fixed.

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

n_kept <- sum(enet_coef != 0)
c(predictors = length(enet_coef), kept_by_model = n_kept)
##    predictors kept_by_model 
##            56            13
coef_df <- data.frame(Predictor = names(enet_coef), Coefficient = enet_coef) |>
  filter(Coefficient != 0) |>
  mutate(Importance = round(100 * abs(Coefficient) / max(abs(Coefficient)), 1),
         Coefficient = round(Coefficient, 4),
         Type = if_else(grepl("^Biological", Predictor), "Biological", "Process")) |>
  arrange(desc(Importance))

importance_df <- head(coef_df, 15)
rownames(importance_df) <- NULL
importance_df
##                 Predictor Coefficient Importance       Type
## 1  ManufacturingProcess32      0.8332      100.0    Process
## 2  ManufacturingProcess09      0.4355       52.3    Process
## 3  ManufacturingProcess17     -0.1723       20.7    Process
## 4  ManufacturingProcess06      0.1486       17.8    Process
## 5    BiologicalMaterial06      0.1372       16.5 Biological
## 6  ManufacturingProcess34      0.0835       10.0    Process
## 7  ManufacturingProcess37     -0.0697        8.4    Process
## 8  ManufacturingProcess39      0.0438        5.3    Process
## 9    BiologicalMaterial03      0.0350        4.2 Biological
## 10 ManufacturingProcess13     -0.0227        2.7    Process
## 11 ManufacturingProcess11      0.0012        0.1    Process
## 12 ManufacturingProcess43      0.0011        0.1    Process
## 13 ManufacturingProcess45      0.0012        0.1    Process
importance_df |>
  ggplot(aes(x = reorder(Predictor, Importance), y = Importance, fill = Type)) +
  geom_col() +
  coord_flip() +
  labs(title = "Most important predictors in the elastic net",
       subtitle = "Importance = absolute coefficient, scaled so the largest is 100",
       x = NULL, y = "Importance") +
  scale_fill_manual(values = c(Biological = "darkseagreen4", Process = "steelblue"))

table(importance_df$Type)
## 
## Biological    Process 
##          2         11
n_top     <- nrow(importance_df)
n_proc    <- sum(importance_df$Type == "Process")
n_bio     <- sum(importance_df$Type == "Biological")
bio_total <- sum(grepl("^Biological", names(chem_X_train)))
bio_share_data <- bio_total / ncol(chem_X_train)

The elastic net keeps 13 of the 56 predictors and sets the rest exactly to zero, so it is a fairly sparse model: it has decided most predictors add nothing once the others are known.

Do either dominate? Yes, the process predictors do. Of the top 13, 11 are manufacturing process measurements and 2 are biological. The single most important predictor is ManufacturingProcess32, a process predictor.

It is worth reading the counts against what is available. Biological predictors make up 20 percent of the predictors in the data but 15 percent of this list, so they appear no more often than their share of the data would suggest. The process predictors lead, and they are the ones that can actually be changed, which is the more useful outcome for the company.

(f) Relationships between the top predictors and yield

My thought process. Knowing which predictors matter is not enough to act on. The company needs to know the direction of each relationship, because raising a setting that is negatively related to yield would make things worse. I compare two views of direction: the simple correlation of each predictor with yield on its own, and the sign of its coefficient in the model, which is its effect with the other predictors held fixed.

top_predictors <- head(importance_df$Predictor, 8)

correlations <- data.frame(
  Predictor   = top_predictors,
  Coefficient = importance_df$Coefficient[match(top_predictors, importance_df$Predictor)],
  Correlation = round(sapply(top_predictors,
                             \(v) cor(predictors_imputed[[v]], yield)), 3)
) |>
  mutate(Type  = if_else(grepl("^Biological", Predictor), "Biological", "Process"),
         Agree = sign(Coefficient) == sign(Correlation))
rownames(correlations) <- NULL

correlations
##                Predictor Coefficient Correlation       Type Agree
## 1 ManufacturingProcess32      0.8332       0.608    Process  TRUE
## 2 ManufacturingProcess09      0.4355       0.503    Process  TRUE
## 3 ManufacturingProcess17     -0.1723      -0.426    Process  TRUE
## 4 ManufacturingProcess06      0.1486       0.392    Process  TRUE
## 5   BiologicalMaterial06      0.1372       0.478 Biological  TRUE
## 6 ManufacturingProcess34      0.0835       0.177    Process  TRUE
## 7 ManufacturingProcess37     -0.0697      -0.159    Process  TRUE
## 8 ManufacturingProcess39      0.0438       0.036    Process  TRUE
predictors_imputed |>
  select(all_of(top_predictors)) |>
  mutate(Yield = yield) |>
  pivot_longer(-Yield, names_to = "Predictor", values_to = "Value") |>
  ggplot(aes(x = Value, y = Yield)) +
  geom_point(alpha = 0.5, colour = "steelblue") +
  geom_smooth(method = "lm", se = FALSE, colour = "darkred") +
  facet_wrap(~ Predictor, scales = "free_x", ncol = 3) +
  labs(title = "Top predictors against yield",
       subtitle = "Values are centred and scaled by the imputation step",
       x = "Predictor value (scaled)", y = "Yield (percent)")

How this could help improve yield. The practical value comes from separating the two kinds of predictors, which is what the question was pointing at from the start.

For the process predictors the company has a direct lever. The model says that raising ManufacturingProcess32, ManufacturingProcess09, ManufacturingProcess06, ManufacturingProcess34 and ManufacturingProcess39, while lowering ManufacturingProcess17 and ManufacturingProcess37 is associated with higher yield, with the other predictors held fixed. Those are settings that can be adjusted, so they are the ones worth testing. Given that one percent of yield is worth about a hundred thousand dollars a batch, even a small reliable gain pays for the work quickly.

For the biological predictors the company cannot change the raw material, but it can measure it on arrival. BiologicalMaterial06 is positively related to yield, so these measurements can be used to grade incoming batches, to predict which will run well, and possibly to negotiate with suppliers on material quality.

For every one of the top eight, the model coefficient and the simple correlation point the same way, so the directions above are consistent whichever view is taken.

An important caution. Everything above is association, not proof of cause. The model shows that these predictors move together with yield, not that changing a setting will cause yield to rise. It is entirely possible that some third factor drives both.

So the right next step is not to change the process on the strength of this model. It is to run controlled experiments on the handful of process settings the model points to. The real value of this analysis is that it narrows 56 predictors down to the 13 the model keeps, and to a short list at the top of that, which makes the experiments far cheaper than testing everything.