# Simulate data for 2 x 2 ANOVA with interaction

What is the best way to simulate ANOVA data for a 2 x 2 design with interaction using a regression approach? I want to generate the data so that I know the true regression coefficients of the model when running lm() in R.

lm(A * B, data=df)


You can simulate this by running your regression equation backwards into R.

You equation is

$$y = \alpha + \beta_1 A + \beta_2 B + \beta_3 A B + \epsilon$$

where $y$ is the dependant variable, $\alpha$ the intercept, $\beta_n$ your coefficients for your $n$ predictors, in this case, $A$, $B$, and $A*B$ (the interaction), and $\epsilon$ the normally distributed noise, or error.

In R:

# Set coefficients
alpha = 10
beta1 = .3
beta2 = -.5
beta3 = -1.1

# Generate 200 trials
A = c(rep(c(0), 100), rep(c(1), 100)) # '0' 100 times, '1' 100 times
B = rep(c(rep(c(0), 50), rep(c(1), 50)), 2) # '0'x50, '1'x50, '0'x50, '1'x50
e = rnorm(200, 0, sd=1) # Random noise, with standard deviation of 1

# Generate your data using the regression equation
y = alpha + beta1*A + beta2*B + beta3*A*B + e

# Join the variables in a data frame
data = data.frame(cbind(A, B, y))

# Fit an ANOVA
model = aov(y ~ A*B, data=data)
summary(model)


Output:

             Df Sum Sq Mean Sq F value   Pr(>F)
A             1   0.97    0.97    0.79 0.375336
B             1  63.76   63.76   51.97  1.2e-11 ***
A:B           1  15.02   15.02   12.24 0.000577 ***
Residuals   196 240.48    1.23
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Because this is a $2\times2$ ANOVA, we can also fit it as a regression model, revealing that the coefficients estimated from the data are close to those we set when simulating it.

lin.model = lm(y ~ A*B, data = data)
summary(lin.model)

Call:
lm(formula = y ~ A * B, data = data)

Residuals:
Min       1Q   Median       3Q      Max
-2.69066 -0.78089  0.00864  0.72583  2.93863

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   9.9081     0.1566  63.251  < 2e-16 ***
A             0.4090     0.2215   1.846 0.066386 .
B            -0.5811     0.2215  -2.623 0.009404 **
A:B          -1.0963     0.3133  -3.499 0.000577 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.108 on 196 degrees of freedom
Multiple R-squared:  0.249, Adjusted R-squared:  0.2375
F-statistic: 21.67 on 3 and 196 DF,  p-value: 3.67e-12


EDIT

In response to your comment, you could generate more complicated designs by specifying them at minimal length (i.e. what the data would look like with one trial per condition), and then use rep on that.

# 3x2 balanced design
A = c(0,0,0,0,1,1,1,1)
B = c(0,0,1,1,0,0,1,1)
C = c(0,1,0,1,0,1,0,1)

# Repeat 50 times, for N = 200
A = rep(A, 50)
B = rep(B, 50)
C = rep(C, 50)

• thanks, this seems promising. I wonder if the procedure of generating trials could be simplified, e.g. using contrast matrices. If not, more complicated designs are pretty error prone due to all the rep(). Any ideas? – beginneR Sep 18 '14 at 10:42
• Off the top of my head (but mind you, I'm no expert, so there may be a better way of doing this), you could try something like what I've added to the bottom of the answer (it doesn't format nicely in a comment). – Eoin Sep 18 '14 at 10:52
• Eoin, I'm wondering if your method perfectly takes into account all assumptions of a two-way fixed-effects ANOVA? Also, I was wondering how your formulation would change if one of your factors (e.g., $A$) could be a random factor? – rnorouzian Oct 11 '17 at 17:32

In addition to the answers already given, note that one can easily take a fitted model and simulate from that by using the simulate command.

It takes a fitted model object and returns one or more simulated sets of data (specified by nsim). See ?simulate.

> df
y A B
1  5.808711 0 0
2  4.551636 0 1
3  6.347596 0 0
4  4.770539 0 1
5  5.694007 0 0
6  7.110004 1 1
7  8.766641 1 0
8  7.245431 1 1
9  8.022540 1 0
10 7.358717 1 1
> abintfit <- lm(y ~ A * B, data=df)
> simulate(abintfit, nsim=2)
sim_1    sim_2
1  5.967064 6.455663
2  5.120627 4.398695
3  5.801663 5.899567
4  4.912768 4.431236
5  5.788909 6.099048
6  7.607631 7.281159
7  8.377021 8.549396
8  7.473821 7.254118
9  8.580709 8.393214
10 7.391657 6.820862


In a 2 x 2 ANOVA design you have four different mean values. These are your regression coefficients in this case and the means corresponding to the four different combinations of the two factors. This will give an interaction of the two factors as you allow each of the four combinations to vary freely of each other. You can use ${\tt rnorm }$ to simulate using the mean values (regression coefficients) and some fixed variance. Then you'll have the set-up you talk about.