過分散データへのGLM適用時のp値

つづいて、過分散データにGLM（誤差構造はポアソン分布、リンク関数はlog）を適用したときの検定はどうなるか、をやってみます。

きのうと同様に関数を定義しますが、返り値は切片および係数のp値としています。

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

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

X2の係数は0ですので、危険率5%で、p < 0.05となるのは5%のはずです。

> set.seed(123)
> test2 <- replicate(2000, test2.f())
> sum(test2["X2", ] < 0.05) / 2000
[1] 0.2525

しかし、25%ほど「有意」になってしまいます。

ランダム効果の大きさを小さくすると、5%に近づきます。

> set.seed(123)
> test2 <- replicate(2000, test2.f(sd.B = 1e-6))
> sum(test2["X2", ] < 0.05) / 2000
[1] 0.055

タグ：GLM