如何从 lme4 中提取公式的随机效应部分

Question

假设我已经安装了这样的模型

mymodel <- lmer(Y~X1+(1|fac1)+(1|fac2),mydata)

如何仅提取公式中的随机效应部分 ((1|fac1)+(1|fac2))？

我知道我能做到

formula(mymodel)[-2]

但这只是返回X1 + (1| fac1) + (1| fac2)

我知道我可以用正则表达式做一些事情，但我希望有更简单的方法。

Answer 1

不需要正则表达式，但是，它仍然是字符串操作。

# stringsplit the output of your formula()
# remove the first entry 
# remove spaces with gsub()
# paste it back together

inp <- "X1 + (1| fac1) + (1| fac2)"

paste(gsub(" ", "", unlist(strsplit(inp, "+", fixed = T))[-1], fixed = T), 
      collapse = " + ")

# [1] "(1|fac1) + (1|fac2)"

Answer 2

一个可能无法概括的简单解决方案：

# This model may not make much sense, just for reproducibility
mymodel <- lmer(Petal.Length~Sepal.Width+(1|Species) + (1|Petal.Width),iris)    
stringr::str_extract_all(formula(mymodel),"\\(.*\\)")[3]
    [[1]]
    [1] "(1 | Species) + (1 | Petal.Width)"

“自动”删除所有空元素：

purrr::compact(stringr::str_extract_all(formula(mymodel),"\\(.*\\)"))
[[1]]
[1] "(1 | Species) + (1 | Petal.Width)"

Answer 3

查找栏

lme4 包提供findbars:

library(lme4)

fo <- Y~X1+(1|fac1)+(1|fac2)

findbars(fo)
## [[1]]
## 1 | fac1
##
## [[2]]
## 1 | fac2

如果需要字符串，我们可以使用以下。 deparse1 将处理 deparse 失败的某些罕见情况，但如果有必要在 R 4.0.0 之前的 R 版本中工作，deparse 将主要作为替代方案。

sapply(findbars(fo), deparse1)
## [1] "1 | fac1" "1 | fac2"

如果所需的结果是公式的 RHS 但没有固定效应项，那么我们可以通过添加括号并使用 reformulate 来重构上述结果。如果需要公式对象，请省略 [[2]]。以上关于deparse1 的讨论也适用于此处。

reformulate(sprintf("(%s)", sapply(findbars(fo), deparse1)))[[2]]
## (1 | fac1) + (1 | fac2)

术语/标签

获取字符结果的另一种方法是使用labels，它将从terms 中提取它们。如果需要公式，请使用reformulate，如上所述。这不使用任何包。

X <- grep("|", labels(terms(fo)), fixed = TRUE, value = TRUE)
X
## [1] "1 | fac1" "1 | fac2"

如上，公式及其右侧可以从X 生成，如下所示：

reformulate(sprintf("(%s)", X))
reformulate(sprintf("(%s)", X))[[2]]

getTerms

另一种方法是使用Terms of a sum in a R expression 中的getTerms 这个简短的函数递归地遍历公式以提取术语。它不使用任何包。

XX <- grep("|", sapply(getTerms(fo[[3]]), deparse1), fixed = TRUE, value = TRUE)
XX
## [1] "(1 | fac1)" "(1 | fac2)"

公式及其右边可以这样生成：

reformulate(XX)
reformulate(XX)[[2]]

Answer 4

看到 G Grthendieck 的回答后，我意识到我可能正在重新发明轮子，但这里有一种方法可以在不使用正则表达式的情况下从模型中取出随机效果部分。它使用递归来检查公式的 AST 中的每个调用，并只保留括号中的调用，然后将其重新构建为表达式。我可能错了，但这感觉比在字符串和语言对象之间切换更安全。可以修改为只退出| 调用。

get_random_effects <- function(mod)
{
  rip_formula <- function(form) 
  {
    
    if(rlang::is_formula(form)) form <- as.list(form)[-c(1:2)][[1]]
    if(is.call(form)) {
      call_list <- as.list(form)
      if(as.character(call_list[[1]]) == "+") 
        return(unlist(lapply(call_list[-1], rip_formula)))
      if(as.character(call_list[[1]]) == "(") 
        return(form)
     } 
  }
  
  re_list <- rip_formula(formula(mod))
  while(length(re_list) > 2) 
    re_list <- c(as.call(list(bquote(`+`), re_list[1:2])), re_list[-(1:2)])
  as.call(list(bquote(`+`), re_list[[1]], re_list[[2]]))
}

所以现在很简单：

get_random_effects(mymodel)
#> (1 | fac1) + (1 | fac2)

Answer 5

您可以使用insight-package 访问各种模型信息，例如公式、预测变量、数据等。insight 提供适用于许多不同模型的类型安全“泛型”。在这种情况下，您可以使用find_formula() 或find_random()：

library(insight)
library(lme4)
data(sleepstudy)
sleepstudy$mygrp <- sample(1:5, size = 180, replace = TRUE)
sleepstudy$mysubgrp <- NA
for (i in 1:5) {
  filter_group <- sleepstudy$mygrp == i
  sleepstudy$mysubgrp[filter_group] <-
    sample(1:30, size = sum(filter_group), replace = TRUE)
}

m <- lmer(
  Reaction ~ Days + (1 | mygrp / mysubgrp) + (1 | Subject),
  data = sleepstudy
)
#> boundary (singular) fit: see ?isSingular

find_formula(m)
#> $conditional
#> Reaction ~ Days
#> 
#> $random
#> $random[[1]]
#> ~1 | mysubgrp:mygrp
#> 
#> $random[[2]]
#> ~1 | mygrp
#> 
#> $random[[3]]
#> ~1 | Subject

find_random(m)
#> $random
#> [1] "mysubgrp:mygrp" "mygrp"          "Subject"       
#> 

find_random(m, split_nested = TRUE)
#> $random
#> [1] "mysubgrp" "mygrp"    "Subject" 

find_random(m, split_nested = TRUE, flatten = TRUE)
#> [1] "mysubgrp" "mygrp"    "Subject" 

find_formula() 和 find_random() 也适用于具有随机效应的零通胀部分的模型，例如对于 glmmTMB 或 brms 包中的模型。寻找可能的随机斜率的“对应方”是find_random_slopes()。