[R-sig-ME] Different random effects variances for outcomes and groups
Sarah Barry
sarah at stats.gla.ac.uk
Wed Sep 19 14:28:03 CEST 2007
Dear all,
I wonder if someone could give me some insight on coding lmer. I have
facial shape data on children in two groups at four time points (3,6,12
and 24 months). Each child has a set of coordinate positions measured
on their face at each time point (the set of coordinates is the same
across individuals and times). Take coordinates 1 and 2 only for now
(reproducible code at the bottom of this email for simulated data).
If I plot the trends for coordinates 1 and 2 for each individual over
time, there is a different amount of variance amongst the individuals
(at least in the intercept, and maybe in the slope) for the two
coordinates and also within the two groups, with group 1 (cleft) having
higher variation than group 0 (control). I want to allow for these
sources of variation in the model. The other thing is that I would
expect coordinate positions within an individual to be correlated so I
also want to allow for this. The model, therefore, would be (for
coordinate r=1,2 measured on individual i at time t, group_i an
indicator variable taking value one for group 1 and zero otherwise):
y_{ir}(t) = \beta_{0r} + \beta_{1r} group_i + \beta_{2r} t + \beta_{3r}
group_i : t + b_{ir0} + group_i * b_{ir1} + \epsilon_{ir}(t),
where \epsilon_{ir}(t) ~ N(0, \sigma2) and the random effect b_{irp} ~
N(0, \sigma_{rp}2), for p=0,1. I think the following code is
appropriate (model 1):
lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+time:group)
+ (0+coord1+coord2+coord1:group+coord2:group|ID), data=simdata)
where coord1 and coord2 are indicator variables for coordinates 1 and 2,
respectively, time is continuous and group is an indicator variable
taking value one for the cleft group and zero for the controls. Does
the random effects part make sense? I'm especially unsure about
allowing correlations between all of the random effects terms, although
I think that's it's appropriate under this parameterisation because each
person has a value for both coordinates 1 and 2, and the group effect is
additional.
An alternative parameterisation is (model 2):
lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+time:group)
+ (0+gp0:coord1+gp0:coord2|ID)+(0+coord1:group+coord2:group|ID),
data=simdata),
where gp0 is an indicator variable taking value one if the individual is
in group 0 and zero otherwise. It seems to me that this should be
equivalent to model 1, but it doesn't appear to be (perhaps this just
comes down to fewer correlations estimated in model 2).
If a correlation between random effects is estimated to be 1 or -1, is
this generally because the model is over-parameterised?
Reproducible code is below.
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(rep(1:n.subj, each=n.times),2)
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$gp0 <- 1-simdata$group
simdata$y <- rep(NA, dim(simdata)[1])
randeff <- c(rep(rnorm(n.subj/2, 0, 15),each=n.times),
rep(rnorm(n.subj/2, 0, 10), each=n.times), rep(rnorm(n.subj/2, 0, 25),
each=n.times), rep(rnorm(n.subj/2, 0, 20), each=3))
for (i in 1:dim(simdata)[1]) simdata$y[i] <- rnorm(1,
randeff[i]+simdata$time[i]+simdata$coord[i]+10*simdata$group[i], 16)
lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+time:group)
+ (0+coord1+coord2+coord1:group+coord2:group|ID), data=simdata)
lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+time:group)
+ (0+gp0:coord1+gp0:coord2|ID)+(0+coord1:group+coord2:group|ID),
data=simdata),
Many thanks for any help.
Best regards,
Sarah
--
Sarah Barry, MSc
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