后续:从 lme4 中的 VarCorr 对象中提取名称并将其粘贴为列名称

Follow-up: Extracting names from a VarCorr object in lme4 and pasting it as column names

我正在跟进 。下面的函数 foo 获取 VarCorr(fit) 输出的 Name 列,并使它们成为 summary(rePCA(fit)) 调用的列名称。

当我们输入fm1fm2时可以正常工作,但我想知道为什么输入fm3会失败?有修复吗?

library(lme4) 
dat <- read.csv('https://raw.githubusercontent.com/rnorouzian/e/master/sng.csv')
fm1 <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin) 
fm2 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
fm3 <- lmer(y ~ A * B * C + (A + B | group) + (C|group), data = dat)


foo <- function(fit) {
  
  obj <- summary(rePCA(fit))
  model <- VarCorr(fit)
  
  Map(function(x, z) {
    colnames(x$importance) <- paste(z, unique(sapply(model, colnames)), sep = '_')
    x
  }, obj, names(obj))
}

#EXAMPLE OF USE:
foo(fm1) ###     OK !
foo(fm2) ###     OK !
foo(fm3) ###     :-( Error in dimnames(x) <- dn

objmodel 的长度不同时函数失败。这是使其适用于 fm3.

的技巧
foo <- function(fit) {
  
  obj <- summary(rePCA(fit))
  model <- VarCorr(fit)
  if(length(obj) == length(model)) {
   obj <- Map(function(x, z) {
    colnames(x$importance) <- paste(z, unique(sapply(model, colnames)), sep = '_')
    x
  }, obj, names(obj))
  }
  else if(length(obj) == 1) {
    colnames(obj[[1]]$importance) <- unlist(mapply(paste, names(model), sapply(model, colnames), MoreArgs = list(sep = '_')))
  }
  return(obj)
}

这returns输出如下:

foo(fm1) 

#$plate
#Importance of components:
#                       plate_(Intercept)
#Standard deviation                  1.54
#Proportion of Variance              1.00
#Cumulative Proportion               1.00

#$sample
#Importance of components:
#                       sample_(Intercept)
#Standard deviation                   3.51
#Proportion of Variance               1.00
#Cumulative Proportion                1.00

foo(fm2) 
#$Subject
#Importance of components:
#                       Subject_(Intercept) Subject_Days
#Standard deviation                   0.967       0.2309
#Proportion of Variance               0.946       0.0539
#Cumulative Proportion                0.946       1.0000

foo(fm3) 
#$group
#Importance of components:
#                       group_(Intercept) group_A group_B group.1_(Intercept) group.1_C
#Standard deviation                 1.418   1.291  0.5129              0.4542 0.0000497
#Proportion of Variance             0.485   0.402  0.0634              0.0498 0.0000000
#Cumulative Proportion              0.485   0.887  0.9502              1.0000 1.0000000

您可以从 'fit@cnms' 中获取列名,这样就省去了使用 'VarCorr' 的麻烦。魔鬼似乎是 fm1 的情况,它给出了一个列表输出,我们可能想要强制 as.data.frame。那么我们可以只用`colnames<-`,不需要Map战斗。

foo2 <- function(fit) {
  obj <- summary(rePCA(fit))
  obj <- as.data.frame(lapply(obj, `[`, "importance"))
  `colnames<-`(obj, paste(names(fit@cnms), unlist(fit@cnms), sep="_"))
}

foo2(fm1)
#                        plate_(Intercept) sample_(Intercept)
# Standard deviation              1.539676           3.512519
# Proportion of Variance          1.000000           1.000000
# Cumulative Proportion           1.000000           1.000000

foo2(fm2)
#                        Subject_(Intercept) Subject_Days
# Standard deviation                0.966868    0.2308798
# Proportion of Variance            0.946050    0.0539500
# Cumulative Proportion             0.946050    1.0000000

foo2(fm3)
#                        group_(Intercept)  group_A   group_B group_(Intercept)    group_C
# Standard deviation              1.385987 1.322335 0.5128262         0.4547251 0.08892506
# Proportion of Variance          0.463190 0.421630 0.0634100         0.0498600 0.00191000
# Cumulative Proportion           0.463190 0.884820 0.9482300         0.9980900 1.00000000

我发现的另一个优点是数字未四舍五入。您可以在函数中构建 rounding 或之后执行此操作:

round(foo2(fm1), 2)
#                        plate_(Intercept) sample_(Intercept)
# Standard deviation                  1.54               3.51
# Proportion of Variance              1.00               1.00
# Cumulative Proportion               1.00               1.00