Firth の方法を用いた場合のモデル選択規準について

この文章について

本ページは、Nagashima and Sato (Stat Med 2017) を概説したものです。
詳細な内容は論文をご参照ください。

はじめに

ロジスティック回帰モデルの推定において、すべての $\mathbf{x}_i$ について、 \[ \begin{cases} \beta' \mathbf{x}_i \geq 0, & \mathrm{event} = 1 \\ \beta' \mathbf{x}_i \leq 0, & \mathrm{event} = 0 \end{cases} \] または、 \[ \begin{cases} \beta' \mathbf{x}_i \leq 0, & \mathrm{event} = 1 \\ \beta' \mathbf{x}_i \geq 0, & \mathrm{event} = 0 \end{cases} \] をみたすとき、(凖)完全分離 ((quasi-)complete separation) の問題が起こることが知られている (Albert and Anderson, 1984)。 ただし、$\mathbf{x}_i$ は個体 $i$ ($=1,\ldots,n$) の $p$ 次元共変量ベクトル、$\beta=(\beta_1,\ldots,\beta_p)$ は $p$ 次元パラメータベクトルとする。 また、Cox 回帰モデルにおいても、同様に単調尤度 (monotone likelihood) の問題が起こりうる。 この問題が起こった場合、尤度関数があるパラメータの軸方向について極限値を持ってしまうため、最尤推定量が一意に定まらなくなり、標準誤差が無限大になる。 共変量の分布のアンバランスなデータ、強力な効果を持つ共変量が存在するデータ、多数の打ち切りが存在 (Cox 回帰の場合) するデータなどで、この問題が起こりうる事が知られている。

ここで、Cox 回帰での事例を説明する。

	require("survival")
	require("ggplot2")
	time  <- c(1, 2, 3)
	event <- c(1, 1, 1)
	x     <- c(1, 1, 0)
	sim   <- data.frame(time, event, x)
	logL  <- 1:100
	beta  <- 1:100/10
	for (i in 1:100) {
	res     <- coxph(sim, formula = Surv(time, event) ~ x,
						iter.max = 0, init = beta[i])
	logL[i] <- res$loglik
	}
	ggplot(data.frame(beta, logL), aes(beta, logL)) +
	geom_point() +
	theme_classic(base_size = 20)
	

このデータの尤度関数は \[ L(\beta_1) = \frac{e^{\beta_1}}{2 e^{\beta_1} + 1} \cdot \frac{e^{\beta_1}}{e^{\beta_1} + 1} \cdot \frac{1}{1} = \frac{1}{2 + 1/e^{\beta_1}} \cdot \frac{1}{1 + 1/e^{\beta_1}}, \] であり、 \[ L(\beta_1) \to \frac{1}{2}, \quad \mathrm{as} \,\, \beta_1 \to \infty, \] となる ($\beta_1$ が正の無限大のときに極限値を持つ)。
対数尤度関数を図示すると以下となる (関数値は $\log(0.5)=-0.6931472\cdots$ に近づく)。

ちなみに、Cox回帰の場合を一般化すると、 \[ L(\beta)= \prod_{i=1}^n \left[ \frac{e^{\beta' \mathbf{x}_i}}{\sum_{k \in R(t_i)} e^{\beta' \mathbf{x}_k}} \right]^{\delta_i}= \prod_{i=1}^n \left[ \frac{1}{\sum_{k \in R(t_i)} e^{\beta' (\mathbf{x}_k - \mathbf{x}_i)}} \right]^{\delta_i}, \] となる。 $j$ 番目のパラメータに着目し、$\beta_j$ 以外を $\hat{\beta}$ なる定数に全て置き換えたとき、尤度関数は \[ \prod_{i=1}^n \left[ \frac{1}{\sum_{k \in R(t_i)} e^{\beta_j (x_{kj} - x_{ij})}} \right]^{\delta_i}, \] と書け、イベント発現あり ($\delta_i=1$) の全ての時点について、 \[ x_{kj} - x_{ij} \geq 0, \] または、 \[ x_{kj} - x_{ij} \leq 0, \] を満たしている時に、単調尤度の問題が生じる事が分かる。 ただし、添え字 $k \in R(t_i)$ である。
例えば、ある離散共変量について、特定の層内のイベント数が $0$ 個の場合には、すべての $i$ について $x_{kj} = x_{ij}$ が成り立ち、単調尤度の問題が起こる。 また、Cox 回帰の単調尤度の問題は、ロジスティック回帰の場合と若干異なり、イベントの順序にも関連して起こることが分かる。

単調尤度・(凖)完全分離の問題を解決する方法に焦点を当てた代表的な研究として Heinze and Schemper (2001, 2002) があり、Firth の方法 (Firth, 1993) が単調尤度・(凖)完全分離のもとでも良い性能を持つことが示されている。 Firth の方法は元々最尤推定量の漸近バイアスを補正する目的で提案され、この方法は対数尤度関数にペナルティ項を加えることで調整を行うものである。 Firth 先生のブログの記事 (Firth, 2023) では、Firth の方法は元々 separation の問題 (原典では最尤推定量の existance の問題として言及されている) をペナルティ付き推定法により解決することをモチベーションとした研究であったが、最終的には漸近バイアスを補正する形の推定方程式が導出され、 一般化された結果が得られてあの形にまとまったと紹介されている (Firth の論文は 1993/3/1 に出版され、昨日でめでたく 30 周年を迎えました)。 Firth の方法のペナルティ項は観察情報行列の行列式の対数をとって $1/2$ 倍した $0.5\log|\mathbf{I}(\beta)|$ であり、$\mathbf{I}(\beta)$ が正定値であれば極大値を持つ関数である。 $\mathbf{I}(\beta)$ は、全イベント時点を通じて各時点のリスクセットにおいて線形従属な共変量ベクトルが存在する場合に半正定値になる (Fleming and Harrington, 1991)。 上述の通り、単調尤度・(凖)完全分離は結果変数に関する条件や、イベントの順序によっても生じるため、必ずしも $\mathbf{I}(\beta)$ は半正定値にならない。 したがって、単調尤度・(凖)完全分離の状況でもペナルティ項が極大値を持つ場合があるため、このときは推定値が得られることが分かる。 漸近バイアス補正を目的としたペナルティ項なので、ペナルティ項自身も最尤推定値の近くで極大値をとる。

Firth の方法で、単調尤度・(凖)完全分離の問題を回避することは可能であるが、この状況下でのモデル選択規準は検討されてこなかった。 現状では SAS/STAT の LOGISTIC Procedure と PHREG Procedure では AIC 型・ BIC 型のモデル選択規準 (以下ではこれらを $\mathrm{AIC}^*$, $\mathrm{BIC}^*$ と呼ぶ) が利用できる。 しかしながら、$\mathrm{AIC}^*$ を用いてモデル選択を行うと、$n \to \infty$ においてパラメータ数が最大のモデルが選ばれる事が分かった。 また、$\mathrm{BIC}^*$ も回帰パラメータの数に比例する負のペナルティ項を含むことが分かった。 本研究では、Cox 回帰に着目し単調尤度の状況で妥当なモデル選択規準を提案することを目的とした。 なお、ロジスティック回帰の場合も同様の議論が可能である。

方法

Cox 回帰モデルと Firth の方法

本稿では、Andersen and Gill (1982) の counting process による定式化と正則条件を仮定した Cox 回帰モデルを議論する。 $\mathbf{N}=(N_1,\ldots,N_i,\ldots,N_n)'$ を $n$ 要素の counting process とし、$N_i(t)$ はスケーリングした時間 $t \in [0, 1]$ においてイベント数をカウントする $0$ または $1$ の値をとる関数とする。 $N_i(t)$ がランダムな intensity process $h_i(t)=Y_i(t) h_0(t)\exp\{ \beta_0' \mathbf{Z}_i(t) \}$ を持つとする。 ただし、$h_0(t)$ はベースラインハザード関数、$Y_i(t)$ は predictable process であり $i$ 番目の個体が時点 $t$ でリスク集合に含まれる場合は $1$ そうでない場合は $0$ をとるとし、$\beta_0=(\beta_1,\ldots,\beta_j,\ldots,\beta_p)'$ を $p$ 次元の回帰パラメータの真値のベクトルとし、$\mathbf{Z}_i(t)=(Z_{i1},\ldots,Z_{ip})$ を $p$ 次元の共変量ベクトルで predictable process であるとする。 対数部分尤度は、 \[ l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \beta)=\sum_{i=1}^n \int_0^1 \beta' \mathbf{Z}_i(x) - \log\{ n S^{(0)}(\beta,x) \} \,\mathrm{d}N_i(x), \] スコアベクトルは、 \[ \mathbf{U}(\beta)=\frac{\partial l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \beta)}{\partial \beta}= \sum_{i=1}^n \int_0^1 \left\{ \mathbf{Z}_i(x) - \frac{\mathbf{S}^{(1)}(\beta,x)}{S^{(0)}(\beta,x)} \right\} \,\mathrm{d}N_i(x), \] 観測情報行列は、 \[ \mathbf{I}(\beta)=-\frac{\partial^2 l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \beta)}{\partial \beta \partial \beta'}= \sum_{i=1}^n \int_0^1 \left[ \frac{\mathbf{S}^{(2)}(\beta,x)}{S^{(0)}(\beta,x)} - \left\{ \frac{\mathbf{S}^{(1)}(\beta,x)}{S^{(0)}(\beta,x)} \right\}^{\otimes 2} \right] \,\mathrm{d}N_i(x), \] である。 ただし、$\mathbf{Y}=(Y_1,\ldots,Y_n)'$、$\mathbf{Z}=(\mathbf{Z}_1,\ldots,\mathbf{Z}_n)'$、$S^{(0)}$ はスカラー、$\mathbf{S}^{(1)}$ はベクトル、$\mathbf{S}^{(2)}$ は行列、$\mathbf{S}^{(k)}(\beta,t)=n^{-1}\sum_{i=1}^n Y_i(t) \mathbf{Z}_i(t)^{\otimes k} \exp\{ \beta' \mathbf{Z}_i(t) \}$ ($k=0,1,2$)、ベクトル $\mathbf{b}$ について $\mathbf{b}^{\otimes 0}=1$, $\mathbf{b}^{\otimes 1}=\mathbf{b}$, $\mathbf{b}^{\otimes 2}=\mathbf{b}\mathbf{b}'$ とする。 通常の Cox 回帰の回帰パラメータの最大部分尤度推定量 $\hat{\beta}_{Cox}$ は推定方程式 $\mathbf{U}(\beta)=\mathbf{0}$ の解である。

Firth の方法は、ペナルティ付き対数部分尤度関数 $l^*(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \beta)=l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \beta) + 0.5\log|\mathbf{I}(\beta)|$、およびそれに対応する修正スコア関数 $\mathbf{U}^*(\beta)=\mathbf{U}(\beta)+\mathbf{a}(\beta)$ を用いている。 ただし、$\mathbf{a}(\beta)=(a_1(\beta),\ldots,a_p(\beta))'$、$a_j(\beta)=\mathrm{tr} \left[ \left\{ \mathbf{I}(\beta) \right\}^{-1} \left\{ \partial \mathbf{I}(\beta)/\partial \beta_j \right\} \right]$ である。 バイアス調整した回帰パラメータの最大部分尤度推定量 $\hat{\beta}$ は推定方程式 $\mathbf{U}^*(\beta)=\mathbf{0}$ の解であり、これは $\hat{\beta}_{Cox}$ と異なる。

$\mathrm{AIC}^*$, $\mathrm{BIC}^*$ の問題点

SAS ではペナルティ付き対数部分尤度を用いて、AIC 型のモデル選択規準 $\mathrm{AIC}^*=-2l^*(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta})+2p$ を用いてる。 しかし、 \[ \mathrm{AIC}^*= -2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta}) + (2 - \log n)p -\log\{n^{-p}|\mathbf{I}(\hat{\beta})|\}, \] と変形できてしまう。 第一項 $-2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta})$ は、基本的にはオーバーフィッティングの観点からはパラメータ数が多いほど小さい値をとる。 第二項 $(2 - \log n)p$ は、$n \geq 8$ ならば $\lt 0$ であるため、この条件下ではパラメータ数が多いほど小さい値をとる項であり、かつこの項は $\log n$ の大きさに依存する。 第三項 $-\log\{n^{-p}|\mathbf{I}(\hat{\beta})|\}$ は、$n \to \infty$ において $n^{-p}|\mathbf{I}(\hat{\beta})| \to_P |\Sigma(\beta_0)|$ より定数に近づく。 したがって、$n$ が大きいときはモデル中の回帰パラメータ数が多いものが選択されやすくなる事が分かる (モデル適合度は無視され、望ましくない性質を持つ)。

影響度は大きく異なるが、 \[ \mathrm{BIC}^*=-2l^*(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta})+p\log d=-2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta})+p\log(1-c)-\log\{n^{-p}|\mathbf{I}(\hat{\beta})|\}, \] もパラメータ数に対して負のペナルティ項を持つ。 ただし、$d$ はイベント数、$c=1-d/n \in (0, 1]$ である。 こちらも望ましくない項が含まれているが、定数オーダーのペナルティ項である ($n$ の大きさには依存しない)。

提案法 ($\mathrm{AICF}$, $\mathrm{BICF}$)

本研究では、リスク関数 $\mathrm{E}_{\mathbf{N}} \mathrm{E}_{\tilde{\mathbf{N}}} [-2l(\tilde{\mathbf{N}}, \tilde{\mathbf{Y}}, \tilde{\mathbf{Z}}; \hat{\beta})] $ の推定量に基づいて、$\hat{\beta}$ をもちいた場合の $\mathrm{AIC}$ 型の規準 $\mathrm{AICF}$ を提案した。 ただし、$\mathrm{E}_{\mathbf{N}}[-2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \beta)]=\int_{\Omega} -2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \beta) \mathrm{d}P$ を表わし、$(\tilde{\mathbf{N}}, \tilde{\mathbf{Y}}, \tilde{\mathbf{Z}})$ は現在得られているデータとは別の新しいの観測値の組を表わす。 ほぼ予想された結果ではあるが、リスク関数に基づく導出の結果 \[ \mathrm{AICF}=-2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta})+2p, \] が導かれ、これは通常の Cox 回帰の場合の $\mathrm{AIC}=-2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta}_{Cox})+2p$ において回帰パラメータの推定量を単純にバイアス修正推定量に置き換えたものに等しい事が分かる。 ペナルティ付きの対数部分尤度関数をベースに推定していても、モデル選択規準には罰則項を含めてはならず、同種の問題を扱う場合にも注意が必要である。

同様に、 \[ \mathrm{BICF}=-2l(\mathbf{N}, \mathbf{Y}, \mathbf{Z}; \hat{\beta})+p \log d, \] も提案した。

シミュレーション

詳細は割愛する。 結果の概要を以下に列挙した。

• 「$\mathrm{AIC}^*$, $\mathrm{BIC}^*$の問題点」の節で述べたとおり、$n$ が大きいか、$d$ が大きいとき、$\mathrm{AIC}^*$ はパラメータ数最大のモデルを $100\%$ 選択した。 $\mathrm{AIC}^*$ はバイアスがあり、リスク関数の推定量を構成していないと考えられる。
• $\mathrm{AICF}$ はバイアスがなく、リスク関数の推定量を構成していると考えられる。 $\mathrm{AICF}$ に基づいて選択されたモデルの推定値を用いて別の新しいの観測値の組におけるデビアンス ($-2l(\tilde{\mathbf{N}}, \tilde{\mathbf{Y}}, \tilde{\mathbf{Z}}; \hat{\beta})$) の平均値を求めると、$\mathrm{AIC}^*$ で同様の計算をした場合よりも値が小さくなった (予測性能は $\mathrm{AIC}^*$ より良い)。
• $\mathrm{BIC}^*$ は条件にもよると思われるが、モデル選択の一致性を持たないと考えられる。
• $\mathrm{BICF}$ はモデル選択の一致性を持つと考えられる。ただし、$d$ が小さいときはあまり上手くいかない。

データ解析

JLGK0901 study (Yamamoto et al., 2014) のデータを用いて、$\mathrm{AIC}^*$ および $\mathrm{AICF}$ を用いてモデル選択した結果を検討する。 他の解析結果の詳細は、論文をご参照のこと。

解析に含めた共変量は $age$ (年齢), $sex$ (性別), $kps$ (Karnofsky performance status), $ntumor$ (腫瘍数), $diameter$ (腫瘍の最大径), $volume$ (腫瘍の体積), $ptumor$ (原発腫瘍の種類), $status$ (脳以外の疾病状態), $neuro$ (神経症状の有無) の $9$ 変数であり、髄膜播種を起こすまでの時間を結果変数とした。 本データは $n = 1073$、$d = 145$ である。 以下に上位 $5$ 番目までのモデル選択結果と、各規準での $1$ 番のモデルにおける解析結果を示す。

Table 1 および Table 3 より、$\mathrm{AIC}^*$ は明らかにパラメータ数が多いモデルを選び、$P$ 値が $1.00$ に近い変数 ($age$, $kps$, $status$ など) まで選んでしまっている事が分かる。 Table 2 および Table 4 から、$\mathrm{AICF}$ の方がパラメータ数が少ないモデルを選択しており、Table 3 で $P$ 値が $1.00$ に近かった変数は選択されておらず、比較的 $P$ 値が小さい変数が選ばれた。 よって、$\mathrm{AICF}$ によって選ばれたモデルは、$\mathrm{AIC}^*$ で選ばれていたような明らかに予測に寄与しなそうな変数を含まないと言える。

まとめ

単調尤度・(凖)完全分離の問題の解決方法として Firth の方法があるが、モデル選択規準については研究されてこなかった。 本研究では、この条件下でのモデル選択規準について検討した。 SAS で計算される $\mathrm{AIC}^*$ や $\mathrm{BIC}^*$ には問題点があり、これを解決する方法として $\mathrm{AICF}$ と $\mathrm{BICF}$ を提案した。 シミュレーションによる性能評価や、実データの解析結果から、既存法の問題点はおおむね解決できたと考えられる。

謝辞

本研究においてJLGK0901 studyのデータ使用をご快諾いただいた、勝田病院 水戸ガンマハウスの山本昌昭教授、および築地神経科クリニックの芹澤徹先生に感謝申し上げます。
本研究はJSPS科研費 26870099, 16K16014 の助成を受けたものです。

文献

1. Nagashima K, Sato Y. Information criteria for Firth's penalized partial likelihood approach in Cox regression models. Statistics in Medicine 2017;36(21):3422–3436. DOI: 10.1002/sim.7368. [Postprint]
2. Albert A, Anderson A. On the existence of maximum likelihood estimates in logistic regression models. Biometrika 1984; 71(1): 1–10. DOI: 10.1093/biomet/71.1.1.
3. Firth D. Bias reduction of maximum likelihood estimates. Biometrika 1993; 80(1): 27–38. DOI: 10.1093/biomet/80.1.27.
4. Heinze G, Schemper M. A solution to the problem of monotone likelihood in Cox regression. Biometrics 2001; 57(1): 114–119. DOI: 10.1111/j.0006-341X.2001.00114.x.
5. Heinze G, Schemper M. A solution to the problem of separation in logistic regression. Statistics in Medicine 2002; 21(16): 2409–2419. DOI: 10.1002/sim.1047.
6. Fleming TR, Harrington DP. Counting Processes and Survival Analysis. New York: Wiley, 1991.
7. Andersen PK, Gill RD. Cox's regression model for counting processes: a large sample study. The Annals of Statistics 1982; 10(4): 1100–1120. DOI: 10.1214/aos/1176345976.
8. Yamamoto M, Serizawa T, Shuto T, Akabane A, Higuchi Y, Kawagishi J, Yamanaka K, Sato Y, Jokura H, Yomo S, Nagano O, Kenai H, Moriki A, Suzuki S, Kida Y, Iwai Y, Hayashi M, Onishi H, Gondo M, Sato M, Akimitsu T, Kubo K, Kikuchi Y, Shibasaki T, Goto T, Takanashi M, Mori Y, Takakura K, Saeki N, Kunieda E, Aoyama H, Momoshima S, Tsuchiya K. Stereotactic radiosurgery for patients with multiple brain metastases (JLGK0901): a multi-institutional prospective observational study. The Lancet Oncology 2014; 15(4): 387–395. DOI: 10.1016/S1470-2045(14)70061-0.
9. Firth D. My best known work, fading citations and a bit of history. [2023-01-05] https://davidfirth.github.io/blog/2023/01/05/f93-citations-and-history/.

履歴

• 2023/03/02 Firth の方法の元々のモチベーションに関する記述を修正
• 2017/06/29 公開