4
$\begingroup$

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.

$\endgroup$
3
  • 2
    $\begingroup$ The defitinion of a median is clear for odd vector lengths, but there are different options for even vector length. median gives you a simple way, quantile gives you a bunch of choices. Each of the types returns a valid median. Usually the results are similar, but there would not be any sense in offering choices if they all came to the same results. If you want the same results as SAS you need to figure out, which way SAS computes medians. What exactly is the problem then? $\endgroup$
    – Bernhard
    Commented Jan 11, 2021 at 10:47
  • $\begingroup$ The problem is, that if I want to be compliant with SAS, I can choose appropriate type. But the default, base R median function doesn't respect this and is used everywhere (also in the mean() for trimming), for winsorization in other packages, by many people. I can relate, who wants to write quantile(...)["50%"] rather than just median(). This makes practically all the analyses not only noncompliant, but also incurable. So I was curious was there any reason to not base median() on quantile() making it flexible, rather than implementing the median entirely from scratch and handling only type 7. $\endgroup$
    – StatNovice
    Commented Jan 11, 2021 at 11:55
  • 1
    $\begingroup$ So the problem is not in that I have no choice. The problem is that the function of major importance is not flexible and fixed at a single algorithm. And because it's massively widespread, it affects a lot of code, which causes then a lot of discrepancies, if there is a team of both SAS and R programmers, like in finances, where I work, who need to explain the differences all the times. I was wondering is there is any statistical justification for that? I assumed someone did this deliberately. $\endgroup$
    – StatNovice
    Commented Jan 11, 2021 at 11:57

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ This does not really answer OP's question. Why does the median function use a different type than quantile? One would except both to be equal, median just being shorthand for quantile(0.5), why is this not so, why this inconsistency? $\endgroup$ Commented Jan 12, 2021 at 7:13
  • $\begingroup$ No, it does not answer that legitimate question at all. And I stated that I do not expect a good answer for that to come up. If it comes up, I'll upvote it. I hope that the information I put together will be helpful for people having to deal with quantiles in a mixed R and SAS context. If you disagree and think that my answer is not usefull at all, you can state so by downvoting it which may or may not nudge me to retract it. Or you accept a link to a pdf in which advantages and disadvantages of some median-finding algorithms are listed as a partial answer whilst waiting for a better one. $\endgroup$
    – Bernhard
    Commented Jan 12, 2021 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.