在我的混合模型中使用lme4预测函数时遇到了问题

Question

我在使用lme4预测函数时遇到了一些困难。当进行预测时，我想要设置一些解释变量为指定水平，但是要在其他变量上取平均值。

这里有一些虚构的数据，是我原始数据集的简化版本：

a <-  data.frame(
    TLR4=factor(rep(1:3, each=4, times=4)), 
    repro.state=factor(rep(c("a","j"),each=6,times=8)), 
    month=factor(rep(1:2,each=8,times=6)), 
    sex=factor(rep(1:2, each=4, times=12)), 
    year=factor(rep(1:3, each =32)), 
    mwalkeri=(sample(0:15, 96, replace=TRUE)), 
    AvM=(seq(1:96))
)

AvM编号是水田鼠的识别编号。响应变量（mwalkeri）是每只田鼠身上跳蚤数量的计数。我感兴趣的主要解释变量是Tlr4，这是一个拥有3个不同基因型（编码为1、2和3）的基因。其他解释变量包括生殖状态（成年或幼年）、月份（1或2）、性别（1或2）和年份（1、2或3）。我的模型如下所示（当然，该模型对于虚构数据现在已经不合适，但这并不重要）：

install.packages("lme4")
library(lme4)
mm <- glmer(mwalkeri~TLR4+repro.state+month+sex+year+(1|AvM), data=a, 
    family=poisson,control=glmerControl(optimizer="bobyqa"))`
summary(mm)

我希望可以预测每个不同Tlr4基因型的寄生虫负担，同时考虑到所有其他协变量。为了达到这个目的，我创建了一个新的数据集来指定我想要设置每个解释变量的水平，并使用预测函数：

b <-  data.frame(
    TLR4=factor(1:3), 
    repro.state=factor(c("a","a","a")),
    month=factor(rep(1, times=3)), 
    sex=factor(rep(1, times=3)), 
    year=factor(rep(1, times=3))
)
predict(mm, newdata=b, re.form=NA, type="response")

这种方法虽然可行，但我更希望能够对多年份进行平均，而不是将年份设定为一个特定的级别。然而，每当我尝试对年份进行平均时，就会出现以下错误信息：

Error in model.frame.default(delete.response(Terms), newdata, na.action = na.action, : factor year has new level

我是否可以跨年份进行平均，而不是选择指定的级别？此外，我还没有找到如何获取这些预测的标准误差的方法。目前我唯一能够获取预测标准误差的方式是使用 lsmeans() 函数（来自 lsmeans 包）：

c <- lsmeans(mm, "TLR4", type="response")
summary(c, type="response")

这会自动生成标准误差，但是它是通过对所有其他解释变量进行平均而生成的。我相信可能可以更改这一点，但如果可以的话，我宁愿使用 predict() 函数。我的目标是创建一个图表，其中Tlr4 基因型在 x 轴上，预测寄生虫负荷在 y 轴上，以演示在考虑所有其他显著协变量的情况下，每个基因型的预测寄生虫负荷的差异。

- Martha

---

你还没有包含一个可重现的示例，所以很难提出具体的编码建议并测试它们是否有效。只要你传入一个新的数据框，并且列名与原始数据集完全匹配，那么使用predict()应该非常容易。你可以尝试在代码中包含一些尝试的内容和精确的错误信息。 - MrFlick

我想花些时间在这上面，但还没有机会。对于每个固定效应，您必须指定一些值，但以下三个选项中的任何一个都是合理的：（1）将其设置为平均值或基线值（“条件”）；（2）允许它在合理范围内变化，然后平均预测值（“边际”）；（3）允许它在一定范围内变化，并找到一种方法在最终图中绘制所有值（分组线条、在分面内、通过颜色或其他方式）。 - Ben Bolker

谢谢您的回复！@MrFlick我编辑了我的原始问题，包含一些我的代码，并提供了一个可复现的例子；希望这能让它更清晰一些。@Ben Bolker是否可以将我一些说明变量设置为特定的水平，但根据生物学上最合理的情况来查看其他变量的平均值？如果我将其设置为数据集中已编码的水平，那么当我使用predict()时我没问题，但是如果我尝试做其他事情（例如取年份的平均值），我就无法使其正常工作。一如既往，任何建议都将不胜感激！ - Martha

你确定在这种情况下要将年/月视为分类变量吗？将年份作为因素进行折叠就像折叠性别一样。没有“平均”性别。我猜你可以分别拟合每个，然后平均预测结果，但这真的会影响标准误差计算。 - MrFlick

此外，在?predict.merMod帮助页面中，他们有一条注释说：“没有计算预测标准误差的选项，因为很难定义一种有效的方法来包含方差参数的不确定性；我们建议使用bootMer来完成此任务。”因此，您可能需要考虑引导法来估计误差（除非您有其他计算方法）。 - MrFlick

Answer 1

你可能会对merTools软件包感兴趣，其中包括一些用于创建反事实数据集的函数，然后在该新数据上进行预测以探索变量对结果的实质影响。这个README和软件包vignette提供了一个很好的例子：

假设我们想要探索具有类别和连续预测变量交互项模型的影响。首先，我们拟合一个带有交互作用的模型：

data(VerbAgg)
fmVA <- glmer(r2 ~ (Anger + Gender + btype + situ)^2 +
       (1|id) + (1|item), family = binomial, 
       data = VerbAgg)

现在我们使用merTools中的draw函数准备数据。在这里，我们从模型框架中绘制平均观测值。然后，我们通过扩展数据框来包含重复的相同观测值，但是具有由var参数指定的不同变量值来wiggle数据。在这里，我们将数据集扩展到所有btype、situ和Anger的值。

# Select the average case
newData <- draw(fmVA, type = "average")
newData <- wiggle(newData, var = "btype", values = unique(VerbAgg$btype))
newData <- wiggle(newData, var = "situ", values = unique(VerbAgg$situ))
newData <- wiggle(newData, var = "Anger", values = unique(VerbAgg$Anger))

head(newData, 10)

#>    r2 Anger Gender btype  situ id        item
#> 1   N    20      F curse other  5 S3WantCurse
#> 2   N    20      F scold other  5 S3WantCurse
#> 3   N    20      F shout other  5 S3WantCurse
#> 4   N    20      F curse  self  5 S3WantCurse
#> 5   N    20      F scold  self  5 S3WantCurse
#> 6   N    20      F shout  self  5 S3WantCurse
#> 7   N    11      F curse other  5 S3WantCurse
#> 8   N    11      F scold other  5 S3WantCurse
#> 9   N    11      F shout other  5 S3WantCurse
#> 10  N    11      F curse  self  5 S3WantCurse

现在我们只需将这个新数据集传递给predictInterval，以生成这些反事实情况的预测。然后我们将预测值绘制在连续变量Anger上，并在两个分类变量situ和btype上分别进行分面和分组。

plotdf <- predictInterval(fmVA, newdata = newData, type = "probability", 
        stat = "median", n.sims = 1000)
plotdf <- cbind(plotdf, newData)
ggplot(plotdf, aes(y = fit, x = Anger, color = btype, group = btype)) + 
  geom_point() + geom_smooth(aes(color = btype), method = "lm") + 
  facet_wrap(~situ) + theme_bw() +
  labs(y = "Predicted Probability")

- jknowles