How do I extract random effects for a particular group from an R lme object?

Question

Hello,

I have successfully fit a linear mixed-effects model, and I'm looking to extract the random effect component for individual groups. I know that the full list of random effects can be extracted using

random.effects(model)

Then print(random.effects(model)) gives a two-column list of group names and random effects, even though the data itself appears to have only one column. My question is whether it is possible to "look up" the random effect associated with a particular group by the group name, or, if not, how I can find a list of group names in the same order as the random effects in the data frame that is output by random.effects().

Thanks
Mark Ch.

Answer 1

> library(nlme)
> d <- data.frame(x=rep(letters, each=5),
                z=rep(LETTERS[1:13], each=10),
                y=rep(rnorm(26, sd=2), each=5) + rep(rnorm(13), each=10) + rnorm(26 * 5))
> r <- ranef(d)   # random.effects is a synonym for this
# Look at the structure of r
> str(r)
List of 2
 $ z:'data.frame':  13 obs. of  1 variable:
  ..$ (Intercept): num [1:13] -1.575 -0.365 -1.817 0.235 2.369 ...
  ..- attr(*, "effectNames")= chr "(Intercept)"
 $ x:'data.frame':  26 obs. of  1 variable:
  ..$ (Intercept): num [1:26] -0.8628 0.0536 1.724 -1.9115 -1.1764 ...
  ..- attr(*, "effectNames")= chr "(Intercept)"
 - attr(*, "label")= chr "Random effects"
 - attr(*, "level")= int 2
 - attr(*, "standardized")= logi FALSE
 - attr(*, "grpNames")= chr [1:2] "z" "x %in% z"
 - attr(*, "class")= chr [1:2] "ranef.lme" "list"
> head(r$x)
    (Intercept)
A/a -0.86283867
A/b  0.05360748
B/c  1.72401850
B/d -1.91145501
C/e -1.17643222
C/f  0.24315559
> head(r$z)
  (Intercept)
A  -1.5752441
B  -0.3648627
C  -1.8167101
D   0.2353324
E   2.3685118
F  -1.7544619
> r$z["A", ]
[1] -1.575244
> r$x["A/a", ]
[1] -0.8628387

Answer 2

Is this what you're looking for?

> library(nlme)    
> fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
            data = Loblolly,
            fixed = Asym + R0 + lrc ~ 1,
            random = Asym ~ 1,
            start = c(Asym = 103, R0 = -8.5, lrc = -3.3))

> str(random.effects(fm1))
Classes ‘ranef.lme’ and 'data.frame':   14 obs. of  1 variable:
 $ Asym: num  -5.57 -5.02 -1.69 -2.36 -2.98 ...
 - attr(*, "effectNames")= chr "Asym"
 - attr(*, "label")= chr "Random effects"
 - attr(*, "level")= int 1
 - attr(*, "standardized")= logi FALSE
 - attr(*, "grpNames")= chr "Seed"
> random.effects(fm1)$Asym
 [1] -5.5654676 -5.0168202 -1.6920794 -2.3587798 -2.9814647 -1.4018554
 [7] -0.1100587 -2.3613150  1.1947892  2.0119121  2.9862349  3.5890462
[13]  4.6094776  7.0963810

Answer 3

The problem, it turns out, was my particular way of indexing the groups. ranef(lme) returns a data frame where the row names are the group names. In my data, I used a very long number to differentiate between groups, which R rounded to a handful of decimal places. That meant it was impossible to exactly reference individual groups by name.

I solved the problem by converting every index to a base-62 number. I used digits and the lower and upper case alphabet as the set of characters in the number. (That is, the number matched [a-zA-Z0-9]*) This dramatically reduced the length of the group name and made it impossible for R to round the group name -- the more characters you use, the shorter it gets.

Now, if I do:

M3.ranef <- ranef(M3)
x <- M3.ranef[group_ID,1]

x is the random effect for the group named group_ID, which is how data frames are supposed to work.