Get started: correct machine learning made easy

Preparing the data

Most prediction functions need the predictors as a numeric matrix (a table with only numbers) and the outcome as a vector (a single row or column of elements). Let’s start with the 40 numeric variables:

x = as.matrix(patients[, 1:40])  # columns X1, X2, ..., X40
y = patients$relapse
x_new = as.matrix(new_patients[, 1:40])
y_new = new_patients$relapse

Some steps of this tutorial involve randomness. To obtain exactly the same results as shown here, we fix the “seed” of the random number generator:

set.seed(1)

A standard logistic regression

Many researchers would start with a standard logistic regression including all the variables:

m_glm = glm(y ~ ., data = data.frame(x), family = "binomial")
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

R warns us that the fitted probabilities are 0 or 1. Keep this in mind, it is a first sign of trouble. Let’s see how well the model predicts relapse in our patients. The function predict returns, for each patient, the probability of relapse, and round converts it into a prediction: 1 (relapse) if the probability is higher than 0.5, or 0 (no relapse) otherwise:

y.pred = round(predict(m_glm, type = "response"))
simple.accuracy(y, y.pred, "bac")
#>   bac bac.lwr bac.upr se    se.lwr se.upr sp    sp.lwr sp.upr
#> 1   1  0.9635       1  1 0.9448164      1  1 0.9563268      1

simple.accuracy compares the actual outcomes (y) with the predictions (y.pred). Here, it returns three numbers, each with the lower (.lwr) and upper (.upr) limits of its 95% confidence interval:

• se (sensitivity): the proportion of patients who relapsed that the model correctly predicted would relapse. 1 (i.e., 100%) means that the model correctly predicted all the patients who relapsed.
• sp (specificity): the proportion of patients who did not relapse that the model correctly predicted would not relapse. 1 (i.e., 100%) means that the model correctly predicted all the patients who did not relapse.
• bac (balanced accuracy): the mean of the two.

A perfect model! Too good to be true? Let’s use our crystal ball and apply the model to the new patients:

y_new.pred = round(predict(m_glm, data.frame(x_new), type = "response"))
simple.accuracy(y_new, y_new.pred, "bac")
#>         bac bac.lwr bac.upr        se   se.lwr    se.upr        sp    sp.lwr
#> 1 0.6777762  0.6465  0.7075 0.6770601 0.632779 0.7190669 0.6784922 0.6343522
#>      sp.upr
#> 1 0.7203496

In new patients, the balanced accuracy drops quite a bit. What happened?

The model has memorized the 100 patients instead of learning general rules. With 40 variables and only 100 patients, the model can adjust to the random particularities of most or all patients, which will not be repeated in new patients. This is called overfitting. It is like a student who memorizes the answers of a mock exam: they get a perfect score in that exam, but not in the real one.

A lasso model

A lasso model is a type of regression designed to reduce overfitting: it makes the coefficients smaller (“shrinks” them) and can remove the variables that do not seem useful. glmnet_fit fits a lasso model in one line (Solanes et al., 2022):

m_lasso = glmnet_fit(x, y, family = "binomial")
#> Warning in glmnet_fit(x, y, family = "binomial"): In easy.glmnet 1.2, the
#> default 'type.measure' for "binomial" changed from "class" to glmnet "default"
#> ("deviance"). Use type.measure = "class" to reproduce the results of previous
#> versions. This warning is shown once per session.

(If you see a warning about type.measure, it only informs users of previous versions of the package that a default setting has changed; you can ignore it.)

To apply it to some patients, use glmnet_predict (or just predict), which returns the probabilities of relapse:

y.pred = round(predict(m_lasso, x))
simple.accuracy(y, y.pred, "bac")
#>        bac bac.lwr bac.upr        se    se.lwr    se.upr        sp    sp.lwr
#> 1 0.950487 0.88875   0.978 0.9545455 0.8620847 0.9904354 0.9464286 0.8638995
#>      sp.upr
#> 1 0.9846831

Again, in the same patients used to build the model, the balanced accuracy is very high . Let’s look at the new patients:

y_new.pred = round(predict(m_lasso, x_new))
simple.accuracy(y_new, y_new.pred, "bac")
#>         bac bac.lwr bac.upr        se    se.lwr    se.upr        sp    sp.lwr
#> 1 0.7399049  0.7105 0.76725 0.6971047 0.6534387 0.7382388 0.7827051 0.7429126
#>      sp.upr
#> 1 0.8188828

The lasso model may work some better than the standard logistic regression in new patients, but its accuracy is still much lower than it seemed. Thus, whatever the method (Radua & Koutsouleris, 2024):

Lesson 1: Never evaluate a model with the same patients that you used to build it.

But in a real study, we have no crystal ball. How can we know the accuracy of a model in new patients if we only have the patients used to build it? The answer is cross-validation.

The idea

Cross-validation splits the patients into several groups called folds (usually 10). Then, for each fold:

1. It builds a model using the patients of all the other folds (the training set).
2. It applies this model to the patients of this fold (the test set).

         Fold 1   Fold 2   Fold 3   ...   Fold 10
Round 1:  test    train    train    ...    train
Round 2:  train   test     train    ...    train
Round 3:  train   train    test     ...    train
...
Round 10: train   train    train    ...    test

At the end, each patient has a prediction from a model that never saw them, just like a new patient. Comparing these predictions with the actual outcomes tells us how well the model works in “new” patients.

Doing it with cv

The function cv does all the work (dividing the patients into folds, building the models, applying them, and collecting the predictions), but it needs you to tell it two things, as two small functions:

• fit_fun: how to build a model from the training set. It receives the predictors and outcomes of the training set and must return the model.
• predict_fun: how to apply a model to the test set. It receives the model and the predictors of the test set and must return the predictions.

fit_fun = function (x_training, y_training) {
  glmnet_fit(x_training, y_training, family = "binomial")
}
predict_fun = function (m, x_test) {
  round(predict(m, x_test))
}

Note that these functions do exactly what we did before with glmnet_fit and glmnet_predict (or just predict), but with the training and test sets of each fold. Now we can run the cross-validation:

res = cv(x, y, family = "binomial", fit_fun, predict_fun, verbose = FALSE)

res$predictions contains, for each patient, the fold, the actual outcome and the prediction:

head(res$predictions)
#>   fold y y.pred
#> 1    4 0      0
#> 2    7 0      0
#> 3    5 1      0
#> 4    9 0      0
#> 5    7 0      0
#> 6    9 1      1

Let’s see the accuracy (we just have to print the cv object):

res
#> 10-fold CV of 100 binomial outcomes:
#>     bac bac.lwr bac.upr    se se.lwr se.upr    sp sp.lwr sp.upr
#> 1 0.718   0.624   0.796 0.614  0.466  0.747 0.821  0.706  0.904

The cross-validation estimates a balanced accuracy very close to the one that the lasso model actually achieved in the new patients. Unlike the evaluation with the same patients, the cross-validation tells the truth, even without a crystal ball.

A note on randomness

cv assigns the patients to the folds randomly (though ensuring that each fold has a similar proportion of patients who relapsed). With another seed, the folds would be different, and the estimate would change a bit, especially in small samples like this one. This is why we used set.seed: to make the results reproducible.

Cross-validation can also be optimistic if used wrongly

Imagine that, before the cross-validation, we select the variables with the lasso model that we fitted in Step 1 with all the patients (m_lasso; recall that the lasso sets some coefficients exactly to zero, so it uses only some of the variables), and then we cross-validate a model that uses only these variables. It seems reasonable: we only keep the variables that seem useful. Let’s try:

selected = names(m_lasso$beta)  # variables selected by the lasso of Step 1
res = cv(x[, selected], y, family = "binomial", fit_fun, predict_fun, verbose = FALSE)
res
#> 10-fold CV of 100 binomial outcomes:
#>     bac bac.lwr bac.upr    se se.lwr se.upr    sp sp.lwr sp.upr
#> 1 0.797   0.707   0.864 0.773  0.634  0.877 0.821  0.706  0.904

The cross-validation estimates a higher balanced accuracy than in step 2! But let’s use our crystal ball: we build the model with these variables and all the patients, and apply it to the new patients:

m_leak = glmnet_fit(x[, selected], y, family = "binomial")
y_new.pred = round(predict(m_leak, x_new[, selected]))
simple.accuracy(y_new, c(y_new.pred), "bac")
#>         bac bac.lwr bac.upr        se    se.lwr    se.upr        sp    sp.lwr
#> 1 0.7266209  0.6965 0.75475 0.7060134 0.6626511 0.7467295 0.7472284 0.7056266
#>      sp.upr
#> 1 0.7856688

Quite less. The cross-validation was too optimistic. Why? Because the variables were selected using all the patients, including those that later played the role of test set! In other words, the selection “saw” the test set, so these patients were no longer like new patients. This is called data leakage: information from the test set leaks into the models.

The solution

There is a simple solution (Radua & Koutsouleris, 2024):

Lesson 2: Every step that learns something from the data is part of the model. Thus, it must be inside fit_fun, and must use only the training set.