Ordinal data models

• Introduction
• Data preparation
• ordered logistic regression (logit)
• Ordinal logistic rgeression (probit)
• CART model
• Ordinal Random forst model.
• Continuation Ratio Model
• Compare models
• Conclusion
• Session information

Introduction

This tutorial aims to explore the most popular models used to predict an ordered response variable. We will use the heart disease data uploaded from kaggle website, where our response will be the chest pain cp variable instead of the target variable used usually.

Data preparation

First, we call the data and the libraries that we need along with this illustration as follows.

library(tidyverse)
library(caret)
library(tidymodels)
mydata<-read.csv("../heart.csv",header = TRUE)
names(mydata)[1]<-"age"

The data at hand has the following features:

• age.
• sex: 1=male,0=female
• cp : chest pain type.
• trestbps : resting blood pressure.
• chol: serum cholestoral.
• fbs : fasting blood sugar.
• restecg : resting electrocardiographic results.
• thalach : maximum heart rate achieved
• exang : exercise induced angina.
• oldpeak : ST depression induced by exercise relative to rest.
• slope : the slope of the peak exercise ST segment.
• ca : number of major vessels colored by flourosopy.
• thal : it is not well defined from the data source.
• target: have heart disease or not.

I think the best start to explore the summary of all predictors and missing values is by using the powerful function skim from skimr package.

skimr::skim(mydata)

Variable type: numeric

For our case we will use the chest pain type cp variable as our target variable since it is a categorical variable. However, for pedagogic purposes, we will manipulate it so that it will be an ordered factor with only three levels no pain,moderate pain, severe pain (instead of 4 levels now).

Looking at the above output, we convert the variables that should be of factor type, which are: sex, target, fbs, resecg, exang, slope, ca, thal. For the response variable cp, we drop its less frequently level with all its related rows, then we rename the remaining ones as no pain for the most frequently one, severe pain for the less frequently one, and moderate pain for the last one.

table(mydata$cp)

  0   1   2   3 
143  50  87  23 

we see the level 3 is the less frequently one.

mydata<-mydata %>%
  modify_at(c("cp", "sex", "target", "fbs", "resecg", "exang", "slope", "ca", "thal"),
            as.factor)
mydata<-mydata[mydata$cp!=3,]
mydata$cp<-fct_drop(mydata$cp,only=levels(mydata$cp))
table(mydata$cp)

  0   1   2 
143  50  87 

According to these frequencies we rename and we order the levels as follows.

mydata$cp<-fct_recode(mydata$cp,no="0",sev="1",mod="2")
mydata$cp<-factor(mydata$cp,ordered = TRUE)
mydata$cp<-fct_infreq(mydata$cp)
mydata$cp[1:5]
[1] mod sev sev no  no 
Levels: no < mod < sev

Similar to the logistic regression, the number of cases in each cell from each cross table between the outcome and each factor should be above the threshold of 5 applied in practice.

xtabs(~cp+sex,data=mydata)
     sex
cp      0   1
  no   39 104
  mod  35  52
  sev  18  32
xtabs(~cp+target,data=mydata)
     target
cp      0   1
  no  104  39
  mod  18  69
  sev   9  41
xtabs(~cp+fbs,data=mydata)
     fbs
cp      0   1
  no  125  18
  mod  70  17
  sev  45   5
xtabs(~cp+restecg,data=mydata)
     restecg
cp     0  1  2
  no  78 62  3
  mod 36 50  1
  sev 19 31  0
xtabs(~cp+exang,data=mydata)
     exang
cp     0  1
  no  63 80
  mod 76 11
  sev 46  4
xtabs(~cp+slope,data=mydata)
     slope
cp     0  1  2
  no  11 84 48
  mod  5 33 49
  sev  2 12 36
xtabs(~cp+ca,data=mydata)
     ca
cp     0  1  2  3  4
  no  65 34 29 14  1
  mod 57 20  2  5  3
  sev 37  8  3  1  1
xtabs(~cp+thal,data=mydata)
     thal
cp     0  1  2  3
  no   1 12 52 78
  mod  1  2 62 22
  sev  0  2 39  9

The following variables do not respect this threshold and hence they will be removed from the predictors set: restecg, exang, slope, ca, and thal.

mydata<-mydata[,setdiff(names(mydata), c("restecg", "exang", "slope", "ca",  "thal"))]

The data is ready and we can now split the data between training and testing set.

set.seed(1122)
parts <- initial_split(mydata, prop=0.8, strata = cp)
train <- training(parts)
test <- testing(parts)

The most popular models that we will use are: ordinal logistic model, cart model, ordinal random forest model, Continuation ratio model.

ordered logistic regression (logit)

Before training this type of model let’s show how it works. For simplicity suppose we have data that has an ordered outcome \(y\) with three class labels (“1”,“2”,“3”) and only two features \(x_1\) and \(x_2\).

First we define a latent variable as a linear combination of the features:

\[y_i^*=\beta_1 X_{i1}+\beta_2 X_{i2}\]

Then since we have three classes we define two thresholds for this latent variable \(\alpha_1\) and \(\alpha_2\) such that a particular observation \(y_i\) will be classified as follows:

\[\begin{cases} y_i=1 & \text{if $y_i^* \leq \alpha_1$} \\ y_i=2 & \text{if $\alpha_1 < y_i^* \leq \alpha_2$} \\ y_i=3 & \text{if $y_i^* > \alpha_2$}\end{cases}\]

Now we can obtain the probability of a particular observation to fall into a specific class as follows:

\[\begin{cases} p(y_i=1)=p(y_i^* \leq \alpha_1)=F(\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2}) \\ p(y_i=2)=p(\alpha_1 < y_i^* \leq \alpha_2)=F(\alpha_2-\beta_1 X_{i1}-\beta_2 X_{i2})-F(\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2}) \\ p(y_i=3)=1-p(y_i=2)-p(y_i=1)\end{cases}\]

It remains now to define the suitable distribution function F. There are two commonly used ones for this type of data, the logit function \(F(x)=\frac{1}{1+exp^{-x}}\) and the normal distribution function aka probit.

Note: there exist other functions like loglog, cloglog, and cauchit.

Using the logit function the probabilities will be.

\[\begin{cases} p(y_i=1)=\frac{1}{1+exp^{-(\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2})}} \\ p(y_i=2)=\frac{1}{1+exp^{-(\alpha_2-\beta_1 X_{i1}-\beta_2 X_{i2})}}-p(y_i=1) \\ p(y_i=3)=1-p(y_i=2)-p(y_i=1)\end{cases}\]

The MASS package provides the method polr to perform an ordinal logistic regression.

library(MASS)
set.seed(1234)
model_logistic<-train(cp~., data=train,
                      method="polr",
                      tuneGrid=expand.grid(method="logistic"))

summary(model_logistic)

Coefficients:
              Value Std. Error  t value
age       0.0112236   0.018219  0.61605
sex1      0.2593720   0.316333  0.81993
trestbps -0.0002329   0.009090 -0.02562
chol     -0.0013238   0.002697 -0.49082
fbs1      0.3188826   0.401836  0.79356
thalach   0.0226246   0.008199  2.75933
oldpeak  -0.3360326   0.163547 -2.05465
target1   1.7234740   0.376279  4.58031

Intercepts:
        Value   Std. Error t value
no|mod   4.5786  1.9271     2.3759
mod|sev  6.5004  1.9527     3.3289

Residual Deviance: 376.4697 
AIC: 396.4697 

This table does not provide the p-values. However, it is not a big problem since we can add the p_values by the following script.

prob <- pnorm(abs(summary(model_logistic)$coefficients[,3]),lower.tail = FALSE)*2
cbind(summary(model_logistic)$coefficients,prob)
                 Value  Std. Error     t value         prob
age       0.0112236479 0.018218848  0.61604597 5.378642e-01
sex1      0.2593719567 0.316332564  0.81993442 4.122535e-01
trestbps -0.0002329023 0.009090066 -0.02562163 9.795591e-01
chol     -0.0013237835 0.002697079 -0.49082122 6.235529e-01
fbs1      0.3188825831 0.401836034  0.79356393 4.274493e-01
thalach   0.0226246089 0.008199317  2.75932853 5.792027e-03
oldpeak  -0.3360326371 0.163547467 -2.05464899 3.991292e-02
target1   1.7234739863 0.376278770  4.58031152 4.642839e-06
no|mod    4.5785821473 1.927119568  2.37586822 1.750771e-02
mod|sev   6.5003986218 1.952726089  3.32888399 8.719471e-04

Using the threshold p-value 0.05, we remove the non significant variables. age, trestbps, chol.

set.seed(1234)
model_logistic<-train(cp~.-age-trestbps-chol, data=train,
                      method="polr",tuneGrid=expand.grid(method="logistic"))
prob <- pnorm(abs(summary(model_logistic)$coefficients[,3]),lower.tail = FALSE)*2
cbind(summary(model_logistic)$coefficients,prob)
              Value  Std. Error    t value         prob
sex1     0.25427581 0.308143065  0.8251875 4.092651e-01
fbs1     0.37177505 0.384667871  0.9664832 3.338024e-01
thalach  0.02050951 0.007487511  2.7391620 6.159602e-03
oldpeak -0.33669473 0.161699555 -2.0822242 3.732199e-02
target1  1.71338020 0.369558584  4.6362885 3.547208e-06
no|mod   4.00836398 1.143111953  3.5065367 4.539789e-04
mod|sev  5.92987585 1.185074388  5.0038005 5.621092e-07

Notice that we do not remove the factors sex and fbs even they are not significant due to the significance of the intercepts.

To well understand these coefficients lets restrict the model with only two predictors.

set.seed(1234)
model1<-train(cp~target+thalach, 
              data=train,
              method = "polr",
              tuneGrid=expand.grid(method="logistic"))
summary(model1)

Coefficients:
          Value Std. Error t value
target1 1.87953   0.333153   5.642
thalach 0.02347   0.007372   3.184

Intercepts:
        Value  Std. Error t value
no|mod  4.6457 1.0799     4.3018 
mod|sev 6.5325 1.1271     5.7959 

Residual Deviance: 383.3144 
AIC: 391.3144 

Let’s plug in these coefficients in the above equations we obtain the probability of each class as follows:

\[\begin{cases} p(no)=\frac{1}{1+exp^{-(4.6457-1.87953X_{i1}-0.02347X_{i2})}} \\ p(mod)=\frac{1}{1+exp^{-(6.5325-1.87953X_{i1}-0.02347X_{i2})}}-p(no) \\ p(sev)=1-p(mod)-p(no)\end{cases}\]

Let’s now predict a particular patient, say the third one.

train[3,c("cp","thalach","target")]
   cp thalach target
4 sev     178      1

We plug in the predictor values as follows:

\[\begin{cases} p(no)=\frac{1}{1+exp^{-(4.6457-1.87953*1-0.02347*178)}} \\ p(mod)=\frac{1}{1+exp^{-(6.5325-1.87953*1-0.02347*178)}}-p(no) \\ p(sev)=1-p(mod)-p(no)\end{cases}=\begin{cases} p(no)=0.1959992 \\ p(mod)=0.6166398-0.1959992=0.4206406 \\ p(sev)=1-0.4206406-0.1959992=0.3833602\end{cases}\]

Using the highest probability this patient will be predicted to have mod pain. Now let’s compare these probabilities with those obtained from function predict.

predict(model1, train[1:3,], type = "prob") %>% tail(1)
         no       mod      sev
4 0.1958709 0.4205981 0.383531

Now we go back to our original model and compute the accuracy rate for the training data.

predict(model_logistic, train) %>% 
  bind_cols(train) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.611

with the logistic regression model we get 61% accuracy for the training set, which is quite bad. So let’s test the model using the testing set now.

predict(model_logistic, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.648

Surprisingly, the accuracy rate for the testing set is about 65%, which is larger than that computed from the training data (61%). This is an indication of an underfitting problem (The opposite effect of overfitting problem). Is there any way to improve the model performance? Maybe yes, by going back and tune some hyperparameters, but since we have an underfitting problem we do not have much hyperparameters for this model except the type of function used which is by default the logistic function, but there exist as well other functions like probit, loglog, …ect.

For our case let’s try this model with the probit function

Ordinal logistic rgeression (probit)

set.seed(1234)
model_probit<-train(cp~.-age-trestbps-chol, data=train,                                        method="polr",
                    tuneGrid=expand.grid(method="probit"))

predict(model_probit, train) %>% 
  bind_cols(train) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.606

This rate is slightly worse than that from the previous model. But what about the testing set.

predict(model_probit, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.593

This one also is worse than the previous model. So this means that the logistic function for this data performs better than the probit one.

When we try many things to improve the model performance and we do not gain much, it will be better to think to try different types of models.

CART model

This is a tree-based model. To train this model we make use of rpartScore package, and for simplification, we will include only the significant predictors from the previous model.

library(rpartScore)
set.seed(1234)
model_cart<-train(cp~.-age-trestbps-chol, data=train,
                      method="rpartScore")
model_cart
CART or Ordinal Responses 

226 samples
  8 predictor
  3 classes: 'no', 'mod', 'sev' 

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 226, 226, 226, 226, 226, 226, ... 
Resampling results across tuning parameters:

  cp          split  prune  Accuracy   Kappa    
  0.02702703  abs    mr     0.5748197  0.2845545
  0.02702703  abs    mc     0.5796085  0.3011122
  0.02702703  quad   mr     0.5711605  0.2764466
  0.02702703  quad   mc     0.5805216  0.3020125
  0.04504505  abs    mr     0.5620975  0.2719646
  0.04504505  abs    mc     0.5966801  0.3274893
  0.04504505  quad   mr     0.5592845  0.2608402
  0.04504505  quad   mc     0.5930817  0.3208220
  0.21621622  abs    mr     0.5303342  0.1266324
  0.21621622  abs    mc     0.6004116  0.3343997
  0.21621622  quad   mr     0.5290009  0.1143360
  0.21621622  quad   mc     0.5928132  0.3225686

Accuracy was used to select the optimal model using the largest value.
The final values used for the model were cp = 0.2162162, split = abs and
 prune = mc.

The caret model uses the bootstrapping technique for hyperparameters tuning. In our case, the largest accuracy rate is about 59.59%, with the complexity parameter **cp**=0.2162162, the **split**=abs, and **prune**= **mc**.

The argument split controls the splitting function used to grow the tree by setting the misclassification costs in the generalized Gini impurity function to the absolute abs or squared quad. The argument prune is used to select the performance measure to prune the tree between total misclassification rate mr or misclassification cost mc.

predict(model_cart, train) %>% 
  bind_cols(train) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.615

Surprisingly, we get approximately the same accuracy rate as the logit model. Let’s check the testing set.

predict(model_cart, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.630

Now wit this model we get a lower accuracy rate than that of the logistic model.

Ordinal Random forst model.

This model is a corrected version of random forest model that takes into account the ordinal nature of the response variable. For more detail about this model read this great paper.

To train ordinal random forest model, we need to call the following packages: e1071, ranger, ordinalForest.

library(ordinalForest)
library(ranger)
library(e1071)

Since the create function train use bootstrapping method to perform hyperparameters tuning to choose the best values, this makes the training process very slow, that is why i save the resulted output and load it again

# set.seed(1234)
# model_forest<-train(cp~.-age-trestbps-chol, data=train,
#                       method='ordinalRF')

# saveRDS(model_forest, "C://Users/dell/Documents/new-blog/content/post/ordinal/model_forest.rds")

model_forest <- readRDS("C://Users/dell/Documents/new-blog/content/post/ordinal/model_forest.rds")

model_forest
Random Forest 

226 samples
  8 predictor
  3 classes: 'no', 'mod', 'sev' 

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 226, 226, 226, 226, 226, 226, ... 
Resampling results across tuning parameters:

  nsets  ntreeperdiv  ntreefinal  Accuracy   Kappa    
   50     50          200         0.5878619  0.3173073
   50     50          400         0.5862998  0.3174139
   50     50          600         0.5793897  0.3052328
   50    100          200         0.5856573  0.3157938
   50    100          400         0.5865246  0.3155563
   50    100          600         0.5908431  0.3248941
   50    150          200         0.5929076  0.3277110
   50    150          400         0.5822072  0.3092856
   50    150          600         0.5878854  0.3209507
  100     50          200         0.5903104  0.3277174
  100     50          400         0.5817918  0.3127438
  100     50          600         0.5863685  0.3197145
  100    100          200         0.5841147  0.3136250
  100    100          400         0.5830929  0.3115348
  100    100          600         0.5875075  0.3228354
  100    150          200         0.5837402  0.3167977
  100    150          400         0.5841116  0.3187849
  100    150          600         0.5888038  0.3250622
  150     50          200         0.5872276  0.3193143
  150     50          400         0.5873943  0.3198634
  150     50          600         0.5874236  0.3224713
  150    100          200         0.5870714  0.3218660
  150    100          400         0.5805275  0.3107488
  150    100          600         0.5816353  0.3135750
  150    150          200         0.5884693  0.3247512
  150    150          400         0.5830154  0.3154980
  150    150          600         0.5885156  0.3253794

Accuracy was used to select the optimal model using the largest value.
The final values used for the model were nsets = 50, ntreeperdiv = 150
 and ntreefinal = 200.

We can plot the important predictors as follows.

plot(varImp(model_forest))

Now we can obtain the accuracy rate for the training rate as follows.

predict(model_forest, train) %>% 
  bind_cols(train) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.845

Great!, with this model, the accuracy rate has largely improved to roughly 84%. But wait, what matters is the accuracy of the testing set.

predict(model_forest, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.519

This is exactly what is called the overfitting problem. The model generalizes poorly to new unseen data. We can go back and tune some other hyperparameters like increasing the minimum size of nodes (default is 5) to fight the overfitting problem. we do not, however, do that here since it is not the purpose of this tutorial.

Continuation Ratio Model

This model uses The vector generalized additive models which are available in the VGAM package. for more detail about these models click here.

library(VGAM)
set.seed(1234)
model_vgam<-train(cp~.-age-trestbps-chol, data=train,
                  method="vglmContRatio", trace=FALSE)
model_vgam
Continuation Ratio Model for Ordinal Data 

226 samples
  8 predictor
  3 classes: 'no', 'mod', 'sev' 

No pre-processing
Resampling: Bootstrapped (25 reps) 
Summary of sample sizes: 226, 226, 226, 226, 226, 226, ... 
Resampling results across tuning parameters:

  parallel  link     Accuracy   Kappa    
  FALSE     logit    0.5962581  0.3323075
  FALSE     probit   0.5942637  0.3302998
  FALSE     cloglog  0.5973844  0.3293056
  FALSE     cauchit  0.5967368  0.3316896
  FALSE     logc     0.5945121  0.3152759
   TRUE     logit    0.5758330  0.2961673
   TRUE     probit   0.5738297  0.2924747
   TRUE     cloglog  0.5838764  0.3014038
   TRUE     cauchit  0.5810184  0.3067004
   TRUE     logc     0.5302522  0.1031624

Accuracy was used to select the optimal model using the largest value.
The final values used for the model were parallel = FALSE and link = cloglog.

the best model is obtained when the argument parallel is FALSE and link is cauchit which is the tangent function.

The accuracy rate of the training data is:

predict(model_vgam, train) %>% 
  bind_cols(train) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.659

And the accuracy of the testing set is:

predict(model_vgam, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth)
# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy multiclass     0.630

This the best accuracy rate compared to the other models.

Compare models

We can compare between the above models using resample caret function.

models_eval<-resamples(list(logit=model_logistic,
                            cart=model_cart,
                            forest=model_forest,
                            vgam=model_vgam))
summary(models_eval)

Call:
summary.resamples(object = models_eval)

Models: logit, cart, forest, vgam 
Number of resamples: 25 

Accuracy 
            Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
logit  0.5060241 0.5731707 0.5822785 0.5871083 0.6097561 0.6627907    0
cart   0.3734940 0.5824176 0.6097561 0.6004116 0.6279070 0.6746988    0
forest 0.5357143 0.5569620 0.5853659 0.5929076 0.6117647 0.6860465    0
vgam   0.4936709 0.5760870 0.6046512 0.5973844 0.6202532 0.6626506    0

Kappa 
               Min.   1st Qu.    Median      Mean   3rd Qu.      Max. NA's
logit   0.189086980 0.2792369 0.3144822 0.3100458 0.3437500 0.4512651    0
cart   -0.004889406 0.3185420 0.3474144 0.3343997 0.3775576 0.4526136    0
forest  0.233434988 0.2852665 0.3289346 0.3277110 0.3554862 0.4618772    0
vgam    0.144558744 0.2993406 0.3367647 0.3293056 0.3690791 0.4142980    0

Based on the training set and using the mean of the accuracy rate we can say that cart model is the best model for this data with 60.97% accuracy for the training set. However, things are different when it comes to use the testing set instead.

tibble(models=c("logit", "cart", "forest", "vgam"), 
       accuracy=c(
  predict(model_logistic, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth) %>% 
  .[, ".estimate"],
  predict(model_cart, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth) %>% 
  .[, ".estimate"],
  predict(model_forest, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth) %>% 
  .[, ".estimate"],
  predict(model_vgam, test) %>% 
  bind_cols(test) %>%
  rename(pred="...1", truth=cp) %>% 
  accuracy(pred, truth) %>% 
  .[, ".estimate"])) %>% 
  unnest()
# A tibble: 4 x 2
  models accuracy
  <chr>     <dbl>
1 logit     0.648
2 cart      0.630
3 forest    0.519
4 vgam      0.630

Using the testing set, the logistic model with the link logit is the best model to predict this data.

Conclusion

We have seen so far how to model ordinal data by exploring several models, and it happened that the logistic model is the best on for our data. However, in general the best model depends strongly on the data at hand.

Session information

sessionInfo()
R version 4.0.1 (2020-06-06)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] splines   stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] VGAM_1.1-3          e1071_1.7-3         ranger_0.12.1      
 [4] ordinalForest_2.3-1 rpartScore_1.0-1    rpart_4.1-15       
 [7] MASS_7.3-51.6       yardstick_0.0.6     workflows_0.1.1    
[10] tune_0.1.0          rsample_0.0.7       recipes_0.1.12     
[13] parsnip_0.1.1       infer_0.5.2         dials_0.0.7        
[16] scales_1.1.1        broom_0.5.6         tidymodels_0.1.0   
[19] caret_6.0-86        lattice_0.20-41     forcats_0.5.0      
[22] stringr_1.4.0       dplyr_1.0.0         purrr_0.3.4        
[25] readr_1.3.1         tidyr_1.1.0         tibble_3.0.1       
[28] ggplot2_3.3.1       tidyverse_1.3.0    

loaded via a namespace (and not attached):
  [1] readxl_1.3.1         backports_1.1.7      tidytext_0.2.4      
  [4] plyr_1.8.6           igraph_1.2.5         repr_1.1.0          
  [7] crosstalk_1.1.0.1    listenv_0.8.0        SnowballC_0.7.0     
 [10] rstantools_2.0.0     inline_0.3.15        digest_0.6.25       
 [13] foreach_1.5.0        htmltools_0.4.0      rsconnect_0.8.16    
 [16] fansi_0.4.1          magrittr_1.5         globals_0.12.5      
 [19] modelr_0.1.8         gower_0.2.1          matrixStats_0.56.0  
 [22] RcppParallel_5.0.1   xts_0.12-0           prettyunits_1.1.1   
 [25] colorspace_1.4-1     skimr_2.1.1          blob_1.2.1          
 [28] rvest_0.3.5          haven_2.3.1          xfun_0.14           
 [31] callr_3.4.3          crayon_1.3.4         jsonlite_1.6.1      
 [34] lme4_1.1-23          survival_3.1-12      zoo_1.8-8           
 [37] iterators_1.0.12     glue_1.4.1           gtable_0.3.0        
 [40] ipred_0.9-9          pkgbuild_1.0.8       rstan_2.19.3        
 [43] DBI_1.1.0            miniUI_0.1.1.1       Rcpp_1.0.4.6        
 [46] xtable_1.8-4         GPfit_1.0-8          StanHeaders_2.21.0-5
 [49] lava_1.6.7           prodlim_2019.11.13   DT_0.13             
 [52] htmlwidgets_1.5.1    httr_1.4.1           threejs_0.3.3       
 [55] ellipsis_0.3.1       loo_2.2.0            pkgconfig_2.0.3     
 [58] nnet_7.3-14          dbplyr_1.4.4         utf8_1.1.4          
 [61] tidyselect_1.1.0     rlang_0.4.6          DiceDesign_1.8-1    
 [64] reshape2_1.4.4       later_1.1.0.1        munsell_0.5.0       
 [67] cellranger_1.1.0     tools_4.0.1          cli_2.0.2           
 [70] generics_0.0.2       ggridges_0.5.2       evaluate_0.14       
 [73] fastmap_1.0.1        yaml_2.2.1           ModelMetrics_1.2.2.2
 [76] processx_3.4.2       knitr_1.28           fs_1.4.1            
 [79] future_1.17.0        nlme_3.1-148         mime_0.9            
 [82] rstanarm_2.19.3      xml2_1.3.2           tokenizers_0.2.1    
 [85] compiler_4.0.1       bayesplot_1.7.2      shinythemes_1.1.2   
 [88] rstudioapi_0.11      reprex_0.3.0         tidyposterior_0.0.3 
 [91] lhs_1.0.2            statmod_1.4.34       stringi_1.4.6       
 [94] highr_0.8            ps_1.3.3             blogdown_0.19       
 [97] Matrix_1.2-18        nloptr_1.2.2.1       markdown_1.1        
[100] shinyjs_1.1          vctrs_0.3.1          pillar_1.4.4        
[103] lifecycle_0.2.0      furrr_0.1.0          data.table_1.12.8   
[106] httpuv_1.5.4         R6_2.4.1             bookdown_0.19       
[109] promises_1.1.1       gridExtra_2.3        janeaustenr_0.1.5   
[112] codetools_0.2-16     boot_1.3-25          colourpicker_1.0    
[115] gtools_3.8.2         assertthat_0.2.1     withr_2.2.0         
[118] shinystan_2.5.0      parallel_4.0.1       hms_0.5.3           
[121] grid_4.0.1           timeDate_3043.102    minqa_1.2.4         
[124] class_7.3-17         rmarkdown_2.2        pROC_1.16.2         
[127] tidypredict_0.4.5    shiny_1.4.0.2        lubridate_1.7.9     
[130] base64enc_0.1-3      dygraphs_1.1.1.6