The apply function also works with multi-dimensional arrays, I think this is what you want to achieve using a 3d array: aaa <- array(NA, dim = c(2, dim(a))) aaa[1,,] <- a aaa[2,,] <- a2 apply(aaa, 3, function(x)lm(x[1,]~x[2,]))