dxy logo
首页丁香园病例库全部版块
搜索
登录

用Stata做logistic回归和Cox回归的限制性立方样图

医疗行业从业者 · 最后编辑于 2018-08-24 · IP 湖南湖南
4.4 万 浏览
推荐、丁当奖励 2 项荣誉
这个帖子发布于 6 年零 339 天前,其中的信息可能已发生改变或有所发展。

最近看到有人询问限制性立方样图的问题,就去查资料试着做了一下,在网友的帮助下终于成功了,现将方法归纳如下:感谢网友@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三个变量。

 

img


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

 

img


可以把图做的更精细些:

需要先再运行一遍此命令

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)

 

img

 

 


 

最后生成的图片保存为PDF或者EPS,在PS里保存tif,就能达到300dpi了。




data.xls (3.84 MB)
203 247 41

全部讨论(0)

默认最新
avatar
203
分享帖子
share-weibo分享到微博
share-weibo分享到微信
认证
返回顶部