使用 by/group by 语句运行 lmer？

Question

我试图找到一种运行 lmer 模型的快速方法，但为每个分组变量单独运行它(在 SAS 中可以使用 by= 语句)。我已经尝试使用 dplyr 使用我发现的代码:

t1<- mod1 %>% group_by(c) %>% do(function(df){lmer(m1.formula,data=df)})

但这似乎不起作用。

任何人都知道如何使用 dplyr 或其他方法来做到这一点？

Answer 1

library("lme4")
data(Orthodont,package="nlme")

您可能需要在这里考虑两个基本问题:

统计:如上所述，考虑在数据集中的每个层(分组变量)上单独运行混合模型有点​​奇怪。通常混合模型的全部目的是将模型拟合到组合数据集，尽管我当然可以想象异常(exception)(下面，我按性别拟合单独的混合模型)。您可能正在寻找类似 lmList 的东西。函数(nlme 和 lme4 都有版本)，它分别在每个层上运行(广义)线性模型(不是混合模型)。这更有意义，尤其是作为一种探索性技术。

计算:专门做您在 dplyr 中要求的操作框架有点难，因为基本dplyr范式假设您正在对数据帧(或数据表)进行操作，可能已分组并返回数据帧。这意味着每个操作返回的位必须是数据帧(而不是例如 merMod 模型对象)。 (@docendodismus 指出你可以通过在下面的代码中指定 do(model = ...) 来实现，但我认为结果对象的结构有点奇怪，会鼓励你重新思考你的问题，如下图所示)

在基础 R 中，您可以这样做:

lapply(split(Orthodont,Orthodont$Sex),
       lmer,formula=distance~age+(1|Subject))

或者

by(Orthodont,Orthodont$Sex,
       lmer,formula=distance~age+(1|Subject))

题外话 :如果您想为每个主题拟合线性(非混合)模型，您可以使用

## first remove 'groupedData' attributes from the data, which don't
## work with lme4's version of lmList
Orthodont <- data.frame(Orthodont) 
lmList(distance~age|Subject,Orthodont)
## note: not including Sex, which doesn't vary within subjects

返回主线程 : 在 plyr ( dplyr 的祖先)框架，您可以更紧凑地按性别拟合单独的混合模型:

library("plyr")
dlply(Orthodont,.(Sex),
       lmer,formula=distance~age+(1|Subject))

detach("package:plyr")

如果你想在 plyr 做，你好像需要 do() (我以为我可以不用它，但我还没有找到方法)，并且您需要创建一个函数，将摘要作为数据框返回。

library("dplyr")

Orthodont %>% group_by(Sex) %>%
    do(lmer(.,formula=distance~age+(1|Subject)))

产生

## Error: Results are not data frames at positions: 1, 2

你可以这样做:

 lmer_sum <- function(x,...) {
      m <- lmer(x,...)
      c(fixef(m),unlist(VarCorr(m)))
        data.frame(rbind(c(fixef(m),unlist(VarCorr(m)))),
                   check.names=FALSE)
 }

( unlist(VarCorr(m)) 给出单个标量随机效应的 RE 方差；需要整个 data.frame(rbind(...)) 东西才能将数值向量转换为一行数据框； check.names=FALSE 保留列名 (Intercept) )

 Orthodont %>% group_by(Sex) %>%
    do(lmer_sum(.,formula=distance~age+(1|Subject)))

这给出了合理的结果。