I am using the mice package in R to impute missing data in small study. The study investigates the effect of a behavioral intervention on the frequency of a particular behavior, i.e., count data that can range from anywhere between 0-5 times for most participants and several hundreds for some participants. I have used a package for imputation of count data (https://github.com/kkleinke/countimp) and have successfully imputed the missing data with a reasonable distribution.

Now I want to report the estimated frequency at pre-treatment and post-treatment as well as the pre- to post-change and a significance test of the change, all based on the pooled data. The convention in the field is to report medians and IQR of pre- and post-frequencies as well as the median change in frequency.

Can this be properly done with multiply imputed data? My current idea is to calculate the medians and IQR with standard errors using the quantile regression library for R (quantreg), by employing intercept only regression models for quantiles .25, .50, .75 for each imputation (pre, post, and change score), and then pooling these estimates with standard errors using mice. However, this essentially means that I am taking the mean of medians, which does not feel right.

I'm pasting code that calculates the median (50th percentile) of the pre- and post-frequencies as well as the change in frequency using quantreg and then pools the estimates using mice. After that are calculations of the corresponding means.

# demo.median contains the imputed data

# vcov.rq is called by pool to get the covariance estimate of the
# quantile regression
vcov.rq <- function(fit)summary(fit, se="nid", cov=T)$cov

l <- list()
l$pre <- summary(pool(with(demo.median, rq(x_PRE ~ 1, tau = .5))))
    l$post <- summary(pool(with(demo.median, rq(x_POST ~ 1, tau = .5))))
l$pre_post <- summary(pool(with(demo.median, rq(I(x_PRE - x_POST) ~ 1, tau = .5))))

# Print pooled medians of pre, post and change score
do.call(rbind.data.frame, l)

l <- list()
l$pre <- summary(pool(with(demo.median, lm(x_PRE ~ 1))))
    l$post <- summary(pool(with(demo.median, lm(x_POST ~ 1))))
l$pre_post <- summary(pool(with(demo.median, lm(I(x_PRE - x_POST) ~ 1))))

# Print pooled means of pre, post and change score
do.call(rbind.data.frame, l)


# Medians

           est        se        t       df     Pr(>|t|)     lo 95     hi 95 nmis        fmi    lambda
pre      8.000 1.3250825 6.037360 91.75287 3.286786e-08 5.3681774 10.631823   NA 0.02110754 0.0000000
post     2.875 1.0997478 2.614236 76.09011 1.077545e-02 0.6847043  5.065296   NA 0.16031744 0.1385330
pre_post 2.400 0.8801746 2.726732 52.96378 8.654540e-03 0.6345650  4.165435   NA 0.34977967 0.3256813

# Means

               est       se        t       df     Pr(>|t|)     lo 95    hi 95 nmis        fmi      lambda
pre      18.252632 3.417037 5.341655 92.05265 6.619670e-07 11.466153 25.03911   NA 0.02104097 0.000000000
post     12.019474 2.413928 4.979217 91.32515 2.999330e-06  7.224731 16.81422   NA 0.02889543 0.007858808
pre_post  6.233158 2.316257 2.691047 91.26050 8.471778e-03  1.632374 10.83394   NA 0.02957225 0.008535556

My question is: Is this proper? If not, can I achieve my goal of estimating absolute and change score medians from multiply imputed data in any other way?


This is a great question, one that I share. I am using a method identical to yours in a current project, but I have found no statistical publications that support this practice.

In general, the sampling distribution of the median depends on the distribution from which the median is calculated/sampled. That being said, I did some small simulations, drawing quantiles from a few different distributions, and found that the distributions of the sampled medians looked close to normal for most of the cases. I use this as back-of-the-napkin evidence that using the "mean of the medians" and the SE from the quantile regression is probably good enough for the purpose of summarizing imputed datasets.

That being said, I hope someone smarter than me does a nice simulation study to show whether I am right or not.

Brant Inman


Excellent question, indeed!

Rubin's Rules for constructing a "pooled" estimated value for the parameter of interest from the estimated values provided by each imputation rely on the normality of the underlying estimator used to produce those estimated values.

In some cases, normality can be achieved only after an appropriate transformation (e.g., log-transformation, Fisher's transformation) is applied to that estimator.

If normality (even after a transformation) is in doubt, then one can stay away from using Rubin's rules and report instead summary statistics (e.g., median, range) of the distribution of estimated values of that parameter obtained from the $m$ imputations.

See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2727536/#!po=0.537634 for more details.

So the real question is:

Does the estimator of a quantile produced by quantile regression have a sampling distribution which is normal, at least for large enough sample sizes? (If yes, the sampling distribution of the difference between two such quantiles would also be expected to be normal.) The normality of the sampling distribution of the quantile of interest can be established either by consulting the literature on quantile regression or via simulation.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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