試行数と成功数が与えられたデータに対して二項回帰モデルを当てはめる [R]

アウトカムが2値変数のときに使う回帰モデルについては、以前まとめたことがあります。

necostat.hatenablog.jp

今回はアウトカムが「割合」、というか「試行数と成功数」として与えられているときの回帰モデルについてまとめてみます。

使用するデータ
cbind(成功数, 試行数)で二項モデルを当てはめる
割合を応答変数として試行数で重み付けして二項モデルを当てはめる
おわりに
参考資料

使用するデータ

今回はsegmentedパッケージのdownデータを使います。

library(segmented)
data(down)

このデータは母体年齢と出生児のダウン症（21トリソミー）の発生の関連を示したデータで、以下の3つの変数を含みます。

age：母体の平均年齢
births：出生数
cases：ダウン症を有する児の数

実際のデータはこんな感じ。

> head(down)
   age births cases
1 17.0  13555    16
2 18.5  13675    15
3 19.5  18752    16
4 20.5  22005    22
5 21.5  23896    16
6 22.5  24667    12

cbind(成功数, 試行数)で二項モデルを当てはめる

発生割合が（アウトカムの確率分布の）パラメータになっているモデルといえば、二項-ロジットモデル（ロジスティック回帰モデル）ですね。glm( )を使って、

fit <- glm(cases ~ age, family = binomial(link=logit), data = down)

と書きたいところですが、casesは2値変数ではないのでエラーになります。

こういうときは、cbind( )を使って次のように書きます。

fit.bin <- glm(cbind(cases, births) ~ age, 
               family = binomial(link=logit), data = down)

> summary(fit.bin)

Call:
glm(formula = cbind(cases, births) ~ age, family = binomial(link = logit), 
    data = down)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4053  -1.9359   0.5364   2.1314   4.7481  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.545157   0.213919   -49.3   <2e-16 ***
age           0.136872   0.006456    21.2   <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: 619.96  on 29  degrees of freedom
Residual deviance: 181.27  on 28  degrees of freedom
AIC: 324.15

Number of Fisher Scoring iterations: 5

リンク関数がロジットなので、係数を指数変換すればオッズ比として解釈できます。

割合を応答変数として試行数で重み付けして二項モデルを当てはめる

発生数を出生数で割って、発生割合にしてアウトカム変数とする方法です。

fit.prop <- glm(cases/births ~ age, 
               family = binomial(link=logit), data = down)

とすると「成功数が整数になっていない！」と怒られます。

次のように、引数weightに試行数（ここでは出生数）を指定します。試行数が大きいところの方が信頼度が高いので、全体の結果に寄与する程度も大きいからです。

fit.prop <- glm(cases/births ~ age, weight = births, 
               family = binomial(link=logit), data = down)

> summary(fit.prop)

Call:
glm(formula = cases/births ~ age, family = binomial(link = logit), 
    data = down, weights = births)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-3.412  -1.943   0.546   2.140   4.767  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.562280   0.214342  -49.28   <2e-16 ***
age           0.137524   0.006469   21.26   <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: 625.21  on 29  degrees of freedom
Residual deviance: 183.86  on 28  degrees of freedom
AIC: 326.74

Number of Fisher Scoring iterations: 5

だいたい似たような推定値になりました。

おわりに

(Dispersion parameter for binomial family taken to be 1)というメッセージが出てきているので、過剰分散（overdispersion）の話とか、擬似二項分布（quasibinomial）の話とかをまとめようと思っていますが、それはまた別の機会に。
「いやいや、割り算なんかせずにオフセット項のあるポアソン回帰だろ！」というご意見もあると思いますので、これも別の機会に。