Combining regression models based on missing data patterns

Question

I have a dataset that contains a few patterns of missingness. For this dataset, I have a training set that is complete and contains all input features. My test set has complete observations for the dependent variable, but there are a few patterns of missing data for the input features. Below, I provide an example:

Example dataset

# Generate a fake test with missing data patterns
subset1 <- data.frame(yield = rnorm(500, 1000, 100)) %>%
  mutate(var1 = yield * rnorm(500, 100, 5))

subset2 <- data.frame(yield = rnorm(500, 1000, 100)) %>%
  mutate(var2 = yield * rnorm(500, 100, 10),
         var3 = -yield * rnorm(500, 100, 4))

subset3 <- data.frame(yield = rnorm(500,1000,100)) %>%
  mutate(var1 = yield * rnorm(500, 100, 5), 
         var3 = -yield * rnorm(500, 100, 4))

test <- plyr::rbind.fill(subset1,subset2,subset3)

# train set has complete observations
train <- data.frame(yield = rnorm(500,1000,100)) %>%
  mutate(var1 = yield * rnorm(500, 100, 5),
         var2 = yield * rnorm(500, 100, 10),
         var3 = -yield * rnorm(500, 100, 4))

md.pattern(plyr::rbind.fill(test,train))

Missing Data Pattern

Note the first row is the train set and the other rows will be the test set

I want to perform regression on this dataset, but I am wary of using imputation, unless I can have a strong justification for it. I could, of course, just remove the rows that contain missing data, however this will remove a very large portion of the dataset. I have encountered the concept of reduced modelling, which identifies the different patterns of missing data in the test set, training separate models based on complete subsets. From what I understood of the concept, this is what I tried:

# Look at the missing data patterns on test set
pattern <- mice::md.pattern(test,plot=F)

# initiate a list to put subsets into
model_list <- list()

for (i in 1:(nrow(pattern) - 1)) {
  # get the missing data pattern
  cols_with_no_missing <- names(which(pattern[i, -ncol(pattern)] == 1))

  # subset the dataset based on this pattern
  subset_df <- test %>%
    select(all_of(cols_with_no_missing)) %>% 
    na.omit()

  # train a model using only these columns
  dat <- train %>%
    select(all_of(cols_with_no_missing))
  
  mod <- train(yield ~ ., 
               data = dat,
               method = 'lm',
               trControl = trainControl(method = 'repeatedcv',
                                        number = 10,
                                        repeats = 10))
  
  # get some model statistics
  obspred <- data.frame(obs=subset_df$yield,
                        pred = predict(mod,newdata = subset_df))
  
  n_obs <- nrow(obspred)
  rsq <- summary(mod)$r.squared
  
  # Create a plot
  plt <- obspred %>%
    ggplot(aes(x = obs, y = pred)) + 
    geom_point() + 
    geom_abline(slope = 1) +
    annotate("text", x = -Inf, y = Inf, 
             label = paste(paste('n',n_obs,sep='='),
                           paste('R²', round(rsq,2),sep='='),
                           paste('yield',paste(cols_with_no_missing[-1], collapse = ' + '),sep=' ~ '),
                           sep='\n'), 
             hjust = -0.1, vjust = 1, size = 4)

  model_list[[i]] <- plt
}

# compare the models test performance
library(patchwork)
wrap_plots(model_list) + plot_layout(axes='collect')

Test performance

My Questions

1. I was wanting to know what are the consequences of using such an approach or whether there are more appropriate ways of doing something similar?
2. Can you create some sort of ensemble model that selects the appropriate sub-model based on the available features? If so, how would you incorporate the changing accuracy (which varies depending on the available features) into that?

My actual dataset is much bigger than this example (~90,000 observations) with 25 features, so in reality I will use random forest regression, but (I think??) the principle should be similar.

Answer 1

1. I would use an XGBoost regressor. Section 3.4 of the paper, "XGBoost: A Scalable Tree Boosting System" by Chen & Guestrin (2016), lays out how they do "sparsity-aware split finding." They explicitly note that the matrix of inputs $X$ may be sparse due to the "presence of missing values in the data." What they do is calculate "a default direction in each tree node" and then "when a value is missing in the sparse matrix $X$, the instance is classified into the default direction." Algorithm 3 in this paper shows how XGBoost calculates the default direction of each node.

2. Yes, you could build an ensemble. For example, you could fit a decision tree or LASSO for each predictor on its own. On the training and cross-validation data, you would toss out cases listwise where there is missingness, since there is only one $x$ variable per sub-model. You'd have a sub-model for each of your $p$ predictors. You would also have some estimate of the out-of-sample error on your holdout set. What you could do is average all of the sub-models where you have data for a given case, using the inverse of the errors to weight each model. I am not sure how this would perform; you could simulate data to see. Again, I would stick with XGBoost.

Answer 2

I was wanting to know what are the consequences of using such an approach

As far as I understand this is answered in the paper you reference:

• Reduced modeling is superior to imputation-based approaches (see Section 3) but has a higher computational cost (first paragraph of Section 4), as you need to fit a model for each "missingness pattern" in your test data. You write that you only have "a few patterns of missingness" in your test data, so maybe this is not an issue for you.

• Both reduced modeling and imputation-based methods assume that the data is missing completely at random (MCAR) (see Section 6). That is, the fact that data is missing is independent of the data (missing and not missing). Is this reasonable to assume in your application? Does $P_\mathrm{test}(\mathrm{available \ data} \mid \mathrm{missingness \ pattern} \ P) = P_\mathrm{train}(\mathrm{data \ masked \ according \ to }\ P)$ hold? Does it hold approximately?

there are more appropriate ways of doing something similar?

The approach is reasonable if the data is MCAR.

Can you create some sort of ensemble model that selects the appropriate sub-model based on the available features?

This is reduced modeling. However, if you don't make any assumptions on the possible missingness patterns, your ensemble has to contain $2^{\mathrm{number \ of \ features}}$ models. This quickly gets unrealistic, as is the case for your setting with 25 features.

An alternative is an "online" approach, where a new model based on a missingness pattern $P$ is fit whenever a prediction for a new test observation with a new missingness pattern $P$ is required.

If so, how would you incorporate the changing accuracy (which varies depending on the available features) into that?

I am unsure whether I understand this. The accuracy of your "ensemble model" would be the average accuracy over all test data points.

Answer 3

Since you training data is complete and has no missing data, to apply such an approach would require you to create the missingness yourself. In theory, you could do that (i.e. take the complete dataset and then create versions of it with all the possible missingness patterns you've seen - or in practice don't create copies, just use different lists of non-missing features). However, if you do that, you're making the assumption that missingness is unrelated to the outcome, when normally you'd choose this approach exactly because you're worried that it's not.

Why do I say that? The reason is that the reason you'd want to fit separate models for each missingness pattern (if you really have training data for all of them), is that you think there's really something so fundamentally different between different missingness patterns that you don't want a model that takes any information from other missingness patterns. Alternatively, you could see versions where the missingness pattern is a categorical predictor (and/or missingness in each variable is a predictor), which if you allow for interactions is almost as flexible (depending on your model type) and in the absence of a lot of data on some missingness patterns ends up shrinking towards an average behavior.

Note also that by randomly setting things missing in different patterns, your making very similar (or if you use the outcome in the imputation, which in prediction modeling you usually don't, you even make stronger) assumptions as you would have to make for some imputation approaches such as multiple imputation (sure, you may argue that you're not making distributional assumption). So, one can debate whether one makes any fewer assumptions by not imputing.

How you would select models, surely this is just a matter of "missingness pattern A" => "model trained for missingness pattern A", "missingness pattern B" => "model trained for missingness pattern B" etc., because otherwise we're not talking about training on complete data minus different missingness patterns, but it rather becomes variable selection (which with suitable tuning of some sparsity inducing hyperparameters you might want to do in an appropriately cross-validated fashion, but shrinkage, model averaging and other forms of regularization may also make sense). In principle, of course, a general way of selecting or averaging models is to cross-validate different ways of doing so with some suitable shrinkage parameter (e.g. form a weighted average of different models with regularization based on some penalty parameter towards a simple average).

Finally, as a side remark: there's something you need to be careful about, if you a) impute, b) don't impute but set data to missing different patterns in multiple copies of the dataset, or c) don't impute and mix together lots of copies of the dataset with data set to missing in different patterns (and then use a model that handles missingness in some way such as XGBoost, or some other model with a missingness indicator and imptutation as mean/median/most frequent). You need to make sure that any validation scheme (such as cross-validation splits) respect that an original record should never occur in both the training part and the validation (or test) pat of a split.