生物医药领域三大回归(线性、logistic和Cox)其实本质还是围绕“一根线”,然而很多的时候,自变量x与y是非线性的关系,如何研究呢?三次样条就是近年比较火的一种方法,常见于剂量-反应关系的研究中。 本期松哥借用网络上的数据,带着大家一起实现下: 研究HDL_CHOL(高密度脂蛋白)对心血管事件(CVD)的影响;其中HDL_CHOL为连续性数值变量,CVD为(0-1)变量。采用Logistic回归进行分析。(1)通过菜单-导入数据,选取数据的格式,本例是xls格式。
(2)选择文件所在位置,请选择将第一行当做变量名读入。
(1)由于在stata中做限制性立方样图要求小数型变量必须是float型,对于整数型可以省略这一步,我们导入数据后发现HDL_CHOL带小数点的double,那么我们把变量重新赋值到一个新变量中,generate var=HDL_CHOL此时var 就变成了HDL_CHO的替身float型。(2)分析需要利用xblc函数,因此需要提前安装,否则提示无法识别。输入findit xblc,在弹出的页面中,选择st0215_1,打开后选择安装即可。(1)设置条样变量的节点,我们可以利用四分位数来设置: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 ---------------------------------------------------------------------- (2)根据节点生成样条变量: mkspline vars = var, knots(1.03 1.27 1.53 2.04) cubic 运行此命令后,打开数据集,可以看到产生了vars1-vars3三个变量。 (3)进行回归过程 logit CVD vars1 vars2 vars3 AGE SEX BMI 如果需要调整其他因素就在后面继续加入新的变量,这里我们调整年龄,性别,BMI,其中CVD是因变量: 如果是做生存分析,就把生存分析的stcox过程来一遍即可:
stset SURVIVAL_T, failure(CVD==1) stcox vars1 vars2 vars3 AGE SEX BMI 生存分析结果: failure _d: CVD == 1 结局变量 analysis time _t: SURVIVAL_T 生存时间 (4)产生给定一个具体的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后面的最小值做参考: (5)产生立方样图: 产生每一个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 (6)图形修饰: 可以把图做的更精细些,需要先再运行一遍此命令 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了。
|