NLMIXED プロシジャを用いた生存時間解析伊藤要二アストラゼネカ株式会社臨床統計プログラミンググループグルプ Survival analysis using PROC NLMIXED Yohji Itoh Clinical Statistics & Programming Group, A

Size: px

Start display at page:

Download "NLMIXED プロシジャを用いた生存時間解析伊藤要二アストラゼネカ株式会社臨床統計プログラミンググループグルプ Survival analysis using PROC NLMIXED Yohji Itoh Clinical Statistics & Programming Group, A"

たつぞうはらしない
7 years ago
Views:

1 NLMIXED プロシジャを用いた生存時間解析伊藤要二アストラゼネカ株式会社臨床統計プログラミンググループグルプ Survival analysis using PROC NLMIXED Yohji Itoh Clinical Statistics & Programming Group, AstraZeneca KK

2 要旨 : NLMIXEDプロシジャの最尤推定の機能を用いて指数分布 Weibull 加速モデル Cox 比例ハザードモデルの場合を例として生存時間解析を試みたそして LIFEREGプロシジャやPHREGプロシジャと同等な解析結果を得ることが可能であることを示したこのことを通して生存時間解析の習得及び応用において NLMIXED プロシジャが強力な道具となり得ることを示したキーワード : 最尤法指数分布 Weibull 加速モデル Cox 比例ハザードモデル 2

3 目的生存時間解析のような複雑な統計解析でも SAS が全部やってくれる SAS の計算内容の詳細を必ずしも理解する必要がないその理論や計算アルゴリズムがなかなか身に付かない自分で IML などを用いてプログラムを書くことが望ましい発表の目的 : でも最尤法のアルゴリズムを自分で書くのはかなり面倒 NLMIXED プロシジャが一般的な最尤推定の計算ルーチンを提供幾つかのモデルを例として生存時間解析のプログラムをNLMIXEDプロシジャを用いて書くことができることを示す NLMIXEDプロシジャが生存時間解析を含む統計計算において非常に強力な道具であることを理解してもらうこと 3

4 生存時間解析における最尤推定指数分布 Weibull 分布及びその別の形 Weibull 加速モデル発表内容 Cox 比例ハザードモデルタイのない場合 ( タイのある場合は論文参照 ) まとめ 4

5 生存関数及び確率密度関数生存時間 T が連続の場合のみを考える確率密度関数 Probability density function: 生存関数 Survival function: ここで θ は分布のパラメタ f (t; θ ), t 0 St ( ; θ ) = PT ( > t ; θ ) = 1 fs ( ; θ ) ds 0 t 5

6 データ及び尤度データ生存時間 : t 1, t 2, t 3,,t n ( 打切り例の場合は最終生存確認時間 ) イベント情報 : δ 1, δ 2, δ 3,,δ n, イベントありなら δ i =1 打切りなら δ i = 0 i 番目の症例の尤度への寄与 L i イベント発現例 : 確率密度 : L i = f (t i ; θ ) 打ち切り例 : t i まで生存したという情報を利用 : L i = S (t i ; θ ) 尤度 : これらを全症例分かけ合わせて n n δi 1 δi L( θ) = [ f( t ; θ) S( t ; θ) ] 対数尤度 : これを対数変換して n i i i i= 1 i= 1 n l ( θ ) = log L ( θ ) = [ δ log f ( t ; θ ) + (1 δ )log S ( t θ )] i i i i i i= 1 i= 1 6

7 生存時間の生存関数最尤推定確率密度関数の2つさえ定義されれば対数尤度が定義されあとは最尤法を適用するだけ NLMIXEDプロシジャを利用できる本来は非線形混合効果モデルのためのプロシジャしかし最尤推定のための一般的な機能が利用可能能 NLMIXEDプロシジャのMODELステートメントに MODEL dependent-variable ~ general(ll) を指定ここで LL はユーザー定義の対数尤度関数 7

8 指数分布 Exponential distribution ハザード ( 定数 ): ht ()= λ 生存関数 : St () = e λt 確率密度関数 : ft () = λe λt ( 以下関数の括弧の中のパラメタは省略 ) 8

9 PROC NLMIXED による指数分布の当てはめ df=1e8 は t 分布ではなく正規分布に基づく推測をするために付加 λt ft () = λe St λ () = t 対数尤度 e9

10 指数分布 : PROC NLMIXED からの出力 10

11 Weibull 分布ハザード関数 : ht () = λγ 1 t γ λ: scale parameter, γ : shape parameter 生存関数 : St () = exp( λt γ ) 確率密度関数 : ft t t γ 1 γ ( ) = λγ exp( λ ) γ =1の場合は h (t) = λ となりこれは指数分布となる 11

12 Weibull 分布の当てはめ PROC NLMIXED γ 1 γ ft ( ) = λγ t exp( λt ) St () = exp( λ t γ ) Output 12

13 Weibull 分布の別の形教科書によく載っているWeibull 分布の生存関数 : St ( ) = exp( λ t γ ) (1) この形では拡張性に乏しいので変数およびパラメタの変換を行う : すなわち w = log( t), σ = 1/ γ, μ = (log λ) / σ t = exp( w ), γ = 1/ σ, λ = exp( μ / σ ) これらを (1) に代入し対数生存時間 w を用いた生存関数は w μ SW ( w) = exp{ exp( μ / σ)[exp( w)] 1/ σ } = exp[ exp( )] σ これより w の確率密度関数は f W 1 w w ( w) exp( μ μ = )exp[ exp( )] σ σ σ 13

14 Weibull 分布の別の形の当てはめ PROC NLMIXED 1 w μ w μ fw w = σ σ σ w μ SW ( w) = exp[ exp( )] σ f ( ) exp( )exp[ exp( )] PROC LIFEREG 14

15 Weibull 分布の別の形 : 結果の比較 PROC NLMIXED からの Output PROC LIFEREG からの Output 15

16 加速モデル Accelerated failure time (AFT) model; 一般論 2 つの治療群がありそれぞれの生存関数を S 0 (t) および S 1 (t) と表わすとするこの 2 つの生存関数の間に S () t = S ( t a), a> という関係が成り立つと仮定する a は加速パラメタ通常は a =exp( β ) とおき先ほどの式は S () t = S ( t exp( - β )), - < β < 1 0 となり β が推定すべきパラメタとされる治療群を表わす変数を z とし z =0は標準治療群 z =1 は新規治療群を表わすものとすると zが与えられた場合の生存関数は Stz ( ; ) = S ( t exp( - β z )) = S ( exp( log t - βz )), z= 0,

17 Weibull 分布を仮定した場合 Weibull 加速 (AFT) モデル先ほどは w=log t によって変数変換した w についての生存関数確率密度関数を考えたがその代わりにそれに群の効果を線形の形で付加した w = log t β z, z =0,1 を用いると先述のものと全く同じ生存関数確率密度関数を用いることができる 17

18 Weibull AFT model の数値例 18

19 Weibull AFT model の数値例 PROC NLMIXED PROC LIFEREG 先述のプログラムに先述のプログラムのモデルにが追加されただけ 19

20 Weibull AFT model の結果の比較 PROC NLMIXED からの Output PROC LIFEREG からの Output 20

21 データ Cox 比例ハザード (PH) モデル r 人の死亡例における生存時間を小さい順に並べる生存時間 t (1) t (2) t (3) t (r) 治療群を表わす変数 z (1) z (2) z (2) z (r) ここで z (i) = 0 は標準治療群 z (i) = 1 は新規治療群を表わす生存時間に同順位 ( タイ ) は存在しないものと仮定する 21

22 部分尤度 partial likelihood ( 式の誘導等の詳細な説明は省略しますが ) PL r = = i = 1 j R t exp( β ) ( i ) z() i exp( βz ) ( ) j ここで r は全イベント数 t (i) は (i) 番目のイベント発現時点 z (i) は (i) 番目の時点にイベント発現した症例の群の情報 β は群の効果を表わすパラメタ R(t (i) ) は時点 t (i) におけるリスクセット ( 時点 t (i) においてリスクに曝されている症例の集合 ) 22

23 前出の数値例 Cox PH モデルの数値例リスクセット z i t i δ i (i) z (i) t (i) δ (i) のサイズ t i の順に並べ替え

24 部分尤度の計算 t (1) =2の死亡例の部分尤度への寄与 β 0 (i) t (i) δ (i) z (i) e PL = e + e + e + e + e + e + e + e + e + e t (2) = 3 の死亡例の部分尤度への寄与 β e PL(2) = β 1 β 0 β 1 β 0 β 1 β 0 β 0 β 1 β 1 e + e + e + e + e + e + e + e + e PL (1) β 0 β 1 β 0 β 1 β 0 β 1 β 0 β 0 β 1 β t (3) =5ではイベント発現がないので無視 t (4) =6の死亡例の部分尤度への寄与 β 1 β e PL = e + e + e + e + e + e + e (4) β 1 β 0 β 1 β 0 β 0 β 1 β 1 24

25 事前のデータ加工 25

26 各行は各イベント発現時点に対応加工後のデータセット z1~z10は各時点におけるリスクセットに含まれる各症例の群の情報 26

27 Cox PH モデルの PROC NLMIXED によるあてはめ PL r = i exp( β ) z() i exp( 1 β z ) ( ) = j R t リスクセットに含まれる例数分だけ回す ( i ) 部分尤度の分母部分 j 27

28 Cox PH モデルの結果の比較 PROC NLMIXEDからのOutput PROC PHREG からの Output 28

29 同順位 ( タイ ) のある場合 Cox 比例ハザードモデルにおいて同順位が存在する場合には幾つかの計算法がある Breslowの近似を用いる場合には多少プログラムが複雑になるものの NLMIXEDプロシジャを用いて計算可能である ( 論文参照 ) 29

30 まとめパラメトリックな生存時間解析においては確率密度関数及び生存関数さえ定義すれば NLMIXED プロシジャを用いて最尤推定が可能 SAS のプロシジャを用いての生存時間解析もその計算内容を理解確認することができブラックボックスでなくすることができる LIFEREG プロシジャがサポートしていないような生存時間分布でもその確率密度関数及び生存関数さえ定義できれば推定可能更に混合分布モデルや Piecewise exponential model などのより複雑なモデルにも応用可能 Cox 比例ハザードモデル解析で示されたように事前のデータ加工等の工夫を加えることによってより複雑な計算アルゴリズムの最尤推定であっても問題によっては計算できる可能性が示唆されたこのようにPROC NLMIXEDが生存時間解析において非常に強力な道具であり更なる応用の可能性が示唆された 30

様々なミクロ計量モデル†

様々なミクロ計量モデル† 担当 : 長倉大輔 ( ながくらだいすけ ) この資料は私の講義において使用するために作成した資料です WEB ページ上で公開しており自由に参照して頂いて構いませんただし内容について一応検証してありますがもし間違いがあった場合でもそれによって生じるいかなる損害不利益について責任を負いかねますのでご了承ください間違いは発見次第継続的に直していますがまだ存在する可能性があります 1 カウントデータモデル

NLMIXED プロシジャを用いた生存時間解析 伊藤要二アストラゼネカ株式会社臨床統計 プログラミング グループグルプ Survival analysis using PROC NLMIXED Yohji Itoh Clinical Statistics & Programming Group, A

NLMIXED プロシジャを用いた生存時間解析伊藤要二アストラゼネカ株式会社臨床統計プログラミンググループグルプ Survival analysis using PROC NLMIXED Yohji Itoh Clinical Statistics & Programming Group, A