从lmer对象中提取处理方法并计算错误条

Question

我正在详细介绍我为背景所做的实验--我很清楚lmer的方法，只是不清楚如何提取我需要的一些值/手工计算它们，因此我在SO上发布了这个，而不是简历。我希望这是正确的地方张贴！

数据在这里。

我的实验有一个平分图设计，有层次:块/图/子图。

有六个街区。每块有两个样地，每个块有两个子样地。处理1有两个层次(A和B)，并在地块一级应用:在每个区块有一个小区接受处理，1个A级和1个接受处理1级B。

处理2适用于子小区级，也有两个级别(C和D)：每个小区有一个子小区接收处理，两个A级和一个子小区接收处理2级B。

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

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

TMT1A:TMT2C

TMT1B:TMT2C

TMT1A:TMT2D

TMT1b:TMT2D

我用lmer作为我的模型来解释分裂图的设计。我运行的是跨年模型，但也是每年的模型(因为实验中的复制不允许在跨年模型中测试一年效应-这些模型最终被过度参数化)。

每年的lmer如下所示：

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

对于这些处理方法随时间的变化的图形表示，我希望提取每个处理的每个级别的处理方法(见上面的四个级别)，并为实验的每一年绘制这些处理方法，类似于本文中的例子。

我想知道，是否有可能从lmer对象中提取四种不同的治疗组合(如上面列出的那些)的治疗手段？还是要用手工计算？

我认为这样做的一种方法是创建另一个表示4种治疗组合的因素(参见粘贴数据中的"TMT1x2“栏)。然后，我可以为每年运行以下模型：

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

并提取治疗手段的每一个四个水平的方式。然而，我不确定这种方法是否能适当地控制分裂的地块设计，因为这个新的4层因子忽略了构成它的层次的嵌套性质(尽管随机效应并没有忽略它)。

此外，如果我真的需要手工计算处理方法，是否有人知道如何做到这一点，说明我的实验中筑巢的程度？

我也想在每种处理方法周围计算错误栏.

如果任何人有任何洞察力，这将是非常感谢！

Answer 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()

Answer 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年6月8日编辑：

您可能会在开发版本的bootMer()中考虑lme4用于(参数引导的)置信区间的预测，它应该包含随机效应的变化中的不确定性，并且将与GLMMs一起工作(例如，参见这条线)。

其思想是从兴趣模型中进行模拟，用模拟值进行修正，并从模型中计算利息统计量。您可以自己完成这些步骤，使用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)))

Answer 3

现在使用lme4，我认为bootMer()可能是最好的方法，因为它可以解释模型中的各种不确定性。但是，对于某些类型的问题，bootMer()无法工作，因为每个模型匹配可能需要多长时间。对于这些较大的问题，有一个称为merTools的R包，它提供了一种predictInterval方法来利用arm::sim来考虑固定效应和随机效应中的不确定性以及模型的残差。这是相当容易使用，并提供更快的预测在情况下，模型需要很长的时间来适应。对于随机效应之间的方差关系定义较好的问题，bootMer()给出的预测区间具有很好的覆盖范围。

要使用它，您只需：

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

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