Is there any equivalent median values to Estimated Marginal Means?

Question

I am wondering if there is any sort of measure/R package/statistical tool to compute "balanced" median values for groups? I would like to do a pairwise comparison of these median values and plot them like one would with the emmeans package.

Thanks

Answer 1

Technically, I think the answer is no. However, you could do something that approaches that. For starters, you could specify fac.reduce = median in the emmeans() call. That will cause it to use the marginal medians rather than the marginal means.

However, that is not going far enough, because emmeans() summarizes a model, and most models predict the mean. So it would be desirable that the model try to predict the medians, or at least some robust measure of location, instead of the means. For example, fit the model using MASS::lms().

It is still touchy because the covariance structure of the marginal medians is not well estimated by these methods. So I wouldn't be too literal in interpreting tests and CIs. If you can fit a Bayesian robust model, then the posterior sample of the marginal medians would be more trustworthy.

Addendum

The more I think about this, the more I think this is fraught with problems, especially in terms of using the emmeans package. The package is quite heavily dependent on linearity of the estimation procedure; and the median is not a linear function of the data.

Just for sake of discussion, suppose we have some kind of unbalanced dataset with two factors $A$ and $B$ in a completely randomized design. For starters, just computing the marginal medians for $A$ and $B$ would be misleading for the same reason that the marginal means are -- which is where your question comes in.

So what can we do? Well, one thing that comes to mind -- parallel to what is done with EMMs -- is to compute the $A\times B$ (cell) medians, then compute the marginal medians of those. But note that even in a balanced experiment, that procedure will not yield the marginal medians of the data.

So another approach would be to use the cell counts as weights and compute weighted marginal medians: those would be a value such that half the weight goes to data above that value and the other half of the weight goes below. That procedure would yield the marginal medians if the data were balanced. But it is also completely outside of the algorithms used in the emmeans package. I do not think adding that capability to this package would be a good idea -- not that it wouldn't be a good thing, just that it'd be trying to fit a square peg in a round hole. Some other package -- hypothetically emmedians -- would be great to have.

Something akin to the first approach could be implemented in emmeans, but it would require a lot of care. First, we fit a model that fits medians or at least some sort of robust estimates of cell location. Second, use regrid() to remove the operability of any linear functions involved (regrid() will make the linfct slot equal to the identity matrix). Third, use the posterior estimates to obtain the desired marginal medians. (If not a Bayesian model, use N.sim in the regrid() call to create something akin to it.). Fourth, use emmobj() to wrap the results in an emmGrid object (or use emmeans() and then hack the resulting object to have the needed posterior estimates). Then you will have some kind of median-of-median estimates, and associated HPD intervals.