11.6 パラメトリック生命表解析

(1) パラメトリックモデル

第4節で「ハザード関数λ(t)の具体的な姿を規定するのは困難な場合が多い」と説明しましたが、あえてハザード関数の具体的な姿を規定して生命表解析を行う手法があります。 その手法はパラメトリック生命表解析と呼ばれ、その手法で用いられるモデルはパラメトリックモデルと呼ばれます。 これに対して第1節と第2節で説明した生存率の計算方法と比較方法はノンパラメトリック生命表解析と呼ばれ、第4節で説明した比例ハザードモデルはセミパラメトリックモデルまたはセミノンパラメトリックモデルと呼ばれます。

最も単純なパラメトリックモデルは第3節で例として説明したハザード関数が常に一定と仮定するモデルです。 このモデルは死亡率関数f(t)つまり生存時間の分布が指数分布になるので指数分布モデルまたは標的モデルと呼ばれます。

このモデルでは死亡数と観察期間からλの最尤推定値を求めることができます。 そしてλは単位時間あたりの死亡率なので、その逆数は1人あたりの生存時間つまり平均生存時間になります。 このモデルではλは時間によらず一定ですが、λが時間によって変化してもt = 0〜∞のλの平均値つまり平均ハザードの逆数は平均生存時間になり、それは生存時間の分布つまり死亡関数f(t)の平均値に相当します。 生存時間の原理面ではハザードは重要な指標です。 しかしハザードをそのまま解釈するよりも、それを逆数にして平均生存時間にした方が解釈しやすくて実際的です。

例えば表11.1.1のデータに指数分布モデルを当てはめると次のようになります。 そして最尤法を利用して2群のλの比つまりハザード比が1かどうかの検定を行うことができます。 その手法をパラメトリック・ハザード比検定(PHR-test：parametric hazard ratio test)と名付けることにしました。 _(注1)(注2)

○A群

S(t) = exp(-0.0194t)
平均生存時間 =
50％生存時間(MST) =
○B群

S(t) = exp(-0.0533t)
平均生存時間 =
50％生存時間(MST) =
○パラメトリック・ハザード比検定
ハザード比(B群／A群)：
検定：χ_β1² = 3.523(p = 0.0605) ＜ χ²(1,0.05) = 3.841 … 有意水準5％で有意ではない
95％信頼区間　下限：HR_L = 0.9561　上限：HR_U = 7.9417

この方法で求めたハザード比2.756は第2節のコックス・マンテル検定で求めたハザード比3.697とは少し異なります。 この方法ではハザードが常に一定と仮定して群ごとのハザードを最尤法によって推定し、その最尤推定値の比からハザード比を計算しています。 それに対してコックス・マンテル検定では死亡時間という重要な情報を用いず、死亡例が発生するたびにその時点のハザードを群ごとに計算して、それらの平均的なハザード比を計算しています。 この計算方法の違いが両者のハザード比の違いの原因です。

参考までに第2節で行ったように表11.1.1の時間間隔を全て1にして累積生存率曲線を描き、これに上記の指数分布モデルを当てはめると図11.6.2のようになります。 第2節で説明したように、ノンパラメトリック手法を用いると図11.6.1も図1.6.2も全く同じ検定結果と推定結果になります。 しかしパラメトリック手法では次のように結果が変わります。 このことから死亡時間という重要な情報を用いるパラメトリック手法の方が正確かつ合理的であることがわかると思います。

図11.6.1のB群では35ヶ月に最後の1例が死亡しているので、カプラン・マイヤー法で計算した累積生存率曲線は35ヶ月で累積生存率がストンと落ちて0になっています。 しかし母集団は例数がもっと多い——理論的には無限例——ので、累積生存率が35ヶ月でいきなり0になるとは考えにくく、理論的生存関数のように35ヶ月以後も徐々に低くなっていくと考えられます。 このことから理論的生存関数の方がより普遍性があり、将来のことをある程度は予測することができるつまり外挿が可能であることがわかります。

検量線に例えればカプラン・マイヤー法で求めた累積生存率曲線は実際のデータを折れ線で結んだ検量線に相当し、理論的生存曲線はデータを直線で回帰した理論的検量線に相当します。 折れ線の検量線よりも理論的検量線の方が精度が高く普遍性があるのと同様に、累積生存率曲線よりも理論的生存曲線の方が精度が高く普遍性があります。

疫学分野では人時間に基づく罹患率(person-time incidence rate)IR_P-Tという指標が利用されます。 これは罹患密度(incidence density)または罹患力(force of morbidity)または単に罹患率(incidence rate)とも呼ばれ、次のように定義されています。

人時間単位としてよく使われるのは人年(person-years)であり、kは5が使われます。 そして例えば「日本における2000〜2010年のガンの人時間に基づく罹患率は10万人年あたり300人」などと表現されます。 この定義式から指数分布モデルにおけるハザードλは人時間に基づく死亡率に相当することがわかると思います。 したがって人時間に基づく罹患率は、特定の期間中は罹患率が一定であり、その間の累積罹患率曲線を指数関数によって近似できるという暗黙の前提を含んだ指標であることがわかります。

(2) 各種のパラメトリックモデル

指数分布モデル以外にも色々なパラメトリックモデルが提唱されているので、代表的なものを紹介しましょう。

(a) ワイブル分布モデル

標的モデルが直列にa個並んだモデルであり、直列モデルとも呼ばれます。 このモデルはf(t)がワイブル分布(Weibull distribution)になり、a = 1の時は指数分布モデルになります。 例えば生命維持にとって大切な臓器がa個あり、そのどれか1つでも損傷すると死亡する場合はこのモデルが当てはまります。 _(注3)

(b) ガンマ分布モデル

標的モデルが並列にb個並んだモデルであり、並列モデルとも呼ばれます。 このモデルはf(t)がガンマ分布(Gamma distribution)になり、b = 1の時は指数分布モデルになります。 例えば生命維持にとって大切な臓器がb個あり、それら全てが損傷すると死亡する場合はこのモデルが当てはまります。 (→付録1　各種の確率分布)

　　
　　 … ガンマ分布　(0 ＜ λ ＜ 1、0 ＜ b)
Γ(b)：ガンマ関数　　Γ(b;λt)：不完全ガンマ関数

(c) 対数正規分布モデル

f(t)が対数正規分布になるモデルです。 例えば理論的には一定時間後に必ず死亡するものの、その時間には誤差があり、その誤差が近似的に対数正規分布する、つまり時間を対数変換した時に誤差が正規分布するような場合はこのモデルが当てはまります。

　　　　
… 対数正規分布
… 標準正規分布の確率分布関数
… 標準正規分布の確率密度関数

これらのモデルを実際の累積生存率曲線に当てはめた時、どのモデルが最もフィットするかを検討することによって疾患による死亡の生理的メカニズムを推測する時の参考になります。 そしてそれにより色々な治療方法の効果を比較するだけでなく特徴を比較することができ、単なる勝ち負けではなく使い分けの検討が可能になります。

単純な勝ち負けの検討の時または生理的メカニズムが不明の時は、指数分布モデルが最もよく用いられます。 このモデルはパラメーターがλだけであり、その推定値を比較的簡単に計算することができるからです。 図11.6.1を見ればわかるように、指数分布モデルは実際の累積生存率曲線にかなりうまく近似します。 医学分野で用いられる普通の回帰直線はこれほどうまくは近似しないでしょう。

(3) 境界内平均生存時間

生存関数S(t)は生存時間tの発生確率を表す関数なので、tについて積分すると平均生存時間になります。 例えば指数分布モデルの生存関数をt = 0〜∞について積分すると、次のように平均生存時間になります。 _(注5)

ところが実際のデータは全例が死亡するまで観察せず、たいていは全例が死亡する前に観察を打ち切ることが多いと思います。 そのため生存関数をt = ∞まで外挿して求めた平均生存時間は、どうしても信頼性が低くなってしまいます。 そこで生存関数を実際のデータが存在する時間t = τまで積分し、その時間までの平均生存時間を求めれば信頼性が低くならないと考えられます。 これを境界内平均生存時間(RMST：Restricted Mean Survival Time)といい、指数分布モデルでは次のようになります。

境界内平均生存時間は生存関数を積分した値なので、累積生存率曲線の曲線下面積AUC(Area Under the Curve)になります。 そのためカプラン・マイヤー法による累積生存率曲線の曲線下面積を計算することによって近似値を求められます。 例えば表11.1.1のデータについて指数分布モデルを用いた時と、カプラン・マイヤー法による累積生存率曲線を用いた時の境界内平均生存時間を求めると次のようになります。

上記のように、指数分布モデルによる境界内平均生存時間とカプラン・マイヤー法による境界内平均生存時間はまずまず近似しています。 指数分布モデルに限らず、各種のパラメトリックモデルの生存関数はカプラン・マイヤー法による累積生存率曲線の中央付近を通ります。 そのためカプラン・マイヤー法による境界内平均生存時間は、どんなパラメトリックモデルについてもまずまずの近似値になります。 したがって群によって分布モデルが異なる時でも、また分布モデルがわからない時でも実用的な近似値として用いることができます。

またA群ではt = 0から35ヶ月までの境界内平均生存時間と、t = 0から56ヶ月までの境界内平均生存時間の2種類を計算しました。 これはB群のデータが35ヶ月までしかないので、それと比較するためです。 このように2群の境界内平均生存時間を比較する時は同じ時間までの値を比較する必要があります。 そのため2群の最終観察時間が異なっている時は、最終観察時間が多い群のデータを少し無駄にしてしまいます。

そのような欠点はあるものの、境界内平均生存時間は群によって分布モデルが異なる時でも群の平均生存時間を比較できるという長所を持っています。 そこでこの長所に着目して、2群の累積生存率曲線が非平行の時つまり比例ハザード性が成り立たない時にハザード比の代替え指標として用いる時があります。

しかし境界内平均生存時間には少々難点があります。 (2)で説明したように、比例ハザード性が成り立たない時は単純な勝ち負けではなく薬剤の特徴や使い分けを検討することが大切です。 ところが境界内平均生存時間はその検討が難しいのです。

例えば表11.6.1のようなデータがあったとします。 これは表11.1.1のB群のデータを変更してC群のデータにしたものです。 そして図11.6.6に示したように、C群のデータにはワイブル分布モデルがうまく当てはまります。

このデータにパラメトリック・ハザード比検定とコックス・マンテル検定を適用し、境界内平均生存時間を求めると次のようになります。

上記のようにパラメトリック・ハザード比検定でもコックス・マンテル検定でもハザード比は1.5くらいであり、全体としてはC群の方がハザードが少し高い(生存率が少し低い)と解釈できます。 そして交互作用のχ²値が少し大きいので2群の累積生存率曲線が非平行で比例ハザード性が成り立っていない可能性があります。 これらのことは図11.6.6のグラフを見ると何となくわかると思います。

ただし例数が少ないので、どの検定結果も有意水準5％で有意ではありません。 第2節で説明したように交互作用の検定は2群の累積生存率曲線の非平行性の検定であり、比例ハザード性が成り立っているかどうかをチェックするためのものです。 しかしこれは有意性検定になりがちなので、検定結果が有意かどうかよりも検定統計量であるχ²値が大きいかどうかをチェックする方が実際的です。

また境界内平均生存時間はA群の方が少し大きく、ハザード比と似た傾向です。 ただし指数分布モデルによって求めた値よりもカプラン・マイヤー法によって求めた値の方が2群の差が小さくなっています。 カプラン・マイヤー法によって求めた値は近似値なので、この程度の誤差は致し方ないでしょう。

図11.6.6を見ると、2本の理論的生存率曲線は30ヶ月後あたりで交わっています。 この交点は指数分布モデルの生存関数とワイブル分布モデルの生存関数から理論的に求めることができ、正確には28ヶ月後になります。 _(注6)

このことからA群は指数分布モデルが当てはまり、指数関数的に生存率が低下していくのに対して、C群はワイブル分布モデルが当てはまり、最初のうちは生存率があまり低下しないものの、20ヶ月後くらいから急激に低下し、28ヶ月以後はA群よりも生存率が低くなることがわかります。 そしてA群とC群がそれぞれ薬剤Aと薬剤Cを用いていたのなら、28ヶ月までは薬剤Cを用い、28ヶ月以後は薬剤Aを用いると効果的な治療ができる可能性が高いと考えられます。

このように比例ハザード性が成り立たない時は、単純な勝ち負けではなく薬剤の特徴や使い分けを検討することが大切です。 その検討をするには累積生存率曲線のグラフを描いて特徴をチェックし、うまく当てはまるパラメトリックモデルを調べるのが有用です。 そしてこのような検討にはハザード比も境界内平均生存時間も適していません。 そもそも境界内平均生存時間は原理的にハザードに反比例するので、ハザード比と同じ特徴を持っているのは当然です。 したがって比例ハザード性が成り立たない時に境界内平均生存時間を代替え指標にするのは、あまり良い方法とはいえません。

(4) パラメトリック多変量生命表解析

比例ハザードモデルと同じように、パラメトリックモデルでもハザードに影響を与える因子を共変数にし、リンク関数を用いた重回帰型モデルを想定して多変量生命表解析を行うことができます。 例えば指数分布モデルでは次のような重回帰型モデルを想定することができます。

比例ハザードモデルの場合、ハザード関数は共変数によっても時間によっても形を変えず、ハザード関数と基準ハザード関数の比を対数変換した対数ハザード比と共変数の間に線形関係があると仮定しました。 この仮定は現実にはほとんど有り得ず、最尤法の適用方法も苦肉の策です。 でも、その代わりハザード関数の具体的な姿を規定する必要がないというメリットがあります。

それに対して指数分布モデルを利用した重回帰型モデルでは、ハザード関数が時間とは無関係に一定であり、対数ハザードと共変数の間に線形関係があると仮定します。 この仮定も現実にはほとんど有り得ません。 しかし数学的な取り扱いが簡単で最尤法を正確に適用することができ、結果が解釈しやすく、生存率の外挿が可能であるという大きなメリットがあります。

このモデルにおける各変数の偏回帰係数は、他の変数が一定という条件下で各変数が「1」増加した時に、対数ハザードがいくつ変化するかを表す値つまり対数ハザードの変化量になります。 したがって偏回帰係数を指数変換した値は、比例ハザードモデルと同じように補正ハザード比(調整ハザード比)になります。 そして最尤法による解が漸近的に正規分布するという性質を利用して、偏回帰係数が0かどうかの検定つまりハザード比が1かどうかの検定と推定を行うことができます。

x_j = 0の時の対数ハザード：η_j0 = ln(λ_j0)
x_j = 1の時の対数ハザード：η_j1 = ln(λ_j1)
対数ハザードの変化量：
∴ 　　HR_j：x_jのハザード比

例えば表11.3.1のデータに指数分布モデルを当てはめ、最尤法を利用して解を求めると次のようになります。 _(注2)

パラメトリックモデルでは生存関数S(t)の具体的な内容を求めることができるので、S(t)を利用して共変数が任意の値の時の理論的生存関数を求めることができます。 例えば重症度が軽症で、治療が無い時と有る時の理論的生存関数は次のようになり、それらをグラフ化すると図11.6.7のようになります。

図11.6.7には、比較のために図11.4.1と同じ比例ハザードモデルによる基準生存関数S₀(t)と理論的生存関数も描いてあります。 全ての共変数が平均値の時の理論的生存関数は共変数を無視した時の理論的生存関数になり、それは基準生存関数に指数分布モデルを当てはめた時の理論的生存関数になります。 上の計算結果と図11.4.1を見ると、指数分布モデルの結果と比例ハザードモデルの結果はよく似ていることがわかると思います。

多変量生命表解析に限らず生命表解析全般で、現在はパラメトリックモデルよりもセミパラメトリックモデルやノンパラメトリック生命表解析の方が多用されています。 しかしセミパラメトリックモデルやノンパラメトリック生命表解析は折れ線で描いた検量線のようなものであり、結果がデータに左右されやすくて普遍性がなく、理論的な取り扱いが困難で、実際のデータよりも先を外挿によって予測することはできません。

それに対してパラメトリック生命表解析は結果がデータに左右されにくくて普遍性があり、理論的な取り扱いが簡単で結果が解釈しやすく、しかもある程度は外挿が可能という特徴を持っています。 そしてパラメトリック生命表解析は、疾患による死亡の生理的なメカニズムを推測する時の参考にすることもできます。 さらにパラメトリック生命表解析は試験の必要例数を比較的簡単に計算できるので、正確な試験計画を立てることができます。 そのためパラメトリック生命表解析はもっと利用されてしかるべきです。 _(注4)

(注1)　指数分布モデルの生存関数に最尤法を適用し、λの最尤推定値を求めると次のようになります。 そしてこれらの値を用いて各種の検定や推定を行うことも可能です。

死亡例の確率：f(t_i) = λ(t_i)S(t_i) = λe^-λt_i 　　脱落例の確率：S(t_i) = e^-λt_i
全体の尤度：
対数尤度：
　　∴
　　∴ 　　I：フィッシャー(Fihser)情報量
※脱落例がある時は死亡例合計の代わりに、1例ごとの累積死亡確率の合計を用いた方が近似が良いという考え方もある。

○2群のハザードの差の検定＝パラメトリック・ハザード差検定と区間推定
δ_λ = λ₁ - λ₂ 　　V(δ_λ) = V(λ₁) + V(λ₂) 　　　(λ₁とλ₂は独立なので共分散 C(λ₁,λ₂) = 0)
＞ t(∞,α)の時、有意水準100α％で有意
100(1-α)％信頼区間：
→ 下限：δ_λL = δ_λ - t(∞,α)SE_δλ 　上限：δ_λU = δ_λ + t(∞,α)SE_δλ

脱落例がある時の近似分散は1例ごとに観測終了時での理論的累積死亡確率{1 - exp(-λt_i)}を求め、それを合計したものを死亡例数の代わりに使うという考え方に基づいたものです。 これは死亡例が全て無限時間後に発生し、脱落例は全てその前に発生した時に、死亡例合計を用いた最初の近似分散と一致します。 そして死亡例が無限時間より前に発生すると分散の分母を小さくし、脱落例が死亡例よりも後に発生すると分散の分母を大きくするので2種類の近似分散は似た値になります。 そこで(注2)で説明するパラメトリック・ハザード比検定との整合性を考慮して、ここでは死亡例合計を用いた近似分散を用いることにします。

指数分布の期待値はλの逆数になり、中央値は期待値にln(2)を掛けた値になります。 そのためλから平均生存時間と50％生存時間(MST)を求めることができます。 ノンパラメトリック生命表解析では全例が死亡していないと平均生存時間を求められず、半数以上が死亡していないと50％生存時間を求められません。 それに対してパラメトリック生命表解析では、死亡例が半数未満でもそれらの推定値を求めることができます。

平均生存時間： 　　
50％生存時間(MST)：

表11.1.1のA群について実際に計算すると次のようになります。 第1節の(注2)で求めたカプラン・マイヤー法によるMSTは35.11でしたから、この方法で求めたMSTとよく近似しています。

　　

(注2)　比例ハザードモデルと同様に、表11.4.2に指数分布モデルを当てはめると次のようになります。 (→11.4 比例ハザードモデル (注1))

表11.4.2 共変数がp個の一般的データ

共変数					観測期間	転帰
x11	…	x1j	…	x1p	t1	d1
:		:		:	:	:
xi1	…	xij	…	xip	ti	di
:		:		:	:	:
xn1	…	xnj	…	xnp	tn	dn

※d_i:死亡例は1、脱落例は0のダミー変数

指数分布モデル：λ(t) = λ　　S(t) = exp(-λt)

η = β₀１_n + β₁ｘ₁ + … + β_jｘ_j + … + β_pｘ_p + ε = Ｘβ + ε = ｙ + ε 　　ｙ = Ｘβ
　　 　　 　　
死亡例の確率：f(t_i|ｘ_i)=λ(t_i|ｘ_i)S(t_i|ｘ_i) 　　脱落例の確率：S(t_i|ｘ_i)
全体の尤度：
対数尤度：
※λ_i = exp(β'ｘ_i)

この場合はハザード関数が規定されているので比例ハザードモデルのように部分尤度関数という苦し紛れの手を用いる必要はありません。 この対数尤度関数にニュートン・ラプソン法を適用し、最尤解を求めると次のようになります。

Ｈ_k = ∇ｇ_k

ｂ_k+1 = ｂ_k - Ｈ_k^-1ｇ_k

比例ハザードモデルと同様に、偏回帰係数が0かどうかの検定つまりハザード比が1かどうかの検定と推定をワルドの検定によって行うことができます。

V(ｂ) = -Ｈ_k^-1 ≒ Ｉ_f^-1 　　E(-Ｈ_k) = Ｉ_f：情報行列 　　[Ｉ_f]_jj^-1：Ｉ_f^-1の第j対角要素
検定： ＞ χ²(1,α)の時、有意水準100α％で有意
推定：100(1-α)％信頼区間
→ 下限: 　上限：

また偏回帰係数が全て0の時の尤度つまり共変数がなくて定数項だけのモデルの尤度と、定数項と共変数がp個のモデルの尤度の比を利用した尤度比検定によって、共変数全体の回帰の検定を行うことができます。 定数項だけのモデルのハザードλ₀は(注1)で導いたように全死亡数を観察期間合計で割った値になります。 そのためλ₀を用いて定数項だけのモデルの尤度を簡単に計算することができます。

対数尤度差：
尤度比検定：-2(L(0) - L(β)) = χ_β² ＞ χ²(p,α)の時、有意水準100α％で有意

また飽和モデルの対数尤度とそれを利用した異質性の検定、そして擬似寄与率とAIC(赤池の情報量基準)を次のようにして計算することができます。 疑似寄与率とAICはモデルがデータにどの程度適合しているか検討する時の指標になります。 (→10.3 ロジスティック回帰分析の計算方法 (注2))

飽和モデルの対数尤度： 　　尤度の自由度 = (p + 1) - 1 = p
　　m：共変数ｘ_iが等しい個体の組数、1組はn_k ≧ 1例
対数尤度差： 　　尤度の自由度 = m - (p + 1) = m - p - 1
尤度比検定：D = -2(L(β) - L_f)=χ_LOF² ＞ χ²(m-p-1,α)の時、有意水準100α％で有意
D：デビアンス(deviance) … モデルとデータのズレの大きさを表す指標
擬似寄与率：
AIC(赤池の情報量基準) = -2L(β) + 2(p + 1) = 2{(p+1) - L(β)}

偏回帰係数の初期値は、比例ハザードモデルと同様に対数変換した生存時間を目的変数にした重回帰分析の偏回帰係数を求め、その符号を反対にしたものを用いることができます。

表11.3.1のデータに指数分布モデルを当てはめ、実際に計算してみましょう。

対数生存期間を目的変数にした重回帰分析の結果：y = 2.9135 + 0.77911x₁ - 0.60483x₂
y：ln(観察期間)　　x₁：治療(0：無 1：有) 　　x₂：重症度(0：症状無 1：軽症 2：重症)
ln(0.036246) = b₀₀ - 0.77911×0.485714 + 0.60483×1.01429
∴b₀₀ = -3.31742 + 0.77911×0.485714 - 0.60483×1.01429 = -3.55247
初期値：
　　

更新されたｂ₁を用いて同様の計算を繰り返すと、4回目で値が収束します。

　　 　　 　　
ｂ₅ = ｂ₄ - Ｈ₄^-1ｇ₄ ≒ ｂ₄
指数分布モデル：y = ln(λ) = -3.639 - 0.667x₁ + 0.708x₂ 　　S(t) = exp(-λt)
○偏回帰係数の検定と推定
・
　ハザード比：HR₁ = exp(-0.667) = 0.513
　ハザード比の95％信頼区間　下限：HR_1L = exp(-1.199) = 0.302 　上限：HR_1U = exp(-0.135) = 0.874
・
　ハザード比：HR₂ = exp(0.708) = 2.030
　ハザード比の95％信頼区間　下限：HR_2L = exp(0.369) = 1.447 　上限：HR_2U = exp(1.047) = 2.849
○尤度比検定
　このモデルの尤度：L(β₀,β₁,β₂) = -231.78
　定数項だけのモデルの尤度：L(β₀) = 56×{ln(0.036246) - 1} = 56×(-4.317426) = -241.7759
　飽和のモデルの尤度：L_f = -231.183
　回帰：χ_β² = -2×(-241.7759 + 231.78) = 19.991(p = 0.00005) ＞ χ²(2,0.05) = 5.991
　ズレ：χ_LOF² = -2×(-231.78 + 231.183) = 1.194(p = 0.754444) ＜ χ²(3,0.05) = 7.815
擬似寄与率：
AIC(赤池の情報量基準) = 469.561

このモデルでもワルドの検定に用いるχ²値をスコア統計量にして、比例ハザードモデルと全く同じ手順で変数選択を行うことができます。 計算過程の説明は省きますが、同じデータについて変数選択を行うと2つの共変数がどちらも選択されて上と同じ結果になります。 (→11.5 変数選択法 (注1))

指数分布モデルにおいて説明変数が1つだけで、しかもそれが0または1の値を取るダミー変数の時は次のようになります。

指数分布モデル：η_i = ln{λ(t|ｘ_i)} = β₀ + β₁x_i1 + ε_i = y_i + ε_i　(x_i1：群1の時は0、群2の時は1になるダミー変数)
対数尤度：

∴
n₁：群1の例数　　n₂：群2の例数　　d：全死亡数 　　：群1の観測期間合計 　　：群2の観測期間合計


∴ 　　
d₁：群1の死亡数　　d₂：群2の死亡数 　　λ₁：群1のハザード　　λ₂：群2のハザード



　　 　　 　　

ちなみに(注1)で求めたλの近似分散にデルタ法を適用すると、対数ハザード差の分散を次のように近似できます。

　　
：デルタ法による近似分散

※ln(λ₁)とln(λ₂)は独立のため共分散 C(ln(λ₁)、ln(λ₂)) = 0

このようにデルタ法による近似分散と最尤法による近似分散は一致します。 ただし脱落例がある時の近似分散を用いると、デルタ法による近似分散と最尤法による近似分散は少し異なります。 しかし、いずれにせよどれも近似分散なので多変量の場合と整合性を取るために最尤法で求めた近似分散を用いることにします。 (→3.4 2標本の計数値　(2)名義尺度(分類データ) (注5))

表11.1.1のデータに指数分布モデルを当てはめ、実際に計算してみましょう。

(注1)より A群：d_A = 6　λ_A = 0.0194 　　B群：d_B = 8　λ_B = 0.0533
b₀=ln(λ_B) = -3.9448 　　b₁ = ln(λ_B/λ_A) = 1.0136
　　
b₁の検定：
b₁の95％信頼区間：b_1L^U = 1.0136 ± 1.96×√0.2917 = 1.0136 ± 1.0585 　下限：b_1L = -0.0449　上限：b_1U = 2.0721
ハザード比：HR = 2.7556
ハザード比の95％信頼区間　下限：HR_L = exp(-0.0449)=0.9561 　上限：HR_U = exp(2.0721) = 7.9417

このように2群比較に指数分布モデルを当てはめると、回帰係数の検定と推定を簡単に計算することができます。 この回帰係数の検定と推定はハザード比の検定と推定に相当するので、パラメトリック・ハザード比検定(PHR-test：parametric hazard ratio test)と名付けることにしました。 2群の累積生存率曲線を比較するには第2節で説明したノンパラメトリック手法よりも、このパラメトリック・ハザード比検定の方が合理的かつ効率的です。

(注3)　表11.4.2にワイブル分布モデルを当てはめ、対数尤度関数にニュートン・ラプソン法を適用して最尤解を求めると次のようになります。

ワイブル分布モデル：λ(t) = λa(λt)^a-1 = λ^aat^a-1 　　S(t) = exp{-(λt)^a}
ln(a) = h_a　　a = exp(h_a) 　　ln(λ) = β₀+β₁x_i1+…+β_jx_ij+…+β_px_ip = β'ｘ_i 　　λ = exp(β'ｘ_i)
　　
対数尤度：







ｈ_k+1 = ｈ_k - Ｈ_k^-1ｇ_k

ln(a)の初期値は0〜0.1程度にし、偏回帰係数の初期値は、指数分布モデルと同様に対数変換した生存時間を目的変数にした重回帰分析の偏回帰係数を求め、その符号を反対にしたものを用いると良いでしょう。 計算過程は省略しますが、表11.3.1のデータにワイブル分布モデルを当てはめると次のような結果になります。

ワイブル分布モデル：ln(a) = 0.131　　a = 1.398
ln(λ) = -3.629 - 0.616x₁ + 0.666x₂ 　　S(t) = exp{-(λt)^1.398}
○パラメータの検定と推定
・a:χ² = 1.482(p = 0.2235) ＜ χ²(1,0.05) = 3.841
　aの95％信頼区間　下限：a_L = exp(-0.0799) = 0.923 　上限：a_U = exp(0.342) = 1.407
・b₀:χ² = 263.152(p = 5.773×10^-15) ＞ χ²(1,0.05) = 3.841
　b₀の95％信頼区間　下限：b_0L = -4.068 　上限：b_0U = -3.191
・b₁:χ² = 6.502(p = 0.0108) ＞ χ²(1,0.05) = 3.841
　ハザード比：HR₁ = exp(-0.616) = 0.540
　ハザード比の95％信頼区間　下限：HR_1L = exp(-1.089) = 0.336 　上限：HR_1U = exp(-0.142) = 0.867
・b₂:χ² = 18.385(p = 0.00002) ＞ χ²(1,0.05) = 3.841
　ハザード比：HR₂ = exp(0.666) = 1.946
　ハザード比の95％信頼区間　下限：HR_2L = exp(0.362) = 1.436 　上限：HR_2U = exp(0.970) = 2.639
○尤度比検定
　このモデルの尤度：L(ｈ) = -231.082
　定数項だけのモデルの尤度：L(a,β₀) = -241.776
　回帰：χ_β² = 21.387(p = 0.00002) ＞ χ²(2,0.05) = 5.991
AIC(赤池の情報量基準) = 468.164

上記のようにaがほぼ1であり、偏回帰係数の値も指数分布モデルとよく似ています。 このことは、図11.6.8の指数分布モデルとワイブル分布モデルの生存関数S(t)がよく似ていることからもわかります。 そしてAICは指数分布の方がわずかに大きな値です。 これらのことから、このデータには指数分布モデルの方がより適していると考えられます。

なお共変数ｘが同じケースについてaと切片b₀だけのワイブル分布モデルを最尤法で求め、その対数尤度を合計することによって飽和モデルの対数尤度を求めることができます。 しかしこれは非常に面倒なので通常は計算しません。 そのため異質性の検定と疑似寄与率も通常は計算しません。 したがって普通は指数分布モデルを用いた方が何かと便利です。

(注4)　生命表解析の必要例数は、パラメトリック・ハザード比検定の原理を利用すると比較的簡単に求めることができます。

　　n = n₁ + n₂ 　　n₁：群1の例数　　n₂：群1の例数
λ₁、λ₂：群1と群2の母ハザード(瞬間死亡率つまり単位時間当たりの死亡率)
∵パラメトリック・ハザード比検定の計算式より

より 　
∴
時点tまで観測する時、死亡例数は d ≒ n{1 - exp(-λ t)} になる。 そこで時点tまで観測する時の全体の例数は次のようになる。

全例が死亡するまで観測する時(t = ∞)：n₁ = d₁
※{1 - exp(-λ₁t)} ≒ {1 - exp(-λ₂t)} ≒ {1 - exp(-λ_mt)} と近似すると少し簡便な式になる。
　　ただし
※λが不明の時は平均生存時間μ_tや50％生存時間μ'_t、あるいは観察期間tと観察終了時の生存率S(t)から求める。
　　 　　S(t) = exp(-λt) →
ただし観察終了時の生存率は誤差が大きくなりがちなので、平均生存時間や50％生存時間を利用した方が正確。

実際の試験では観測時間tを指定するのが難しい時もあると思います。 その時は便宜的にt = ∞として計算し、求められた必要例数を死亡例数として設定します。 そして実際の死亡例数がその例数に達するまで試験を続けます。 ただし試験計画を確実なものにするためには、やはり観測期間を指定する方が安全です。 (→1.8 科学的研究の種類とデザイン (注1))

有意水準を5％、検出力を80％、群1の母ハザード推定値を0.2(年単位)、群2の母ハザード推定値を0.1(年単位)、観測期間を5年、群1の例数と群2の例数の比率を1として、実際に計算してみましょう。

s = 1　　

n₂ = 68　　2群合計：n = 136　(2群合計死亡例数 = 68)

ちなみに全例が死亡するまで観測する場合は、次のように必要例数が約半分になります。 このように生命表解析では観測期間が長くなるほど必要例数が少なくなります。

n₂ = 33　　2群合計 = 2群合計死亡例数：n = 66

また探索型研究において、得られたハザードについてある程度の精度を確保したい時は信頼区間の原理に基づいて必要例数を計算します。

n：例数　　λ：母ハザード　　e：信頼度(1-α)での絶対精度(信頼区間幅の半分の値)

信頼係数を95％、母ハザード推定値を0.1(年単位)、絶対精度を0.05として実際に計算してみましょう。

・観測期間が5年の場合
　
・全例が死亡するまで観測する場合
　

ノンパラメトリック生命表解析の場合は必要例数の計算が非常に複雑になるので、上記のようにパラメトリック・ハザード比検定やハザードの信頼区間の原理に基づいて求めた値を流用すると便利です。

(注5)　死亡率関数f(t)は確率分布であり、t = 0〜∞まで積分すると1になります。 そこでf(t)に時間tを掛けて積分するとtの期待値E(t)つまり平均生存時間になります。

∴

ここでt = 0〜∞まで積分するとt・S(t) → 0になり、期待値E(t)は生存関数S(t)を積分した値になります。 そしてS(t)をt = 0〜tまで積分した値からt・S(t)を引いた値は、tの直前の時間をτとすると、S(t)をt = 0〜τまで積分した値になります。 tは死亡時間なのでτはtの直前までの生存時間になり、S(t)をt = 0〜τまで積分した値は死亡時間tの直前までの平均生存時間になります。 これが境界内平均生存時間です。

指数分布モデルの場合、境界内平均生存時間とその分散は次のようになります。 この分散を用いて境界内平均生存時間の信頼区間を求めることができます。

RMST(t=0-τ) =

また生存関数S(t)をカプラン・マイヤー法による累積生存率曲線で近似すると、その曲線下面積AUCが境界内平均生存時間の近似値になります。 第1節の(注1)と(注2)で説明したカプラン・マイヤー法による生存率の計算法とグリーンウッドの近似式を利用すると、境界内平均生存時間とその分散は次のようになります。

… 死亡例だけ足し合わせる

(注6)　指数分布モデルの生存関数とワイブル分布モデルの生存関数の交点は次のようにして求めることができます。

指数分布モデルの生存関数：S(t) = exp(-λ_et)
ワイブル分布モデルの生存関数：S(t) = exp{-(λ_wt)^a}
exp(-λ_et) = exp{-(λ_wt)^a} → λ_et = (λ_wt)^a → (a-1)ln(t) = ln(λ_e) - a・ln(λ_w)
∴

表11.6.1のデータについて求めると次のようになります。

A群の指数分布モデルの生存関数：S(t) = exp(-0.0194・t)
C群のワイブル分布モデルの生存関数：S(t) = exp{-(0.0280602・t)^2.46149}

最終更新日：2021年12月22日 　第5節へ 第12章へ

　webmaster@snap-tck.com 　Copyleft (C) 2000 SNAP(Sugimoto Norio Art Production)