[R] CI on median

[Simon Fear]

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}}