可視化で理解しない「負の二項分布」
匿名知的集団であるホクソエムの親分が可視化で理解する「負の二項分布」 - ほくそ笑むで言及している負の二項分布ですが、シンクロニシティというか、同じ匿名知的集団ホクソエムなので当たり前と言えば当たり前なのですが、私も負の二項分布、特に負の二項分布のパラメータを推定するための関数であるglm.nb関数の挙動を調べる必要があったのでまとめる。数式変形も大体理解しているが、参考テキストや資料によって表記のギリシア文字のぶれがあるので、その辺どれを使うのがいいのか決めてからまとめたい。
ダミーデータの作成
glm.nb関数に食わせるためのダミーデータをrnegbin関数(負の二項分布に従う乱数を生成)作成する。もちろんシードは71だ。
ここで
- 負の二項分布の平均はexp(3.7*x+2)、形状パラメータは10.7
と適当に設定している。このあたりのパラメータの呼称が混乱の原因になるように思うが、ここでは
でいう`shape parameter`を用いて書いた形で議論する(そのほうがおいおいGamma–Poisson mixtureとしての負の二項分布としての見通しがよくなる)。
set.seed(71) size <- 10^4 x <- -(1:size)/size shape <- 10.7 mu <- exp(3.7*x+2) y1 <- sapply(mu, function(m)rnegbin(1, mu=m, theta=shape))
モデルの推定
負の二項分布に従うモデルのフィッティングを行うためにglm.nb関数で推定を実施する。表記法は(g)lm関数などと同様formula式でOK。
> library("MASS") > model_glmnb <- glm.nb(y1 ~ x) > summary(model_glmnb) Call: glm.nb(formula = y1 ~ x, init.theta = 9.775593004, link = log) Deviance Residuals: Min 1Q Median 3Q Max -3.1997 -0.9142 -0.4375 0.5292 3.6430 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.02512 0.01309 154.7 <2e-16 *** x 3.76445 0.03668 102.6 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for Negative Binomial(9.7756) family taken to be 1) Null deviance: 24389 on 9999 degrees of freedom Residual deviance: 10194 on 9998 degrees of freedom AIC: 29351 Number of Fisher Scoring iterations: 1 Theta: 9.776 Std. Err.: 0.746 2 x log-likelihood: -29345.396
ここまでの使い方は適当にググるとでてくるのだが、問題はこの結果の解釈で、それについてはあまり言及されていない気がするので詳細に書くと。。。
- 切片(Intercept)のはおおよそ: 2.02512
- 特徴量(説明変数)xの係数: 3.76445
- Theta(shape変数の値): 9.776
と読む。この係数(2.02512、3.76445、9.776)は元のダミーデータを設定した際の値がほぼ再現されており、良い推定を与えていることがわかる。
また、ここで推定されているもの(predict, fitted関数でとれるもの）は、（lm関数などの結果から類推できるが、）上でいう負の二項分布の平均値(mu)の推定値になる。実際、モデルオブジェクトに入っている値は以下のようにexp(係数×データ)で計算されるもので、上で書いたmuの定義exp(3.7*x+2)そのものになっている。以下では同じ答え（推定値）を出す三通りの方法＋厳密な値（私が指定）を算出している。推定値と厳密な値も当然近しいものになる。
> head(fitted(model_glmnb, x)) 1 2 3 4 5 6 7.574199 7.571348 7.568499 7.565650 7.562803 7.559956 > head(model_glmnb$fitted.values) 1 2 3 4 5 6 7.574199 7.571348 7.568499 7.565650 7.562803 7.559956 > head(as.numeric(exp(cbind(1, x) %*% beta))) [1] 7.574199 7.571348 7.568499 7.565650 7.562803 7.559956 > head(mu) [1] 7.386323 7.383590 7.380859 7.378128 7.375399 7.372671
元のデータの復元
ここから先は上述の親分のブログにあるように、元のデータと同じようにサンプリングされたデータを
- ガンマ分布からサンプリングする
- ポアソン分布からサンプリングする
という手順で再現する。そして、その頻度分布が近しいものですねと、ただしくサンプリングされていますねということをチェックする。
まずは、実際に設定したパラメータを用いてサンプリングする。ここで、scale=m/shapeと設定しているが、これは数式変形をがんばると見えてくるところだ。
lambda <- sapply(mu, function(m){rgamma(1, shape=shape, scale=m/shape)}) y2 <- sapply(lambda, function(l)rpois(1, l)) plot(density(y1), col="red") lines(density(y2), col="blue")
となり、ほぼ一致する。
同様の作業を推定したモデルで実行すると
beta <- model_glmnb$coefficients theta <- model_glmnb$theta lambda <- sapply(model_glmnb$fitted.values, function(m){rgamma(1, shape=theta, scale=m/theta)}) y3 <- sapply(lambda, function(l)rpois(1, l)) hist(y1, breaks=0:max(y1, y3), col="red") hist(y3, breaks=0:max(y1, y3), col="blue", add=TRUE)
となり、こちらもほぼ一致する。
