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

...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))

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:


# 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))

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


1 Answer 1



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.


fit0 <- rq(mpg ~ wt, data = mtcars, tau = 0.5)
# 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

# 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

#           [,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.).


mtcars[1, "mpg"] <- NA

imp <- mice(mtcars, m = 10)

implist <- mids2mitml.list(imp)

fit1 <- with(implist, rq(mpg ~ wt, tau = 0.5))
# 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. 

#                 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))

#           2.5%     97.5%
# [1,] 31.559696 41.518307
# [2,] -7.103394 -3.875969
  • 1
    $\begingroup$ 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? $\endgroup$
    – E James
    Aug 10, 2019 at 12:40
  • 1
    $\begingroup$ 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. $\endgroup$
    – SimonG
    Aug 13, 2019 at 9:18

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.