如何从包装在 dlply 中的 lme4 中提取随机效应和方差分量

Question

这篇帖子How can I extract elements from lists of lists in R? 回答了我的一些问题，但这对我来说仍然不太有效，而且我需要做的事情超出了我的 R 知识范围。

我有 2 个环境（=试验）、2 年和 5 个感兴趣的特征（由 trait_id 定义）的现场试验数据。 GID 是唯一的线路标识符。我在 lme4 中的模型是：

mods <- dlply(data,.(trial,trait_id), 
 function(d)
  lmer(phenotype_value ~(1|GID)+(1|year)+(1|year:GID)+(1|year:rep), 
       na.action = na.omit,data=d))

运行它会返回一个包含 10 个元素的大型列表，我想将每次试验的所有特征的 GID 随机效应存储在数据框中。我尝试了几件事：

blup=lapply(mods,ranef, drop = FALSE)
blup1=blup[[1]]
blup2=blup1$GID  

会给我一个 df，其中包含每次试验一个特征的随机效应，我希望有更精简的东西，可以在列名中保留一些信息，如 $irrigation.GRYLD。

这是一个只有两个特征（GRYLD、PTHT）、2 年（11OBR、12OBR）和两个代表的可重现示例：

structure(list(GID = structure(c(1L, 2L, 3L, 4L, 5L, 5L, 1L, 
2L, 4L, 3L, 1L, 2L, 3L, 4L, 5L, 5L, 1L, 2L, 4L, 3L, 1L, 2L, 3L, 
4L, 5L, 5L, 2L, 1L, 4L, 3L, 1L, 2L, 3L, 4L, 5L, 5L, 2L, 1L, 4L, 
3L, 1L, 2L, 3L, 4L, 5L, 5L, 1L, 2L, 4L, 3L, 1L, 2L, 3L, 4L, 5L, 
5L, 1L, 2L, 4L, 3L, 1L, 2L, 3L, 4L, 5L, 5L, 2L, 1L, 4L, 3L, 1L, 
2L, 3L, 4L, 5L, 5L, 2L, 1L, 4L, 3L), .Label = c("A", "B", "C", 
"D", "E"), class = "factor"), year = structure(c(1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("11OBR", 
"12OBR"), class = "factor"), trial = structure(c(1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("heat", 
"irrigation"), class = "factor"), rep = c(1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L), trait_id = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("GRYLD", 
"PTHT"), class = "factor"), phenotype_value = c(3.93, 3.38, 1.65, 
4.33, 2.45, 2.48, 3.98, 3.3, 4.96, 1.53, 87.5, 69.5, 65.5, 84.5, 
77, 81, 94.5, 84.5, 89, 81, 6.56, 4.3, 5.76, 7.3, 5.73, 4.14, 
5.93, 6.96, 8.43, 5.81, 114.5, 100, 104.5, 110, 110, 106, 99, 
97.5, 105, 100, 0.119, 0.131, 0.681, 0.963, 0.738, 1.144, 0.194, 
0.731, 0.895, 0.648, 35, 50, 45, 50, 45, 50, 55, 45, 50, 55, 
2.79, 3.73, 3.96, 4.64, 5.03, 2.94, 3.78, 4.14, 3.89, 3.21, 90, 
95, 105, 100, 105, 85, 95, 100, 100, 95)), .Names = c("GID", 
"year", "trial", "rep", "trait_id", "phenotype_value"), class = "data.frame", row.names = c(NA, 
-80L))

Answer 1

我不太很确定你想要什么作为输出格式，但是怎么样：

all_ranef <- function(object) {
    rr <- ranef(object)
    ldply(rr,function(x) data.frame(group=rownames(x),x,check.names=FALSE))
}
ldply(mods,all_ranef)

##         trial trait_id      .id   group   (Intercept)
## 1        heat    GRYLD year:GID 11OBR:A  7.935352e-01
## 2        heat    GRYLD year:GID 11OBR:B  1.960487e-01
## 3        heat    GRYLD year:GID 11OBR:C -1.504116e+00
## ...
## 82 irrigation     PTHT year:rep 12OBR:2 -1.595022e+00
## 83 irrigation     PTHT     year   11OBR  2.915033e+00
## 84 irrigation     PTHT     year   12OBR -2.915033e+00

这工作得相当好，因为您的所有随机效果都是仅拦截的。如果模型中有一些随机斜率项，您可能想要reshape2:::melt() 单个随机效应，或者使用rbind.fill() 将数据帧与不同的随机效应列组合起来。

library("ggplot2"); theme_set(theme_bw())
ggplot(vals, aes(y=group,x=`(Intercept)`))+
    geom_point(aes(colour=interaction(trial,trait_id)))+
    facet_wrap(~.id,scale="free")

顺便说一句，通常不建议使用只有 2 个级别 (YEAR) 的因子作为分组变量...