用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了。