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

[David Afshartous]

Hi Sarah,
A few comments below.
David


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

> Hi David,
> 
> Thanks for your reply.  I comment on your various questions below.
> 
> David Afshartous wrote:
> 
>> 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).
>> 
>>  
>> 
> Yes, you are right here (apart from a typo - there are 100 x 2 x 3 x 2 =
> 1200 observations).  The groups do contain different children.
> 
>> 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.
>> 
>>  
>> 
> I think you're right here and this corresponds with what you say below.
> I think my model should be instead written as
> 
> y_{ir}(t) = \beta_{0r} + \beta_{1r} group_i + \beta_{2r} t + \beta_{3r}
> group_i : t  + group_i * b_{irp} + \epsilon_{ir}(t),
> 
> where p=0,1 represents the value of 'group_i' and is distributed as you
> suggest.  This seems to make more sense than everyone having a random
> intercept from one distribution and group 1 having an additional one.
> This way they are completely separate.
> 
>> 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?).
>> 
>>  
>> 
> Yes, I should have considered this more carefully.  Certainly the random
> effects for coordinates 1 and 2 for a particular individual should have
> been simulated from a multivariate normal distribution rather than from
> two separate normals.  I don't think that there is so much a group
> component and a coordinate component, as a coordinate component that
> differs depending on which group the individual comes from.  So, for
> example, instead of the randeff statement below we could have (repeated
> across times and inserted in the correct positions for each group):
> 
> randeff.gp0 <- mvrnorm(n.subj/2, c(0,0), matrix(c(100,50,50,400),
> nc=2))  # coordinates 1 and 2 (respectively) for individuals in group 0 #
> randeff.gp1 <- mvrnorm(n.subj/2, c(0,0), matrix(c(225,70,70,625),
> nc=2))  # coordinates 1 and 2 (respectively) for individuals in group 1 #
> 
> I think the group itself should only be a fixed effect because, while
> individuals are randomly sampled from two populations of cleft-lip
> patients and healthy controls, the groups are fixed (not sampled from a
> population of groups).  What do you think?  This has been one point that
> I have wondered long and hard about.
> 

I agree that group should be a fixed effect based on your description. In
your current model, group is a fixed effect, and the random effect on the
intercept is for child ID, although this depends on both the group and
coordinate for the child.  You do not have a random effect for group per se.
This can be contrasted w/ instead employing separate random effects for the
nesting levels (Pinheiro & Bates p.40 is a good example of a nesting
structure where random effects are used for each level of the nesting.)  It
would be interesting to observe the difference in the results.


>> 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
>>  
>> 
> I do agree and I'm not sure now what model the first lmer2 model
> statement is actually fitting, because I would have expected at least
> some of the coord1:group and coord2:group random effects to be zero.
> The second lmer2 statement gives more sensible answers and corresponds
> with the rewritten model above.  Also when I check the estimated random
> effects variances and covariances against the actual values, there is
> good correspondence so I think I will proceed with this parameterisation.
> 
> Many thanks for your help,
> Sarah
> 
>> 
>> 
>> 
>> 
>> 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+t
>>> im
>>> 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+t
>>> im
>>> 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+t
>>> im
>>> 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+t
>>> im
>>> e:group) 
>>> + (0+gp0:coord1+gp0:coord2|ID)+(0+coord1:group+coord2:group|ID),
>>> data=simdata),
>>> 
>>> Many thanks for any help.
>>> 
>>> Best regards,
>>> Sarah
>>>    
>>> 
>> 
>>  
>>