如果画一个NEJM 期刊式的生存曲线

一般期刊的生存曲线是这样子的：


而 NEJM 的生存曲线是这样子的：

仔细观察，其实不难发现差别在于：

• 图中是否有No. at risk 的表格

• 图中是否内嵌了一个局部放大的K-M曲线

第一点，No. at risk 表格，利用 SAS 自带的生存分析过程步 proc lifetest中plot选项的atrisk子选项，可以轻松实现。但自带的过程步，效果确实不敢恭维。不过，好在我们可以用ods output抓取生存曲线的信息，然后用proc sgplot过程中的scatter及xaxistable语句完成。具体可见此前的推文:

• 如何绘制带No.at risk 表格的生存曲线？

第二点，内嵌一个局部放大的 K-M 曲线，利用 SAS 实现起来确实有一定的难度，原因是SAS 无法通过 sgplot过程步直接实现，因此需要借助GTL。

如何利用GTL绘制这种内嵌局部放大的K-M曲线呢? 这里，给大家提供一种思路。

• 利用自带过程proc lifetest及ods output抓取 K-M 曲线数据

• 利用 GTL 及抓取数据集绘制一个原尺寸的,带No. at risk的K-M 曲线

• 利用 GTL 及抓取数据集绘制一个局部放大的不带No. at risk的K-M 曲线

• 利用 GTL 定义一个绘图模板，将局部放大图内嵌入原始尺寸图

• 渲染绘图模板，获得嵌套的 K-M 图

由于 GTL 的语法比较繁杂，对大家来说，比较陌生，不过中间的两个GTL代码，我们可以利用proc sgplot的tmplout=选项，让 SAS 帮我们自动生成。

具体如何做呢?  请看下面的SAS代码及讲解。

*===1.抓取数据集;
ods output Survivalplot=SurvivalPlotData;
proc lifetest data=sashelp.BMT plots=survival(atrisk=0 to 2500 by 500);
time T * Status(0);
strata Group / test=logrank adjust=sidak;
run;

*===2.获得原始尺寸带no.at risk的GLT代码;
proc sgplot data=SurvivalPlotData noborder noautolegend tmplout="d:\p1.sas";
step x=time y=survival / group=stratum name='s';
scatter x=time y=censored / markerattrs=(symbol=plus) name='c';
scatter x=time y=censored / markerattrs=(symbol=plus) GROUP=stratum;
xaxistable atrisk / x=tatrisk class=stratum colorgroup=stratum;
run;

*===3.获得局部放大尺寸不带no.at risk的GLT代码;
proc sgplot data=SurvivalPlotData noborder tmplout="d:\p2.sas";
step x=time y=survival / group=stratum name='s';
scatter x=time y=censored / markerattrs=(symbol=plus) name='c';
scatter x=time y=censored / markerattrs=(symbol=plus) GROUP=stratum;
xaxis values=(0 to 500 by 100);
yaxis label=" ";
run;

*===4.将局部放大图内嵌与原始尺寸图;
proc template;
define statgraph cars_inset_bar_bar;
begingraph;
    layout lattice / columnweights=preferred rowweights=preferred columndatarange=union rowdatarange=union columns=1;
    layout overlay / walldisplay=(fill) xaxisopts=(Label="Days from randomization");
        StepPlot X=Time Y=Survival / primary=true Group=Stratum LegendLabel="Survival Probability" NAME="s";
        ScatterPlot X=Time Y=Censored / subpixel=off Markerattrs=( Symbol=PLUS) LegendLabel="Censored" NAME="c";
        ScatterPlot X=Time Y=Censored / subpixel=off Group=Stratum Markerattrs=( Symbol=PLUS) LegendLabel="Censored" NAME="SCATTER";
                layout overlay / width=30pct height=35pct halign=0.85 valign=0.85
walldisplay=(fill) opaque=false
yaxisopts=( Label=" " type=linear ) x2axisopts=(labelFitPolicy=Split) xaxisopts=(Label=" " labelFitPolicy=Split type=linear linearopts=( tickvaluelist=( 0 100 200 300 400 500 ) viewmin=0 viewmax=500 ) ) x2axisopts=(labelFitPolicy=Split);
            StepPlot X=Time Y=Survival / primary=true Group=Stratum LegendLabel="Survival Probability" NAME="s";
                    ScatterPlot X=Time Y=Censored / subpixel=off Markerattrs=( Symbol=PLUS) LegendLabel="Censored" NAME="c";
                    ScatterPlot X=Time Y=Censored / subpixel=off Group=Stratum Markerattrs=( Symbol=PLUS) LegendLabel="Censored" NAME="SCATTER";
                endlayout;
        endlayout;
        Layout Overlay /  walldisplay=none xaxisopts=(display=none griddisplay=off displaySecondary=none) x2axisopts=(display=none griddisplay=off displaySecondary=none);
                AxisTable Value=AtRisk X=tAtRisk / labelPosition=min Class=Stratum colorGroup=Stratum Display=(Label);
        endlayout;
    endlayout;
    endgraph;
end;
run;

*===5.渲染GTL;
ods html style=gghtml2;
proc sgrender data=SurvivalPlotData template=cars_inset_bar_bar;
run;

具体效果如下图：

我开通了视频号，也许以后分享视频教程就更方便了，欢迎围观。