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, AstraZeneca KK

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

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

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

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

データ及び尤度 データ 生存時間 : 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

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

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

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

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

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

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

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

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

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

加速モデル Accelerated failure time (AFT) model; 一般論 2 つの治療群があり それぞれの生存関数を S 0 (t) および S 1 (t) と表わすとする この 2 つの生存関数の間に S () t = S ( t a), a> 0 1 0 という関係が成り立つと仮定する 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, 1 0 0 16

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

Weibull AFT model の数値例 18

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

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

データ 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

部分尤度 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

前出の数値例 Cox PH モデルの数値例 リスクセット z i t i δ i (i) z (i) t (i) δ (i) のサイズ 0 2 1 1 0 2 1 10 0 5 0 2 1 3 1 9 0 11 1 3 0 5 0 8 0 15 1 4 1 6 1 7 0 21 1 5 0 11 1 6 t 1 3 1 i の順に並べ替え 6 1 12 0 5 1 6 1 7 0 15 1 4 1 12 0 8 0 21 1 3 1 24 1 9 1 24 1 2 1 48 1 10 1 48 1 1 23

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

事前のデータ加工 25

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

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

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

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

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