从 lmer 对象中提取处理均值并计算误差线答案

【问题标题】：Extract treatment means from an lmer object and calculate error bars从 lmer 对象中提取处理均值并计算误差线
【发布时间】：2013-08-27 13:57:14
【问题描述】：

[我正在详细说明我的背景实验 - 我很清楚 lmers 的方法，只是不清楚如何提取我需要的一些值/手动计算它们，因此我在 SO 上发布了这个而不是简历。我希望这是发帖的正确位置！]

data are here。

我的实验采用裂区设计，级别为：块/图/子图。

有 6 个区块。每个区块有 2 个地块，每个地块有两个子地块。处理 1 有两个级别（A 和 B），并应用于地块级别：在每个区块中，有一个地块接受处理 1 级别 A，一个地块接受处理 1 级别 B。

处理 2 应用于子小区级别，也有两个级别（C 和 D）：每个小区有一个子小区接受处理 2 级别 A，一个子小区接受处理 2 级别 B。

实验进行了 3 年。我很感兴趣这两种治疗方法的每种组合如何影响我的因变量 (DV)。

因此，我有 4 种治疗组合：

TMT1A:TMT2C

TMT1B:TMT2C

TMT1A:TMT2D

TMT1b:TMT2D

我在我的模型中使用 lmer 来解释裂区设计。我正在运行一个跨年模型，但也依次运行一个模型（因为实验中的复制不允许在跨年模型中测试年份效应 - 模型最终被过度参数化）。

每年的lmers 如下所示：

m2011<- lmer (DV2011~ TMT1*TMT2 + (1|Block/TMT1))
m2012<- lmer (DV2012~ TMT1*TMT2 + (1|Block/TMT1))
m2013<- lmer (DV2013~ TMT1*TMT2 + (1|Block/TMT1))

对于这些处理均值随时间变化的图形表示，我想提取每年每个处理的每个级别（请参阅上面的四个级别）的处理均值，并为实验的每一年绘制这些图，类似于the example in this post

我想知道，是否可以从lmer 对象中提取四种不同治疗组合（如上面列出的那些）的治疗手段？还是必须手动计算？

我认为这样做的一种方法是实际创建另一个代表 4 种治疗组合的因子（请参阅粘贴数据中的“TMT1x2”列）。然后我可以每年运行以下模型：

m2011<- lmer (DV2011~ TMT1x2 + (1|Block/TMT1))

并以这种方式提取 4 个级别中的每个级别的处理方法。但是我不确定这种方法是否适合控制裂区设计，因为这个新的 4 水平因子忽略了构成它的水平的嵌套性质（尽管随机效应不会忽略它）...

此外，如果我确实需要手动计算处理均值，有谁知道如何考虑到我的实验中的嵌套级别？

我还想计算每个处理方法周围的误差线...

如果有人对此有任何见解，将不胜感激！

【问题讨论】：

您可能会发现 plotLMER.fnc 包中的 languageR 很有用。帮助页面上有一个示例“绘制两个因素之间的交互”。
谢谢亨利克。这对我的高斯模型很有用。您知道如何提取 mcmc 模拟提供的 HPD 值吗？我想在 ggplot 中绘制处理均值和误差值，因为它看起来更好:)。不幸的是，我也有一些带有二项式误差分布的模型。您知道可以为具有二项式误差的模型创建误差线的工具吗？ plotLMER.fnc 不能，因为它使用 mcmc sim 估计这些...
我认为我的答案太长了，所以我将其作为答案发布
嗨 Henrik - 谢谢 - 我期待着阅读它！

标签： r ggplot2 lme4 mixed-models lmer

【解决方案1】：

使用包languageR 中的函数的替代方案。我称你的数据集为df。

library(lme4)
library(languageR)
library(ggplot2)

# fit model
# n.b. I don't claim that this is a sensible model
# It is just used to demonstrate the plot
mod <- lmer(DV ~ TMT1 * TMT2 + (1|Block), data = df)

# create MCMC matrix
mcmc <- pvals.fnc(mod, nsim = 1000, withMCMC = TRUE)
# pval.fnc also calculates MCMC-based p-values and HPD confidence intervals,
# and plot the posterior distributions of the parameters

# plot using plotLMER.fnc 
# in addition, set withList = TRUE to create a list of data frames with plot data
# which can be used for a (possibly prettier) plot in ggplot
ll <- plotLMER.fnc(mod, withList = TRUE, pred = "TMT1", 
               intr = list(
                 "TMT2",
                 c("C", "D"),
                 "end",
                 list(c("red",  "blue"), rep(1, 2))),
               addlines = TRUE,
               mcmcMat = mcmc$mcmc)

 # here follows additional steps to plot using ggplot 

 # convert list to data frame
 df <- do.call(rbind, ll$TMT1)

 # rename 
 names(df)[names(df) == "Levels"] <- "TMT1"

 # add TMT2
 df$TMT2 <- rep(c("C", "D"), each = 2)

# plot using ggplot
dodge <- position_dodge(width = 0.1)
ggplot(data = df, aes(x = TMT1, y = Y, col = TMT2, group = TMT2)) +
   geom_point(position = dodge, size = 3) +
   geom_errorbar(aes(ymax = upper, ymin = lower, width = 0.1), position = dodge) +
   geom_line(position = dodge) +
   ylab("DV") +
   theme_classic()

【讨论】：

谢谢 Henrik - 这太棒了！你知道这个 HPD 值的 mcmc 模拟是否考虑了随机效应方差以及固定效应？另外，有没有一种方法可以与非高斯误差分布结合使用？
关于 mcmc 的详细信息，我认为您最好查阅 R 中相关函数的帮助文本。您还可以在作者的 this paper 中阅读有关 mcmcsamp 和 pvals.fnc 的更多信息languageR 和“lme4-Bates”。 Chapter 7 here 也是相关的。
关于来自 glmer 的估计 CI，我曾经在 arm 包中使用过 sim。不过，我从未将它用于不同因素组合的预测值。我尝试在 r-sig-mixed-models 上进行快速搜索，发现 this 和 an answer from "lme4-Bolker"。 “languageR 方法应该是这里的黄金标准”......“考虑随机效应参数的不确定性”。

【解决方案2】：

我认为您要的是某种形式的predict()，在lme4 中没有类mer 的默认方法（至少是CRAN 上的版本）。但是，您可以使用ez::ezPredict。

library(ez)
library(ggplot2)
to_predict <- expand.grid(TMT1=c("A","B"), TMT2=c("C","D"))
t_means <- rbind(ezPredict(m2011, to_predict=to_predict, boot=F), ezPredict(m2012, to_predict=to_predict, boot=F), ezPredict(m2013, to_predict=to_predict, boot=F) )
t_means$YEAR = rep(2011:2013, each = 4)
ggplot(t_means, aes(x=YEAR, y=value, color=TMT1:TMT2)) + geom_point() + geom_line()

此函数具有一些可能被证明有用的附加功能，例如提供引导值。

如果您想要的只是处理均值的点估计，手动计算同样容易，尤其是所有三个模型都具有相同的设计矩阵：

mm = unique(model.matrix(m2011))
Y_bar <- c(mm%*%fixef(m2011), mm%*%fixef(m2012), mm%*%fixef(m2013))
ggplot(t_means, aes(x=YEAR, y=Y_bar, color=TMT1:TMT2)) + geom_point() + geom_line()

我不太确定您所说的 “计算处理方法……考虑到我的实验中的嵌套级别”是什么意思。混合模型中的随机效应是结构化的、与总体水平效应（固定效应）的正态分布偏差。查看随机效应估计值ranef(m2011) 和相关的设计矩阵m2011@Zt 可能是有益的。

因此，如果您只想绘制总体水平的处理均值，您可以简单地使用上述固定效应fixef(m2011) 和固定效应设计矩阵model.matrix(m2011) 的估计值。如果您想在总体水平预测中包含一些不确定性度量，或者想要对每个区块/地块/子地块进行预测，则需要同时使用随机效应和固定效应。我建议您首先查看标题“预测的预测和/或置信度（或预测）区间”下的http://glmm.wikidot.com/faq。

2013 年 8 月 26 日编辑：

您可以在lme4 的开发版本中考虑bootMer() 用于预测周围的（参数自举）置信区间，它应该包含随机效应方差中的不确定性，并且适用于 GLMM（例如参见this thread )。

这个想法是从感兴趣的模型进行模拟，用模拟值重新拟合，并从重新拟合的模型中计算感兴趣的统计量。您可以自己完成这些步骤，simulate() 和 refit()：

t_sim <- apply(simulate(model, 999), 2, function(x) combn(unique(model.matrix(model))%*%fixef(refit(model, x)), 2, diff) )

它会生成 999 个处理均值之间成对差异的引导表示，您可以在上面使用 quantile()（或您希望的任何引导置信区间）：

apply(t_sim, 1, function(.) quantile(., c(0.975, 0.025)))

【讨论】：

谢谢内特。阅读该链接，他们似乎无法拟合考虑随机效应不确定性的误差线。使用 plotLMER.fnc 我可以使用 mcmc 模拟高斯模型来估计误差，但是我有一些二项式模型，这不适用于...这是一种耻辱，因为我想测试 4 tmt 平均值之间的差异，我不知道怎么做。我的模型只告诉我每种治疗的主要效果以及它们之间的相互作用。他们不测试 4 种治疗水平是否存在差异。有没有办法在 lme4 中设置对比度来实现这一点？
@Sarah 你可能会考虑参数引导。在lme4的开发版中有这个程序的功能，也可以用simulate()和refit()自己做。我在上面的答案中添加了后者的示例。

【解决方案3】：

现在有了lme4，我认为bootMer() 可能是最好的选择，因为它考虑了模型中的各种不确定性。但是，对于某些类别的问题，bootMer() 无法工作，因为每个模型拟合可能需要多长时间。对于这些较大的问题，有一个名为 merTools 的 R 包，它提供了一个 predictInterval 方法来使用 arm::sim 来解释固定和随机效应的不确定性以及模型的残差。在模型需要很长时间才能拟合的情况下，它相当容易使用并且提供预测的速度要快得多。它很好地覆盖了bootMer() 产生的预测区间，用于解决随机效应之间关系的方差相当明确的问题。

要使用它，您只需：

library(merTools)
preds <- predictInterval(m2011, newdata = myData, level = 0.95, n.sims = 1000)

还有其他几个用户可配置的选项，但结果是一个预测对象，类似于从 lm 请求预测间隔时产生的预测对象——一个三列 data.frame，列有 fit、lwr 和 @ 987654332@。

【讨论】：