カウント（すなわち、0以上の整数）のモデリングに関しては、我々はしばしばポアソン回帰から始めます。 これは一般化線形モデルで、応答が予測変数の重み付き合計の条件付きでポアソン分布を持っていると仮定されるものである。 たとえば、NFLのクォーターバックの脳震盪の記録数を、プレーしたスナップ数とオフェンスラインの合計経験年数の関数としてモデル化することができます。 しかし、ポアソン回帰の潜在的な欠点は、回数のばらつきを正確に表現できない可能性があることである。

ポアソン分布は、その平均と分散の両方を持つ \(\lambda) によってパラメータ化されます。 覚えておくと便利ですが、現実的ではないことが多いです。 カウントの分布は、通常、平均と等しくない分散を持ちます。 ポアソン分布と仮定した（あるいは仮定していた）データでこのようなことが起こった場合、分散が平均より小さいか大きいかによって、過小分散または過大分散と呼ばれます。 このような挙動を示すカウントデータに対してポアソン回帰を実行すると、モデルがうまく適合しないことになります。

この問題に対処する1つのアプローチとして、負の二項回帰がある。 負の二項分布はポアソン分布と同様に、0以上の整数の出現確率を記述する。 このことは、平均値とは異なるばらつきを持つ計数をモデル化するのに有効な近似となる可能性を示唆しています。 負の二項分布の分散は、その平均の関数であり、さらに分散パラメータと呼ばれるパラメータkを持ちます。 カウントを負の二項分布の確率変数Yとすると、Yの分散は

$ var(Y) = \mu + \mu^{2}/k $$

分散パラメータが大きくなると、分散は平均と同じ値に収束し、負の二項分布はポアソン分布に変化する。

負の二項分布を説明するために、Alan Agresti (2002)の著書『Categorical Data Analysis』から、いくつかのデータを使ってみましょう。 データは13.4.3節の表13.6に示されています。 このデータは1308人を対象にしたアンケート調査で、殺人事件の被害者を何人知っているかという質問である。 変数は、respが知っている被害者の数、raceが回答者の人種（黒人か白人か）である。 人種は、その人が何人の殺人被害者を知っているかを説明するのに役立つのでしょうか？

> # Table 13.6> # Agresti, p. 561> black <- c(119,16,12,7,3,2,0)> white <- c(1070,60,14,4,0,0,1)> resp <- c(rep(0:6,times=black), rep(0:6,times=white))> race <- factor(c(rep("black", sum(black)), rep("white", sum(white))),+ levels = c("white","black"))> victim <- data.frame(resp, race)

モデリングに入る前に、データを調査してみましょう。 まず、ほとんどの回答者が白人であることに気づきます。

> table(race)racewhite black 1149 159 

黒人は白人よりも平均カウントが高いです。

> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

それぞれの人種について、サンプルの分散は平均のおよそ2倍である。

> with(victim, tapply(resp, race, var)) white black 0.1552448 1.1498288 

最後に、人種別のカウントの分布を見てみましょう。

> table(resp, race) raceresp white black 0 1070 119 1 60 16 2 14 12 3 4 7 4 0 3 5 0 2 6 1 0

モデルフィッティングに移ります。

> # Poisson model> pGLM <- glm(resp ~ race, data=victim, family = poisson)> summary(pGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.38321 0.09713 -24.54 <2e-16 ***raceblack 1.73314 0.14657 11.82 <2e-16 ***

人種は非常に有意であることがわかります。 黒人は殺人の被害者を知っている可能性が非常に高いようだ。 しかし、係数1.73は何を意味するのだろうか。 この2値予測変数の単純なモデルでは，対数期待カウントの差である．

> exp(coef(pGLM))raceblack 5.658419 > # same thing> mean(victim$resp)/mean(victim$resp) 5.658419

実際にこのモデルで予測し、その結果を指数化すると、標本平均が得られます：

> exp(predict(pGLM, newdata = data.frame(race=c("white","black")))) 1 2 0.09225413 0.52201258 > # same thing> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

これは白人の既知の犠牲者のカウントは平均と分散が0に等しくポアソンとして分布すると言うことです。09、一方黒人の既知の犠牲者の数は、平均と分散が0.52に等しいポアソンとして分布している。 我々は、すでに探索的分析から、観測された分散がもっと大きいことを知っているので、モデルの推定された分散をあまり喜んではいけません。

> # fitted counts for Poisson GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans 1 2 0.09225413 0.52201258 > fittedW <- dpois(0:6,lambda = fmeans) * sum(victim$race=="white") > fittedB <- dpois(0:6,lambda = fmeans) * sum(victim$race=="black") > data.frame(Response=0:6,BlackObs=black, BlackFit=round(fittedB,1), + WhiteObs=white, WhiteFit=round(fittedW,1)) Response BlackObs BlackFit WhiteObs WhiteFit1 0 119 94.3 1070 1047.72 1 16 49.2 60 96.73 2 12 12.9 14 4.54 3 7 2.2 4 0.15 4 3 0.3 0 0.06 5 2 0.0 0 0.07 6 0 0.0 1 0.0

上記では、まず予測された平均を fmeans というオブジェクトに保存しました。 次に、0から6を得る確率を予測するために、推定された平均とともにdpois関数を使用して適合カウントを生成しました。 そして、これらの確率に回答者数を掛けて、適合カウントを得ました。 最後に、観測値と適合値を簡単に比較できるように、すべてをデータフレームにまとめました。 より顕著な点として、0カウントはアンダーフィット、1カウントはオーバーフィットになることが挙げられます。

ルートグラムを使用して、カウント回帰モデルの適合度を視覚化することができます。 countregパッケージのrootogram()関数を使うと簡単にできます。

> countreg::rootogram(pGLM)

赤い曲線が理論上のポアソンフィットです。 赤い線上の各点から “垂れ下がる “棒は、その高さが期待値と観測値の差を表している。 0以下の棒はアンダーフィットを示す。 0より大きい場合はオーバーフィットである。 カウントは、小さなカウントが大きなカウントに埋もれてしまわないように、平方根変換されている。 カウント 2 以上はアンダーフィット、カウント 1 はオーバーフィットが大きいことがわかります。 カウントの変動はどちらの人種でも大きいことに気がついた。

Rで負の二項モデルを適合させるために、MASSパッケージ（Rにインストールされているパッケージ）のglm.nb()関数に注目します。

> # Fit a negative binomial model> library(MASS)> nbGLM <- glm.nb(resp ~ race, data=victim)> summary(nbGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.3832 0.1172 -20.335 < 2e-16 ***raceblack 1.7331 0.2385 7.268 3.66e-13 *** Theta: 0.2023 

最初に、係数が以前と同じであることに注意してください。 もう一度、レース係数を指数化して標本平均の比率を求め、予測をして元の標本平均を得ることができます。

> # fitted counts for Negative Binomial GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans # same as pGLM 1 2 0.09225413 0.52201258 

しかし、人種係数の標準誤差が大きく、我々の推定値の不確実性が高いことに注目してください（0.15に対して0.24）。 これは、カウントで観察されたばらつきを考えると、理にかなっています。 また、θの推定値にも注目してください。 これは分散パラメータです。

> nbGLM$theta 0.2023119

そして、カウントの推定分散を得るためにそれを使用できます。

> fmeans + fmeans^2 * (1/nbGLM$theta) 1 2 0.134322 1.868928 

これらはポアソンモデルによって与えられるものよりもはるかに観察された分散に近いです。

> countreg::rootogram(nbGLM)

これはポアソンモデルのルートグラムよりはるかに良いようです。 カウント1から3まではわずかにアンダーフィット/オーバーフィットしていますが、それ以外はかなり良さそうです。

我々の負の二項モデルに対するさらなる洞察を得るために、そのパラメータを使ってデータをシミュレートし、シミュレートされたデータを観察されたデータと比較してみましょう。 以下では、関数の「連鎖」を可能にする%>%演算子にアクセスするために、magrittrパッケージをロードします。

まず、白人と黒人の参加者それぞれについて、カウントを集計し、棒グラフを作成します。 次に、負の二項分布からデータをシミュレートするために、モデル・パラメータを使用します。 2つのパラメータはmuとsize（つまり、分散パラメータ）です。 各人種の適切な係数を抽出するために、coef()関数を使用していることに注意してください。 白人の場合は切片だけで、黒人の場合は切片と傾きです（したがって、それらを合計します）。 そして、対数スケールから元のスケールに変換するために、指数関数で計算します。

library(magrittr) # for %>% operatorop <- par(mfrow=c(1,2))set.seed(1)victim$resp %>% `)) %>% table() %>% barplot(main = "Simulated White")victim$resp %>% `[`(victim$race=="black") %>% table() %>% barplot(main = "Observed Black")rnbinom(n = 159, size = nbGLM$theta, mu = exp(sum(coef(nbGLM)))) %>% table() %>% barplot(main = "Simulated Black")par(op)

シミュレーションデータは、観測データと非常に似ており、このデータをモデルするために負の二項回帰を選ぶことに自信を与えてくれます。 これらのプロットはまた、我々のモデルの条件付の性質を示しています。 カウントの負の二項分布は、人種に依存する、あるいは条件付きです。 各人種の平均は異なるが、分散パラメータは共通である。

• Agresti, Alan (2002), Categorical Data Analysis, Wiley.
• Kleiber, Christian & Zeileis, Achim (2016): Visualizing Count Data Regressions Using Rootograms, The American Statistician, DOI: 10.1080/00031305.2016.1173590

この記事に関する質問や不明点はUVA Library StatLabにお問い合わせください：[email protected]

UVA図書館StatLab記事のコレクション全体を見る

Clay Ford
統計研究コンサルタント
University of Virginia Library