[R-sig-ME] Different random effects variances for outcomes and groups

[David Afshartous]

Hi Sarah,

Did you receive any replies to your post?  Some comments below, perhaps
tangential to you main question but possibly of interest.

First let me make sure I understand you data structure.  You have 100
children in each of 2 groups, and for each child you take 3 measurements at
coordinate 1 and coordinate 2.
Hence there are 100 x 2 x 3 x 3 = 1200 observations.  Moreover, the children
in the two groups are different children, hence the different IDs in the
data (as opposed to the same children being in both groups e.g. when each
child receives both treatment and placebo; this affects the number of random
effects).

In your model equation:
> 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),
> 

You have a random intercept model, where the intercept is broken down
according to fixed effects for coordinate and group, and similarly for the
slope.  I assume you want the random effect for the intercept stratified
according to both group and coordinate.  I'm not sure how the terms above
reflect this; perhaps all you need is the term b_{irp} ~ > N(0,
\sigma_{rp}2), for p=0,1; r = 1,2.

RE your simulation code to generate the random effects, I assume you have
broken "randeff" into blocks of 300 ( = 1200/4) because the group x
coordinate stratification yields 4 distinct combinations.
However, I have a question RE the way in which these random effects are
simulated.  For instance, consider patient 1 in group 1.  According to your
simulation code, two separate random normals will be generated to reflect
the random effect of this patient at coordinate 1 and coordinate 2, with
random normal variance equal to the sum of the group random effect variance
and the coordinate random effect variance.  However, I don't think this
reflects the nesting appropriately.  Perhaps the group component should only
be generated once as it is the same child in the same group; and the
coordinate component is the part that needs to be generated twice. This of
course will increase the correlation between the two realized values.
(BTW, did you chose the standard deviations values of 15, 10, 25, and 20 to
reflect the aforementioned sub-component structure, and if so what were the
standard deviations of the sub-components?).

With respect to random effects, I assume your model will generate 400 unique
random effects estimates, i.e., two (for each coordinate) for each of the
200 children.  And each of these may be viewed as the sum of the
sub-components of coordinate and group.  Running your first lmer2 model
statement yields a 200 x 4 matrix for the estimated random effects, w/ each
row being a patient and the columns corresponding apparently to the
aforementioned subcomponents:

An object of class ³ranef.lmer²
[[1]]
           coord1       coord2  coord1:group coord2:group
1    -0.502182860   7.98888012  -4.590717867    5.6973871
2    -0.190673717   3.38674017  -1.849503513    2.4098210
3    -0.981561080  12.80952815  -8.127985669    9.1788197
Etc

However, I would think some of these cells need to be 0, e.g., each patient
is only in 1 group and thus shouldn't have a random effect estimate from
both groups?  Or am I reading the table completely wrong?

Now, when I ran your second lmer2 model statement and checked the estimated
random effects (too messy to copy here), I got two lists of 2 random effects
per child (1 for each coordinate), where it appears that the two lists
correspond to the two groups, and apparently there are 0's for children that
were not in the respective group. Based on the estimated random effects
produced by the two model statements, I think that the second more
representative of what you're trying to do.  Have a look at the random
effects for the second model statement and let me know if you agree.

Cheers,
David






On 9/19/07 8:28 AM, "Sarah Barry" <sarah at stats.gla.ac.uk> wrote:

> 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+tim
> e: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+tim
> e: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+tim
> e:group) 
> + (0+coord1+coord2+coord1:group+coord2:group|ID), data=simdata)
> lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+tim
> e:group) 
> + (0+gp0:coord1+gp0:coord2|ID)+(0+coord1:group+coord2:group|ID),
> data=simdata),
> 
> Many thanks for any help.
> 
> Best regards,
> Sarah