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)
Thanks in advance
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)
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.
