負の二項回帰モデル R で行う方法

負の二項回帰モデルを R で行う方法

↑1万人以上の医療従事者が購読中

負の二項モデルを R で行う方法
まとめ

負の二項モデルを R で行う方法

まず、データを読み込む

X <- c(rep(1, 913), rep(2, 1345), rep(3, 519), rep(9, 863))
Y <- c(rep(0, 910), rep(1, 2), rep(2, 1), rep(4, 0),
       rep(0, 1333), rep(1, 11), rep(2, 1), rep(4, 0),
       rep(0, 509), rep(1, 9), rep(2, 1), rep(4, 0),
       rep(0, 849), rep(1, 13), rep(2, 0), rep(4, 1))
dat <- data.frame(X, Y)

Yの分布を確認

> table(dat$Y)

   0    1    2    4 
3601   35    3    1

Yの分布図はすごいことに、なっている

plot(table(dat$Y))

XとYの関係を見ていこう

まずクロス集計表

> table(dat$X, dat$Y)
   
       0    1    2    4
  1  910    2    1    0
  2 1333   11    1    0
  3  509    9    1    0
  9  849   13    0    1

Y がゼロばかりなのが目立つ

負の二項回帰モデルは、MASS パッケージの glm.nb() を使う

library(MASS)
nb.res <- glm.nb(Y ~ factor(X), data=dat)
summary(nb.res)

結果は、Xが１に比べて、３とか、９とかはYに関連がある

> summary(nb.res)

Call:
glm.nb(formula = Y ~ factor(X), data = dat, init.theta = 0.05009536409, 
    link = log)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.18802  -0.18228  -0.13295  -0.09165   3.10151  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -5.4304     0.5214 -10.415   <2e-16 ***
factor(X)2    0.7912     0.6030   1.312   0.1895    
factor(X)3    1.5764     0.6334   2.489   0.0128 *  
factor(X)9    1.5032     0.5948   2.527   0.0115 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.0501) family taken to be 1)

    Null deviance: 210.18  on 3639  degrees of freedom
Residual deviance: 199.65  on 3636  degrees of freedom
AIC: 465.67

Number of Fisher Scoring iterations: 1


              Theta:  0.0501 
          Std. Err.:  0.0236 

 2 x log-likelihood:  -455.6660