Bootstrap Anova and posthoc comparisons for non-normal data

Question

I have a dataset like the following (n=1400):

register    country PC1
CMT BD  0.528902409041985
CMT IN  0.659599336404661
CMT LK  0.424746884028921
CMT PK  0.617481735398022
CMT UK  0.432241778651171
CMT US  0.520006978931032
TWT BD  -0.120412754435259
TWT IN  -0.775416939396557
TWT LK  -0.331060813776788
TWT PK  -0.0476004644598422
TWT UK  -0.751168065821314
TWT US  -0.861747850448701
TXM BD  -0.899207300872416
TXM IN  -1.90230790510253
TXM LK  0.257287440181
TXM PK  -1.3102770881823
WBF BD  -0.38312607807368
WBF IN  -1.4048106311512
WBF LK  -0.238559559698407
WBF PK  0.0249239934526432
WBF UK  -0.467017637887557
WBF US  -0.423802534509881
WBS BD  1.53739431443881
WBS IN  0.275786018712733
WBS LK  1.32988601584956
WBS PK  1.68224760320901
WBS UK  1.6017172088108
WBS US  1.34625059689434

I am interested in ANOVA and if significant groups comparisons using emmeans package in R. afex::check_homogeneiety throws unequal variance warning for PC1. The residuals are not normally distributed as per afex::check_noarmality. See also qqplot below):

Which means that I cannot use anova() and emmeans in one go like this:

library(emmeans)
m_dims <- lm(PC1 ~ register*country, data = dims)
m_dims
anova(m_dims)
em_dims <- emmeans(m_dims, pairwise ~ country | register)

See the sample output:

Analysis of Variance Table

Response: PC1
                   Df Sum Sq Mean Sq  F value    Pr(>F)    
register            4 776.63 194.157 468.6266 < 2.2e-16 ***
country             5  20.55   4.111   9.9222 2.452e-09 ***
register:country   18  33.39   1.855   4.4769 1.411e-09 ***
Residuals        1372 568.43   0.414                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

And for emmeans()$contrasts (here only first two levels of register)

register = CMT:
 contrast estimate    SE   df t.ratio p.value
 BD - IN  -0.29759 0.129 1372  -2.312  0.1899
 BD - LK   0.42673 0.129 1372   3.315  0.0121
 BD - PK  -0.05777 0.129 1372  -0.449  0.9977
 BD - UK   0.07512 0.129 1372   0.584  0.9921
 BD - US   0.13296 0.129 1372   1.033  0.9069
 IN - LK   0.72432 0.129 1372   5.626  <.0001
 IN - PK   0.23981 0.129 1372   1.863  0.4257
 IN - UK   0.37271 0.129 1372   2.895  0.0445
 IN - US   0.43055 0.129 1372   3.344  0.0109
 LK - PK  -0.48451 0.129 1372  -3.764  0.0024
 LK - UK  -0.35161 0.129 1372  -2.731  0.0698
 LK - US  -0.29377 0.129 1372  -2.282  0.2021
 PK - UK   0.13290 0.129 1372   1.032  0.9071
 PK - US   0.19073 0.129 1372   1.482  0.6762
 UK - US   0.05783 0.129 1372   0.449  0.9977

register = TWT:
 contrast estimate    SE   df t.ratio p.value
 BD - IN  -0.38951 0.129 1372  -3.026  0.0303
 BD - LK  -0.13149 0.129 1372  -1.021  0.9109
 BD - PK  -0.12868 0.129 1372  -1.000  0.9182
 BD - UK   0.10248 0.129 1372   0.796  0.9682
 BD - US   0.34901 0.129 1372   2.711  0.0737
 IN - LK   0.25802 0.129 1372   2.004  0.3403
 IN - PK   0.26083 0.129 1372   2.026  0.3279
 IN - UK   0.49199 0.129 1372   3.822  0.0019
 IN - US   0.73852 0.129 1372   5.737  <.0001
 LK - PK   0.00281 0.129 1372   0.022  1.0000
 LK - UK   0.23397 0.129 1372   1.817  0.4547
 LK - US   0.48050 0.129 1372   3.733  0.0027
 PK - UK   0.23116 0.129 1372   1.796  0.4689
 PK - US   0.47769 0.129 1372   3.711  0.0029
 UK - US   0.24653 0.129 1372   1.915  0.3932

So I decided to use bootstrapping to resample my data, apply anova() and emmeans() on each sample and calculate the usual statistics for register, country, register*country: p-values, F statistic, degrees of freedom, R-sq etc. from anova() output, and pair wise comparisons of each country level (PK, UK, US, LK, IN, BD) within each register level (CMT, TWT, TXM, WBF, WBS). As per my very limited understanding of bootstrapping, I thought of averaging the resulting distributions to get each statistic, e.g. median p-value from all 1000 or more anova() outputs from my data samples. My questions:

1. Am I correct to assume that the p-value or any statistic obtained this way is a robust alternative to the one time output of anova() (or the subsequent emmeans()) as I showed above?
2. If my assumption is not correct, how should I proceed to apply bootstrap in this scenario? Before writing this post, I have consulted various blog posts and searched for ready-made solutions/functions in R but could not find anything suitable or convincing. Some references An R script for bootstrap ANOVA and post hoc comparisons. (I changed lsmeans to emmeans but it outputs same p-value for each post-hoc comparison which I do not understand why, so I left it). Bootstrap resampling with tidymodels Bootstrap Anova. Bootstrap followup contrasts (no ANOVA bootstrapping).

Answer 1

You can use emmeans() directly with bootstrap results by inputting them into emmobj() as posterior samples. Here is a simple example with three independent samples of somewhat unequal size generated from Cauchy distributions.

library(emmeans)

set.seed(08.28)

# Cauchy samples with medians 12, 15, 14 resp.
dat = list(rt(50, 1) + 12,  rt(60, 1) + 15, rt(50, 1) + 14)

boot = replicate (200, sapply(dat, function(x) {
    y = sample(x, length(x), replace = TRUE)
    median(y)
}))

# each column is a set of 3 medians. Need them to be the rows
boot = t(boot)

# frequentist estimates -- won't actually be used
means = apply(boot, 2, mean)
V = cov(boot)

# marginal means
EMM = emmobj(means, V, levels = LETTERS[1:3], post.beta = boot)

summary(EMM)
##  level estimate lower.HPD upper.HPD
##  A         12.1      11.6      12.5
##  B         15.0      14.7      15.3
##  C         13.9      13.7      14.6
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

summary(pairs(EMM))
##  contrast estimate lower.HPD upper.HPD
##  A - B       -2.96    -3.445     -2.38
##  A - C       -1.92    -2.730     -1.35
##  B - C        1.05     0.369      1.54
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

^{Created on 2022-08-28 by the reprex package (v2.0.1)}

Note the medians of the medians do come close to the true values used in the simulation; and of course so do the differences.

Note I used the median as the estimator in the bootstrap. I strongly encourage you to not use the mean there, which is a poor estimator when you have heavy tails.