Thanks to Martin, Marc, Dennis, Rashid and Bill for comment. The short answer is there is no factorial function within the base package, so use the result that gamm(x + 1) = x!. This what the factorial() function in the package gregmisc does: factorial <- function (x) gamma(1 + x) Bill V. also writes: "Your best move is usually to use the lgamma() function (for log gamma) and do any cancellation in the additive log scale first. e.g. for binomial coefficients you would use something like nCr <- function(n, r) round(lgamma(n+1) - lgamma(r+1) - lgamma(n-r+1)) This is vectorized, of course, so works for vectors of values of n and r." Martin M. notes " ... something which we have thought to be too trivial to provide an extra function name for -- particularly in light of the use of the term "factorial" in something like "factorial design" in statistics. If this question continues to pop up with such regularity, we probably should reconsider..." On that note, I think it would be nice when teaching students [with limited stats background] to be able to avoid introducing the gamma function to heighten their confusion! cheers Peter ********************************************************************* Dr Peter Caley CSIRO Entomology GPO Box 1700, Canberra, ACT 2601 Email: peter.caley@csiro.au Ph: +61 (0)2 6246 4076 Fax: +61 (0)2 6246 4000 ********************************************************************* [[alternate HTML version deleted]]