[R-sig-ME] problem with estimation of individual and interaction term with small residual variance

[Florian Klinglmueller]

dear list,

i have a problem with estimating the variance components of the 
following model:

i have 10 animals, from each of which 4 different tissue preparations 
are measured in replicates of 5. i'm interested in testing if there are 
any differences between the 4 preparations.
so i have a fixed effect for the preparation, then a random effect for 
the individual and i'd like to include a random interaction term since 
there might be slight differences in how the preparations are done from 
each animal. this leads me to the model:

y(ijk) = 1 + x(i) + z(j) + z(ij) + e(ijk)

where y(ijk) is the measurement from each of the samples, x(i) is the 
fixed effect for the 4 preparations, z(j) is the individual random 
effect, z(ij) is the interaction random effect of animal and preparation 
and finally e(ijk) the residual.

i want to estimate the variances of z(j), z(ij) and e(ijk).

i did a little simulation study (see code below) over several scenarios 
for the different variance components. if i set the residual standard 
error to a very small value (~0.01) keeping individual and interaction 
variance at one, i observe is that the individual effect is 
underestimated as being close to zero (mean sd for z(ij) estimate over 
1000 simulations: .01), although it should be 1. for se=0 i get false 
convergence, however not for small values. consequently the estimates 
are negatively biased, with the amount of bias dependent on the amount 
of residual standard deviation. This also doesn't get better if i 
increase the sample sizes.

so my question is: is this a numerical problem or do i have misspecified 
something?

thanks and regards
florian
 
--- Simulation Code -----

sim = function(rsd){
    sb = 1 # sd of individual random effect
    sg = 1 # sd of interaction random effect
    se = rsd  # residual standard deviation
    a = 1:4 # preparation fixed effecs
    ia = factor(rep(1:4,each=10*5)) # dummy variable for preparation
    ib = factor(rep(1:10,each=5,4))  # dummy variable for individual
    y = a[ia]+rnorm(10,0,sb)[ib] +rnorm(4*10,0,sg)[ia:ib]+rnorm(4*10*5,0,se)

    lmer(y~1+ia+(1|ib)+(1|ia:ib))
}

## 1000 runs of the simulation
foo <- replicate(1000,attr(VarCorr(sim(2))[[2]],'stddev'))