Dear guRus, is there a version of boot() that deals with array data and array-valued statistics? For example: foo <- array(rnorm(120),dim=c(3,5,8)) means <- apply(foo, MARGIN=c(2,3), FUN=mean) means contains the means over the first dimension of foo: a 5x8 array. Now I would like to bootstrap this array, perhaps with something along the lines of some array.boot(), with an additional MARGIN parameter: means.boot <- array.boot(foo, statistic=mean, MARGIN=c(2,3), R=1000) I would like means.boot to contain (e.g., in analogy to boot(), in a component $t) an array of dimensions 1000x5x8, containing componentwise means of sampled slices of foo, e.g., to get confidence intervals like this: apply(means.boot$t,MARGIN=c(2,3),FUN=quantile,probs=c(.975,.025)) I have been playing around with plyr and various flavors of apply, and searching only yielded lots of hits for boot.array(), which is something completely different... Any help would be greatly appreciated! Cheers Stephan