4.A 過大分散カウントデータのモデル
Ph.D. Biochemistsによって作られた出版物の数に関するLong (1990) のデータを使って、ポアソンモデル、過剰分散ポアソン、負の二項モデル、ゼロインフレートポアソンモデルの適用を説明する。
データセットの変数は
| Variable | Description |
|---|---|
art |
博士課程の過去3年間の論文である。 |
fem |
女性の場合は1つコーディング |
mar |
既婚の場合は1つコーディング |
kid5 |
6歳未満の子供の数 |
phd |
博士号の威光 |
kid5 |
博士号の威光 |
| 6歳未満の子供の数4487D. program | |
ment |
過去3年間のメンターによる論文 |
これらのデータは Long and Freese (2001) でも分析されており、Stata ウェブサイトから入手できます：
> library(foreign)> ab names(ab) "art" "fem" "mar" "kid5" "phd" "ment"> r c(mean=r, var=r, ratio=r/r) mean var ratio 1.692896 3.709742 2.191358
論文数の平均値は 1.69, 分散値は 3.71 で平均値の少し2倍の数値になっていることがわかりました。 データは過分散ですが、もちろん共変量はまだ考慮していません。
A Poisson Model
ここで、Long and Freese(2001)が用いたモデル、5つの予測変数すべてを使った単純な加法モデルを当てはめてみましょう。
> mp summary(mp)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = poisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.304617 0.102981 2.958 0.0031 ** femWomen -0.224594 0.054613 -4.112 3.92e-05 ***marMarried 0.155243 0.061374 2.529 0.0114 * kid5 -0.184883 0.040127 -4.607 4.08e-06 ***phd 0.012823 0.026397 0.486 0.6271 ment 0.025543 0.002006 12.733We see that the model obviously doesn't fit the data. The five-percent critical value for a chi-squared with 909 d.f. is
>qchisq(0.95, df.residual(mp)) 980.2518> deviance(mp) 1634.371> pr sum(pr^2) 1662.547で、devianceとPearsonのカイ二乗はともに1600台である。
4.A.1.のようになる。
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524ここで分散が平均に等しいのではなく、比例していると仮定し、ピアソンのカイ二乗をそのd.f.で割ったスケールパラメータφを推定してみる。
Rは
quasipoissonfamily:> mq summary(mq)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = quasipoisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.304617 0.139273 2.187 0.028983 * femWomen -0.224594 0.073860 -3.041 0.002427 ** marMarried 0.155243 0.083003 1.870 0.061759 . kid5 -0.184883 0.054268 -3.407 0.000686 ***phd 0.012823 0.035700 0.359 0.719544 ment 0.025543 0.002713 9.415The estimates are exactly the same as before, but the standard errorsare about 35% larger. We can verify this fact easily. First we writea useful function to extract standard errors and then use it onour fits:
> se round(data.frame(p=coef(mp), q=coef(mq),+ se.p=se(mp), se.q=se(mq), ratio=se(mq)/se(mp)), 4) p q se.p se.q ratio(Intercept) 0.3046 0.3046 0.1030 0.1393 1.3524femWomen -0.2246 -0.2246 0.0546 0.0739 1.3524marMarried 0.1552 0.1552 0.0614 0.0830 1.3524kid5 -0.1849 -0.1849 0.0401 0.0543 1.3524phd 0.0128 0.0128 0.0264 0.0357 1.3524ment 0.0255 0.0255 0.0020 0.0027 1.3524を使用すると、この計算をしてくれます。別のアプローチとしては、ポアソンモデルをあてはめ、標準誤差のロバスト推定量またはサンドイッチ推定量を使用することです。
4.A.2 Negative Binomial Regression
ここで同じ予測変数で負の二項モデルを当てはめてみる。
> library(MASS)> mnb summary(mnb)Call:glm.nb(formula = art ~ fem + mar + kid5 + phd + ment, data = ab, init.theta = 2.264387695, link = log)Deviance Residuals: Min 1Q Median 3Q Max -2.1678 -1.3617 -0.2806 0.4476 3.4524 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.256144 0.137348 1.865 0.062191 . femWomen -0.216418 0.072636 -2.979 0.002887 ** marMarried 0.150489 0.082097 1.833 0.066791 . kid5 -0.176415 0.052813 -3.340 0.000837 ***phd 0.015271 0.035873 0.426 0.670326 ment 0.029082 0.003214 9.048 1/mnb$theta 0.4416205> -2*(logLik(mp)-logLik(mnb))'log Lik.' 180.196 (df=6)Rの
thetaは乗法的なランダム効果の精度であり、表中の1/σ2に相当する。このパラメータの有意性を検定するには、このモデルとポアソンモデルの対数尤度の差の2倍、180.2を計算し、それを1dfのカイ2乗として扱うことを考えればよいでしょう。 しかし、帰無仮説がパラメータ空間の境界上にあるので、通常の漸近法は適用されない。
より良い近似は、統計量を0と1つのd.f.を持つカイ2乗の50：50混合として扱うことを示すいくつかの仕事がある。
回帰係数に関する仮説の検定には、Wald検定または尤度比検定のいずれかを使用することができ、これは完全な分布の仮定をしているので可能です。
glmにはnegative.binomial族もあり、これはパラメータthetaが与えられていれば使用することができます。 (これは負の二項が固定分散のglm族であるという結果に基づいています)。> mnbg all(abs(coef(mnb)==coef(mnbg))The standard errors would differ, however, because allows for the fact that was estimated, whereas does not.
Unobserved Heterogeneity
R has a function to compute the density of a gamma distribution with given shape and scale (or its reciprocal the rate). In particular, when the random effect hasvariance the density is. This can beused to draw the density.
> v = 1/mnb$theta> gd = data.frame(x = seq(0.0, 3, 0.1), g = dgamma(x, shape = 1/v, scale = v))> ggplot(gd, aes(x, g)) + geom_line()> ggsave("gamdenr.png", width=500/72, height=400/72, dpi=72)
また、
qgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE)を使っても分位数が計算できます。 特に、ランダム効果が分散vの場合、四分位はqgamma((1:3)/4, shap e= 1/v, scale = v)となる。> qgamma((1:3)/4, shape = 1/v, scale = v) 0.5114167 0.8572697 1.3347651観察されない異質性の分布のQ1にいる生化学者の論文出版数は、彼らの観察特性からの予想より49%少なく、中央値の人は14%少なく、Q3の人は予想より33%多く出版している。
推定値と標準誤差の比較
ポアソンモデル、過分散ポアソンモデル、負の二項モデルでのパラメータ推定値と標準誤差を比較してみましょう：
> round(data.frame(+ p=coef(mp),q=coef(mq),nb=coef(mnb),+ se.p=se(mp),se.q=se(mq),se.nb=se(mnb)),4) p q nb se.p se.q se.nb(Intercept) 0.3046 0.3046 0.2561 0.1030 0.1393 0.1373femWomen -0.2246 -0.2246 -0.2164 0.0546 0.0739 0.0726marMarried 0.1552 0.1552 0.1505 0.0614 0.0830 0.0821kid5 -0.1849 -0.1849 -0.1764 0.0401 0.0543 0.0528phd 0.0128 0.0128 0.0153 0.0264 0.0357 0.0359ment 0.0255 0.0255 0.0291 0.0020 0.0027 0.0032負の二項推定値はポアソンモデルに基づくものと大きな違いはなく、両方のセットが同じ結論を導くことになります。
標準誤差を見ると、過大分散に対するどちらのアプローチも非常によく似た標準誤差を推定しており、通常のポアソン回帰は標準誤差を過小評価していることがわかる。
適合度
負の2項モデルの適合度は、デビアンス
> deviance(mnbg) 1004.281を用いて評価することができる。
分散関数
過分散のポワソンモデルと負の2項モデルは、異なる分散関数を持っている。
ここに負の2項線形予測変数に基づくグループがあり、5(5)95%で区切り、ほぼ同じ大きさの20のグループを作るために
cutを使って作成されています。predict()は、特に指定されない限り、線形予測器を計算することに注意してください。 応答のスケールで予測するには、オプションtype="response"を追加してください。> xb g m v png("c4afig1.png",width=500,height=400)> plot(m, v, xlab="Mean", ylab="Variance", + main="Mean-Variance Relationship")> mtext("Articles Published by Ph.D. Biochemists",padj=-0.5)> x lines(x, x*phi, lty="dashed")> lines(x, x*(1+x/mnb$theta))> legend("topleft", lty=c("dashed","solid"), + legend=c("Q. Poisson","Neg. Binom."), inset=0.05)> dev.off();null device 1
グラフは平均と分散をプロットし、分散がφμである過分散ポアソンモデルと分散がμ（1+μσ2）である負の二項モデルに対応するカーブを重ねて表示しています。 負の二項分散関数は、あまり違いはないが、二次関数であるため、より速く上昇することができ、ハイエンドでより良い仕事をする。
4.A.3 ゼロ膨張ポアソンモデル
回数データでよく見られるのは、ポアソンモデルで予想されるものよりゼロが多くなることである。 これは実際に我々のデータで問題になっている。
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071サンプル中の科学者の30.0%が博士号取得後の3年間に論文を発表していないことがわかりますが、ポアソンモデルは20.9%だけが論文なしであることを予測します。
このような状況をモデル化する1つの方法は、データが2つの集団の混合から来たと仮定することです。1つはカウントが常にゼロであり、もう1つはカウントが平均μのポアソン分布を持っているところです。
博士号を持つ生化学者の出版物の文脈では、出版物が重要でない仕事を考えていた人もいれば、出版物の記録が期待される学術的な仕事を目指していた人もいたと想像できます。
そして結果の分布は2つのパラメータ、πは「常にゼロ」の確率、μは「常にゼロ」グループでない人の出版数の平均でモデル化することができる。 共変量を導入する自然な方法は、常にゼロの確率πと常にゼロクラスでないものの平均μのlogをモデル化することです。
このタイプのモデルは、
psclパッケージのzeroinfl()関数を使用してRでフィットすることができます。 同じ変数が両方の式に含まれる場合、モデル式は通常通り指定することができる。ここに、両方の方程式にすべての共変数を含むゼロインフレートモデルがあります。
> library(pscl)> mzip summary(mzip)Call:zeroinfl(formula = art ~ fem + mar + kid5 + phd + ment, data = ab)Pearson residuals: Min 1Q Median 3Q Max -2.3253 -0.8652 -0.2826 0.5404 7.2976 Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.640838 0.121307 5.283 1.27e-07 ***femWomen -0.209145 0.063405 -3.299 0.000972 ***marMarried 0.103751 0.071111 1.459 0.144565 kid5 -0.143320 0.047429 -3.022 0.002513 ** phd -0.006166 0.031008 -0.199 0.842378 ment 0.018098 0.002294 7.888 3.07e-15 ***Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.577059 0.509386 -1.133 0.25728 femWomen 0.109746 0.280082 0.392 0.69518 marMarried -0.354014 0.317611 -1.115 0.26502 kid5 0.217097 0.196482 1.105 0.26919 phd 0.001274 0.145263 0.009 0.99300 ment -0.134114 0.045243 -2.964 0.00303 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 21 Log-likelihood: -1605 on 12 Dfインフレート方程式を見ると、「常にゼロ」クラスであることの唯一の有意な予測因子は、メンターによって出版された記事の数であり、メンターによる各記事は出版しないことの12.6%低い確率と関連していることが分かります。
常にゼロのクラスでないものの平均論文数または論文の方程式を見ると、女性および5歳未満の子供を持つ科学者の著しい不利、およびメンターによる出版数の大きな正の効果があり、各論文は出版数の期待値で1.8％の増加と関連する。
モデルが過剰ゼロの問題を解決することを検証するため、πおよびμを予測し、出版なしの結合確率を計算する。 これらを得るために、
predict()関数に"zero"と"count"というオプションがあります。> pr mu zip mean(zip) 0.2985679そこでこのモデルは過剰なゼロの問題を解決し、生化学者の29.9％が論文を発表しないことを予測し、観察値の30.0％にはるかに近いものになりました。 このモデルは
zeroinfl()にdist="negbin"パラメータをつけて適合させることができる。 インフレータブル方程式の代替リンクとしては、link="probit"を使用して指定できるプロビットがある。Model Comparison with AIC
As it happens, for this data, the negative binomial is solves the problem too! 負の二項式で記事が0になる確率はこちらです。 このモデルは、生化学者の30.4%が博士課程の最後の3年間に論文を発表しないことを予測し、30.0%という観測値に非常に近い値です。 異なる数のパラメータを持つモデルを比較する非常に簡単な方法は、
AIC = -2logL + 2p
ここでpはモデル中のパラメータの数であると定義するAkaikeの情報基準（AIC）を計算することである。 第一項は基本的にデビアンス、第二項はパラメータ数に対するペナルティで、
AIC関数を用いてほとんどのモデルで計算することができます。 また、計算を検証できるように「手で」求めています。 我々のデータでは> c(nb=AIC(mnb), zip=AIC(mzip)) nb zip 3135.917 3233.546 > mzip$rank aic sapply(list(mnb,mzip), aic) 3133.917 3233.546このデータセットでは、負の二項モデルがパーシモンと適合度の点で明らかに勝者であった。
我々が見ることができる他の診断基準は、予測されたカウントと観測されたカウントの周辺分布と分散関数である。
ゼロ切り捨てモデルとハードルモデル
我々がカバーしていない他のモデルは、ゼロを含まないデータ用に設計されたゼロ切り捨てポアソンと負の二項モデルである。 一般的な例としては、少なくとも1日以上入院しているような場合です。 賢明なアプローチは、ゼロを除外し、他の確率の和が1になるように再スケールしたポアソンまたは負の2項モデルを適合させることです。 μは期待される結果ではなく、ゼロを含む基本的な分布の平均であるため、これらのモデルの解釈には注意が必要である。 これらのモデルは
VGAMパッケージのvglm()関数でpospoissonとposnegbinomialのファミリーを使って実装されている。ゼロの過剰（または不足）に対する別のアプローチは、ゼロと正の数を区別するためにロジットモデル、次に正の数に対してゼロ切断ポワソンまたは負の二項モデルで2段階プロセスを使用することである。 私たちの例では、出版する人としない人を区別するためにロジット・モデルを使用し、そして、少なくとも1つを出版した人の論文数に対して、切断されたポアソンまたは負の2項モデルを使用することができます。 これらのモデルは、しばしばハードル・モデルと呼ばれます。 これらはRで個別のロジットとゼロカットされたポアソンまたは負の2項モデルを使い、単に対数尤度を加えるか、
psclパッケージのhurdle()関数を使ってフィットさせることができます。ハードル モデルとゼロインフレート モデルを比較すると、ゼロと 1 つ以上の区別はハードル モデルのほうが明確ですが、平均の解釈はゼロインフレート モデルのほうが明確です。 Rにおけるカウントデータの回帰モデル」
。
