ＪＭＰ V4 による生存時間分析

(1)

1

ＪＭＰ

V4 による生存時間分析

∼生存時間分析の実行∼

於：SAS東京

2000.11.18

（

株）

リコーＣＳ・

品質本部

廣野元久

(2)

１．はじめに

n

今回の生存時間分析のデータでは，

主に工業製品での信頼性試験データを扱う

n

例題は主に実験データ，要因の水準(値)

(3)

Motohisa HIRONO 3/186

n

生存時間分析の基礎

n

生存関数，ハザード関数，故障（死亡）の分類，

Weibull分布，打切りデータ，競合リスク

n

生存時間の分析（一変量の解析）

n

カプラン・マイヤー法，確率プロット

n

非線形回帰分析（要因解析）

n

比例ハザードモデル，加速モデル(Weibull回帰)

反応論モデル

本セミナーの内容

(4)

(5)

生存時間分析とは

n

ある基準となる時刻から，

目的の反応が得られるまでの

(生存)時間データ(Survival Time)を

対象とした解析手法の総称

n

工業分野では信頼性データ解析と呼ばれる

(6)

目的の反応とは

n

観測する個体に対して

1度だけ

非再起的に発生する事象

(Event)

n

例）

n

死亡

癌や循環器系の臨床研究での患者の再発・死亡

n

故障

システムや機器の信頼性研究での故障

(7)

生存時間データの特徴

n

時間データである（

非負、長期の観測）

n

観測打ち切りデータがある

n

故障や死亡の原因が複数ある

n

加速性や比例ハザード性がある

n

どの分布を仮定するかがあらかじめ不明

n

正規分布の理屈が使いにくい

n

データ収集に時間とお金がかかる

(8)

Time

Start of Study

End

Lost

Died

Time

Censor

Died

Start of Observation

Died

生存時間データ

(9)

Start of Study

End

Time

Censor

Died

Start of Observation

Died

生存時間データ

(10)

補足；

完全データ

(11)

補足；

定時打切試験のデータ

時点ｔ

ｃ

までは何個故障しても試験をする。

(12)

補足；

定数打切試験のデータ

Ｘ印が故障時点、０印が試験打切り時点

＃ｒとそれ以降とはｔの値は同じだが＃ｒは

完全データ、その他は打切りデータである

(13)

(14)

(15)

補足；

定時観測データ

（図中の破線のどこかでアイテムは故障したが

定時観測のため時刻が特定できない）

(16)

生存時間分析の目的

n

(目的の)の反応が，ある生存時間区間に

n

集中して発生するか

n

時間に依存しないでランダムに発生するか

n

生存している確率（信頼度）はいくらか

n

制御因子や環境因子の影響により生存時間

に違いはあるか

n

設計パラメータの決定

n

加速性，比例ハザード性による予測

(17)

Survival Time Analysis

生存時間の分布のご利益

注目する時点まで生存する

確率（Proportion）と

生存関数（Survival Function）の記述

製品の信頼性の記述，

生存時間データの分布の仮定

•例）

変圧器が

5000時間まで生存する確率

•例）

ベアリングの10％が

故障するまでの時間

(18)

生存時間データの分布の特徴

n

時間データである

非負である

n

中には非常に長生きな個体がある

分布の裾が右に長い

n

左右対称の単峯分布ではない

正規分布の理論が使い辛い

n

対象によって生存時間の分布が異なる

対数正規分布，指数分布

,Weibull分布など

(19)

生存関数

(Survival Function)

f(t)

t

故障

生存

密度関数 f(t)

生存関数S(t)とは時点tまで生存している

個体の累積比率(Proportion)を表す関数

分布関数F(t)とは

S(t)=1-F(t)

(20)

生存関数

(Survival Function)

S(t)

1-S(t)=F(t)

t

1

0 故障

故障

生存

1−分布関数を生存関数 S(t)

(21)

ハザード関数の導入

n

ハザード関数

(hazard function)とは

時点ｔ

まで生存したという条件付きで，

時点ｔ

の瞬間にイベントが発生する

確率

(rate)を表す関数

_{h(t)＝f(t)/S(t)}

n

例）

40歳死亡率

40歳で死亡するためには

40歳まで生きている必要があり，

その条件の人の中で次の

1年間で死亡する確率

(22)

ハザード関数と生存関数

S(t)

ﾊｻﾞｰﾄﾞは

時点ｔのS(t)を

１としたときの

S(t)の

傾きの絶対値

S(t)

( )

dS t

dt

(23)

分布を表現する関数の関係

( )

(

)

( )

( ) ( )

( )

(

)

( ) (

)

( )

(

( )

)

0 t 0 t 0

Pr

1 ,

Pr

|

lim

log

1 lim

t

S t

ob T t

F t

S t

f t F t

f u du

ob t T t

t T t

h t

t

d

S t

t

dS t

t S t

dt

S t

dt

∆ → ∆ →

=

≥

= −

=

≤ < +∆

≥

=

∆

−

+∆

=

=−

= −

∆ •

∫

生存関数；

分布関数；

密度関数；

ハザード関数；

(24)

分布を表現する関数の関係

( )

(

( )

)

( )

{

( )

}

0 log

log

exp

t

d

S t

h t

dt

H t

h u du

S t

H t

=−

=

=−

=

−

∫

ハザード関数；

累積ハザード関数；

(25)

n

生存関数、ハザード関数、累積ハザード関

数は数学的に等価

n

データ解析的には等価でない

n

ハザード関数は確率的な誤差の影響を受け

やすい

n

自動車の瞬間速度の測定が困難であるのと

同じ

n

臨床では生存関数、信頼性では累積ハザード

関数が好まれる

分布を表現する関数の関係

(26)

n

ハザード関数が大きな値をとるほど故障の

リスクは高く、生存関数は早く０に近づく

分布を表現する関数の関係

( )

( ) ( )

_{( )}

( )

0

0 ,

exp

a

t

a

t

a

h t

a h t

h t

a h t

a

h t

S

t

a h u du

h u du

S t

• = •

=









=

_

− •

_

=

_

−

_











_

_



=

_

_

−

_

_

=



_

_







∫

a倍

a乗

(27)

故障（

死亡）の分類

n

ハザード関数は瞬間死亡率であるからハ

ザード関数の傾向により

3つに分類できる

n

初期故障型；ハザード関数が単調減少

n

例）

乳幼時期の死亡率

n

偶発故障型；ハザード関数が一定

n

例）

青壮年期の死亡率

n

磨耗故障型；ハザード関数が単調増加

n

例）

老年期の死亡率

(28)

補足；

初期故障型

（DFR：decreasing failure rate)

h(t)

(29)

補足；

偶発故障型

（CFR：constant failure rate）

h(t)

(30)

補足；

摩耗故障型

（Ｉ

FR：increasing failure rate)

h(t)

(31)

補足；

バスタブ曲線

(Bath-Tab Curve）

h(t)

t

初期故障型

摩耗故障型

偶発故障型

合成曲線

部品･ユニットの故障率は実際には前述の各故障率のパターン

が合成された形で時間の関数として表されることが多い

(32)

Weibull分布

n

工業では主に

Weibull分布を仮定

( )

t

1

exp

t

f t

β β

β

α α

α

−

_

_



 

=

_{ }

• _

−

_{ }

_

 

_

_

 

_

_

( )

1 ( )

exp

t

S t

F t

β

α







 



= −

=

_

−

_{ }

_

 







( )

f t

_{( )}

( )

t

1 h t

S t

β

α α

−

 

=

_{=  }

 

密度関数

生存関数

ハザード関数

(33)

Weibull分布のパラメータ

Weibull分布のパラメータはα，β，γの3つ

通常は位置パラメータγを0と仮定する

形状パラメータβと尺度パラメータαを推定する

形状パラメータは分布の形を決めるものであり

Weibull分布では

β ＜1 のとき初期故障型

β ＝1 のとき偶発故障型（指数分布）

β ＞1 のとき磨耗故障型

に対応する

(34)

補足；

Weibull分布

0

0.2

0.4

0.6

0.8

1

0

2

4 6 t

f(t)

β=.5,α=2 β=1,α=2 β=2,α=2 β=4,α=2

0

0.5

1

1.5

2

2.5

0

0.5

1

1.5 _{2 t}

h(t)

β=.5,α=2 β=1,α=2 β=2,α=2 β=4,α=2

α =2 で βを変えたワイブル分布の f(t),h(t)

(35)

Weibull分布のパラメータ

尺度パラメータαは

累積の故障割合が

63.2％に達するときの時点

注）形状パラメータの表記は

JMPではβ(Beta) 日本の信頼性データ解析ではｍ

尺度パラメータの表記は

JMPではα(Alpha) 日本の信頼性データ解析ではη

(36)

n

データを読み込む

n

Weibull分布の確認

n

オーバーレイプロットを描く

n

1）密度関数を描く

n

2）生存関数と分布関数を描く

n

3）3つのハザード関数を描く

ハンズオン

1 Weibull.JMPの解析

(37)

変数の意味など

n

変数

time が生存時間

n

α＝

10000時間，β＝2.5のWeibull分布

n

time，f(t),F(t),S(t)

n

βの値を変えてみよう

n

h1(t)…β＝2.5

n

h2(t)…β＝1.0

n

h3(t)…β＝0.5

n

データ件数

40件

ﾊｻﾞｰﾄﾞ関数はどう変わる

(38)

操作

1.1 ｵｰﾊﾞｰレイﾌﾟﾛｯﾄを描く

(39)

操作

1.2 変数の役割を指定

1.Timeを選択し

2.Xをｸﾘｯｸする

3.f(t)を選択し

4.Yをｸﾘｯｸする

5.OKを選択

(40)

操作1.3 密度関数の描画

1.Connect

(41)

操作1.4 生存関数と

分布関数の描画

同様な操作で，分布関数，生存関数を描画する

同様な操作で，ﾊｻﾞｰﾄﾞ関数を描画する

(42)

ハンズオン

1 Weibull.JMPのまとめ

グラフより，ワイブル分布の

ハザード関数h(t)は

形状パラメータβの値によって，

単調減少 β＜1

一定 β＝1

単調増加 β＞1

することが分かる

(43)

ＪＭＰ

Calculator

n

メニューの

ColsでColumn Info.を選ぶ

n

Current Properties がFormulaを確認

n

Edit Formula をクリック

n

計算式を表示

,編集

n

2つのリストボックスの役割を理解

n

左；現在登録されている

Col Name

n

右；利用する関数群

n

クリックすると関数を表示

(44)

練習問題

n

LnNormal.JMPを読み込み,ハザード関数

と生存関数を求め

,グラフにしてみよう

( )

1 exp

1 log

2

2 t

f t

t

µ

σ

π



₋











=

_

−

_

_

_











対数正規分布の密度関数

(45)

補足；

正規分布

(

)

2 2

1 ( )

exp

,

2

2 x

f x

µ

x

σ

σ π



₋





=

_

−

_

− ∞ < < +∞







(

)

2 2

1 ( )

exp

,

2

x

t

F x

µ

dt

x

σ

σ π

−∞



₋





=

_

−

_

− ∞ < < +∞







∫

(46)

補足；

指数分布

故障の発生が時間に依存しない

( )

exp

( )

0,

0 f t

=

λ

−

λ

t t

≥

λ

>

( )

1 exp

( )

0,

0 F t

= −

−

λ

t t

≥

λ

>

指数型累積ハザードプロット

f(t) t

(47)

補足；

対数正規分布

n

ＥＭによる寿命に適用される

(

)

2 2

ln

1

1 ( )

exp

0

2

2 t

f t

t

µ

σ

σ π



₋





 

=

_{ }

_

−

_

>

 

_

_

_

_

( )

(

)

2 2 0

ln

1

1 exp

0

2

t

x

F t

dxt

x

µ

σ

π



₋





=

_

−

_

>







∫

0 0.05 0.1 0.15 0.2 0 3 6 9 12 t f(t) _μ=2,σ=.5 μ=2,σ=1 μ=2,σ=2 0 0.05 0.1 0.15 0.2 0.25 0 3 6 9 12 t h(t) μ=2,σ=.5 μ=2,σ=1 μ=2,σ=2

対数正規分布の f(t)と h(t)

(48)

補足；

寿命に用いられる主な分布

n

ゴンぺルツ分布

ワイブル分布

n

正規分布

対数正規分布

n

ロジスティック分布

対数ロジスティック分布

・指数分布＿偶発故障

n

ガンマ分布，極値分布，一般化レーリー分布

など

(49)

49

(50)

補足；

確率プロット

確率紙

多くの分布に対応したものが開発されている

確率紙はデータを昇順に並べたとき,

データと仮定する分布関数との関係が

直線になるように工夫された用紙である

(51)

生存時間の分布を調べる方法

n

確率プロットの利用

n

V4からヒストグラムに各種の確率プロットが追加

された

n

正規分布

,対数正規分布

n

weibull分布,3パラメータweibull分布

n

極値分布

,指数分布

n

ガンマ分布

,ベータ分布

n

Weibullプロットを例にとる

(52)

Weibull プロットの原理

両辺に対数をとると

再度，両辺に対数をとると

( )

1 ( )

exp

t

S t

F t

β

α







 



= −

=

_

−

_{ }

_

 







( )

log

S t

t

β

α







 



= −

_

_{ }

_

 







log

( )

t

S t

β

α

 

−

_{=  }

 

( )

{

}

log

S t

β

log

t

α

 

−

=

_{ }

 

( )

{

}

(53)

ハンズオン

2 Weibull.JMPの解析2

n

データを読み込む

n

Weibullプロットの原理を理解する

n

Columnを1つ追加する

n

ＪＭＰ

Calculatorを使い,-lnS(t)を計算する

n

AnalyzeのFit Y ｂｙX で,Xにtime,Yに-lnS(t)を

指定する

n

縦軸と横軸を対数にする

n

プロットが直線になることを確認する

(54)

ハンズオン

2 Weibull.JMPのまとめ

n

Weibullプロットの原理は

(55)

練習問題

n

Wiring．JMPを読み込み、Timeのヒストグ

ラムを作ってみよう

n

AnalyzeからDistributionを選ぶ

n

Fit Distributionを使いWeibull分布を当て

はめてみよう

n

他の分布を当てはめたときはプロットはど

うなるか調べてみよう

(56)

補足；

プロットが直線でないとき

n

プロットが直線的でない場合がある

n

寿命分布が

Weibull分布に従わない

n

位置パラメータγが予想される

n

打切りデータを考慮していないなど

n

位置パラメータγが存在する場合の解析

n

競合リスクモデルの解析

n

混合分布の解析

n

複合

Weibull分布の解析

(57)

57

(58)

補足；

信頼性について

劣化状態

故障時点

；ｙ（

測定

）

初期特性

or

製造条件

ストレス

加速寿命試験

市場データ

(59)

補足；

製品などの寿命予測の流れ

スタート 信頼性試験の 計画

・設計パラメータの決定

・反応速度論モデルなどを仮定して環境因子と水準の決定

・寿命分布の仮定

・実験計画法を利用した割付と試験計画の立案

・故障の定義など

(60)

補足；

製品などの寿命予測の流れ

信頼性試験の

計画

試験の実施

・個々の試料の故障解析

・寿命の測定

故障モード決定 データ解析開始 変数変換

(61)

補足；

製品などの寿命予測の流れ

・寿命分布の確認 多変量解析 ・残差の検討 ・加速係数や活性化エネルギーの評価 ・実寿命予測の妥当性 データ解析開始 1変量解析 変数変換 例） 温度因子 摂氏 → 1/（kBT） 寿命予測モデル の作成 ＯＫ 2変量解析 ・曲線関係 ・加速性のﾁｪｯｸ ・時間依存性

(62)

補足；

工程データ解析との比較

解析ステップ

工程データ

生存時間ﾃﾞｰﾀ

1 １変量解析

基本統計量

ヒストグラム

箱髯図

Kaplan-Meier推定量

生存関数プロット

確率紙・累積ハザード紙

2 １因子(質)

分散分析

ログランク検定

_{一般化Wilcoxon検定}

3 1因子(量)

単回帰分析

_{直交多項式}

反応論による加速モデル

_{比例ハザードモデル}

4 多因子(質量)

実験計画法

_{重回帰分析}

反応論による複合加速モデル

_{複合比例ハザードモデル}

(63)

Survival メニュー

n

Survival Distribution

n

ノンパラメトリックの

Kaplan-Meier法

n

仮定した分布のパラメータを求める確率プロット

n

右側打ち切り，グループの比較，競合リスク

n

Parametric Regression

n

生存時間に分布を仮定した非線形回帰

n

加速モデル

n

Proportional Hazards

n

生存時間に分布を仮定しない

Cox回帰

n

比例ハザードモデル（

セミパラメトリック）

n

Recurrence 今回は対象外

(64)

４．打ち切りデータの処理と生存

時間分布（

Weibull分布）の確認

(65)

打ち切り

(censor)データ

n

観測される個体について，正確な生存時

間が測定できるとは限らないので打ち切り

(censor)が生じる場合がある

1.時間の原点（

測定の開始時点）

が不明確な場合

例）製品やシステムなどの信頼性研究では

納品時点は分かっても，実際のユーザーの

使用開始時点は正確には分からない

例）デバイスなどでは信頼性試験終了時点に，

故障していない個体がある

2.反応の発生時点が分からない場合

左側打ち切り

右側打ち切り

(66)

Censor変数の意味

n

ＪＭＰでは，

Censor変数にルールがある

n

0；

打ち切りのない完全なデータを意味する

n

1，2，…（0以外）；

打ち切りが生じたデータ

(67)

ハンズオン

3 Wiring.JMPの解析

n

データを読み込む

n

一変量の解析

,分布を確認する

n

Survival Disuribution を使う

n

Survival プロットを描く

n

Weibullプロットを描く

n

Weibullプロットに参照線を追加する

n

Weibull分布のパラメータを推定する

n

他の分布（対数正規分布，指数分布）を試す

(68)

変数の意味など

n

電子デバイス

(Al配線)の加速寿命試験における

寿命分布を推定する

.

n

変数の意味

n

Time(加速寿命試験による電子デバイスが故

障に至るまでの時間

)

n

Censor(打切り)

打切りがあるので、Distributionは

使ってはいけない

(69)

操作

3.1 打切り情報を

無視した解析

1.Survival Distribution をクリックして

2.Timeをクリックして

3.Yをクリック

4.OKをクリック

(70)

操作

3.2 Weibullプロットの実行

Weibull Plot

Weibull Fit

(71)

操作

3.3 パラメータの推定

Weibull プロットが

描画される

打ち切りが考慮されて

いないので生存時間が短めに推定される

(72)

操作

3.4 打切り情報を

考慮した解析

1.Survival Distribution をクリックして

2.Timeをクリックして

3.Yをクリック

6.OKをクリック

4.Censorをクリックして

5.Censorをクリック

(73)

操作

3.5 信頼区間の追加

1.Plot Optionsから

Show Confid Intervalを選択

2. Kaplan-Meier 法による

平均値,標準偏差の推定値

95％信頼区間が描画される

(74)

操作

3.6 Weibullプロットの実行

2.Weibull プロットを描画する

3.Save Estimatesを選択

1.Weibull Plot,Weibull Fit を選択

(75)

操作

3.7 推定値の保存と単回帰

(76)

ハンズオン

3 Wiring.JMPのまとめ

n

グラフの下に

Extreme-value Parameter Estimates と

Weibull Parameter Estimates が表示される

n

表には，

Weibull分布のパラメータの推定値の他に，

推定値の信頼水準

95％の上下限値，および反応数が

表示される

n

DeltaはWeibullプロットの傾きでBeta=1/Delta

n

LambdaはWeibullプロットの63.2％点で

Alpha=e

Lambda

n

Weibull分布の他にExponential PlotとLogNormal

Plotができる

(77)

Kaplan･Meier 推定量

n

故障

(イベント)があった時点を

とする

n

t

_i

時点の故障数をd

_i

とする

n

t

_i

時点の直前までのリスク集合の大きさをn

_i

とする

n

リスク集合の大きさとは

,その直前まで生き残っていた個体数

n

故障が発生するリスクにさらされた個体数

n

リスク集合の大きさは

,故障と打切りにより時点ごとに異なる

(時間の経過とともにリスク集合の大きさは小さくなる)

{ }

1 n i i

t

₌

( )

1

2

1

2 ˆ

₁

i

n

i

t

n

i

d

S t

n

_<

n











 



= −

_

_{ }

× −

_

× × −

_

_

=

_

−

_



 



L





∏





(78)

Kaplan･Meier推定量

n

カプラン･マイヤー推定量は

,各時点の生存

率の積になっている

n

対数を取ると

n

各項は近似的に独立なので

( )

1

2

1

2 ˆ

₁

i

n

i

t

n

i

d

S t

n

_<

n











 



= −

_

_{ }

× −

_

× × −

_

_

=

_

−

_



 



L





∏





( )

1

ˆ

log

log 1

n i i _i

d

S t

n

=





=

_

−

_





∑

1

ˆ

log

log 1

n i i i

d

V

S

V

n

=











_{ =}

_

_

₋

_

_





∑

_

_

_

_

(79)

Kaplan-Meier推定量

n

テーラー展開の近似により

n

二項分散により

n

したがって

,

2

1 ˆ

ˆ

log

ˆ

V

S

V S

S









_≈

_

_

_•

 





_

_

 

(

)

log 1

i i i i i i

d

V

n

n n

d









−

=









₋









(

)

1

ˆ

log

n i i _i _i _i

d

V

S

n n

d

=



_{ =}





∑

₋

(

)

2 1

ˆ

n i i i i i

d

V S

S

n n

d

=

  =

 

∑

₋

グリーンウッドの公式

(80)

Kaplan-Meier推定量

n

打切りがない場合は

n

率の分散の公式に一致する

1

1 ˆ

n

i

n

d

S

n

=

−

=

∑

( )

1 ˆ

ˆ

₁

ˆ

_/

V S

  =

_{ }

S

−

S

n

(81)

補足；

累積ハザード法

n

時点

t

_i

におけるハザード率は､ハザード関

数による表現により､

( )

(

,

)

,

i

t t

t

h t

t

n

∆

+

× =

における故障数

時間以上の寿命を持つ個数

Δｔ

_{は微小な値で､限りなく0に近づけると}

( )

1 i

k

i

k

t

H t

=

∑

k

における故障数,r

t 時間以上の寿命を持つ個数n

タイがないときは

1

1 ( )

t

i

k

H t

=

∑

逆順位

(82)

補足；

データ例

n

39､79，90，115 66，96 （ｈ)

完全データ

打切データ

から分布関数を推定する

(83)

補足；

計算例

6 観測値打切り順位逆順位故障数 λ(ｔ

)

H(t)

F(t)

39

0

1

6 1 0.1667 0.1667 0.15352

66

1

2

5

0 0 0.1667

73

0

3

4

1 0.25 0.4167 0.34076

90

0

4

3 1 0.3333

0.75 0.52763

96

1

5

2

0

0.75

115

0

6

1

1 1.75 0.82623

=count(データ範囲)

データ､打切り有無､順位を入力

＝A$1-C3+1

＝E3/D3

If文を使う

累積する

If文を使い完全データだけF(t)

を推定する

(84)

補足；

ジョンソン法

n

打切りを考慮した平均順位を計算して

F(t)

を推定する

n

打切りがなければ､故障順位と順位は一致

1

1 1 (

1)

j

n

MON

n

i

−

+ −

=

+

+ − −

打切りがあっても順位を推定できる方法

(85)

補足；

ジョンソン法

6 順位(i)

打切り故障順位(j) k=i-1 MONj

1

0

1

0

1

2

0

2

1

2

3

0

3

2

3

4

0

4

3

4

5

0

5

4

5

6

0

6

5

6 順位(i)

打切り故障順位(j) k=i-1 MONj' MONj

1

0

1

0

1

2

1

3

0

2

2.2

4

0

3

3.4

5

1

3

4

3.4

6

0

4

5

5.2

(86)

補足；

推定方法の比較

ジョンソン法と累積ハザード法とK-M法

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

(87)

87

(88)

1因子（質的変数）の解析

n

処理の違うグループ間の生存時間の比較

n

Kaplan-Meier 法の実行

n

グループ毎の生存関数の表示

n

ノンパラメトリック検定

n

ログランク検定

n

一般化

Wilcoxon検定

n

パラメトリック検定

n

生存時間分布の指定

n

パラメトリックモデルの実行

(89)

ノンパラメトリック検定

n

ログランク検定

n

全ての時点で同じウエイトがかかる

n

一般化

Wilcoxon検定

n

ウエイトがその時点のリスク集合の大きさ

n

ｔ

が小さい初期のウエイトが大きい

(90)

補足；

ログランク検定

封止樹脂 A(群 1)：150，170，220，280

封止樹脂 B(群 2)：221，290，340，390

(値は人工データ；

データ数は説明を簡便にするためのものであり,

このような小さなサンプルサイズを薦めている

のではない)

(91)

補足；

ログランク検定

第

j

時点

(

j

=

1,2,

L

,

k

)

における 2×2 の分割表を

以下のように表記する.

1 1 1 1 1 1 1 1 11 11 11 11 2 2 2 2 2 2 2 2 21 21 21 21 1 1 1 1

1 ₁

1

2

j j j j j j j k k k k j j j j k k k k j j j j k k k k

t

_t

t

d

n

d

n

_d

_n

_d

_n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

d

n

−

₋

−

₋

−

L

時点故障

生存

計

時点故障

生存

計

時点故障

生存

計

第群

計

_計

計

両群の生存時間が同じという条件（期待度数）では,

故障数は平均的に

1 j

j

1 j

/

j

,

2 j

j

2 j

/

j

m

= ×

d

n

m

= ×

d

n

(92)

補足；

ログランク検定

分散は

(

)

(

)

(

)

(

)

1 2 1 1 1 ₂ 1 2 2 2 2 ₂

1

j j j ---j j j

---n -d

,

n -1

n -d

n -1

j j j j j j j j j j j j j j j j j j j j j j j j j j

n n d

n

d

n

v

d

n

n n d

n

d

n

v

d

n

−



 



=

_

−

_{ }

_

=

−



 



−



 



=

_

−

_{ }

_

=

−



 



1

,

/

j j j

n

=

d

p

=

n

の二項分布を考え

下線部分の有限母集団修正を加えたもの

と考えれば良い.

(93)

補足；

ログランク検定

故障数とその平均との差

1j

d

1j

m

1j

,

2 j

d

2 j

m

2 j

δ

=

−

δ

=

−

が大きくなければ,生存関数に有意な差は認められない

(

)

(

)

1j

d

1j

m

1j

,

2j

d

2j

m

2j

δ

=

−

δ

=

−

∑

を調べれば良いだろう.

この差の二乗を分散で割った値が

近似的にχ2 分布に従うことを使って検定する.

(94)

補足；

ログランク検定

同時点での故障がなければ

1 1 1 1 1 1 1 1 2 1 ₂

1

j j j j j j j j j j j j j j j

n

d

n

v

n

δ



₋











_

=



_

−



_

=



_

−



_

=













₋



×

=

∑

群の故障 1

群2の故障

(95)

補足；

一般化ウィルコクスン検定

ログランク検定を少しかき方を変える

(

)

(

)

(

)

1 1 1 1 1 1 2 1

,

1

j j j j j j j j j j j j j j j j j j

n

w d

m

w

d

n

n n d

n

d

v

w

n

δ

=

−

=



_

−



_





−

=

−

∑

(

)

(

)

(

)

1 1 1 1 1 1 2 1

,

1

j j j j j j j j j j j j j j j j j j j

n

w d

m

w

d

n

n n d

n

d

v

w

n

n n

δ

=

−

=



_

−



_





−

=

−

∑

一般化ウィルコクスン検定では

(96)

補足；

一般化ウィルコクスン検定

同時点での故障がなければ

(

)

(

)

(

)

1 1 1 2 1 1 1 1 1 2 1 1 2

1

2

j j j j j j j j j j j j j j j j j j j j j j j

n

w

d

n d

n

n n d

n

d

v

w

n n

n

δ

=



_

−

_



_

=

−

=

Σ



_

−

=

−







−

=

−

∑

群の故障

(97)

(98)

ハンズオン

4 Rats.JMPの解析

n

データを読み込む

n

一因子（質的変数）の解析

,処理の比較

n

Survival Time Modelingを使う

n

2群のSurvival プロットを描く

n

２群の生存時間の違いを検定する

n

推定量を保存する

n

２群をまとめて生存時間を推定する

(99)

変数の意味など

n

発癌剤の投与量を変えた

2つのグループ

のラットの生存時間の比較

n

解析に用いる変数

n

day （ラットの生存日数）

n

Group （発癌剤の投与量の違い)

n

Censor （打ち切りの有無)

(100)

操作

4.1 役割の指定

1.Survival Distribution をクリックして

2.daysをクリックして

3.Yをクリック

8.OKをクリック

4.Censorをクリックして

5.Censorをクリック

6.Groupをクリックして

(101)

操作

4.2 ノンパラ検定の実行

5％有意でない

2群のラットの

生存日の違いは

このデータからでは

分からない

(102)

操作

4.3 カテゴリの併合

カテゴリー(処理の違い)を

併合してみる

Plot Optionsから

Show Conbinedを

クリックする

併合されたときの統計量が出力される

(103)

ハンズオン

4 Rats.JMPのまとめ

n

Tests Between Groupsに表示されているよう

に，2群の生存関数が等しいという帰無仮説の

下での検定結果（Prob>ChiSq）から，いずれの

検定でも5％で有意でないことが分かる.

n

ウインドウの左下隅にあるﾁｪｯｸボタンをクリッ

クして，Show Combinedを選ぶと，2群を合併

した生存関数がグラフに追加される.

n

グラフの軸をダブルクリックするとグラフにグリット

を表示できるなどのオプション機能がある

(104)

生存関数の推定値の出力

n

左端の

daysには，時間の順に反応，打ち切りが

あった日が出力される

.

n

2番目のSurvivalのカラムは，生存関数の推定

値を示している

n

Survival=(At Risk - N Failed)/At Risk

n

3番目のFailureのカラムは，累積のイベントの

割合（分布関数）を示しており，

Survival+Failure=1

という関係がある

n

SurvStdErr(Survival Standard Error)の

カラムは，得られた生存関数の値を母集団に対

する推定値と考えた場合の標準誤差を表示

(105)

生存関数の推定値の出力

n

N Failedのカラムは1番目のカラムの時点で発生

した反応の数である

n

N Censoredのカラムは 1番目のカラムの時点で

打ち切りのあった数である

n

右端の

At Riskのカラムは，リスク集合の大きさ（ス

タート時の標本数から，その時点までの反応と打ち

切りを受けた個体を除いた数）を示している

n

Quantilesの表にある，MeanとStdDev

（

Standard Deviation）は生存時間の平均の推

定値と標準偏差の推定値である

(106)

練習問題

n

PWB．JMPを読み込み周辺温度の違いが

寿命にどのような影響を与えるか考察して

みよう．

n

各条件の寿命は

Weibull分布にしたがってい

ると考えられるか

n

周辺温度の違いが寿命に影響を与えると考え

られるか

n

周辺温度は温度にどのくらい影響を与えるか

(107)

補足；

PWBのデータ

温度条件

生存時間

65℃

180 270 370 430 500 570 620 710 790 960

75℃

90 170 210 220 270 310 320 360 430 540

85℃

30 50 60 70 90 100 110 120 130 150

(108)

補足；Ｉ

ＣのＥＭ寿命の

ロットによる層別

ロット1 生存時間 143 156 178 180 192 206 209 213 216 打切り 0 0 0 0 0 0 0 0 0 生存時間 216 220 224 227 230 234 246 265 304 打切り 1 0 1 0 0 0 0 0 0 ロット2 生存時間 163 172 185 198 204 205 222 225 230 233 打切り 0 0 0 0 1 0 0 0 0 0 生存時間 233 235 239 240 261 280 280 296 323 344 打切り 1 0 0 0 0 0 0 0 0 1

ロット別の EM 寿命

(109)

χ

2

_{値自由度ｐ値}

ログランク検定 2.9414 1 0.0863

一般化 Wilcoxon 検定 2.4823 1 0.1151

補足；Ｉ

ＣのＥＭ寿命の

ロットによる層別

(110)

μ

σ

L95%μ U95%μ L95%σ U95%σ 故障数

ロット1 5.36851 0.18645 5.27724 5.46453 0.13642 0.27642

16 ロット2 5.48055 0.20769 5.38403 5.58331 0.15307 0.30443

17 モデル

-対数尤度

χ

2

値

自由度

ｐ値

要因

推定値

標準誤差

モデルの差

1.34425232

2.6885

1 0.1011

ロット1

5.2596613

0.1052976

フルモデル

-2.8245958

ロット2

5.3695364

0.0657343

縮約モデル

-1.4803435

σ

0.1976291

0.0245172

対数正規分布を仮定したときのパラメータの推定値

対数正規分布を仮定したときのロット間差の検定

補足；Ｉ

ＣのＥＭ寿命の

ロットによる層別

(111)

111

(112)

競合リスクモデル

n

イベント

(故障,死亡)の原因が複数考えら

れる場合に発生原因別に解析したい

n

得られたデータはイベントの生起時間のみ

n

イベントは，原因Ａ，Ｂ，Ｃ…でおきる

n

原因Ａ，Ｂ，Ｃ，…は背反で生存時間に依

存しない

n

原因Ａでイベントがおきたとき，

原因Ｂ，Ｃ，…でイベントは起きない

(113)

補足；

ジャーナル部分の故障

T

4 :シャフトの変形

T

3 :疲れ破壊

T

２ :ジャーナル部の摩耗

T

１ :スプラインの変形

原因別シャフトの故障

(114)

競合リスクモデルの解析の仕方

n

着目した原因で故障，死亡が発生した場合はイ

ベントが発生したと考える

n

着目した原因以外で故障，死亡が発生した場合

は，観測が打ち切られたと考えて打ち切りデータ

として扱う

n

故障の原因が同じようなものは

1つにまとめる（群

間の検定をおこなう）

(115)

補足；

競合リスクモデル

n

競合リスクモデル

n

リスクとは故障原因となりうるもの

n

リスクが同時に複数存在しているモデル

n

競合リスクモデルの下でのデータ

n

最小値とその値が観測される原因となったリ

スク対ﾃﾞｰﾀで

,ﾗﾝﾀﾞﾑ打切りデータとなる

n

競合リスクモデルの目的

n

リスクが除去された場合の効果の把握

n

故障原因ごとの寿命特性の解析

(116)

補足；

故障を層別して解析

n

故障モード

(リスク)を層別して解析する

n

市場の製品は複数のサブシステムから構成

n

故障のモードも複数ある

n

それらは競合リスクモデルであり

n

ドレニックの定理により

,m=1(指数分布)として

推定される

n

そこから得られる情報は僅かであり

,故障のモー

ドにより層別することが重要である

(117)

補足；

競合リスクモデルによる

打切りデータ

■は故障モードＡ ▲は故障モードＢ 1 2 3 4 5 6 7 0 50 100 150 200 動作時間アイテム番号 ■はモードＡの故障時刻 ○はモードＡの中途打切り 1 2 3 4 5 6 7 0 50 100 150 200 動作時間アイテム番号

ＪＭＰ V4 による生存時間分析

ＪＭＰ

V4 による生存時間分析

∼生存時間分析の実行∼

於：SAS東京

2000.11.18

（

株）

リコー ＣＳ・

品質本部

廣野 元久

１．はじめに

n

今回の生存時間分析のデータでは，

主に工業製品での信頼性試験データを扱う

例題は主に実験データ，要因の水準(値)

n

生存時間分析の基礎

生存関数，ハザード関数，故障（死亡）の分類，

Weibull分布，打切りデータ，競合リスク

n

生存時間の分析（一変量の解析）

カプラン・マイヤー法，確率プロット

n

非線形回帰分析（要因解析）

比例ハザードモデル，加速モデル(Weibull回帰)

反応論モデル

本セミナーの内容

生存時間分析とは

n

ある基準となる時刻から，

目的の反応が得られるまでの

(生存)時間データ(Survival Time)を

対象とした解析手法の総称

n

工業分野では信頼性データ解析と呼ばれる

目的の反応とは

n

観測する個体に対して

1度だけ

非再起的に発生する事象

(Event)

n

例）

死亡

癌や循環器系の臨床研究での患者の再発・死亡

故障

システムや機器の信頼性研究での故障

生存時間データの特徴

n

時間データである（

非負、長期の観測）

n

観測打ち切りデータがある

n

故障や死亡の原因が複数ある

n

加速性や比例ハザード性がある

n

どの分布を仮定するかがあらかじめ不明

n

正規分布の理屈が使いにくい

n

データ収集に時間とお金がかかる

Time

Start of Study

End

Lost

Died

Time

Censor

Censor

Died

Start of Observation

Died

Died

Died

Died

生存時間データ

Start of Study

リコーＣＳ・

廣野元久