玄関雑学の部屋雑学コーナー統計学入門

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

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

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

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

λ(t)=λ(定数)  S(t)=exp(-λt)  F(t)=1 - exp(-λt)   f(t)=λ exp(-λt) … 指数分布 (0<λ<1)
標的モデル
図11.3.1 各種の関数

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

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

○A群
λ A = 死亡数 観察期間合計 = 6 310 0.0194
S(t)=exp(-0.0194t)
平均生存時間=
50%生存時間(MST)=
○B群
λ B = 死亡数 観察期間合計 = 8 150 0.0533
S(t)=exp(-0.0533t)
平均生存時間=
50%生存時間(MST)=
○パラメトリック・ハザード比検定
ハザード比(B群/A群):
検定:χβ12=3.523(p=0.0605)<χ2(1,0.05)=3.841 … 有意水準5%で有意ではない
95%信頼区間 下限:HRL=0.9561 上限:HRU=7.9417
図11.6.1 指数分布モデルによる理論的生存曲線 図11.6.2 時間間隔を1にした時の理論的生存曲線

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

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

○A群
λ A = 死亡数 観察期間合計 = 6 138 0.0435
S(t)=exp(-0.0435t)
○B群
λ B = 死亡数 観察期間合計 = 8 77 0.1039
S(t)=exp(-0.1039t)
○パラメトリック・ハザード比検定
ハザード比 (B 群/ A ) HR 0.1039 = 0.0435 2.3896
検定:χβ12=2.602(p=0.1067)<χ2(1,0.05)=3.841 … 有意水準5%で有意ではない
95%信頼区間 下限:HRL=0.8291 上限:HRU=6.8870

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

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

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

IR P-T = 特定の期間中に新規発生した事象数 総危険人時間単位 × 10 k
総危険人時間単位:観測期間合計  k:2〜6

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

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

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

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

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

λ(t)=λa(λt)a-1aa ta-1  S(t)=exp{-(λt)a}  F(t)=1 - exp{-(λt)a}
f(t)=λa(λt)a-1exp{-(λt)a}=λaa ta-1exp{-(λt)a} …ワイブル分布 (0<λ<1、0<a)
直列モデル
図11.6.3 ワイブル分布モデル

(b) ガンマ分布モデル

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

  
   … ガンマ分布 (0<λ<1、0<b)
Γ(b):ガンマ関数  Γ(b;λt):不完全ガンマ関数
並列モデル
図11.6.4 ガンマ分布モデル

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

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

    
… 対数正規分布
… 標準正規分布の確率分布関数
… 標準正規分布の確率密度関数
図11.6.5 対数正規分布モデル

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

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

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

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

η=ln{λ(x1,…,xp)}=β0 + β1x1 + … + βjxj + … + βpxp + ε
λ(x1,…,xp)=exp(β0 + β1x1 + … + βjxj + … + βpxp + ε)
S(t)=e-λt=e-exp(β0 + β1x1 + … + βjxj + … + βpxp + ε)t
η:対数ハザード  β0:定数  βj:偏回帰係数(j=1,…,p)   ε:回帰誤差

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

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

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

xj=0の時の対数ハザード:ηj0=ln(λj0)
xj=1の時の対数ハザード:ηj1=ln(λj1)
対数ハザードの変化量:
  HRj:xjのハザード比

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

y=ln{λ(x1,x2)}=b0 + b1x1 + b2x2=-3.639-0.667x1 + 0.708x2
λ=exp(-3.639 - 0.667x1 + 0.708x2)   S(t|x1,x2)=e-λt=e-exp(-3.639 - 0.667x1 + 0.708x2)t
○x1:HR1=exp(b1)=exp(-0.667)=0.513
95%信頼区間 下限:HR1L=exp(-1.199)=0.302  上限:HR1U=exp(-0.135)=0.874
χβ12=6.030(p=0.0141)>χ2(1,0.05)=3.841 … 有意水準5%で有意
○x2:HR2=exp(b2)=exp(0.708)=2.030
95%信頼区間 下限:HR2L=exp(0.369)=1.447 上限:HR2U=exp(1.047)=2.849
χβ22=16.792(p=0.00004)>χ2(1,0.05)=3.841 … 有意水準5%で有意
x1:治療(0:無 1:有)   x2:重症度(0:症状無 1:軽症 2:重症)

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

○治療無(x1=0)で軽症(x2=1)の時の理論的生存関数:図11.6.5の赤色の曲線
y=-3.639 - 0.667×0 + 0.708×1=-2.931  λ=exp(-2.931)=0.0533   S(t|x1=0,x2=1)=e-0.0533t
○治療有(x1=1)で軽症(x2=1)の時の理論的生存関数:図11.6.5の青色の曲線
y=-3.639 - 0.667×1 + 0.708×1=-3.593  λ=exp(-3.593)=0.0274   S(t|x1=1,x2=1)=e-0.0274t
図11.6.6 指数分布モデルによる理論的生存関数

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

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

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


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

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

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

脱落例がある時の近似分散は、1例ごとに観測終了時での理論的累積死亡確率{1-exp(-λti)}を求め、それを合計したものを死亡例数の代わりに使うという考え方に基づいたものです。 これは死亡例が全て無限時間後に発生し、脱落例は全てその前に発生した時に死亡例合計を用いた最初の近似分散と一致します。 そして死亡例が無限時間より前に発生すると分散の分母を小さくし、脱落例が死亡例よりも後に発生すると分散の分母を大きくするので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個の一般的データ
共変数観測期間転帰
x11x1jx1pt1d1
:::::
xi1xijxiptidi
:::::
xn1xnjxnptndn
※di:死亡例は1、脱落例は0のダミー変数
指数分布モデル:λ(t)=λ  S(t)=exp(-λt)

η0n + β11 + … + βjj + … + βpp + ε=β + ε= + ε   =β
        
死亡例の確率:f(ti|i)=λ(ti|i)S(ti|i)   脱落例の確率:S(ti|i)
全体の尤度:
対数尤度:
※λi=exp(β'i)

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



k=k

k+1=k - k-1k

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

V()=-k-1f-1   E(-k)=f:情報行列   [f]jj-1f-1の第j対角要素
検定:≧χ2(1,α)の時、有意水準100α%で有意
推定:100(1-α)%信頼区間
→ 下限:  上限:

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

対数尤度差:
尤度比検定:-2(L(0) - L(β))=χβ2≧χ2(p,α)の時、有意水準100α%で有意

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

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

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

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

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

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

        
5=4 - 4-144
指数分布モデル:y=ln(λ)=-3.639 - 0.667x1 + 0.708x2   S(t)=exp(-λt)
○偏回帰係数の検定と推定

 ハザード比:HR1=exp(-0.667)=0.513
 ハザード比の95%信頼区間 下限:HR1L=exp(-1.199)=0.302  上限:HR1U=exp(-0.135)=0.874

 ハザード比:HR2=exp(0.708)=2.030
 ハザード比の95%信頼区間 下限:HR2L=exp(0.369)=1.447  上限:HR2U=exp(1.047)=2.849
○尤度比検定
 このモデルの尤度:L(β012)=-231.78
 定数項だけのモデルの尤度:L(β0)=56×{ln(0.036246) - 1}=56×(-4.317426)=-241.7759
 飽和のモデルの尤度:Lf=-231.183
 回帰:χβ2=-2×(-241.7759 + 231.78)=19.991(p=0.00005)>χ2(2,0.05)=5.991
 ズレ:χLOF2=-2×(-231.78 + 231.183)=1.194(p=0.754444)<χ2(3,0.05)=7.815
擬似寄与率:
AIC(赤池の情報量基準)=469.561

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

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

指数分布モデル:ηi=ln{λ(t|i)}=β0 + β1xi1 + εi=yi + εi (xi1:群1の時は0、群2の時は1となるダミー変数)
対数尤度:


n1:群1の例数  n2:群2の例数  d:全死亡数   :群1の観測期間合計   :群2の観測期間合計


  
d1:群1の死亡数  d2:群2の死亡数   λ1:群1のハザード  λ2:群2のハザード



        

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

  
デルタ法による近似分散

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

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

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

(注1)より A群:dA=6 λA=0.0194   B群:dB=8 λB=0.0533
b0=ln(λB)=-3.9448   b1=ln(λAB)=1.0136
  
b1の検定:
b1の95%信頼区間=1.0136±1.96×√0.2917=1.0136±1.0585  下限:b1L=-0.0449 上限:b1U=2.0721
ハザード比:HR=2.7556
ハザード比の95%信頼区間 下限:HRL=exp(-0.0449)=0.9561  上限:HRU=exp(2.0721)=7.9417

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

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

ワイブル分布モデル:λ(t)=λa(λt)a-1aata-1  S(t)=exp{-(λt)a}
ln(a)=ha  a=exp(ha)   ln(λ)=β01xi1+…+βjxij+…+βpxip=β'i   λ=exp(β'i)
  
対数尤度:







k+1=k - k-1k

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

ワイブル分布モデル:ln(a)=0.131  a=1.398
ln(λ)=-3.629 - 0.616x1 + 0.666x2   S(t)=exp{-(λt)1.398}
○パラメータの検定と推定
・a:χ2=1.482(p=0.2235)<χ2(1,0.05)=3.841
 aの95%信頼区間 下限:aL=exp(-0.0799)=0.923  上限:aU=exp(0.342)=1.407
・b02=263.152(p=5.773×10-15)>χ2(1,0.05)=3.841
 b0の95%信頼区間 下限:b0L=-4.068  上限:b0U=-3.191
・b12=6.502(p=0.0108)>χ2(1,0.05)=3.841
 ハザード比:HR1=exp(-0.616)=0.540
 ハザード比の95%信頼区間 下限:HR1L=exp(-1.089)=0.336  上限:HR1U=exp(-0.142)=0.867
・b22=18.385(p=0.00002)>χ2(1,0.05)=3.841
 ハザード比:HR2=exp(0.666)=1.946
 ハザード比の95%信頼区間 下限:HR2L=exp(0.362)=1.436  上限:HR2U=exp(0.970)=2.639
○尤度比検定
 このモデルの尤度:L()=-231.082
 定数項だけのモデルの尤度:L(a,β0)=-241.776
 回帰:χβ2=21.387(p=0.00002)>χ2(2,0.05)=5.991
AIC(赤池の情報量基準)=468.164
図11.6.7 指数分布モデルとワイブル分布モデル

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

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

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


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

より  

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

全例が死亡するまで観測する時(t=∞):n1=d1
※{1-exp(-λ1t)}≒{1-exp(-λ2t)}≒{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  

n2=68  2群合計:n=136 (2群合計死亡例数=68)

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


n2=33  2群合計=2群合計死亡例数:n=66

ノンパラメトリック生命表解析の場合は必要例数の計算が非常に複雑になるので、このパラメトリック・ハザード比検定の計算式から求めた値を流用すると便利です。