Why does median() differ from quantile()["50%"] in R and doesn't honor the type? Is this an estimator or sample quantile?

Question

I just found, that when I want a specific type of quantile in R, say the number 3, it doesn't agree with median calculated in the "traditional way" in R.

> set.seed(100); x <- rnorm(100)
> sapply(1:9, function(t) setNames(round(quantile(x, type=t)["50%"], 7), paste("type:", t)))
   type: 1    type: 2    type: 3    type: 4    type: 5    type: 6    type: 7    type: 8    type: 9 
-0.0688440 -0.0594199 -0.0688440 -0.0688440 -0.0594199 -0.0594199 -0.0594199 -0.0594199 -0.0594199 

> median(x)
[1] -0.0594199  # doesn't agree with types 1, 3, 4

I was told, that median() in R is sample median and uses the textbook definition, while quantile() doesn't support the textbook definition in neither way, but rather uses the 1-9 approach described in some book.

Yet the documentation of quantiles also says it's "sample quantile".

So - shouldn't "sample 50% quantile" be equal to "sample median"?

I found that the median() is used also in mean() for case we want the 50% trimmed mean (=median).

So, if I wanted to replicate calculations done by others in other statistical software, like SAS, I have to stick to quantile(), where I can set the type, but I cannot use the median(), which is very often used by others, who don't care of that, so my calculations may diverge somewhere.

Is there any reason behind that in R?

Apart from R, is there any reason for making the sample median different than sample quantile(50%)?

---

Part of the mean() source code in R:

mean.default <- function(x, trim = 0, na.rm = FALSE, ...)
{
...
    if(trim > 0 && n) {
...
    if(trim >= 0.5) return(stats::median(x, na.rm=FALSE))
...
}

The median() code:

median.default <- function(x, na.rm = FALSE, ...)
{
...
    half <- (n + 1L) %/% 2L
    if(n %% 2L == 1L) sort(x, partial = half)[half]
    else mean(sort(x, partial = half + 0L:1L)[half + 0L:1L])
}

Source of the quantile(): https://github.com/wch/r-source/blob/trunk/src/library/stats/R/quantile.R

Note to the Moderators: I don't ask "how to do something in R" and the manuals don't answer that general question about any potential reasons for doing that.

Answer 1

Whenever there is an even number of unique values, there is an infinite number of correctly calculated medians. If we consider the values 1 to 4, every value between 2 and 3 (excluding the limits) will split the values into to equally sized halves.

Usually the exact choice of algorithm is of minor importance as sample medians are random numbers and more often then not, different algorithms will yield very similar results. Obviously exceptions can be found and it would indeed be nice to have some universally accepted standard so each statistic program could produce identical results.

Unfortunately, differing results from different software packages are a common problem we have to deal with.

Bergman, Ludbrook and Sporen: Different Outcomes of the Wilcoxon—Mann—Whitney Test from Different Statistics Packages, February 2000 The American Statistician 54(1):72-77 ( DOI: 10.1080/00031305.2000.10474513 ) found largely differing p values from 11 statistical software packages in Mann-Whitney-U-Tests performed on data with ties and many internet pages are concerned with Type I, II and III sums of squares because people try to get identical ANOVA results from SPSS and R.

A worthwhile read in the context of this question is Jorine Putter, Liza Faber: Quartiles within SAS, 2012 downloadable from https://www.lexjansen.com/phuse/2012/pp/PP16.pdf

The authors describe the 5 different ways to compute quantiles that SAS provides and acknowledge the problems possibly associated with this choice:

Many times during the reporting of a study, programmers blindly report whichever statistics are generated by default by the specific SAS® procedure ...

Unbeknownst to the programmer, the reported statistics may not accurately reflect what the statistician is expecting to see.

The purpose of this paper is to educate programmers on the different methods for calculating Q1 and Q3, and ensuring the statistician has clearly documented the appropriate method to use.

Besides ample explanation of the SAS ways of doing quantiles the authors do not only mention the 9 different ways to compute quantiles that R provides. In their table 3 we find which parameter in SAS equals which type in R. From that we can take that

type = 1 equals QNTLDEF = 3
type = 2 equals QNTLDEF = 5
type = 3 equals QNTLDEF = 2
type = 4 equals QNTLDEF = 1
type = 6 equals QNTLDEF = 4

R's standard for the quantile function is type = 7, which apparantely has no perfext match in SAS. From the linked paper:

The default method used by SAS, is QNTLDEF=5.

I do not believe that it will be feasible to ask the programmers of SAS's or R's original function why they opted for a particular definition of a quantile for which function. The linked paper states that R's predefined type = 7 equals the definition Excel uses. That is not a great statistical reason but might be a very practical one.

A first helpful step towards equal results in both packages might be to redefine R's median function along the lines of

median <- function(x, na.rm = FALSE, type = 2, ...) 
              quantile(x, probs = .5, na.rm = na.rm, type = type, ...)

This will of course only change the results of median within that scope/environment.