過分散データに負の二項分布のモデルをあてはめ

過分散データに負の二項分布のモデルをあてはめて、p値をみてみました。MASSパッケージのglm.nb関数を使用。

以下のような関数を定義しました。

library(MASS)

test4.f <- function(N = 200, NB = 10,
                    m1 = 0, m2 = 0,
                    sd1 = 0.2, sd2 = 0.3, sd.B = 0.8,
                    intercept = 1, coef.X1 = 1,
                    summary = FALSE) {
    X1 <- rnorm(N, 0, sd1)
    X2 <- rnorm(N, 0, sd2)
    B <- rep(rnorm(NB, 0, sd.B), each = N / NB)
    Y <- rpois(N, exp(intercept + coef.X1 * X1 + B))

    g1 <- glm.nb(Y ~ X1 + X2)
    if (summary) print(summary(g1))
    return(summary(g1)$coefficients[, "Pr(>|z|)"])
}

結果です。

> set.seed(1234)
> test4 <- replicate(2000, test4.f())
 警告メッセージ: 
 glm.nb(Y ~ X1 + X2) で:  alternation limit reached
> sum(test4["X2", ] < 0.05) / 2000
[1] 0.057

2000回の繰り返しのなかには収束しないデータも生成されたようですが、p<0.05となる割合はだいたいよい値でした。

タグ：負の二項分布 GLM