Continuous and survival outcomes, regression to the mean, and multisite data

The data

Our patients completed a symptom scale twice: at baseline (score_pre) and at follow-up (score_post). We want to predict how much each patient’s symptoms will change:

change = patients$score_post - patients$score_pre
change_new = new_patients$score_post - new_patients$score_pre
summary(change)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  -29.00  -11.25   -0.50   -1.10    6.25   27.00

Negative values mean improvement, and positive values mean worsening.

Predicting the change

For a continuous outcome, use family = "gaussian". Everything else is as in the previous tutorials:

set.seed(1)
fit_fun = function (x_training, y_training) {
  glmnet_fit(x_training, y_training, family = "gaussian")
}
predict_fun = function (m, x_test) {
  predict(m, x_test)
}
res_x = cv(x, change, family = "gaussian", fit_fun, predict_fun, verbose = FALSE)
res_x
#> 10-fold CV of 100 gaussian outcomes:
#>     cor cor.lwr cor.upr mse.pred_div_mse.mean mse.mean mse.mean.lwr
#> 1 0.285  0.0942   0.456                 0.943      158          122
#>   mse.mean.upr mse.pred mse.pred.lwr mse.pred.upr calibr.large calibr.large.lwr
#> 1          195      149          108          191       -0.188            -2.62
#>   calibr.large.upr calibr.slope calibr.slope.lwr calibr.slope.upr
#> 1             2.25        0.648            0.212             1.08

For continuous outcomes, res shows:

• cor: the correlation between the predicted and the actual values, from -1 to 1. The higher, the better; its square is roughly the proportion of the variance of the outcome that the model explains.
• mse.pred_div_mse.mean: the mean squared error of the predictions divided by the mean squared error of simply predicting the mean for everyone. Values lower than 1 mean that the model predicts better than the mean (e.g., 0.8 means that the model reduces the error by 20%).
• The calibration-in-the-large (ideally 0) and the calibration slope (ideally 1), as explained in the second tutorial.

A correlation of 0.29 seems promising. But before celebrating, we must check something.

The prophecy: “the sicker shall improve, the healthier shall worsen”

Let’s build a model with only one predictor, the baseline score:

x_pre = cbind(score_pre = patients$score_pre)
x_pre_new = cbind(score_pre = new_patients$score_pre)
res_pre = cv(x_pre, change, family = "gaussian", fit_fun, predict_fun, verbose = FALSE)
res_pre
#> 10-fold CV of 100 gaussian outcomes:
#>     cor cor.lwr cor.upr mse.pred_div_mse.mean mse.mean mse.mean.lwr
#> 1 0.561    0.41   0.683                 0.685      158          122
#>   mse.mean.upr mse.pred mse.pred.lwr mse.pred.upr calibr.large calibr.large.lwr
#> 1          195      108         77.3          140     -0.00573            -2.08
#>   calibr.large.upr calibr.slope calibr.slope.lwr calibr.slope.upr
#> 1             2.07        0.976            0.688             1.26

The baseline score alone predicts the change even better, with a correlation of 0.56. The model simply predicts that the patients with high baseline scores will improve, and those with low baseline scores will worsen:

plot(patients$score_pre, change, xlab = "Baseline score", ylab = "Change")
abline(lm(change ~ patients$score_pre))

Is the baseline score really telling us something about the patients’ future? No. The scores at baseline and at follow-up are only weakly correlated (cor(patients$score_pre, patients$score_post) = 0.25 in our data), so a patient with a very high baseline score will probably have a less extreme score at follow-up, simply because extreme scores are rare. For instance, a patient who scored very high at baseline may have been having a bad week, and at follow-up, that bad week is likely over. The opposite happens with the patients who scored very low at baseline. This statistical phenomenon is called regression to the mean, and it makes the change look predictable even when nothing that we measured can tell us how a patient will evolve (Radua et al., 2026).

Why did the predictors seem to predict the change?

In our data, the predictors are strongly associated with the baseline score. Thus, the model with the predictors could estimate the baseline score of each patient and, from it, “predict” the change exactly as the baseline-only model did. Indeed, its predictions are strongly correlated with the baseline score:

cor(res_x$predictions$y.pred, patients$score_pre)
#> [1] -0.6980462

So, the predictors seemed to predict the change, but only because of regression to the mean.

Does the model add anything beyond regression to the mean?

Because of regression to the mean, a model predicting change is only interesting if it predicts better than the baseline score alone. Thus, we must build a full model, including the predictors and the baseline score, and compare it with the baseline-only model:

x_full = cbind(x, score_pre = patients$score_pre)
x_full_new = cbind(x_new, score_pre = new_patients$score_pre)
res_full = cv(x_full, change, family = "gaussian", fit_fun, predict_fun, verbose = FALSE)
simple.accuracy.test(change, res_full$predictions$y.pred, res_pre$predictions$y.pred,
                     "cor", alternative = "greater")
#>        cor1      cor2   p.value                           test.method
#> 1 0.4891928 0.5613898 0.9899259 cocor.dep.groups.overlap/williams1959

simple.accuracy.test compares the accuracy of two sets of predictions of the same patients. Here, it tests whether the correlation of the full model (the first set of predictions) is greater than that of the baseline-only model (the second set). It is not: the predictors add nothing beyond regression to the mean.

Let’s check it with the crystal ball, building the models with our 100 patients and applying them to the new patients:

m_pre = fit_fun(x_pre, change)
m_full = fit_fun(x_full, change)
c(baseline_only = cor(predict(m_pre, x_pre_new), change_new),
  full = cor(predict(m_full, x_full_new), change_new))
#> baseline_only          full 
#>     0.6105832     0.5860801

In the new patients, the full model does not predict the change better than the baseline-only model either.

Lesson 3: When predicting change, always compare your model with a model that uses only the baseline score. Otherwise, you may be reporting regression to the mean as if it was a prediction, even if your model does not include the baseline score.

Note that our example is extreme on purpose: in real data, the scores at baseline and at follow-up are usually more correlated, and thus regression to the mean is smaller. But even with moderate correlations, it can explain a large part of the apparent accuracy of a model predicting change (Radua et al., 2026).

The data

Sometimes we want to predict not only whether an event happens, but when. The column surv contains the number of years until relapse or last follower up, and whether the patient relapsed (1) or was followed without relapse until that time (0, “censored”, i.e., we only know that the patient did not relapse before that time). R stores this kind of data as a Surv object (from the package survival). Censors are marked with + when printing this object; for example, 1+ represents a patient followed for one year without relapse, whereas 2 represents a patient who relapsed at two years.

y_surv = patients$surv
y_surv_new = new_patients$surv
head(y_surv)
#> [1] 1+ 2  1  2  3  1

Predicting survival

For survival outcomes, use family = "cox", which fits a lasso Cox model (Solanes et al., 2022; Palau et al., 2022):

set.seed(1)
fit_fun = function (x_training, y_training) {
  glmnet_fit(x_training, y_training, family = "cox")
}
predict_fun = function (m, x_test) {
  predict(m, x_test)
}
res_surv = cv(x, y_surv, family = "cox", fit_fun, predict_fun, verbose = FALSE)
#> Warning in glmnet_fit(x_training, y_training, family = "cox"): In easy.glmnet
#> 1.2, an artificial intercept is used for "cox" to keep the mean of linear
#> predictors constant across folds. Use cox.a0 = FALSE to reproduce the results
#> of previous versions. This warning is shown once per session.
res_surv
#> 10-fold CV of 100 cox outcomes:
#>     hr hr.lwr hr.upr hr.per.sd hr.per.sd.lwr hr.per.sd.upr calibr.slope
#> 1 2.29   1.34   3.91      1.51          1.16          1.96        0.829
#>   calibr.slope.lwr calibr.slope.upr
#> 1            0.295             1.36

(If you see a warning about cox.a0, it only informs users of previous versions of the package that, for Cox models, glmnet_fit now adds an artificial intercept so that the risk scores of the different folds are on the same scale; you can ignore it.)

For Cox models, predict does not return a probability but a risk score (the “linear predictor”): the higher, the sooner the patient is expected to relapse. res_surv shows:

• hr: the hazard ratio associated with one unit of the risk score.
• hr.per.sd: the hazard ratio associated with one standard deviation of the risk score. Unlike hr, it does not depend on the scale of the risk score, so it is the one to use to compare models. For instance, a value of 2 means that the hazard of relapse doubles with each standard deviation of the risk score.
• The calibration slope (ideally 1).

Discrimination at different times

We may also want to know how well the model discriminates the patients who relapse before a given time from those who do not. surv2binary converts the survival outcome into a binary outcome at each time:

y_binary = surv2binary(y_surv)
y_binary[[1]]$time
#> [1] 1
head(y_binary[[1]]$status)
#> [1] 0 0 1 0 0 1

Then we can estimate the AUC at each time. Note that, at a given time, the status of the patients censored before that time is unknown (NA), so we must exclude them:

for (i in 1:4) {
  status = y_binary[[i]]$status
  known = which(!is.na(status))
  auc = simple.accuracy(status[known], res_surv$predictions$y.pred[known], "auc")
  cat("Year", y_binary[[i]]$time, ": AUC =", round(auc$auc, 2),
      "in", length(known), "patients\n")
}
#> Year 1 : AUC = 0.72 in 100 patients
#> Year 2 : AUC = 0.65 in 84 patients
#> Year 3 : AUC = 0.63 in 65 patients
#> Year 4 : AUC = 0.5 in 57 patients

The AUC decreases with time because, in our data, the relapses after the first year are unrelated to the predictors (we know it because we simulated the data this way).

The data

In this part, we want to predict whether each patient will need surgery (needs_surgery: 1 = yes, 0 = no). Our 100 patients come from four hospitals, which have very different proportions of patients requiring surgery:

site = patients$site
table(site, patients$needs_surgery)
#>             
#> site          0  1
#>   Hospital A  3 22
#>   Hospital B  7 13
#>   Hospital C 15  9
#>   Hospital D 24  7

In addition, our patients have three laboratory measures (L1, L2, and L3). These measures are not related to the need of surgery, but each hospital’s laboratory is calibrated slightly differently, so their values differ between hospitals:

round(tapply(patients$L1, site, mean), 1)
#> Hospital A Hospital B Hospital C Hospital D 
#>       33.0       27.8       22.4       18.7

The effects of the site

Let’s see if the laboratory measures predict the need of surgery:

x_lab = as.matrix(patients[, c("L1", "L2", "L3")])
x_lab_new = as.matrix(new_patients[, c("L1", "L2", "L3")])
y = patients$needs_surgery
y_new = new_patients$needs_surgery

set.seed(1)
fit_fun = function (x_training, y_training, ...) {
  glmnet_fit(x_training, y_training, family = "binomial")
}
predict_fun = function (m, x_test, ...) {
  predict(m, x_test)
}
res_lab = cv(x_lab, y, family = "binomial", fit_fun, predict_fun, verbose = FALSE)
res_lab
#> 10-fold CV of 100 binomial outcomes:
#>     auc auc.lwr auc.upr calibr.large calibr.large.lwr calibr.large.upr
#> 1 0.796   0.708   0.884     -0.00359           -0.455            0.448
#>   calibr.slope calibr.slope.lwr calibr.slope.upr
#> 1        0.961            0.551             1.44

An AUC of 0.8! Did we discover a new biomarker? (We hid the warning about type.measure, which you can ignore as explained in the first tutorial.) Let’s apply the model to the new patients, who come from a fifth hospital:

m_lab = fit_fun(x_lab, y)
y_new.pred = predict(m_lab, x_lab_new)
simple.accuracy(y_new, y_new.pred, "auc")
#>         auc  auc.lwr   auc.upr
#> 1 0.5435303 0.505887 0.5811737

The model does not work in the new hospital. What happened is that the model learned to recognize the hospital from the laboratory measures, and the hospitals differ in their proportion of patients requiring surgery. Thus, it predicted a high risk for the patients who seemed to come from the hospital with most surgeries, and a low risk for those who seemed to come from the hospital with fewest. This does not work in another hospital (Radua & Koutsouleris, 2024, rule 5).

Considering the site

To detect this problem, give the sites to cv with the argument site:

set.seed(1)
res_lab = cv(x_lab, y, family = "binomial", fit_fun, predict_fun, site = site,
             verbose = FALSE)
res_lab
#> 10-fold CV of 100 binomial multisite outcomes:
#>     auc site.method      auc.method
#> 1 0.556       covar adapted AROC.sp
#> Note: 'multisite.accuracy' called with default parameters. Use 'multisite.accuracy' directly to use different parameters.

When you give the sites:

• cv assigns the patients to the folds separately for each site, so that all folds have patients from all sites.
• cv gives the site of the training and test patients to fit_fun and predict_fun, respectively. Our functions ignore them (this is why we added ... to their arguments), but you could use them, e.g., to remove the differences between sites in the laboratory measures. Following Lesson 2, these differences must then be estimated using only the training patients.
• res estimates the accuracy with the function multisite.accuracy of the package multisite.accuracy, which controls for the effects of the site (Solanes et al., 2021). Here, the AUC drops to around 0.5 or lower, i.e., the laboratory measures do not predict the need of surgery within each hospital.

Removing the differences between sites

In multisite studies, it is also convenient to try to remove the differences between sites in the predictors, e.g., with the ComBat method of the package combat.enigma (Radua et al., 2020). Following Lesson 2, the differences between sites must be estimated in the training set (combat_fit) and then removed from both the training and the test sets (combat_apply). This is precisely what the sites given to fit_fun and predict_fun are for:

library(combat.enigma)
fit_fun = function (x_training, y_training, site_training) {
  combat = combat_fit(x_training, site_training, verbose = FALSE)
  x_training = combat_apply(combat, x_training, site_training, verbose = FALSE)$dat.combat
  list(
    combat = combat,
    lasso = glmnet_fit(x_training, y_training, family = "binomial")
  )
}
predict_fun = function (m, x_test, site_test) {
  x_test = combat_apply(m$combat, x_test, site_test, verbose = FALSE)$dat.combat
  predict(m$lasso, x_test)
}
res = cv(x_lab, y, family = "binomial", fit_fun, predict_fun, site = site,
         verbose = FALSE)
res

Fortunately, a benchmark study found that ComBat is enough to remove the effects of the site for lasso models, but not for other techniques such as random forests, support vector machines, or Gaussian processes (Solanes et al., 2023; see also the MAREoS benchmark at https://www.imardgroup.com/mareos-benchmark/benchmark.html). Note, however, that ComBat cannot be applied to patients from a site that was not in the training set unless we have some data from that site to estimate its differences.

Of course, the best test of a model is to apply it to patients from other sites, as we did with the crystal ball.