SSブログ

過分散データへの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


タグ:glmm GLM
nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:日記・雑感

nice! 1

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0