JCF,

Você pode usar o que já está disponível na agricolae, no caso a função design.lsd() para gerar o delineamento. São funções eu recomendo os usuários aplicarem para fazer o sorteio dos níveis dos fatores às unidades experimentais. A partir dele você obtém a matriz do delineamento e daí o resto todo,

require(agricolae)
help(design.lsd, help_type="html")

varieties <- c("perricholi","yungay","maria bonita","tomasa")
lsd <- design.lsd(varieties, seed=23)
lsd
str(lsd)

X <- model.matrix(~row+col+varieties, data=lsd)
dim(X)
colnames(X)
betas <- runif(ncol(X)) # use os valores paramétricos aqui
lsd$Ey <- X%*%betas     # valor esperado de y
lsd$y <- rnorm(nrow(X), mean=lsd$Ey, sd=1) # realização de y
plot(y~varieties, lsd)

À disposição.
Walmes.

==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================