Two-Sided Randomized P-Value

For a randomization test with observed statistic $T_{\text{obs}}$, a two-sided randomized p-value measures how often the randomized statistics are at least as extreme as the observed one in absolute value:

$$p = P\left(|T| \ge |T_{\text{obs}}| \mid H_0\right).$$

This is especially useful when the randomization distribution is not symmetric.

R Template

obs <- mean(yA) - mean(yB)
y <- c(yA, yB)
nA <- length(yA)
B <- 10000

rand_stat <- replicate(B, {
  idx <- sample(seq_along(y), size = nA, replace = FALSE)
  mean(y[idx]) - mean(y[-idx])
})

p_value_two_sided <- mean(abs(rand_stat) >= abs(obs))

To switch from difference in means to difference in medians, replace mean() by median().