Evaluating and interpreting your model

Probabilities and thresholds

In the first tutorial, we converted the probabilities into predictions of 0 or 1 with round, i.e., we predicted a relapse if the probability was higher than 0.5. With these predictions, we could estimate the balanced accuracy:

y.pred = round(res$predictions$y.pred)
simple.accuracy(res$predictions$y, y.pred, "bac")
#>         bac bac.lwr bac.upr        se    se.lwr    se.upr        sp    sp.lwr
#> 1 0.7581169   0.667  0.8315 0.6590909 0.5124696 0.7855432 0.8571429 0.7484684
#>      sp.upr
#> 1 0.9300705

However, 0.5 is an arbitrary threshold. For instance, if relapses are very harmful and the preventive treatment is cheap and safe, a clinician might prefer to treat all patients with a probability of relapse higher than 0.2. With another threshold, the sensitivity and the specificity would be different, and so would the balanced accuracy.

The area under the ROC curve

The area under the ROC curve (AUC) measures how well the model discriminates, without choosing any threshold. It is the probability that, if we take at random a patient who relapsed and a patient who did not, the model estimates a higher probability of relapse for the first one. Thus, an AUC of 0.5 means that the model discriminates no better than tossing a coin, and an AUC of 1 means that it discriminates perfectly.

The AUC was the first number shown by res, followed by the lower (auc.lwr) and upper (auc.upr) limits of its 95% confidence interval. You can also obtain them with simple.accuracy:

simple.accuracy(res$predictions$y, res$predictions$y.pred, "auc")
#>         auc   auc.lwr   auc.upr
#> 1 0.8011364 0.7142033 0.8880695

Are the probabilities correct?

A model may discriminate well but estimate wrong probabilities. For instance, it could correctly estimate that one patient has a higher risk than another, but estimate probabilities of 60% and 40% when the actual probabilities are 30% and 20%. If a clinician started a preventive treatment in all patients with a probability higher than 50%, they would unnecessarily treat the first patient. The agreement between the estimated probabilities and what actually happens is called calibration.

Calibration-in-the-large and calibration slope

simple.calibration summarizes the calibration with two numbers, each with its 95% confidence interval:

simple.calibration(res$predictions$y, res$predictions$y.pred, "binomial")
#>   calibr.large calibr.large.lwr calibr.large.upr calibr.slope calibr.slope.lwr
#> 1  -0.06905186       -0.5616115        0.4129425    0.8131361        0.4807838
#>   calibr.slope.upr
#> 1         1.215632

• The calibration-in-the-large (calibr.large in the output) tells whether the estimated probabilities are, on average, too high (negative values) or too low (positive values). Its ideal value is 0.
• The calibration slope (calibr.slope in the output) tells whether the estimated probabilities are too extreme (values lower than 1, e.g., probabilities too close to 0 or 1) or too modest (values higher than 1, e.g., probabilities too close to the mean). Its ideal value is 1.

There are no universally accepted cutoffs, but a useful guide is to check whether the confidence intervals include the ideal values. Also, some authors consider a calibration slope of at least 0.9 acceptable, i.e., predictions no more than 10% too extreme (Riley et al., 2020).

These numbers were also shown by res. Here, the calibration-in-the-large is close to zero and the slope is close to one, i.e., the probabilities are, on average, about right, and only slightly too extreme.

Two notes of caution. First, with small samples like this one, these estimates are quite imprecise, as their wide confidence intervals show. Second, the calibration-in-the-large depends on the proportion of patients who relapse in the sample. If this proportion differs from the proportion in the population where the model will be used, the model will be miscalibrated there, however well calibrated it seems in the cross-validation. We will come back to this below.

How often is each variable selected?

The lasso selects variables, but which ones it selects may change a bit with the specific patients used to build the model. The cross-validation built 10 models, each with a slightly different group of patients, and they are stored in res$models. glmnet_get.items.count tells in which proportion of these models each variable was selected (Radua et al., submitted). We may look at those selected in at least 90% of the folds:

count = glmnet_get.items.count(res$models)
count[which(count >= 0.9)]
#> X11 X14 X17 X21 X22 X28 X34 X38 X39  X5  X7 X33 
#> 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.9

Variables selected in all or almost all models are stable: the lasso finds them whatever the patients used. Variables selected in few models are unstable, and we should be careful interpreting them.

How much does each variable weigh?

Being selected does not mean being important: a variable may be selected in all models but have a tiny coefficient. glmnet_get.items.relevance estimates the relevance of each variable, i.e., its share of the total weight of the model, from the size of its standardized coefficient (standardized so that variables measured in different units can be compared) (Solanes et al., 2022; Salazar de Pablo et al., 2025). It averages the relevance across the 10 models, and the relevances of all variables add up to 1. We may look at those with at least 5% relevance:

relevance = glmnet_get.items.relevance(res$models)
round(relevance[which(relevance >= 0.05)], 3)
#>   X21   X39   X22   X14    X5    X7   X11 
#> 0.139 0.131 0.128 0.101 0.077 0.063 0.059

The two functions answer different questions, so it is useful to look at both: a relevant variable selected in all models is much more convincing than a variable that is relevant in a few models only.

Two options

The cross-validation estimates how well the models work, but with the cv(), we built 10 models (one for each fold), so which model should you use with future patients? There are two reasonable options.

The first option is to build the model again, using the same fit_fun with all the patients:

m_all = glmnet_fit(x, y, family = "binomial")

The second option is to choose, among the 10 models built during the cross-validation, the most representative one. glmnet_get.main.model selects the model whose variables are most similar to the variables of the other models (according to the Dice coefficient), and if there are several equally representative models, it averages them (Sobregrau et al., 2024; Salazar de Pablo et al., 2025):

m_main = glmnet_get.main.model(res$models, verbose = FALSE)

The first option uses all the patients, so it may be slightly more accurate. The second option uses one of the models whose accuracy was actually estimated in the cross-validation, and it is useful when, as we will see in another tutorial, the cross-validation involves several models per fold (e.g., one for each imputation of the missing values). Let’s use our crystal ball to see how both work in the new patients:

y.pred_all = predict(m_all, x_new)
y.pred_main = predict(m_main, x_new)
simple.accuracy(y_new, y.pred_all, "auc")
#>         auc   auc.lwr   auc.upr
#> 1 0.8273424 0.8011382 0.8535465
simple.accuracy(y_new, y.pred_main, "auc")
#>         auc  auc.lwr   auc.upr
#> 1 0.8160386 0.788921 0.8431562
simple.calibration(y_new, y.pred_all, "binomial")
#>   calibr.large calibr.large.lwr calibr.large.upr calibr.slope calibr.slope.lwr
#> 1    0.2572726       0.09140591         0.423091    0.8667944        0.7511321
#>   calibr.slope.upr
#> 1        0.9902875
simple.calibration(y_new, y.pred_main, "binomial")
#>   calibr.large calibr.large.lwr calibr.large.upr calibr.slope calibr.slope.lwr
#> 1     0.253726       0.09716647        0.4102683     1.003944        0.8670087
#>   calibr.slope.upr
#> 1          1.14981

Both models discriminate similarly in the new patients, and their AUCs are in line with the one estimated by the cross-validation.

However, the calibration-in-the-large is now a bit positive: the models slightly underestimate the probability of relapse. This might be because fewer of our 100 patients relapsed compared to the new patients. This is one reason why models should be validated in independent samples, and why sharing your model is so useful (Radua & Koutsouleris, 2024): others can test it, and recalibrate it if needed, in their own patients.

Estimating the probability of relapse of a new patient

Finally, to estimate the probability of relapse of a single new patient, use predict with the data of this patient. For instance, for the first new patient:

new_patient = x_new[1, , drop = FALSE]
predict(m_main, new_patient)
#> [1] 0.9158238

The model estimates that this patient has a high probability of relapse. Note the drop = FALSE: x_new[1, ] would return a vector, whereas predict needs a matrix, even if it only has one row.