two.way = function(n = 8, mu = 5, alpha = c(2, 3), beta = c(.5, .7),
               ab = c(.2, .3, .4, .5) ){

m1 = mu + alpha[1] + beta[1] + ab[1]
m2 = mu + alpha[1] + beta[2] + ab[2]
m3 = mu + alpha[2] + beta[1] + ab[3]
m4 = mu + alpha[2] + beta[2] + ab[4]

 y = c(rnorm(n, m1),
       rnorm(n, m2),
       rnorm(n, m3),
       rnorm(n, m4))

 A = rep(1:2, each = 2*n)
 B = rep(rep(1:2, each = n), 2)

anova(aov(y ~ A * B))
}
# Example of use:
two.way()