Hello,

I simulate a multilevel model, with 2 levels,  of the form:

y_ij=beta_0+beta_1*x_ij+alfa_i+eta_ij

Individual-level units: i=1,...,150
Group-level units: j=1,...,10
The error therms are assumed to follow N(0,1)


I want to compute the magnitude of the bias of the estimators when there is
correlation between x_ij and alfa_i and between x_ij and eta_ij.

I assign the following true values for the beta parameters: beta_0=10,
beta_1=2. I run 250 replications.

I study 2 cases, depending on the magnitude of the correlation:

i)  corr(x_ij, alfa_i) = 0 , corr(x_ij, eta_ij) = 0.3
ii)  corr(x_ij, alfa_i) = 0.3 , corr(x_ij, eta_ij) = 0.3


What I am doing is actually a replications of a paper by Ebbes et al.
(2004). The problem encountered is that I do not get the same magnitude of
the bias as the authors. I get considerably smaller bias.

Below are the results of the authors across 250 replications (means and
std. deviations):


                                                     Case


i                       ii

Fixed Effects          beta_0             -                       -
                               beta_1       2.43 (0.04)          2.42 (0.04)


Random Effects     beta_0        7.88 (0.20)          5.79 (0.30)
                               beta_1       2.42 (0.04)          2.84 (0.06)
                              sigma_alfa  0.99 (0.13)          0.00 (0.00)
                              sigma_eta   0.90 (0.03)         0.91 (0.04)


My code:


require(lme4)
require(MASS)
k <- 10            # total number of firms
n <- 15            # number of employees within each firm
j <- k*n        # total number of employees
beta0=10
beta1=2
rho1=0.0001     # this is Case iii
rho2=0.3
beta_zero <- rep(NA,250)
beta_one <- rep(NA, 250)
for (i in 1:250) {
c <- cbind(FID = sort(rep(1:k,n)), Fint =rnorm(j,0,1),Z = rnorm(j,0,1))
            #generate level 2 data
alfa_i <- rnorm(j, mean=0, 1)
eta_ij <- rnorm(j, mean=0, 1)
b <- cbind(EID = 1:j, Eint=rnorm(j,0,1), X = rnorm(j,mean =0, sd
=1)+eta_ij+alfa_i)  # generate level 1 data
m1dat <- data.frame(c,b)
m1dat <- within(m1dat, {
   EID <- factor(EID)
   FID <- factor(FID)}
   )
outcome <- data.frame(y_ij =
beta0+beta1*(m1dat[,'X'])+rho1*alfa_i+rho2*eta_ij )
m1dat <- cbind(m1dat,outcome)
rand_eff <- lmer(y_ij ~1+X+(1|FID), data=m1dat)
beta_zero[i] <- fixed.effects(rand_eff)[[1]]
beta_one[i] <- fixed.effects(rand_eff)[[2]]
}
 means <- c(mean(beta_zero),mean(beta_one))



Thanks in advance
-- 


*Raluca Ioana Gui*


<http://www.socialnetworks.uzh.ch/>

	[[alternative HTML version deleted]]

