そんなに頻繁に更新するつもりはありません。
小ネタより、長めの実用的なプログラムの紹介が多いかも。
第1回はオンコロジー領域でもよく使用する、Kaplan-Meier Curve（KM曲線）を作成してみようと思います。
KM曲線では、シリーズとしていくつか段階を経て自分なりにアレンジし、論文発表のFigureに使えそうなグラフを最終的に載せようと思っています。
早速ですが、①では治験総括報告書や論文でも使用できそうな基本的なモノクロのKM曲線を作成しようと思います。
私の場合は上記以外に、被験者数設計のSimulationの際の生存曲線のイメージを作成することにも使用したりします。
発生させた乱数が合っているかの確認にも使えるので（特に区分指数分布とか）。
今回のSGPLOTで重要なステートメントは、基本的ではありますがSTEP、SCATTER、XAXISTABLEあたりでしょうか。
まず最初に適当なデータから作成します。
手元にすでにTime to Eventのデータセットをお持ちの方はそちらを使用してみてください。
data work;
do k = 1 to 2;
do i=1 to 250;
y = rand('UNIFORM'); output;
end;
end;
run;
data adtte;
set work;
length TRTP $7;
/*---------------Treatment group---------------*/
if k = 1 then do ; TRTP = 'active' ; TRTPN = 1; end;
else if k = 2 then do ; TRTP = 'control'; TRTPN = 2; end;
/*---------------Hazard---------------*/
if TRTPN = 1 then hazard = -0.58*(log(0.5))/6 ;
else if TRTPN = 2 then hazard = -log(0.5)/6 ;
/*---------------Time to event---------------*/
AVAL = -log(y) / hazard;
/*---------------Censor---------------*/
CNSR = rand('BINOMIAL',0.8,1);
/*---------------Follow up time---------------*/
if AVAL >= 24 then do;
AVAL = 24; CNSR = 0;
end;
run;
とこんな感じに作成しました。24カ月で試験を打ち切りにしています。
作成した変数のうちPROC LIFETESTに使用する変数は、TRTPN（投与群の数値変数）、AVAL（Time to Event）、CNSR（打ち切り情報）になっています。ADaMに合わせてます。
このデータセットをPROC LIFETESTにかけてあげて、ods outputでAt risk情報も同時に出力します。
ods output Survivalplot = sgplot ;
proc lifetest data=adtte METHOD= KM atrisk maxtime = 24
plots = survival (outside atrisk = 0 to 24 by 3 atrisktick );
time AVAL*CNSR(0);
strata TRTPN / test = logrank;
run;
ここからデータ加工をしていきますが、データ加工せずにこのままSGPLOTでKM曲線を作成することもできます！
なので、飛ばしたい方はそのままSGPLOTのステップに飛んでいただいてもOKです。
ではなぜデータ加工をするのか...といいますと、ods outputの出力では、At riskで0になったら、それ以降のオブザベーションを作成してくれないんです（便利なオプションがあるならすいません）。
わかりにくいので説明します。
例えば今回の例ではAt risk数を24カ月まで3カ月ごとに出力するようになっています。
もしこの治験で片方の群が17カ月までに生存率が0になったとします。
そうすると、18カ月のAt riskは0でそれ以降のオブザベーションは全く発生しません。
KM曲線に、18、21、24のAt risk数を0で表示させたい場合、加工しないと21、24ではブランクで出力されてしまいます。
見栄えの問題で、私はブランクがないほうが良いと感じるので、加工していきます（今回のSample datasetではこの問題は起きません）。
data form ;
do StratumNum = 1 to 2 ;
do time = 0 to 24 by 3;
tAtRisk = time; output;
end;
end;
run;
data atrisk0;
merge form sgplot;
by StratumNum time ;
run;
data sgplot2;
set atrisk0;
if AtRisk = . then do;
AtRisk = 0;
tAtRisk = time;
end;
run;
proc sort data = sgplot2;
by StratumNum time ;
run;
これで、オブザベーションが発生しなかったAt riskを0で追加する加工を行っています。
SGPLOTをかける前に、StratumNumのFormatを作成しておきます。
proc format;
value treat 1 = "Active" 2 = "Control";
run;
それではSGPLOTにかけていきます。
ods graphics / width=1100px height=700px border=off
reset=index imagename= "KMcurve" imagefmt=tiff ;
proc sgplot data = sgplot2 noautolegend;
step x=Time y=Survival / group=StratumNum
lineattrs=(thickness=2 color=black ) name="a" ;
scatter x=Time y=Censored / group=StratumNum
markerattrs=(size=15 color=black symbol=plus) ;
xaxistable AtRisk / x=tAtRisk class=StratumNum
valueattrs=(size=13) labelattrs=(size=13)
title="No. at Risk" titleattrs=(size=13);
yaxis values=(0 to 1 by 0.1)
label="Probability of Overall Survival"
valueattrs=(size=13) labelattrs=(size=13);
xaxis values=(0 to 24 by 3) offsetmin=0.03 offsetmax=0.02
label="Time(Months)" valueattrs=(size=13)
labelattrs=(size=13);
keylegend "a" / location = inside position = topright
titleattrs=(size=13)valueattrs=(size=13)
acros=1 noborder ;
format StratumNum treat. ;
run;
完成図は以下になります。
画質悪くて申し訳ないですが、こんな出力になります。
細かいATTRSの設定はお好みで行ってください。
＊STEPで生存曲線をプロット
＊SCATTERで打ち切り記号をプロット
＊XAXISTABLEでAt risk数をプロット
＊YAXISとXAXISで軸の設定
＊KEYLEGENDで凡例の設定
＊FORMATでFormatの指定
となります。
次回以降、このプログラムをベースに追加を行っていきたいと思います。
STEPやSCATTERのATTRSオプションにおいて、色やシンボルや線種の指定は、1種類しかできません。
これは他のグラフを作成する場合でも同じです。
今回線種は指定していないので、デフォルトの2種類で出力されています。
色やシンボルも、ATTRSで指定しない場合、デフォルトの設定で2種類出力されます。
なので群ごとに色やシンボルを自分自身で好きに決めたい場合はどうすればいいか、
また、好きな色を指定する簡単な方法などを次回紹介してみようと思います。
と言っても2018年のユーザー総会で発表した内容ですが...
おわり
0 件のコメント:
コメントを投稿