# Bootstrap hypothesis test for median of differences

I have not found any examples of bootstrap hypothesis tests for the median of differences. Hence, I'd like to suggest my approach. Question: Do you a agree that the reproducible example below would be the correct way of testing the null hypothesis that the median of differences is 0 (against the alternative hypothesis that it is larger than 0)?

In addition, I am trying to relate this to the paper Two Guidelines for Bootstrap Hypothesis Testing. This paper is different to my approach here because instead of computing p-values, it finds critical t-values corresponding to certain significance levels. Nevertheless, it seems that my approach fulfills the first guideline: Resample from $$\hat{\theta}^*-\hat{\theta}$$ (because of my transformation of the differences d to d - median(d) before doing bootstrap samples). However, I don't understand how to incorporate the second guideline: Base the test on the bootstrap distribution of $$(\hat{\theta}^*-\hat{\theta}) / \hat{\sigma}^*$$. I'd be glad about any hints.

Hypotheses

H0: median(d) = 0

H1: median(d) > 0,

where d = x1 - x2 and the values are assumed to be paired. For illustration, the data sample might look as follows, where for each id, the corresponding values of x1 and x2 represent a pair.

id     x1      x2      d
1   -0.58   -0.62   0.04
2    0.23    0.04   0.19
3   -0.79   -0.91   0.12
4    1.65    0.16   1.49
5    0.38   -0.65   1.03



Explanation of Approach

Transformation: To sample under the H0, I first transform the values of d by subtracting their median. This ensures that among the transformed values d_H0 = d - median(d) the H0: median(d) = 0 is true.

Bootstrap sampling: Then, I draw R bootstrap samples: I sample from d_H0 with replacement and compute the median for each sample, obtaining R medians of differences.

Computing p-value: The p-value is computed as percentage of cases where the R medians are larger than median(d), the median of the differences in the 1 given data sample. There is a normalization constant added (hence +1 in the numerator and the denominator).

Reproducable Example (in R)

# -------------------------------------------------
# Function to get bootstrapped statistics t_star
# -------------------------------------------------
my_boot = function(d_H0, R){

N = length(d_H0)
t_star = numeric(R)

for (i in 1:R){
t_star[i] = median(sample(d_H0, size = N, replace = TRUE))
}

return(t_star)

}

# -------------------------------------------------
# Generate sample
# -------------------------------------------------
set.seed(1)
x1 = rnorm(100) + 0.05
x2 = rnorm(100)
d = x1 - x2
t = median(d)

# -------------------------------------------------
# Adjust sample to fulfill H0: median(d) = 0
# -------------------------------------------------
d_H0 = d - t

# -------------------------------------------------
# Conduct bootstrap sampling
# -------------------------------------------------
R = 5000
t_star = my_boot(d_H0, R)

# -------------------------------------------------
# Compute p-value
# -------------------------------------------------
p = (sum(t_star > t) + 1) / (R + 1)
p # 0.03

• Are you looking specifically for a bootstrap solution? I have code for a closed-form solution to test difference in medians. Commented Apr 6, 2020 at 17:25
• Which closed form solution are you referring to? I have so far only come across / been recommended to use the Wilcoxon signed-rank test. I'm open to change my mind but I generally prefer bootstrap tests. Commented Apr 6, 2020 at 17:31
• ncbi.nlm.nih.gov/pubmed/12243307 I have R code for it I can show in an answer below. Commented Apr 6, 2020 at 17:40
• I actually did come across this paper. Yes, an answer would be highly appreciated! Commented Apr 6, 2020 at 17:43

Bonett & Price (2002) proposed a closed-form solution for this, so not requiring bootstrapping. At the very least, you could compare this to your bootstrapping to see how often that the two methods agree.

Forgive me for the verbose code; I did a lot of comments and did it step-by-step for the purposes of being readable.

# test from table 3 of b&p 2002
x1 <- c(77, 87, 88, 114, 151, 210, 219, 246, 253, 262, 296, 299, 306,
376, 428, 515, 666, 1310, 2611)
x2 <- c(59, 106, 174, 207, 219, 237, 313, 365, 458, 497, 515, 529,
557, 615, 625, 645, 973, 1065, 3215)

# sort vectors
x1 <- sort(x1)
x2 <- sort(x2)

# get medians
x1_mdn <- median(x1)
x2_mdn <- median(x2)

# stuff to calculate variance of medians
x1_n <- length(x1)
x2_n <- length(x2)

x1_aj <- round((x1_n + 1) / 2 - x1_n ^ (1 / 2))
x2_aj <- round((x2_n + 1) / 2 - x2_n ^ (1 / 2))

z <- 1.855 # from table 1 of b&p 2002, see p. 376

# calculate variance
x1_var <- ((x1[x1_n - x1_aj + 1] - x1[x1_aj]) / (2 * z)) ^ 2
x2_var <- ((x2[x2_n - x2_aj + 1] - x2[x2_aj]) / (2 * z)) ^ 2

# contrast coefficients, such that its median(d) - median(dg)
x1_cj <- 1
x2_cj <- -1

# median difference
mdn_diff <- x1_mdn * x1_cj + x2_mdn * x2_cj

# standard error
mdn_diff_se <- (((x1_cj ^ 2) * x1_var) + ((x2_cj ^ 2) * x2_var)) ^ (1 / 2)

# 95% confidence interval
lb <- mdn_diff - 1.96 * mdn_diff_se
ub <- mdn_diff + 1.96 * mdn_diff_se

# within roundng error of p. 376 of b&p 2002
paste0(mdn_diff, " [", round(lb), ", ", round(ub), "]")


Reference

Bonett, D. G., & Price, R. M. (2002). Statistical inference for a linear function of medians: Confidence intervals, hypothesis testing, and sample size requirements. Psychological Methods, 7(3), 370–383. doi:10.1037/1082-989x.7.3.370

• Thanks a lot, this is very helpful as a reference. Commented Apr 6, 2020 at 18:17
• @jollycat No worries. I may try to do a quick simulation study, comparing that above equation and comparing to yours in the post. If I get around to it tonight, I'll update question and ping you. Commented Apr 6, 2020 at 19:53
• one thing to note is that I have a directed hypothesis. So, the code from my post tests this directed (one-sided) hypothesis specifically. With the confidence intervals, I believe you can only test undirected (two-sided) hypotheses. So, the two are not directly comparable. Commented Apr 6, 2020 at 20:08