Dear list,

My understanding is that a 1.0 version of lme4 is coming out soon, but
currently I'm getting some disconcerting results from a simple simulation.
lme4.0, in contrast, generates sensible results. I was investigating the
properties of mixed models fitted to very small numbers of level-2 units,
and got the results below. Maybe there's some reason for this, but I'd be
interested to hear... Or alternatively, I thought people might be
interested to know.

Cheers,
Malcolm



library(lme4) # lme4_0.99999911-1
nsims <- 200 # number of simulations
set.seed(080813)

dgp <- function(N=5,n=6) { within(data.frame(grp = gl(N,n), x1 =
rnorm(N*n), x2 = runif(N)[rep(1:N, each=n)]), y <- 1 + x1 + x2 +
rnorm(N)[grp] + rnorm(N*n)) }

c.mer <- function(mod) { c(fixef(mod), c(unlist(lapply(VarCorr(mod),
diag)), attr(VarCorr(mod), "sc")^2)) } # function to extract FEs, REs

sim <- function(N, n) { c.mer(lmer(y ~ x1 + x2 + (1 | grp), dgp())) }

apply(sapply(1:nsims, sim), 1, function(xx) c(mean(xx),sd(xx)))
# the means, in row 1, should be 1s across the board... but aren't

     (Intercept)          x1       x2 grp.(Intercept)
[1,]    1.418810 -0.02859056 0.143737        1.151423 2.0353251
[2,]    1.474223  0.28381662 2.815086        1.138676 0.6091787

detach(package:lme4)
library(lme4.0) # lme4.0_0.9999-1

set.seed(080813)

apply(sapply(1:nsims, sim), 1, function(xx) c(mean(xx),sd(xx)))
# everything's much closer to 1...

     (Intercept)        x1       x2 grp.(Intercept)
[1,]    1.055962 1.0008875 1.050862        1.124547 1.0029512
[2,]    1.335184 0.2055197 2.288807        1.155740 0.2802394

	[[alternative HTML version deleted]]

