Thank you all for your help so far.
Having as a starting point Ben Colker's Code (http://rpubs.com/bbolker/11703), I have put together the following code that will allow me to generate data for a LMM with two fixed effects (Treatment and Time), their interaction, and and random intercepts and slopes for each participant.

expdat <- expand.grid(kid = factor(1:500), Time = factor(1:4), Treat = c("XTx", "BAU"))
expdat$obs <- factor(seq(nrow(expdat)))
set.seed(101)
nsim <- 20
beta <- c(100, -7, 8, 15, 20)
theta <- c(15, 15, 15, 15, 1, 0, 0, 0, 0, 0)
ss <- simulate(~Treat + Time + (Time | kid), nsim = nsim, family = gaussian, weights = rep(25, nrow(expdat)), newdata = expdat, newparams = list(theta = theta, beta = beta, sigma = 1))
expdat$Outcome <- ss[, 1]
fit1 <- lmer(Outcome ~Treat + Time + Treat*Time + (Time | kid), data = expdat)

I am puzzled now with who to specify "theta" to give specify correlational matrices for the time points. How is theta transformed to the correlational table that I see in the model output? I am assuming that theta corresponds to ΘD to define the variance covariance matrix for the random effects. 
A second question is how I should establish the magnitude of sigma as it appears to play a significant role in the simulation.

Thank you,
Gerasimos 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Gerasimos Fergadiotis, PhD., CCC-SLP
Assistant Professor
Speech & Hearing Sciences
724 SW Harrison St., RM 84B
Portland State University
Portland, Oregon 97201
Phone: 503.725.2217
Fax: 503.725.9171
Web page: http://aaldresearch.wix.com/aald
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



	[[alternative HTML version deleted]]

