過分散データへのGLM適用時のパラメータ推定値 [統計]
過分散データをむりやりGLMにあてはめたとき、パラメータの推定値がどうなるか実験してみました。
Rで、以下のような関数を定義します。ポアソン分布でデータを生成しますが、ブロックによるランダム効果を含めているので、過分散となります。
test1.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, 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(coef(g1)) }
1回だけ実行してみます。
> set.seed(123) > test1.f(summary = TRUE) Call: glm(formula = Y ~ X1 + X2, family = poisson) Deviance Residuals: Min 1Q Median 3Q Max -2.5625 -1.0869 -0.2506 0.7244 3.2706 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.90249 0.04558 19.799 < 2e-16 *** X1 1.14140 0.22602 5.050 4.42e-07 *** X2 0.19364 0.14879 1.301 0.193 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 373.63 on 199 degrees of freedom Residual deviance: 347.18 on 197 degrees of freedom AIC: 825.38 Number of Fisher Scoring iterations: 5 (Intercept) X1 X2 0.9024888 1.1414028 0.1936370
残差の自由度が197で、残差逸脱度が347.18ですから、しっかり過分散になっています。
2000回繰り返します。
> set.seed(123) > test1 <- replicate(2000, test1.f()) > apply(test1, 1, summary) (Intercept) X1 X2 Min. 0.248 -0.6141 -1.036000 1st Qu. 1.073 0.7673 -0.153800 Median 1.267 0.9855 -0.003162 Mean 1.269 0.9892 -0.003051 3rd Qu. 1.455 1.2090 0.144500 Max. 2.472 3.0440 0.964900
X1, X2の係数の平均値はそれぞれデータを生成したモデルの値とだいたい一致していますが、切片の値が大きくなるようです。
では、ランダム効果を大きくして再度ためしてみます。
> set.seed(123) > test1 <- replicate(2000, test1.f(sd.B = 1.2)) > apply(test1, 1, summary) (Intercept) X1 X2 Min. 0.1888 -1.4340 -1.361000 1st Qu. 1.2860 0.7123 -0.177100 Median 1.5790 1.0140 -0.000645 Mean 1.6050 0.9926 0.008430 3rd Qu. 1.9030 1.2800 0.196300 Max. 4.0300 3.3150 1.289000
切片の平均値はさらに大きくなりました。ランダム効果が大きくなると、大きな値が出やすくなる一方、小さい値は0までしかありませんから、このように偏りが生じるのではないかと思われます。
GLMMをつかうと、このような偏りは生じませんでした。test1.fのglmの部分を以下のように変更します。
g1 <- glmmML(Y ~ X1 + X2, cluster = B, family = poisson)
GLMMの結果です。
> set.seed(123) > test1 <- replicate(2000, test1.f(sd.B = 1.2)) > apply(test1, 1, summary) (Intercept) X1 X2 Min. -0.09565 0.2740 -0.4342000 1st Qu. 0.74440 0.8817 -0.0740600 Median 0.99770 1.0020 -0.0003356 Mean 1.00200 1.0020 0.0001078 3rd Qu. 1.26900 1.1190 0.0717100 Max. 2.17700 1.6680 0.4227000
コメント 0