Mais algumas opções.

dados <- expand.grid(A=factor(paste("A", 1:3, sep="")),
                     B=factor(paste("B", 1:4, sep="")),
                     rep=1:4)
dados$y <- rgamma(n=48, shape=2, scale=2)

m0 <- by(data=dados, INDICES=dados$B,
         FUN=lm, formula=y~rep)

lapply(m0, summary)
lapply(m0, anova)
sapply(m0, coef)

require(plyr)

m0 <- dlply(.data=dados, "B", .fun=lm, formula=y~rep)

lapply(m0, summary)
lapply(m0, anova)
sapply(m0, coef)

À disposição.
Walmes.