通常のポアソン回帰

一般に、ポアソン回帰は0以上の整数値を取るデータに対して適用される回帰モデルで、特に「稀なイベント」に使われています。 例えば、交差点での交通事故発生件数や製品製造ラインでの故障品などが代表的でしょうか。 また、故障品のモデリングなど「全体でいくつ作っているか」も重要になる場合は、 オフセット項を導入することで、割合もモデリングすることが可能になります。FYI: J-Stageの解説論文

また、ポアソン回帰は0/1のbinary outcomeに対しても適用することができます。 イメージする上では、Petersen, et al. 2008のフレーズがわかりやすいように感じます。

It is well known that when the prevalence is low and the sample size is large, probabilities from the Poisson distribution can often be used to approximate probabilities from the binomial distribution. Similarly, one can think of an existing sample of binomial data (0 or 1) as being approximately Poisson, where the probability of a value of 2 or greater is low enough that no values greater than 1 occurred in the obtained sample.

もし、ポアソン回帰をbinary outcomeに対して用いる場合は、Zou, 2004を引用すると良さげです (鉄板論文らしい)。 それでは、実際にすすめていきます。まずは通常のポアソン回帰から。

poisson <- glm(outcome ~ exposure + sex, data = df, family = poisson(link = "log"))
summary(poisson)

## 
## Call:
## glm(formula = outcome ~ exposure + sex, family = poisson(link = "log"), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2286  -0.8966  -0.6704   0.7917   1.1979  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.91140    0.05074 -17.960  < 2e-16 ***
## exposure     0.62997    0.07837   8.038 9.14e-16 ***
## sex         -0.58137    0.08460  -6.872 6.34e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1469.0  on 1999  degrees of freedom
## Residual deviance: 1330.6  on 1997  degrees of freedom
## AIC: 2894.6
## 
## Number of Fisher Scoring iterations: 5

exp(coef(poisson)) # risk ratio

## (Intercept)    exposure         sex 
##   0.4019608   1.8775456   0.5591328

exp(confint(poisson)) # 95%CI, basic poisson regression

##                 2.5 %    97.5 %
## (Intercept) 0.3633534 0.4433375
## exposure    1.6080557 2.1867087
## sex         0.4726644 0.6586717

ここで注目していただきたいポイントは信頼区間です。 今回のポアソン回帰では、RR = 1.88 [95%CI: 1.61-2.19] となりました。

一方、2×2表で得たリスク比は 1.88 [95%CI: 1.70-2.07] であり、ポアソン回帰で得た信頼区間の幅の方が大きくなっています。

ここが、Zou (2004) の言うところの “On the other hand, use of Poisson regression tends to provide conservative results” にあたるのでしょう。

ロバスト分散を使った修正ポアソン回帰

そこで、Zouの手法による修正ポアソン回帰のロバスト分散の計算を実行します。
参考1: Zou’s Modified Poisson Regression, 参考2: いろいろ情報盛りだくさん

pacman::p_load(lmtest, sandwich)

modified_poiss <- coeftest(poisson, vcov = sandwich)

係数は同じですね。一方、標準誤差は随分違います。修正ポアソン回帰の方が、小さな標準誤差を得ることに成功しました。

modified_poiss # 修正ポアソン回帰

## 
## z test of coefficients:
## 
##              Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept) -0.911401   0.039205 -23.2469 < 2.2e-16 ***
## exposure     0.629965   0.050511  12.4718 < 2.2e-16 ***
## sex         -0.581368   0.070023  -8.3025 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

では、修正ポアソン回帰に基づく信頼区間を計算します。 面倒なことに、自分で計算をしなくてはならないようです。

GetConfint <- function(res_coeftest, siglevel=0.05, digits=4){
  temp <- exp(cbind(
    RiskRatio = res_coeftest[,1], 
    LowerCI = res_coeftest[,1] + qnorm(siglevel/2)*res_coeftest[,2],
    UpperCI = res_coeftest[,1] - qnorm(siglevel/2)*res_coeftest[,2]
  ))
  p_value <- res_coeftest[,4]
  
  result <- cbind(temp, p_value = p_value)
  result <- round(result, digits)
  
  return(result)
}

GetConfint(modified_poiss)

##             RiskRatio LowerCI UpperCI p_value
## (Intercept)    0.4020  0.3722  0.4341       0
## exposure       1.8775  1.7006  2.0729       0
## sex            0.5591  0.4874  0.6414       0

2×2表で得たリスク比は 1.88 [95%CI: 1.70-2.07] ですから、通常のポアソン回帰で得た1.88 [95%CI: 1.61-2.19]よりも狭い信頼区間を得ることができました。

※ 2×2表に戻るなら→ データのチェック

モデル式

GLM族なので、a) 線形予測子、b) リンク関数、c) 確率分布の3つで記述することができます。
やっていきましょう。

log-binomial regressionのモデルは以下の通りです（説明変数が2つとします）。
参考: Chen, et al. 2018

\[\log{(p_i)} = b_0 +b_1x_{i1}+b_2x_{i2} \] \[y_i \sim Binomial(n, p_i)\]

ただ、今回は個人が集計単位でy_iは0/1のみがアウトカムなので、n=1となります。
(二項分布なので、コイン投げを1回するイメージですね)

なぜ”log”?

ロジスティック回帰はとても良くできたモデルだと思います。

• 解釈が比較的簡単
• 確率のシグモイド的な振る舞いを考えることができる（「ある閾値を超えるとイベントが起こる」といった考え方）
• 0≤p≤1の確率pに対し、0≤ p/(1-p) <∞であり、log(p/(1-p))は負の無限大から無限大を取るので、被説明変数として適切

そんな中log-binomialモデルでは、なぜlog(p)を考えるのでしょうか？ オッズではなく確率をモデルできることは強みです。 ただし、0≤p≤1の確率pに対し、log(p)は負の無限大から0を取るので、実数全体を取りうるわけではありません。ここがどうにも妙なところです。

この妙な点は、収束が難しいことも引き起こしてしまうようです。

However, the Log link function in Log-Binomial models restricts the probabilities of an outcome to be greater than or equal to zero, that is, to fall within the bounds [0, ∞). Due to this mismatch between the bounds of the model and the allowable outcome, in practice, the Log-Binomial model will routinely fail to converge and will not provide the parameter estimates (Localio et al. 2007).

(from: SAS)

なんだか難しそうですが、とりあえずやってみましょう。

log-binomial regの実践

logbin <- glm(outcome ~ exposure + sex, data = df, family = binomial(link = "log"))

summary(logbin)

## 
## Call:
## glm(formula = outcome ~ exposure + sex, family = binomial(link = "log"), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6735  -1.0139  -0.7141   1.3502   1.7271  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.91148    0.03919 -23.257   <2e-16 ***
## exposure     0.62844    0.04974  12.635   <2e-16 ***
## sex         -0.57996    0.06956  -8.338   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2674.1  on 1999  degrees of freedom
## Residual deviance: 2423.8  on 1997  degrees of freedom
## AIC: 2429.8
## 
## Number of Fisher Scoring iterations: 6

exp(coef(logbin)) # risk ratio

## (Intercept)    exposure         sex 
##   0.4019306   1.8746794   0.5599207

exp(confint(logbin)) # 95%CI

##                 2.5 %    97.5 %
## (Intercept) 0.3714111 0.4331468
## exposure    1.6988379 2.0656133
## sex         0.4868184 0.6397534

2×2表で得たリスク比は 1.88 [95%CI: 1.70-2.07] でしたから、普通にできてますね。

収束の問題が起こるのは、データが十分にない（サンプルサイズ・イベント数など）場合のようです。Williamson, et al. 2013

Recently there was a paper published in Stroke [4], where in the statistical methods section the authors indicated that: “As a first approach to the multivariable analysis, we used a log-binomial model, but owing to the sparseness of data, this failed to converge.

デフォルトの設定では、大きな問題にはなりませんでした。

# startの引数設定
## 説明変数が切片と2つの因子 (曝露と性別) 
## したがって、切片用にlog(mean(y)), 投入した説明変数用にrep(0, 2)となる
logbin2 <- glm(outcome ~ exposure + sex, data = df, family = binomial(link = "log"),
               start = c(log(mean(df$outcome)), rep(0, 2))) 


summary(logbin2)

## 
## Call:
## glm(formula = outcome ~ exposure + sex, family = binomial(link = "log"), 
##     data = df, start = c(log(mean(df$outcome)), rep(0, 2)))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6735  -1.0139  -0.7141   1.3502   1.7271  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.91148    0.03919 -23.257   <2e-16 ***
## exposure     0.62844    0.04974  12.635   <2e-16 ***
## sex         -0.57996    0.06956  -8.338   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2674.1  on 1999  degrees of freedom
## Residual deviance: 2423.8  on 1997  degrees of freedom
## AIC: 2429.8
## 
## Number of Fisher Scoring iterations: 5