2014-01-29 75 views
5

我試圖預測在我的二項式數據上運行的glmer模型隨時間推移的值(x軸的天數)。 Total Alive和Total Dead是計數數據。這是我的模型,以及下面的相應步驟。glmer - 用二項式數據預測(cbind計數數據)

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:指數)。

然後,我們使用疏浚命令來確定主效應(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軸上的總體積和X軸上的天數,以及取決於模型輸出的擬合線。我理解,這是棘手,因爲該模型級聯的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" 

回答

12

您問題在於你使用的lme4之前的版本爲lme4,它沒有實現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 REForm=NA`1.0.5) :

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)) 

enter image description here