用Stata做logistic回归和Cox回归的限制性立方样图
最近看到有人询问限制性立方样图的问题，就去查资料试着做了一下，在网友的帮助下终于成功了，现将方法归纳如下：感谢网友@arikjin的知道，感谢@wangyiyisheng 提供的数据。下面以stata15进行如下分析：
有什么问题，尽管留言：
Logistic回归型:这段代码和命令均来自stata的xblc包：我们做的是关于HDL_CHOL对于心血管事件(CVD)的影响的
1载入数据：use "C:\Users\ALIENWARE\Desktop\data.dta" (路径根据自己的文件位置改变)。
2，由于在stata中做限制性立方样图要求小数型变量必须是float型，对于整数型可以省略这一步，我们导入数据后发现HDL_CHOL带小数点的double，那么我们把变量重新赋值到一个新变量中，generate var=HDL_CHOL此时var 就变成了HDL_CHO的替身float型。
3.设置条样变量的节点，我们可以利用四分位数来设置：
generate all=1
table all, c(freq p25 var p50 var p75 var p95 var)
产生如下结果：
----------------------------------------------------------------------
all | Freq. p25(var) med(var) p75(var) p95(var)
----------+-----------------------------------------------------------
1 | 14,757 1.03 1.27 1.53 2.04
----------------------------------------------------------------------
4.根据节点生成样条变量：
mkspline vars = var, knots(1.03 1.27 1.53 2.04) cubic
运行此命令后，打开数据集，可以看到产生了vars1-vars3三个变量。
5接下来就是进行回归过程，
logit CVD vars1 vars2 vars3 AGE SEX
如果需要调整其他因素就在后面继续加入新的变量，这里我们调整年龄和性别，其中CVD是因变量:
------------------------------------------------------------------------------
CVD | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
vars1 | -.5354012 .2771065 -1.93 0.053 -1.07852 .0077175
vars2 | 3.200325 1.901794 1.68 0.092 -.5271229 6.927772
vars3 | -6.922344 4.595648 -1.51 0.132 -15.92965 2.08496
AGE | .960004 .0288162 33.31 0.000 .9035253 1.016483
SEX | -.4745218 .0775497 -6.12 0.000 -.6265164 -.3225272
_cons | -7.379377 .3509652 -21.03 0.000 -8.067256 -6.691497
如果是做生存分析，就把生存分析的stcox过程来一遍即可：
stset SURVIVAL_T, failure(CVD==1)
stcox vars1 vars2 vars3 AGE SEX
生存分析结果：
failure _d: CVD == 1 结局变量
analysis time _t: SURVIVAL_T 生存时间
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
vars1 | .5047087 .1237369 -2.79 0.005 .3121446 .8160667
vars2 | 34.36153 57.55235 2.11 0.035 1.289405 915.7047
vars3 | .0005944 .0024037 -1.84 0.066 2.15e-07 1.645798
AGE | 2.669694 .0712306 36.80 0.000 2.533672 2.813018
SEX | .629417 .0432078 -6.74 0.000 .5501812 .7200642
6 产生给定一个具体的var值的OR或HR图：这一步可有可无
对于logistic回归和cox回归命令一样
xblc vars1-vars3, covname(var) at(0.21 0.93 1.27 1.53 2.04) reference(0.21) eform scatte
at：任意设置值，这个值一定要在变量var出现过，要升序排列，reference（）要设置at后面的最小值做参考：
7.产生立方样图：
产生美一个var值对应的OR或HR及可信区间，0.21-5.07是变量var最小值和最大值
levelsof var if inrange(var, 0.21,5.07)
xblc vars1-vars3, covname(var) at(`r(levels)') eform reference(0.21) line
可以把图做的更精细些：
需要先再运行一遍此命令
levelsof var if inrange(var, 0.21,5.07)
保存每个值对应的P、OR/HR lb ub: capture drop pa or lb ub
xblc vars1-vars3, covname(var) at(`r(levels)') eform reference(0.21) line generate(pa or lb ub)
增加标题：注意cox回归要把里面的or换成hr
twoway (rcap lb ub pa, sort) (scatter or pa, sort), legend(off) scheme(s1mono) ytitle("OR或HR可改") xtitle("横坐标可改") name(f1, replace)
最后生成的图片保存为PDF或者EPS，在PS里保存tif，就能达到300dpi了。
