[R-sig-ME] Zero estimated random effects

[Sarah Barry]

Hi,

I have longitudinal facial shape data on two treatment groups, which I 
break down into individual coordinates and fit linear mixed effects 
models to each of the pairwise combinations of these coordinates (and 
subsequently aggregate the results).  I would like to fit a random 
intercept with a different variance for each of the two coordinates 
within any model, along with various fixed effects (for coordinate r=1,2 
measured on individual i at time t):

y_{ir}(t) = \beta_{0r} + b_{ir} + \beta_{1r} group_i + \beta_{2r} t + 
\beta_{3r} group_i : t + \epsilon_{ir}(t), 

where \epsilon_{ir}(t) ~ N(0, \sigma^2) and the random effect b_{ir} ~ 
N(0, \sigma_r^2).

The problem is that for a small selection of coordinate pairs, the 
responses for one coordinate are much more variable for the other and 
when I try to fit the model using lmer2, it estimates all the random 
effects for the coordinate with less variation as zero.  This results in 
the variance of that random effect also being estimated as zero.  This 
doesn't seem to happen with lme, even though the variance is still 
estimated as very small.  Using lmer, a warning is given and the 
variance is estimated as very small, although still not exactly zero. 

I don't think it's unreasonable to estimate a different random effects 
variance for each coordinate, especially since the variation is clearly 
so different.  But am I overparameterising, is my syntax wrong or is 
there a problem with lmer2?  I enclose a reproducible example below, 
including code for a plot of the responses.

set.seed(100)
n.subj <- 200
n.times <- 3
n.coords <- 2
simdata <- data.frame(coord1=c(rep(1,n.subj*n.times), 
rep(0,n.subj*n.times)))
simdata$coord2 <- c(rep(0,n.subj*n.times), rep(1,n.subj*n.times))
simdata$coord <- ifelse(simdata$coord1==1, 1, 2)
simdata$ID <- rep(1:n.subj, each=n.times*n.coords)
simdata$time <- rep(1:n.times, n.subj*n.coords)
simdata$group <- rep(c(1,0,1,0), each=n.subj*n.times/2)
simdata$y <- rep(NA, dim(simdata)[1])
for (i in 1:dim(simdata)[1]) simdata$y[i] <- rnorm(1, 
simdata$group[i]*(simdata$coord[i]-1)*40+simdata$time[i], 
1+simdata$coord[i]*10)
plot(simdata$y[simdata$coord==1 & simdata$group==0], 
simdata$y[simdata$coord==2 & simdata$group==0], xlim=range(simdata$y), 
ylim=range(simdata$y), pch=20, xlab="Coordinate 1", ylab="Coordinate 2")
points(simdata$y[simdata$coord==1 & simdata$group==1], 
simdata$y[simdata$coord==2 & simdata$group==1], pch=20, col="red")

require(lme4)
lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+time:group) 
+ (0+coord1|ID)+(0+coord2|ID), data=simdata)
lmer(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+time:group) 
+ (0+coord1|ID)+(0+coord2|ID), data=simdata)
require(nlme)
summary(lme(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+time:group), 
random=list(ID=pdDiag(~-1+coord1+coord2)), data=simdata))

Thanks in advance for any help.

Regards,
Sarah

-- 
Sarah Barry, MSc
PhD student
Department of Statistics
University of Glasgow
Tel: +44 (0)141 330 2474
Fax: +44 (0)141 330 4814
www.stats.gla.ac.uk/~sarah