stata两种方法制作限制立方条图

我们前面已经讲了怎么使用R语言制作限制立方条图，今天来讲讲怎么使用stata制作限制立方条图，参考了网上常见的两种做法，顺便比较一下做图出来的差异，然后在和R语言比较一下。
首先要导入stata的xblc包，我们可以打入命令

search xblc, net

安装后以后就可以制作了限制立方条图了，我们继续使用我们的乳腺癌数据来制作先把数据导入

我们先来看看数据：
age表示年龄，pathsize表示病理肿瘤大小（厘米），lnpos表示腋窝淋巴结阳性，histgrad表示病理组织学等级，er表示雌激素受体状态，pr表示孕激素受体状态，status结局事件是否死亡，pathscat表示病理肿瘤大小类别（分组变量），ln_yesno表示是否有淋巴结肿大，time是生存时间，后面的agec是我们自己设定的，不用管它。
假设我们想看年龄和生存率的关系：
第一种这种是参考丁香园版主统小计的做法，
首先使用使用四分位数来设置样条的节点

generate all=1
table all, c(freq p25 age p50 age p75 age p95 age)//生成四分位数
mkspline age = age, knots(46 56 67 78) cubic//根据上面的四分位进行设置节点，46 56 67 78)就是这么来的

由此得到age1-age3等3个节点

然后建立COX回归模型，如果你是其他模型的话就在这一步建立其他模型

通过菜单和命令都可以建模，我就不详细说了，可以看我的其他文章有详细说明

stset time, failure(status==1)
stcox age1 age2 age3 ln_yesno
stcox age1 age2 age3 i.ln_yesno

然后通过下面命令得到一个年龄和HR的表

xblc age1-age3, covname(age) at(35 40 51 60 79) reference(50) eform scatte

这一步可有可无，主要是要产生一个年龄和HR的表，at可以选任意值，我主要是针对年龄50上下自由选的，但这个一定是要在年龄变量里面有的的年龄。

接下来做出立方条图

levelsof age if inrange(age, 35,79)
xblc age1-age3, covname(age) at(`r(levels)') eform reference(50) line

由上图可知，我们的参考点是50岁，为什么要选50岁呢，由下面表可知50岁的时候HR刚好为1，曲线上和曲线下对它有参考作用，这和我们使用R语言做出来的参考点是一致的。

我们还可以美化一下图形，先要生成pa hr lb ub这几个变量

levelsof age if inrange(age, 35,79)
xblc age1-age3, covname(age) at(`r(levels)') eform reference(50) line generate(pa hr lb ub)

然后做图，这个命令我要好好解释一下，twoway (line hr pa, sort) (line lb ub pa, sort lc(black black)这里是做出三条线，black black的意思是两条边线是黑色的，lp(- -)代表两条边线使用虚线表示

twoway (line hr pa, sort) (line lb ub pa, sort lc(black black) lp(- -)),legend(off) yscale(log) ylab(0.25 .5 1 2 5) xtitle(Age) ytitle(Hazard ratio) yline(1) title(Effect of age) subtitle (Adjusted for treatment)

接下来讲另一种方法，其实主要是在生存样条函数这里不同，是通过mkspline函数生成了样条函数

mkspline spl_age=age, cubic nknots(4)

这样就生成了spl_age1-- spl_age3等3个样条函数，也可以生成4个，我这里是为了好比较只生成3个

建模的过程也是一样的，我就简单点使用代码带过了
建模

stset time, failure(status==1)
stcox spl_age1 spl_age2 spl_age3 i.ln_yesno

做图，reference(50)也是和上面一样，从预测表里面选出

quietly levelsof age, local(ages)
levelsof age, local(ages)
xblc spl_age*, covname(spl_age1) at(`ages') reference(50) eform gen(ptage hr lb ub)
twoway (line hr ptage, sort) (line lb ub ptage, sort lc(black black) lp(- -)),legend(off) yscale(log) ylab(0.25 .5 1 2 5) xtitle(Age) ytitle(Hazard ratio) yline(1) title(Effect of age) subtitle(Adjusted for treatment)

OK，两个图都做出来了，我们进行比较一下

可以看出两种stata做出来的拐点基本一致，都是60多这样，50岁前的置信区间稍有不同，基本上两个图很接近，。我们再来比较一下和R做出来的图有上面不同

首先stata和R做出来的图型都是U型曲线，只是R做出来的U弯的没有这么明显，R的置信区间明显宽过stata,应该是算法不同，其次拐点基本一致，都是60多岁这样，因此两个软件做出来的结论应该基本一致，怎么取舍就看个人了。但是论图形的美观度，R就完爆stata了，ggplot2做图能力太强了。
更多精彩文章请关注公众号：零基础说科研