Is R's trimmed means function biased?

I've been trying to understand how R's trimmed mean function works. I suspect it might be biased, but would like to get feedback here before I file a bug-report (if this is an inappropriate forum for such a question please let me know).

Consider the following example (sorted for clarity):

> x <- c(1, 2, 3, 45, 55, 56, 57, 58, 403, 900)
> length(x)
[1] 10
> percent.to.trim <- 0.25


One would expect that the trimmed mean would be based on length(x) - 2 * length(x) * percent.to.trim = 5 observations. If this is the case, there are three possible methods of calculation:

1. To remove the two smallest observations and the three largest observations. In this case, mean(x[3:7]) = 43.2.
2. To remove the three smallest observations and the two largest observations. In this case, mean(x[4:8]) = 54.2.
3. To remove the two smallest observations, the two largest observations and to average the third smallest and third largest observation (this is to be preferred for obvious reasons). In this case mean(c(x[4:7], mean(c(x[3], x[8])))) = 48.7.

However, R gives the following result:

> mean(x, trim=percent.to.trim)
[1] 45.66667


This must result from the trimmed mean being based on length(x) - 2 * floor(length(x) * percent.to.trim) = 6 observations, thus removing only the first two and last two observations:

> mean(x[3:8])
[1] 45.66667


Is this a bug?

I'm using R version 2.13.1 (2011-07-08), svn rev 56322.

• It's not a bug. If you read the help file, you'll see that an equal number of observations are trimmed from each end. Oct 5, 2011 at 7:44
• @mark999 Are you getting that from the term symmetrically trimmed mean? Or is there more explicit wording that I am somehow missing? Method (3) in my question could also be considered symmetrical IMHO. Oct 5, 2011 at 7:49
• Yes, there's that, but also the fraction (0 to 0.5) of observations to be trimmed from each end of x before the mean is computed. Oct 5, 2011 at 7:57
• @mark999 Thanks. The trimmed from each end is telling, although I guess it is still ambiguous as rounding behaviour is not specified. As the answers below point out, that is what source-code is for. Oct 5, 2011 at 8:08
• The thread at stats.stackexchange.com/questions/4252/… discusses why neither of the first two options is appropriate. (The third one is an intriguing generalization.) In short, trimming is done by equal amounts at both ends to minimize the bias of the estimator.
– whuber
Oct 5, 2011 at 15:46

In such cases it is a good idea to look at the source code. In this case the source code of mean.default

> mean.default
function (x, trim = 0, na.rm = FALSE, ...)
{
if (!is.numeric(x) && !is.complex(x) && !is.logical(x)) {
warning("argument is not numeric or logical: returning NA")
return(NA_real_)
}
if (na.rm)
x <- x[!is.na(x)]
if (!is.numeric(trim) || length(trim) != 1L)
stop("'trim' must be numeric of length one")
n <- length(x)
if (trim > 0 && n) {
if (is.complex(x))
stop("trimmed means are not defined for complex data")
if (any(is.na(x)))
return(NA_real_)
if (trim >= 0.5)
return(stats::median(x, na.rm = FALSE))
lo <- floor(n * trim) + 1
hi <- n + 1 - lo
x <- sort.int(x, partial = unique(c(lo, hi)))[lo:hi]
}
.Internal(mean(x))
}


From this the trimming is based on the lower and upper index bounds lo <- floor(n * trim) + 1 and hi <- n + 1 - lo, which are the design decisions implemented in R. Now, while writing this I see that Nick Sabbe posted more or less the same answer, but I leave the entire code posted here as well.

P.S. I think the clue to reading the help page is that it says a symmetrically trimmed mean, not the symmetrically trimmed mean.

Apart from how it should work, here's how it does work:

You can simply get to the implementation of the relevant part of mean by typing it at the prompt. Then you'll notice that it is a generic, so for your case, mean.default will be called.

Next you call up its implementation by once again typing mean.default at the R prompt. There, you can see that the relevant part of the implementation reads:

lo <- floor(n * trim) + 1
hi <- n + 1 - lo
x <- sort.int(x, partial = unique(c(lo, hi)))[lo:hi]


So, for 10 items and trim=0.25, the first 2 and the last 2 are dropped.

HTH.

• Thanks Nick, I still haven't gotten my head around how S works enough to have figured out the mean.default trick and was stuck on the mean generic error message. Oct 5, 2011 at 8:05