[R-sig-ME] Xj correlated with Uj?

[Malcolm Fairbrother]

Dear list,

Code below simulates a simple two-level dataset, but where a level-2 covariate (Xj) is correlated with the level-2 random effects (Uj). In fitting the model, using the call below, the estimate for the fixed effect of Xj is substantially biased.

Is there a way to recover the correct parameter (1)? Presumably it's possible if one knows the variance-covariance matrix "s", but could it be done without that knowledge?

Any thoughts (particularly though not only suggested code) would be much appreciated.

- Malcolm


library(lme4); library(multicore); library(mnormt)
N <- 25
T <- 30
s <- matrix(c(2,0.7,0.7,3), ncol=2)
dgp <- function(N, T, s) {
	dat <- data.frame(group=1:N, rmnorm(n=N, varcov=s))[rep(1:N,each=T),]
	names(dat)[2:3] <- c("Xj", "Uj")
	dat$y <- 1 + dat$Xj + dat$Uj + rnorm(nrow(dat))
	dat
	}
res1 <- do.call("rbind", mclapply(1:100, function(yy) fixef(lmer(y ~ Xj + (1 | group), dgp(N=N, T=T, s=s)))))
colMeans(res1)