data624_hw07

6.2

a.

Loading the packages:

Loading required package: ggplot2

Loading required package: lattice

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

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

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

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

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

The data has the following shape:

dim(fingerprints)

[1]  165 1107

paste('The fingerprints matrix has', ncol(fingerprints), 'columns.')

[1] "The fingerprints matrix has 1107 columns."

paste('The fingerprints matrix has', nrow(fingerprints), 'rows.')

[1] "The fingerprints matrix has 165 rows."

b.

fingerprints <- as.data.frame(fingerprints)

zero_var_cols <- nearZeroVar(fingerprints)
zero_var_names <- names(fingerprints)[zero_var_cols]

df <- fingerprints |> 
  select(-all_of(zero_var_cols))

paste('After applying nearZeroVar() the fingerprints data has', ncol(df), 'features left as predictors, and', nrow(df), 'rows. ')

[1] "After applying nearZeroVar() the fingerprints data has 388 features left as predictors, and 165 rows. "

After applying nearZeroVar() the data went from 1107 to 388 features.

c.

The data is split into a training and testing set.

set.seed(1234)

inTrain <- createDataPartition(
  y = permeability,
  p = 0.75,
  list = FALSE
)

training <- df[inTrain,]
testing <- df[-inTrain,]

train_y <- permeability[inTrain]
test_y <- permeability[-inTrain]

Then a Partial Least Squares model is trained.

plsFit <- train(
  x = training,
  y = train_y,
  method = "pls",
  preProc = c("zv", "center", "scale"),
  tuneGrid = data.frame(ncomp = 1:20)
)

summary(plsFit)

Data:   X dimension: 125 388 
    Y dimension: 125 1
Fit method: oscorespls
Number of components considered: 2
TRAINING: % variance explained
          1 comps  2 comps
X           22.87    34.97
.outcome    27.94    44.42

plsFit

Partial Least Squares 

125 samples
388 predictors

Pre-processing: centered (388), scaled (388) 
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 125, 125, 125, 125, 125, 125, ... 
Resampling results across tuning parameters:

  ncomp  RMSE      Rsquared   MAE      
   1     13.42334  0.2672258  10.072437
   2     13.29547  0.3138451   9.386395
   3     13.80270  0.2783473   9.989636
   4     13.77598  0.2957272   9.980013
   5     13.51692  0.3298663   9.822628
   6     13.45904  0.3406477   9.899885
   7     13.41837  0.3546579  10.037745
   8     13.32789  0.3691099   9.997005
   9     13.57084  0.3644617  10.261974
  10     13.91231  0.3569511  10.578343
  11     14.00905  0.3560731  10.711118
  12     14.35028  0.3435015  10.968892
  13     14.51800  0.3417688  11.141786
  14     14.60902  0.3447673  11.224760
  15     14.73226  0.3468174  11.328549
  16     14.82189  0.3451504  11.374250
  17     14.93592  0.3444535  11.434104
  18     15.14040  0.3416878  11.592096
  19     15.26297  0.3409856  11.705482
  20     15.40570  0.3412974  11.834812

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

The optimal number of latent variables is 2, and the corresponding resampled R² value is 0.3138451.

d.

The held out set is predicted below using the fitted model.

predictions <- predict(plsFit, newdata = testing)
pls_best <- postResample(pred = predictions, obs = test_y)
pls_best

      RMSE   Rsquared        MAE 
10.2050816  0.6372984  7.1543535 

The R² value is improved at 0.6372984 for the test set, compared to 0.3138451 for the training set.

e.

Ordinary Linear Regression

The data has more predictors than observations. As a result Ordinary Least Squares is not a practical choice for regression without reducing the features.

Ridge Regression

Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

Loaded glmnet 4.1-10

A ridge regression using 10 fold cross validation is trained below.

ctrl <- trainControl(
  method = 'cv',
  number = 10
)

ridgeGrid <- expand.grid(
  lambda = 10^seq(-4, 1, length.out = 50)
)

set.seed(1234)

ridge_caret_model <- train(
  x = training,
  y = train_y,
  method = "ridge",
  preProcess = c('zv', 'center', 'scale'),
  trControl = ctrl,
  tuneGrid = ridgeGrid
  )

ridge_caret_model$bestTune

      lambda
32 0.1456348

plot(ridge_caret_model)

The lowest RMSE and lambda value from model training is shown below.

ridge_caret_model$results[
  which.min(ridge_caret_model$results$RMSE),
]

      lambda     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
32 0.1456348 12.33407 0.4591478 9.774303 3.242247  0.2663755 2.797893

predictions <- predict(ridge_caret_model, newdata = testing)
postResample(pred = predictions, obs = test_y)

      RMSE   Rsquared        MAE 
11.3256612  0.5824795  8.2003048 

Below is a second version of a ridge regression using the “glmnet” method instead of the “ridge” method.

ctrl <- trainControl(method = 'cv', number = 10)

lambdas <- 10^seq(3, -3, length.out = 100)

set.seed(1234)

ridge_model <- train(
  x = training,
  y = train_y,
  method = "glmnet",
  trControl = ctrl,
  preProcess = c('zv', 'center', 'scale'),
  tuneGrid = expand.grid(alpha = 0, lambda = lambdas),
  lambda = lambdas
)

ridge_model 

glmnet 

125 samples
388 predictors

Pre-processing: centered (388), scaled (388) 
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 112, 112, 112, 113, 113, 113, ... 
Resampling results across tuning parameters:

  lambda        RMSE      Rsquared   MAE      
  1.000000e-03  17.91421  0.3018878  13.476126
  1.149757e-03  17.90179  0.3020367  13.468748
  1.321941e-03  17.88791  0.3022109  13.460431
  1.519911e-03  17.87108  0.3024286  13.450340
  1.747528e-03  17.85231  0.3026612  13.438968
  2.009233e-03  17.83072  0.3029642  13.425896
  2.310130e-03  17.80765  0.3032585  13.411947
  2.656088e-03  17.78195  0.3036264  13.396631
  3.053856e-03  17.75286  0.3040345  13.379124
  3.511192e-03  17.71927  0.3045230  13.358738
  4.037017e-03  17.68296  0.3050378  13.336848
  4.641589e-03  17.64170  0.3056475  13.311815
  5.336699e-03  17.59532  0.3063371  13.283437
  6.135907e-03  17.54360  0.3071230  13.251785
  7.054802e-03  17.48522  0.3080244  13.215712
  8.111308e-03  17.42025  0.3090539  13.175299
  9.326033e-03  17.34941  0.3102003  13.131918
  1.072267e-02  17.26934  0.3115431  13.082200
  1.232847e-02  17.18330  0.3130109  13.028430
  1.417474e-02  17.08801  0.3146876  12.968052
  1.629751e-02  16.98357  0.3165758  12.901175
  1.873817e-02  16.86977  0.3187080  12.827609
  2.154435e-02  16.74735  0.3210744  12.756433
  2.477076e-02  16.61794  0.3236745  12.679241
  2.848036e-02  16.48186  0.3264847  12.595901
  3.274549e-02  16.34017  0.3294482  12.508143
  3.764936e-02  16.18935  0.3326848  12.415271
  4.328761e-02  16.02822  0.3362687  12.318108
  4.977024e-02  15.85907  0.3401884  12.218362
  5.722368e-02  15.68202  0.3444386  12.113269
  6.579332e-02  15.49776  0.3490143  12.000774
  7.564633e-02  15.30682  0.3539308  11.880103
  8.697490e-02  15.11271  0.3590494  11.754052
  1.000000e-01  14.91474  0.3644514  11.627436
  1.149757e-01  14.71501  0.3700247  11.493573
  1.321941e-01  14.51550  0.3757234  11.352915
  1.519911e-01  14.31650  0.3815163  11.204916
  1.747528e-01  14.12090  0.3872997  11.053684
  2.009233e-01  13.92811  0.3930788  10.898360
  2.310130e-01  13.74163  0.3987483  10.747456
  2.656088e-01  13.56012  0.4043313  10.607974
  3.053856e-01  13.38308  0.4097719  10.464839
  3.511192e-01  13.21257  0.4150208  10.325749
  4.037017e-01  13.04986  0.4200969  10.199979
  4.641589e-01  12.89496  0.4248070  10.088683
  5.336699e-01  12.74765  0.4293836   9.976139
  6.135907e-01  12.61009  0.4335097   9.875478
  7.054802e-01  12.48145  0.4375468   9.791077
  8.111308e-01  12.35943  0.4411163   9.702813
  9.326033e-01  12.24743  0.4446508   9.622018
  1.072267e+00  12.14504  0.4474889   9.536351
  1.232847e+00  12.04744  0.4504906   9.455762
  1.417474e+00  11.95964  0.4529845   9.377367
  1.629751e+00  11.87839  0.4552719   9.304873
  1.873817e+00  11.80471  0.4572512   9.242389
  2.154435e+00  11.73831  0.4588912   9.181484
  2.477076e+00  11.67848  0.4602584   9.123768
  2.848036e+00  11.62556  0.4612874   9.065519
  3.274549e+00  11.57874  0.4620206   9.006663
  3.764936e+00  11.53697  0.4625275   8.952320
  4.328761e+00  11.50132  0.4627039   8.893705
  4.977024e+00  11.47050  0.4626578   8.841265
  5.722368e+00  11.44412  0.4624053   8.794204
  6.579332e+00  11.42419  0.4618377   8.762348
  7.564633e+00  11.40795  0.4611197   8.725180
  8.697490e+00  11.39773  0.4601553   8.691106
  1.000000e+01  11.39029  0.4591872   8.656272
  1.149757e+01  11.38992  0.4579119   8.621826
  1.321941e+01  11.39198  0.4566779   8.585684
  1.519911e+01  11.40111  0.4551777   8.553736
  1.747528e+01  11.41260  0.4537726   8.522714
  2.009233e+01  11.42966  0.4522275   8.512529
  2.310130e+01  11.45179  0.4505775   8.510735
  2.656088e+01  11.47766  0.4489023   8.506474
  3.053856e+01  11.50745  0.4472585   8.509234
  3.511192e+01  11.54118  0.4455975   8.513689
  4.037017e+01  11.57854  0.4440084   8.519866
  4.641589e+01  11.61924  0.4424580   8.533368
  5.336699e+01  11.66250  0.4410304   8.560304
  6.135907e+01  11.70963  0.4396114   8.593711
  7.054802e+01  11.75980  0.4383094   8.625096
  8.111308e+01  11.81362  0.4370853   8.663317
  9.326033e+01  11.87102  0.4359475   8.718228
  1.072267e+02  11.93211  0.4348856   8.782535
  1.232847e+02  11.99836  0.4338682   8.862434
  1.417474e+02  12.06879  0.4328758   8.946798
  1.629751e+02  12.14418  0.4318705   9.036693
  1.873817e+02  12.22474  0.4308128   9.131354
  2.154435e+02  12.31057  0.4296660   9.232948
  2.477076e+02  12.40177  0.4284000   9.338132
  2.848036e+02  12.49835  0.4269891   9.451584
  3.274549e+02  12.60024  0.4254099   9.564515
  3.764936e+02  12.70717  0.4236369   9.683812
  4.328761e+02  12.81867  0.4216411   9.804000
  4.977024e+02  12.93407  0.4193957   9.922333
  5.722368e+02  13.05249  0.4168852  10.038184
  6.579332e+02  13.17295  0.4141119  10.150957
  7.564633e+02  13.29523  0.4112025  10.261096
  8.697490e+02  13.41612  0.4079405  10.366884
  1.000000e+03  13.53567  0.4043057  10.468823

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 = 11.49757.

plot(ridge_model)

ridge_model$results[
  which.min(ridge_model$results$RMSE),
]

   alpha   lambda     RMSE  Rsquared      MAE   RMSESD RsquaredSD    MAESD
68     0 11.49757 11.38992 0.4579119 8.621826 3.003354  0.2797527 2.335832

predictions <- predict(ridge_model, newdata = testing)
ridge_best <- postResample(pred = predictions, obs = test_y)
ridge_best

     RMSE  Rsquared       MAE 
9.9850654 0.6338777 7.2552879 

This second ridge version performs slightly better than the first version.

Lasso Model

cv_fit <- cv.glmnet(x = as.matrix(training), y = train_y, alpha = 1)
opt_lambda <- cv_fit$lambda.min

print(opt_lambda)

[1] 0.7796381

cv_fit

Call:  cv.glmnet(x = as.matrix(training), y = train_y, alpha = 1) 

Measure: Mean-Squared Error 

    Lambda Index Measure    SE Nonzero
min 0.7796    53   142.3 21.87      37
1se 2.8678    25   163.1 26.94       7

cv_rmse <- sqrt(
  cv_fit$cvm[which.min(cv_fit$cvm)]
)

cv_rmse

[1] 11.92765

lasso_pred <- predict(
  cv_fit,
  newx = as.matrix(testing),
  s = "lambda.min"
)

lasso_best <- postResample(
  pred = as.numeric(lasso_pred),
  obs = test_y 
)

lasso_best

      RMSE   Rsquared        MAE 
11.9434642  0.4788445  8.7477749 

Elastic Net

An elastic net model is fit below.

enetGrid <- expand.grid(
  lambda = c(0.001, 0.01, 0.1, 1),
  fraction = seq(0.05, 1, length.out = 20)
)

set.seed(1234)

enet_model <- train(
  x = training,
  y = train_y,
  method = "enet",
  trControl = ctrl,
  preProcess = c("zv", "center", "scale"),
  tuneGrid = enetGrid
)

enet_model$bestTune

   fraction lambda
43     0.15    0.1

best_fraction <- enet_model$bestTune$fraction
best_lambda <- enet_model$bestTune$lambda
plot(enet_model)

enet_pred <- predict(enet_model, newdata = testing)

best_enet <- postResample(
  pred = enet_pred,
  obs = test_y
)

best_enet

      RMSE   Rsquared        MAE 
11.6355832  0.5054059  8.4412912 

df <- bind_rows(
  "Elastic_net" = best_enet, 
  "Lasso" = lasso_best, 
  "Ridge" = ridge_best, 
  "PLS" = pls_best,
  .id = "Model") |> 
  arrange(RMSE)

df

# A tibble: 4 × 4
  Model        RMSE Rsquared   MAE
  <chr>       <dbl>    <dbl> <dbl>
1 Ridge        9.99    0.634  7.26
2 PLS         10.2     0.637  7.15
3 Elastic_net 11.6     0.505  8.44
4 Lasso       11.9     0.479  8.75

Inspecting the results from the different models, the Ridge model performed the best in terms of the root mean square error, outperforming the PLS model by a small margin.

f.

All of the models had RMSE values within a range of 2. These values are relatively close. In particular the Ridge and PLS are very close in both RMSE and R² values. The Ridge value could be a useful substitute, however it would not likely make a notable difference.

6.3

a.

library(AppliedPredictiveModeling)
library(impute)
data(ChemicalManufacturingProcess)

chem <- ChemicalManufacturingProcess

b.

The following code imputes missing values with the median value for each feature. However, a different method will be used when training the model below as part of the preProcess.

chem_impute <- ChemicalManufacturingProcess |>
  mutate(
    across(
      -Yield,
      ~ ifelse(is.na(.x), median(.x, na.rm = TRUE), .x)
    )
  )

c.

The following code splits the data and trains a ridge model on the training data set.

set.seed(1234)

split <- createDataPartition(
  y = chem$Yield,
  p = 0.8,
  list = FALSE
)

train_chem <- chem[split,]
test_chem <- chem[-split,]

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

lambdas <- 10^seq(3, -3, length.out = 100)

set.seed(1234)

chem_ridge_model <- train(
  x = train_chem |> select(-Yield),
  y = train_chem$Yield,
  method = "glmnet",
  trControl = ctrl,
  preProcess = c('zv', 'knnImpute'),
  tuneGrid = expand.grid(alpha = 0, lambda = lambdas),
  lambda = lambdas
)

chem_ridge_model$bestTune

   alpha   lambda
53     0 1.417474

plot(chem_ridge_model)

chem_ridge_model$results[
  which.min(chem_ridge_model$results$RMSE),
]

   alpha   lambda     RMSE  Rsquared       MAE    RMSESD RsquaredSD     MAESD
53     0 1.417474 1.216674 0.6069249 0.9952091 0.2286583  0.1696863 0.1815832

The model selected a lambda value of 1.417474 and an RMSE of 1.246949.

d.

The model is then used to predict Yield based on the held out testing data.

chem_pred <- predict(
  chem_ridge_model,
  newdata = test_chem |> select(-Yield)
)

postResample(pred = as.numeric(chem_pred), obs = test_chem$Yield)

     RMSE  Rsquared       MAE 
2.1062984 0.1448504 1.2953941 

When comparing the results of the RMSE on the training data and the testing data, the model had a 1.247 RMSE on the training data, and a 2.1063 RMSE on the testing data. Additionally the R² on training data was 0.5866 and on the testing data it was 0.14485. The test results are notably weaker than on the training set indicating limited generalization of this training/test split.

e.

The following code creates a dataframe of the coefficients and sorts them in descending order by the absolute value of the coefficient.

chem_coef <- coef(
  chem_ridge_model$finalModel,
  s = chem_ridge_model$bestTune$lambda
)

coef_matrix <- as.matrix(chem_coef)

coef_df <- data.frame(
  Feature = rownames(coef_matrix),
  Coefficient = coef_matrix[,1]
) |> 
  filter(Feature != "(Intercept)") |> 
  arrange(desc(abs(Coefficient)))

head(coef_df, 10)

                                      Feature Coefficient
ManufacturingProcess32 ManufacturingProcess32  0.27528220
ManufacturingProcess09 ManufacturingProcess09  0.22267584
ManufacturingProcess36 ManufacturingProcess36 -0.19881229
ManufacturingProcess13 ManufacturingProcess13 -0.18956699
ManufacturingProcess17 ManufacturingProcess17 -0.18518187
ManufacturingProcess37 ManufacturingProcess37 -0.14380561
ManufacturingProcess06 ManufacturingProcess06  0.13584095
ManufacturingProcess34 ManufacturingProcess34  0.13290285
ManufacturingProcess11 ManufacturingProcess11  0.11042892
ManufacturingProcess15 ManufacturingProcess15  0.09976618

Reviewing the list of coefficients, the ten largest values are all associated with a Manufacturing process, after which there appears to be a mix of Biological and Manufacturing.

f.

top_features <- head(coef_df$Feature, 10)

train_chem |>
  select(Yield, all_of(top_features)) |>
  pivot_longer(
    cols = -Yield,
    names_to = "Predictor",
    values_to = "Value"
  ) |>
  ggplot(aes(x = Value, y = Yield)) +
  geom_point(alpha = 0.6, na.rm = TRUE) +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = FALSE,
    na.rm = TRUE
  ) +
  facet_wrap(~ Predictor, scales = "free_x", ncol = 2) +
  labs(
    x = "Predictor value",
    y = "Yield",
    title = "Relationships between top predictors and Yield"
  ) +
  theme_minimal()

In exploring the relationship between yield and the coefficients for the top ten coefficients by absolute value. In examining the processes, improvements in ones that have the most extreme absolute values will have mean the best opportunity for improvements. By focusing attention or resources to improving those manufacturing processes there is the most potential for improving yield.