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

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

第1節と第2節で説明した生存率の計算方法と比較方法はノンパラメトリック生命表解析と呼ばれ、第4節で説明した比例ハザードモデルによる重回帰型生命表解析はセミパラメトリックモデルまたはセミノンパラメトリックモデルと呼ばれることがあります。 比例ハザードモテルによる重回帰型生命表解析がそのように呼ばれる理由は、この手法は特定の被験者のハザード比を予想できる点はパラメトリック的ですが、特定の被験者の生存率や生存時間を予想できず、ハザード比に基いた死亡順序しか予想できない点がノンパラメトリック的だからです。

この手法が生存率や生存時間を予想できないのは、ハザード関数λ(t)をブラックボックスにしたまま強引に生命表解析を行うからです。 生命表解析の大きな目的である生存率や生存時間を予想するためには、やはりハザード関数λ(t)の具体的な内容をきちんと定義して生命表解析を行う必要があります。 その正攻法の手法をパラメトリック生命表解析といい、その手法で用いられるモデルをパラメトリックモデルといいます。

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

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

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

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

第4節で説明したように、表11.1.1の時間間隔を全て1にした表11.2.3と時間間隔を間延びさせた表11.4.3に上記の指数分布モデルを当てはめると、図11.6.1と図11.6.2と図11.4.4のようになります。 そしてノンパラメトリック手法を用いると図11.6.1も図11.6.2も図11.4.4も全く同じハザード比3.697になり、その推定結果と検定結果も全く同じになります。 また比例ハザードモデルによる重回帰型生命表解析を用いても図11.6.1と図11.6.2と図11.4.4の結果は全く同じハザード比3.377になり、その推定結果と検定結果は全く同じになります。 これは明らかに非合理です。

しかしパラメトリック手法では、ハザード比は図11.6.1が2.756、図11.6.2が2.390、図11.4.4が4.391になり、それに応じて推定結果と検定結果も変わります。 このことからハザード関数λ(t)をきちんと定義して解を求めるパラメトリック手法の方が正確かつ合理的であることがわかると思います。

図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人」などと表現されます。 また暴露群と非暴露群の罹患率の比をIRR(Incidence Rate Ratio：罹患率比 または 発生率比)といい、例えば「IRR=2 だから暴露群の罹患率は非暴露群の罹患率の2倍ある」などと表現されます。

上記の定義式から指数分布モデルにおけるハザードλは人時間に基づく死亡率に相当し、ハザード比HRは人時間に基づく死亡率比に相当することがわかると思います。 したがって人時間に基づく罹患率と罹患率比は、特定の期間中は罹患率が一定であり、その間の累積罹患率曲線を指数関数によって近似できるという暗黙の前提を含んだ指標であることがわかります。

(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) パラメトリック多変量生命表解析

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

比例ハザードモデルの場合、ハザード関数は共変数によっても時間によっても形を変えず、ハザード関数と基準ハザード関数の比を対数変換した対数ハザード比と共変数の間に線形関係があると仮定します。 そしてハザード関数をブラックボックスにしたまま生命表解析を行うために、部分尤度法というトリックを用いて精度の悪い近似解を求めます。

それに対して指数分布モデルを用いた重回帰型モデルでは、ハザード関数が共変数とも時間とも無関係に一定であり、対数ハザードと共変数の間に線形関係があると仮定します。 そしてハザード関数を具体的に定義しているので、通常の最尤法を用いて精度の良い解を求めることができます。 しかもハザード関数に基いて理論的な生存関数を求めることができるので、特定の被験者の生存率や生存時間を予測することができます。 被験者の予後を予測することは生命表解析の大きな目的のひとつですから、これは非常に大きな利点です。

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

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

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

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

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

ハザード関数は対数ハザードηと共変数の重回帰式ですから、重回帰分析と同様に特定の共変数とηの偏回帰式を求めることができます。 そしてその式の両辺を指数変換すれば特定の共変数とハザードλの偏回帰曲線を表す式になります。 例えば治療が無い時と有る時について、重症度とハザードの偏回帰曲線を描くと図11.6.9のようになります。 (→7.2 重回帰分析結果の解釈)

上図からわかるように、ハザード比HR₁が同じ0.513でも、重症度が症状なし(x₂=0)と重症(x₂=2)の時では治療有と治療無のハザード差はかなり違います。 またハザード比HR₂が同じ2.030でも、症状なし(x₂=0)と軽症(x₂=1)のハザード差と軽症(x₂=1)と重症(x₂=2)のハザード差もけっこう違います。 このように、このグラフはハザードの具体的な値とその変化がわかるのでハザード比の意味を感覚的に理解するのに役立つと思います。

多変量生命表解析に限らず生命表解析全般で、現在はパラメトリックモデルよりも比例ハザードモデルやノンパラメトリック生命表解析の方が多用されています。 しかし比例ハザードモデルやノンパラメトリック生命表解析は結果の精度が悪くて普遍性がなく、理論的生存関数を求められないので予後を予測することはできません。

それに対してパラメトリック生命表解析は結果の精度が良くて普遍性があり、理論的生存関数を求められるので予後を予測することができます。 そしてパラメトリック生命表解析は、疾患による死亡の生理的なメカニズムを推測する時の参考にすることもできます。 さらにパラメトリック生命表解析は試験の必要例数を比較的簡単に計算できるので、正確な試験計画を立てられます。 そのためパラメトリック生命表解析はもっと利用されてしかるべきです。 _(注4)

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

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

○2群のハザードの差の検定：パラメトリック・ハザード差検定と区間推定
δ_λ = λ₁ - λ₂ 　　V(δ_λ) = V(λ₁) + V(λ₂) 　　　(λ₁とλ₂は独立なので共分散 C(λ₁,λ₂) = 0)

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対角要素

推定：100(1-α)％信頼区間
→ 下限： 　上限：

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

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

なお重回帰分析では、全て説明変数の値が平均値の時の目的変数yの値はy切片と一致します。 そのためこの場合も、全ての共変数の値が平均値の時のハザードは切片だけのモデルのハザードλ₀と一致すると考えがちだと思います。 でも重回帰分析のような最小2乗解と違って、この場合は最尤解なので一致するとは限りません。

また飽和モデルの対数尤度とそれを利用した異質性の検定、そして擬似寄与率と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₁の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 = λ^aa t^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.398：χ² = 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₀ = -3.629：χ² = 263.152(p = 5.773×10^-15) ＞ χ²(1,0.05) = 3.841
　b₀の95％信頼区間　下限：b_0L = -4.068 　上限：b_0U = -3.191
・b₁ = -0.616：χ² = 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₂ = 0.666：χ² = 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年の場合
　
・全例が死亡するまで観測する場合
　

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

最終更新日：2021年12月22日 　第5節へ 第7節へ
　Copyleft (C) 2000 SNAP(Sugimoto Norio Art Production)