使用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")

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

• 统计:如上所述，考虑在数据集内的每个层(分组变量)上单独运行混合模型有点奇怪。通常，混合模型的全部要点是将模型拟合到组合数据集，尽管我当然可以想象例外情况(下面，我按性别拟合单独的混合模型)。您可能正在寻找类似于lmList函数( nlme和lme4都有版本)的东西，它分别在每个层上运行(广义的)线性模型(而不是混合模型)。这更有意义，特别是当探索性technique.
• computational:在dplyr框架中具体执行您所要求的操作时，这有点困难，因为基本的dplyr范例假设您正在对数据帧(或数据表)进行操作，可能对数据帧进行分组并返回数据帧。这意味着每个操作返回的位必须是数据帧(而不是merMod模型对象)。(@docendodismus指出，您可以通过在下面的代码中指定do(model = ...)来做到这一点，但我认为结果对象的结构有点奇怪，这会鼓励您重新考虑您的问题，如下所示)

在base R中，您可以这样做：

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

或

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

Digression：如果您想要将线性(非混合)模型适合于每个主题，您可以使用

## 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)))

这给出了合理的结果。

Answer 2

问题是您错误地调用了do() -它不适用于这样的匿名函数。do()中的参数是在数据上下文中计算的，因此当您说function(df)时，do将尝试使用数据的df列。它没有该列，所以它失败了(带有一条隐秘的消息)。

您可以使用.引用分组中的整个数据框，并且不需要匿名函数。您只需使用.变量直接调用(嵌套的)函数：

t1 <- mod1 %>% group_by(c) %>% do(lmer(m1.formula, .))

未测试，因为您没有提供可重现的示例。