Missing data

Removing the patients with missing values

We could use only the patients with complete data:

complete = complete.cases(x)
sum(complete)
#> [1] 51

We would lose 49 of our 100 patients, and with them, a lot of information. Moreover, the model could not predict the future patients with missing values, i.e., 51% of our new patients.

Removing the variables with missing values

We could also remove X21 and X39. But these were the two most relevant variables in the previous tutorial! Removing them would waste all the information that they have in the 70% of patients in whom they were recorded.

Imputing the missing values

The better solution is to impute the missing values, i.e., to replace them with plausible values estimated from the other variables. For instance, if X21 is associated with other variables, we can estimate the value of X21 of a patient from the values of these other variables.

Note, however, that simply replacing each missing value with its prediction (e.g., with predict) is not enough. A prediction is the expected value for a patient, so all the imputed values would lie exactly on the regression line: they would vary less than real values, and they would be more strongly associated with the other variables than real values are. Thus, a good imputation must add random noise of the right size to the predictions, so that the imputed values vary as much as real values would.

How it works

impute.glmnet.matrix_fit fits, for each variable, lasso models that estimate it from the other variables (Solanes et al., 2022; Palau et al., 2023; Salazar de Pablo et al., 2025). It does it even for the variables without missing values, because future patients may have missing values in any variable. Then, impute.glmnet.matrix uses these models to replace the missing values.

As explained above, impute.glmnet.matrix does not simply use the predictions of these models: it adds random noise with the size of the prediction error of each imputation model (for binary variables, it randomly assigns 0 or 1 according to the predicted probability). Moreover, an imputed value is only an estimate, and pretending that we know it exactly would make us overconfident. For this reason, impute.glmnet.matrix imputes the values several times, creating several complete datasets, each with slightly different imputed values. This is called multiple imputation. Then, we will fit a model with each imputed dataset, and average their predictions.

Doing it

As with data.frame2glmnet.matrix in the first tutorial, there are two functions: one that learns how to impute (impute.glmnet.matrix_fit) and one that applies the imputation (impute.glmnet.matrix). Let’s see what they do with all our patients:

m_impute = impute.glmnet.matrix_fit(x, verbose = FALSE)
x_imputed = impute.glmnet.matrix(m_impute, x, nimp = 5, verbose = FALSE)

nimp is the number of imputations. x_imputed is a list of 5 complete matrices. For example, these are the first rows of the first two imputed datasets:

x_imputed[[1]][1:6, c("X20", "X21", "X22", "X39")]
#>        X20      X21      X22      X39
#> 1 14.61144 30.09261 28.00546 22.44802
#> 2 22.72770 26.22413 11.16442 20.28340
#> 3 24.17220 21.78907 25.90762 38.00132
#> 4 29.48826 15.32596 36.30944 18.24790
#> 5 24.85255 30.19348 28.55986 31.20224
#> 6 34.29194 23.58225 30.78636 31.80661
x_imputed[[2]][1:6, c("X20", "X21", "X22", "X39")]
#>        X20      X21      X22      X39
#> 1 14.61144 25.71262 28.00546 22.44802
#> 2 22.72770 26.22413 11.16442 20.28340
#> 3 24.17220 21.78907 25.90762 38.00132
#> 4 29.48826 15.32596 36.30944 18.24790
#> 5 24.85255 30.19348 28.55986 31.20224
#> 6 34.29194 23.58225 30.78636 31.80661

The observed values are the same in all the imputed datasets, whereas the imputed values (those that were NA) differ slightly. In this tutorial, we use 5 imputations so that the code runs faster, but in a real study, we recommend using more (the default is 20).

Note that we only imputed the predictors, never the outcome: y was not included in x. In the new patients, the outcome is what we want to predict, so it cannot be used to impute anything. Similarly, patients whose outcome is missing cannot be used to build the model, so they must be excluded (with survival outcomes, patients lost to follow-up are instead censored, as you will see in the last tutorial).

Lesson 2 again

Imputing learns from the data: impute.glmnet.matrix_fit finds the models to estimate the missing values, and it learns these models from the patients we give it. Thus, following Lesson 2 of the first tutorial, the imputation must be inside fit_fun, using only the training patients. Then, predict_fun imputes the test patients with the models learned from the training patients, exactly as we would do with a future patient.

The code

Now, the model has two parts: the imputation models, and one lasso model for each imputed dataset. fit_fun returns both in a list:

nimp = 5
fit_fun = function (x_training, y_training) {
  m = list(
    impute = impute.glmnet.matrix_fit(x_training, verbose = FALSE),
    lasso = list()
  )
  x_imputed = impute.glmnet.matrix(m$impute, x_training, nimp = nimp, verbose = FALSE)
  for (i in 1:nimp) {
    m$lasso[[i]] = glmnet_fit(x_imputed[[i]], y_training, family = "binomial")
  }
  m
}

predict_fun imputes the test patients, applies each lasso model to its imputed dataset, and averages the predictions:

predict_fun = function (m, x_test) {
  x_imputed = impute.glmnet.matrix(m$impute, x_test, nimp = nimp, verbose = FALSE)
  y.pred = NULL
  for (i in 1:nimp) {
    y.pred = cbind(y.pred, predict(m$lasso[[i]], x_imputed[[i]]))
  }
  rowMeans(y.pred)
}

(cbind puts the predictions of the different imputations side by side, as columns of a matrix, and rowMeans averages each row, i.e., each patient.)

Now we can run the cross-validation as usual. This is the step that takes a few minutes:

res = cv(x, y, family = "binomial", fit_fun, predict_fun, verbose = FALSE)
res
#> 10-fold CV of 100 binomial outcomes:
#>     auc auc.lwr auc.upr calibr.large calibr.large.lwr calibr.large.upr calibr.slope
#> 1 0.828   0.746    0.91     -0.00995           -0.464            0.438         1.35
#>   calibr.slope.lwr calibr.slope.upr
#> 1            0.838             1.97

(We hid the warning about type.measure, which you can ignore as explained in the first tutorial.)

Despite the missing values, the cross-validation estimates an AUC similar to the one we obtained in the previous tutorial with complete data.

The calibration slope, however, is higher than 1, i.e., the estimated probabilities are somewhat too modest. This may be because the imputed values of X21 and X39 are estimated from the other variables only, so they are less associated with relapse than the real values would be, and the model relies a bit less on these variables than it should. Keep in mind, though, that with 100 patients, the calibration slope is quite imprecise (see the second tutorial for how to interpret it).

Several models per fold

In the previous tutorial, each fold had one lasso model. Now, each fold has 5 lasso models (one per imputation), stored in the lasso element of the list returned by fit_fun. The functions to interpret the models accept this structure: we only have to tell them the name of this element with childname:

count = glmnet_get.items.count(res$models, "lasso")
count[which(count >= 0.9)]
#>  X11  X14  X22  X38  X39   X5   X7  X21  X28 
#> 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.98 0.94
relevance = glmnet_get.items.relevance(res$models, "lasso")
round(relevance[which(relevance >= 0.05)], 3)
#>   X39   X22   X21   X14   X38    X5    X7   X11 
#> 0.129 0.128 0.120 0.115 0.079 0.078 0.078 0.057

Here, the count and the relevance are calculated across the 10 folds × 5 imputations = 50 models.

The model to use with new patients

As in the previous tutorial, there are two options. The first option is to build the model with all the patients, using fit_fun. The result includes the imputation and the lasso models, so predict_fun can apply it to new patients, including those with missing values:

m_all = fit_fun(x, y)
y.pred_all = predict_fun(m_all, x_new)
cbind(simple.accuracy(y_new, y.pred_all, "auc"),
      simple.calibration(y_new, y.pred_all, "binomial"))
#>         auc   auc.lwr   auc.upr calibr.large calibr.large.lwr calibr.large.upr
#> 1 0.8078015 0.7799928 0.8356103    0.2020289       0.05262644        0.3514366
#>   calibr.slope calibr.slope.lwr calibr.slope.upr
#> 1     1.185639         1.021262         1.360574

The AUC in the new patients is in line with the cross-validation estimate, and the calibration slope is also higher than 1. The calibration-in-the-large is positive for the reason explained in the previous tutorial: fewer of our 100 patients relapsed than of the new patients.

The second option is to choose the most representative lasso model with glmnet_get.main.model:

m_main = glmnet_get.main.model(res$models, "lasso", verbose = FALSE)
m_main
#> 'glmnet_fit' of binomial outcomes:
#>  (Intercept)           X5           X7          X11          X14          X17 
#> -6.398284081  0.064851256  0.043309458  0.056137214  0.085388405  0.020602076 
#>          X19          X21          X22          X23          X27          X28 
#> -0.015577887 -0.078470447  0.076552895  0.007648375  0.006067221 -0.020975668 
#>          X32          X33          X34          X38          X39          X40 
#> -0.044776576  0.007032242  0.017432266 -0.074625196  0.112609397 -0.017468938

This model is useful to report and interpret the model (a single formula instead of 50 lasso models). However, to apply it to new patients with missing values, we would also need imputation models, so use the first option to predict new patients.