开发者

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

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 on开发者_StackOverflowe 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.


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.


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


> 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
0

上一篇:

下一篇:

精彩评论

暂无评论...
验证码 换一张
取 消

最新问答

问答排行榜