[R-sig-ME] Different random effects variances for outcomes and groups
David Afshartous
dafshartous at med.miami.edu
Wed Sep 26 20:43:11 CEST 2007
Hi Sarah,
comments below; please excuse the e-mail length, I'm trying to tie up a
previously related thread.
Consider: 1) random term for each child, stratified for some partitioning.
2) random term for each child, PLUS random terms for other grouping factors,
e.g., group and/or coordinate (whether this is a good ideas is a separate
issue). In terms of the mathematical equation, 1) involves only 1 random
term (aside from the lowest level residual error term), whose variance needs
subscripts for the desired stratification. On the other hand, 2) will
involve multiple random terms (w/o subscripts for their variance; one could
attempt this, but let's assume not for now). Moreover, 2) does not
accomplish stratification, i.e., differential variability for some defined
partition.
Now, with respect to your three models (below I'll number them from the
beginning to avoid e-mail thread confusion), I don't think any of them has a
random effect for group. This would require a "(... | group)" term. All of
the models have the random effect defined by "(...| ID)"; thus the random
term is defined per child (ID), and the "(...|" part defines a variance
stratification of this random term. However, for this to make sense, I think
that the stratification must uniquely partition the data (comments below for
each model on what I mean by this):
Model 1:
Barryfm1
=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)
Data: simdata
AIC BIC logLik MLdeviance REMLdeviance
10702 10794 -5333 10687 10666
For this model the random statement is:
(0+coord1+coord2+coord1:group+coord2:group|ID)
One way to assess what this does is to sum the left hand side of "(...| ID)"
for a particular observation. For the partition and hence variance
stratification to work, I think that this sum should always be equal to 1.
Otherwise, I think the stratification will be redundant. Here, the sum will
equal 1 or 2, e.g., all observations in the second group will have sum equal
to 2. Now, regarding the random effects,
> ranef(Barryfm1)
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
...
there are 4 components for each child. Recall that there is really only 1
random intercept term and this is being stratified for a partition, and it
looks like this partition is redundant. For instance, one might be tempted
to say that the 4 components above correspond to the 4=2x2 cells defined by
coord and group, and that for child1 we have an estimated intercept random
effect of -0.502182860 under coord1 and group 0, and 7.98888012 + 5.6973871
when under coord2 and group 1; however, it doesn't seem that the indicator
variables are defined that way, since coord1 is not defined as specific or
"baseline" to group 0. Or is perhaps lmer imposing this? (I calculated the
variance of these estimated random effects and compared them w/ the
estimated model random effect variances and they didn't seem to "match" ...
To be sure, they shouldn't be the same, but one would except some type of
correspondence.)
Model 2:
Barryfm2 =
lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+t
ime:group)
+ (0+gp0:coord1+gp0:coord2|ID)+(0+coord1:group+coord2:group|ID),
data=simdata)
Data: simdata
AIC BIC logLik MLdeviance REMLdeviance
10694 10765 -5333 10687 10666
This model statement seems to accomplish the desired partitioning, viz., the
two "(...| ID)" components will have the left hand side sum to 1 for all
observations and partition the observations to your desired stratification.
Moreover, the way you have it gives you the covariance for the coordinates
within the same group that you wanted. Finally, as mentioned before, the
estimated random effects exhibit variability that seems to correspond to the
estimated model random effects variances.
Model 3:
Barryfm3=
lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+t
ime:group)
+ (0+coord1+coord2+group|ID), data=simdata)
Data: simdata
AIC BIC logLik MLdeviance REMLdeviance
10697 10768 -5335 10690 10669
For this model the random statement is "(0+coord1+coord2+group|ID)".
Similar to comments for model 1, this doesn't provide a clean partitioning
to yield a variance stratification, nor does it have a random effect for
group. Random effects are:
> ranef(Barryfm3)
An object of class ³ranef.lmer²
[[1]]
coord1 coord2 group
1 -1.596120e+00 16.48042299 -3.145684833
2 -6.500204e-01 6.89408428 -1.243928305
...
As with 1), I'm not sure how to interpret this stratification. Once again,
there a random intercept defined by the "(...| ID)" or child grouping, but
this is stratified by "...|" in a way that doesn't seem to make sense, yet
we have estimates of random effects. As w/ 1), based on the way the
indicators are defined, I don't think we can say that the estimated random
effect for child 1 in coord1 and group 0 is -1.596120e+00, and 16.48042299
-3.145684833 for coord2 and group 1.
In summary, it would be interesting to get closure on what is going on in 1)
and 3), and why the logLik is so similar to 2), and whether lmer should be
producing any estimates at all if the random intercept variance
stratification is really redundant as I suggest. Perhaps Prof. Bates could
shed light on the latter issue.
Cheers,
David
On 9/26/07 5:57 AM, "Sarah Barry" <sarah at stats.gla.ac.uk> wrote:
> David,
>
> Perhaps that's actually what I was fitting before:
>
> (2)
> 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)
>
> with a separate group effect for each coordinate. An alternative with
> the same group effect across coordinates would be:
>
> (3)
> lmer2(y~-1+coord1+coord2+coord1:(time+group+time:group)+coord2:(time+group+tim
> e:group)
> + (0+coord1+coord2+group|ID), data=simdata)
>
> Model (2) gives the same loglik as model (1) (the desired model with
> stratified variances), but a higher AIC. Model (3) gives a slightly
> lower loglik, with the same AIC as model (2). BIC is a bit more
> varying, being lowest for model (1) and highest for model (2).
>
> Sarah
>
> David Afshartous wrote:
>
>> One additional comment below ...
>>
>>
>> On 9/25/07 10:03 AM, "David Afshartous" <dafshartous at med.miami.edu> wrote:
>>
>>
>>
>>> 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.
>>>
>>>
>>>
>>
>> Moreover, allowing a random effect for group is definitely not what you want
>> as this would not model the differential variability in your data. Each
>> child would have a single group random effect (distributed the same for each
>> group), and there would be a child random effect that depends on the
>> coordinate, but they would not be distributed differently across the groups
>> (you would get 1 list instead of two lists for the random effects in your
>> second model statement). Having said that, I estimated both cases for some
>> simple repeated measures data that has a treatment factor. Specifically:
>> 1) Stratified RE variance: fixed effects for treatment, RE for each subject,
>> stratified per treatment group.
>> 2) RE for treatment group: fixed effects for treatment, RE for each subject
>> (not stratified per treatment), plus RE for treatment group.
>> Although I really don't want option 2), I was curious about the differences
>> in model fit. Both models were extremely similar for AIC/BIC/LogLik, which
>> is surprising since 2) doesn't model at all the differential growth curve
>> variability per treatment in my data. I didn't try this comparison for your
>> data as I was unsure of how to modify your lmer2 statement (perhaps mine was
>> much simpler due to the single factor of treatment). If you are able to do
>> the appropriate modification I'd be interested 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
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>>
>>
>>
>>
More information about the R-sig-mixed-models
mailing list