Solved – How to simulate an unreplicated factorial design

experiment-designrrandom-generationregressionsimulation

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.

An example from Montgomery for $2^4$ unreplicated factorial design

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

\begin{equation}
y=\sum\limits_{i=0}^{p}x_{i}\beta_{i}+\varepsilon_{i}
\end{equation}

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).
enter image description here

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$.
enter image description here

And then, I run the code below

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

y=read.csv("beta2.txt", header=TRUE)
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.

Best Answer

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