0

我有一个相当复杂的 ZINB 模型。我试图复制我正在尝试做的基本结构：

MyDat<-cbind.data.frame(fac1 = rep(c("A","B","C","D"),10), 
       fac2=c(rep("X",20),rep("Y",20)), 
       offset=c(runif(20, 50,60),runif(20,150,165)), 
       fac3=rep(c(rep("a1",4),rep("a2",4),rep("a3",4),rep("a4",4),rep("a5",4)),2),
       Y=c(0,0,0,1,0,0,11,10,0,0,0,5,0,0,0,35,60,0,0,0,0,2,0,0,16,0,0,0,0,0,3,88,0,0,0,0,0,0,27,0))

f<-formula(Y~fac1+ offset(log(offset))|fac3+ fac2) 
ZINB <-zeroinfl(f, dist = "negbin",link = "logit", data = MyDat)
summary(ZINB)

该模型的主要目标是查看 fac1 在四个级别上的影响。其他变量只是采样过程的产物。

这是输出：

Call:
zeroinfl(formula = f, data = MyDat, dist = "negbin", link = "logit")

Pearson residuals:
      Min        1Q    Median        3Q       Max 
-0.418748 -0.338875 -0.265109 -0.001566  2.682920 

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.7192     0.9220  -1.865 0.062239 .  
fac1B        -4.4161     1.4700  -3.004 0.002663 ** 
fac1C        -1.2008     1.2896  -0.931 0.351778    
fac1D         0.1928     1.3003   0.148 0.882157    
Log(theta)   -1.7349     0.4558  -3.806 0.000141 ***

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -11.5899   210.8434  -0.055    0.956
fac3a2       -0.4775     2.4608  -0.194    0.846
fac3a3      -11.2284   427.5200  -0.026    0.979
fac3a4       10.7771   210.8056   0.051    0.959
fac3a5       -0.3135     2.3358  -0.134    0.893
fac2Y        11.8292   210.8298   0.056    0.955
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 0.1764 
Number of iterations in BFGS optimization: 76 
Log-likelihood: -63.82 on 11 Df

我已经查阅了论文和统计书籍和论坛，但我仍然不确定如何呈现这些信息。我真正想要的是一个条形图，显示 Y 轴上的效果和 X 上的 4 个级别。

如果我理解正确，fac1 的 A 级当前设置为 0，是我的参考级别（如果我在这里错了，请纠正我）。因此，我可以绘制 4 个级别（包括 A 级为零）的图。这似乎并不理想。我真的很想在所有级别都有 95% 的 CI。

我也可以使用 predict 函数，但是 predict.zeroinfl 没有给出错误估计，我不确定如何解释偏移量的影响。

类似的论文只是将原始数据的箱线图放在预测的箱线图旁边并进行比较。我觉得我应该可以做得更好。

下面是创建预测值的代码和绘图：

MyDat$phat<-predict(ZINB, type="response")
MyDat$phat_os<-MyDat$phat/MyDat$offset

plot(phat~fac1, MyDat)

预测图

引导是要走的路吗？我已经尝试过这个并且遇到了各种各样的麻烦，因为我不确定是否有必要。

提前谢谢你，如果我在做一个愚蠢的疏忽/假设，请放轻松。我还在学习，但这些统计数据感觉有点遥不可及。

4

1 回答 1

1

对于初学者，您可以绘制模型系数及其置信区间。该arm包具有该coefplot功能，但它没有zeroinfl模型的方法，因此我在下面使用创建了一个简单的系数图ggplot2。模型的predict方法zeroinfl不提供预测的置信区间，但这个对 CrossValidated 问题的回答显示了如何为zeroinfl模型构建自举置信区间。

关于fac1:的级别A是参考级别，因此其他级别的系数与 相关fac1 = "A"。

library(pscl)
library(ggplot2)

MyDat<-cbind.data.frame(fac1 = rep(c("A","B","C","D"),10), 
                       fac2=c(rep("X",20),rep("Y",20)), 
                       offset=c(runif(20, 50,60),runif(20,150,165)), 
                       fac3=rep(c(rep("a1",4),rep("a2",4),rep("a3",4),rep("a4",4),rep("a5",4)),2),
                       Y=c(0,0,0,1,0,0,11,10,0,0,0,5,0,0,0,35,60,0,0,0,0,2,0,0,16,0,0,0,0,0,3,88,0,0,0,0,0,0,27,0))

f<-formula(Y ~ fac1 + offset(log(offset))|fac3 + fac2) 
ZINB <-zeroinfl(f, dist = "negbin",link = "logit", data = MyDat)

# Extract coefficients and standard errors from model summary
coefs = as.data.frame(summary(ZINB)$coefficients$count[,1:2])
names(coefs)[2] = "se" 
coefs$vars = rownames(coefs)

# Coefficient plot
ggplot(coefs, aes(vars, Estimate)) + 
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  geom_errorbar(aes(ymin=Estimate - 1.96*se, ymax=Estimate + 1.96*se), 
                lwd=1, colour="red", width=0) +
  geom_errorbar(aes(ymin=Estimate - se, ymax=Estimate + se), 
                lwd=2.5, colour="blue", width=0) +
  geom_point(size=4, pch=21, fill="yellow") +
  theme_bw()

这就是情节的样子。

在此处输入图像描述

于 2016-04-13T20:55:25.807 回答