2

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

数据在这里。

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

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

lmer每年的s 如下所示：

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

此外，如果我确实需要手动计算处理方法，有没有人知道如何计算我的实验中的嵌套水平？

我还想计算这些治疗手段中的每一个周围的误差线......

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

4

3 回答 3

2

我认为您要的是某种形式的predict()，其中没有类mer的默认方法lme4（至少是 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 日：

您可能会考虑bootMer()在开发版本的lme4（参数自举）预测置信区间中，它应该包含随机效应方差中的不确定性，并且可以与 GLMM 一起使用（例如，参见这个线程）。

这个想法是从感兴趣的模型进行模拟，用模拟值重新拟合，并从重新拟合的模型中计算感兴趣的统计量。您可以自己完成这些步骤，使用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)))
于 2013-08-25T20:51:23.523 回答
2

使用 package 中的函数的替代方法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()
于 2013-08-26T17:54:20.210 回答
0

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

要使用它，您只需执行以下操作：

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

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

于 2015-08-13T15:51:55.797 回答