R のロジスティック回帰: glm() と family = binomial

結果が yes/no のとき

線形回帰は数値を予測します。しかし、モデル化する価値のある問いの多くは二値です。顧客は解約するか、患者は回復するか、メールはクリックされるか。0/1 の結果に直線を当てはめるとたちまち破綻します。直線は平然と −0.3 や 1.4 という確率を予測しますが、それは意味をなしません。

ロジスティック回帰は、対数オッズ（ロジット）変換を通じて結果の 確率 をモデル化することでこれを解決します: log(p / (1 − p)) = intercept + slope × x。対数オッズのスケールは数直線全体にわたるので、そこでは線形の方程式が自然に当てはまります。そして元に戻すと、すべての予測がおなじみの S 字曲線に沿って (0, 1) に押し込められます。この技の代償は、係数が対数オッズのスケール上に存在することです。ロジスティック回帰を 読む という仕事のすべては、それを人間が理解できる何かに翻訳し直すことにあります。

当てはめ: glm() と family = binomial

glm()（一般化線形モデル）は lm() の兄貴分であり、family = binomial がロジスティック回帰を選択します。mtcars では am がトランスミッションの種類（1 = マニュアル、0 = オートマチック）を記録しています。燃費のよい車はマニュアルである傾向があるでしょうか?

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

出力を読む前に 2 点あります。第一に、family = binomial は省略可能ではありません。省くと glm() は黙って通常の最小二乗を当てはめます。第二に、結果変数は二値でなければなりません。0/1、論理値、または 2 水準の因子です（R は 第 2 水準の確率をモデル化します）。

では summary をブロックごとに見ていきましょう。

• Coefficients — mpg の Estimate は約 0.31 で、これは対数オッズの傾きです。mpg が 1 増えるごとに、マニュアルである対数オッズが 0.31 加算されます。正なら「確率を上げる」、負なら「下げる」という意味ですが、符号を越えると誰の直感もこのスケールでは働きません。だからこそ次の節が存在します。
• z value と Pr(>|z|) — 回帰の t 検定と同じ論理（Estimate ÷ Std. Error、そして「これはゼロでありうるか?」の p 値）ですが、正規近似を使うため t ではなく z になります。ここでは p ≈ 0.011 で、mpg とトランスミッションの種類の関連はノイズである可能性が低いということです。
• Null deviance と Residual deviance — 逸脱度は glm の世界での当てはまりの悪さです（小さいほど良い）。Null deviance（自由度 31 で 43.2）は切片のみのモデル、Residual deviance（自由度 30 で 29.7）があなたのモデルです。説明変数 1 個分の自由度に対して約 13.6 減少したことが、「決定係数が上がった」に相当する glm 版の指標です。
• AIC — 当てはまりと複雑さのバランスを取るモデル比較のスコアで、低い方が勝ちます。単独では意味がなく、同じデータ上の候補モデル間で有用です。

対数オッズからオッズ比へ: exp(coef())

指数変換により、係数は加法的な対数オッズのスケールから乗法的な オッズ のスケールへ移ります。

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

exp(0.307) ≈ 1.36 です。そして覚えるべき誠実な文型がこれです。「mpg が 1 増えるごとに、マニュアルトランスミッションであるオッズが約 1.36 倍になる。」 オッズ比が 1 より大きければオッズを上げ、1 より小さければ下げ、ちょうど 1 なら効果なしです。だからこそ、オッズ比に対する信頼区間の判断基準は「区間が 1 を含まないか?」になります（ゼロではありません。ゼロが境界だったのは対数オッズのスケール上での話です）。

言葉づかいに注意してください。オッズは確率ではありません。オッズ = p / (1 − p) なので、確率 0.75 はオッズ 3 です。オッズ を 1.36 倍することは、確率 を 1.36 倍することとは違います。そして結果が頻繁に起こる場合、その差は大きくなります。まれな結果に対するオッズ比 2 は「リスクがおよそ 2 倍」のように振る舞いますが、50% で起こる結果に対してはまったくそうなりません。結果がまれでない限り、オッズ比をリスク比の言い回し（「1.36 倍起こりやすい」）で報告してはいけません。

予測確率: type = "response" の落とし穴

実務で最もよく見かけるロジスティック回帰のバグがこれです。

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

最初の呼び出しは 既定の type = "link" を返します。負の値も含む、対数オッズのスケール上の予測です。2 つ目は実際の確率を返します。「確率」が負になったり 1 を超えたりするのは、これが理由です。ブロックを実行してみましょう。15 mpg の車はマニュアルである見込みがほぼなく、30 mpg の車は非常にマニュアルらしく、S 字曲線がその中間を曲がっていきます。

分類: 閾値と混同行列

確率は閾値を選ぶことで予測クラスになります。既定の選択は 0.5 です。そして誠実な成績表は、予測と実測の表です。

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

対角のセルが正しい判定、対角以外の 2 つのセルが 2 種類の 異なる 誤りです（オートマチックをマニュアルと予測する誤りと、その逆）。全体の正解率だけではモデルをひどく美化しかねません。顧客の 95% が解約しないなら、「誰も解約しないと予測する」は 95% を得点しつつ、解約者を 1 人も捕まえません。ですから常に両方の誤りの種類を見ましょう。そして 0.5 は慣習であって法律ではありません。2 種類の誤りのコストが異なるなら、それに応じて閾値を動かしましょう。

誠実さのための但し書きがひとつ。この表は、モデルを当てはめたのと同じデータでモデルを採点しており、その分だけ甘い評価になります。本当の評価は、モデルが一度も見ていないデータを取り置いて行います。

複数の説明変数

lm() とまったく同じで、+ で項を追加します。そしてすべての解釈に「他を一定に保ったうえで」という但し書きが付きます。

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

指数変換した各係数は、いまや 他の説明変数について同様である車の中で その説明変数が 1 単位増えたときのオッズの倍率です。仕組みはそのまま拡張されますが、線形回帰での但し書きも同様に拡張されます。相関する説明変数は互いの係数を組み替えます。

注意点

• 完全分離。 ある説明変数が結果を完全に分割してしまう場合（ある mpg より上の車はすべてマニュアル、下はすべてオートマチック）、最尤法の係数は無限大になろうとします。R は glm.fit: fitted probabilities numerically 0 or 1 occurred と警告し、途方もない標準誤差を伴う巨大な係数を報告します。その数値を世に出してはいけません。モデルを単純にするか、データを増やすか、罰則付きの手法（brglm2 や logistf パッケージ）を使いましょう。
• 十分なイベント数。 制約となるのは総行数ではなく、まれな方 の結果の件数です。古い経験則では説明変数 1 つあたり 10〜15 件のイベントが求められます。ここでの 32 台の例は仕組みを教えるためのものであり、論文に載せられる標本サイズの手本ではありません。
• オッズ比はリスク比ではありません — 結果が頻繁に起こる場合の話で、上で扱いましたが、あなたが見落としても査読者は見つけるので繰り返しておきます。

この章のまとめ

• 二値の結果 → glm(y ~ x, data = df, family = binomial)。family を決して忘れないこと。
• 生の係数は対数オッズです。exp(coef(fit)) がオッズ比を返し、その区間にとっての帰無値は 1 です。
• 文型はこうです。「x が 1 単位増えるごとに、結果のオッズが exp(b) 倍になる。」
• 確率には predict(..., type = "response") を使います。既定は対数オッズを返し、これが最大の混乱の元です。
• 閾値で分類し、混同行列で評価しましょう。正解率だけでは嘘をつくことがあります。
• 分離の警告に注意し、イベント数を数え、オッズ比をリスク比に見せかけないこと。

次は、ここまで見てきたすべての区間の裏側にある仕組み — t.test()、confint()、prop.test() による信頼区間です。