[R-sig-ME] Simulating linear mixed models - the Venables approach
Douglas Bates
bates at stat.wisc.edu
Wed Apr 9 19:17:28 CEST 2008
On 4/9/08, Hank Stevens <HStevens at muohio.edu> wrote:
> 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))
More information about the R-sig-mixed-models
mailing list