R help - Sep 2003 - CI on median

On Sat, 27 Sep 2003 13:25:57 +0200
Troels Ring <tring at gvdnet.dk> wrote:
> 
> Dear friends, I'm probably wrong but is there anything better than 
> bootstrap to get a confidence interval of the median from a population with
> unspecified distribution ?
> Best wishes
> Troels Ring, Aalborg, Denmark
The bootstrap doesn't work perfectly for the median.  Standard nonparametric
statistics texts provide a distribution-free CI for the median, like:

medianCI <- function(y) {
    y <- sort(y[!is.na(y)])
    n <- length(y)
    if(n < 4) return(c(NA,NA))
    r <- pmin(qbinom(c(.025,.975), n, .5) + 1, n)  ## Exact 0.95 C.L.
    c(y[r[1]],y[r[2]])
}

In an upcoming release of the Hmisc package you can also get the Harrell-Davis
distribution-free estimator of the median which is slightly more efficient than
the traditional estimator, for small samples.  A new hdquantile function
provides an estimate of the standard error for this estimator.

---
Frank E Harrell Jr    Professor and Chair            School of Medicine
                      Department of Biostatistics    Vanderbilt University

Dear friends, I'm probably wrong but is there anything better than 
bootstrap to get a confidence interval of the median from a population with 
unspecified distribution ?
Best wishes
Troels Ring, Aalborg, Denmark

You can also get CIs for medians as a side effect of boxplot.

If for some reason you need to get the same result as you
would get from SAS, you could use the following (non-optimised
but hopefully intelligible) instead.

HTH

###############################################
# 95% confidence interval function for the median
# prints out CI, returns invisible list of median.hat, lowlim, highlim
# matches the output of SAS proc FREQ CIPCTLDF option
###############################################
#
"ci.median.sas" <- function(x, do.print=TRUE)
{
x <- x[!is.na(x)]
median.hat <- median(x)
nobs <- length(x)
x.ranked <- sort(x)
if (nobs < 2 ) {
  if (do.print) print(paste(median.hat, " (---, ---)", sep =
""))
  invisible(return(median.hat, lowlim=NA, highlim=NA))
}
# in general, rank of low cutoff point is the position of q in 0:nobs,
so add
one
# to give the position in 1:nobs
lowrank <- qbinom(0.025, nobs, 0.5) + 1
# but make sure that this choice does not use up too much probability,
making
a
# 95% interval impossible (except for very small samples, low will
always be
one):
plow <- pbinom(lowrank-1, nobs, 0.5) 
if (plow > 0.05 && lowrank > 1) {
  lowrank <- lowrank - 1
  plow <- pbinom(lowrank-1, nobs, 0.5)
}
# then find the corresponding high cutoff giving (better than) 95%
coverage
overall:
# nb for very small samples plow will still exceed 0.05, so we use
min(1,p)
to
# prevent potential eror
highrank <- qbinom(min(1, 0.95+plow), nobs, 0.5) + 1
# again +1 for index in 1:nobs rather than 0:nobs
lowlim <- x.ranked[lowrank]
highlim <- x.ranked[min(highrank, nobs)]
if (do.print) print(paste(median.hat, " (", lowlim, ", ",
highlim, ")",
sep ""))
invisible(return(median.hat, lowlim, highlim))
}
 

Simon Fear
Senior Statistician
Syne qua non Ltd
Tel: +44 (0) 1379 644449
Fax: +44 (0) 1379 644445
email: Simon.Fear at synequanon.com
web: http://www.synequanon.com
 
Number of attachments included with this message: 0
 
This message (and any associated files) is confidential and\...{{dropped}}