Pooling across quantile regression analyses from multiple imputed datasets (quantreg, MICE)

Question

We have a dataset looking at predictors of reading comprehension ability, with a few missing data points here and there. After lots of going round in circles I think that multiple imputation is the best option for dealing with the missing data, and have been testing this out in a basic regression model. E.g., ...

mult.imp <- mice(raw_data)
mult.mod <- with(mult.imp, lm(comp ~ ageMonths + nonverbal + vocab))
summary(pool(mult.mod)

...and everything seems to be working fine and as expected.

However, we need to conduct quantile regression models on these analyses, and I can't seem to make the pool() function play ball with the quantile regression output. For example... (just at a single quantile for simplicity)

mult.rq <- with(mult.imp, rq(comp ~ ageMonths + nonverbal + vocab, tau = 0.5))
summary(pool(mult.rq)

This gives me the error:

Error in rq.fit.br(x, y, tau = tau, ci = TRUE, ...) : unused arguments (effects = "fixed", exponentiate = FALSE)

And sometimes also this:

Vectorizing 'logLik' elements may not preserve their attributes

(Note the quantile regression code on the single non-imputed dataset is fine)

Does anyone know if/how this can be worked around? Sorry, I'm pretty new to this but here is my attempt at a reproducible example:

library(mice)
library(quantreg)

# Introduce NAs to engel dataset from quantreg package
dataNA <- as.data.frame(lapply(engel, function(cc) cc[ sample(c(TRUE, NA), prob = c(0.90, 0.10), size = length(cc), replace = TRUE) ]))

# Create datasets using multiple imputations
imp <- mice(dataNA)

# Run quantile regression model
rq.mod <- with(imp, rq(foodexp ~ income, tau = 0.5))
summary(pool(rq.mod))

Any help would be extremely gratefully received, as I think we might have to lose data if I can't make it work (alternative suggestions welcome...). Many thanks in advance!!

Answer 1

Problem

For quantile regression, there is no agreed upon method to calculate standard errors (SEs), which are usually required to pool results under MI. For this reason, the vcov() method, which normally gives the variance-covariance matrix of the model parameters, is undefined for quantile regression.

library(quantreg)
data(mtcars)

fit0 <- rq(mpg ~ wt, data = mtcars, tau = 0.5)
summary(fit0)
# 
# Call: rq(formula = mpg ~ wt, tau = 0.5, data = mtcars)
# 
# tau: [1] 0.5
# 
# Coefficients:
#             coefficients lower bd upper bd
# (Intercept) 34.23224     32.25029 39.74085
# wt          -4.53947     -6.47553 -4.16390

vcov(fit0)
# Error (because undefined).

Provisional solution

In case you have a preferred method for calculating SEs, you can define your own vcov() method and extract the variance-covariance matrix from the summary().

vcov.rq <- function(object, ...) summary(object, se = "nid", covariance = TRUE)$cov

vcov(fit0)
#           [,1]       [,2]
# [1,]  4.909645 -1.5705424
# [2,] -1.570542  0.5499756

Notice that this will hard-code the method for calculating SEs, so it should be used with great care! I also wouldn't recommend doing this for bootstrap SEs, because bootstrap results are best pooled directly and not via SEs!

Once you have that, any package that uses the vcov() method to pool results under MI should work out-of-the-box (e.g., mitml, mitools, miceadds etc.).

library(mice)
library(mitml)

mtcars[1, "mpg"] <- NA

imp <- mice(mtcars, m = 10)

implist <- mids2mitml.list(imp)

fit1 <- with(implist, rq(mpg ~ wt, tau = 0.5))
testEstimates(fit1)
# 
# Call:
# 
# testEstimates(model = fit1)
# 
# Final parameter estimates and inferences obtained from 20 imputed data sets.
# 
#              Estimate Std.Error   t.value        df   P(>|t|)       RIV       FMI 
# (Intercept)    35.160     2.712    12.964   228.857     0.000     0.405     0.294 
# wt             -4.788     0.849    -5.641   425.155     0.000     0.268     0.215 
# 
# Unadjusted hypothesis test as appropriate in larger samples. 
# 

confint(testEstimates(fit1))
#                 2.5 %    97.5 %
# (Intercept) 29.816379 40.503897
# wt          -6.456713 -3.119769

Other options

One other option would be to conduct a bootstrap with each of the imputed data sets and pool the bootstrap estimates across imputations to obtain a (naïve) percentile-based confidence interval. This method is a bit more difficult to apply, because there is currently no R function that automatizes this (that I am aware of).

Note that this method requires more imputations than pooling based on SEs: Link.

bs1 <- with(implist, boot.rq(y = mpg, x = model.matrix(mpg ~ wt), R = 5000)$B)

bs1.pooled <- do.call(rbind, bs1)
bs1.ci <- apply(bs1.pooled, 2, quantile, probs=c(.025, .975))

t(bs1.ci)
#           2.5%     97.5%
# [1,] 31.559696 41.518307
# [2,] -7.103394 -3.875969