Linear and Bayesian Regression Models with tidymodels package

Make random forest model

We specify the model using the parsnip package (Kuhn and Vaughan 2020a). This package provides a tidy, unified interface to models for a range of models without getting bogged down in the syntactical minutiae of the underlying packages. The parsnip package help us to specify ;

• the type of model e.g random forest,
• the mode of the model whether is regression or classification
• the computational engine is the name of the R package.

Based on the information above, we can use parsnip package to build our model as;

rf.model = rand_forest() %>%
  set_engine(engine = "ranger") %>%
  set_mode(mode = "regression")

rf.model

Random Forest Model Specification (regression)

Computational engine: ranger 

Train random forest model

Once we have specified the model type, engine and mode, the model can be trained with the fit function. We simply parse into the fit the specified model and the transformed training set extracted from the prepped recipe.

rf.fit = rf.model %>%
  fit(chl ~ ., data = crtr.training)

predict with random forest

To get our predicted results, we use the predict() function to find the estimated chlorophyll-a. First, let’s generate the estimated chlorophyll-a values by simply parse the random forest model rf.model we specified and the transformed testing set we created from a prepped recipe. We also stitch the predicted values to the transformed train set with bind_cols function;

rf.predict = rf.fit %>%
  predict(crtr.testing) %>%
  bind_cols(crtr.testing) 

rf.predict

# A tibble: 12 x 8
     .pred temperature      do      ph ammonia phosphate nitrate     chl
     <dbl>       <dbl>   <dbl>   <dbl>   <dbl>     <dbl>   <dbl>   <dbl>
 1  0.0278       0.750  0.884  -0.947    0.258   -0.560   -0.625 -1.50  
 2 -0.0202       0.370  1.44   -0.870    0.291   -0.560   -0.677  0.504 
 3  0.0453      -0.392  0.977  -2.33     0.291   -0.560   -0.521 -0.0237
 4 -0.173       -0.265  0.861  -1.41     0.588    0.0415   1.66  -0.235 
 5 -0.0443      -0.201  1.07   -1.02     1.35     2.01     2.49   0.346 
 6 -0.274       -1.15   0.0720  0.284    0.291   -0.961   -0.573  0.399 
 7 -0.140       -0.265 -0.624   0.515    0.158   -1.12    -0.625 -0.499 
 8  0.0790       0.814 -0.392  -0.0239   2.80    -0.480   -0.625 -0.763 
 9 -0.193        0.814 -1.39   -0.485   -1.10    -0.720   -0.573  0.0819
10 -0.231        1.00  -1.41   -0.101    0.621    0.482    0.206  1.35  
11 -0.417        0.370 -0.856   0.438   -1.10    -0.640   -0.573  1.61  
12  1.30         1.38  -0.949   2.44     0.456    0.362   -0.677 -0.657 

When making predictions, the tidymodels convention is to always produce a tibble of results with standardized column names. This makes it easy to combine the original data and the predictions in a usable format:

Evaluate the rf model

So far, we have built a model and preprocessed data with a recipe. We predicted new data as a way to bundle a parsnip model and recipe together. The next step is to assess or evaluate the accuracy of the model. We use a metrics function from yardstick package (Kuhn and Vaughan 2020b)to assess the accuracy of the model by comparing the predicted versus the original outcome variable. Note that we use the predicted dataset we just computed using predict function.

rf.predict %>%
  metrics(truth = chl, estimate = .pred)

# A tibble: 3 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.10 
2 rsq     standard       0.221
3 mae     standard       0.857

Make linear model

The good of tidymodels is that when we change the model, we do not need to start over again from the beginning but rather change the engine. For instance, we replace the ranger engine with lm to specify the linear model using the parsnip package (Kuhn and Vaughan 2020a) .

lm.model = linear_reg() %>%
  set_engine(engine = "lm") %>%
  set_mode(mode = "regression")

Train Linear model

Once we have specified the model type, engine and mode, the model can be trained with the fit function;

lm.fit = lm.model %>%
  fit(chl ~ ., data = crtr.training)

Predict with linear model

Once the model is fitted, This fitted object lm_fit has the lm model output built-in, which you can access with lm_fit$fit, but there are some benefits to using the fitted parsnip model object when it comes to predicting. To predict the values we use predict function and parse the model and standardized testing values we computed from the recipe (R Core Team 2018). Note that here we use the transformed test set and not the original from the split object. In this case we use the model to predict the value and stitch the testing values using the bind_cols function;

lm.predict = lm.fit %>%
  predict(crtr.testing) %>%
  bind_cols(crtr.testing) 

lm.predict

# A tibble: 12 x 8
     .pred temperature      do      ph ammonia phosphate nitrate     chl
     <dbl>       <dbl>   <dbl>   <dbl>   <dbl>     <dbl>   <dbl>   <dbl>
 1 -0.236        0.750  0.884  -0.947    0.258   -0.560   -0.625 -1.50  
 2 -0.0649       0.370  1.44   -0.870    0.291   -0.560   -0.677  0.504 
 3 -0.728       -0.392  0.977  -2.33     0.291   -0.560   -0.521 -0.0237
 4 -0.764       -0.265  0.861  -1.41     0.588    0.0415   1.66  -0.235 
 5  0.0740      -0.201  1.07   -1.02     1.35     2.01     2.49   0.346 
 6  0.0151      -1.15   0.0720  0.284    0.291   -0.961   -0.573  0.399 
 7 -0.135       -0.265 -0.624   0.515    0.158   -1.12    -0.625 -0.499 
 8  0.323        0.814 -0.392  -0.0239   2.80    -0.480   -0.625 -0.763 
 9 -0.772        0.814 -1.39   -0.485   -1.10    -0.720   -0.573  0.0819
10 -0.117        1.00  -1.41   -0.101    0.621    0.482    0.206  1.35  
11 -0.263        0.370 -0.856   0.438   -1.10    -0.640   -0.573  1.61  
12  1.13         1.38  -0.949   2.44     0.456    0.362   -0.677 -0.657 

Evaluate linear model

Once we have our lm.predict dataset that contains the predicted and test values, we can now use the metrics fiction to evaluate the accuracy of the model.

lm.predict%>%
  metrics(truth = chl, estimate = .pred)

# A tibble: 3 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard      1.07  
2 rsq     standard      0.0337
3 mae     standard      0.929 

Estimate stats

Sometimes you may wish to plot predicted values with different predictors. To do that you need to create a new tidied data from the the model with tidy function and parse interval = TRUE argument as shown in the code below. This create a tibble shown below and the same data is plotted in figure 3.

lm.fit.stats = lm.fit %>% 
  tidy(interval = TRUE)

lm.fit.stats

# A tibble: 7 x 5
  term         estimate std.error statistic p.value
  <chr>           <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)  1.25e-15     0.172  7.25e-15  1.00  
2 temperature -4.82e- 2     0.258 -1.87e- 1  0.854 
3 do           1.83e- 1     0.214  8.58e- 1  0.400 
4 ph           3.88e- 1     0.307  1.26e+ 0  0.220 
5 ammonia      1.59e- 1     0.209  7.64e- 1  0.453 
6 phosphate    3.98e- 1     0.225  1.77e+ 0  0.0913
7 nitrate     -3.00e- 1     0.274 -1.09e+ 0  0.286 

lm.fit.stats %>%
  slice(-1) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(size = 4)+
  geom_errorbar(aes(ymin = estimate-std.error, ymax = estimate+std.error), width = .2)+
  scale_y_continuous(breaks = seq(-1.5,1.5,0.5))+
  ggpubr::theme_pubclean()+
  theme(axis.text = element_text(size = 10))+
  labs(x = "", y = "Estimated chl")

Figure 3: Estimated value of chlorophyll concentration at different predictors

Make Bayes Model

For Bayesian, we also change the engine and specified are called prior and prior_intercept. It turns out that linear_reg() has a stan engine. Since these prior distribution arguments are specific to the Stan software, they are passed as arguments to parsnip::set_engine().

prior.dist = rstanarm::student_t(df = 1)

set.seed(401)

## make the parsnip model
bayes.model = linear_reg() %>%
  set_engine(engine = "stan", 
             prior_intercept = prior.dist, 
             prior = prior.dist) %>%
  set_mode(mode = "regression")

This kind of Bayesian analysis (like many models) involves randomly generated numbers in its fitting procedure. We can use set.seed() to ensure that the same (pseudo-)random numbers are generated each time we run this code. The number 123 isn’t special or related to our data; it is just a “seed” used to choose random numbers.

Train Bayes model

Once we have defined the Bayesian model, we train it using a transformed testing set;

## train the bayes model
bayes.fit = bayes.model%>%
  fit(chl ~ ., data = crtr.testing)

bayes.fit

parsnip model object

Fit time:  1.7s 
stan_glm
 family:       gaussian [identity]
 formula:      chl ~ .
 observations: 12
 predictors:   7
------
            Median MAD_SD
(Intercept)  0.4    0.4  
temperature -0.5    0.6  
do          -0.4    0.4  
ph          -0.2    0.3  
ammonia     -0.3    0.3  
phosphate    0.6    0.8  
nitrate     -0.2    0.6  

Auxiliary parameter(s):
      Median MAD_SD
sigma 1.0    0.3   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

Predict Bayes fit

bayes.predict = bayes.fit %>%
  predict(crtr.testing) %>%
  bind_cols(crtr.testing)

bayes.predict

# A tibble: 12 x 8
     .pred temperature      do      ph ammonia phosphate nitrate     chl
     <dbl>       <dbl>   <dbl>   <dbl>   <dbl>     <dbl>   <dbl>   <dbl>
 1 -0.500        0.750  0.884  -0.947    0.258   -0.560   -0.625 -1.50  
 2 -0.570        0.370  1.44   -0.870    0.291   -0.560   -0.677  0.504 
 3  0.211       -0.392  0.977  -2.33     0.291   -0.560   -0.521 -0.0237
 4 -0.0976      -0.265  0.861  -1.41     0.588    0.0415   1.66  -0.235 
 5  0.475       -0.201  1.07   -1.02     1.35     2.01     2.49   0.346 
 6  0.280       -1.15   0.0720  0.284    0.291   -0.961   -0.573  0.399 
 7  0.0727      -0.265 -0.624   0.515    0.158   -1.12    -0.625 -0.499 
 8 -0.902        0.814 -0.392  -0.0239   2.80    -0.480   -0.625 -0.763 
 9  0.681        0.814 -1.39   -0.485   -1.10    -0.720   -0.573  0.0819
10  0.554        1.00  -1.41   -0.101    0.621    0.482    0.206  1.35  
11  0.552        0.370 -0.856   0.438   -1.10    -0.640   -0.573  1.61  
12 -0.119        1.38  -0.949   2.44     0.456    0.362   -0.677 -0.657 

Evaluate Bayes model

bayes.predict %>%
  metrics(truth = chl, estimate = .pred)

# A tibble: 3 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       0.646
2 rsq     standard       0.415
3 mae     standard       0.533

Bayes.fit.stats

To update the parameter table, the tidy() method is once again used:

bayes.stats = bayes.fit %>% 
  tidy(intervals = TRUE)

bayes.stats

# A tibble: 7 x 5
  term        estimate std.error  lower upper
  <chr>          <dbl>     <dbl>  <dbl> <dbl>
1 (Intercept)    0.373     0.425 -0.394 1.07 
2 temperature   -0.472     0.564 -1.44  0.559
3 do            -0.428     0.416 -1.14  0.340
4 ph            -0.170     0.341 -0.772 0.461
5 ammonia       -0.318     0.312 -0.864 0.244
6 phosphate      0.600     0.757 -0.812 1.89 
7 nitrate       -0.206     0.625 -1.22  0.942

bayes.stats %>% 
  slice(-1) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(size = 4)+
  geom_errorbar(aes(ymin = lower, ymax = upper), width = .2)+
  scale_y_continuous(breaks = seq(-1.5,1.5,0.5))+
  ggpubr::theme_pubclean()+
  theme(axis.text = element_text(size = 10))+
  labs(x = "", y = "Estimated chl")