判断自变量共线性,在SPSS 中可套用线性回归的共线性检验,看VIF<10即可。
打开D盘中“20190810北京临床预测模型’中的“03_logistic_regression‘’,点开“1126 training-218.csv”。此为上次SPSS中选出的70%作为预测模型,注意将变量名都转换为英文。
将所要用到的数据在桌面上建立一个文件夹:Stata,将其放入。更改保存的路径名:C:\Users\Administrator\Desktop。导入文件。
可进性相关的统计操作。
描述性统计:
连续型变量:sum X if Y=1,detail
正态性检验:sktest X;qladder X;gladder X;
卡方检验:病例对照:tab X Y,column chi e
队列研究:tab X Y,row chi e(为什么使用Fisher检验?不应该 满足T<1或者n<40?)
U检验:ranksum X,by(Y)
table1的快速制作
其操作所需的代码见桌面上的文件夹Stata中的"1126 training model",若出现Stata 16许可证过期的提示,将系统时间调回2018年便可打开。
ssc install table1:告诉软件给table1留个位子。
pwd:查看储存位置,可以修改保存自自己建立的文件夹。如:cd"C:\Users\Administrator\Desktop\Stata"
table1,by(yqol) vars(age cate\bmi cate\mouthopen conts\rmorph cate\rcurve cate\rnum cate\gender cate\presymp cate\smoke cate\pgl cate\time conts\pgv cate\winter cate\flap cate\procedure cate\experience cate) format(%8.1f) onecol missing cmissing test pdp(2)**saving(C:\Users\Administrator\Desktop\Stata/table1.xls)**或者直接saving(“table1,xls”)保存至默认路径。
代码说明:-format(%8.1f),cformat(%8.1f) 可分别设定对连续型变量及分类变量保留一位小数,括号中“%”是格式化代码,
“8”是包含整数、小数及符号在内的总宽度,“.”代表小数点,“1”是小数点后位数,“f”是“fixed”的缩写,
表示指定总宽度及保留小数位数(注意:小数点后位数一定要小于总宽度,此例中小数点后位数为“1” < 总宽度“8”,
平时可以根据研究需要更改这两个数字);
-onecol可让多分类变量的分组等级显示在此变量名下方,而非单独成列;
-missing可将分类变量的缺失例数及比例显示出来,cmissing可报告连续型变量非缺失例数;
-test可将我们对每个变量所使用的特定分组比较方法显示在表格中;
**-pdp(#)**可设定分组比较p值保留小数位数,如pdp(2)为保留两位小数。
建模
logit Y X
enter :加入全部变量,b3.pgl指以3为对照,1、2和3比(pgl水平分为1,2,3)
enter1:加入全部变量,i.pgl指以1为对照,3、2和1比
enter2 :留下自己感兴趣的变量,剔除一些无临床探究价值的变量。
逐步法:
stepwise, pe(0.2): logit yqol age bmi mouthopen rmorph rcurve rnum gender presymp smoke pgl time pgv winter flap procedure experience
est store f1
**pe代表向前:一个一个进,(0.2)指变量的p<0.2时变量纳入模型。pr代表向后:变量一开始全部进入模型,然后挨个往外踢,(0.02)指P>0.02时,将变量踢出。界值设定可以随意,0.1-0.9都可以尝试,然后比较AIC ,众多模型中最后AIC区域稳定的即为理想的模型。AIC值与变量个数的关系要看不同变量的AIC值的差异有无统计学差异,进而选择要AIC值小的还是变量个数少的。
est store 模型名 (储存模型) est stats 模型名(显示模型的AIC)
lrtest 模型1 模型2 比较两个模型的AIC值差异有无统计学意义。
注意,不能直接拿危险因素来做预测模型,危险因素更加关注P值大小,而预测模型注重β值,可称之为预测因子。预测模型文章的讨论不似危险因素逐个分析,不研究为什么这个变量有预测价值,而是从整体上辩证,建模后的验证。 (训练集,验证集)预测因子与危险因素截然不同。
最好样本量大,不然容易过拟合。过拟合即在训练集结果很好,拿到验证集不好。如ROC在训练集0.9,到了验证集0.6。
ROC曲线
主要语句:1.产生预测概率:gen pred0=1/(1+exp(-(β0 +X1β1+…)))
2.生成ROC表格(仅一个模型):roctab Y(结局变量) pred0
3.计算约登指数和临界值:roctab yqol pred0,detail,将结果输出为表格,在Excel中计算临界值即可。
4.生成其他文献构建的模型:gen pred1=1/(1+exp(-(-5.384181 +time .0395471 +presymp * 0.8765776 + pgl 0.6712864+smoke0.6849501 + age0 .5099984 + gender0.4002588)))
5.将自己构建的模型与已有文献的模型进行计较(两个模型):roccompY pred1 pred0,还可直观看到曲线:roccomp Y pred1 pred0,graph
6.显示全部或部分的预测概率:list yqol pred1;list yqol pred1 in 1/10(展示前十个)
7.更换已有模型:如replace pred0**=1/(1+exp(-(β0 +X1β1+…)))
8.若有三个模型需要互相比较:rocgold Y pred0 pred1 pred2。这条命令是指pred1和pred2均与pred0(金标准即参照)比较。
9.若有相关代码的疑惑,可用谷歌搜索,如roc? in stata,即可出现很多视频。
Hosmer-lemeshow检验
estat gof, group(10) table
三个地方愿意P值>0.05:正态性检验,方差齐性检验,Hosmer-lemeshow检验。
注:视频中讲解的hl命令目前还不能安装,但使用estat gof, group(10) table(注意要紧跟在logistics模型之后使用)结果相同。命令的来源参考
所有的命令也可通过窗口菜单操作。
用R进一步进行calibiration。将结果导出,保存为.csv格式。
本次用到的命令汇总:
描述性统计:
1.连续型变量:sum X if Y=1,detail
2.正态性检验:sktest X;qladder X;gladder X;
3.卡方检验:病例对照:tab X Y,column chi e
队列研究:tab X Y,row chi e(为什么使用Fisher检验?不应该 满足T<1或者n<40?)
4.U检验:ranksum X,by(Y)
ssc install table1*:告诉软件给table1留个位子。
5.pwd:查看储存位置,可以修改保存自自己建立的文件夹。如:cd"C:\Users\Administrator\Desktop\Stata"
table1,by(yqol) vars(age cate\bmi cate\mouthopen conts\rmorph cate\rcurve cate\rnum cate\gender cate\presymp cate\smoke cate\pgl cate\time conts\pgv cate\winter cate\flap cate\procedure cate\experience cate) format(%8.1f) onecol missing cmissing test pdp(2)**saving(C:\Users\Administrator\Desktop\Stata/table1.xls)**或者直接saving(“table1,xls”)保存至默认路径。
代码说明:-format(%8.1f),cformat(%8.1f) 可分别设定对连续型变量及分类变量保留一位小数,括号中“%”是格式化代码,
“8”是包含整数、小数及符号在内的总宽度,“.”代表小数点,“1”是小数点后位数,“f”是“fixed”的缩写,
表示指定总宽度及保留小数位数(注意:小数点后位数一定要小于总宽度,此例中小数点后位数为“1” < 总宽度“8”,
平时可以根据研究需要更改这两个数字);
-onecol可让多分类变量的分组等级显示在此变量名下方,而非单独成列;
-missing可将分类变量的缺失例数及比例显示出来,cmissing可报告连续型变量非缺失例数;
-test可将我们对每个变量所使用的特定分组比较方法显示在表格中;
**-pdp(#)**可设定分组比较p值保留小数位数,如pdp(2)为保留两位小数。
6.logit Y X
enter :加入全部变量,b3.pgl指以3为对照,1、2和3比(pgl水平分为1,2,3)
enter1:加入全部变量,i.pgl指以1为对照,3、2和1比
enter2 :留下自己感兴趣的变量,剔除一些无临床探究价值的变量。
逐步法:
stepwise, pe(0.2): logit yqol age bmi mouthopen rmorph rcurve rnum gender presymp smoke pgl time pgv winter flap procedure experience
est store f1
**pe代表向前:一个一个进,(0.2)指变量的p<0.2时变量纳入模型。pr代表向后:变量一开始全部进入模型,然后挨个往外踢,(0.02)指P>0.02时,将变量踢出。界值设定可以随意,0.1-0.9都可以尝试,然后比较AIC ,众多模型中最后AIC区域稳定的即为理想的模型。AIC值与变量个数的关系要看不同变量的AIC值的差异有无统计学差异,进而选择要AIC值小的还是变量个数少的。
7.est store 模型名 (储存模型) est stats 模型名(显示模型的AIC)
lrtest 模型1 模型2 比较两个模型的AIC值差异有无统计学意义。
8.ROC:主要语句:
8.1.产生预测概率:gen pred0=1/(1+exp(-(β0 +X1β1+…)))
8. 2.生成ROC表格(仅一个模型):roctab Y(结局变量) pred0
8. 3.计算约登指数和临界值:roctab yqol pred0,detail,将结果输出为表格,在Excel中计算临界值即可。
8.4.生成其他文献构建的模型:gen pred1=1/(1+exp(-(-5.384181 +time .0395471 +presymp * 0.8765776 + pgl 0.6712864+smoke0.6849501 + age0 .5099984 + gender0.4002588)))
8.5.将自己构建的模型与已有文献的模型进行计较(两个模型):roccompY pred1 pred0,还可直观看到曲线:roccomp Y pred1 pred0,graph
8. 6.显示全部或部分的预测概率:list yqol pred1;list yqol pred1 in 1/10(展示前十个)
8.7.更换已有模型:如replace pred0**=1/(1+exp(-(β0 +X1β1+…)))
8. 8.若有三个模型需要互相比较:rocgold Y pred0 pred1 pred2。这条命令是指pred1和pred2均与pred0(金标准即参照)比较。
8.9.若有相关代码的疑惑,可用谷歌搜索,如roc? in stata,即可出现很多视频。
9.Hosmer-lemeshow检验:estat gof, group(10) table