DATA 624 HW 7
Zineb Tamnat
2026-11-08
Exercise 6.2
A
data(permeability)
dim(fingerprints)## [1] 165 1107length(permeability)## [1] 165head(fingerprints[, 1:6])## X1 X2 X3 X4 X5 X6
## 1 0 0 0 0 0 1
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 1
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0head(permeability)## permeability
## 1 12.520
## 2 1.120
## 3 19.405
## 4 1.730
## 5 1.680
## 6 0.510B
set.seed(624)
# Splitting data into training and test sets first
y <- as.numeric(permeability)
train_index <- createDataPartition(
y,
p = 0.80,
list = FALSE
)
# Separating training and test predictors
x_train_original <- fingerprints[train_index, , drop = FALSE]
x_test_original <- fingerprints[-train_index, , drop = FALSE]
y_train <- y[train_index]
y_test <- y[-train_index]
# Identifying near zero variance predictors using training data only
nzv <- nearZeroVar(x_train_original)
# Removing the same predictors from both sets
x_train <- x_train_original[, -nzv, drop = FALSE]
x_test <- x_test_original[, -nzv, drop = FALSE]
cat("Original number of predictors:", ncol(fingerprints), "\n")## Original number of predictors: 1107cat("Number of predictors removed:", length(nzv), "\n")## Number of predictors removed: 700cat("Number of predictors remaining:", ncol(x_train), "\n")## Number of predictors remaining: 407The data was first split into training and test sets to avoid using test set information during predictor filtering. Using the nearZeroVar function on the training data, 700 of the original 1,107 molecular fingerprint predictors were removed due to near zero variance. This left 407 predictors for modeling. The same predictors were removed from the test set to keep consistency.
C
# Setting up 10 fold cross validation
train_control <- trainControl(
method = "cv",
number = 10
)
# Training PLS model
set.seed(624)
pls_model <- train(
x = x_train,
y = y_train,
method = "pls",
preProcess = c("center", "scale"),
tuneGrid = expand.grid(ncomp = 1:20),
trControl = train_control,
metric = "Rsquared"
)
# Displaying results
print(pls_model)## Partial Least Squares
##
## 133 samples
## 407 predictors
##
## Pre-processing: centered (407), scaled (407)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 121, 118, 119, 121, 119, 120, ...
## Resampling results across tuning parameters:
##
## ncomp RMSE Rsquared MAE
## 1 13.61193 0.2886746 10.387959
## 2 12.57515 0.3958010 9.234593
## 3 12.64539 0.4264515 9.693109
## 4 12.80453 0.4322974 9.835668
## 5 13.05880 0.4072410 9.683955
## 6 13.04766 0.4290457 9.719455
## 7 12.71229 0.4190529 9.666694
## 8 12.58213 0.4025824 9.652278
## 9 12.36276 0.4105045 9.585815
## 10 12.50241 0.3968828 9.663708
## 11 12.68950 0.3819727 9.788914
## 12 12.54813 0.3962018 9.613786
## 13 12.48248 0.4044788 9.291697
## 14 12.77969 0.3799871 9.569038
## 15 13.06601 0.3644893 9.608370
## 16 13.23343 0.3640904 9.717218
## 17 13.49368 0.3600821 9.933227
## 18 13.39359 0.3665454 9.878643
## 19 13.28448 0.3774457 9.749059
## 20 13.26420 0.3818645 9.596228
##
## Rsquared was used to select the optimal model using the largest value.
## The final value used for the model was ncomp = 4.cat("Optimal number of latent variables:",
pls_model$bestTune$ncomp, "\n")## Optimal number of latent variables: 4cat("Resampled R-squared:",
max(pls_model$results$Rsquared), "\n")## Resampled R-squared: 0.4322974plot(pls_model)
The PLS model was tuned using 10 fold cross validation with 1 to 20 latent variables. Predictors were centered and scaled before modeling. The optimal model used 4 latent variables, achieving a resampled R² of 0.4323, suggesting moderate predictive performance.
D
# Predicting permeability for the test set
pls_predictions <- predict(pls_model, newdata = x_test)
# Calculatng the test set performance
test_results <- postResample(
pred = pls_predictions,
obs = y_test
)
print(test_results)## RMSE Rsquared MAE
## 10.6166826 0.5068919 8.2109498cat("Test set R-squared:",
test_results["Rsquared"], "\n")## Test set R-squared: 0.5068919The PLS model achieved a test set R² of 0.5069, compared to the cross-validated training R² of 0.4323. The slightly higher test R² suggests that the model performed somewhat better on the held out test data. The test RMSE was 10.6167 and the MAE was 8.2109, indicating moderate predictive accuracy.
E
# Using the same cross validation settings
set.seed(624)
# Ridge regression
ridge_model <- train(
x = x_train,
y = y_train,
method = "glmnet",
preProcess = c("center", "scale"),
tuneGrid = expand.grid(
alpha = 0,
lambda = 10^seq(-3, 3, length = 50)
),
trControl = train_control,
metric = "RMSE"
)
# Lasso regression
set.seed(624)
lasso_model <- train(
x = x_train,
y = y_train,
method = "glmnet",
preProcess = c("center", "scale"),
tuneGrid = expand.grid(
alpha = 1,
lambda = 10^seq(-3, 3, length = 50)
),
trControl = train_control,
metric = "RMSE"
)
# Elastic net
set.seed(624)
elastic_model <- train(
x = x_train,
y = y_train,
method = "glmnet",
preProcess = c("center", "scale"),
tuneGrid = expand.grid(
alpha = seq(0.1, 0.9, by = 0.2),
lambda = 10^seq(-3, 3, length = 50)
),
trControl = train_control,
metric = "RMSE"
)
# Predicting on the test set
ridge_pred <- predict(ridge_model, newdata = x_test)
lasso_pred <- predict(lasso_model, newdata = x_test)
elastic_pred <- predict(elastic_model, newdata = x_test)
# Comparing test performance
model_comparison <- data.frame(
Model = c("PLS", "Ridge", "Lasso", "Elastic Net"),
RMSE = c(
RMSE(pls_predictions, y_test),
RMSE(ridge_pred, y_test),
RMSE(lasso_pred, y_test),
RMSE(elastic_pred, y_test)
),
Rsquared = c(
R2(pls_predictions, y_test),
R2(ridge_pred, y_test),
R2(lasso_pred, y_test),
R2(elastic_pred, y_test)
)
)
print(model_comparison)## Model RMSE Rsquared
## 1 PLS 10.616683 0.5068919
## 2 Ridge 9.847412 0.6311740
## 3 Lasso 11.258977 0.4669111
## 4 Elastic Net 10.555135 0.5206896Ridge regression, lasso regression and elastic net were compared with the PLS model using the same training and test sets. Ridge regression achieved the best predictive performance on the test set with the lowest RMSE of 9.8474 and the highest R² of 0.6312. Elastic net also slightly outperformed PLS, while lasso performed worse. The regularized regression models were tuned using cross validated RMSE. Overall, ridge regression was the strongest model among those evaluated on the test set.
F
Although ridge regression performed best, I would not recommend replacing laboratory permeability experiments with these models yet. The ridge model achieved a test R² of 0.6312, indicating that its predictions still have considerable uncertainty. However, the model could be useful for screening compounds and identifying promising candidates before laboratory testing. Further validation on additional compounds would be needed before considering it as a replacement.
Exercise 6.3
A
# Loading the dataset
data(ChemicalManufacturingProcess)
# Checking dimensions
dim(ChemicalManufacturingProcess)## [1] 176 58# Viewing the first six columns
head(ChemicalManufacturingProcess[, 1:6])## Yield BiologicalMaterial01 BiologicalMaterial02 BiologicalMaterial03
## 1 38.00 6.25 49.58 56.97
## 2 42.44 8.01 60.97 67.48
## 3 42.03 8.01 60.97 67.48
## 4 41.42 8.01 60.97 67.48
## 5 42.49 7.47 63.33 72.25
## 6 43.57 6.12 58.36 65.31
## BiologicalMaterial04 BiologicalMaterial05
## 1 12.74 19.51
## 2 14.65 19.36
## 3 14.65 19.36
## 4 14.65 19.36
## 5 14.02 17.91
## 6 15.17 21.79# Viewing the response variable
head(ChemicalManufacturingProcess$Yield)## [1] 38.00 42.44 42.03 41.42 42.49 43.57# Checking predictor and response counts
cat("Number of manufacturing runs:",
nrow(ChemicalManufacturingProcess), "\n")## Number of manufacturing runs: 176cat("Number of predictors:",
ncol(ChemicalManufacturingProcess) - 1, "\n")## Number of predictors: 57B
# Separating predictors and response
process_x <- ChemicalManufacturingProcess[,
setdiff(names(ChemicalManufacturingProcess), "Yield")
]
process_y <- ChemicalManufacturingProcess$Yield
# Counting missing values before imputation
cat("Missing values before imputation:",
sum(is.na(process_x)), "\n")## Missing values before imputation: 106# Imputing missing values using the median
impute_model <- preProcess(
process_x,
method = "medianImpute"
)
process_x_imputed <- predict(impute_model, process_x)
# Counting missing values after imputation
cat("Missing values after imputation:",
sum(is.na(process_x_imputed)), "\n")## Missing values after imputation: 0# Checking dimensions
dim(process_x_imputed)## [1] 176 57The predictor dataset contained 106 missing values. Median imputation was used to replace missing values with the median of each respective predictor. After imputation, no missing values remained and all 57 predictors were retained.
For model training and evaluation, imputation was performed within the modeling workflow using the training data rather than the fully imputed dataset to avoid introducing information from the test set.
C
set.seed(624)
# Splitting into training and test sets
process_train_index <- createDataPartition(
process_y,
p = 0.80,
list = FALSE
)
process_x_train <- process_x[
process_train_index, , drop = FALSE
]
process_x_test <- process_x[
-process_train_index, , drop = FALSE
]
process_y_train <- process_y[process_train_index]
process_y_test <- process_y[-process_train_index]
# Setting up 10 fold cross validation
process_control <- trainControl(
method = "cv",
number = 10
)
# Training ridge regression with preprocessing
set.seed(624)
process_ridge_model <- train(
x = process_x_train,
y = process_y_train,
method = "glmnet",
preProcess = c("medianImpute", "center", "scale"),
tuneGrid = expand.grid(
alpha = 0,
lambda = 10^seq(-3, 3, length.out = 50)
),
trControl = process_control,
metric = "RMSE"
)
# Displaying model results
print(process_ridge_model)## glmnet
##
## 144 samples
## 57 predictor
##
## Pre-processing: median imputation (57), centered (57), scaled (57)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 130, 130, 130, 130, 129, 130, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared MAE
## 1.000000e-03 2.058464 0.4926793 1.212355
## 1.325711e-03 2.058464 0.4926793 1.212355
## 1.757511e-03 2.058464 0.4926793 1.212355
## 2.329952e-03 2.058464 0.4926793 1.212355
## 3.088844e-03 2.058464 0.4926793 1.212355
## 4.094915e-03 2.058464 0.4926793 1.212355
## 5.428675e-03 2.058464 0.4926793 1.212355
## 7.196857e-03 2.058464 0.4926793 1.212355
## 9.540955e-03 2.058464 0.4926793 1.212355
## 1.264855e-02 2.058464 0.4926793 1.212355
## 1.676833e-02 2.058464 0.4926793 1.212355
## 2.222996e-02 2.058464 0.4926793 1.212355
## 2.947052e-02 2.058464 0.4926793 1.212355
## 3.906940e-02 2.058464 0.4926793 1.212355
## 5.179475e-02 2.058464 0.4926793 1.212355
## 6.866488e-02 2.058464 0.4926793 1.212355
## 9.102982e-02 2.058464 0.4926793 1.212355
## 1.206793e-01 2.014197 0.4947619 1.198042
## 1.599859e-01 1.910358 0.4982735 1.167054
## 2.120951e-01 1.805335 0.5008074 1.140723
## 2.811769e-01 1.702731 0.5022961 1.114853
## 3.727594e-01 1.606130 0.5027752 1.089659
## 4.941713e-01 1.518812 0.5026634 1.067737
## 6.551286e-01 1.443972 0.5026800 1.049265
## 8.685114e-01 1.383586 0.5043035 1.036697
## 1.151395e+00 1.339829 0.5096568 1.031175
## 1.526418e+00 1.317197 0.5203688 1.030800
## 2.023590e+00 1.318298 0.5334747 1.040196
## 2.682696e+00 1.338678 0.5373449 1.053640
## 3.556480e+00 1.366675 0.5273740 1.081975
## 4.714866e+00 1.394308 0.5116532 1.112114
## 6.250552e+00 1.418874 0.4968682 1.140397
## 8.286428e+00 1.441438 0.4846370 1.169374
## 1.098541e+01 1.464119 0.4746166 1.199496
## 1.456348e+01 1.488825 0.4658162 1.228336
## 1.930698e+01 1.516505 0.4571843 1.256775
## 2.559548e+01 1.546823 0.4480554 1.283172
## 3.393222e+01 1.578416 0.4387056 1.309155
## 4.498433e+01 1.609526 0.4299462 1.332114
## 5.963623e+01 1.638583 0.4224297 1.352228
## 7.906043e+01 1.664526 0.4163660 1.370009
## 1.048113e+02 1.686852 0.4116555 1.385121
## 1.389495e+02 1.705513 0.4080617 1.400895
## 1.842070e+02 1.720754 0.4053458 1.413597
## 2.442053e+02 1.732978 0.4033021 1.423593
## 3.237458e+02 1.742644 0.4017718 1.431484
## 4.291934e+02 1.750200 0.4006162 1.437768
## 5.689866e+02 1.756059 0.3997483 1.442937
## 7.543120e+02 1.760572 0.3990958 1.446926
## 1.000000e+03 1.768884 0.3570241 1.454233
##
## Tuning parameter 'alpha' was held constant at a value of 0
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0 and lambda = 1.526418.cat("Optimal lambda:",
process_ridge_model$bestTune$lambda, "\n")## Optimal lambda: 1.526418cat("Optimal cross-validated RMSE:",
min(process_ridge_model$results$RMSE), "\n")## Optimal cross-validated RMSE: 1.317197plot(process_ridge_model)
A ridge regression model was trained using an 80/20 training test split and 10 fold cross validation. Missing values were handled using median imputation and predictors were centered and scaled. The optimal regularization parameter was λ = 1.5264 which produced a cross validated RMSE of 1.3172. The cross validated RMSE of 1.3172 indicates a typical prediction error of approximately 1.32 percentage points in manufacturing yield as measured by root mean squared error.
D
# Predicting yield for the test set
process_predictions <- predict(
process_ridge_model,
newdata = process_x_test
)
# Calculating test set performance
process_test_results <- postResample(
pred = process_predictions,
obs = process_y_test
)
print(process_test_results)## RMSE Rsquared MAE
## 1.4867793 0.5058082 1.0652821# Comparing cross validation and test RMSE
cat("Cross-validated RMSE:",
min(process_ridge_model$results$RMSE), "\n")## Cross-validated RMSE: 1.317197cat("Test set RMSE:",
process_test_results["RMSE"], "\n")## Test set RMSE: 1.486779The ridge regression model achieved a test set RMSE of 1.4868 compared to the cross validated training RMSE of 1.3172. The slightly higher test RMSE indicates that the model performed somewhat worse on unseen data although the results were relatively close. The test R² was 0.5058 suggesting moderate predictive performance.
E
# Calculating variable importance
process_importance <- varImp(process_ridge_model, scale = TRUE)
# Extracting and sort importance scores
importance_df <- data.frame(
Predictor = rownames(process_importance$importance),
Importance = process_importance$importance$Overall
)
importance_df <- importance_df[
order(-importance_df$Importance),
]
# Identifying predictor type
importance_df$Type <- ifelse(
grepl("^BiologicalMaterial", importance_df$Predictor),
"Biological",
"Process"
)
# Displaying top 10 predictors
top_10 <- head(importance_df, 10)
print(top_10)## Predictor Importance Type
## 44 ManufacturingProcess32 100.00000 Process
## 21 ManufacturingProcess09 71.22841 Process
## 48 ManufacturingProcess36 67.09188 Process
## 46 ManufacturingProcess34 63.60742 Process
## 29 ManufacturingProcess17 53.75954 Process
## 25 ManufacturingProcess13 53.65121 Process
## 49 ManufacturingProcess37 45.97006 Process
## 18 ManufacturingProcess06 45.11969 Process
## 3 BiologicalMaterial03 42.03703 Biological
## 6 BiologicalMaterial06 39.03873 Biological# Counting predictor types in top 10
table(top_10$Type)##
## Biological Process
## 2 8# Plotting top 10 predictors
ggplot(top_10, aes(
x = reorder(Predictor, Importance),
y = Importance,
fill = Type
)) +
geom_col() +
coord_flip() +
scale_fill_manual(
values = c(
"Biological" = "#8DAA91",
"Process" = "#633B70"
)
) +
labs(
title = "Top 10 Most Important Predictors",
x = "Predictor",
y = "Variable Importance"
) +
theme_minimal()
The variable importance analysis showed that manufacturing process
predictors dominated the model accounting for 8 of the top 10 most
important variables. ManufacturingProcess32 was the most important
predictor followed by ManufacturingProcess09 and ManufacturingProcess36.
Only two biological predictors, BiologicalMaterial03 and
BiologicalMaterial06 appeared in the top 10. This suggests that
manufacturing process variables may be particularly useful for
predicting product yield.
F
# Selecting top 10 predictors
top_predictors <- top_10$Predictor
# Create dataset for plotting
plot_data <- ChemicalManufacturingProcess[
, c("Yield", top_predictors)
]
# Convert to long format
library(tidyr)
plot_long <- pivot_longer(
plot_data,
cols = all_of(top_predictors),
names_to = "Predictor",
values_to = "Value"
)
# Plot predictor relationships with yield
ggplot(plot_long, aes(x = Value, y = Yield)) +
geom_point(alpha = 0.6, color = "darkgreen") +
geom_smooth(
method = "lm",
se = FALSE,
color = "black",
na.rm = TRUE
) +
facet_wrap(~ Predictor, scales = "free_x", ncol = 2) +
labs(
title = "Relationships Between Top Predictors and Yield",
x = "Predictor Value",
y = "Product Yield (%)"
) +
theme_minimal()
The scatterplots show that several important predictors have relationships with product yield. ManufacturingProcess32, ManufacturingProcess09 and ManufacturingProcess06 show positive relationships, while ManufacturingProcess13, ManufacturingProcess17 and ManufacturingProcess36 show negative relationships. ManufacturingProcess34 and ManufacturingProcess37 have relatively weak relationships with yield. Both biological predictors also show positive trends.
These findings could help improve future manufacturing yields by identifying process variables worth investigating. For example, increasing ManufacturingProcess32 or decreasing ManufacturingProcess36 is worth testing. Since manufacturing process variables can be adjusted, they provide opportunities for improving yield. However, further experiments would be needed to determine whether changing these variables actually causes an increase in yield.