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

【问题标题】：running lmer with a by/group by statement?使用 by/group by 语句运行 lmer？
【发布时间】：2015-01-19 10:45:05
【问题描述】：

我正在尝试找到一种快速运行 lmer 模型的方法，但要为每个分组变量单独运行它（在 SAS 中，可以使用 by= 语句）。我已经尝试使用 dplyr 和我找到的代码：

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

但这似乎不起作用。

有人知道如何使用 dplyr 或其他方法做到这一点吗？

【问题讨论】：

您为什么要这样做？通过像这样对数据进行分组，您忽略数据的结构，而混合模型是为了在您的模型中包含此结构。
请重现示例...？
大家好，感谢您的回答。我想要的是为每个组（“站”这里）运行单独的混合模型。 Ben：我不确定你所说的可重现是什么意思（对不起，我是一个新人）。你的意思是输出还是DF本身？这是我得到的错误：> t1% group_by(stn) %>% do(function(df){lmer(m1.formula,data=df)}) 错误：不知道怎么做使用最新的 dplyr 1.4 thx Z 处理类型对列表 im
请参阅here，了解如何提供可重现的示例。它允许我们尝试一些事情以提供帮助。

标签： r lme4

【解决方案1】：

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

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

统计：如上所述，考虑在数据集中的每个层（分组变量）上单独运行混合模型有点奇怪。通常混合模型的全部意义在于将模型拟合到组合数据集，尽管我当然可以想象例外情况（下面，我按性别拟合单独的混合模型）。您可能正在寻找类似 @987654322@ 函数（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)))

给出合理的结果。

【讨论】：

我目前无法对其进行测试，但我认为这将在 dplyr 中工作：Orthodont %>% group_by(Sex) %>% do(model = lmer(.,formula=distance~age+(1|Subject))) 尽管您稍后需要提取列表列中的条目（我所做的只是包括 @987654346 @里面do())
@docendodiscimus：它确实有效，但它产生了一个非常不寻常（或至少对我来说非常熟悉）的对象......

【解决方案2】：

问题是您错误地调用了do() - 它不适用于这样的匿名函数。 do() 中的参数是在数据的上下文中评估的，所以当你说function(df) 时，do 将尝试使用数据的df 列。它没有那一列，所以它失败了（带有一个神秘的消息）。

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

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

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

【讨论】：