Confidence interval for median - which is more appropriate bootstrap or binom/exact/SAS method? Why are results different?

Question

I am unsure which method to use for calculating confidence interval for median values. I know the data set is small (n = 30). I've read these discussions which suggest both "bootstrap" and "binomial/SAS/exact" methods are feasible:

• Confidence intervals for median
• Confidence interval for median
• https://statisticsbyjim.com/hypothesis-testing/bootstrapping/

Can anyone provide insights on which one is more appropriate for the example data set and also why the CI ranges for the median is so different between the two methods for these two data sets I have below using the "DescTools" library in R? I have other data sets (example data3) where the results between the two methods are similar.

library(DescTools)
data1 = c(8,    7,  8,  9.5,    1,  20, 8,  7.5,    3,  20.5,   2.5,    5.5,    15.5,   2,  4,  1,
          17,   2,  3.5,    8.5,    8.5,    2.5,    11, 4,  10.5,   7.5,    12, 5,  16.5,   8.5)
data2 = c(7.1,  32.0,   3.8,    1.6,    19.6,   6.0,    7.2,    14.9,   0,  2.0,    5.7,
          19.4, 13.1,   15.5,   11.3,   9.6,    13.9,   5.6,    12.6,   1.0,    1.9,
          8.1,  15.9,   0.8,    6.1,    8.1,    18.0,   4.6,    5.5,    15.6)

data3 = c(16.1, 10.4,   0.5,    12.2,   7.2,    1.7,    21.6,   6.3,    0.8,    3.2,    12.6,   20.0,   3.4, 7.3,   3.5,
          7.5,  15.8, 4.7, 8.3, 11.9,   1.6,    9.0, 8.6,   11.7,   8.1, 5.8, 3.3,  7.9,    7.0,    8.5)


medianCI_Bootstrap_dF1 = MedianCI(data1, na.rm = TRUE, method = "boot")
medianCI_Binom_dF1 = MedianCI(data1, na.rm = TRUE, method = "exact")
medianCI_Bootstrap_dF1 
medianCI_Binom_dF1 

medianCI_Bootstrap_dF2 = MedianCI(data2, na.rm = TRUE, method = "boot")
medianCI_Binom_dF2 = MedianCI(data2, na.rm = TRUE, method = "exact")
medianCI_Bootstrap_dF2 
medianCI_Binom_dF2

medianCI_Bootstrap_dF3 = MedianCI(data3, na.rm = TRUE, method = "boot")
medianCI_Binom_dF3 = MedianCI(data3, na.rm = TRUE, method = "exact")
medianCI_Bootstrap_dF3 
medianCI_Binom_dF3

Answer 1

To begin, I confess I'm not familiar with either the SAS/binomial method or (just from your bit of R code) the exact type of bootstrap CI you are using.

I notice from your plot that the median is not near the middle of either of these types of CIs.

There are various correct methods for making CIs, based on different assumptions and criteria, so both kinds of results might be 'correct'. [The ultimate goal is to get a style of confidence limits that bracket the true value of the population median $\eta$ 95% of the time over the long run. Any style that does that is OK.]

A very elementary quantile bootstrap of your first sample in R gives the result $(4.5, 8.5):$

x1 = c(8, 7, 8, 9.5, 1, 20, 8, 7.5, 3, 20.5, 2.5, 5.5, 
       15.5, 2, 4, 1, 17, 2, 3.5, 8.5, 8.5, 2.5, 11, 4, 
       10.5, 7.5, 12, 5, 16.5, 8.5)
n = length(x1)
set.seed(2020)
h.re = replicate(10000, median(sample(x1, n, rep=T)))
quantile(h.re, c(.025,.975))
 2.5% 97.5% 
  4.5   8.5 

From a boxplot and a stripchart of your data the CI $(4.5,8.5)$ [or $(4,8.5)]$ seems more reasonable (to me) than the interval $(7,11).$

par(mfrow=c(2,1))
 boxplot(x1, horizontal=T, col="skyblue2")
  abline(v = c(4.5, 8.5), col="red")
 stripchart(x1, meth="stack", pch=20)
  abline(v = c(4.5, 8.5), col="red")
par(mfrow=c(1,1))

Getting any kind of nonparametric CI from these data may be problematic because you have only $n = 30$ observations with only 21 uniquely different values.

length(x1);  length(unique(x1))
[1] 30
[1] 21

Notes: Three other possibilities come to mind.

(1) Another style of 95% bootstrap CI, gives the same result $(7,11)$ as your bootstrap CI.

h.obs = median(x1)
d.re = replicate(5000, median(sample(x1, n, rep=T)) - h.obs)
UL = quantile(d.re, c(.975,.025))
h.obs - UL
97.5%  2.5% 
    7    11 

The rationale is that if we knew the distribution of the the difference $D$ between sample median $H$ and the population median $\eta,$ we could find values $L$ and $U$ that cut probability $0.025$ from the lower and upper tails, respectively of that distribution. Thus, $P(L \le D = H - \eta \le U) = 0.95.$ Then we could pivot to get the 95% CI $(H-U, H-L).$

Not knowing actual values of $H$ and $L,$ we bootstrap $D,$ temporarily using the observed median $H$ as a proxy for unknown $\eta.$ Finally, take quantiles of the bootstrap distribution of $D$ to get approximate values $L^*$ and $U^*$ of $L$ and $U,$ respectively.

(2) Contemplating that your data might be exponential with mean $\mu,$ so that $\bar X/\mu \sim \mathsf{Gamma}(30,30),$ a 95% CI for $\mu$ would be $(5.8, 11.9).$ Then a 95% CI for the median $\eta$ would be $(4.0,8.2)$ [not much different from your binomial method.]

mean(x1)/qgamma(c(.975,.025), 30,30)
[1]  5.762466 11.857195
log(2)*mean(x1)/qgamma(c(.975,.025), 30,30)
[1] 3.994237 8.218782

(3) Jittering the data slightly to break ties, so that the nonparametric 95% CI from Wilcoxon's signed rank procedure will work, I got the interval $(5.5, 9.8).$

wilcox.test(x1+runif(30,-.1,.1), conf.int=T)$conf.int
[1] 5.479651 9.751046
attr(,"conf.level")
[1] 0.95

Two more runs:

wilcox.test(x1+runif(30,-.1,.1), conf.int=T)$conf.int
[1] 5.479582 9.763138
attr(,"conf.level")
[1] 0.95

wilcox.test(x1+runif(30,-.1,.1), conf.int=T)$conf.int
[1] 5.496640 9.752084
attr(,"conf.level")
[1] 0.95