GLMER - 用于二项数据（cbind计数数据）的预测

Question

我正在尝试预测随时间变化（x轴上的天数）我的二项数据上运行的glmer模型的值。Total Alive和Total Dead是计数数据。这是我的模型以及相应的步骤：

full.model.dredge<-glmer(cbind(Total.Alive,Total.Dead)~(CO2.Treatment+Lime.Treatment+Day)^3+(Day|Container)+(1|index),
                         data=Survival.data,family="binomial")

我们已经在代码中解决了过度离散问题 (1:index)。

然后，我们使用 dredge 命令确定最佳拟合模型及其主效应 (CO2.Treatment、Lime.Treatment、Day) 和相应的交互作用。

dredge.models<-dredge(full.model.dredge,trace=FALSE,rank="AICc")

然后为它们创建了一个工作区变量。

my.dredge.models<-get.models(dredge.models)

我们随后进行了模型平均，以平均最佳拟合模型的系数。

silly<-model.avg(my.dredge.models,subset=delta<10)

但是现在我想创建一个图表，Y轴上是Total Alive，X轴上是Days，还有一条根据模型输出的拟合线。我知道这很棘手，因为该模型连接了Total.Alive和Total.Dead（请参见模型中的cbind(Total.Alive,Total.Dead)）。

当我尝试运行预测命令时，会出现错误。

# 9: In UseMethod("predict") :
#   no applicable method for 'predict' applied to an object of class "mer"

- user3250522

Answer 1

你的主要问题是使用了一个版本低于1.0的lme4，该版本没有实现方法。(更新最简单，但我相信如果由于某些原因你不能更新，那么在http://glmm.wikidot.com/faq有一种方法可以手动进行预测，通过提取固定效应设计矩阵和系数...)实际上，预测结果并没有问题，默认情况下预测对数几率或概率(type="response")；如果你想预测数字，你需要适当地乘以N。

虽然你没有给出一个示例，但这里有一个可重现(尽管有点琐碎)的示例，使用内置的cbpp数据集(我确实会得到一些警告信息--no non-missing arguments to max; returning -Inf--但我认为这可能是因为模型中只有一个非平凡固定效应参数?)

library(lme4)
packageVersion("lme4")  ## 1.1.4, but this should work as long as >1.0.0
library(MuMIn)

方便后续使用(ggplot)，添加一个比例变量：

cbpp <- transform(cbpp,prop=incidence/size)

拟合模型（您也可以使用glmer(prop~..., weights=size, ...)）

gm0 <- glmer(cbind(incidence, size - incidence) ~ period+(1|herd),
           family = binomial, data = cbpp)
dredge.models<-dredge(gm0,trace=FALSE,rank="AICc")
my.dredge.models<-get.models(dredge.models)
silly<-model.avg(my.dredge.models,subset=delta<10)

预测确实有效：

predict(silly,type="response")

创建绘图：

library(ggplot2)
theme_set(theme_bw())  ## cosmetic
g0 <- ggplot(cbpp,aes(period,prop))+
    geom_point(alpha=0.5,aes(size=size))

建立一个预测框架：

predframe <- data.frame(period=levels(cbpp$period))

在人口层面上进行预测 (ReForm=NA -- 在lme4 `1.0.5中可能需要改为 REForm=NA)：

predframe$prop <- predict(gm0,newdata=predframe,type="response",ReForm=NA)

将其添加到图表中：

g0 + geom_point(data=predframe,colour="red")+
    geom_line(data=predframe,colour="red",aes(group=1))

- Ben Bolker