Skip to content

random effects specification

16 messages · Douglas Bates, Jesse Whittington, Ken Beath +2 more

Original
Sébastien Bihorel
# 
Hi,

In the past I've used lme to fit simple mixed models to longitudinal
data (growth), but now I'm trying to learn lmer and its different syntax
to fit a more complex model.  I have a structure with subjects (id,
random factor) exposed to 4 different treatments and a continuous
response variable is measured (n).  The subjects come from 2 different
communities, so it's a nested design very much like the Oats data in the
nlme package.  The interest is in the fixed effects of community and
treatment, and their interaction, so I thought this could be modelled in
lmer with this call:

lmer(n ~ treatment + community + (1 | id/treatment), mydata)

but got this error:

Error in lmerFactorList(formula, mf, fltype) : 
  number of levels in grouping factor(s) ?treatment:id? is too large


Am I using the right formula here?  Thanks.
Douglas Bates
#
On Thu, Apr 3, 2008 at 3:45 PM, Sebastian P. Luque <spluque at gmail.com> wrote:
It seems that the observations are indexed by subject and treatment so
the number of levels in the factor treatment:id equals the number of
observations.  You can't estimate a variance for such a term and also
estimate a residual variance.

I would start with

n ~ treatment * community +(1|id)
Sébastien Bihorel
# 
On Thu, 3 Apr 2008 17:32:40 -0500,
"Douglas Bates" <bates at stat.wisc.edu> wrote:
[...]
Yes, the observations are indexed by subject and treatment in the sense
that id levels are the same within treatments of the same community, but
are different among communities.  This is a subset of the data:

---<---------------cut here---------------start-------------->---
id  community treatment     n
 1          A         1 13.93
 2          A         1 14.42
 3          A         1 13.56
 1          A         2 14.61
 2          A         2 14.74
 3          A         2 15.59
 1          A         3 13.95
 2          A         3 15.21
 3          A         3 14.51
 1          A         4 13.61
 2          A         4 14.99
 3          A         4 15.13
 4          B         1 14.79
 5          B         1 13.41
 6          B         1 14.71
 4          B         2 14.69
 5          B         2 13.46
 6          B         2 14.28
 4          B         3 14.30
 5          B         3 13.18
 6          B         3 13.58
 4          B         4 14.54
 5          B         4 13.25
 6          B         4 14.09
---<---------------cut here---------------end---------------->---

Of course, there are many more individuals, but the levels of id differ
among communities, and are the same among treatments.  lmer did converge
rapidly with your suggested formula though:

---<---------------cut here---------------start-------------->---
Linear mixed-effects model fit by REML 
Formula: n ~ treatment * community + (1 | id) 
   Data: isotope.m.ph 
 AIC BIC logLik MLdeviance REMLdeviance
 450 481   -216        410          432
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.193    0.439   
 Residual             0.232    0.481   
number of obs: 240, groups: id, 61

Fixed effects:
                               Estimate Std. Error t value
(Intercept)                     14.9748     0.1170   128.0
treatment2                       0.0884     0.1222     0.7
treatment3                      -0.2829     0.1222    -2.3
treatment4                       0.3568     0.1222     2.9
communitysanikiluaq             -0.5749     0.1678    -3.4
treatment2:communitysanikiluaq  -0.2471     0.1763    -1.4
treatment3:communitysanikiluaq  -0.7479     0.1763    -4.2
treatment4:communitysanikiluaq  -0.6169     0.1763    -3.5

Correlation of Fixed Effects:
            (Intr) trtmn2 trtmn3 trtmn4 cmmnty trtm2: trtm3:
treatment2  -0.522                                          
treatment3  -0.522  0.500                                   
treatment4  -0.522  0.500  0.500                            
cmmntysnklq -0.697  0.364  0.364  0.364                     
trtmnt2:cmm  0.362 -0.693 -0.347 -0.347 -0.524              
trtmnt3:cmm  0.362 -0.347 -0.693 -0.347 -0.524  0.503       
trtmnt4:cmm  0.362 -0.347 -0.347 -0.693 -0.524  0.503  0.503
---<---------------cut here---------------end---------------->---


However, I don't understand how (1 | id) accounts for treatment being
nested within community.  Maybe it's time for me to re-read some more
examples from "Mixed-effects models in S and S-plus".  Thanks.


Cheers,
Douglas Bates
#
On Thu, Apr 3, 2008 at 6:32 PM, Sebastian P. Luque <spluque at gmail.com> wrote:
I'm not sure that I understand what you mean by "treatment being
nested within community".  Does this mean that there are really 8
different treatments because treatment 1 in community A is different
from treatment 1 in community B?  If so, then it would make sense to
me to simply create a new factor that is the interaction of treatment
and community.

Perhaps I am approaching the community factor incorrectly.  In your
data there are two communities so, even if it would be reasonable to
model community effects as random effects, that would be difficult.
With only two levels I think it is best modeled as a fixed effect,
which would mean that questions about treatment and community are
related to the fixed effects.
Sébastien Bihorel
# 
On Fri, 4 Apr 2008 07:17:36 -0500,
"Douglas Bates" <bates at stat.wisc.edu> wrote:
[...]
I was not employing the term "nested" properly.  The number of levels
for both community and treatment are 2 and 4, respectively, just as in
the example.  The same 4 treatments were used in both communties, so in
fact, treatment is crossed with community, not nested.  However,
subjects are nested within communities because each subject belongs to
one community only, yet received all 4 treatments.  Sorry for this
confusion.
Could you please show a formula for the case where each individual is
seen at both communities (community and treatment still being fixed)?
This would help me understand the syntax better.

Thank you so much for your help.
1 day later
Ken Beath
#
On 05/04/2008, at 12:05 AM, Sebastian P. Luque wrote:
Once they are considered fixed effects, concepts of crossing and  
nesting are irrelevant. They are simply covariates. So a model of the  
form n ~ treatment + community +(1|id) or if the treatment effect is  
allowed to vary between communities n ~ treatment *community +(1|id)  
is appropriate. The main problem is your subject id are not unique.  
You will need to define a new id. The easiest way is to add a  
different large number to id depending on community.
Same model as previous, provided a subject only receives a treatment  
once. If a subject receives the same treatment more than once then  
there needs to be a random effect that models the correlation between  
repeated measurements of the same treatment, so the model is  
y~treatment+community+(1|id/treatment) One problem that may have  
occurred in your original attempts is that id and treatment need to be  
factors.

Ken
Douglas Bates
#
On Sat, Apr 5, 2008 at 10:14 PM, Ken Beath <kjbeath at kagi.com> wrote:
I agree with everything up to here.
That approach contradicts your later advice to represent a factor
variable as a factor in R.  If id is a factor (as it should be) you
can't add  a large number to it.

The specification (1|treatment:id) generates unique id's.

To me the convention that different experimental units should be given
the same level of 'id' is just another nonsensical aspect of the
traditional approaches to random-effects models using observed and
expected mean squares, for which it makes sense to index the
observations by group and by unit within group.

If we could manage to unlearn old habits and just give each subject a
unique id at the start it would make life easier.
Yes, that is one way of expressing an interaction between a random
effect for id and a fixed effect for treatment.

It expands to two random effects terms (1|id) + (1|id:treatment).  The
first is the effect for person and the second is the effect of
different individuals having different responses to the levels of
treatment.

A more general model (and consequently more difficult to estimate on
occasion) has possible correlations of the random effects for
different levels of treatment within individual.  The term is written
(treatment|id).
Ken Beath
#
On 07/04/2008, at 1:04 AM, Douglas Bates wrote:
I was thinking of this in terms of the original data, and should have  
mentioned that a conversion may be required in R, as described in  
help(factor).
This does give problems fitting although I get sensible results.



The need for id to be a factor is not consistent but depends on the  
model (this is using the version of lme4 on CRAN). In the following  
simulation (for the model with subjects treated in both communities)   
id as a numeric works for (1|id) but needs to be a factor for (1|id/ 
treatment).


library(lme4)

nsubjects <- 100

id <- rep(1:nsubjects,each=8)
treatment <- rep(1:4,times=nsubjects*2)
community <- rep(1:2,each=4,times=nsubjects)

thedata <- data.frame(id,treatment,community)

subject <- data.frame(id=1:nsubjects,subject=rnorm(nsubjects))

thedata <- merge(thedata,subject)

treatment <- data.frame(id=rep(1:nsubjects, 
4),treatment=rep(1:4,each=nsubjects),subtreat=rnorm(4*nsubjects))

thedata <- merge(thedata,treatment)

thedata$y <- thedata$subject+thedata$subtreat+thedata$community+thedata 
$treatment+rnorm(8*nsubjects)/10

thedata$treatment <- factor(thedata$treatment)
thedata$community <- factor(thedata$community)

# ignoring  treatment within subject correlations
lmer0 <- lmer(y~treatment+community+(1|id),data=thedata)

print(lmer0)

thedata$id <- factor(thedata$id)

lmer1 <- lmer(y~treatment+community+(1|id/treatment),data=thedata)

print(lmer1)


Ken
Sébastien Bihorel
# 
Hi,

Fortunately, I have learnt to use unique factor levels for subjects,
only as a lucky accident though, but am glad it's a good convention to
use in mixed models.

Thanks everyone for your useful feedback, this has been very instructive!


Cheers,
1 day later
Sébastien Bihorel
# 
Hi,

I think I'm misunderstanding something because I expected a different
result below (building up a reproducible example).

---<---------------cut here---------------start-------------->---
## Simulating data
library(lattice)
library(lme4)
set.seed(1000)
rCom <- rnorm(2, mean=5, sd=0.5)
rTr <- rep(rCom / 1.1, 2)
nbase <- rnorm(240, 10, 0.1)
dta <- within(expand.grid(community=LETTERS[1:2], treatment=letters[1:4],
                          id=factor(1:30)), {
                              n <- rCom[as.numeric(community)] +
                                  rTr[as.numeric(treatment)] + nbase
                          })
dta <- dta[order(dta$community, dta$treatment), ]
## Simulate an interaction
dta$n[dta$community == "A"] <- rev(dta$n[dta$community == "A"])
## Have a look
xyplot(n ~ treatment | community, data=dta, groups=id,
       type="b", pch=19, cex=0.3)

## We fit LME with community and treatment fixed effects, their
## interactions, and use random effects for subject.
n.lmer1 <- lmer(n ~ community * treatment + (1 | id), dta)
## Remove the interaction, for comparison.
n.lmer2 <- lmer(n ~ community + treatment + (1 | id), dta)
## Compare these models.
anova(n.lmer1, n.lmer2)                 # interaction term needed

## So let's test the treatment effect at each community separately.
n.lmerA <- lmer(n ~ treatment + (1 | id), dta,
                subset=community == "A")
## Here I expected some terms to be significantly different
## from zero, but that's not the case:
summary(n.lmerA)
---<---------------cut here---------------end---------------->---

Looking at the data, there seems to be an obvious treatment effect, yet
this doesn't show in the lmer summary.  What am I misunderstanding here?
Thanks.


Cheers,
Reinhold Kliegl
# 
Looking at your data, the lmer tells you very nicely the pattern of
means shown in your plot:
Linear mixed model fit by REML
Formula: n ~ treatment + (1 | id)
   Data: dta
 Subset: community == "A"
    AIC    BIC logLik deviance REMLdev
 -181.7 -164.9  96.83   -218.5  -193.7
Random effects:
 Groups   Name        Variance   Std.Dev.
          (Intercept) 5.5211e-13 7.4304e-07
 Residual             9.8079e-03 9.9035e-02
Number of obs: 120, groups: id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 18.78058    0.01808  1038.7
treatmentb   0.33300    0.02557    13.0
treatmentc   0.01150    0.02557     0.4
treatmentd   0.32939    0.02557    12.9

Treatment b is significanlty higher, c is at the same level, and d is
again significantly higher than treatment a. (With t-values of 13 we
may say that even without HPD intervals.)

The same is true for community B, except that here treatment b and d
are significantly lower than a and that what the estimates tell you.
If you learn how to combine estimates, you can derive the pattern of
means directly from your first analyses.

Reinhold
On Tue, Apr 8, 2008 at 6:22 PM, Sebastian P. Luque <spluque at gmail.com> wrote:
Sébastien Bihorel
# 
On Tue, 8 Apr 2008 19:31:08 +0200,
"Reinhold Kliegl" <reinhold.kliegl at gmail.com> wrote:
[...]
Of course, what was I thinking!  Thanks.
You mean using approaches as those used in estimable() from the gmodels
package?  If so, is that equivalent to doing a separate analysis for
each level of community with mixed models?  I think it isn't for fixed
effects only models.


Cheers,
Reinhold Kliegl
#
On Tue, Apr 8, 2008 at 8:13 PM, Sebastian P. Luque <spluque at gmail.com> wrote:
A            B
a 18.78058 18.74912
b 19.11358 18.38756
c 18.79208 18.73067
d 19.10998 18.40533
...
Fixed effects:
                      Estimate Std. Error t value
(Intercept)           18.78058    0.01784  1052.7
communityB            -0.03146    0.02523    -1.2
treatmentb             0.33300    0.02523    13.2
treatmentc             0.01150    0.02523     0.5
treatmentd             0.32939    0.02523    13.1
communityB:treatmentb -0.69456    0.03568   -19.5
communityB:treatmentc -0.02995    0.03568    -0.8
communityB:treatmentd -0.67318    0.03568   -18.9

You can reconstruct the table of means from the coefficients:
                      A                                                   B
a :  intercept                               intercept + communityB
b:   intercept + treatmentb          intercept + treatmentb +
communityB + communityB:treatmentb
c:   intercept + treatmentc          intercept + treatmentc +
communityB + communityB:treatmentc
d:   intercept + treatmentd          intercept + treatmentd +
communityB + communityB:treatmentd

This is a consequence of the default treatment contrasts associated
with factors as attributes.
b c d
a 0 0 0
b 1 0 0
c 0 1 0
d 0 0 1
B
A 0
B 1

You can change the default, that is specify your own contrasts as
factor attributes most flexibly via C (see ?C). lmer uses the factor
attribute (like, e.g.,  lm and other programs) and provides test
statistics for each of the contrasts (and their products).

Best
Reinhold
Sébastien Bihorel
# 
[sorry for the line wrap mangling]

On Tue, 8 Apr 2008 21:47:04 +0200,
"Reinhold Kliegl" <reinhold.kliegl at gmail.com> wrote:
[...]
Yes, but IIUC, this is not necessarily equivalent to fitting separate
models with subsets of the terms/levels.  Correct me if I'm wrong, but I
think the estimates will differ when the design is unbalanced (e.g. some
subjects not receiving all treatments) because a straight calculation
from the full model assumes the coefficients have equal weight.  For
this example data, the design is balanced, so this might not be a
problem.

Thanks,
Reinhold Kliegl
# 
That is correct, but what is the issue? In general, I think it is best
to stay with one model. This gives you more precise estimates, because
you have more observations and fewer model parameters, for example,
you estimate only one residual variance. (Of course, there are also
situations where it is best to run separate analyses for different
parts of the data, for reasons of ease of communication, reviewer
requests, etc.).
This is correct. Note, however, that for unbalanced designs the model
estimates correct for differences in reliability (e.g., due to
differences in the number of observations). So, for prediction, the
estimates may be better than the observed means.

Reinhold