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 の効果を 4 つのレベルにわたって調べることです。他の変数は、サンプリング プロセスの単なるアーティファクトです。

出力は次のとおりです。

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.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 に答える