Simulating linear mixed models - the Venables approach

In case you missed it on the R-help list, I urge readers of this list
to consider the understated elegance of the code Bill Venables posted
for simulating data from a simple random effects model.
set.seed(7658943)

fph <- 0.4
Sigh <- sqrt(0.0002)
Sigi <- sqrt(0.04)

reH <- rnorm(90, fph, Sigh)  ## hospid effects
dta <- within(expand.grid(hospid = 1:90, empid = 1:80),
fpi1 <- reH[hospid] + rnorm(7200, fph, Sigi))

One is reminded of John Keats

 'Beauty is truth, truth beauty,?that is all	
 Ye know on earth, and all ye need to know.'
Hi Doug,
I couldn't find this email on the r-help list. Would you mind  
elaborating **briefly** on what is elegant about this?
Eager to learn,
Hank

In case you missed it on the R-help list, I urge readers of this list
to consider the understated elegance of the code Bill Venables posted
for simulating data from a simple random effects model.

set.seed(7658943)

fph <- 0.4
Sigh <- sqrt(0.0002)
Sigi <- sqrt(0.04)

reH <- rnorm(90, fph, Sigh)  ## hospid effects
dta <- within(expand.grid(hospid = 1:90, empid = 1:80),
          fpi1 <- reH[hospid] + rnorm(7200, fph, Sigi))

One is reminded of John Keats

'Beauty is truth, truth beauty,?that is all
Ye know on earth, and all ye need to know.'

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Dr. Hank Stevens, Associate Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/~stevenmh/
http://www.cas.muohio.edu/ecology
http://www.muohio.edu/botany/

"If the stars should appear one night in a thousand years, how would men
believe and adore." -Ralph Waldo Emerson, writer and philosopher  
(1803-1882)
Hi Doug,
 I couldn't find this email on the r-help list. Would you mind elaborating
**briefly** on what is elegant about this?
 Eager to learn,
 Hank
 On Apr 7, 2008, at 12:29 PM, Douglas Bates wrote:
In case you missed it on the R-help list, I urge readers of this list
to consider the understated elegance of the code Bill Venables posted
for simulating data from a simple random effects model.
set.seed(7658943)

fph <- 0.4
Sigh <- sqrt(0.0002)
Sigi <- sqrt(0.04)

reH <- rnorm(90, fph, Sigh)  ## hospid effects
dta <- within(expand.grid(hospid = 1:90, empid = 1:80),
         fpi1 <- reH[hospid] + rnorm(7200, fph, Sigi))
I find the use of within, expand.grid and indexing on the random
effects to be elegant.  Actually, on reexamining the code I think that
Bill should change the first call to rnorm to have mean = 0, not mean
= fph.

The code for the simulation, data display and fit could be collapsed to

library(lme4)
set.seed(7658943)
reH <- rnorm(90, sd = sqrt(0.0002))
dta <- within(data.frame(hospid = gl(90,80)),
              fpi1 <- reH[hospid] + rnorm(length(hospid), 0.4, sqrt(0.04)))
dotplot(reorder(hospid, fpi1) ~ fpi1, dta)
(fm1 <- lmer(fpi1 ~ 1|hospid, dta, verb = TRUE))