[R-sig-ME] Power Analysis for Linear Mixed Model with Covariates

[Ben Bolker]

Gerasimos Fergadiotis <gf3 at ...> writes:

>  Thank you all for your help so far.  Having as a starting point Ben
> Bolker'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 ?? to define the variance
> covariance matrix for the random effects.

theta is the "lower cholesky factor" of the (SCALED) variance-covariance
matrix (i.e. the random-effects vcov matrix divided by the (scalar) residual
variance): that is, if you have (t1,t2,t3) for the theta-vector of
a 2x2 var-cov matrix, the var-cov matrix will be

  ( t1  0  )  (  t1  t2  )   =    (  t1^2  t1*t2 )
  ( t2 t3  )  (   0  t3  )        (  t1*t2 t3^2  )

the results are similar (but much more complex ...) for the 4-by-4
var-cov matrix you're going to use here.   If you're feeling lazy
you can feed this to Wolfram Alpha:

{{t1, 0, 0, 0 }, {t2, t3, 0, 0}, {t4, t5, t6, 0},{t7,t8,t9,t10}}  
 {{t1, t2, t4, t7 }, {0, t3, t5, t8}, {0, 0, t6, t9},{0,0,0,t10}} 

You can also specify the variance-covariance matrix you want
and use cholesky decomp, e.g.

m <- matrix(1,nrow=4,ncol=4)
diag(m) <- 3
lc <- t(chol(m))
lc[lower.tri(lc,diag=TRUE)]  ## desired theta vector
lc %*% t(lc)

Note that 

> A second question is how I should establish
> the magnitude of sigma as it appears to play
> a significant role in > the simulation.

  You just have to figure out what's a sensible value,
or a range of sensible values, for your system.