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

这确实奏效，但我真的更希望将多年平均，而不是将年份设置为一个特定水平。但是，每当我尝试平均年份时，都会收到以下错误消息：

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轴上具有预测的寄生虫负担的图表，以证明每种基因型在寄生虫负担方面的预测差异，同时考虑了所有其他重要的协变量。

Answer 1

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

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

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

现在我们使用 merTools 中的函数准备数据draw。在这里，我们从模型框架中得出平均观察值。然后，我们wiggle通过扩展数据框来包含重复的相同观察结果，但参数指定的变量值不同var。在这里，我们将数据集扩展到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")