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

我在尝试在混合模型上使用 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")

这确实有效，但我真的更喜欢对年份进行平均，而不是将年份设置为一个特定水平。但是，每当我尝试平均年份时，我都会收到此错误消息：

model.frame.default(delete.response(Terms), newdata, na.action = na.action, 中的错误：因子年份有新级别

我是否可以对年份进行平均而不是选择指定的级别？另外，我还没有弄清楚如何获得与这些预测相关的标准误差。我能够获得预测标准误差的唯一方法是使用lsmeans()函数（来自 lsmeans 包）：

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

它会自动生成标准错误。然而，这是通过对所有其他解释变量进行平均而生成的。我确信可能可以改变这一点，但我宁愿使用predict()如果可以的话。我的目标是创建一个图表，其中 x 轴为 Tlr4 基因型，y 轴为预测的寄生虫负担，以证明每种基因型的寄生虫负担的预测差异，同时考虑所有其他重要的协变体。

您可能感兴趣merTools包，其中包括几个函数，用于创建反事实数据集，然后对新数据进行预测，以探索变量对结果的实质性影响。自述文件和包 vignette 就是一个很好的例子：

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

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

现在我们使用以下方法准备数据drawmerTools 中的函数。在这里，我们从模型框架中得出平均观察值。然后我们wiggle通过扩展数据框以包含重复但具有由指定变量的不同值的相同观察来获取数据var范围。在这里，我们将数据集扩展到所有值btype, situ, and 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 and 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")