[R-sig-ME] Zero estimated random effects
Sarah Barry
sarah at stats.gla.ac.uk
Tue Aug 14 11:40:26 CEST 2007
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
More information about the R-sig-mixed-models
mailing list