CRAN availability of tidymodels packages:

package	CRAN	description
broom		Convert statistical analysis objects from R into tidy format
rsample		Classes and functions to create and summarize different types of resampling objects
dials		Tools for creating tuning parameter values
textrecipes		Extra recipes for Text Processing
yardstick		Tidy methods for measuring model performance
parsnip	NR	A tidy unified interface to models
probably	NR	Tools for post-processing class probability estimates
recipes		A preprocessing engine to generate design matrices
embed		Extra recipes for categorical predictor embeddings
infer		An R package for tidyverse-friendly statistical inference
tidyposterior		Bayesian comparisons of models using resampled statistics

Unified Modelling Syntax

The declared goal of the tidymodels metapackage is to provide a unified modelling synthax similar to scikit-learn in the python domain or an improved version of caret but adhering to the tidyverse principles. parsnip is going to be the core package while dials will provide suitable objects and functions for parameter tuning. The amount of supported models is still a bit meager so we will not explore these packages any further for the moment.

Statistical Tests and Model Selection

The regular statistical test supported by R have the same problem as the modelling implementations, they lack a uniform tidyverse compatible synthax. Further traditional statistical tests have lately gotten a bit out of fashion. The following criticism keeps popping up:

• Specific statistical requirements for each test. The strategies for selecting the right statistical tests are a bit convoluted and a certain set of statistical requirements need to be full-filled for each of them.

• Interpretation of P Values. There is a pickiness when it comes to interpreting P Values, the perfect definition eludes me and is completly useless to a none-statistician. Allen Downey has a refreshing practical approach to P values in which he uses a bayesian approach to show that indeed from small p values (<= 0.01) one can conlude that the observed effect has a low probability to be the result of chance (post)

• Disregard of Effect Size. If we have a large sample even irrelevant effects will result in low p-values and if we have a small sample only very large effects will result in low p-values. If we detect a relevant effect with a low p-value we cannot be sure that the magnitude of the effect is reproducible. Typically the effect size will decrease the larger the sample. The Null hypothesis does not incorporate a minimum effect size.

As a remedy for the issue of the convoluted statisical requirements for each test a workaround has again been proposed by Allen Downey. He proposes to simulate data that assumes that there is no connection between two hypothetical sets of data that we want to compare (the null hypothesis is true). ( post1, post2 ). Similar to bootstrapping this method is none-parametric and we can use the simulated data to calculate a set of summary statistics. Then we can compare the distribution of these statistics against the actual value. infer allows us to do just that and on-top offers a tidy synthax to the conventional R implementations of standard statistical tests.

However even the simulation technique does not really help us to judge the effect size properly. This is something that can be adressed using bayesian modelling techniques, which will provide you with a posterior distribution of your response variable which allows you to sufficiently judge the effect size.

When using any k-fold cross-validation strategy for model training and validation we can apply statistical tests on each set of k performance metrics to select the best performing model. In general we run into the same issues as discussed above. In order to adress them we can either use the simulation technique of the inferpackage or use tidyposterior which uses Bayesian modelling to compare performance metrics which allows us to define a relevant effect size to test against.

In general I think tidyposterior is probably best practise, however to reduce complexity I am personally quite happy with the 1 SE rule. Simply plotting the mean value with the SE and then picking the simplest model that is within 1SE of the model with the highest performance. Thus I will not include these packages in my modelling workflow for the moment.

Resampling, Feature Engineering and Performance Metrics

rsample, recipes and yardstick are packages that give an overall complete impression and can be used with caret. rsample allows us to create cross validation pairs by indexing an existing dataframe and facilitates the use of modelling dataframes. If supports a variety of resampling methods such as not only limited to k-fold cross validation but also bootstrapping and nested cross validation. recipes allows straight forward feature engineering and preprocessing and yardstick allows us to easily calculate performance metrics from model predictions.

Data

data('BostonHousing')
df = as_tibble( BostonHousing )
summary(df)

##       crim                zn             indus       chas         nox        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   0:471   Min.   :0.3850  
##  1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   1: 35   1st Qu.:0.4490  
##  Median : 0.25651   Median :  0.00   Median : 9.69           Median :0.5380  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14           Mean   :0.5547  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10           3rd Qu.:0.6240  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74           Max.   :0.8710  
##        rm             age              dis              rad        
##  Min.   :3.561   Min.   :  2.90   Min.   : 1.130   Min.   : 1.000  
##  1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100   1st Qu.: 4.000  
##  Median :6.208   Median : 77.50   Median : 3.207   Median : 5.000  
##  Mean   :6.285   Mean   : 68.57   Mean   : 3.795   Mean   : 9.549  
##  3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188   3rd Qu.:24.000  
##  Max.   :8.780   Max.   :100.00   Max.   :12.127   Max.   :24.000  
##       tax           ptratio            b              lstat      
##  Min.   :187.0   Min.   :12.60   Min.   :  0.32   Min.   : 1.73  
##  1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38   1st Qu.: 6.95  
##  Median :330.0   Median :19.05   Median :391.44   Median :11.36  
##  Mean   :408.2   Mean   :18.46   Mean   :356.67   Mean   :12.65  
##  3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23   3rd Qu.:16.95  
##  Max.   :711.0   Max.   :22.00   Max.   :396.90   Max.   :37.97  
##       medv      
##  Min.   : 5.00  
##  1st Qu.:17.02  
##  Median :21.20  
##  Mean   :22.53  
##  3rd Qu.:25.00  
##  Max.   :50.00

p_hist = ggplot(df) +
  geom_histogram( aes(lstat) ) +
  lims( x = c(0,40) )

p_ecdf = ggplot(df) +
  stat_ecdf(aes(lstat) ) +
  lims( x = c(0,40) )

gridExtra::grid.arrange( p_hist, p_ecdf )

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

df_cor = df %>%
  select_if( is.numeric ) %>%
  gather( key = 'variable', value = 'value', - lstat) %>%
  group_by(variable) %>%
  nest() %>%
  mutate( cor = map_dbl(data, function(x) cor(x$lstat, x$value) ) ) %>%
  unnest(cols = c(data)) %>%
  ungroup() %>%
  mutate( variable = fct_reorder(variable, cor)
          , cor = round(cor,2) )

df_label = df_cor %>%
  group_by( variable, cor) %>%
  summarise( pos = max(value) *.9 )

ggplot( df_cor, aes(lstat, value) ) +
  geom_point( alpha = 0.2 ) +
  geom_smooth( method = 'lm') +
  geom_label( aes( x = 5, y = pos, label = cor)
             , df_label
             , color = 'pink') +
  facet_wrap(~variable, scales = 'free_y')

## `geom_smooth()` using formula 'y ~ x'

df %>%
  select_if( is.factor ) %>%
  bind_cols( df['lstat'] ) %>%
  gather( key = 'variable', value = 'value', - lstat) %>%
  ggplot( aes( x = value, y = lstat) ) +
  geom_violin() +
  geom_boxplot( alpha = 0.5 ) +
  ggpubr::stat_compare_means() +
  facet_wrap( ~ variable )

Preprocessing with recipe

Note we are intentionally standardizing the response variable since the unit of lstat is irrelevant for this demo

We will

• Yeo Johnson Transform
• Scale
• Center
• remove co-correlating variables (threshold 0.5)
• dummy encode

rec = recipe(df, lstat ~ . )

summary(rec)

## # A tibble: 14 x 4
##    variable type    role      source  
##    <chr>    <chr>   <chr>     <chr>   
##  1 crim     numeric predictor original
##  2 zn       numeric predictor original
##  3 indus    numeric predictor original
##  4 chas     nominal predictor original
##  5 nox      numeric predictor original
##  6 rm       numeric predictor original
##  7 age      numeric predictor original
##  8 dis      numeric predictor original
##  9 rad      numeric predictor original
## 10 tax      numeric predictor original
## 11 ptratio  numeric predictor original
## 12 b        numeric predictor original
## 13 medv     numeric predictor original
## 14 lstat    numeric outcome   original

rec = rec %>%
  step_scale( all_numeric() ) %>%
  step_center( all_numeric() ) %>%
  step_YeoJohnson( all_numeric() ) %>%
  step_corr( all_numeric(), - all_outcomes(), threshold = 0.5 ) %>%
  step_dummy( all_nominal() )

prep_rec = prep(rec, df)
prep_rec

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         13
## 
## Training data contained 506 data points and no missing data.
## 
## Operations:
## 
## Scaling for crim, zn, indus, nox, rm, age, dis, rad, ... [trained]
## Centering for crim, zn, indus, nox, rm, age, dis, rad, ... [trained]
## Yeo-Johnson transformation on crim, zn, indus, nox, rm, age, ... [trained]
## Correlation filter removed indus, nox, age, dis, tax, medv, crim [trained]
## Dummy variables from chas [trained]

df_prep = bake(prep_rec, df )

Resampling with rsample

rs = rsample::vfold_cv(df, v = 5)

rsample::pretty.vfold_cv(rs)

## [1] "5-fold cross-validation"

Convert to caret-compatible object

rs_caret = rsample::rsample2caret(rs)

Modelling with caret

car = function( method, recipe, rsample, data){
  

  car = caret::train( recipe
                      , data
                      , method = method
                      , trControl = caret::trainControl(index = rsample$index
                                                        , indexOut = rsample$indexOut
                                                        , method = 'cv'
                                                        , verboseIter = T
                                                        , savePredictions = T
                                                        , search = 'random')
                      , metric = 'RMSE'
                      , tuneLength = 100
                       )
  
  return( car )
}

# c = car( 'lm', rec, rs_caret, df)

df_m = tibble( methods = c('lm', 'rpart', 'cubist', 'parRF', 'earth', 'xgbTree') )

df_m = df_m %>%
  mutate( c = map(methods, car, rec, rs_caret, df ) )

Assess Performance with yardstick

df_pred = df_m %>%
  mutate( pred = map(c, 'pred' )
          , pred = map(pred, as_tibble )
          , best_tune = map(c, 'bestTune') )
  
df_pred  

## # A tibble: 6 x 4
##   methods c          pred                   best_tune       
##   <chr>   <list>     <list>                 <list>          
## 1 lm      <tran.rcp> <tibble [506 × 5]>     <df[,1] [1 × 1]>
## 2 rpart   <tran.rcp> <tibble [18,216 × 5]>  <df[,1] [1 × 1]>
## 3 cubist  <tran.rcp> <tibble [48,070 × 6]>  <df[,2] [1 × 2]>
## 4 parRF   <tran.rcp> <tibble [3,036 × 5]>   <df[,1] [1 × 1]>
## 5 earth   <tran.rcp> <tibble [13,156 × 6]>  <df[,2] [1 × 2]>
## 6 xgbTree <tran.rcp> <tibble [50,600 × 11]> <df[,7] [1 × 7]>

filter(df_pred, methods == 'cubist') %>%
  .$pred 

## [[1]]
## # A tibble: 48,070 x 6
##        obs rowIndex    pred committees neighbors Resample
##      <dbl>    <int>   <dbl>      <int>     <int> <chr>   
##  1 -1.87          4 -1.32            3         5 Fold1   
##  2  0.740         8 -0.0712          3         5 Fold1   
##  3  1.55          9  0.616           3         5 Fold1   
##  4  0.855        11 -0.240           3         5 Fold1   
##  5 -0.351        15 -0.0309          3         5 Fold1   
##  6 -0.123        19  1.15            3         5 Fold1   
##  7  0.481        26  1.06            3         5 Fold1   
##  8  0.292        27  0.225           3         5 Fold1   
##  9  0.560        28  0.608           3         5 Fold1   
## 10 -0.0793       30 -0.0862          3         5 Fold1   
## # … with 48,060 more rows

params_as_str = function(df, params){
  
  symbols = map( names(params), as.name )
  
  df %>%
    mutate( desc_values = pmap_chr( list( !!! symbols), paste )
            , desc_params = paste( names(params), collapse = ' ' ) )
}

# params_as_str(df_pred$pred[[6]], df_pred$best_tune[[6]] )

df_pred = df_pred %>%
  mutate( pred = map2(pred, best_tune, params_as_str )
          , pred = map(pred, select, Resample, desc_params, desc_values, rowIndex, obs, pred)
          ) %>%
  unnest(pred)

df_pred

## # A tibble: 133,584 x 9
##    methods c     Resample desc_params desc_values rowIndex     obs    pred
##    <chr>   <lis> <chr>    <chr>       <chr>          <int>   <dbl>   <dbl>
##  1 lm      <tra… Fold1    intercept   TRUE               4 -1.87   -0.815 
##  2 lm      <tra… Fold1    intercept   TRUE               8  0.740  -0.627 
##  3 lm      <tra… Fold1    intercept   TRUE               9  1.55   -0.177 
##  4 lm      <tra… Fold1    intercept   TRUE              11  0.855  -0.726 
##  5 lm      <tra… Fold1    intercept   TRUE              15 -0.351   0.0395
##  6 lm      <tra… Fold1    intercept   TRUE              19 -0.123   0.644 
##  7 lm      <tra… Fold1    intercept   TRUE              26  0.481   0.534 
##  8 lm      <tra… Fold1    intercept   TRUE              27  0.292   0.250 
##  9 lm      <tra… Fold1    intercept   TRUE              28  0.560   0.218 
## 10 lm      <tra… Fold1    intercept   TRUE              30 -0.0793 -0.331 
## # … with 133,574 more rows, and 1 more variable: best_tune <list>

df_best_models = df_pred %>%
  group_by( methods, desc_params, desc_values) %>%
  yardstick::rmse(obs, pred) %>%
  group_by( methods ) %>%
  mutate( rnk = rank(.estimate, ties.method = 'first' ) ) %>%
  filter( rnk == 1 ) %>%
  select( - rnk ) %>%
  arrange(.estimate) %>%
  ungroup() %>%
  mutate( methods = fct_reorder(methods, .estimate) )

df_best_models

## # A tibble: 6 x 6
##   methods desc_params          desc_values          .metric .estimator .estimate
##   <fct>   <chr>                <chr>                <chr>   <chr>          <dbl>
## 1 parRF   mtry                 2                    rmse    standard       0.512
## 2 xgbTree nrounds max_depth e… 638 4 0.02251800388… rmse    standard       0.540
## 3 cubist  committees neighbors 36 8                 rmse    standard       0.540
## 4 earth   nprune degree        8 2                  rmse    standard       0.563
## 5 rpart   cp                   0.00450606062062921  rmse    standard       0.592
## 6 lm      intercept            TRUE                 rmse    standard       0.624

performance = yardstick::metric_set( yardstick::rmse, yardstick::rsq, yardstick::mae, yardstick::mape )

df_perf = df_best_models %>%
  select(methods, desc_params, desc_values) %>%
  left_join(df_pred) %>%
  group_by( methods, Resample) %>%
  performance(obs, pred) %>%
  mutate( methods = as.factor(methods)
          , methods = fct_relevel(methods, levels(df_best_models$methods) )) %>%
  group_by(methods, .metric) %>%
  mutate( me = mean(.estimate)
             , se = sd(.estimate)/sqrt(n()) )

## Joining, by = c("methods", "desc_params", "desc_values")

df_1se = df_perf %>%
  group_by(methods, .metric, me, se) %>%
  summarise() %>%
  mutate(  ymin = me - se
          , ymax = me + se ) %>%
  group_by(.metric) %>%
  mutate( rnk = rank(me, ties.method = 'first')
          , rnk_desc = rank( desc(me), ties.method = 'first')
          ) %>%
  rename( best_method = methods ) %>%
  filter( (rnk == 1 & .metric != 'rsq') | (.metric == 'rsq' & rnk_desc == 1) )

df_1se

## # A tibble: 4 x 8
## # Groups:   .metric [4]
##   best_method .metric      me      se    ymin    ymax   rnk rnk_desc
##   <fct>       <chr>     <dbl>   <dbl>   <dbl>   <dbl> <int>    <int>
## 1 parRF       mae       0.394  0.0149   0.379   0.409     1        6
## 2 parRF       rmse      0.510  0.0221   0.488   0.532     1        6
## 3 parRF       rsq       0.703  0.0372   0.666   0.740     6        1
## 4 lm          mape    167.    32.0    135.    199.        1        6

len = levels(df_perf$methods) %>%
  length()

col_class = RColorBrewer::brewer.pal('Greys', n = 9) %>% rev()
col_folds = RColorBrewer::brewer.pal('Dark2', n = 8) %>%
  easyalluvial::palette_filter(greens = F, greys = F)

pal = c( col_class[2], col_folds[1:5], col_class[4], col_class[6] )

df_perf %>%
  left_join( select(df_1se, .metric, ymin, ymax, best_method) ) %>%
  mutate( classification = case_when( best_method == methods ~ 'best'
                                      , me >= ymin & me <= ymax ~ 'in'
                                      , T ~ 'out' )
          ) %>%
  ggplot( aes(methods, .estimate) ) +
    geom_rect( aes(ymin = ymin, ymax = ymax )
               , xmin = 0, xmax = len
               , fill = col_class[7]
               , alpha = 0.05 ) +
    geom_line( aes( group = Resample, color = Resample)
               , size = .5, alpha = 1 ) +
    stat_summary( aes( color = classification)
                  , size = 1
                  , fun.data = function(x) mean_se( x, mult = 1) ) +
    scale_color_manual( values = pal ) +
    theme( legend.position = 'none') +
    labs( y = '', x = '', caption = 'grey area: 1SE range of best model'
        , title = 'CV Performance Metrics + SE') +
    facet_wrap(~.metric, scales = 'free_y', ncol = 1) 

## Joining, by = ".metric"