Inconsistent results with median test (Coin-package) R

Question

I am trying to calculate the difference between medians with a 95%CI for continuous variables of independent samples. I need to do this for group A versus group B and group A versus group C. I made two datasets, sample 1 and sample 2:

My descriptives are:

• Sample 1

  • Group A, var1: median 1.5 (20 subjects, #10 ="1" and #11="2")
  • Group B, var1: median 3 (35 subjects, #18 ="3")

• Sample 2

  • Group A, var1: median 1.5 (20 subjects, #10 ="1" and #11="2")
  • Group C, var1: median 2 (18 subjects, #9 ="2" and #10="2")

Observations:

• Group A, var1: 1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,6,8,8
• Group B, var1: 1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,4,4,4,4,5,5,6,6,7,7,7,8,9,9,18,20
• Group C, var1: 1,1,1,1,1,1,2,2,2,2,3,4,5,6,7,7,7,20

I tried to use the median_test from the "Coin" package, but it returns a "difference in location" of 0.999 for sample 1, and 1.0 for sample 2. The median difference is 1.5 (sample 1) and 2.5(sample 2) though, and I don't know how to make these results match.

For sample 1

install.packages("coin"') ; library("coin")

median_test(formula = var1 ~ Group, data = dat2, ties.method = "average-scores", mid.score = "0.5",conf.int = TRUE, conf.level = 0.95)

Output

Asymptotic Two-Sample Brown-Mood Median Test

data: var1 by Group(B, A)

Z = 2.3951, p-value = 0.01662

alternative hypothesis: true mu is not equal to 0

95 percent confidence interval: -5.996509e-05 2.000034e+00

sample estimates: difference in location 0.9999581

For sample 2

median_test(formula = Opnameduur ~ Sedatie, data = dat, ties.method = "average-scores", mid.score = "0.5", conf.int = TRUE, conf.level = 0.95)

Output

Asymptotic Two-Sample Brown-Mood Median Test

data: Opnameduur by Sedatie (morfine, ketamine)

Z = 1.2606, p-value = 0.2074

alternative hypothesis: true mu is not equal to 0

95 percent confidence interval: 2.775156e-07 3.999948e+00

sample estimates: difference in location 1.000018

I have tried adjusting the ties.method argument to "average-scores", and adjusted mid.score to "0.5", but this does not change the different estimate.

My guess is that it has to do with the mid-ranks vs average-score calculation of the medians.

My questions are:

1. Is the "difference in location" is really the difference between the medians?
2. Is there a way to adjust this test so it uses the correct median? Or a different test that works better?
3. Or is there a way to adjust the way the median is calculated in the descriptives (with some kind of mid-ranks argument)?

Thank you!

Answer 1

Because you are interested in comparing medians with data that is essentially discrete with many ties, you are going to have some difficulties constructing pretty confidence intervals.

One approach to calculating confidence intervals around median-like differences in location is to use the Hodges-Lehmann approaches. Something like this approach is followed by wilcox.test, which is commonly used by R users. This approach is appropriate even if it is not reporting differences in medians per se. Consult the documentation of the function for details.

The wilcox.exact function also produces confidence intervals about median-like differences in location.

For some specifics on the Hodges-Lehmann approach, see: Link

Another approach is to use bootstrapping to construct confidence intervals about the actual difference in medians. With discrete values, this can come out a little clunky, but directly addresses the difference in medians.

An example of bootstrapped confidence intervals is shown below.

Explicit answers to questions:

1. Is the "difference in location" is really the difference between the medians? As you have determined, no, that's not what the median_test function reports.

2. Is there a way to adjust this test so it uses the correct median? I would suspect not, as it is not what the test does.

2.b. Or a different test that works better? The code for bootstrapped confidence interval for the difference in medians, using the boot package below, reports confidence intervals for the simple difference in medians. I'm not sure this is a "better" approach than that used by median_test, wilcox.test, or wilcox.exact, but it might be simpler to explain or present. An alternative is to use the confidence intervals from wilcox.exact, as the documentation gives an explicit reference for the procedure (Bauer 1972). Code below.

3. Or is there a way to adjust the way the median is calculated in the descriptives (with some kind of mid-ranks argument)? I would imagine not.

Code:

Group.A=c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,6,8,8)  
Group.B=c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,4,4,4,
          4,5,5,6,6,7,7,7,8,9,9,18,20)
Group.C=c(1,1,1,1,1,1,2,2,2,2,3,4,5,6,7,7,7,20)

var1 = c(Group.A, Group.B, Group.C)

Group = c(rep("A", length(Group.A)),
          rep("B", length(Group.B)),
          rep("C", length(Group.C)))

data = data.frame(var1, Group)

dat2 = data[data$Group=="A" | data$Group=="B",]
dat2$Group=factor(dat2$Group)

dat3 = data[data$Group=="A" | data$Group=="C",]
dat3$Group=factor(dat3$Group)

dat4 = data[data$Group=="B" | data$Group=="C",]
dat4$Group=factor(dat4$Group)

##########################################

if(!require(exactRankTests)){install.packages("exactRankTests")}
library(exactRankTests)
wilcox.exact(var1 ~ Group, data=dat2, conf.int=T)
wilcox.exact(var1 ~ Group, data=dat3, conf.int=T)
wilcox.exact(var1 ~ Group, data=dat4, conf.int=T)

###########################################

### Adopted from http://rcompanion.org/handbook/G_10.html

if(!require(boot)){install.packages("boot")}
library(boot)

Function = function(input, index){
                    Input = input[index,]
                    M1 = median(Input$var1[Input$Group==levels(Input$Group)[1]])
                 M2 = median(Input$var1[Input$Group==levels(Input$Group)[2]])
                 Stat = M2-M1
                 return(Stat)}

Boot = boot(dat2,
             Function,
             R=5000)
mean(Boot$t[,1])
hist(Boot$t[,1])
boot.ci(Boot,
        conf = 0.95,
        type = c("perc"))

Boot = boot(dat3,
             Function,
             R=5000)
mean(Boot$t[,1])
hist(Boot$t[,1])
boot.ci(Boot,
        conf = 0.95,
        type = c("perc"))

Boot = boot(dat4,
             Function,
             R=5000)
mean(Boot$t[,1])
hist(Boot$t[,1])
boot.ci(Boot,
        conf = 0.95,
        type = c("perc"))