スライド 1

Transcription

1 SAS による生存時間分布の予測 Death Note の統計学 東京理科大学 浜田知久馬, 魚住龍史 Death Note Predcton of Survval Dstrbuton usng SAS Statstcs for Death Note Chkuma Hamada Toko Unverst of Scence

2 要旨 分位点回帰を行うプロシジャ QUANTREG プロシジャを打ち切りがある場合に拡張した QUANTLIFE プロシジャがバージョン 9.4 から新たに追加された. このプロシジャを利用するとメディアン等の分位点と共変量の関係をモデル化することが出来る. 本発表では, チュートリアルとして LIFETEST プロシジャ,PHREG プロシジャ, LIFEREG プロシジャ,QUANTLIFE プロシジャを用いた生存関数の予測について解説する. キーワード : 生存時間分布, 比例ハザードモデル, 加速モデル, 分位点回帰 QUANTREG( 分位点回帰 )+LIFETEST( ノンパラ ) =QUANTLIFE( 分位点回帰の生存時間への拡張 ) 2

3 内容 Death Note 生存時間の予測ノンパラメトリックな予測 KM(LIFETEST) セミパラメトリックな予測 Co(PHREG) パラメトリックな予測加速 (LIFEREG) 分位点の予測 QUANTLIFE 3

4 生存時間解析の SAS の教科書 生存時間解析応用編 4

5 第 章はじめに. 生存時間解析とSAS.. LIFETEST プロシジャ ( 第 2 章 )..2 PHREG プロシジャ ( 第 3 章 )..3 POWER プロシジャ ( 第 4 章 )..4 LIFEREG プロシジャ..5 SURVEYPHREG プロシジャ..6 QUANTLIFE プロシジャ..7 ICLIFETEST プロシジャ..8 ICPHREG プロシジャ 5

6 第 2 章生存関数のノンパラメトリックな 推定と検定 (LIFETEST プロシジャ ) 2. ノンパラメトリックな生存関数と ハザード関数の推定 2.. カプラン マイヤー推定量 2..2 カプラン マイヤープロット 2..3 生存関数の信頼区間 2..4 生存関数の信頼バンド 2..5 ハザードの推定と可視化 群間の比較 多重比較法による解析 6

7 第 3 章コックス回帰によるハザード比の推定とその拡張 (PHREGプロシジャ) 3. PHREG プロシジャによる様々な線型仮説に対する検討 3.2 共変量及び多重性の調整 3.3 最大対比法の適用 3.4 モデルの評価 3.4. 残差統計量 累積ショーンフェルド残差プロットによる比例ハザード性の評価 累積マルチンゲール残差プロットによる線型モデルの仮定の評価 3.5 フレイルティモデルと周辺コックスモデルによる クラスター生存時間データの解析 7

8 第 4 章生存時間解析における 例数設計 (POWER プロシジャ ) 4. 生存時間解析における例数設計の概要 : 4.2 フリードマンの方法とショーンフェルドの方法 4.3 POWER プロシジャによる生存時間解析の 例数設計 4.3. ラカトスの方法 POWER プロシジャによる実行例 フリードマンの方法とショーンフェルドの 方法による例数との比較 区分直線モデルによる例数設計 8

9 応用編に関連したチュートリアル 浜田知久馬 藤井陽介 (2003). 生存時間解析の症例数設計. 浜田知久馬 安藤英一 (2005). POWER プロシジャによる症例数設計. 浜田知久馬 (20). 生存時間解析再入門 生存時間解析のミステリーをひも解く. 浜田知久馬 (203). SAS 生存時間解析プロシジャの最新の機能拡張. 9

10 共変量 (covarate) 予後因子 (prognostc factor) 背景因子 (baselne factor) 生存時間に対する影響を調べたい変数 ) 治療法の違いを表す共変量 2) 患者固有の性質を表す共変量 : 年齢, 検査値 3) 外性的な共変量 ( 環境要因 ): 施設, 時代 時間非依存性共変量 : 性, 人種 時間依存性共変量 : 検査値, 治療法 ( tme dependent covarate ) 0

11 生存時間に影響を与える因子 性別 : 男性 80.5 歳 (8 位 ) 女性 86.8 歳 ( 位 ) (204) 疾患 : 糖尿病, 高血圧食生活 : 赤ワイン, 果物地域 : 沖縄県は長寿だった生活習慣 : 喫煙, 運動職業 : 医師, タクシードライバー, 関取遺伝子 ( ゲノム疫学 )

12 生存時間解析の方針 ノンパラメトリック (LIFETEST) 特定の分布を仮定しない推定 Kaplan-Meer 法, ログランク検定, 一般化ウイルコクソン検定 セミパラメトリック (PHREG) 特定の分布を仮定しないハザード比の推定 Co の比例ハザードモデル パラメトリック (LIFEREG) 特定の分布を仮定した生存時間解析 加速モデル 2

13 3 Kaplan-Meer 法死亡時点で階段が低下 t t 2 t 3 t 4 /7 /

14 Kaplan-Meer 法 瞬間の単位で生き残る確率を積算 時点 t まで生き残る :t <t の全てを生き残る. 死亡が起きてない時点 : 生き残る確率 = 死亡が起きた時点 : 生き残る確率 = n : 時点 のリスク集合の大きさ d : 時点 の死亡数死亡がおきてない時点は水平で, 死亡が起きた時点で生存関数 S(t) は階段状に低下 (d = で打ち切りがなければ,/n づつ低下 ) d n

15 5 生存割合の推移時点 t まで : ( 死亡がおきてない ) t ~t 2 : t 2 ~t 3 : t 3 ~t 4 : t 4 ~ :

16 死亡 : 階段低下打ち切り : ヒゲ t t 2 t 3 t 4

17 加速 加速 ( 時間に対するモデル化 ) 比例ハザード ( 時間の逆数のモデル化 ) T T ep( β T z 0 β>0:z T ( リスクは減少 ) β<0:z T ( リスクは増大 ) 比例 h( t) h0 ( t)ep( β T z ) β>0:z T ( リスクは増大 ) β<0:z T ( リスクは減少 ) 2つのモデル間でβの符号は反対 ) 7

18 8 加速モデルと比例ハザードモデル 加速モデル 比例ハザードモデル )) ( log log( )) ( log log( ) ( ) ( ) )ep( ( ) ( 0 ) ep( 0 0 t S t S t S t S t h t h T T T z β z β z β log log ) ep( ) ( ) ep( T T t S t S T T T T T z β z β z β

19 加速モデル : 横軸が定数倍,S(t)=S 0 (t/2) :T 0 :T ep(β T z)=2 S( ) S 0 S 任意の % 点が定数倍 9

20 比例ハザードモデル : 縦軸がべき乗 S(t)=S 0 (t) 2 ep(β T z)=2のとき

21 二重対数プロット 加速モデル log(-log(s)) 水平方向に 平行移動 log T z log T 0 β a z=0, log(tme) 2

22 二重対数プロット 比例ハザードモデル log(-log(s)) β ph 垂直方向に 平行移動 log( log S ( t )) log(tme) z log( log z=0, S 0( t )) 22

23 ワイブル確率プロット ( 二重対数プロット ) ワイブル分布の生存関数 log( log S( t)) S( t) log S( t) ep( t t log( log S( t)) log log t log t logtを横軸,log(-logs(t)) を縦軸にプロットすると傾きγ の直線になる. ) 切片 :log 傾き : 23

24 二重対数プロット ワイブル分布 log(-log(s)) 傾き :γ β a 切片 :log β ph β ph = γβ a log( log S( t)) log log t log(tme) 24

25 二重対数プロット ) 横軸 :log(tme) 縦軸 :log(-log(s)) 2) 比例ハザードモデルのとき, 上下に平行移動 3) 加速モデルのとき, 左右に平行移動 4) ワイブル分布のときは直線 ( 上下, 左右に同時に平行移動 ) 比例ハザードかつ加速モデル -β ( 加速 ) γ =β ( 比例 ) 25

26 例題データ : 加速モデル data test; do drug=0,; do = to 20; tme=*(drug+);output;end;end;run; proc prnt;run; proc sgplot; scatter =tme =drug ;run; 26

27 薬剤群 T 2T 0 例題データ : 人工的な加速モデル T 40 生存時間が 2 倍 対照群 T 生存時間 27

28 LIFETEST による解析 ( ノンパラ ) proc lfetest data= test plots=(s,lls); tme tme; strata drug; run; 28

29 生存関数の予測 ( ノンパラ ) Kaplan-Meer プロット 29

30 二重対数プロット ( ノンパラ ) T 2 T 0 logt log 2 logt 0 log(2)=0.693 対数生存時間比 log(tme) 30

31 PHREG による解析 ( セミパラ ) data dumm; do drug=0 to ;_group_=;output; end; proc phreg data=test plots(overla=group)=(s); model tme= drug; baselne out=out covarates=dumm survval=s loglogs=lls/nomean;run; data out;set out;logw=log(tme); label s='survval' lls='lls'; proc gplot data=out; plot s*tme=drug lls*logw=drug;;run; 3

32 生存関数の予測 ( セミパラ ) 32

33 二重対数プロット ( セミパラ ) log(-log(s)) 尤推定値の分析 パラメータ パラメータ推定値 ハザード比 drug 対数ハザード比 -.369=log(0.254) 共変量の効果を誤特定 33

34 LIFEREG による解析 ( パラ ) ワイブル分布 proc lfereg data=test; model tme=drug/dst=webull; output out=out cdf=cdf; data out;set out;s=-cdf; lls=log(-log(s));logw=log(tme); proc gplot; plot s*tme=drug lls*logw=drug; smbol =splne c=red w=3 l=; smbol2 =splne c=green w=3 l=;run; 34

35 生存関数の予測 ( パラ ) T 2T 0 35

36 log(-log(s)) 二重対数プロット ( パラ ) 生存時間分布を誤特定 =.262 log2=0.693 Analss of Mamum Lkelhood Parameter Estmates Parameter Estmate Intercept drug Scale Webull Shape

37 比例ハザードモデルと加速モデル ) ノンパラ共変量をモデル化できない. 既存のモデルは仮定が厳しすぎる. 2) 比例ハザードモデルハザードが時点によらず定数倍二重対数プロットが縦軸方向に平行移動 3) 加速モデル ( ワイブル分布 ) 時間の % 点が定数倍二重対数プロットが横軸方向に平行移動特定の生存時間分布を仮定する必要がある. 37

38 logt ( ) QUANTLIFE による解析 0( ) ( ) proc quantlfe data= test method=na log plot=quantplot; model tme=drug /quantle=0.00 to 0.95 b 0.05; baselne out=qlout covarates=dumm survval=qlsurv quantle=ql;run; proc gplot data=qlout;plot qlsurv*ql=drug; smbol =steplj w=3 l= c=red; z τ 分位点に対する回帰分析 加速 ( 対数線形 ) モデル smbol2 =steplj w=3 l=2 c=blue;run; 38

39 Parameter Estmates 加速モデル logt Quantle Parameter DF Estmate Standard Error 0( ) ( ) 95% Confdence Lmts t Value Pr > t Intercept drug Intercept drug Intercept <.000 drug Intercept <.000 drug Intercept <.000 drug Intercept <.000 drug Intercept <.000 drug Intercept <.000 drug < Intercept <.000 ( ) drug <.000 z 39

40 パラメータ推定値 Quantle logt ( ) 0( ) ( ) z ( ) log(2) log() = to 20 0( ) log(20 ) log(2)=0.693 % 点によらず係数は一定加速モデルが成立 plot=quantplot オプション F( t) S( t) 40

41 生存関数の予測 4

42 N=9 で対数変換なし ( 加法モデル ) data test; do drug=0,; do = to 9;tme=*(drug+);output;end;end;run; proc quantlfe data= test method=na plot=quantplot; model tme=drug / quantle=0.00 to 0.95 b 0.05; baselne out=qlout covarates=dumm survval=qlsurv quantle=ql; run; T 0( ) ( ) proc gplot data=qlout;plot qlsurv*ql=drug; smbol =steplj w=3 l= c=red; smbol2 =steplj w=3 l=2 c=blue;run; 42 z

43 生存関数の予測 43

44 加法モデル Parameter Estmates Quantle Paramete DF Estmate Standard 95% Confdence t Value Pr > t r Error Lmts Intercept drug Intercept drug Intercept drug Intercept drug Intercept <.000 drug Intercept <.000 drug Intercept <.000 drug Intercept <.000 drug Intercept <.000 T 0( ) ( ) drug <.000 z 44

45 パラメータ推定値 Quantle T 20 0( ) ( ) 20 0( ) ( ) z F( t) S( t) F( t) S( t) 45

46 QUANTLIFE の文法 PROC QUANTLIFE < optons > ; BASELINE < optons > ;* 予測する共変量の指定 ; CLASS varables ; EFFECT name = effect-tpe ( varables < / optons > ) ; MODEL response < censor(lst) > = < effects > < / QUANTILE= > ; OUTPUT <OUT=SAS-data-set > < keword=name... keword=name > ; TEST effects < / optons > ;

47 生存時間解析のプロシジャの比較 プロシジャ焦点モデル LIFETEST 生存関数 アプローチ 共変量打ち切り ノンパラ層別右側 PHREG ハザード関数 比例ハザードモデル セミパラ モデル調整 右側 LIFEREG 生存時間 加速モデル パラ モデル調整 右側, 左側, 区間 QUANTLIFE 生存時間の分位点 加速モデル 加法モデル セミパラ * モデル調整 右側 *QUANTLIFE: 生存時間分布に特定の分布を仮定しないが, 共変量で % 点が変化 47

48 最小 2 乗法と絶対値の和の最小 48 5 : mn ) 2( ) (9 ) (6 ) (5 ) (3 ) ( mn ) ( medan A A S A na A da ds A A A A A A S,,,, :

49 算術平均 :4.8 中央値 :5.0 2 乗和を最小絶対値の和を最小 A 49

50 最小 2 乗の確認 S=Σ(Y -4.8) 2 =(-4.8) 2 + (3-4.8) 2 + (5-4.8) 2 + (6-4.8) 2 + (9-4.8) 2 =(-3.8) 2 + (-.8) 2 + (0.2) 2 + (.2) 2 + (4.2) 2 =36.8 S =Σ(Y -5.0) 2 =(-5.0) 2 + (3-5.0) 2 + (5-5.0) 2 + (6-5.0) 2 + (9-5.0) 2 =(-4.0) 2 + (-2.0) 2 + (0) 2 + () 2 + (4) 2 =

51 最小絶対値の確認 S=Σ Y -4.8 = = =.2 S =Σ Y -5.0 = = =.0 5

52 絶対値の和の変化 A 52

53 最小 2 乗法の模式図 単回帰分析 Y no X X X 53

54 最小 2 乗推定直線による予測式 : 残差 : 残差平方和 : b 0 b の関数偏微分して 0 正規方程式 54 b b 0 e mn ) ( ) ( e b b e S

55 最小 2 乗推定量 55 e e e S S b b e e b db ds b b b S b b b b b nb e b b db ds ) )( ( ) )( ( ) )( ( ) )( ( 0 2 ) ( 2 )) ( )( ( 2 )) ( ( ) ( ) (

56 単回帰分析の結果 b 0.6 b 2.4

57 % 点 ( 分位点 ) 回帰の指標 57 n T n T n T n T n T n T T T T T S S S S β β β β β β β β β β : : : : : 95% ) ( : % 0.5 : : 点回帰点回帰メディアン回帰回帰分析

58 proc quantreg; model = /quantle= ; % 点 50% 点 2 5% 点

59 QUANTREG プロシジャの結果 2 切片 傾き 分位点 分位点

60 QUANTREG プロシジャの結果 50% 点 メディアン回帰 : S 0.5 β n T 2

61 5% 点 : 負の残差が出ないようにする S : n T β 0.05 T β : n T β 0.95 T β

62 95% 点 : 正の残差が出ないようにする n T n T T T S β β β β : :

63 分位点回帰の特徴 ノンパラメトリック回帰 ) 外れ値に対して頑健 ( 平均値に対する中央値の頑健性.) 2) 等分散性が成り立たない場合でも対応可 ( 期待値が大きくなると分散が大きくなる.) 3) % 点についてモデル化 % 点ごとに異なった回帰係数をあてはめることができる. 63

64 203 年プロ野球打者のデータ 年俸 ( 万円 ) 阿部 回帰分析 90% 予測区間 通算安打 64

65 年俸万円 203 年プロ野球打者のデータ 分位点回帰 ) 不等分散性に対応可 2) 外れ値に対して頑健 分位点回帰 95% 点 50% 点 5% 点 通算安打 65

66 % 点回帰 (quantlfe の記載 ) n T n T n T T n T T T I r u I u u β β β β β β β β : : ) ( 0) ) (( ) ( 0)) ( ( ) (

67 分位点回帰の 打ち切りデータに対する拡張方針 打ち切りを起こした個体 : 打ち切り時点から の どこかで死亡が起きる. ) 打ち切り時点で死亡 2) ( 最終観察時点 ) の時点で死亡 に重みを付けて分割 後の時点の累積死亡率を推定するときは,) の重みを大きくする. 重みは逐次的に更新 67

68 Kaplan-Meer-Tpe Estmator イベント, 打ち切りの寄与 T T T w or w w ) ( : : )) ( ( ) ( * * 最終観察時点打ち切り時点生存時間打ち切り : イベント : β β β : 推定したい分位点個体の累積分布打ち切り時点の :

69 打ち切りデータの重み 69 の死亡の重み時点 : 打ち切り時点の死亡の重み : 推定したい分位点打ち切り時点の個体の累積分布 : ) ( ) ( : w w

70 Kaplan-Meer-Tpe Estmator 打ち切りデータの重み : 0.4, : 0.4, : 0.4, : 0.4, で死亡最終観察時点で打ち切り : * 累積分布が 0.3 の時点では生存累積分布が.0 の時点では死亡

71 重みの模式図 時点 での死亡の重み 累積死亡分布関数 時点 t での死亡の重み t 時間 最終観察時点 7

72 重みの持つ意味 KM 推定量 : * * * 72 0.,0.2,0.3,0.4,0.52,0,76, ) : ( 0.9,0.8,0.7,0.6,0.48,0.24,0 ) : ( : t F t S t

73 data data; nput t cards; LIFETEST のプログラム ; proc lfetest plots=s(atrsk=0 to 0 b ); tme t*c();run; 73

74

75 * * * 時点 6の累積分布関数 τ 時点 9の累積分布関数 τ の重みで死亡 この式をτ について解くとKM 推定量

76 EXAMPLE: PRIMARY BILIARY CIRRHOSIS STUDY 原発性胆汁性肝硬変 (PBC) N=48 死亡 :6, 打ち切り :257 Tme: 生存時間 ( 年 ) Status: 打ち切り変数 ( : 死亡,0: 打ち切り ) Age: 試験登録時の年齢 Albumn: アルブミン濃度 (g/dl) Blrubn: ビリルビン (mg/dl) Edema: 浮腫の有無 Protme: プロトロンビン時間 ( 秒 )

77 PHREG による解析 proc phreg data=pbc model Tme*Status(0)=logBlrubn logprotme logalbumn Age Edema; 最尤推定値の分析 パラメータ 自由度 パラメータ推定値 標準誤差 カイ 2 乗値 Pr > Ch Sq ハザード比 logblrubn < logprotme logalbumn < Age < Edema

78 PHREG による生存時間分布の予測 プロットする参照共変量 ( 平均値 ) logblrubn logprotme logalbumn Age Edema

79 LIFEREG による解析 ( ワイブル ) proc lfereg data=pbc; model Tme*Status(0)=logBlrubn logprotme logalbumn Age Edema; Analss of Mamum Lkelhood Parameter Estmates Parameter DF Estmate Standard Error Ch- Square Pr > ChSq Intercept <.000 logblrubn <.000 logprotme logalbumn <.000 Age <.000 Edema Scale Webull Shape

80 QUANTLIFE のプログラム ods graphcs on; proc quantlfe data=pbc log method=na plot=(quantplot survval) seed=268; model Tme*Status(0)=logBlrubn logprotme logalbumn Age Edema /quantle=( ); run;

81 分位点に対するパラメータ推定値

82 QUANTLIFE のパラメータ推定値 LIFEREG の推定値 有意差なし.67 生存時間を有意に低下 82

83 QUANTLIFE のパラメータ推定値 分位点によらず, 回帰係数はほぼ一定 加速モデルの妥当性が示唆 83

84 QUANTLIFEによる生存時間分布の予測 plot=(survval) オプション プロットする参照共変量 ( 平均値 ) logblrubn logprotme logalbumn Age Edema

85 QUANTLIFE の利点 任意の % 点に対して, 特定の生存時間分布を仮定しないで共変量の影響を評価, 生存時間分布を予測 任意の % 点に対する係数の一様性を評価することで加速モデルの妥当性を評価 ( 加速モデルでは任意の % 点に対して 係数が一定 ) 85

86 参考文献 Koenker, R. and Bassett, G. W. (978), Regresson Quantles, Econometrca, 46, SAS Global Forum 203 Statstcs and Data Analss Koenker, R. and Gelng, O. (200), Reapprasng Medfl Longevt: A Quantle Regresson Survval Analss, Journal of the Amercan Statstcal Assocaton, 96, Peng, L. and Huang, Y. (2008), Survval Analss wth Quantle Regresson Models, Journal of the Amercan Statstcal Assocaton, 03, Portno, S. (2003), Censored Regresson Quantles, Journal of the Amercan Statstcal Assocaton, 98, Ln, G. and Rodrguez,N.R. (203), Usng the QUANTLIFE Procedure for Survval Analss SAS Global Forum 203,Paper 42 86