I will try to check Rand Wilcox's book and see if this is addressed. There may be appropriate functions that weren't included in the WRS2 package. There's another way to download his other functions (https://dornsife.usc.edu/labs/rwilcox/software/ ), but there's no documentation, so you need his book to understand what function you should use for what.

For a trimmed means analysis, I'm wondering if it might make sense to use bootstrapping to determine the confidence interval for the trimmed mean for the individual groups or interaction groups, and use that as a post-hoc.

I wonder, though, if trimmed means are what you are interested in. You might consider a non-parametric approach with aligned ranks anova (although the post-hoc tests for interaction are a little different with this.)? Or transforming your original data? Or using a generalized linear model that would be appropriate for the kind of data you have?

This question was already on ResearchGate:

https://www.researchgate.net/post/3-way-Robust-ANOVA-post-hoc-test-in-R

Since you have a 2x2x2 design, I was wondering why you would need a post hoc test.

Salvatore S. Mangiafico Thank you for your reply and tips.

I opted to go with trimmed means over an aligned ranks anova because of the many outliers (1.5xIQR) my data has, besides being heteroskedastic and non-normal. Transforming data is often recommended against in my field (I'm researching stereotypes on female leaders), as it makes it harder to draw conclusions from the data. Would you agree with these choices?

I'm always wary about removing data without good reason, which is essentially what determining the trimmed mean does. The whole idea of outliers is often misused. Many things in the world are not normally or symmetrically distributed. If these measurements are naturally, say, right skewed, then we might call the fewer, higher values "outliers", even though they are legitimate observations that fit perfectly well with the real distribution of the measurements.

A good approach in these cases is to use a generalized linear model that matches the kind of data being measured.

Back to WRS2: there is a post-hoc test for 2-way trimmed means, mcp2atm(). If you don't have a three-way interaction, then looking at the post-hoc analyses for the 2-way interactions would work. Using this function is a little unintuitive: you have to print the contrasts from the result to understand what the results refer to.

I'll put this here for convenience: There is a function in the original WRS package to compute the confidence interval of a trimmed mean. The code is at the following link. [The package doesn't need to be installed, but it will need to pull some functions from WRS2 to work in that case]. https://rdrr.io/rforge/WRS/man/trimpb.html .

I checked in the Rand Wilcox book, Robust Estimation and Hypothesis Testing, 4th edition. There is some discussion about post-hoc comparisons after a 3-way analysis with trimmed mean. Unfortunately, I can't quite figure out how to use the functions. Maybe it's just me. The functions can all be downloaded with the R call source("https://dornsife.usc.edu/assets/sites/239/docs/Rallfun-v38.txt") . From there, the relevant functions are bbbmcp(), which might rely on the output from con3way() (?), and bbbtrim(). I'd be willing to send relevant pages from the text in private message, if that's helpful.

Otherwise, for a trimmed mean analysis, I'd be tempted to simply generate the confidence intervals for the groups of interest and consider non-overlapping confidence intervals to be statistically different. This isn't the most savvy approach, but it is easy to interpret and easy to plot. I shared a function in a previous response that can be used to do this, or I can share a easier method to do this for multiple groups more succinctly in private message.