r - 用by/group by语句运行lmer？

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

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

但这似乎不起作用。

有人知道如何使用dplyr或其他方法执行此操作吗？

最佳答案

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

您可能需要在此处考虑两个基本问题:

statistics:如上所述，考虑在数据集中的每个层次(分组变量)上分别运行混合模型有点​​奇怪。通常，混合模型的全部目的是使模型适合组合的数据集，尽管我当然可以想象到例外(下面，我按性别拟合了单独的混合模型)。您可能正在寻找类似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)))

给出合理的结果。