Se é que eu entendi sua dúvida... vars=list() for(i in 1:30) vars[[i]] <- rnorm(48) names(vars) <- paste('var',1:30,sep='') vars varis <- do.call('cbind',vars) da <- expand.grid(FatA=rep(LETTERS[1:6],4), FatB=letters[1:2]) dad <- data.frame(da,varis) tabs <- apply(dad[,-c(1:2)], 2, function(x) aov(x ~ FatA*FatB,dad)) lapply(tabs,summary)