Reference or book on simulation of experimental design data in R

Question

Can someone tell a reference and/or book that explain how to use R for simulation of experimental design data? Thanks

Answer 1

Statistical models in S, by Chambers and Hastie (Chapmann and Hall, 1991; or the so-called White Book), and to a lesser extent Modern Applied Statistics with S, by Venables and Ripley (Springer, 2002, 4th ed.), include some material about DoE and the analysis of common designs in S and R. Vikneswaran wrote An R companion to "Experimental Design", although it is not very complete (IMHO), but there are a lot other textbooks in the Contributed section on CRAN that might help you get started.

Apart from textbook, the CRAN Task View on Design of Experiments (DoE) & Analysis of Experimental Data has some good packages that ease the creation and analysis of various experimental designs; I can think of dae, agricolae, or AlgDesign (which comes with a nice vignette), to name a few.

Answer 2

I have the feeling that pretty recent "Introduction to Scientific Programming and Scientific Simulation Using R" by Owen Jones, Robert Maillardet, and Andrew Robinson (2009) could be what you are looking for.

There is also a very positive review for it in the Journal of Statistical Software..

Although this book is not specifically targeted at simulating experimental data it will probably get you in the direction you want to go.

Answer 3

Here is an example of some code that I wrote for this purpose. The experimental design is: there are four levels of nitrogen and six replicates at each level. These data could be tested using a one-way ANOVA, but as the levels are continuous, I tested the fit of different curves.

set.seed(1)
library(nlme)
library(ggplot2)
### Below is a set of practice data 

## 1. four levels of Nitrogen:
##    0,1,4,10
N <- c(rep(0,6),rep(1,6),rep(4,6),rep(10,6))

## 2. variance 
s <- 2

## 3. Data simulated to provide examples of the
##    Various hypothesized responses of Y to N 

## 3.1 asymptotic incr Y = 10*N/(1+N) + 10
asym <- c(rnorm(6,10,s),rnorm(6,15,s),rnorm(6,18,s),rnorm(6,19,s))
## 3.2 Y = 0*N + 10 the Null model
m0 <- c(rnorm(24,10,s))
## 3.3 Y =0.2*N+10 a shallow slope
m1  <- c(rnorm(6,10,s),rnorm(6,10.2,s),rnorm(6,10.8,s),rnorm(6,12,s))
## 3.4 Y = 1*N + 10 a more steep slope
m4 <- c(rnorm(6,10,s),rnorm(6,14,s),rnorm(6,26,s),rnorm(6,50,s))
## 3.5 Y = 4*log10(N)+ 10 an log-linear response
lm4 <- c(rnorm(6,10,s),rnorm(6,12.4,s),rnorm(6,15.6,s),rnorm(6,18.2,s))
## 3.6 'Hump' with max at N=1 g m-2 yr 
hump <- c(rnorm(6,10,s),rnorm(6,20,s),rnorm(6,9,s),rnorm(6,8,s))

## A function to compare the fit of five models:

fn.BIC.lmnls <- function (x,y,shape){
  foo.null <- lm ( y~1)
  foo.poly1 <- lm(y~x)
  foo.poly2 <- lm(y ~ x + I(x^2))
  foo.poly3 <- lm(y ~ x + I(x^2) + I(x^3) )
  foo.mm <- nls ( y~ (a*x)/(b + x),start=list(a=1,b=1))
  bic <- BIC(foo.null,foo.poly1,foo.poly2,foo.poly3,foo.mm)
  print(bic)
  return(bic)
}

### now, plot data and print BIC values for each of the data sets
par(mfrow = c(3,2))
fn.BIC.lmnls    (N,m0,"Y = 0*N +10")
fn.BIC.lmnls    (N,m1,"Y = 0.2*N + 10")
fn.BIC.lmnls (N,m4,"Y = 1*N + 10") 
fn.BIC.lmnls (N,lm4,"Y = 0.4*log10(N) + 10")
fn.BIC.lmnls (N,asym,"Y = 10+10*N/(1+N)")#Y = 20*N/(5+N)
fn.BIC.lmnls (N,hump,"Y =  [10,20,9,8]")