One small problem with apply is that it is not idempotent - it always puts the result vector in the first dimension of the result array. For example: a <- array(1:24, c(2,3,4)) dim(apply(a, c(1,2), I)) dim(apply(a, c(1,2), diff)) I've written a wrapper around apply that rearranges the dimensions into the original order. # Idempotent apply iapply <- function (X, MARGIN, FUN, ...) { dims <- c((1:length(dim(X)))[!(1:length(dim(X)) %in% MARGIN)], MARGIN) res <- apply(X, MARGIN, FUN) if (length(dim(res)) == length(dims)) { aperm(res, order(dims)) } else { res } } dim(iapply(a, c(1,2), I)) dim(iapply(a, c(1,2), diff)) Hopefully this may be of use to someone else. The function isn't completely satisfying as apply doesn't deal nicely with functions that return higher-dimensional arrays (eg. apply(a, 1, I)) and reduces 1-D arrays to vectors (eg. apply(a,3,sum)). Any comments on the code would be greatfully appreciated. Hadley