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

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!!

# 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

• This is enormously helpful, thank you so much for taking the time to explain the issue to me! I'm really pleased to have it working following your example. I've been playing around with various methods to see if I can get to the point I need to be able to (comparing models across different quantiles). Just to check - am I right in saying that the comparing models in the mitml formats only work for nested models? Aug 10, 2019 at 12:40
• I assume you are referring to testModels() in mitml. Yes, only nested models can be compared in this way. The added vcov method will make it possible to compare nested models using the D1 method, which is essentially a pooled Wald test over multiple parameters. Aug 13, 2019 at 9:18