# How to simulate an unreplicated factorial design?

I am dealing with a unreplicated factorial design. I have some illustrative examples but I need to simulate some unreplicated factorial designs. I do not how and what to use. Can $R$ handle this?

For example, I would like to analyse a $2^{4}$ factorial design (factors are A, B, C and D) with only one run and 15 contrasts. I have a single column for response. I would like to compare some methods in the literature to see which method detects active effects better. Thus, I set the active effects to have the same magnitude of $1.5\sigma$ and I would like to generate $100$ response vectors using errors that are i.i.d. with $\mathcal N(0 ,1)$. My true model has four active effects and I would like to simulate $100$ response vectors using this true model $y=3+1.5A+1.5B+1.5C+1.5BC$. But I do not know how to generate data like this using R.

Thanks gung for your reply. I just wrote a simple code before I saw your answer here. I think, I need to build up a bit more R knowledge. Anyway, here it is:

For the analysis of unreplicated factorial designs with $k$ factors and $p=2^{k}-1$ factorial effects (the main effects and interactions), the following model is generally used

$$y=\sum\limits_{i=0}^{p}x_{i}\beta_{i}+\varepsilon_{i}$$

So, Firstly I introduced my sign table for $2^{4}$ and $\beta$ coefficients of so-called active effects.

Sign table consists of rows (runs) and columns (contrasts with general mean).

And then, I created my regression equation with magnitudes of active effects and zeros of remaining inactive effects. My simulated model, for example, was $y=3+1.5A+1.5B+1.5C+1.5BC$.

And then, I run the code below

x=read.csv("sign2.txt", header=TRUE)
sign= as.matrix(x)
is.matrix(sign)

beta= as.matrix(y)
is.matrix(beta)

signt=t(sign)

bs=t(beta %*% signt)

epsilon=matrix( rnorm(16*1,mean=0,sd=1), 16, 1)

response=bs+epsilon


However, unfortunately, it's for one simulation. I will put a loop command to run the simulation n-times.

• Note that there are some differences between your approach & what you had posted originally / I responded to. Eg, the slopes were $1.5$ in the question, whereas this uses $2$. Also, I used reference level coding for the categorical variables, but you are using effect coding. – gung - Reinstate Monica Aug 20 '13 at 20:25
• I just changed the magnitudes in the regression equation. It is my mistake. Sorry. By the way, how you coded factors is not correct, I believe. You need to use (-1, 1) not (0, 1). Because in order to get predicted values or response values for each run (in this case, it's 16) you need to multiply every row in the sign table with the corresponding β coefficients of active effects in the regression equation and then, to add error that is normally distributed. – ARAT Aug 20 '13 at 21:18

The key to generating random data like yours in R is to use ?rnorm. You may also want to set the random seed to a fixed value, via ?set.seed, so that your simulation can be exactly replicated in the future. For convenience, you may prefer to use ?expand.grid to create your factor combinations, although you can do it manually as well (and could conceivably prefer that for the clarity it affords). Once you have set up your simulation, you will want to run it many times to get a sense of its long run behavior; this can be done using ?replicate, I believe, but I usually just nest it in a for loop because I specialize in writing comically inefficient code.

Here is an example:

> set.seed(1)
> A <- c(0,1)
> B <- c(0,1)
> C <- c(0,1)
> D <- c(0,1)
> Xmat <- expand.grid(A=A, B=B, C=C, D=D)
> Xmat
A B C D
1  0 0 0 0
2  1 0 0 0
3  0 1 0 0
4  1 1 0 0
5  0 0 1 0
6  1 0 1 0
7  0 1 1 0
8  1 1 1 0
9  0 0 0 1
10 1 0 0 1
11 0 1 0 1
12 1 1 0 1
13 0 0 1 1
14 1 0 1 1
15 0 1 1 1
16 1 1 1 1
> y <- 3 + 1.5*Xmat$A + 1.5*Xmat$B + 1.5*Xmat$C + 1.5*Xmat$B*Xmat\$C +
+      rnorm(n=16, mean=0, sd=1)
> y
[1] 2.373546 4.683643 3.664371 7.595281 4.829508 5.179532 7.987429 9.738325
[9] 3.575781 4.194612 6.011781 6.389843 3.878759 3.785300 8.624931 8.955066