Standard error of the median

Question

Is the following formula right if I want to measure the standard error of the median in case of a small sample with non normal distribution (I'm using python)?

 sigma=np.std(data)
 n=len(data)
 sigma_median=1.253*sigma/np.sqrt(n)

Answer 1

Based on some of @mary's comments I think the following is appropriate. She seems to be selecting the median because the sample is small.

If you were selecting median because it's a small sample that's not a good justification. You select median because the median is an important value. It says something different from the mean. You might also select it for some statistical calculations because it's robust against certain problems like outliers or skew. However, small sample size isn't one of those problems it's robust against. For example, when sample size gets smaller it's actually much more sensitive to skew than the mean.

Answer 2

The magic number 1.253 comes from the asymptotic variance formula: $$ {\rm As. Var.}[\hat m] = \frac1{4f(m)^2 n} $$ where $m$ is the true median, and $f(m)$ is the true density at that point. The magic number 1.253 is $\sqrt{\pi/2}$ from the normal distribution so... you still are assuming normality with that.

For any distribution other than the normal (and mary admits that this is doubtful in her data), you would have a different factor. If you had a Laplace/double exponential distribution, the density at the median is $1/2b$ and the variance is $2b^2$, so the factor should be $1/\sqrt{2} = 1.414$ -- the median is the maximum likelihood estimate of the shift parameter, and is more efficient than the mean. So you can start picking your magic numbers in different ways...

Getting the median estimate $\hat m$ is not such a big deal, although you can start agonizing about the middle values for the even number of observations vs. inverting the cdf or something like that. More importantly, the relevant density value can be estimated by kernel density estimators, if needed. Overall, this of course is relatively dubious as three approximations are being taken:

1. That the asymptotic formula for variance works for the small sample;
2. That the estimated median is close enough to the true median;
3. That the kernel density estimator gives an accurate value.

The lower the sample size, the more dubious it gets.

Answer 3

Sokal and Rohlf give this formula in their book Biometry (page 139). Under "Comments on applicability" they write: Large samples from normal populations. Thus, I am afraid that the answer to your question is no. See also here.

One way to obtain the standard error and confidence intervals for the median in small samples with non-normal distributions would be bootstrapping. This post provides links to Python packages for bootstrapping.

Warning

@whuber pointed out that bootstrapping the median in small samples isn't very informative as the justifications of the bootstrap are asymptotic (see comments below).

Answer 4

Not a solution here, but perhaps helpful:

Suppose your data distribution is $p(x)$, and let $P(x) = \int_{-\infty}^x p$ be the cumulative density function. So the median of the distribution is the number m such that P(m) = 1/2.

Following this helpful page we can compute the distribution of a number $x$ being the median of $n$ samples of this distribution. I think it is $q(x) = c_n p(x) (P(x)(1-P(x)))^{(n-1)/2}$. Here $c_n$ is the appropriate constant to make this a probability distribution, and I think it is n-1 choose (n-1)/2 if n is odd (unsure on that).

Finally, you would like to know the variance of q(x), which you may be able to reason about with this formula.

Answer 5

There is an empirical procedure for obtaining a confidence interval for the sample median. The procedure is non-parametric and relies on the binomial distribution. It can be found in Ott and Longnecker, 2015 in the section named ‘Inferences about the median’. Stata implements the procedure as the ‘centiles’ command and the Stata doc provides a mathematical justification with references.

Here is the procedure for a 95% CI and a python script. The standard error of the median is determined from the CI. The results are the same as the results from the Stata 'centiles' command.

import numpy as np
from scipy.stats import binom
# n = 25
data = [1.1, 1.2, 2.1, 2.6, 2.7, 2.9, 3.6, 3.9, 4.2, 4.3, 4.5, 4.7, 5.3,
        5.6, 5.8, 6.5, 6.7, 6.7, 7.8, 7.8, 14.2, 25.9, 29.5, 34.8, 43.8]
median = np.median(data)
print(f'median: {median}')
# the distribution
n = 25
p = .5
rv = binom(n, p)
# 95% critical value
q = .05
binom_critical = rv.ppf(q=q)
print(f'binom 95% critical value: {binom_critical}')
# the 95% CI for the median
L_q = int(binom_critical)
U_q = int(n - binom_critical)
print(f'L_q: {L_q} U_q: {U_q}')
lower_ci = data[L_q - 1]
upper_ci = data[U_q - 1]
print(f'lower_ci: {lower_ci} upper_ci: {upper_ci}')
median: 5.3
binom 95% critical value: 8.0
L_q: 8 U_q: 17
lower_ci: 3.9 upper_ci: 6.7

The standard error is CI/2. In this case 6.7 – 3.9 / 2 = 1.4. This is analogous to the normal case, where if given the CI for the mean, the standard error is calculated as:

For the non-parametric method for the median, there is no t-statistic because the confidence level is embodied in the critical values for the binomial distribution.