How to extract PLSR coefficients as for glmnet using tidymodels

Question

I tuned a glmnet regression model and extracted the coefficients as described here. That works wonderfully. However, when I use the same form of coefficient extraction for PLSR with mixOmics engine, I obtain single values per term and component as demonstrated here. For further external use I need the coefficients of PLSR in the first form. I can achieve this by using the optimal hyperparamterset with the plsr() function from the pls package and then extracting it with coef() as shown at the end of the code below. However, I would like to avoid this extra step because I cannot pass parameters like predictor_prop to plsr and thus results may vary.

Is there a more elegant way to extract the overall model coefficients of the PLSR as for glmnet or can I calculate them from the component values?

library(tidymodels)
library(plsmod)

data(Chicago)
Chicago <- Chicago %>% select(ridership, Clark_Lake, Austin, Harlem)

# create cross-validation dataset
folds <- vfold_cv(Chicago)

# create recipe
rec <- recipe(ridership ~ ., Chicago) %>% 
  step_normalize(all_predictors()) %>%
  prep(training = Chicago)

# define model
mod <- parsnip::pls(mode = "regression", 
                    num_comp = tune(), 
                    predictor_prop = tune()) %>% 
  set_engine("mixOmics")

# define workflow
wf <- workflow() %>%
  add_recipe(rec) %>%
  add_model(mod)

# run grid tuning
set.seed(123)
res <- tune_grid(wf, resamples = folds, grid = 5)

# get best model
res_best <- res %>% select_best("rmse")

# fit best model and extract coefficients
wf %>% 
  finalize_workflow(res_best) %>%
  fit(Chicago) %>%
  extract_fit_parsnip() %>% 
  tidy()

# extracting coefficients using plsr from pls package and coef function
p <- pls::plsr(ridership ~ ., data = Chicago, scale = T, center = T, ncomp = 3)
coef(p, intercept = T)

Thank you for the awesome tidymodels framework and everyone who makes it what it is!

Answer 1

As far as I can tell you are doing the correct thing.

The coef() functions only shows you the result for 3 comps but you can get the same result by adding filter(component == 3) in the following code

wf %>% 
  finalize_workflow(res_best) %>%
  fit(Chicago) %>%
  extract_fit_parsnip() %>% 
  tidy() %>%
  filter(component == 3)
#> # A tibble: 4 × 4
#>   term       value type       component
#>   <chr>      <dbl> <chr>          <dbl>
#> 1 Clark_Lake     0 predictors         3
#> 2 Austin         1 predictors         3
#> 3 Harlem         0 predictors         3
#> 4 Y              1 outcomes           3

The reason why you are getting 0s and 1s is because the hyper parameter tuned values of predictor_prop is quite low, giving you are sparse representation

library(tidymodels)
library(plsmod)

data(Chicago)
Chicago <- Chicago %>% select(ridership, Clark_Lake, Austin, Harlem)

# create cross-validation dataset
folds <- vfold_cv(Chicago)

# create recipe
rec <- recipe(ridership ~ ., Chicago) %>% 
  step_normalize(all_predictors()) %>%
  prep(training = Chicago)

# define model
mod <- parsnip::pls(mode = "regression", 
                    num_comp = tune(), 
                    predictor_prop = tune()) %>% 
  set_engine("mixOmics")

# define workflow
wf <- workflow() %>%
  add_recipe(rec) %>%
  add_model(mod)

# run grid tuning
set.seed(123)
res <- tune_grid(wf, resamples = folds, grid = 5)

# get best model
res_best <- res %>% select_best("rmse")

res_best
#> # A tibble: 1 × 3
#>   predictor_prop num_comp .config             
#>            <dbl>    <int> <chr>               
#> 1         0.0869        3 Preprocessor1_Model1

# fit best model and extract coefficients
wf %>% 
  finalize_workflow(res_best) %>%
  fit(Chicago) %>%
  extract_fit_parsnip() %>%
  tidy()
#> # A tibble: 12 × 4
#>    term       value type       component
#>    <chr>      <dbl> <chr>          <dbl>
#>  1 Clark_Lake     1 predictors         1
#>  2 Clark_Lake     0 predictors         2
#>  3 Clark_Lake     0 predictors         3
#>  4 Austin         0 predictors         1
#>  5 Austin         0 predictors         2
#>  6 Austin         1 predictors         3
#>  7 Harlem         0 predictors         1
#>  8 Harlem        -1 predictors         2
#>  9 Harlem         0 predictors         3
#> 10 Y              1 outcomes           1
#> 11 Y              1 outcomes           2
#> 12 Y              1 outcomes           3

wf %>% 
  finalize_workflow(
    tibble(predictor_prop = 0, num_comp = 3)
  ) %>%
  fit(Chicago) %>%
  extract_fit_parsnip() %>%
  tidy()
#> # A tibble: 12 × 4
#>    term       value type       component
#>    <chr>      <dbl> <chr>          <dbl>
#>  1 Clark_Lake     1 predictors         1
#>  2 Clark_Lake     0 predictors         2
#>  3 Clark_Lake     0 predictors         3
#>  4 Austin         0 predictors         1
#>  5 Austin         0 predictors         2
#>  6 Austin         1 predictors         3
#>  7 Harlem         0 predictors         1
#>  8 Harlem        -1 predictors         2
#>  9 Harlem         0 predictors         3
#> 10 Y              1 outcomes           1
#> 11 Y              1 outcomes           2
#> 12 Y              1 outcomes           3

wf %>% 
  finalize_workflow(
    tibble(predictor_prop = 0.5, num_comp = 3)
  ) %>%
  fit(Chicago) %>%
  extract_fit_parsnip() %>%
  tidy()
#> # A tibble: 12 × 4
#>    term        value type       component
#>    <chr>       <dbl> <chr>          <dbl>
#>  1 Clark_Lake  0.908 predictors         1
#>  2 Clark_Lake  0     predictors         2
#>  3 Clark_Lake  0     predictors         3
#>  4 Austin      0.419 predictors         1
#>  5 Austin     -0.859 predictors         2
#>  6 Austin     -0.406 predictors         3
#>  7 Harlem      0     predictors         1
#>  8 Harlem     -0.513 predictors         2
#>  9 Harlem      0.914 predictors         3
#> 10 Y           1     outcomes           1
#> 11 Y           1     outcomes           2
#> 12 Y           1     outcomes           3

wf %>% 
  finalize_workflow(
    tibble(predictor_prop = 1, num_comp = 3)
  ) %>%
  fit(Chicago) %>%
  extract_fit_parsnip() %>%
  tidy()
#> # A tibble: 12 × 4
#>    term        value type       component
#>    <chr>       <dbl> <chr>          <dbl>
#>  1 Clark_Lake  0.593 predictors         1
#>  2 Clark_Lake  0.738 predictors         2
#>  3 Clark_Lake  0.321 predictors         3
#>  4 Austin      0.576 predictors         1
#>  5 Austin     -0.111 predictors         2
#>  6 Austin     -0.810 predictors         3
#>  7 Harlem      0.562 predictors         1
#>  8 Harlem     -0.665 predictors         2
#>  9 Harlem      0.491 predictors         3
#> 10 Y           1     outcomes           1
#> 11 Y           1     outcomes           2
#> 12 Y           1     outcomes           3

^{Created on 2022-09-08 by the reprex package (v2.0.1)}