ゼロ膨張モデルを使用してすべての因子レベルについてSEを推定する

Question

私はかなり複雑な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) 

このモデルの主な目的は、4つのレベル間で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 

私は論文や統計書籍やフォーラムを相談しましたが、私はまだ、この情報を提示するかどうかはわかりません。私が本当に欲しいのは、X軸のY軸と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) 

Predictions plot

移動するための方法をブートストラップされていますか？私はこれを試して、私が必要であるかどうかわからない何かのためにあらゆる種類のトラブルに遭遇する。

私が愚かな見落とし/仮定をしている場合は、事前にありがとうございます。私はまだ学んでいますが、これらの統計は私の手の届かないところから少し外れています。

Answer 1

最初に、モデル係数を信頼区間でプロットすることができます。 armパッケージにはcoefplotという機能がありますが、zeroinflモデルのメソッドはありませんので、ggplot2を使用して簡単な係数プロットを作成しました。 zeroinflモデルのpredictメソッドは、予測の信頼区間を提供しませんが、zeroinflモデルのブートストラップ信頼区間を構成する方法を示しています（this answer to a question on CrossValidated）。

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

次に、プロットの外観を示します。