Top 5 This Week

関連記事

7. パラメトリック生存モデル:指数・ワイブル・対数正規分布

- 本サイト運営者のサービスの紹介 -

パラメトリックアプローチの概要と位置づけ

生存時間解析において、観測対象がイベントを経験するまでの時間 $T$ を確率変数として扱います。基本的な量として生存関数

$$S(t) = P(T > t)$$

が定義され、時点 $t$ までイベントが発生しない確率を表します。ハザード関数は

$$h(t) = \frac{f(t)}{S(t)}$$

と定義され、$f(t)$ は生存時間の確率密度関数です。ハザード関数は時点 $t$ まで生存した個体が次の瞬間にイベントを経験する瞬間的な強度を表します。

生存時間解析の手法はノンパラメトリック法、半パラメトリック法、パラメトリック法に大別されます。カプラン・マイヤー推定量に代表されるノンパラメトリック法は分布形状に仮定を置かず、観測データのみから生存曲線を推定します。Cox比例ハザードモデルに代表される半パラメトリック法はベースラインハザードを非パラメトリックに扱いつつ共変量効果をパラメトリックに推定します。パラメトリック法は生存時間が特定の確率分布族に従うと仮定したうえでパラメータを推定します。

分布仮定を置くことにより推定の統計効率が向上し、観測期間外への外挿が可能となります。たとえば、多くの観測が打ち切られた試験終了時点を超えた将来の生存確率を推定できます。ただし、仮定した分布族がデータの真の生成過程と一致しない場合(分布誤特定)、推定量に系統的なバイアスが生じます。各分布はハザード関数の形状により特徴付けられ、指数分布は定数ハザード、ワイブル分布は単調ハザード、対数正規分布および対数ロジスティック分布は山型の非単調ハザードを表現できます。

仮定として、打ち切りメカニズムが生存時間と独立であること(非情報的打ち切り仮定)が必要です。また、生存時間が特定の確率分布族に従うことを仮定します。限界として、分布誤特定によるバイアスリスクが存在し、ノンパラメトリック法と比較してモデルの適合に関する追加的な診断が必要となります。

指数分布モデル

指数分布は最も単純なパラメトリック生存モデルであり、率パラメータ $\lambda_E > 0$ を持ちます。確率密度関数は

$$f(t) = \lambda_E \exp(-\lambda_E t), \quad t \geq 0$$

生存関数は

$$S(t) = \exp(-\lambda_E t)$$

ハザード関数は

$$h(t) = \lambda_E$$

であり、ハザードが時間に依存しない定数となります。これが指数分布の本質的な性質です。

指数分布の重要な数学的性質として記憶なし性があります。すでに時間 $s$ が経過した個体がさらに $t$ 以上生存する条件付き確率は

$$P(T > s + t \mid T > s) = P(T > t)$$

と表され、過去の経過時間が将来のイベント発生確率に影響を与えないことを示します。これは電子部品の偶発故障期(バスタブ曲線の中央部)のように、部品が摩耗せず劣化も進行しない状況の近似として用いられます。平均生存時間は $E[T] = 1/\lambda_E$ となります。

打ち切りデータを含む場合の最尤推定量は解析的に求まり、

$$\hat{\lambda}_E = \frac{d}{\displaystyle\sum_{i=1}^{n} t_i}$$

となります。ここで $d$ は観測期間内にイベントが発生した個体数、$\sum_{i=1}^{n} t_i$ は全観測時間の総和です。

仮定として、ハザードが時間に依存しない一定ハザード仮定が成立する必要があります。限界として、実際の故障・死亡過程ではハザードが経過時間とともに変化することが多く、適用範囲は限られます。また、後述するワイブル分布の形状パラメータ $\gamma = 1$ の特殊ケースに相当するため、単独モデルとしての表現力は最小限です。

ワイブル分布モデル

ワイブル分布は形状パラメータ $\gamma > 0$ と尺度パラメータ $\lambda > 0$ の2つのパラメータを持ち、信頼性工学で広く用いられます。生存関数は

$$S(t) = \exp\!\left(-\!\left(\frac{t}{\lambda}\right)^{\!\gamma}\right)$$

ハザード関数は

$$h(t) = \frac{\gamma}{\lambda}\left(\frac{t}{\lambda}\right)^{\!\gamma-1}$$

累積ハザード関数は

$$H(t) = \left(\frac{t}{\lambda}\right)^{\!\gamma}$$

です。形状パラメータ $\gamma$ によりハザードの形状が決まります。$\gamma < 1$ のとき単調減少ハザード(初期故障期)、$\gamma = 1$ のとき定数ハザード(偶発故障期)、$\gamma > 1$ のとき単調増加ハザード(摩耗故障期)となります。

$\gamma = 1$ のとき、ワイブル分布は指数分布と一致します。ただし、この対応における記号の関係に注意が必要です。ワイブル分布の尺度パラメータ $\lambda$ は、指数分布の節で用いた率パラメータ $\lambda_E$ とは異なる量であり、$\lambda = 1/\lambda_E$ の関係にあります。$\gamma = 1$ を代入するとワイブル分布のハザード関数は $h(t) = 1/\lambda = \lambda_E$ となり、指数分布の定数ハザードと一致します。

尺度パラメータ $\lambda$ は特性寿命とも呼ばれ、$S(\lambda) = \exp(-1) \approx 0.368$、すなわち全個体の約63.2%が故障する時点を表します。

ワイブルプロットは分布の適合を視覚的に診断する手法です。累積ハザード関数から $\log(-\log S(t)) = \gamma \log t – \gamma \log \lambda$ という線形関係が導かれるため、横軸に $\log t$、縦軸に $\log(-\log \hat{S}(t))$ を取ったプロットが直線状に並ぶかを確認します。傾きから形状パラメータ $\gamma$ をグラフィカルに推定することも可能です。

仮定として、ハザードが単調(単調増加または単調減少)であることが必要です。限界として、バスタブ型のような非単調ハザードは1つのワイブル分布では表現できません。また、小標本では形状パラメータの推定量に大きな標準誤差が生じます。

対数正規・対数ロジスティック分布モデル

指数分布とワイブル分布は単調ハザードしか表現できませんが、術後死亡リスクや感染症の急性期と回復期のように、時間とともに増加した後に減少する山型ハザードを持つ現象も存在します。対数正規分布と対数ロジスティック分布はこのような非単調ハザードを表現できます。

対数正規分布では、$\log T$ が正規分布 $N(\mu, \sigma^2)$ に従うと仮定します。生存関数は

$$S(t) = 1 – \Phi\!\left(\frac{\log t – \mu}{\sigma}\right)$$

ここで $\Phi$ は標準正規分布の累積分布関数です。

対数ロジスティック分布では、$\log T$ がロジスティック分布に従うと仮定します。生存関数は

$$S(t) = \frac{1}{1 + (t/\lambda)^{\gamma}}$$

ハザード関数は

$$h(t) = \frac{(\gamma/\lambda)(t/\lambda)^{\gamma-1}}{1 + (t/\lambda)^{\gamma}}$$

であり、$\gamma > 1$ のとき山型のハザード形状を持ちます。

両モデルはAFT(加速故障時間)モデルの枠組みで統一的に記述できます。$W$ を標準正規またはロジスティック分布に従う誤差項とすると、

$$\log T = \mu + \sigma W$$

と表されます。対数ロジスティック分布はさらに比例オッズ構造を持ちます。時点 $t$ における故障オッズは $(t/\lambda)^\gamma$ であり、共変量を導入することで比例オッズモデルとして解釈できます。

仮定として、対数変換した生存時間 $\log T$ が正規分布またはロジスティック分布に従うことが前提です。限界として、対数正規分布のハザードは $t \to \infty$ で $0$ に収束するため、長期生存者の多い集団では現実妥当性が低くなる場合があります。また、いずれの分布も最大1つの山しか持てず、複雑なハザード形状には対応できません。

各分布モデルのハザード関数形状の比較

(Fig1. 各分布モデルのハザード関数形状の比較(指数・ワイブル・対数正規・対数ロジスティック))

分布 生存関数 $S(t)$ ハザード関数 $h(t)$ ハザード形状 典型的適用場面
指数分布 $\exp(-\lambda_E t)$ $\lambda_E$(定数) 一定 偶発故障期・電子部品
ワイブル($\gamma < 1$) $\exp(-(t/\lambda)^\gamma)$ $({\gamma}/{\lambda})(t/\lambda)^{\gamma-1}$(減少) 単調減少 初期故障期・品質選別後部品
ワイブル($\gamma = 1$) $\exp(-t/\lambda)$ $1/\lambda = \lambda_E$(定数) 一定(指数分布と等価) 偶発故障期(ただし $\lambda=1/\lambda_E$)
ワイブル($\gamma > 1$) $\exp(-(t/\lambda)^\gamma)$ $({\gamma}/{\lambda})(t/\lambda)^{\gamma-1}$(増加) 単調増加 摩耗故障期・疲労破壊
対数正規 $1 – \Phi((\log t – \mu)/\sigma)$ 山型($t\to\infty$ で $0$) 非単調(山型) 術後回復・一部感染症
対数ロジスティック $1/(1+(t/\lambda)^\gamma)$ 山型($\gamma > 1$ のとき) 非単調(山型) 急性疾患・癌治療後経過

打ち切りデータの最尤推定

パラメトリック生存モデルのパラメータは、打ち切り情報を適切に反映した対数尤度関数の最大化により推定されます。観測値 $(t_i, \delta_i)$($\delta_i = 1$ は事象発生、$\delta_i = 0$ は打ち切り)に対する対数尤度関数は

$$\ell(\theta) = \sum_{i=1}^{n}\!\left[\delta_i \log f(t_i;\theta) + (1-\delta_i)\log S(t_i;\theta)\right]$$

と構成されます。完全観測($\delta_i = 1$)は確率密度 $f(t_i;\theta)$ で寄与し、「ちょうど $t_i$ に事象が起きた確率密度」を意味します。打ち切り観測($\delta_i = 0$)は生存関数 $S(t_i;\theta)$ で寄与し、「少なくとも $t_i$ まで生存した確率」を表します。

最尤推定量はスコア方程式

$$\frac{\partial \ell}{\partial \theta} = 0$$

の解として定義されます。指数分布では解析解が存在しますが、ワイブル分布や対数正規分布では一般に数値最適化が必要です。推定量の標準誤差は観測情報行列(フィッシャー情報行列の推定量)

$$\hat{I}(\theta) = -\left.\frac{\partial^2 \ell}{\partial \theta\, \partial \theta^{\top}}\right|_{\theta=\hat{\theta}}$$

から計算されます。Wald信頼区間は $\hat{\theta} \pm z_{\alpha/2} \cdot \mathrm{SE}(\hat{\theta})$ で構成されますが、パラメータが正の実数に制約される場合この対称区間が下限として負値を返すことがあり、対数スケールでの信頼区間構成が推奨されます。

仮定として、独立打ち切り仮定(打ち切りメカニズムが生存時間と独立であること)および正しい分布族の特定(正則条件の成立)が必要です。限界として、分布誤特定の場合、最尤推定量はバイアスを持ち一致性を失います。また、複数パラメータを持つモデルでは対数尤度曲面が複数の局所最大解を持つ場合があり、初期値の選択が推定結果に影響することがあります。

モデル比較と適合度診断

複数のパラメトリックモデルを比較する際には、AIC(赤池情報量規準)およびBIC(ベイズ情報量規準)が用いられます。

$$\mathrm{AIC} = -2\ell(\hat{\theta}) + 2k$$
$$\mathrm{BIC} = -2\ell(\hat{\theta}) + k\log n$$

ここで $k$ はパラメータ数、$n$ は標本サイズです。AICはパラメータ1つにつき2のペナルティを課し、BICは標本サイズに依存するより厳しいペナルティを課します。小標本ではAICが相対的に複雑なモデルを選択する傾向があり、大標本ではBICのペナルティが支配的となります。

視覚的な適合診断として確率プロット(ワイブルプロット・対数正規プロット)があります。ワイブルプロットでは $\log(-\log \hat{S}(t))$ を $\log t$ に対してプロットし、点列の直線性を確認します。Cox-Snell残差は

$$r_i = H(t_i;\hat{\theta}) = -\log S(t_i;\hat{\theta})$$

と定義されます。モデルが正しければ $r_i$ は近似的に標準指数分布(率パラメータ1)に従うため、カプラン・マイヤー推定量による $r_i$ の累積ハザードプロットが傾き1の直線から乖離していないかを確認します。入れ子モデル(指数分布はワイブル分布の特殊ケース)の比較には尤度比検定も有効です。

KM推定量と複数のパラメトリックモデルの生存曲線適合比較

(Fig2. KM推定量と複数のパラメトリックモデルの生存曲線適合比較(シミュレーションデータ))

ワイブル確率プロット:線形性による分布適合診断と形状パラメータのグラフィカル推定

(Fig3. ワイブル確率プロット:線形性による分布適合診断と形状パラメータのグラフィカル推定)

仮定として、比較するモデルが同一データセットに適用されていること、および漸近理論が適用できる十分な標本サイズが確保されていることが必要です。限界として、AICおよびBICは相対的な指標であり、比較対象のすべてのモデルが誤特定であっても最良とされるモデルが選択されてしまいます。確率プロットによる判断は主観的要素を含み、定量的な確信を単独では与えません。

信頼性工学への応用

産業機器の故障時間解析はパラメトリック生存モデルの主要な適用領域です。半導体デバイスや工業部品の故障時間データにワイブル分布を当てはめ、MTTF(平均故障時間)やB10寿命を推定して保証期間設計と予防保全スケジュールの最適化を行います。

MTTFは期待生存時間であり、各分布について次のように表されます。指数分布では

$$\mathrm{MTTF} = \frac{1}{\lambda_E}$$

ワイブル分布では

$$\mathrm{MTTF} = \lambda\,\Gamma\!\left(1 + \frac{1}{\gamma}\right)$$

となります。ここで $\Gamma(\cdot)$ はガンマ関数です。B10寿命(10%故障寿命)は

$$S(t_{B10}) = 0.9$$

を満たす時点 $t_{B10}$ として定義され、全個体の10%が故障するまでの時間を意味します。保証期間設計では $t_{B10}$ を基準として設定されることが多いです。信頼性関数は

$$R(t) = S(t)$$

と表され、時点 $t$ における個体の生存確率を工学的な文脈で表したものです。

バスタブ曲線は初期故障期(高ハザード・減少型)、偶発故障期(一定ハザード)、摩耗故障期(増加ハザード)の3つの段階からなります。それぞれワイブル分布の $\gamma < 1$、$\gamma = 1$($\lambda = 1/\lambda_E$)、$\gamma > 1$ に対応します。加速寿命試験では、高温・高電圧などのストレス条件下で取得したデータにパラメトリックモデルを当てはめ、推定したパラメータをアレニウス則などの加速モデルを介して通常使用条件へ外挿します。

仮定として、個々の機器の寿命がIID(独立同分布)に従うこと、および加速劣化条件から通常使用条件への外挿モデル(アレニウス則等)が成立することが必要です。限界として、実際の多部品システムでは部品間に故障依存性が生じることが多く、独立性仮定が成立しないケースがあります(直列・並列システム構造)。また、環境変動や経年劣化による非定常性は定常パラメトリックモデルでは捉えきれません。加速寿命試験で推定した分布パラメータを通常使用条件へ外挿する際、加速モデルの仮定が成立しない場合に信頼性予測が大きく外れるリスクがあります。

Popular Articles