S 中的统计模型，由 Chambers 和 Hastie（Chapmann 和 Hall，1991 年；或所谓的白皮书），以及在较小程度上与 S 的现代应用统计，由 Venables 和 Ripley （Springer，2002 年，第 4 版），包括一些关于 DoE 的材料以及 S 和 R 中常见设计的分析。 Vikneswaran为“实验设计”写了一个 R 伴侣，虽然它不是很完整（恕我直言），但在 CRAN的投稿部分还有很多其他教科书这可能会帮助您入门。

除了教科书外，CRAN Task View on Experiment of Experiments (DoE) & Analysis of Experimental Data还有一些很好的软件包，可以简化各种实验设计的创建和分析；我可以想到dae、agricolae或AlgDesign（带有漂亮的小插图），仅举几例。

我感觉最近由 Owen Jones、Robert Maillardet 和 Andrew Robinson（2009 年）撰写的“使用 R 的科学编程和科学模拟简介”可能是您正在寻找的内容。

Journal of Statistical Software对此也有非常积极的评价。.

虽然这本书不是专门针对模拟实验数据的，但它可能会让你朝着你想要的方向前进。

这是我为此目的编写的一些代码的示例。实验设计为：有四个氮水平，每个水平六个重复。这些数据可以使用单向方差分析进行测试，但由于水平是连续的，我测试了不同曲线的拟合。

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]")