Using projpred with brm_multiple and mice

Hi,

I tried to use projpred after multiple imputation using the mice package. It did not produce results as expected (see here for plots of the output. I suspected that this was due to a combination of lack of predictive power and because of multiple imputation outside of the statistical model with mice may not be supported. @avehtari confirmed to me that multiple imputation might not yet be supported and to open an issue here.

In case it is in fact supported already: Below code reproduces the output. I fit random noise with and without using multiple imputation. I would expect that RMSE is ~equal between models thanks to the horseshoe prior. This is the case without imputations but when using multiple imputations RMSE increases when more variables are included.

One more thing: Note that models have divergent transitions. In my real dataset I circumvented that by using different priors that are less regularizing (normal(0, 0.2) for all betas). I am aware that the horseshoe prior is important to prevent overfitting but I did not find a way yet to fit the models without divergent transitions when using it in my dataset (any hints are welcome if there are better ways to circumvent that).

  library(tidyverse)
  library(brms)
  library(projpred)
  library(mice)


  # simulate data where all 25 beta are 0
  set.seed(1)
  n <- 75
  n_feat <- 25
  X <- map_dfc(1:n_feat, function(x) {
    df <- tibble(rnorm(n, 0, 2))
    colnames(df) <- glue::glue("var_{x}")
    df
  })
  y <- rnorm(n, 0, 1)
  df <- bind_cols(tibble(y = y), X)

  # projpred 
  n <- nrow(df) # 75
  D <- ncol(df[, -1]) # 25
  p0 <- 3 # prior guess for the number of relevant variables
  tau0 <- p0/(D-p0) * 1/sqrt(n) 
  m1 <- brm(
    y ~ ., 
    family = gaussian(), 
    data = df,
    prior = prior(horseshoe(scale_global = tau0, scale_slab = 1), class = b),
    seed = 1, 
    chains = 2, 
    iter = 2000,
    # control = list(adapt_delta = 0.8)
  )
  refmodel <- get_refmodel(m1)
  vs <- cv_varsel(refmodel, method = "forward", cv_method = "LOO")
  plot(vs, stats = c("elpd", "rmse"))



  # create ~ 30% random missingness per var
  for (row in seq_along(df$var_1)) {
    for (col in seq_along(colnames(df))) {
      # ~ 30% missingness
      if (runif(0, 1, n = 1) > 0.7) {
        df[row, col] <- NA
      }
    }
  }

  # 10 imputations (predictive mean matching)
  dfimp <- mice(df, m = 10)

  # use multple imputed datasets with projpred 
  m2 <- brm_multiple(
    y ~ ., 
    family = gaussian(), 
    data = dfimp,
    prior = prior(horseshoe(scale_global = tau0, scale_slab = 1), class = b),
    seed = 1, 
    chains = 2, 
    iter = 2000,
    # control = list(adapt_delta = 0.8)
  )

  refmodel2 <- get_refmodel(m2)
  vs2 <- cv_varsel(refmodel2, method = "forward", cv_method = "LOO")
  plot(vs2, stats = c("elpd", "rmse"))

}

Run the supplied R reproduction comparing brm() with brm_multiple(), then inspect how get_refmodel() and cv_varsel() handle the resulting model. Compare the RMSE and divergent-transition behavior with and without mice imputations; done means establishing whether this workflow is supported and correcting or documenting the observed result.