統計学入門－第11章

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

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

このモデルの場合、死亡数と観察期間からλの最尤推定値を求めることができ、さらにλの逆数が平均生存時間になります。例えば、表11.1のデータに指数分布モデルを当てはめると次のようになります。 _(注1)

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

また図11.10のB群では、35ヶ月に最後の1例が死亡しているため、カプラン・マイヤー法で計算した累積生存率曲線は35ヶ月で累積生存率がストンと落ちて0になっています。しかし母集団は例数がもっと多い(理論的には無限例)ため、累積生存率が35ヶ月でいきなり0になるとは考えにくく、理論的生存関数のように35ヶ月以後も徐々に低くなっていくと考えられます。

このことから理論的生存率関数の方がより普遍性があり、将来のことをある程度は予測することができる、つまり外挿が可能であることがわかります。

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

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

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

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

(b) ガンマ分布モデル

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

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

これらのパラメトリックモデルのうち、最もよく用いられるのは指数分布モデルです。このモデルはパラメーターがλだけであり、その推定値を比較的簡単に計算することができます。一般的な相関分析や回帰分析では、2つの項目の関係を最も単純な直線関係で近似してしまうのが普通ですから、生命表解析でも最も単純な指数分布モデルで近似してしまってかまわないでしょう。

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

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

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

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

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

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

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

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

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

それに対してパラメトリック生命表解析は、結果がデータに左右されにくく、理論的な取り扱いが簡単で結果が解釈しやすく、しかもある程度は外挿が可能という特徴を持っています。したがって、パラメトリック生命表解析はもっと利用されてしかるべきだと思います。

(注1)　指数分布モデルの生存関数に最尤法を適用し、λの最尤推定値を求めると次のようになります。これらの値を用いて各種の検定や推定を行うことが可能ですが、残念ながら現在はまだあまり一般的ではありません。

死亡例の確率：f(t_i)=λ(t_i)S(t_i)=λe^-λt_i
脱落例の確率：S(t_i)=e^-λt_i
全体の尤度：

対数尤度：

∴

指数分布の期待値はλの逆数になり、中央値は期待値にln(2)を掛けた値になるため、λから平均生存時間と50％生存時間(MST)を求めることができます。

平均生存時間：

50％生存時間(MST)：

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

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

・指数分布モデル

死亡例の確率：f(t_i|ｘ_i)=λ(t_i|ｘ_i)S(t_i|ｘ_i)
脱落例の確率：S(t_i|ｘ_i)
全体の尤度：

対数尤度：

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

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

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

対数尤度差：

尤度比検定：

の時有意水準100α％で有意

偏回帰係数の初期値は、比例ハザードモデルと同様に、対数変換した生存時間を目的変数にした重回帰分析の偏回帰係数を求め、その符号を反対にしたものを用いることができます。定数項については、全ての共変数に平均値を代入した時の対数ハザードが、(注1)の計算式で求めたハザードλ₀の対数変換値と一致するように調整します。

表11.8のデータに指数分布モデルを当てはめ、実際に計算すると次のようになります。

・対数生存期間を目的変数にした重回帰分析の結果
　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
初期値：