Firthのロジスティック回帰とは ― 完全分離・希少イベントで推定が破綻する問題をRで解決する ―
この記事でわかること
- 完全分離・準完全分離が起きると、通常のロジスティック回帰の最尤推定（MLE）はなぜ破綻するのか
- 推定破綻はエラーで教えてくれるとは限らず、どのような症状（標準誤差・p値・オッズ比・信頼区間）で気づけばよいか
- Firth（1993）のペナルティ付き最尤法が、なぜ分離下でも有限の推定値を返せるのか
- ペナルティ項がJeffreysの事前分布に対応するというベイズ的な解釈
- 希少イベント（rare events）で生じる小標本バイアスを、Firth法がどのように補正するのか
- RのlogistfパッケージでFirthのロジスティック回帰を実装し、通常のglmと比較する具体的な方法
はじめに
臨床試験や市販後調査では、ある有害事象（AE）がごく少数の被験者にしか発生しない、あるいは特定の群では一件も発生しないという状況が珍しくありません。このような「希少イベント」を目的変数としてロジスティック回帰を当てはめようとすると、通常の最尤推定（MLE）が正しく機能しなくなることがあります。これが本記事のテーマである「完全分離（complete separation）」および「準完全分離（quasi-complete separation）」の問題です。
厄介なのは、この推定破綻が必ずしもエラーとして明示的に知らされるとは限らない点です。Rのglm()は何のエラーも出さずに、標準誤差が異常に大きく、オッズ比が桁外れの値になった結果を平然と返すことがあります。解析者がその症状に気づかなければ、意味のない推定結果をそのまま報告してしまいかねません。
この問題に対する実務的な解決策が、Firth（1993）が提案したペナルティ付き最尤法、いわゆる「Firthのロジスティック回帰」です。尤度関数にペナルティ項を加えることで、分離が生じている状況でも必ず有限の推定値を得ることができ、さらに希少イベントで生じやすい小標本バイアスも同時に低減できます。
本記事は、製薬企業の生物統計担当者、R・SASを用いる統計解析実務者、そして統計学を学ぶ大学院生・学生を主な対象としています。小規模な臨床試験や稀な有害事象の解析でロジスティック回帰がうまく収束しない、あるいは信じがたいオッズ比が出てしまったという経験がある方に、特に読んでいただきたい内容です。前半では完全分離とは何か、Firth法がなぜ機能するのかという数理的な背景を整理し、後半ではRのlogistfパッケージを使った具体的な実装方法と結果の解釈を、実測データに基づいて解説します。
Firthのロジスティック回帰とは
通常のロジスティック回帰では、目的変数（イベントの有無）を最もよく説明するパラメータを、最尤法（Maximum Likelihood Estimation, MLE）によって求めます。この推定が破綻する典型的な状況が「完全分離」です。完全分離とは、ある説明変数（あるいは説明変数の組み合わせ）の値によって、目的変数の実現値（0か1か）が完全に予測できてしまう状態を指します。たとえば、ある群では有害事象が一件も発生せず、別の群では発生している、という状況がこれに当たります。
このとき、尤度関数はその群に対応する係数を無限大（または無限小）に近づけるほど大きくなり続けてしまうため、最尤推定量は数学的に存在しません。推定アルゴリズムは反復計算の中で係数を発散させ続け、最終的に非常に大きい（しかし本来は無限大であるべき）値で計算が打ち切られます。目的変数が完全にではなく、ごく一部の重複を除いてほぼ予測できてしまう場合は「準完全分離」と呼ばれ、これも同様に推定を不安定にします。
分離が起きても、Rは必ずしもエラーや警告で教えてくれるとは限りません。実務では、①標準誤差が数百〜数千と桁外れに大きい、②z値がほぼ0でp値が1に近い、③オッズ比が天文学的な数値になる、④信頼区間が(0, ∞)に近づく、といった症状から自分で気づく必要があります。後述のRでの実測例で、この症状を具体的に確認します。
希少イベント（rare events）、つまり目的変数の1（イベントあり）が少ないデータでは、この完全分離・準完全分離が起きやすくなります。さらに、完全に分離していなくても、イベント数が少ないだけでMLEには小標本バイアスが生じ、オッズ比を過大に見積もる傾向があることが知られています。製薬領域では、稀な有害事象や小規模な第II相試験など、まさにこの状況に直面する場面が少なくありません。
Firth（1993）が提案したペナルティ付き最尤法は、この2つの問題――分離による推定破綻と、希少イベントによる小標本バイアス――の両方に対して有効です。Firth法は尤度関数にペナルティ項を加えることで、分離が生じている状況でも必ず有限の推定値を返し、さらにMLEが持つ一次のバイアスを除去します。次章では、このペナルティがどのような数理的根拠に基づいているのかを見ていきます。
Firth法の数理的背景
通常のロジスティック回帰のMLEは、対数尤度 \( l(\beta) \) を最大化するパラメータ \( \beta \) を求めます。分離が生じている状況では、この最大化問題に有限の解が存在しないことは前章で述べたとおりです。Firth法はこの対数尤度に補正項を加えた、次の「ペナルティ付き対数尤度」を最大化します。
\[ l^{*}(\beta) = l(\beta) + \tfrac{1}{2}\log\bigl|\,I(\beta)\,\bigr| \]
ここで \( I(\beta) \) はFisher情報行列（パラメータの推定精度に関する情報量を表す行列）、\( |I(\beta)| \) はその行列式です。通常の対数尤度 \( l(\beta) \) に \( \tfrac12\log|I(\beta)| \) という項を加えることで、係数が発散しようとするとペナルティが働き、尤度の増加が頭打ちになるように調整されています。
このペナルティ項には、ベイズ統計の観点からも意味があります。\( \tfrac12\log|I(\beta)| \) はJeffreysの不変事前分布（パラメータの変換に対して不変な、無情報に近い事前分布）の対数密度に一致することが知られています。つまりFirth推定量は、「Jeffreys事前分布を置いたときの事後分布の最頻値（MAP推定量）」というベイズ的な解釈を持っています。頻度論の枠組みで導出された手法でありながら、ベイズ的な正則化とも自然に結びついている点がFirth法の特徴です。
実際の推定では、対数尤度を微分して得られる通常のスコア方程式 \( U(\beta) = 0 \) の代わりに、次の修正スコア方程式を解きます。
\[ U(\beta)_j^{*} = U(\beta)_j + \tfrac12\,\mathrm{tr}\{I(\beta)^{-1}\,\partial I(\beta)/\partial\beta_j\} = 0 \]
右辺第2項がFirthによる補正分です。2×2分割表の場合、この補正は「観測されたセル度数それぞれに0.5を足してからオッズ比を計算する」という、疫学で古くから知られるHaldaneの連続性補正とほぼ同じ効果を持つことが知られています。直感的に言えば、Firth法はどのセルにも「半分だけの仮想的なイベント」を足すことで、実現値が0であるセルからオッズが0または無限大になってしまうことを防いでいます。
この補正の効果は大きく2つあります。第一に、MLEが本来持つ \( O(n^{-1}) \) のオーダーの一次バイアスを解消し、特に希少イベントでオッズ比が過大評価される問題を緩和します。第二に、完全分離・準完全分離が生じている状況でも、パラメータの推定値が必ず有限の値に収束します。なお、信頼区間や検定については、通常のWald法（標準誤差から正規近似で構成する方法）は分離下では信頼できないため、Firth法ではペナルティ付きプロファイル尤度に基づく信頼区間・検定を用いるのが既定の方法です（Rのlogistfパッケージの既定設定でもこの方法が使われます）。次章では、実際にRでこれらを確認していきます。
RでFirthのロジスティック回帰を実装する
ここからは実際にRを使って、完全分離が起きるデータに対して通常のロジスティック回帰（glm）がどのように破綻し、Firthの方法（logistfパッケージ）がどう解決するのかを、実測結果で確認していきます。環境はR 4.5.1、logistf 1.26.1です。
ステップ1：完全分離が起きるデータを用意する
小規模な第II相試験を想定します。対照群45例でAE（有害事象）が1件も起きず（0/45）、実薬群45例で7件発生した（7/45）という状況です。対照群のイベントが0であることが、まさに完全分離を生みます。
d1 <- data.frame(
group = factor(rep(c("control","active"), each = 45),
levels = c("control","active")),
ae = c(rep(0L, 45), c(rep(1L, 7), rep(0L, 38)))
)
table(group = d1$group, AE = d1$ae)
tapply(d1$ae, d1$group, mean)
> AE
> group 0 1
> control 45 0
> active 38 7
>
> control active
> 0.000 0.156
対照群はイベント0/45（発生率0.0%）、実薬群は7/45（15.6%）です。対照群のセルが0であるため「対照群のオッズ」が0となり、実薬 vs 対照のオッズ比が理論上∞に発散します。これが完全分離です。この状態ではMLE（最尤推定）が正しく推定できません。
ステップ2：通常のロジスティック回帰（glm）が破綻する
まず普通のglm(family = binomial)を当てはめてみます。
f0 <- glm(ae ~ group, family = binomial, data = d1)
summary(f0)$coefficients
exp(coef(f0)) # オッズ比
exp(confint.default(f0)["groupactive", ]) # Wald 95%CI
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -20.56607 2643.087 -0.007781077 0.9937917
> groupactive 18.87439 2643.088 0.007141040 0.9943023
>
> (Intercept) groupactive
> 1.170226e-09 1.574144e+08
>
> 2.5 % 97.5 %
> 0 Inf
実薬群の係数は18.87、標準誤差は2643という異常な大きさです。z値はほぼ0、p値は0.994でまったく有意になりません。オッズ比は1.57×10⁸（1億超）、Wald信頼区間は(0, ∞)と、いずれも解釈不能な値になっています。
ここで重要なのは、Rがエラーも警告も出さずにこの無意味な結果を返したという点です。この例では”fitted probabilities numerically 0 or 1 occurred”の警告すら表示されません。完全分離は「巨大な標準誤差・p値≒1・桁外れのオッズ比・信頼区間が(0, ∞)」という症状から、自分で気づく必要があります。
ステップ3：Firth（logistf）で推定する
glm()をlogistf()に置き換えるだけです。
library(logistf)
f1 <- logistf(ae ~ group, data = d1)
summary(f1)
> coef se(coef) lower 0.95 upper 0.95 Chisq p method
> (Intercept) -4.510860 1.421963 -9.3497907 -2.569464 58.249234 2.309264e-14 2
> groupactive 2.875104 1.476918 0.7163346 7.752860 7.933699 4.852251e-03 2
>
> Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
> Likelihood ratio test=7.933699 on 1 df, p=0.004852251, n=90
実薬群の係数は2.875と有限の値に収まりました。信頼区間・p値はmethod=2、すなわちペナルティ付きプロファイル尤度で計算されています。尤度比検定はχ²=7.93（df=1）、p=0.00485で有意です。glmでは検出できなかった群間差を、Firthは適切に検出できています。
ステップ4：オッズ比と95%信頼区間（Firth）
係数を指数変換してオッズ比の形にします。
data.frame(
OR = round(exp(coef(f1)), 3),
L95 = round(exp(f1$ci.lower), 3),
U95 = round(exp(f1$ci.upper), 3),
p = signif(f1$prob, 3)
)
> OR L95 U95 p
> (Intercept) 0.011 0.000 0.077 2.31e-14
> groupactive 17.727 2.047 2328.221 4.85e-03
実薬群のオッズ比は17.73、95%信頼区間は(2.05, 2328.22)です。区間の幅は非常に広い（イベント数が少ないため当然です）ものの、下限が2.05と1を上回るため「実薬群でAEのオッズが有意に高い」と結論できます。glmが返した1.57×10⁸のような無意味な値ではなく、解釈可能な有限の推定値になった点が決定的な違いです。
多変量・希少イベントへの応用
Firthの方法は2×2の完全分離だけでなく、多変量モデルや「分離はしていないが希少イベント」の状況でも役立ちます。ここでは、希少AE（全160例中11件＝6.9%）を治療群・バイオマーカー高/低・年齢で説明する多変量モデルを考えます。バイオマーカー低群ではAEがほぼ発生しておらず（1/115）、準分離に近い状況です。
set.seed(2024)
n <- 160
biomarker <- factor(sample(c("low","high"), n, TRUE, prob = c(0.7,0.3)),
levels = c("low","high"))
group <- factor(sample(c("control","active"), n, TRUE),
levels = c("control","active"))
age <- round(rnorm(n, 62, 9))
lp <- -5.0 + 1.2*(group=="active") + 3.2*(biomarker=="high") + 0.03*(age-62)
ae <- rbinom(n, 1, plogis(lp))
d2 <- data.frame(ae, group, biomarker, age)
table(biomarker, ae) # 低群はイベントがほぼ無い
> total events: 11 / 160 (6.9%)
> ae
> 0 1
> low 114 1
> high 35 10
通常のglm（MLE）とlogistf（Firth）を並べて比較します。
g0 <- glm(ae ~ group + biomarker + age, family = binomial, data = d2) # MLE
g1 <- logistf(ae ~ group + biomarker + age, data = d2) # Firth
> term OR_MLE OR_Firth Firth_L95 Firth_U95 Firth_p
> (Intercept) 0.000 0.001 0.000 0.198 1.07e-02
> groupactive 2.568 2.372 0.654 9.530 1.90e-01
> biomarkerhigh 33.371 21.567 4.762 205.809 1.22e-05
> age 1.041 1.038 0.954 1.130 3.77e-01
| 項目 | OR（MLE） | OR（Firth） | Firth 95%CI下限 | 上限 | Firth p |
|---|---|---|---|---|---|
| (Intercept) | 0.000 | 0.001 | 0.000 | 0.198 | 0.0107 |
| group（active vs control） | 2.568 | 2.372 | 0.654 | 9.530 | 0.190 |
| biomarker（high vs low） | 33.371 | 21.567 | 4.762 | 205.809 | 1.22×10⁻⁵ |
| age | 1.041 | 1.038 | 0.954 | 1.130 | 0.377 |
バイオマーカー高のオッズ比は、通常のMLEでは33.37と過大ですが、Firthでは21.57（95%CI 4.76–205.81, p=1.22×10⁻⁵）に縮小しました。これは希少イベントで生じる小標本バイアスをFirthが補正した結果です。治療群（p=0.19）・年齢（p=0.38）は有意でなく、この点はMLEとFirthで結論が一致しています。完全分離が起きていなくても、イベント数が少ないときはFirthがオッズ比の過大評価を抑える点が実務上重要です。
なお、分離が無く十分なイベント数がある場合、FirthはMLEにほぼ一致します。実測では、n=400・173イベントというイベント数が十分な設定で、係数はMLE（切片 −0.373、x 0.984）とFirth（切片 −0.371、x 0.974）でほとんど同一の値になりました。この事実から、Firthは分離の有無にかかわらず既定手法として安全に運用できるといえます。
実務でのポイント
・小規模試験・希少AE・稀な曝露カテゴリでロジスティック回帰が分離したら、まずFirth（logistf）を検討します。
・分離の見分け方：オッズ比が桁外れ・標準誤差が数百〜数千・p値≒1・信頼区間が(0, ∞)。エラーが出なくても疑うことが大切です。
・信頼区間・検定はWaldでなく、ペナルティ付きプロファイル尤度（logistfの既定）を使います。Wald法は分離下では信頼できません。
・分離が無く十分なイベント数があるときはMLEにほぼ一致するため、既定手法として安全に使えます。
・代替手法として正確ロジスティック回帰（exact logistic）、罰則化（Lasso/Ridge、glmnet）、ベイズ（弱情報事前分布）があります。ただしFirthは追加設定が少なく実装が容易です。
・SASでは
PROC LOGISTICのFIRTHオプション、Rではlogistfパッケージで実装できます。参考書籍
Firthのロジスティック回帰とその周辺を深く理解するために役立つ書籍を3冊紹介します。
関連記事
Firthのロジスティック回帰とあわせて読むと理解が深まる記事を紹介します。
オッズ比・相対リスク・ハザード比の違い
本記事で頻出するオッズ比の解釈を、相対リスク・ハザード比と比較しながら整理した記事です。
順序ロジスティック回帰（比例オッズモデル）とは
2値ではなく順序カテゴリを扱う、ロジスティック回帰の姉妹モデルです。
ポアソン回帰とは
同じGLM一族の手法で、件数・発生率データを扱う際に用います。希少事象の解析という点でFirthと関心が重なります。
p値とは何か？
Firthの出力にも登場するp値（プロファイル尤度検定）を正しく理解するための基礎記事です。
まとめ
本記事では、ロジスティック回帰で完全分離・準完全分離が起きるとMLEが破綻する現象を、実際にRのglm()で再現しました。対照群のイベントが0件という状況では、オッズ比が1.57×10⁸、標準誤差が2643という無意味な値が、エラーも警告もなく返ってくることを確認しました。
これに対してFirthのペナルティ付き最尤法（logistfパッケージ）は、同じデータに対して有限で解釈可能な推定値（オッズ比17.73、95%CI 2.05–2328.22）を返し、群間差を正しく検出できました。さらに多変量モデルでは、完全分離が起きていない場合でも、希少イベントによる小標本バイアスをFirthが縮小してくれることを確認しています。一方で分離もなくイベント数も十分なときはMLEとほぼ一致するため、Firthは特別な場面だけでなく既定の手法として安全に運用できます。
小規模な臨床試験や希少な有害事象を扱う生物統計の実務では、こうした「MLEが黙って壊れる」場面に遭遇することが少なくありません。今回取り上げたFirthの方法とペナルティ付きプロファイル尤度による信頼区間・検定は、その対処法として押さえておく価値があると思います。あわせて、ポアソン回帰とはの記事もご覧いただくと、希少事象データの扱い方への理解がさらに広がるかと思います。
