如何從lme4中提取隨機效應和方差組件dlply

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 

這工作相當好，因爲所有的隨機效應是截距。如果在模型中有一些隨機斜率項，您可能想單獨使用隨機效應，或使用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）...