DATA 624: HOMEWORK 7

Author

Pascal Hermann Kouogang Tafo

Problem 6.2

a) Let’s Load the data

Before loading the AppliedPredictiveModeling library, I installed the package in the R console so that it is not being installed every time that we run our code.

library(AppliedPredictiveModeling)

Warning: package 'AppliedPredictiveModeling' was built under R version 4.5.3

data(permeability)

# let's verify the matrix of the fingerprints and permeability

dim(fingerprints)

[1]  165 1107

length(permeability)

[1] 165

head(permeability)

  permeability
1       12.520
2        1.120
3       19.405
4        1.730
5        1.680
6        0.510

str(permeability)

 num [1:165, 1] 12.52 1.12 19.41 1.73 1.68 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:165] "1" "2" "3" "4" ...
  ..$ : chr "permeability"

The matrix fingerprints contains 165 compounds and 1107 binary molecular predictors for the compounds. The permeability data contains permeability response. Before diving deeper in the analysis our data, let’s generate a general descriptive statistics to have an overall initial interpretation of our data.

summary(permeability[, 1])

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.06    1.55    4.91   12.24   15.47   55.60

This result tells us that we have a data that is right skewed because the mean is way greater than the Median. Moreover, the minimum and maximum values are respectively 0.06 and 55.60.

##b) Filter out predictors with near-zero variance

Let’s load the useful package and check how many predictors are left for modeling

# Load the caret library

library(caret)

Warning: package 'caret' was built under R version 4.5.3

Loading required package: ggplot2

Warning: package 'ggplot2' was built under R version 4.5.3

Loading required package: lattice

Warning: package 'lattice' was built under R version 4.5.2

# Let's apply the nearZeroVar function to the fingerprint predictors

nzv <- nearZeroVar(fingerprints)

# Let's filter out predictors with near zero variance from the fingerprints matrix

fingerprints_filtering <- fingerprints[, -nzv]

# Let's check the number of predictors left for modeling

remaining_pred <- ncol(fingerprints_filtering)
remaining_pred

[1] 388

The result shows us that out of the initial 1107 predictors, there are only 388 predictors remaining for further modeling.

c) Split the data , Pre-process and Tune a PLS(Partial Least Squares) model

Splitting data into training and testing set

We are going to split the 165 compounds data into a 80/20 percent ratio with 80% for training set and 20% for the testing set.

# Let's load the necessary library after installing the package in the console

library(pls)

Warning: package 'pls' was built under R version 4.5.3


Attaching package: 'pls'

The following object is masked from 'package:caret':

    R2

The following object is masked from 'package:stats':

    loadings

# Let's split the data into training and testing set

set.seed(123)

training_index <- createDataPartition(permeability, p = 0.8, list = FALSE)

fingerprints_train <- fingerprints_filtering[training_index, ]
permeability_train <- permeability[training_index]

fingerprints_test  <- fingerprints_filtering[-training_index, ]
permeability_test  <- permeability[-training_index]

Let’s define the 10-fold cross-validation which will be useful when training our partial least squares model.

# configuration of the resampling and tuning protocol in a 50 distinc training and validation iterations

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

# Training the PLS model with 10-fold cross validation

set.seed(123)

pls_model <- train(
  x = fingerprints_train,
  y = permeability_train,
  method = "pls",
  tuneLength = 20,  # evaluating the model using between 1 to 20 latent component
  trControl = ctrl
)

# Let's view the tuning model table result

pls_model

Partial Least Squares 

133 samples
388 predictors

No pre-processing
Resampling: Cross-Validated (10 fold, repeated 5 times) 
Summary of sample sizes: 120, 119, 118, 120, 121, 119, ... 
Resampling results across tuning parameters:

  ncomp  RMSE      Rsquared   MAE      
   1     14.00543  0.2984506  10.761458
   2     11.72971  0.5008207   8.508109
   3     11.49422  0.5244486   8.678784
   4     11.70128  0.5121349   8.964770
   5     11.66094  0.5232784   8.797789
   6     11.55216  0.5301295   8.557093
   7     11.49226  0.5319143   8.598236
   8     11.49660  0.5317352   8.808766
   9     11.70230  0.5263586   9.020912
  10     12.12583  0.5063114   9.215251
  11     12.44899  0.4867215   9.447173
  12     12.49900  0.4754215   9.446703
  13     12.57708  0.4703572   9.441098
  14     12.70370  0.4625174   9.586998
  15     12.89016  0.4520250   9.746643
  16     13.12010  0.4416082   9.873428
  17     13.28203  0.4385247   9.935144
  18     13.36355  0.4349297  10.036145
  19     13.54903  0.4290731  10.233540
  20     13.51941  0.4294793  10.256468

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was ncomp = 7.

Result

The PLS model result shows that the optimal number of latent variables values represented by ncomp is equal to 7. The corresponding resampled estimate of R^2 is 0.530.

Prediction of the test Set response

# let's predict the response of the permeability for the test set

pls_test_predictions <- predict(pls_model, newdata = fingerprints_test)

# Let's find the value of the test set estimate of R^2

pls_test_Rsquare <- cor(pls_test_predictions,permeability_test)^2
pls_test_Rsquare

[1] 0.3497131

We compare the predicted permeability values to the actual permeability test value and we find that the test set estimate of R^2 value is 0.3497.

e) Let’s build other models discussed in this chapter.

I will build the ridge regression, Lasso regression, Elastic net models. It is important to remember that x_train represents fingerprints_train, y_train represents permeability_train above.

pp   <- c("zv", "center", "scale")

# Ridge Regression

set.seed(123)

ridge_fit <- train(x = fingerprints_train, y = permeability_train,
                   method     = "ridge",
                   tuneGrid   = expand.grid(lambda = seq(0.01, 0.4, length = 15)),
                   preProcess = pp,
                   trControl  = ctrl)
ridge_fit$bestTune

     lambda
6 0.1492857

# Elastic Net

library(elasticnet)

Warning: package 'elasticnet' was built under R version 4.5.2

Loading required package: lars

Warning: package 'lars' was built under R version 4.5.2

Loaded lars 1.3

enet_grid <- expand.grid(lambda   = c(0.01, 0.05, 0.1, 0.2),
                         fraction = seq(0.05, 1, length = 20))
set.seed(100)
enet_fit <- train(x = fingerprints_train, y = permeability_train,
                  method     = "enet",
                  tuneGrid   = enet_grid,
                  preProcess = pp,
                  trControl  = ctrl)
enet_fit$bestTune

   fraction lambda
45     0.25    0.1

# Comparison

models <- list(PLS = pls_model, Ridge = ridge_fit,  `Elastic net` = enet_fit)

best_cv <- function(fit) {
  r <- merge(fit$bestTune, fit$results)
  c(CV_RMSE = r$RMSE, CV_Rsq = r$Rsquared)
}

test_perf <- function(fit) {
  p <- postResample(predict(fit, fingerprints_test), permeability_test)
  c(Test_RMSE = unname(p["RMSE"]), Test_Rsq = unname(p["Rsquared"]))
}

comparison <- do.call(rbind, lapply(models, function(m) c(best_cv(m), test_perf(m))))
comparison <- data.frame(
  Model  = names(models),
  Tuning = c(paste("ncomp =", pls_model$bestTune$ncomp),
             paste("lambda =", round(ridge_fit$bestTune$lambda, 3)),
             paste0("lambda = ", enet_fit$bestTune$lambda,
                    ", fraction = ", round(enet_fit$bestTune$fraction, 2))),
  round(comparison, 3), row.names = NULL)
comparison

        Model                        Tuning CV_RMSE CV_Rsq Test_RMSE Test_Rsq
1         PLS                     ncomp = 7  11.492  0.532    11.976    0.350
2       Ridge                lambda = 0.149  12.330  0.505    12.980    0.401
3 Elastic net lambda = 0.1, fraction = 0.25  11.579  0.540    11.490    0.394

Problem 6.3

##a) Loading the data

library(AppliedPredictiveModeling)

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.5.3

Warning: package 'tibble' was built under R version 4.5.2

Warning: package 'tidyr' was built under R version 4.5.2

Warning: package 'readr' was built under R version 4.5.2

Warning: package 'purrr' was built under R version 4.5.3

Warning: package 'dplyr' was built under R version 4.5.3

Warning: package 'stringr' was built under R version 4.5.2

Warning: package 'forcats' was built under R version 4.5.2

Warning: package 'lubridate' was built under R version 4.5.2

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ lubridate 1.9.5     ✔ tibble    3.3.1
✔ purrr     1.2.2     ✔ 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 errors

# Load data with fallback for environment compatibility

if ("ChemicalManufacturingProcess" %in% data(package = "AppliedPredictiveModeling")$results[, "Item"]) {
  data("ChemicalManufacturingProcess")
  raw_df <- ChemicalManufacturingProcess
} else {
  data("chemicalManufacturing")
  raw_df <- data.frame(Yield = yield, processPredictors)
}

dim(raw_df)

[1] 176  58

str(raw_df)

'data.frame':   176 obs. of  58 variables:
 $ Yield                 : num  38 42.4 42 41.4 42.5 ...
 $ BiologicalMaterial01  : num  6.25 8.01 8.01 8.01 7.47 6.12 7.48 6.94 6.94 6.94 ...
 $ BiologicalMaterial02  : num  49.6 61 61 61 63.3 ...
 $ BiologicalMaterial03  : num  57 67.5 67.5 67.5 72.2 ...
 $ BiologicalMaterial04  : num  12.7 14.6 14.6 14.6 14 ...
 $ BiologicalMaterial05  : num  19.5 19.4 19.4 19.4 17.9 ...
 $ BiologicalMaterial06  : num  43.7 53.1 53.1 53.1 54.7 ...
 $ BiologicalMaterial07  : num  100 100 100 100 100 100 100 100 100 100 ...
 $ BiologicalMaterial08  : num  16.7 19 19 19 18.2 ...
 $ BiologicalMaterial09  : num  11.4 12.6 12.6 12.6 12.8 ...
 $ BiologicalMaterial10  : num  3.46 3.46 3.46 3.46 3.05 3.78 3.04 3.85 3.85 3.85 ...
 $ BiologicalMaterial11  : num  138 154 154 154 148 ...
 $ BiologicalMaterial12  : num  18.8 21.1 21.1 21.1 21.1 ...
 $ ManufacturingProcess01: num  NA 0 0 0 10.7 12 11.5 12 12 12 ...
 $ ManufacturingProcess02: num  NA 0 0 0 0 0 0 0 0 0 ...
 $ ManufacturingProcess03: num  NA NA NA NA NA NA 1.56 1.55 1.56 1.55 ...
 $ ManufacturingProcess04: num  NA 917 912 911 918 924 933 929 928 938 ...
 $ ManufacturingProcess05: num  NA 1032 1004 1015 1028 ...
 $ ManufacturingProcess06: num  NA 210 207 213 206 ...
 $ ManufacturingProcess07: num  NA 177 178 177 178 178 177 178 177 177 ...
 $ ManufacturingProcess08: num  NA 178 178 177 178 178 178 178 177 177 ...
 $ ManufacturingProcess09: num  43 46.6 45.1 44.9 45 ...
 $ ManufacturingProcess10: num  NA NA NA NA NA NA 11.6 10.2 9.7 10.1 ...
 $ ManufacturingProcess11: num  NA NA NA NA NA NA 11.5 11.3 11.1 10.2 ...
 $ ManufacturingProcess12: num  NA 0 0 0 0 0 0 0 0 0 ...
 $ ManufacturingProcess13: num  35.5 34 34.8 34.8 34.6 34 32.4 33.6 33.9 34.3 ...
 $ ManufacturingProcess14: num  4898 4869 4878 4897 4992 ...
 $ ManufacturingProcess15: num  6108 6095 6087 6102 6233 ...
 $ ManufacturingProcess16: num  4682 4617 4617 4635 4733 ...
 $ ManufacturingProcess17: num  35.5 34 34.8 34.8 33.9 33.4 33.8 33.6 33.9 35.3 ...
 $ ManufacturingProcess18: num  4865 4867 4877 4872 4886 ...
 $ ManufacturingProcess19: num  6049 6097 6078 6073 6102 ...
 $ ManufacturingProcess20: num  4665 4621 4621 4611 4659 ...
 $ ManufacturingProcess21: num  0 0 0 0 -0.7 -0.6 1.4 0 0 1 ...
 $ ManufacturingProcess22: num  NA 3 4 5 8 9 1 2 3 4 ...
 $ ManufacturingProcess23: num  NA 0 1 2 4 1 1 2 3 1 ...
 $ ManufacturingProcess24: num  NA 3 4 5 18 1 1 2 3 4 ...
 $ ManufacturingProcess25: num  4873 4869 4897 4892 4930 ...
 $ ManufacturingProcess26: num  6074 6107 6116 6111 6151 ...
 $ ManufacturingProcess27: num  4685 4630 4637 4630 4684 ...
 $ ManufacturingProcess28: num  10.7 11.2 11.1 11.1 11.3 11.4 11.2 11.1 11.3 11.4 ...
 $ ManufacturingProcess29: num  21 21.4 21.3 21.3 21.6 21.7 21.2 21.2 21.5 21.7 ...
 $ ManufacturingProcess30: num  9.9 9.9 9.4 9.4 9 10.1 11.2 10.9 10.5 9.8 ...
 $ ManufacturingProcess31: num  69.1 68.7 69.3 69.3 69.4 68.2 67.6 67.9 68 68.5 ...
 $ ManufacturingProcess32: num  156 169 173 171 171 173 159 161 160 164 ...
 $ ManufacturingProcess33: num  66 66 66 68 70 70 65 65 65 66 ...
 $ ManufacturingProcess34: num  2.4 2.6 2.6 2.5 2.5 2.5 2.5 2.5 2.5 2.5 ...
 $ ManufacturingProcess35: num  486 508 509 496 468 490 475 478 491 488 ...
 $ ManufacturingProcess36: num  0.019 0.019 0.018 0.018 0.017 0.018 0.019 0.019 0.019 0.019 ...
 $ ManufacturingProcess37: num  0.5 2 0.7 1.2 0.2 0.4 0.8 1 1.2 1.8 ...
 $ ManufacturingProcess38: num  3 2 2 2 2 2 2 2 3 3 ...
 $ ManufacturingProcess39: num  7.2 7.2 7.2 7.2 7.3 7.2 7.3 7.3 7.4 7.1 ...
 $ ManufacturingProcess40: num  NA 0.1 0 0 0 0 0 0 0 0 ...
 $ ManufacturingProcess41: num  NA 0.15 0 0 0 0 0 0 0 0 ...
 $ ManufacturingProcess42: num  11.6 11.1 12 10.6 11 11.5 11.7 11.4 11.4 11.3 ...
 $ ManufacturingProcess43: num  3 0.9 1 1.1 1.1 2.2 0.7 0.8 0.9 0.8 ...
 $ ManufacturingProcess44: num  1.8 1.9 1.8 1.8 1.7 1.8 2 2 1.9 1.9 ...
 $ ManufacturingProcess45: num  2.4 2.2 2.3 2.1 2.1 2 2.2 2.2 2.1 2.4 ...

b) Using an imputation function to fill in missing values

We will use the preProcess as our imputation function

# total count of missing values across the entire predictor matrix
total_na <- sum(is.na(raw_df[, -1]))
total_na

[1] 106

# Let's generate the overall missingness rate of the predictor space
pct_na   <- mean(is.na(raw_df[, -1])) * 100
pct_na

[1] 1.056619

# Let's produce a filtered, named vector containing only columns that have missing data, alongside their exact missing value counts.

na_cols  <- colSums(is.na(raw_df[, -1]))[colSums(is.na(raw_df[, -1])) > 0]
na_cols

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) Splitting, preProcess and tune a model of our choice

We will split our data in 80/20 percent ratio where 80% training set and 20% test set.

# Splitting

set.seed(123)

train_idx <- createDataPartition(raw_df$Yield, p = 0.80, list = FALSE)
train_raw <- raw_df[train_idx, ]
test_raw  <- raw_df[-train_idx, ]

# Let's tune a PLS model using the 10-fold cross validation

# Establish 10-fold repeated Cross validation

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

# Fit and tune PLS model

set.seed(624)

pls_fit <- train(
  Yield ~ .,
  data = train_raw,
  method = "pls",
  tuneLength = 20,
  trControl = ctrl,
  preProc = c("zv", "nzv", "knnImpute"),
  na.action = na.pass
)

# Let's find the Optimal tuning parameter and resampled metric

best_ncomp <- pls_fit$bestTune$ncomp
best_ncomp

[1] 3

pls_resampled_rmse <- getTrainPerf(pls_fit)$TrainRMSE
pls_resampled_rmse

[1] 1.23011

pls_resampled_r2   <- getTrainPerf(pls_fit)$TrainRsquared
pls_resampled_r2

[1] 0.6031408

pls_resampled_mae  <- getTrainPerf(pls_fit)$TrainMAE
pls_resampled_mae

[1] 0.9769446

The optimal value of the performance metric is 3.

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?

# Let's create the prediction test set

test_preds <- predict(pls_fit, newdata = test_raw)


# Calculate out-of-sample performance metrics

test_metrics <- postResample(pred = test_preds, obs = test_raw$Yield)

Warning in pred - obs: longer object length is not a multiple of shorter object
length
Warning in pred - obs: longer object length is not a multiple of shorter object
length

test_metrics

    RMSE Rsquared      MAE 
1.830727       NA 1.483320

Interpretation

e) Identify the most important predictors

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_fit)$importance # Extract variable importance as a data frame
pls_imp_df <- data.frame(
  Predictor = rownames(pls_imp),  # Add row names as a column for predictors
  Importance = pls_imp$Overall,
  Type = ifelse(grepl("^Biological", rownames(pls_imp)), "Biological Material", "Manufacturing Process")
) %>% arrange(desc(Importance))

# Top 10 predictors

head(pls_imp_df, 10)

                Predictor Importance                  Type
1  ManufacturingProcess32  100.00000 Manufacturing Process
2  ManufacturingProcess17   87.59217 Manufacturing Process
3  ManufacturingProcess13   86.03757 Manufacturing Process
4  ManufacturingProcess09   85.72915 Manufacturing Process
5  ManufacturingProcess36   82.88148 Manufacturing Process
6  ManufacturingProcess06   68.67342 Manufacturing Process
7  ManufacturingProcess33   62.55021 Manufacturing Process
8    BiologicalMaterial06   62.01526   Biological Material
9    BiologicalMaterial03   61.32828   Biological Material
10   BiologicalMaterial08   60.43679   Biological Material