glmer - 使用二项式数据进行预测（cbind 计数数据）答案

【问题标题】：glmer - predict with binomial data (cbind count data)glmer - 使用二项式数据进行预测（cbind 计数数据）
【发布时间】：2014-02-21 21:21:15
【问题描述】：

我正在尝试预测在我的二项式数据上运行的 glmer 模型随时间变化的值（x 轴上的天数）。 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），我们已经考虑了过度分散。

然后我们使用疏浚命令来确定具有主要影响（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"

【问题讨论】：

标签： r glm predict lme4

【解决方案1】：

您的大部分问题是您使用的是 lme4 的 1.0 之前的版本，它没有实现 predict 方法。（更新将是最简单的，但我相信如果你因为某种原因不能，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))

【讨论】：