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

NLMIXEDプロシジャを用いた生存時間解析

NLMIXEDプロシジャを用いた生存時間解析

伊藤要二 アストラゼネカ株式会社 臨床統計 プログラミング グル プ 臨床統計・プログラミング・グループ

S

i

l

l

i

i

Survival analysis using

PROC NLMIXED

PROC NLMIXED

Yohji Itoh

Clinical Statistics & Programming Group, AstraZeneca KK

要旨：

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

目的

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

発表内容

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

生存関数及び確率密度関数

生存時間Tが連続の場合のみを考える。

• 確率密度関数 Probability density function: • 確率密度関数 Probability density function:

f (t;θ ), t  0 • 生存関数 Survival function: = > =



( ; ) ( ; ) 1 ( ; ) t S t θ P T t θ f s θ ds ここでθは分布のパラメタ = > = −



0 ( ; ) ( ; ) 1 ( ; ) S t θ P T t θ f s θ ds ここでθは分布のパラメタ

データ及び尤度

デ タ • データ • 生存時間 : t₁, t₂, t₃, …, t_n (打切り例の場合は最終生存確認時間) イベ ト情報 δ δ δ δ • イベント情報: δ₁, δ₂, δ₃, …, δ_n , イベントありなら、 δ_i = 1 打切りなら δ 0 打切りなら、 δ_i = 0 • i 番目の症例の尤度への寄与 L_i • イベント発現例： 確率密度: L_i = f (t_i ; θ ) • 打ち切り例： t_i まで生存したという情報を利用: L_i = S (t_i ; θ ) • 尤度： これらを全症例分かけ合わせて i i n n i i i L ( ) = [ ( ; )f t S t( ; )1− ]

∏

θ

∏

θ δ θ δ • 対数尤度： これを対数変換して i=1 i=1 n n l( )θ



l L_i( )θ



[ lδ_i f t( _i θ) (1 δ_i )l S t( _iθ)]

最尤推定

• 生存時間の • 生存関数数 • 確率密度関数 の２つさえ定義されれば、対数尤度が定義され、あとは最尤法を適用する だけ • NLMIXEDプロシジャを利用できる • 本来は非線形混合効果モデルのためのプロシジャ。 か 最尤推定 ため 般的な機能が利 能 しかし、最尤推定のための一般的な機能が利用可能 • NLMIXEDプロシジャのＭＯＤＥＬステートメントに

MODEL dependent-variable ~ general(ＬＬ)

指数分布 Exponential distribution

• ハザード（定数）:

h t

( )

=

λ

• ハザード（定数）:

生存関数

h t

( )

=

λ

• 生存関数:

t

S t

( )

=

e

−λ

• 確率密度関数:

t

f t

( )

=

λ

e

−λ

（

以下、関数の括弧の中のパラメタは省略）

PROC NLMIXEDによる指数分布の当てはめ

data ex1; input t event; cards; 3 1 6 0 6 0 12 1 24 1 「df=1e8」は、ｔ分布ではなく、正規分布 に基づく推測をするために付加 48 1 ;

proc nlmixed data=ex1 df=1e8; proc nlmixed data=ex1 df=1e8;

parms lambda=1.;

if event=1 then lh=lambda*exp(-lambda*t); ← f t( ) = λe−λt

else lh= exp(-lambda*t); ll=log(lh); model t general(ll); ← ← 対数尤度 ( ) = −λ t S t e model t general(ll);

指数分布： PROC NLMIXEDからの出力

Parameter Estimates Parameter Estimates

Standard

Parameter Estimate Error DF t Value Pr > ￨t￨ Parameter Estimate Error DF t Value Pr > ￨t￨ lambda 0.04301 0.02151 1E8 2.00 0.0455 lambda 0.04301 0.02151 1E8 2.00 0.0455

Weibull分布

• ハザード関数: λ: scale parameter, 1

( )

h t

=

λγ

t

γ− p , γ : shape parameter • 生存関数:

S t

( )

=

exp(

−

λ

t

γ

)

• 確率密度関数:

f t

( )

=

λγ

t

γ −1

exp(

−

λ

t

γ

)

γ = 1の場合は、h (t) = λとなり、これは指数分布となる。 γ 場合 、 ( ) り、 指数分布 。

Weibull分布の当てはめ

PROC NLMIXED

proc nlmixed data=ex1 df=1e8; parms lambda=1 gamma=1;

PROC NLMIXED

parms lambda=1 gamma=1;

if event=1 then lh=lambda*gamma*t**(gamma-1) *exp(-lambda*t**gamma); else lh=exp(-lambda*t**gamma); S t( ) ( λtγ ) ← f t( ) = λγtγ−1 exp(−λtγ ) else lh=exp(-lambda*t**gamma); ll=log(lh); model t general(ll); S t( ) = exp(−λtγ ) ← run; Parameter Estimates Output Parameter Estimates Standard

Parameter Estimate Error DF t Value Pr > ￨t￨ Parameter Estimate Error DF t Value Pr > ￨t￨ lambda 0.01612 0.02856 1E8 0.56 0.5724 gamma 1.2916 0.4950 1E8 2.61 0.0091 gamma 1.2916 0.4950 1E8 2.61 0.0091

Weibull分布の別の形

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

S ( )w = exp{ exp(− −μ σ/ )[exp( )]w 1/σ } = exp[ exp(− − μ )]

σ

• これよりw の確率密度関数は

W

w w

f ( )w = 1 exp( − μ )exp[ exp(− − μ )]

σ σ σ

Weibull分布の別の形の当てはめ

proc nlmixed data=ex1 df=1e8;

parms mu=1 sigma=1; f ( ) 1 (w − μ ) [ (w − μ )]

PROC NLMIXED

parms mu=1 sigma=1;

w=log(t);

if event=1 then lh=1/sigma*exp((w-mu)/sigma) *exp(-exp((w-mu)/sigma));

1

W

w w

f ( )w = exp( μ )exp[ exp(− μ )]

σ σ σ *exp(-exp((w-mu)/sigma)); else lh=exp(-exp((w-mu)/sigma)); ll=log(lh); ll=log(lh); model t general(ll); run; W w S ( )w = exp[ exp(− − μ)] σ

proc lifereg data=ex1; PROC LIFEREG

proc lifereg data=ex1;

model t*event(0)= /d=weibull; run;

Weibull分布の別の形： 結果の比較

Parameter Estimates

PROC NLMIXEDからのOutput

Standard

Parameter Estimate Error DF t Value Pr > ￨t￨ mu 3.1956 0.3976 1E8 8.04 <.0001 sigma 0.7742 0.2967 1E8 2.61 0.0091

Analysis of Parameter Estimates

PROC LIFEREGからのOutput

Standard 95% Confidence

Chi-Parameter DF Estimate Error Limits Square Pr > ChiSq Intercept 1 3.1956 0.3976 2.4163 3.9749 64.60 <.0001 Intercept 1 3.1956 0.3976 2.4163 3.9749 64.60 <.0001 Scale 1 0.7742 0.2967 0.3653 1.6407

加速モデル

Accelerated failure time (AFT) model; 一般論

• 2つの治療群があり、それぞれの生存関数をS₀(t) およびS₁(t) と表わすと する。この２つの生存関数の間に という関係が成り立つと仮定する。a は加速パラメタ。 = ⋅ > 1( ) 0( ), 0 S t S t a a • 通常は、 a = exp(−β ) とおき、先ほどの式は ( ) ( exp( )) S t S t β ¥ < < ¥β となり、β が推定すべきパラメタとされる。 治 変 1( ) 0( exp( )), S t = S t ⋅ -β - ¥ < < ¥β • 治療群を表わす変数をz とし、 z = 0は標準治療群 z =1は新規治療群 z =1は新規治療群 を表わすものとすると、zが与えられた場合の生存関数は 0 0

( ; ) ( exp( )) (exp(logt )), 0, 1

S t z( ; ) = S t₀( exp(⋅ -βz)) = S₀(exp(logt - βz)), z = 0, 1

Weibull加速（AFT）モデル

• Weibull分布を仮定した場合

先ほどは、 w = log t によって変数変換したw についての生存関数、確率密度関数を考えたが、 その代わりに、それに群の効果を線形の形で付加した その代わりに、それに群の効果を線形の形で付加した w = log t −β z, z = 0,1 を用いると、先述のものと全く同じ生存関数、確率密度関数を用いることが できる できる。

Weibull AFT modelの数値例

data ex2;

input group t event; cards; cards; 0 2 1 0 5 0 0 5 0 0 11 1 0 15 1 0 21 1 0 21 1 1 3 1 1 6 1 1 6 1 1 12 0 1 24 1 1 48 1 1 48 1 ;

Weibull AFT modelの数値例

proc nlmixed data=ex2 df=1e8; PROC NLMIXED

先述のプログラムに

proc nlmixed data=ex2 df=1e8; parms mu=1 sigma=1 beta=0; w=log(t) - beta * group;

if event=1 then lh=1/sigma*exp((w-mu)/sigma)

先述のプログラムに 「_{- beta * group」} 及びbetaの初期値が 追加されているだけ

if event=1 then lh=1/sigma*exp((w-mu)/sigma) *exp(-exp((w-mu)/sigma)); else lh=exp(-exp((w-mu)/sigma)); 追加されているだけ で、他は全く同じ else lh=exp(-exp((w-mu)/sigma)); ll=log(lh); model t general(ll); run; run; PROC LIFEREG 先述のプログラムの

proc lifereg data=ex2;

model t*event(0)=group /d=weibull; run;

先述のプログラムの モデルに_groupが追 加されただけ

Weibull AFT modelの結果の比較

Standard

PROC NLMIXEDからのOutput

Parameter Estimate Error DF t Value Pr > ￨t￨ mu 2.6155 0.3665 1E8 7.14 <.0001 sigma 0.7315 0.2065 1E8 3.54 0.0004 beta 0.5940 0.5196 1E8 1.14 0.2530

Standard 95% Confidence

Chi-Parameter DF Estimate Error Limits Square Pr > ChiSq

PROC LIFEREGからのOutput

Parameter DF Estimate Error Limits Square Pr > ChiSq Intercept 1 2.6155 0.3665 1.8973 3.3338 50.94 <.0001 group 1 0.5940 0.5196 -0.4244 1.6123 1.31 0.2530 group 1 0.5940 0.5196 -0.4244 1.6123 1.31 0.2530 Scale 1 0.7315 0.2065 0.4206 1.2722

Cox比例ハザード(PH)モデル

データ

r 人の死亡例における生存時間を小さい順に並べる

r 人の死亡例における生存時間を小さい順に並べる

生存時間

t

₍₁₎

t

₍₂₎

t

₍₃₎



t

_{( )}

生存時間

t

₍₁₎

t

₍₂₎

t

₍₃₎

t

_(r)

治療群を表わす変数

z

₍₁₎

z

₍₂₎

z

₍₂₎



z

_(r)

ここで、z

_(i)

= 0は標準治療群、 z

_(i)

= 1は新規治療群を表わす。

ここで、z

_(i)

0は標準治療群、 z

_(i)

1は新規治療群を表わす。

生存時間に同順位（タイ）は存在しないものと仮定する。

部分尤度 partial likelihood

（式の誘導等の詳細な説明は省略しますが） ( ) exp( ) r i z PL =

∏

β

ここで r は全イベント数 1 ( ) ( _i )exp( ) i _{j R t j} PL z

β

= _∈ =

∏ 

ここで、r は全イベント数 ｔ_(i)は(i)番目のイベント発現時点 ｚ_(i)は(i)番目の時点にイベント発現した症例の群の情報 β は群の効果を表わすパラメタ、 β R(t_(i))は時点t_(i)におけるリスクセット (時点t においてリスクに曝されている症例の集合) (時点t_(i)においてリスクに曝されている症例の集合)

Cox PHモデルの数値例

前出の数値例

z_i t_i δ_{i (i) z}

(i) t(i) δ(i)

リスクセット のサイズ

i i i

0 2 1

0 5 0

(i) z_(i) t_(i) δ_(i)

1 0 2 1 10 2 1 3 1 9 のサイズ 0 11 1 0 15 1 2 1 3 1 9 3 0 5 0 8 4 1 6 1 7 0 21 1 1 3 1 4 1 6 1 7 5 0 11 1 6 6 1 12 0 5 t_i の順に並べ替え 1 6 1 1 12 0 5 7 0 15 1 4 8 0 21 1 3 1 24 1 1 48 1 9 1 24 1 2 10 1 48 1 1

部分尤度の計算

• t₍₁₎ = 2の死亡例の部分尤度への寄与 (i) t(i) δ(i) z(i)

0

e PL

β⋅

• t₍₂₎ = 3の死亡例の部分尤度への寄与 (i) t(i) δ(i) z(i)

1 2 1 0 2 3 1 1 (1) 0 1 0 1 0 1 0 0 1 1 PL eβ⋅ eβ⋅ eβ⋅ eβ⋅ eβ⋅ eβ⋅ eβ⋅ eβ⋅ eβ⋅ eβ⋅ = + + + + + + + + + t₍₂₎ = 3の死亡例の部分尤度への寄与 2 3 1 1 3 5 0 0 4 6 1 1 1 (2) 1 0 1 0 1 0 0 1 1 e PL e e e e e e e e e β β β β β β β β β β ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ = + + + + + + + + • t₍₃₎ = 5ではイベント発現がないので無視 4 6 1 1 5 11 1 0 6 12 0 1 • t₍₄₎ = 6の死亡例の部分尤度への寄与 6 12 0 1 7 15 1 0 8 21 1 0 eβ⋅1  8 21 1 0 9 24 1 1 1 ( 4) 1 0 1 0 0 1 1 e PL e e e e e e e β β⋅ β⋅ β⋅ β⋅ β⋅ β⋅ β⋅ = + + + + + +  10 48 1 1

事前のデータ加工

data ex2;

input group t event;

proc iml;

sort ex2 out=s by t; use s; event; cards; 0 2 1 0 5 0 use s; read all; nrec=nrow(t); vname="i"￨￨"time"￨￨(repeat("z",nrec)` 0 5 0 0 11 1 0 15 1 vname="i"￨￨"time"￨￨(repeat("z",nrec)` +left(char(1:nrec,2,0))); vec=j(1,nrec+2,0);

create phdata from vec [colname=vname]; 0 15 1

0 21 1 1 3 1 1 6 1

create phdata from vec [colname=vname]; do i=1 to nrec;

if event[i]=1 then do; 1 6 1

1 12 0 1 24 1 1 48 1

if event[i]=1 then do; vec=i￨￨t[i]￨￨group`; append from vec;

end; 1 48 1 ; run; end; group[i]=.; end; close phdata;

加工後のデータセット

• 各行は各イベント発現時点に対応 • z1～z10は各時点におけるリスクセットに含まれる各症例の群の情報情 i time z1 z2 z3 z4 z5 z6 z7 z8 z9 z10 1 2 0 1 0 1 0 1 0 0 1 1 2 3 . 1 0 1 0 1 0 0 1 1 2 3 . 1 0 1 0 1 0 0 1 1 4 6 . . . 1 0 1 0 0 1 1 5 11 . . . . 0 1 0 0 1 1 5 11 . . . . 0 1 0 0 1 1 7 15 . . . 0 0 1 1 8 21 . . . 0 1 1 9 24 . . . 1 1 10 48 . . . 1

Cox PHモデルのPROC NLMIXEDによるあてはめ

proc nlmixed data=phdata df=1e8;

parms beta 0; ( )

exp(

)

(

)

r i

z

PL

β

β

=

∏ 

parms beta 0; array z[10] z1-z10; denominator=0; ( ) 1 _{( )}

exp(

)

i i= _{j R t∈}

β

z

_j

∏ 

denominator=0; do j=i to 10; denominator=denominator+exp(beta*z[j]); ← リスクセットに含まれる例数分だけ回す ←部分尤度の分母部分 denominator=denominator+exp(beta*z[j]); end; logpl=log(exp(beta*z[i])/denominator); 部分尤度の分母部分 logpl=log(exp(beta*z[i])/denominator); model time general(logpl);

run; run;

Cox PHモデルの結果の比較

Parameter Estimates

PROC NLMIXEDからのOutput

Standard

Parameter Estimate Error DF t Value Pr > ￨t￨ beta -0.9031 0.8724 1E8 -1.04 0.3006

Analysis of Maximum Likelihood Estimates

PROC PHREGからのOutput

Analysis of Maximum Likelihood Estimates

Parameter Standard Hazard Variable DF Estimate Error Chi-Square Pr > ChiSq Ratio Variable DF Estimate Error Chi-Square Pr > ChiSq Ratio group 1 -0.90307 0.87240 1.0715 0.3006 0.405

同順位(タイ）のある場合

• Cox比例ハザードモデルにおいて同順位が存在する場合には幾つかの 計算法がある。

• Breslowの近似を用いる場合には、多少プログラムが複雑になるものの、 NLMIXEDプロシジャを用いて計算可能である（論文参照）。

まとめ

• パラメトリックな生存時間解析においては、確率密度関数及び生存関数さ え定義すれば、NLMIXEDプロシジャを用いて最尤推定が可能 • SASのプロシジャを用いての生存時間解析もその計算内容を理解・ 確認することができ、ブラック・ボックスでなくすることができる。 • LIFEREGプロシジャがサポートしていないような生存時間分布でも、 その確率密度関数及び生存関数さえ定義できれば、推定可能

• 更に、混合分布モデルやPiecewise exponential modelなどのより複 雑なモデルにも応用可能 • Cox比例ハザードモデル解析で示されたように、事前のデータ加工等の 工夫を加えることによって、より複雑な計算アルゴリズムの最尤推定で あっても 問題によっては計算できる可能性が示唆された あっても、問題によっては計算できる可能性が示唆された。 • このようにPROC NLMIXEDが生存時間解析において非常に強力な道具 であり 更なる応用の可能性が示唆された であり、更なる応用の可能性が示唆された。