log-rank検定と層別log-rank検定の計算過程をSASで考える

2024年10月18日 log-rank検定と層別log-rank検定の計算過程を SASで考える臨床開発事業本部データサイエンスセンター統計解析２部飯田絢子

目次 1. 用語,添え字について 2. log-rank検定と層別log-rank検定の概説 3. log-rank(通常/層別)検定とCMHの関係 4. 通常のlog-rank検定の計算イメージ 5. 層別log-rank検定の計算イメージ 6. 通常のlog-rank検定の再現 7. 層別log-rank検定の再現 8. まとめ Copyright©EPS All rights reserved. 2

3.

用語，添え字について • 用語 ➢ リスク人数：「number at risk」や「リスク集合」とも呼ばれるが，本発表では「リスク人数」とする ➢ 通常のlog-rank：本発表において，層別log-rank検定に対して，層別でない場合のlog-rankを指す. ➢ 層別(予後)因子：アウトカムに影響のある因子.「交絡因子」と記載されることもある • 使う添え字 ➢ i：層番号 ➢ j：時点番号 Copyright©EPS All rights reserved. 3

4.

5.

log-rank検定とは • 帰無仮説H0：S1(t)=S2(t). 群間で生存時間関数が同じである • 対立仮説H1：S1(t)≠S2(t). 群間で生存時間関数が異なる ෠S t = ςtj<t(1 − dj ) nj dj :時点tjでのイベント発生人数 nj：number at risk. 全体の例数 - 時点tjより前に死亡または打ち切りを起こした人数 ※添え字ｊ：よく「i」が用いられますが，以降で層番号「i」を用いるため，時点を「j」としています • S1(t)=S2(t)か，どのように調べる？死亡，打ち切り：カテゴリカルデータ →分割表を用いて調べる Outcome Drug Outcome 合計反応あり反応なし a ｂリスク人数イベント発生イベント未発生 m1 Drug O1j N1j – O1j N1j Placebo O2j N2j – O2j N2j (Nj-N1j) Oj Nj-Oj Nj Factor Factor Placebo Total ｃｄ m2 n1 n2 N 生存時間解析以外で見る分割表 Copyright©EPS All rights reserved. Total 生存時間解析での分割表(j時点) 以降，log-rankの分割表と呼ぶ． 5

6.

7.

log-rank(通常/層別)検定とCMHの関係通常のlog-rankおよび層別log-rank CMH(Cochran-Mantel-Haenszel検定) i 番目の層におけるj 時点目※通常のlog-rankは層「i」なし石本りさ. Cochran-Mantel-Haenszel検定の概要とSASでの実装方法. SASユーザー総会2024. https://sas-user2024.ywstat.jp/download.html?n=10&key=paper，Accessed Oct 10, 2024 ※j：イベントが発生した時点のみが対象. イベント発生最終時点はt とする(j=1~t) 層はk 個あるとする(i=1~k) Outcome Drug Factor Placebo リスク人数イベント発生イベント未発生 O1ij N1ij – O1ij (event_d) O2ij (event_p) N2ij – O2ij Oij Nij-Oij Total N1ij (ret_lag_risk_d) N2ij (Nij-N1ij) (ret_lag_risk_p) Nij 上記表のDrug ×イベント発生セル位置(赤枠で囲った位置) において， E1ij=N1ij*Oij/Nij V1ij=Oij*N1ij*(Nij-Oij)*N2ij / ( (Nij**2)(Nij-1) ) 2 σk σtj=1 O1ij−E1ij i=1 χ2L = Z12= σk σt i=1 j=1 V1ij 2 t σk σ O1ij−N1ij∗Oij/Nij i=1 j=1 = Oij∗N1ij∗(Nij−Oij)∗N2ij σk σt i=1 j=1 Nij2 Nij−1 文字が違っていて見づらいが，分割表に対応する位置を見ると，同じ式を使っていることが分かる ※層別log-rankは時点jを考慮するので，式としてはσtj=1 が追加される点のみCMHの式と異なる ()：プログラムでの変数名 ※層別log-rankのプログラム上ではiは出てこない (sex=F, Mで層別．i=1は”F”, i=2は”M”のように置き換えて考える) Copyright©EPS All rights reserved. 7

https://sas-user2024.ywstat.jp/download.html?n=10&key=paper

8.

通常のlog-rankの計算イメージ① データから，時点ごとの各種発生数やリスク人数の情報を作成するこのような，両群イベント0のデータ行については分割表にまとめないイベント発生時点ごとに分割表に落とし込む Outcome Drug Factor Placebo Total リスク人数 Outcome イベント発生イベント未発生リスク人数 j=t O1ij N1ij Drug (event_d) イベント未発生 (ret_lag_risk_d) イベント発生 … Outcome リスク人数 Factor O1ijO2ij N1ijN2ij Drug (Nij-N1ij) Placebo (event_d) (ret_lag_risk_d) (event_p) イベント未発生イベント発生 Factor Outcome N2ij(ret_lag_risk_p) j=4(Day 179) リスク人数 O2ij O1ij N1j Total Oij Nij-Oij Nij (Nij-N1ij) Placebo Drug (event_p) (event_d) (ret_lag_risk_d) (ret_lag_risk_p) イベント発生イベント未発生 Factor j=3(Day 171) Total Nij-Oij Nij N2j O2ijOij O1j N1j (Nj-N1j) Placebo N1j(event_p) – O1j (event_d) (ret_lag_risk_d) (ret_lag_risk_p) j=2(Day 157) N2j Total O2j Oij Nij-Oij Nj (Nj-N1j) N2j – O2j (event_p) (ret_lag_risk_p) j=1(Day 156) Oj Nj-Oj Nj 時点ごとの結果の統合は，どのように計算しているのか？ Copyright©EPS All rights reserved. 計算過程を次ページで紹介する 8

9.

通常のlog-rankの計算イメージ② イベント発生時点ごとに分割表に落とし込む Outcome リスク人数 j=t (O1t-E1t)の計算結果 + … + V1tの計算結果 + … + ＝＝ Outcome イベント発生イベント未発生 … リスク人数 O1ij N1ij Drug (event_d) イベント未発生 (ret_lag_risk_d) イベント発生 Outcome リスク人数 Factor N2ij O1ijO2ij N1ij j=4(Day 179) (O14-E14)の計算結果 V14の計算結果 Drug (Nij-N1ij) Placebo (event_d) (ret_lag_risk_d) + + (event_p)イベント未発生イベント発生 Factor N2ij(ret_lag_risk_p) j=3(Day 171) (O13-E13)の計算結果 V13の計算結果 O2ijO1j N1j Total Oij Nij-Oij Nij (Nij-N1ij) Placebo Drug N1j – O1j (event_p) (event_d) (ret_lag_risk_d) + + (ret_lag_risk_p) Factor N2j Total Oij O2j Nij-Oij Nij j=2(Day 157) (O12-E12)の計算結果 V12の計算結果 (Nj-N1j) Placebo N2j – O2j (event_p) + + (ret_lag_risk_p) Total Oj Nj-Oj Nj j=1(Day 156) (O11-E11)の計算結果 V11の計算結果 E1ij=N1ij*Oij/Nij V1ij=Oij*N1ij*(Nij-Oij)*N2ij / ((Nij**2)(Nij-1)) χL2＝Z12＝ t ෍ O1j−E1j j=1 2 ÷ t ෍ V1j j=1 ・イベント発生時点ごとに(O1j-E1j)が計算される 2 t →t時点まで足して二乗 σj=1 O1j−E1j ・イベント発生時点ごとにV1jが計算される →時点tまで足す 2 ・χL2＝z12＝ σtj=1 O1j−E1j ÷σtj=1 V1j ・prob = 1 - probchi(z1, 1)． probchi(確率変数の値z1,自由度)：指定した自由度のカイ二乗分布において確率変数が統計量z1以下となる範囲の確率 Copyright©EPS All rights reserved. 9

10.

層別log-rankの計算イメージ① データを層別し，時点ごとの各種発生数やリスク人数の情報を作成する男性(i=1) 女性(i=2) d群(Drug X群) 時点男性層イベント (i=1) 発生数 ↓ ↓ ↓ 打ち切りリスク人数発生数 ↓ ↓ Outcome リスク人数 p群(Placebo群) イベント打ち切りリスク人数発生数発生数 ↓ ↓ ↓ イベン Outcome イベント未発ト発生リスク生イベン人数イベン O1ij ト未発 N1ij Outcome (event_ (ret_lag_ Drugト発生リスク生 Fa 人数 O1ijd)イベン N1ijrisk_d) イベン j=t ト未発 cto N2ij (event_ (ret_lag_ Drug ト発生 O2ij … 生 (Nij-N1ij) d) (event_ risk_d) Fa r Placeb o O1ij N1ijN2ij(ret_lag_ cto O2ijp) N1ij – (ret_lag_ (event_ risk_p) Drug (Nij-N1ij) r Placeb j=4(Day 253) O1ij (event_ d) risk_d) Fa Total Oij Nij-Oij(ret_lag_ Nij o p) cto N2ijrisk_p) j=3(Day 249) O2ij (Nij-N1ij) r Total Placeb N2ij – Oij Nij-Oij Nij (event_ o O2ij (ret_lag_ j=2(Day 242) p) risk_p) Total Oij Nij-Oij Nij j=1(Day 237) Copyright©EPS All rights reserved. d群(Drug X群) 時点女性層イベント打ち切りリスク人数 (i=2) 発生数発生数 ↓ ↓ ↓ ↓ ↓ p群(Placebo群) イベント打ち切りリスク人数発生数発生数 ↓ ↓ ↓ Outcome リスクイベン人数 Outcome イベント未発ト発生リスク生イベン人数イベン O1ij ト未発 N1ij Outcome ト発生 (event_ (ret_lag_ Drug リスク生 d)イベン risk_d) 人数 O1ij N1ij Fac イベン j=t ト未発 N2ij (event_ (ret_lag_ Drugト発生 tor O2ij … 生 (Nij-N1ij) d) (event_ risk_d) Fa Placebo (ret_lag_ O1ij N1ijN2ij cto p) N1ij – O2ij (event_ (ret_lag_ risk_p) Drug (Nij-N1ij) r Placeb j=4(Day 179) O1ij (event_ d) risk_d) Fa Total Oij Nij-Oij (ret_lag_ Nij o p) cto N2ijrisk_p) j=3(Day 171) O2ij (Nij-N1ij) r Total Placeb N2ij – Oij Nij-Oij Nij (event_ o O2ij (ret_lag_ j=2(Day 157) p) risk_p) Total Oij Nij-Oij Nij j=1(Day 156) 10

11.

層別log-rankの計算イメージ② 女性(i=2) 男性(i=1) Outcome リスク人 j=t (O1t-E1t)の計算結果イベンイベント数ト発生未発生 Outcome + リスク人 … O1ij N1ij … イベンイベント数 (event_ (ret_lag_r Outcome Drug ト発生未発生リスク人 + isk_d) イベン O1ijd) イベント N1ij数 Fac (O14-E14)の計算結果 j=4 未発生 N2ij (event_ (ret_lag_r Drug ト発生 tor O2ij O1ij N1ij (Nij-N1ij) d) (event_ N1ij – isk_d) FacDrugPlacebo (event_ (ret_lag_r N2ij(ret_lag_r j=3 (O13-E13)の計算結果 p) O1ij tor isk_p) d)O2ij isk_d) (Nij-N1ij) Fac Placebo (event_ + Total Oij Nij-OijN2ij Nij (ret_lag_r tor O2ij p) N2ij – (Nij-N1ij) (O12-E12)の計算結果 j=2 isk_p) Placebo (event_ O2ij (ret_lag_r Total Nij-Oij Nij p)Oij isk_p) (O11-E11)の計算結果 j=1 Total Oij Nij-Oij Nij V1tの計算結果 j=t + … … ＝ t ෍ O1j−E1j j=1 V14の計算結果 χL2＝ t ෍ O1j−E1j j=1 女性＋ j=4 (O14-E14)の計算結果 + V1tの計算結果 … j=3 + (O13-E13)の計算結果 V13の計算結果 + + V12の計算結果 j=2 + (O12-E12)の計算結果 + + V11の計算結果 j=1 ＝ t ෍ V1j j=1 t ෍ O1j−E1j j=1 ÷ + 男性 V12の計算結果 + (O11-E11)の計算結果＝ t ෍ V1j j=1 + + V14の計算結果 + V13の計算結果 2 t ෍ O1j−E1j j=1 男性 + + … + + + (O1t-E1t)の計算結果 V11の計算結果＝ t ෍ V1j j=1 t ෍ V1j j=1 女性層別データに関して，・各時点×層別に(O1j-E1j)が計算される→層別にt時点まで足す→男性の σtj=1 O1j−E1j と女性の σtj=1 O1j−E1j を足してから二乗・各時点×層別にV1jが計算される→層別に時点tまで足す→男性の σtj=1 Vj と女性のσtj=1 Vj を足す t σ ・χL2＝Z12＝ σk i=1 j=1 O1ij−E1ij ・prob = 1 - probchi(z1, 1) Copyright©EPS All rights reserved. 2 t σ ÷σk i=1 j=1 Vij 11

12.

13.

使用データセット SAS Help Center: Getting Started: LIFETEST Procedure の例を用いる data exposed; input Days Status Treatment Sex $ @@; format Treatment Rx.; datalines; 179 1 1 F 378 0 1 M 256 1 1 F 355 1 1 M 262 1 1 M 319 1 1 M 256 1 1 F 256 1 1 M 255 1 1 M 171 1 1 F 224 0 1 F 325 1 1 M 225 1 1 F 325 1 1 M 287 1 1 M 217 1 1 F 319 1 1 M 255 1 1 F 264 1 1 M 256 1 1 F 237 0 0 F 291 1 0 M 156 1 0 F 323 1 0 M 270 1 0 M 253 1 0 M 257 1 0 M 206 1 0 F 242 1 0 M 206 1 0 F 157 1 0 F 237 1 0 M 249 1 0 M 211 1 0 F 180 1 0 F 229 1 0 F 226 1 0 F 234 1 0 F 268 0 0 M 209 1 0 F ; Copyright©EPS All rights reserved. Days：時間(日数)を表す変数 Status：0=打ち切り１=イベント Treatment：Drug X, Placebo Sex：F, M(層別因子) Drug X群の割付人数：20人 Placebo群の割付人数：20人 Sex＝Fにおける Drug X群の割付人数：9人 Placebo群の割付人数：11人 sex=Mにおける Drug X群の割付人数：11人 Placebo群の割付人数：9人 13

https://go.documentation.sas.com/doc/ja/pgmsascdc/9.4_3.4/statug/statug_lifetest_gettingstarted.htm

14.

15.

通常のlog-rank検定をプロシジャで行うと proc lifetest data=Exposed time Days*Status(0); strata Treatment/ test=LOGRANK; run; 構文： proc lifetest data=[入力データセット]; time [時間を表す変数] * [イベント/打ち切りを表す変数(打ち切り側の水準)]; strata [比較したい群変数] / test=LOGRANK; 例）変数名Status run; 0=打ち切り１=イベントの場合，「Status(0)」この計算過程をSASで考えていく Copyright©EPS All rights reserved. 15

16.

通常のlog-rankの分割表を完成させるために目指すデータセットの形プログラム1_通常のlogrank検定_iida.sasのデータセットWk4に該当 d群(Drug X群) 時点 ↓ イベント発生数打ち切り発生数 ↓ ↓ Copyright©EPS All rights reserved. リスク人数 ↓ p群(Placebo群) イベント発生数打ち切り発生数 ↓ ↓ リスク人数 ↓ Censr_dとCensr_pはlog-rankの検定の計算に直接は使用しないが， ret_lag_risk_d, ret_lag_risk_pを算出するために必要 16

17.

通常のlog-rank計算のためのデータセット作成① 群×Daysごとのイベント発生数と打ち切り発生数をカウント /*群×Daysごとのイベント,打ち切りの件数をカウント*/ proc sort data=Exposed;by Treatment Days;run; proc freq data=Exposed(where=(status^=.)) noprint; tables Status/out=wk1; by Treatment Days; run; /*転置*/ proc sort data=wk1;by Days;run; proc transpose data=wk1 out=wk1_2 delimiter = _; id Treatment Status; var COUNT; by Days; run; /*変数名変換*/ data wk1_3; set wk1_2; rename Drug_X_0=Censr_d 変数名の接尾語，接頭語について. Drug_X_1=Event_d XXXXX_d：Drug_X群情報, Placebo_0=Censr_p XXXXX_p：Placebo群情報, Placebo_1=Event_p Event_X：イベント発生件数, ; Censr_X：打ち切り発生件数 run; ※「X」は何かしら文字が入る箇所 Copyright©EPS All rights reserved. 17

18.

通常のlog-rank計算のためのデータセット作成② atrisk(リスク人数)の作成 /*転置してDrug_Xのatrisk及びPlaceboのatrisk列を作成*/ proc sort data=o_survival out=NDK_atrisk nodupkey; by Time atrisk; run; proc transpose data=NDK_atrisk out=wk1_atrisk; id Stratum;/*Treatment*/ by Time ;/*Days*/ var atrisk; run; ods output SURVIVALPLOT = o_survival; proc lifetest data=Exposed; time Days*Status(0); strata Treatment/ test=LOGRANK; run; SURVIVALPLOTによるatriskの計算再現については, プログラム0_リスク人数の計算_iida.sasを参照 /*Time=0の行のatriskを埋める(両群とも割付20人のため，20を入れる)*/ data wk2_atrisk; set wk1_atrisk; if Time=0 and Drug_X=. then atrisk_d=20; else atrisk_d=Drug_X; if Time=0 and Placebo=. then atrisk_p=20; else atrisk_p=Placebo; rename Time =Days; run; Copyright©EPS All rights reserved. 18

19.

通常のlog-rank計算のためのデータセット作成③ wk2_atrisk, wk1_3データセットのマージ data wk3; merge wk2_atrisk wk1_3; by Days; run; d群(Drug X群) イベント発生数 ↓ 打ち切り発生数 ↓ p群(Placebo群) リスク人数 ↓ イベント発生数 ↓ 打ち切り発生数 ↓ リスク人数 ↓ 次ページからは，欠測になっている箇所を埋めていくプログラム．例えば， Placebo群でイベントか打ち切りが発生した(例ではイベントのみ発生)が， Drug_X群では何も発生していないDay156~157については，SURVIVALPLOT ( o_survivalデータセット)ではリスク人数の計算がされていないため，手動でリスク人数を算出していく．Placebo群も同様に欠測しているリスク人数箇所を計算． Copyright©EPS All rights reserved. 19

20.

通常のlog-rank計算のためのデータセット作成④ リスク人数の計算 data Wk4; set Wk3; /*各群のイベントと打ち切りの変数について，欠測値は0件に直す*/ if event_d=. then event_d=0; if Censr_d=. then Censr_d=0; if event_p=. then event_p=0; if Censr_p=. then Censr_p=0; /*＝＝Drug_X群のリスク人数(ret_lag_risk_d)の計算＝＝*/ /*①各時点における，atrisk - イベント件数 - 打ち切り件数を計算*/ if atrisk_d^=. then risk_d=atrisk_d - Event_d - Censr_d; /*②前の行の①の値を引き継ぐ*/ lag_risk_d=lag(risk_d); /*③②の操作でもまだlag_risk_dが欠測の場合(欠測レコードが2連続以上)は，欠測している時点より前の，値を持つ最終行のlag_risk_dをretainする*/ retain ret_lag_risk_d; if lag_risk_d^=. then ret_lag_risk_d=lag_risk_d; /*④1行目データ(Day 0)のret_lag_risk_dは割付人数(20人)にする*/ if _N_=1 then ret_lag_risk_d=20; run; 同様に，Placebo群のリスク人数(ret_lag_risk_p)を計算詳細はプログラム1_通常のlogrank検定_iida.sasの Wk4参照 Copyright©EPS All rights reserved. 20

21.

リスク人数の求め方とリスク人数の欠測が続く場合の埋め方 Day 156 のリスク人数は 20-0=20 Day 179のリスク人数は 19-1=18 (20-1-1 =18) ・リスク人数の計算：リスク人数＝割付人数－(リスク人数を求めたい行の前の行までのイベントと打ち切り累計数) →前の行にリスク人数があれば，リスク人数がある行の，「リスク人数 - 同行のイベント数 - 同行の打ち切り数」が次の行のリスク人数となる・リスク人数の欠測レコードが連続する場合，前の欠測行を埋めてからリスク人数を順に計算していく →プログラムが少し複雑？ Copyright©EPS All rights reserved. 21

22.

通常のlog-rank計算のためのデータセット作成⑤ log-rank計算用にrename data Wk5; set Wk4; rename ret_lag_risk_d=N1j ret_lag_risk_p=N2j event_d=O1j event_p=O2j ; run; 時点 j の分割表 Outcome Factor Total Copyright©EPS All rights reserved. リスク人数 (死亡する可能性のある人数) イベント発生イベント未発生 O1j N1j N1j – O1j Drug (event_d) (ret_lag_risk_d) N2j O2j N2j – O2j Placebo (Nj-N1j) (event_p) (ret_lag_risk_p) Oj Nj-Oj Nj 22

23.

通常のlog-rank検定の値(統計量,p値) の算出① data Wk6; set Wk5 END=EOF; /*少なくとも一方の群でイベントが発生したデータに絞る*/ where Days^=0 and ^(O1j=0 and O2j=0); Nj=N1j+N2j; Oj=O1j+O2j; E1j=N1j*Oj/Nj; 時点 j の分割表リスク人数 (死亡する可能性のある人数) Outcome Factor イベント発生イベント未発生 O1j N1j N1j – O1j Drug (event_d) (ret_lag_risk_d) N2j O2j N2j – O2j Placebo (Nj-N1j) (event_p) (ret_lag_risk_p) Total Oj Nj-Oj Nj E1j=N1j*Oj/Nj /*＝＝z12の分母＝＝*/ V1j /*各時点*/ Oj∗N1j∗(Nj−Oj)∗(Nj−N1j) V1j=divide (Oj*N1j*(Nj-Oj)*(Nj-N1j),(Nj**2)*(Nj-1)); = Nj2 Nj−1 時点和 /*V1jの時点和*/ V1+V1j; /*z12の分母*/ if EOF=1 then do; denom=V1; end; V1=σtj=1 V1j χL2＝Z12＝ t ෍ O1j−E1j j=1 2 ÷ t ෍ V1j j=1 次ページに続く Copyright©EPS All rights reserved. 23

24.

通常のlog-rank検定の値(統計量,p値) の算出② /*data Wk6の続き*/ 時点 j の分割表 /*＝＝z12の分子＝＝*/ /*各時間における，期待値と実測値の差*/ OE1j=(O1j-E1j); /*全時間における，期待値と実測値の差(全時点和)*/ OE1+OE1j; リスク人数 (死亡する可能性のある人数) Outcome Factor イベント発生イベント未発生 O1j N1j N1j – O1j Drug (event_d) (ret_lag_risk_d) N2j O2j N2j – O2j Placebo (Nj-N1j) (event_p) (ret_lag_risk_p) Total Oj Nj-Oj Nj E1j=N1j*Oj/Nj OE1j=(O1j-E1j) 時点和 OE1=σtj=1 OE1j /*z12の分子*/ if EOF=1 then do; Num=OE1**2; end; /*z12の計算 ※プログラム上は変数名「z12」が書けない為「Z1」と表記*/ if EOF=1 then do; z1=Num/denom; end; Copyright©EPS All rights reserved. χL2＝Z12＝ t ෍ O1j−E1j j=1 Num=OE1**2 2 t ÷ ෍ V1j j=1 denom=V 1 (前ページ) z12=Num/denom 24

25.

通常のlog-rank検定の値(統計量,p値) の算出③ /*data Wk6の続き*/ /*p値計算*/ if EOF=1 then do; prob = 1 - probchi(z1, 1); end; run; ・2×2分割表の為，カイ二乗分布の自由度は1 ・probchi(z1,自由度)：指定した自由度のカイ二乗分布において，横軸(確率変数)が統計量z1以下となる範囲の確率青い部分：probchi(z1, 1) 統計量(変数z1)の値白い部分：p値一致．算出した「z1」=SAS出力結果「カイ二乗値」算出した「prob」=SAS出力結果「Pr>Chi-Square」 Copyright©EPS All rights reserved. 25

26.

27.

層別log-rank検定をプロシジャで行うと proc lifetest data=Exposed ; time Days*Status(0); strata Sex/ group=Treatment test=LOGRANK ; run; 構文： proc lifetest data=[入力データセット]; time [時間を表す変数] * [イベント/打ち切りを表す変数(打ち切り側の水準)]; strata [層別変数] / group=[比較したい群変数] test=LOGRANK; 例）変数名Status run; 0=打ち切り１=イベントの場合，「Status(0)」この計算過程をSASで考えていく Copyright©EPS All rights reserved. 27

28.

層別log-rankの分割表を完成させるために目指すデータセットの形プログラム2_層別logrank検定_iida.sasのデータセットWk5_0に該当 d群(Drug X群) 性別時点 ↓ ↓ イベント発生数 ↓ 打ち切り発生数 ↓ p群(Placebo群) リスク人数 ↓ イベント発生数 ↓ 打ち切り発生数 ↓ リスク人数 ↓ 対応プログラムスクショ入れる Copyright©EPS All rights reserved. Censr_dとCensr_pはlog-rankの検定の計算に直接は使用しないが， ret_lag_risk_d, ret_lag_risk_pを算出するために必要 28

29.

層別log-rank計算のためのデータセット作成①② • 通常のlog-rank検定のプログラムとほぼ同じの為，省略. プログラム2_層別logrank検定_iida.sasにおける「data wk2_atrisk;」まで省略 • 通常のlog-rank検定のプログラムと異なる箇所：層別因子であるSEXを入れる例） /*性別×群×Daysごとのイベント発生数と打ち切り発生数をカウント*/ proc sort data=Exposed;by Sex Days;run; proc freq data=Exposed(where=(status^=.)) noprint; tables Treatment*Status/out=wk1; by Sex Days; run; Copyright©EPS All rights reserved. 29

30.

層別log-rank計算のためのデータセット作成③ wk2_atrisk, wk1_3データセットのマージ data wk3; format Sex Days event_d Censr_d atrisk_d event_p Censr_p atrisk_p; keep Sex Days event_d Censr_d atrisk_d event_p Censr_p atrisk_p; merge wk2_atrisk wk1_3; by Sex Days; run; 次ページからは，性別ごとにデータセットを分割した上で，欠測になっている箇所を埋めていくプログラム．女性のPlacebo群でイベントか打ち切りが発生したが，女性のDrug_X群では何も発生していないDay226~237については，SURVIVALPLOT ( o_survival データセット)ではリスク人数の計算がされていないため，手動でリスク人数を算出していく．Placebo群も同様に欠測しているリスク人数箇所を計算．男性についても同様に計算 Copyright©EPS All rights reserved. 30

31.

層別log-rank計算のためのデータセット作成④ リスク人数の計算「log-rank計算のためのデータセット作成④」スライドと同様の為，以降の記載は省略する %macro Wk4(Sex); data Wk4_&Sex.; set Wk3; /*層(性別)別データセットにしたうえで以下計算していく*/ where Sex="&Sex."; /*各群のイベントと打ち切りの変数について，欠測値は0件に直す*/ /*＝＝性別ごとのDrug_X群のリスク人数(ret_lag_risk_d)の計算＝＝*/ /*①各行のリスク人数(ret_lag_risk_d)を計算*/ /*②前の行のリスク人数を引き継ぐ*/ /*③lag_risk_dが欠測の場合は欠測前の行のlag_risk_dをretainする*/ /*④1行目データ(Days=0)のret_lag_risk_dは性別＊群の割付人数にする*/ if Sex="F" and Days=0 then do; ret_lag_risk_d=9; end; if Sex="M" and Days=0 then do; ret_lag_risk_d=11; end; 同様に，Placebo群のリスク人数(ret_lag_risk_p)を計算 %mend Wk4; %Wk4(Sex=F); %Wk4(Sex=M); Copyright©EPS All rights reserved. 詳細はプログラム2_層別logrank検定_iida.sasの%macro Wk4参照 31

32.

層別log-rank計算のためのデータセット作成⑤ log-rank計算用にrename data Wk5_0; set Wk4_F Wk4_M; ; run; data Wk5_1; set Wk5_0; rename ret_lag_risk_d=N1j ret_lag_risk_p=N2j event_d=O1j event_p=O2j ; run; i 番目の層におけるj 時点目 Outcome Factor Total リスク人数イベント発生イベント未発生 O1ij N1ij N1ij – O1ij Drug (event_d) (ret_lag_risk_d) N2ij O2ij N2ij – O2ij Placebo (Nij-N1ij) (event_p) (ret_lag_risk_p) Oij Nij-Oj Nij 以降スライドにおける，プログラムと分割表の添え字記号の差について data Wk5,%macro Wk4(Sex)において，層ごとでデータセット内の変数名を変えていないため，プログラムにおける変数名に「i」はない．プログラム変数名は分割表の「i」を除いた名前と対応している．例）プログラムの変数N1jは分割表のN1ijに対応 Copyright©EPS All rights reserved. 32

33.

層別log-rank検定の値(統計量,p値) の算出① %macro Wk6:層(Sex)ごとに統計量の要素を計算するマクロ i 番目の層におけるj 時点目 Outcome %macro Wk6(Sex); イベント発生イベント未発生 data Wk6_&Sex.; O1ij N1ij – O1ij Drug (event_d) set Wk5_1 END=EOF; Factor O2ij N2ij – O2ij Placebo /*層別データセットに関し， (event_p) Total Oij Nij-Oj 少なくとも一方の群でイベントが発生したデータに絞る*/ where Sex="&Sex." and Days^=0 and ^(O1j=0 and O2j=0); Nj=N1j+N2j; Oj=O1j+O2j; リスク人数 N1ij (ret_lag_risk_d) N2ij (Nij-N1ij) (ret_lag_risk_p) Nij /*＝＝z12の分母の要素(統合前の各層のz12の分母)＝＝*/ V1j=divide (Oj*N1j*(Nj-Oj)*(Nj-N1j),(Nj**2)*(Nj-1));/*各時点*/ V1+V1j;/*時点和*/ if EOF=1 then do; denom=V1; end; 次ページに続く Copyright©EPS All rights reserved. 33

34.

層別log-rank検定の値(統計量,p値) の算出② /*%macro Wk6の続き*/ /*＝＝z12の分子の要素(統合前の各層のz12の分子)＝＝*/ E1j=N1j*Oj/Nj; /*各時間における，期待値と実測値の差*/ OE1j=(O1j-E1j); /*全時間における，期待値と実測値の差(全時点和)*/ OE1+OE1j; if EOF=1 then do; Num=OE1;/*この時点では二乗しない*/ end; run; /*後のmergeのために，変数名をrenameし，分子と分母(最終行)のみ抽出*/ data Wk6_&Sex._2; set Wk6_&Sex. END=EOF; keep Days Num Denom BEF_sum; rename Num=Num_&Sex. Denom=Denom_&Sex. ; if EOF=1 then do;BEF_sum=1; output;end; run; %mend Wk6; Num_M Num_F %Wk6(Sex=F);/*女性側のデータセット作成*/ %Wk6(Sex=M);/*男性側のデータセット作成*/ Copyright©EPS All rights reserved. Denom_M Denom_F 34

35.

層別log-rank検定の値(統計量,p値) の算出③ /*層別で計算したz12の分子,分母要素についてmerge*/ data Wk7; merge Wk6_F_2 Wk6_M_2; by BEF_sum; run; /*層別で計算したz12の分子,分母要素について計算結果の統合*/ data Wk8; set Wk7; total_Num=Num_F + Num_M; total_denom=Denom_F + Denom_M; total_Num /*z12の計算 ※プログラム上は変数名「z12」が書けない為「Z1」と表記*/ z1=total_Num**2/total_denom; /*p値の計算*/ prob = 1 - probchi(z1, 1); run; total_denom 一致．算出した「z1」=SAS出力結果「カイ二乗値」算出した「prob」=SAS出力結果「Pr>Chi-Square」 Copyright©EPS All rights reserved. 35

36.

まとめ • log-rank検定は，生存時間関数が群間で異なるか検定する手法. • 生存時間関数が群間で異なるかは，以下の分割表を用いて統計量を計算 →カイ二乗分布でp値を計算 2 t σk i=1 σj=1 O1ij−E1ij 2 2 χL = Z1 = σk σt i=1 j=1 V1ij 2 t σk σ O1ij−N1ij∗Oij/Nij i=1 j=1 = Oij∗N1ij∗(Nij−Oij)∗N2ij σk σt i=1 j=1 Nij2 Nij−1 E1ij=N1ij*Oij/Nij V1ij=Oij*N1ij*(Nij-Oij)*N2ij / ((Nij**2)(Nij-1) i 番目の層におけるj 時点目 Outcome Drug Factor Placebo Total リスク人数イベント発生 O1ij (event_d) イベント未発生 O2ij (event_p) N2ij – O2ij Oij Nij-Oj N1ij – O1ij N1ij (ret_lag_risk_d) N2ij (Nij-N1ij) (ret_lag_risk_p) Nij • SASプロシジャではproc lifetest で数行で検定結果が出力できるが，上記の計算式を用いたSASプログラムで結果を再現した • 計算のポイント：時点和や層別の結果をいつ統合するか • 統計量Z12の分子と分母は別々に計算して時点和を求めてから計算する. • 層別の場合は，各層の分子要素の時点和結果を足し，各層の分母要素の時点和結果を足したうえで統計量Z12を計算 Copyright©EPS All rights reserved. 36

37.

参考文献 [1]SAS Help Center: Getting Started: LIFETEST Procedure https://go.documentation.sas.com/doc /ja/pgmsascdc /9.4_3.4/statug /statug_lifetest_gettingsta rted.htm，Accessed Oct 10, 2024 [2]中澤秀夫．―基礎科学から医学・医療を見る―生存時間データ解析と比例ハザードモデル. 日医大医会誌 2015; 11（1） https://www.jstage.jst.go.jp/article/manms /11/1/11_29/_pdf，Accessed Oct 10, 2024 [3]リハビリテーション統計学_ログランク検定 https://y2pt.com/13847/,Accessed Oct 10, 2024 [4]折村奈美. データステップを用いた統計検定の再現. SASユーザー総会2024. (2024) https://sas-user2024.ywstat.jp /download.html?n=63&key=paper，Accessed Oct 10, 2024 [5]石本りさ. Cochran-Mantel-Haenszel検定の概要とSASでの実装方法. SASユーザー総会2024. https://sas-user2024.ywstat.jp /download.html?n=10&key=paper，Accessed Oct 10, 2024 [6]層別解析とサブグループ解析の違いは？統合方法や結果の見方もわかりやすく！ | いちばんやさしい、医療統計 https://best-biostatistics.com/review/sub.html,Accessed Oct 10, 2024 [7]サブグループ解析の解釈｜医療統計｜ビーリンサイト https://www.blincyto.jp/documents/medical_statistics09,Accessed Oct 10, 2024 Copyright©EPS All rights reserved. 37

38.

log-rank検定と層別log-rank検定の計算過程をSASで考える

森岡裕[SASユーザー総会世話人]

関連スライド

MMRM入門

t検定の基礎(Studentのt検定とWelchのt検定)

SASのコールルーチンは57個あるねん10分で全部説明するねん

RTF出力を活用した解析報告書の作成

色々な確率分布

重回帰分析における多重共線性について

各ページのテキスト