# Is there any equivalent median values to Estimated Marginal Means?

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

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.

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