2
$\begingroup$

Here is my understanding:

LightGBM by default handles missing values by putting all the values corresponding to a missing value of a feature on one side of a split, either left or right depending on which one maximizes the gain.

It seems that this approach should performa better than simple imputation with mean, median or mode, because in all cases missing values are all grouped together, with the default approach the correspond to NaN and with either of the imputations to a single number. However, in addition to imputed values, let's say median, there are other values that are not imputed that are very close, that might confuse the model.

Of course we can add a binary feature for every feature that has missing values that keeps track where the corresponding value in a feature is imputed or not, but then we get very similar behavior as to what we get the default approach.

We could also train a separate model, that predicts missing values of a feature based on other features in the dataset, but even here if a feature_a is strongly correlated with feature_with_na, then wouldn't the default approach get the same information from feature_a and then it still has not missing values of feature_with_na so it should have the same information. If no features are correlated with feature_with_na, then we can't reliable predict missing values.

I tried building a toy example were any of the imputation strategies perform better than the default missing values approach, but I could not do it, in my case default approach always does better.

Does the default missing values handling better (at least not worse) than any imputation strategy for LightGBM (or any Gradient Boosting algorithm)? If not what am I missing?

$\endgroup$

1 Answer 1

2
$\begingroup$

Firstly, the default is in many prediction situations pretty good. It's an improvement over solely imputing e.g. a mean or median and having a special 0-1 variable for whether the value was missing (because both approaches more or less do the same thing/have the same limitations, but the LightGBM default requires fewer splits in trees to capture the predictive value of missingness).

Secondly, where more sophisticated imputation will outperform the approach is when you understand something about the data generating mechanism. E.g. imagine you are trying to learn a linear relationship like $Y_i = \beta_0 + \beta_1 \times X_i + \beta_2 Z_i + \epsilon_i$, that any value of $X_i$ that is $>x_\max$ gets set to missing and that $X_i$ and $Z_i$ are somewhat weakly or moderately correlated. This is a situation, where you'd expect something like a multiple imputation that captures that the the missing values must be $>x_\max$ + fitting LightGBM to each imputation (and then in the end averaging the predictions across models) to outperform the default handling in LightGBM by a good bit.

EDIT/ADDITION: Illustrative example

In the example that follows what I described above, the default imputation actually does a bit better than a single value imputation to a value above the censoring threshold (that surprised me a bit), but both areoutperformed by multiple imputation using our knowledge of the problem at hand. The RMSE results with 90% CIs are below:

approach     rmse rmse.lcl rmse.ucl
       1 1.741364 1.734353 1.748346
       2 1.769537 1.762528 1.776518
       3 1.727564 1.720581 1.734519

Of course, you could repeat this simulation several times and see what you find. Here's the code for this:

library(tidyverse)
library(lightgbm)

cut_point <- qnorm(0.25)

simulate_data <- function(samples, rho=0.5){
  tibble(epsilon = rnorm(samples, sd=0.5),
         x = rnorm(samples),
         u = rnorm(samples),
         z = rho*x + sqrt(1-rho^2) * u,
         y = 5 + 2.5 * x - 2.5 * z + epsilon) %>%
    mutate(x = ifelse(x>cut_point, NA_real_, x)) %>%
    dplyr::select(y, x, z)
}

set.seed(123)
# Simulate a lot of training and test data
train_data <- simulate_data(samples=250)
test_data <- simulate_data(samples=100000)

# Relying on default imputation
dtrain1 <- lgb.Dataset(train_data %>% 
                         dplyr::select(x,z) %>% 
                         as.matrix(), 
                      label = train_data %>% pull(y))

lgb.cveval1 <- lgb.cv(params = list(objective = "regression", 
                                    metric = "l2",
                                    num_iterations=1000,
                                    learning_rate=0.01,
                                    bagging_fraction=0.7,
                                    feature_fraction=1.0),
                      nfold=5,
                      data = dtrain1)


lgb.model1 <- lgb.train(params = list(objective = "regression", 
                                      metric = "l2",
                                      num_iterations=lgb.cveval1$best_iter,
                                      learning_rate=0.01,
                                      bagging_fraction=0.7,
                                      feature_fraction=1.0),
                        data = dtrain1)

# Manually imputing to be higher than censoring value
dtrain2 <- lgb.Dataset(train_data %>% 
                         dplyr::select(x,z) %>% 
                         mutate(x=ifelse(is.na(x),cut_point+0.1,x)) %>% 
                         as.matrix(), 
                       label = train_data %>% pull(y))
lgb.cveval2 <- lgb.cv(params = list(objective = "regression", 
                                      metric = "l2",
                                      num_iterations=1000,
                                      learning_rate=0.01,
                                      bagging_fraction=0.7,
                                      feature_fraction=1.0),
                     nfold=5,
                     data = dtrain2)
lgb.model2 <- lgb.train(params = list(objective = "regression", 
                                   metric = "l2",
                                   num_iterations=lgb.cveval2$best_iter,
                                   learning_rate=0.01,
                                   bagging_fraction=0.7,
                                   feature_fraction=1.0),
                     data = dtrain2)

# Multiple imputation

library(brms)

imputation_model <- brm(x |  cens(censor) + trunc(lb=lbx) ~ 0 + z,
                        data = train_data %>%
                          dplyr::select(x,z) %>%
                          mutate(censor = ifelse(is.na(x), "right", "none"),
                                 lbx = -Inf,
                                 x = ifelse(is.na(x), cut_point, x)))

imputations <- predict(object = imputation_model, 
                       newdata = train_data %>%
                         filter(is.na(x)) %>%
                         dplyr::select(z) %>%
                         mutate(censor = "none",
                                lbx = cut_point), 
                       summary = F)

test_imputations <- predict(object = imputation_model, 
                       newdata = test_data %>%
                         filter(is.na(x)) %>%
                         dplyr::select(z) %>%
                         mutate(censor = "none",
                                lbx = cut_point), 
                       summary = F)

test_preds <- list()

for (imputation in 1:100){
  tmp_train <- train_data %>% 
    filter(is.na(x)) %>%
    dplyr::select(-x) %>%
    bind_cols(tibble(x=imputations[imputation,])) %>%
    bind_rows(train_data %>%
                filter(!is.na(x)))
  
  dtrain3 <- lgb.Dataset(tmp_train %>% 
                           dplyr::select(x,z) %>%
                           as.matrix(), 
                         label = tmp_train %>% pull(y))
  if (imputation==1){
    lgb.cveval3 <- lgb.cv(params = list(objective = "regression", 
                                        metric = "l2",
                                        num_iterations=1000,
                                        learning_rate=0.01,
                                        bagging_fraction=0.7,
                                        feature_fraction=1.0),
                          nfold=5,
                          data = dtrain3)
    
  }
  
  lgb.model3 <- lgb.train(params = list(objective = "regression", 
                                        metric = "l2",
                                        num_iterations=lgb.cveval3$best_iter,
                                        learning_rate=0.01,
                                        bagging_fraction=0.7,
                                        feature_fraction=1.0),
                          data = dtrain3)
  
  tmp_test <- test_data %>% 
    filter(is.na(x)) %>%
    dplyr::select(-x) %>%
    bind_cols(tibble(x=test_imputations[imputation,])) %>%
    bind_rows(test_data %>%
                filter(!is.na(x)))
  
  test_preds[[imputation]] <- predict(lgb.model3,
                                      tmp_test %>% 
                                        dplyr::select(x,z) %>%
                                        as.matrix())
  
}



avg_test_preds <- map_dbl(1:dim(test_data)[1], 
                          function(x) mean(map_dbl(
                            1:length(test_preds), 
                            function(y) test_preds[[y]][x])))

# Evaluate the different approaches

tibble(approach=1,
       predy = predict(lgb.model1, 
                  test_data %>% 
                    dplyr::select(x,z) %>% 
                    as.matrix()),
       y = test_data$y) %>%
  bind_rows(
    tibble(approach=2,
           predy=predict(lgb.model2, 
                         test_data %>% 
                           dplyr::select(x,z) %>%
                           mutate(x=ifelse(is.na(x), cut_point+0.1, x)) %>%
                           as.matrix()),
           y=test_data$y)) %>%
  bind_rows(
    tibble(approach=3,
           predy=avg_test_preds,
           y=tmp_test$y)) %>%
  group_by(approach) %>%
  summarize(rmse = sqrt(mean( (predy-y)^2)),
            rmse.lcl = sqrt( rmse^2 - sd( (predy-y)^2 ) / sqrt(n()) * qnorm(0.95) ),
            rmse.ucl = sqrt( rmse^2 + sd( (predy-y)^2 ) / sqrt(n()) * qnorm(0.95) )) %>%
  data.frame()
$\endgroup$
6
  • $\begingroup$ The question has a paragraph arguing that the tree can just split on $Z$ then. But I think if you make $X$ a linear combination of more than one other independent variable, then the tree would have to make lots of layers of splits just to recover the masked values of $X$, which is even less likely. $\endgroup$ Commented Feb 15, 2023 at 13:56
  • $\begingroup$ Thank you for the answer. I am a bit confused about the approach on improving the default performance. If every value in X that is greater than $x_\max$ gets set to missing, and we use Z to estimate X. In training data all X values are less or equal to $x_\max$, and therefore LightGBM will not predict anything greater than that, but we also want to make sure that imputed values are not less then $x_\max$, so we essentially fill every missing value with $x_\max$. Perhaps, I just misunderstood. $\endgroup$
    – Akavall
    Commented Feb 16, 2023 at 8:53
  • 1
    $\begingroup$ If you use multiple imputation assuming (in this hypothetical case correctly) that there's a multivariate normal distribution for X and Z jointly, and uses the information on censoring (which you can either explicitly build in or do by discarding samples that are not fitting the constraint), you should be able to impute pretty well. Sure, a tree-base imputation model can't do this. LightGBM would eventually with enough data learn that missing values behave like values at the maximum of its range and how to use the information from Z, but it would be a lot less efficient than imputation. $\endgroup$
    – Björn
    Commented Feb 16, 2023 at 10:17
  • $\begingroup$ Thanks a lot for the code example! I am not too familiar with R, but will go through it. $\endgroup$
    – Akavall
    Commented Feb 17, 2023 at 7:37
  • 2
    $\begingroup$ Sort of. The example I simulated is actually something that LightGBM should cope with on its own given enough data, it's just more data efficient to help it by giving it the extra information we have by doing an appropriate imputation for the particular situation. There's situations that would be much worse that it would not be able to cope with on its own. The scenario of removing the highest 20% of the outcome variable that you describe would certainly be tough on LightGBM (it's a case of missing not at random that all imputation methods should struggle with to some extent). $\endgroup$
    – Björn
    Commented Feb 20, 2023 at 9:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.