Brunner-Munzel Test

This is a generic function that performs a generalized Brunner-Munzel test in a fashion similar to t.test.

Usage

# Default S3 method
brunner_munzel(
  x,
  y,
  paired = FALSE,
  alternative = c("two.sided", "less", "greater", "equivalence", "minimal.effect"),
  mu = 0.5,
  alpha = 0.05,
  scale = c("probability", "difference", "logodds", "odds"),
  test_method = c("t", "logit", "perm", "perm_logit"),
  R = 10000,
  p_method = NULL,
  perm = "deprecated",
  max_n_perm = "deprecated",
  ...
)

# S3 method for class 'formula'
brunner_munzel(formula, data, subset, na.action, ...)

Arguments

x

a (non-empty) numeric vector of data values.

y

an optional (non-empty) numeric vector of data values.

paired

a logical indicating whether you want a paired test.

alternative

a character string specifying the alternative hypothesis, must be one of:

"two.sided" (default): true relative effect is not equal to mu
"less": true relative effect is less than mu
"greater": true relative effect is greater than mu
"equivalence": true relative effect is between the lower and upper bounds specified in mu
"minimal.effect": true relative effect is less than the lower bound or greater than the upper bound specified in mu

mu

a number or vector specifying the null hypothesis value(s):

For standard alternatives ("two.sided", "less", "greater"): a single value representing the hypothesized relative effect (default = 0.5, i.e., stochastic equality).
For "equivalence" or "minimal.effect": two values representing the lower and upper bounds for the relative effect. Values must be between 0 and 1.

Note: mu is always specified on the probability scale regardless of the scale argument. The scale argument only affects how results are reported; it does not change the hypothesis being tested.

alpha

alpha level (default = 0.05)

scale

a character string specifying the scale for the reported estimate, standard error, confidence interval, and null value. The test itself always operates on the probability scale internally; this argument only transforms the output.

"probability": (default): $p = P(X > Y) + 0.5 \cdot P(X = Y)$, range $[0, 1]$, null at stochastic equality = 0.5
"difference": $P(X > Y) - P(X < Y) = 2p - 1$, range $[-1, 1]$, null = 0
"logodds": $\log[p / (1 - p)]$, range $(-\infty, \infty)$, null = 0
"odds": $p / (1 - p)$, range $(0, \infty)$, null = 1

test_method

a character string specifying the test method to use:

"t" (default): approximate t-distribution with Satterthwaite-Welch degrees of freedom
"logit": logit transformation for range-preserving confidence intervals
"perm": studentized permutation test (recommended for the standard null of 0.5 when the smaller group has fewer than 30 observations, and for paired designs)
"perm_logit": studentized permutation test on the logit scale, with range-preserving confidence intervals (recommended for equivalence, minimal effect, and other tests against a null value other than 0.5)

See "Choosing a test method" in Details.

R

the number of permutations for the permutation test (default is 10000). Only used when test_method = "perm" or "perm_logit".

p_method

a character string specifying the method for computing permutation p-values. Only used when test_method = "perm" or "perm_logit":

NULL (default): Automatically selects "exact" for exact permutation tests (when all permutations are enumerated) and "plusone" for randomization tests (when permutations are sampled with replacement).
"exact": Uses b/R where b is the count of permutation statistics at least as extreme as observed. Appropriate when all permutations are enumerated.
"plusone": Uses (b+1)/(R+1), which guarantees p > 0 and provides exact Type I error control for randomization tests where permutations are sampled with replacement (Phipson & Smyth, 2010).

perm

Use test_method = "perm" instead.

max_n_perm

Use R instead.

...

further arguments to be passed to or from methods.

formula

a formula of the form lhs ~ rhs where lhs is a numeric variable giving the data values and rhs either 1 for a one-sample test or a factor with two levels giving the corresponding groups. For paired tests, use the default method with x and y vectors instead of the formula method.

data

an optional matrix or data frame (or similar: see model.frame) containing the variables in the formula formula. By default the variables are taken from environment(formula).

subset

an optional vector specifying a subset of observations to be used.

na.action

a function which indicates what should happen when the data contain NAs. Defaults to getOption("na.action").

Value

A list with class "htest" containing the following components:

"statistic": the value of the test statistic.
"parameter": the degrees of freedom for the test statistic.
"p.value": the p-value for the test.
"conf.int": a confidence interval for the relative effect appropriate to the specified alternative hypothesis.
"estimate": the estimated relative effect.
"null.value": the specified hypothesized value of the relative effect.
"stderr": the standard error of the relative effect.
"alternative": a character string describing the alternative hypothesis.
"method": a character string indicating what type of test was performed.
"data.name": a character string giving the name(s) of the data.

Details

This function is made to provide a test of stochastic equality between two samples (paired or independent), and is referred to as the Brunner-Munzel test.

This tests the hypothesis that the relative effect, discussed below, is equal to the null value (default is mu = 0.5).

The estimate of the relative effect, which can be considered as value similar to the probability of superiority, refers to the following:

$$\hat p = P(X>Y) + \frac{1}{2} \cdot P(X=Y)$$

Note, for paired samples, this does not refer to the probability of an increase/decrease in paired sample but rather the probability that a randomly sampled value of X is greater than a randomly sampled value of Y. This is also referred to as the "relative" effect in the literature. Therefore, the results will differ from the concordance probability provided by the ses_calc function.

The brunner_munzel function is based on the npar.t.test and npar.t.test.paired functions within the nparcomp package (Konietschke et al. 2015).

Test Methods

Four test methods are available:

"t": The default method uses a t-distribution approximation with Satterthwaite-Welch degrees of freedom. This is appropriate for moderate to large sample sizes.
"logit": Uses a logit transformation to produce range-preserving confidence intervals that are guaranteed to stay within [0, 1]. It is conservative in small samples.
"perm_logit": A studentized permutation test carried out on the logit scale. The observed statistic is $(\mathrm{logit}(\hat p) - \mathrm{logit}(p_0)) / SE_{logit}$, with $SE_{logit} = SE / (\hat p (1 - \hat p))$, and the permutation distribution is the same statistic computed for each permuted data set. The confidence interval inverts this test and is back-transformed, so it always stays strictly inside (0, 1) and always agrees with the p-value. Estimates of exactly 0 or 1 are moved half a step inside the unit interval (by $1/(2 n_1 n_2)$ for two samples, $1/(2 n^2)$ for paired samples) before transforming. Null values must be strictly between 0 and 1. The R and p_method arguments work as for "perm".
"perm": A studentized permutation test following Neubert & Brunner (2007). This method is recommended for the standard null of 0.5 in small samples, as it provides better control of Type I error rates in these situations. Note: when permutations are sampled (rather than fully enumerated) and p_method = "exact" (b/R) is requested, the p-value can be 0 if no sampled permuted statistic is as extreme as the observed value. The default plusone method ((b+1)/(R+1)) for sampled permutations avoids this. With full enumeration the observed arrangement is always among the permutations, so b/R is never 0.

Choosing a test method

Recommended defaults:

Hypothesis and design	Recommendation
Null of 0.5, two-sample, smaller group < 30	`"perm"`
Null of 0.5, two-sample, smaller group 30 or more	`"t"` (or `"perm"`; they agree)
Null of 0.5, paired, 5 or fewer pairs	No good option (see below); treat the design as underpowered
Null of 0.5, paired, 6 or more pairs	`"perm"`
Null of 0.5, paired, large sample with approximately uncorrelated pairs	`"t"` is adequate
Equivalence, minimal effect, or any other null value	`"perm_logit"` (most powerful) or `"logit"` (more conservative)
Estimate near 0 or 1 (for the confidence interval)	`"perm_logit"` or `"logit"`

These recommendations come from simulations with normal, heavily tied ordinal, unequal-variance, and skewed data, at $\alpha = 0.05$. Smaller $\alpha$ levels were not examined and may require larger samples for the "t" approximation to be adequate.

Equivalence, minimal effect, and other null values. The advantage of "perm" only holds for the null of 0.5: a permutation test is exact only when the groups are exchangeable, which is impossible when the true relative effect equals a bound other than 0.5. In simulations with bounds of (0.4, 0.6), (0.3, 0.7), and (0.2, 0.8) and 7 to 30 observations per group or pairs:

"t" and "perm" behaved almost identically. Equivalence tests held their level but were overly conservative, with low power (for example about 0.31 with 10 per group and bounds of (0.2, 0.8), compared with 0.53 for "logit" and 0.64 for "perm_logit"). Minimal effect tests had inflated Type I error, up to 15–18% in small samples with bounds of (0.2, 0.8) and still about 9–10% with 30 per group. The function warns when "t" or "perm" is used with a null value other than 0.5.
"logit" held its level for both equivalence (Type I error at most about 5.5%) and minimal effect tests (mostly at or below 5%), with more power than "t" or "perm".
"perm_logit" had the highest power. Its Type I error at the bounds was about 5–7% for equivalence tests and mostly about 5% for minimal effect tests, but reached 8–12% for minimal effect tests in small two-sample designs with unequal spread and wide bounds. Use "logit" when strict Type I error control matters more than power in such designs.

Estimates near 0 or 1. The "t" and "perm" intervals can collapse when the estimate is at or near 0 or 1 (their standard error goes to zero), and their coverage fell as low as 52–82% for a true relative effect of 0.95 with 15 or fewer per group. "perm_logit" kept coverage near 95% in these settings with equal spread and in paired designs, but could undercover (roughly 72–85%) with unequal spread and a true relative effect of 0.9 or more. "logit" never fell below 90% coverage but produces wide intervals.

Two-sample designs. With fewer than about 10 observations in the smaller group, "t" is liberal (Type I error of roughly 7–9% at 5 per group), whereas "perm" stays at or below the nominal level when the two distributions are identical. From about 15 observations in the smaller group, "t" and "perm" were practically equivalent in every scenario: their Type I error rates differed by at most one percentage point, and they reached different decisions for about 1% of data sets or fewer. The cut-off of 30 in the table leaves a safety margin. It applies to the smaller group, because unbalanced designs are where "t" took longest to agree.

Paired designs. The within-pair correlation matters more than the sample size. With positively correlated pairs (the usual paired design), "t" is conservative. For example, with a latent correlation of 0.8 its Type I error was about 2–4% for 7 to 50 pairs, while "perm" stayed near 5%. "perm" is therefore preferable, and more powerful, at every sample size studied (up to 50 pairs). With approximately uncorrelated pairs, "t" and "perm" agree from roughly 15–20 pairs. With 5 or fewer pairs, "perm" cannot reject at two-sided $\alpha = 0.05$: full enumeration of $2^5 = 32$ within-pair swaps gives a smallest attainable p-value of $2/32 = 0.0625$. The level of "t" is also uncertain at such small sample sizes.

The "logit" method was markedly conservative in small samples in both designs (Type I error of roughly 1–3% with 15 or fewer observations per group) and remained below nominal at 50, so it sacrifices considerable power. Its advantage is an interval that always stays within [0, 1], which is useful when the estimated relative effect is close to 0 or 1. It is not a remedy for small samples.

Limitations. Neither "t" nor "perm" is exact when the two distributions differ in spread or shape.

In two-sample designs with very small samples, both can be liberal when the smaller group is the more variable one.
In paired designs with strongly correlated pairs and unequal marginal distributions, "perm" can be slightly liberal (about 6%).
"perm" can be conservative with heavily tied data in very small samples.

Messages and warnings. The function prints messages that follow these recommendations:

For the null of 0.5, it recommends "perm" when a two-sample test with fewer than 30 observations in the smaller group uses "t" or "logit", and whenever a paired test uses "t" or "logit".
For two-sample designs it notes that a permutation test is probably unnecessary only for very large samples (smaller group > 250 and total > 500). This is deliberately conservative. Using a permutation test in large samples is never invalid, only slower.
It warns when "t" or "perm" is used for a minimal effect test or any other test against a null value other than 0.5, and suggests "perm_logit" or "logit".
When the permutation distribution is too coarse for the test to ever reject at the requested alpha (for example, 5 or fewer pairs), the function says so.

Hypothesis Testing

For the standard alternatives, the null hypothesis is that the relative effect equals mu:

"two.sided": H0: p = mu vs H1: p ≠ mu
"less": H0: p ≥ mu vs H1: p < mu
"greater": H0: p ≤ mu vs H1: p > mu

For equivalence and minimal effect testing using the two one-sided tests (TOST) procedure:

"equivalence": H0: p ≤ mu[1] OR p ≥ mu[2] vs H1: mu[1] < p < mu[2]

Tests whether the relative effect falls within the specified bounds. The p-value is the maximum of the two one-sided p-values.
"minimal.effect": H0: mu[1] < p < mu[2] vs H1: p ≤ mu[1] OR p ≥ mu[2]

Tests whether the relative effect falls outside the specified bounds. The p-value is the minimum of the two one-sided p-values.

Test Statistic and P-value Calculation

The test statistic is calculated as:

$$t = \sqrt{N} \cdot \frac{\hat{p} - p_0}{s}$$

where $N$ is the total sample size (or $n$ for paired samples), $\hat{p}$ is the estimated relative effect, $p_0$ is the null hypothesis value, and $s$ is the rank-based standard error.

For equivalence testing, two test statistics are computed:

$$t_{low} = \sqrt{N} \cdot \frac{\hat{p} - p_{low}}{s}$$ $$t_{high} = \sqrt{N} \cdot \frac{\hat{p} - p_{high}}{s}$$

where $p_{low}$ and $p_{high}$ are the lower and upper equivalence bounds. The one-sided p-values are:

$p_1$: p-value for H1: p > $p_{low}$ (from the lower bound test)
$p_2$: p-value for H1: p < $p_{high}$ (from the upper bound test)

For equivalence: $p_{TOST} = \max(p_1, p_2)$

For minimal effect: $p_{MET} = \min(1 - p_1, 1 - p_2)$. For test_method = "perm" and "perm_logit", the minimal effect one-sided p-values are instead computed by counting the opposite tails of the permutation distribution directly (with the same p_method counting rule).

Confidence Intervals for Equivalence Testing

When alternative = "equivalence" or alternative = "minimal.effect", the confidence interval is computed at the $1 - 2\alpha$ level (default: 90% CI when $\alpha = 0.05$). This follows the standard TOST procedure where the $(1 - 2\alpha) \times 100\%$ CI corresponds to two one-sided tests at level $\alpha$.

Permutation Confidence Intervals

When test_method = "perm", the confidence interval is obtained by inverting the same studentized permutation test that produces the p-value: $\hat{p} - SE \cdot q$, where $SE$ is the standard error from the original sample and $q$ is the order statistic of the permutation distribution at which the permutation p-value (with the same p_method counting rule) crosses $\alpha$. The two-sided interval uses the same absolute-value rule as the two-sided p-value. As a result, the confidence interval and p-value always agree, apart from the interval being clamped to [0, 1]. If the number of permutations is too small for the test to ever reject at the requested $\alpha$, the corresponding limit is set to 0 or 1.

When test_method = "perm_logit", the same inversion is carried out on the logit scale, $\mathrm{logit}(\hat{p}) - SE_{logit} \cdot q$, and the limits are back-transformed. No clamping is needed, and the interval and p-value always agree.

Permutation Tests with Non-0.5 Null Values

When test_method = "perm" or "perm_logit" and mu != 0.5, the permutation distribution is constructed by centering the permuted test statistics at 0.5 (the value implied by exchangeability), while the observed test statistic is centered at the hypothesized null value. This approach is justified asymptotically, because the studentized permutation distribution converges to the same limit regardless of the centering (Janssen, 1997; Neubert & Brunner, 2007). It is not exact: when the true relative effect equals a null value other than 0.5, the groups cannot be exchangeable. In small samples, how well this works depends strongly on the scale of the statistic; see "Choosing a test method" above.

Because the equivalence bounds are specified directly on the relative effect scale (i.e., as probabilities between 0 and 1), this avoids the difficulty noted by Arboretti et al. (2021, point IU.7) that arises when margins must be expressed in terms of rank transformations.

These are uncalibrated (naive) procedures. For the IU direction ("equivalence"), the procedure can be conservative when sample sizes are small or when the equivalence bounds are close to 0.5. For the UI direction ("minimal.effect"), the identity-scale "perm" procedure can instead be liberal, because its interval collapses when the estimate approaches 0 or 1 (see "Choosing a test method"). See Arboretti et al. (2021) for a detailed discussion of calibration.

References

Arboretti, R., Pesarin, F. & Salmaso, L. (2021). A unified approach to permutation testing for equivalence. Stat Methods Appl 30, 1033-1052. doi: 10.1007/s10260-020-00548-0

Brunner, E., Munzel, U. (2000). The Nonparametric Behrens-Fisher Problem: Asymptotic Theory and a Small Sample Approximation. Biometrical Journal 42, 17 -25.

Neubert, K., Brunner, E., (2006). A Studentized Permutation Test for the Nonparametric Behrens-Fisher Problem. Computational Statistics and Data Analysis.

Munzel, U., Brunner, E. (2002). An Exact Paired Rank Test. Biometrical Journal 44, 584-593.

Munzel, U., Hauschke, D. (2003). A nonparametric test for proving noninferiority in clinical trials with ordered categorical data. Pharmaceutical Statistics, 2, 31-37.

Konietschke, F., Placzek, M., Schaarschmidt, F., & Hothorn, L. A. (2015). nparcomp: an R software package for nonparametric multiple comparisons and simultaneous confidence intervals. Journal of Statistical Software 64 (2015), Nr. 9, 64(9), 1-17. http://www.jstatsoft.org/v64/i09/

Janssen, A. (1997). Studentized permutation tests for non-i.i.d. hypotheses and the generalized Behrens-Fisher problem. Statistics & Probability Letters, 36(1), 9-21.

Phipson, B., & Smyth, G. K. (2010). Permutation P-values should never be zero: calculating exact P-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology, 9(1), Article 39.

Examples

data(mtcars)
# Standard test of stochastic equality
brunner_munzel(mpg ~ am, data = mtcars)
#> Sample size in at least one group is small (< 30). The permutation test (test_method = 'perm') is recommended.
#> 
#> 	Two-sample Brunner-Munzel test
#> 
#> data:  mpg by am
#> t = -4.2653, df = 20.893, p-value = 0.0003479
#> alternative hypothesis: true relative effect is not equal to 0.5
#> 95 percent confidence interval:
#>  0.009114802 0.330966169
#> sample estimates:
#> P('0'>'1') + .5*P('0'='1') 
#>                  0.1700405 
#> 

# Test using logit transformation for range-preserving CIs
brunner_munzel(mpg ~ am, data = mtcars, test_method = "logit")
#> Sample size in at least one group is small (< 30). The permutation test (test_method = 'perm') is recommended.
#> 
#> 	Two-sample Brunner-Munzel test (logit)
#> 
#> data:  mpg by am
#> t = -2.8922, df = 20.893, p-value = 0.008748
#> alternative hypothesis: true relative effect is not equal to 0.5
#> 95 percent confidence interval:
#>  0.06147768 0.39053873
#> sample estimates:
#> P('0'>'1') + .5*P('0'='1') 
#>                  0.1700405 
#> 

# Permutation test (recommended for the null of 0.5 in small samples)
brunner_munzel(mpg ~ am, data = mtcars, test_method = "perm", R = 999)
#> 
#> 	Two-sample Brunner-Munzel Randomization test
#> 
#> data:  mpg by am
#> t-observed = -4.2653, N-permutations = 999, p-value = 0.003
#> alternative hypothesis: true relative effect is not equal to 0.5
#> 95 percent confidence interval:
#>  0.003426061 0.336654911
#> sample estimates:
#> P('0'>'1') + .5*P('0'='1') 
#>                  0.1700405 
#> 

# Test against a specific null value (logit scale recommended for nulls other than 0.5)
brunner_munzel(mpg ~ am, data = mtcars, mu = 0.3, test_method = "logit")
#> 
#> 	Two-sample Brunner-Munzel test (logit)
#> 
#> data:  mpg by am
#> t = -1.3464, df = 20.893, p-value = 0.1926
#> alternative hypothesis: true relative effect is not equal to 0.3
#> 95 percent confidence interval:
#>  0.06147768 0.39053873
#> sample estimates:
#> P('0'>'1') + .5*P('0'='1') 
#>                  0.1700405 
#> 

# Equivalence test: is the relative effect between 0.35 and 0.65?
# Permutation test on the logit scale (recommended for equivalence testing)
brunner_munzel(mpg ~ am, data = mtcars,
               alternative = "equivalence",
               mu = c(0.35, 0.65),
               test_method = "perm_logit", R = 999)
#> NOTE: Permutation-based TOST for equivalence/minimal.effect testing.
#> 
#> 	Two-sample Brunner-Munzel Randomization test (logit)
#> 
#> data:  mpg by am
#> t-observed = -1.7628, N-permutations = 999, p-value = 0.975
#> alternative hypothesis: equivalence
#> null values:
#> lower bound upper bound 
#>        0.35        0.65 
#> 90 percent confidence interval:
#>  0.07950668 0.32333005
#> sample estimates:
#> P('0'>'1') + .5*P('0'='1') 
#>                  0.1700405 
#> 

# Minimal effect test: is the relative effect outside 0.4 to 0.6?
brunner_munzel(mpg ~ am, data = mtcars,
               alternative = "minimal.effect",
               mu = c(0.4, 0.6),
               test_method = "logit")
#> 
#> 	Two-sample Brunner-Munzel test (logit)
#> 
#> data:  mpg by am
#> t = -2.1525, df = 20.893, p-value = 0.0216
#> alternative hypothesis: minimal.effect
#> null values:
#> lower bound upper bound 
#>         0.4         0.6 
#> 90 percent confidence interval:
#>  0.0738647 0.3448182
#> sample estimates:
#> P('0'>'1') + .5*P('0'='1') 
#>                  0.1700405 
#> 

# Report on the difference scale: P(X>Y) - P(X<Y)
brunner_munzel(mpg ~ am, data = mtcars, scale = "difference")
#> Sample size in at least one group is small (< 30). The permutation test (test_method = 'perm') is recommended.
#> 
#> 	Two-sample Brunner-Munzel test
#> 
#> data:  mpg by am
#> t = -4.2653, df = 20.893, p-value = 0.0003479
#> alternative hypothesis: true relative effect (difference) is not equal to 0
#> 95 percent confidence interval:
#>  -0.9817704 -0.3380677
#> sample estimates:
#> P('0'>'1') - P('0'<'1') 
#>               -0.659919 
#> 

# Report as odds
brunner_munzel(mpg ~ am, data = mtcars, scale = "odds")
#> Sample size in at least one group is small (< 30). The permutation test (test_method = 'perm') is recommended.
#> 
#> 	Two-sample Brunner-Munzel test
#> 
#> data:  mpg by am
#> t = -4.2653, df = 20.893, p-value = 0.0003479
#> alternative hypothesis: true relative effect (odds) is not equal to 1
#> 95 percent confidence interval:
#>  0.009198646 0.494692726
#> sample estimates:
#> odds(P('0'>'1') + .5*P('0'='1')) 
#>                         0.204878 
#>