My next suggestion (I don't have time to work out or test an example at the moment): library(abind) tmparr <- abind(m1,m2,m3,...,along=3) OR tmparr <- do.call(c(matlist,list(along=3))) apply(tmparr,c(1,2),mean,na.rm=TRUE) or something along those lines.