临床预测模型常用统计算法SAS实现

本文为研究者手稿，具体参考 杨凯璇, 谷鸿秋. 临床预测模型常用统计模型及其SAS实现 [J/OL]. 中国卒中杂志, 2024, 19(5): 496.

随着大数据与精准医学时代的到来，临床预测模型在诊断、治疗决策、病人预后管理等方面的应用越来越多。在公共卫生领域，临床预测模型可以识别高危人群，针对患有或患病风险较高的人群进行预防干预；在临床实践中，临床预测模型可以告知患者及医生诊断或预后的概率，以便更加客观的做出医疗决策。开发临床预测模型，依据不同的预测结局类型，选择不同的模型类型，并通过统计软件来实现，对预测模型研究至关重要。本文将对临床预测模型中常用的统计模型及其SAS实现方法进行介绍以便广大研究者开展预测模型类研究。

1 如何选择模型

预测模型研究中，研究的目的（如诊断、预后），研究的设计类型（如横断面、队列）、预测结局的变量类型（连续、分类、生存）都会影响最终的模型类型的选择。其中，研究目的和设计类型决定了预测结局的变量类型，而结局变量的类型是模型类型选择的最终决定因素。因此，本文将按结局变量的类型来介绍如何进行模型的选择。

2常用预测模型及公式

结局类型可以从数据类型的角度出发，分为连续变量、分类变量以及时间-事件变量。不同的结局类型，决定了预测模型适用的模型类型，比如连续型结局变量常采用线性回归模型建模，二分类结局变量则采用logistic回归建模，而时间-事件变量则常用Cox模型建模。对于连续型结局变量，还应考虑是否正态分布，对于分类变量，还需具体区分是否多分类，是否有序；对于时间-事件变量，则需考虑是否存在竞争风险。具体的模型选择可见表1的总结。

2.1 连续变量

对于连续型结局变量，在选择模型类型时，还需要进一步考察其数据分布情况。若服从正态分布，选择线性回归模型即可。若为偏态分布，需要进一步确认其分布类型，如对数正态分布、伽马分布，泊松分布等，以便选择合适的模型。

2.1.1 正态分布

当结局为正态分布的连续变量时，一般选择线性回归模型,如预测体重、血压、生物学年龄等。例如，有案例利用大脑侧脑室周围及脑白质深部高信号体积预测脑神经血管年龄，由于脑神经血管年龄是一个大体符合则正态分布的连续变量，因此，研究者采用了线性回归来构建预测模型[10]。

在SAS中，线性回归可以通过PROC REG，PROC GLM，ROC GENMOD，PROC MIXED，以及PROC GLIMMIX等过程步实现。以sashelp中的“heart”数据集为例，用年龄（AgeAtStart），身高（height）来预测体重(weight)，将预测结果保存到数据集"pred",预测值变量名为"yhat"，SAS的实现代码如何下：

proc reg data=sashelp.heart;
   model weight=AgeAtStart height; 
   output out=pred predicted=yhat; 
quit;
 

proc glm data=sashelp.heart;
   class sex; 
   model weight=AgeAtStart height sex; 
   output out=pred predicted=yhat; 
quit; 

proc genmod data=sashelp.heart;
   class sex; 
   model weight=AgeAtStart height sex; 
  output out=pred predicted=yhat; 
quit; 

proc mixed data=sashelp.heart;
   class sex; 
   model weight=AgeAtStart height sex/outp=pred;
quit;

proc glimmix data=sashelp.heart;
   class sex; 
   model weight=AgeAtStart height sex; 
   output out=pred predicted=yhat; 
quit;

PROC REG的预测变量只能为连续变量。在 PROC REG, PROC GLM, PROC GENMOD和 PROC GLIMMIX中，均可使用 “output”语句中的out=和redicted=P 选项导出预测值；在 PROC MIXED中，需使用“model”语句中的outp=选项指定导出预测结果的数据集名称，默认预测值变量名为“pred”。PROC GENMOD 和PROC GLIMMIX 分别可以通过GEE或者混合效应模型处理群组数据。

2.1.2 对数正态分布

当结局为对数正态分布的连续变量时，一般选择连接函数为log，分布为正态分布的广义线性回归模型。比如手术时长，C反应蛋白、抗体滴度等血液指标等，其中位数小于均值，取值仅为正值，且异方差性（方差随期望增大而增大）也很常见，此类数据常符合对数正态分布。对于此类数据，研究者可能会采用先log 转换而后拟合简单线性回归，但考虑此方法的一些缺陷，比如：结局变量的尺度已经改变，log转换并不能同时改善非线性和异方差性等因素，因此建议选择广义线性模型，而非log 转换后的简单线性回归。

在SAS中，广义线性回归可使用 PROC GENMOD，``PROC GLIMMIX`过程步实现，以 sashelp 中的“heart”数据集为例，用年龄（AgeAtStart），身高（height），性别（sex）来预测胆固醇（Cholesterol），并将预测结果保存到数据集"pred",预测值变量名为"yhat"，SAS的实现代码如何下：

proc genmod data=sashelp.heart;
   class sex;
   model Cholesterol=AgeAtStart height sex/dist=normal link=log; 
   output out=pred predicted=yhat; 
quit;

proc glimmix data=sashelp.heart;
   class sex; 
   model Cholesterol=AgeAtStart height sex/dist=normal link=log; 
  output out=pred predicted(ilink)=yhat; 
quit;

在“model”语句中，通过“dist=normal”，“link=log”选项指定分布和连接函数; PROC GENMOD中可使用“output”语句中的 out= 和 predicted= 选项导出预测值；PROC GLIMMIX中可使用“output”语句中的out=和predicted(ilink)=选项导出预测值。

2.1.3 伽马分布

对于一些右偏态分布的连续变量（如医疗费用，住院时长等），伽马分布也是常用的一种分布类型。某研究在研究医保状态和总的医疗费用的关系时，分别尝试了基于log 转换后的最小二乘算法（OLS）和基于伽马分布及log 连接函数的广义线性回归，发现异方差性可以导致OLS模型中的偏差和低效，因此研究者最终使用了基于伽马分布及log 连接函数的广义线性回归[11]。此外，也有研究者基于伽马分布预测了乳腺癌药物反应[2]。

在SAS中，可使用 PROC GENMOD，``PROC GLIMMIX过程步实现，代码参考对数正态分布，只需将dist=normal改为dist=gamma`。

2.1.4 泊松分布

对于偏态的连续变量，除了上述的对数正态分布，伽马分布外，泊松和负二项分布也是常用的分布类型。泊松回归常用于单位时间/空间内某罕见事件发生次数的预测。对于住院天数，院内并发症的个数结局，也可以考虑采用泊松分布。泊松分布是负二项分布的一个特例，当泊松分布的方差大于均值时，采用负二项分布能更准确地描述数据的变异性。

某研究预测初次全膝关节置换术患者的住院时长时，考虑到数据的分布特性，分别尝试了泊松回归和负二项模型建立模型[3]。

在SAS中，同样可使用PROC GENMOD，PROC GLIMMIX过程步实现，代码参考对数正态分布，需将dist=normal改为dist=poisson。

2.2 分类变量

临床研究中的无明显时间属性的事件类结局为分类变量，如是否患病，患哪种类型的病。具体而言，分类变量可分为二分类、无序多分类以及有序多分类三种情况。

2.2.1 二分类变量

二分类结局变量是临床中最为常见的结局类型。比如，是否患有某种疾病，是否伤残等。针对二分类变量结局，logistic回归模型是拟合此类结果使用最为广泛的统计模型。Logistic 回归模型应用最为经典的预测模型场景便是诊断模型。例如，有研究者就基于多因素logistic回归模型构建了一个完全基于临床变量的肺栓塞诊断模型[4]。此外，对于一些短期内的事件结局和并发症的预测，也可以应用logistic构建预测模型。例如，某研究者，根据格拉斯哥结果量表评分，采用logistic回归模型构建了三个月后患者死亡及功能性残疾预测模型[13]；也有研究者利用住院资料基于logistic回归模型构建了糖尿病酮症的预测模型[12]。

在SAS中，logistic回归可使用 PROC LOGISTIC，PROC GENMOD, PROC GLIMMIX 过程步实现，以sashelp中的“heart”数据集为例，用年龄（AgeAtStart），性别（sex），舒张压（Diastolic），吸烟情况（Smoking）以及胆固醇（Cholesterol）来预测结局是否死亡（status），并将预测结果保存到数据集"pred",预测概率变量名为"phat"，SAS的实现代码如何下：

proc logistic data=sashelp.heart;
  class sex; 
  model status(event="Dead")=AgeAtStart sex Diastolic Smoking Cholesterol; 
  output out=pred predicted=phat; 
run; 

proc logistic data=sashelp.heart;
  class sex; 
  model status(event="Dead")=AgeAtStart sex Diastolic Smoking Cholesterol; 
  output out=pred predprobs=(individual); 
run;

proc genmod data=sashelp.heart;
  class sex;
  model status(event="Dead")=AgeAtStart sex Diastolic Smoking Cholesterol/ dist=binomial link=logit; 
  output out=pred predicted=phat; 
run;
 
proc glimmix data=sashelp.heart;
  class status sex; 
  model status(event="Dead")=AgeAtStart sex Diastolic Smoking Cholesterol/ dist=binary link=logit; 
  output out=pred predicted(ilink)=phat;
run;

PROC LOGISTIC中，在output语句中使用predprobs=(individual)选项可输出预测结局每个分类的预测概率，变量名为“IP_x”, x表示分类的取值；PROC GENMOD中，使用“model”语句，指定dist=binomial，link=logit 选项建模，使用“output”语句中的out=和predicted=选项导出预概率;PROC GLIMMIX中，使用“model”语句，指定 dist=binary，link=logit 选项建模，使用“output”语句中的out=和predicted(ilink)=选项导出预测概率。

2.2.2 无序多分类

无序多分类变量，常见于临床研究中的疾病分型。针对无序多分类变量结局，可拟合无序logistic回归，也称名义 logistic 回归。构建无序logistic 回归模型时，选择结局中的某一个类别作为参照，剩余各类别分别与选定的类别对比，构建模型。例如，有研究者想要对急诊科的发热儿童的诊断做出预测，诊断结局分为三个类别，为肺炎，其他细菌感染和无细菌感染。此结局为无序多分类变量，因此研究者使用了无序logistic回归构建了一个诊断模型[5]。

在SAS中，无序logistic回归可使用PROC LOGISTIC，PROC GLIMMIX过程步实现，以sashelp中的“heart”数据集为例，通过年龄（AgeAtStart），性别（sex），舒张压（Diastolic），吸烟情况（Smoking）以及胆固醇（Cholesterol）来预测死亡原因（DeathCause），并将预测结果保存到数据集"pred",预测概率变量名为"phat"，SAS的实现代码如何下：

proc logistic data=sashelp.heart;
  class sex DeathCause; 
  model DeathCause(event="Cancer")=AgeAtStart sex Diastolic Smoking Cholesterol/link=glogit; 
  output out=pred predicted=phat; 
run; 

proc logistic data=sashelp.heart;
  class sex DeathCause;
  model DeathCause(event="Cancer")=AgeAtStart sex Diastolic Smoking Cholesterol/link=glogit; 
  output out=pred predprobs=(individual); 
run;
 
proc glimmix data=sashelp.heart;
  class sex DeathCause; 
  model DeathCause(event="Cancer")=AgeAtStart sex Diastolic Smoking Cholesterol/ dist=multinomial link=glogit; 
  output out=pred predicted(ilink)=phat;
run;

PROC LOGISITC的“model”语句中，指定link=glogit建模，在“output”语句中，使用out=和predicted=或predprobs=(individual)选项导出预测概率。其中predprobs=(individual)会输出每个分类的预测概率。PROC GENMOD无法拟合无序多分类logistic回归。PROC GLIMMIX的predicted(ilink)= 可以获得每一个水平的预测概率，不过数据集形式是一个水平一条观测。

2.2.3 有序多分类

有序多分类变量，常见于临床研究中疾病的严重程度，比如心功能分级等,可拟合有序logistic回归。有序logistic回归模型是基于累计概率构建的回归模型。若结局变量包括G个等级，则需要构建j=G-1 个模型。比如，在缺血性卒中功能结局预后的研究中，通常采用改良的Rankin评分量表（0-5分）作为结局。大多数研究将mRS以2分为界值将结局转化为二分类结局，但也有研究者将mRS评分视为0-5的等级变量，采用有序logistic回归模型建立预测模型预测缺血性卒中患者出院90天后的mRS评分[6]。

在SAS中，有序logistic回归模型可使用PROC LOGISTIC，PROC GENMOD，PROC GLIMMIX过程步实现，以sashelp中的“heart”数据集为例，通过年龄（AgeAtStart），性别（sex），身高（height），体重（weight）来预测血压分级（bp_status），并将预测结果保存到数据集“pred”,预测概率变量名为“phat”，SAS的实现代码如何下：

proc logistic data=sashelp.heart;
  class sex; 
  model bp_status=sex height weight AgeAtStart; 
  output out=pred predprobs=(individual); 
run;

proc genmod data=sashelp.heart;
  class sex; 
  model bp_status=sex height weight ageatstart/dist=multinomial link=clogit; 
  output out=pred predicted=phat; 
run; 

proc glimmix data=sashelp.heart;
  class sex; 
  model bp_status=sex height weight ageatstart/dist=multinomial link=clogit; 
  output out=pred predicted(ilink)=phat; 
run;

PROC LOGISTIC中，在output语句中使用predprobs=(individual)选项可输出预测结局每个分类的预测概率，变量名为“IP_x”, x表示分类的取值；PROC GENMOD的predicted=和 PROC GLIMMIX的 predicted(ilink)=也可以实现，不过输出的是累积的概率，需要额外的计算方能获得每一个水平的概率。

2.3 时间-事件变量

临床研究中，很多事件类的结局不仅关注其是否出现了某种结局，还关注出现这些结局所经历的时间长短，此类结局多为时间-事件变量,采用生存分析的方法进行分析，另外，根据有无竞争风险事件，生存分析又分为存在和不存在竞争风险两种情况。

2.3.1 不存在竞争风险事件

对于不存在（或者不考虑）竞争风险事件的时间-事件变量，主要应用Cox比例风险模型(proportional hazards regression model)进行建模。不考虑竞争风险的时间-事件变量，在临床研究中非常普遍。例如，有研究者曾基于瑞士缺血性卒中患者信息，相开发建立了缺血性卒中后1年和5年的癫痫发作的预测模型，此预测模型便是基于不考虑竞争风险事件的Cox比例风险模型[7]。以及中国人群动脉粥样硬化性心血管疾病10年发病风险的预测模型，也是基于Cox比例风险模型开发而成[14]。

在SAS中，Cox回归模型可使用PROC PHREG过程步实现，以sashelp中的“heart”数据集为例，首先对“heart”数据集进行整理，“s”为不考虑竞争风险事件的生存结局，“s_c”为考虑竞争风险事件的生存结局。再通过性别（sex），身高（height），体重（weight），舒张压（Diastolic）等来预测生存结局（s），并将预测结果保存到数据集"pred",预测生存概率变量名为“_1_minus_phat”，SAS的实现代码如何下：

data heart;
 set sashelp.heart;
 s=(status="Dead"); /*s=1代表死亡;s=0代表未死亡*/
 if missing(DeathCause) then s_c=0;
 else if DeathCause="Cancer" then s_c=1;
 else if DeathCause="Cerebral Vascular Disease" then s_c=2;
 else if DeathCause="Coronary Heart Disease" then s_c=3;
 else s_c=4;  /*s_c：0=删失;1=结局事件;2,3,4=竞争风险事件)*/
run;

proc phreg data=heart;
  class sex;  
  model AgeAtDeath*s(0)= sex height weight Diastolic Smoking Cholesterol; 
  baseline covariates=heart out=pred survival=_1_minus_phat timelist=70; 
run;

在“baseline”语句中可使用out=选项指定输出结果的数据集名称，survival=选项指定输出结果数据集中生存概率的变量名称，“timelist”选项指定预测何时的生存概率。covariates=指定预测因子数据集，该SAS数据集需包含所有解释变量的集合，数据集中的所有变量都会复制到out=数据集中。

2.3.2 存在竞争风险事件

不考虑竞争风险事件，这其中有技术层面的因素（研究者对考虑竞争风险的统计模型相对陌生），也有现实的因素（竞争风险事件率低，对结果的影响小）。不过，现在越来越多的研究者开始考虑规避方法学上的瑕疵，采用考虑竞争风险事件的方法，目前常用的处理竞争风险事件的模型有两种：原因别Cox模型与Fine-Gray模型。

2.3.2.1 原因别Cox模型

原因别Cox（cause-specific Cox）模型的思路是分别对每一个原因别进行Cox回归，回归时将发生了竞争风险事件的研究对象归为删失[15]。例如，有研究者在开发原发性肺癌患者的10年内第二原发肺癌风险的预测模型时，由于原发性肺癌患者有相当大比例的患者在发展第二原发肺癌之前因合并症而死亡，因此研究者在开发预测模型时，采用了原因别Cox回归以处理竞争风险事件[8]。

在SAS中，原因别Cox模型可使用PROC PHREG过程步实现，以“heart”数据集为例，通过性别（sex），身高（height），体重（weight），舒张压（Diastolic）等来预测生存结局（s_c），并将预测结果保存到数据集"pred",结局累计发生率变量名为“phat”，SAS的实现代码如何下：

proc phreg data=heart;
  class sex;
  model AgeAtDeath*s_c(0)= sex height weight Diastolic Smoking Cholesterol/eventcode(cox)=1; 
   baseline covariates=heart out=pred cif=phat timelist=70; 
run;

在“model”语句中，使用“eventcode(cox)=”选项指定结局变量中关注的结局事件的编码；在“baseline”语句中可使用“out=”选项指定输出结果的数据集名称，“cif=”选项指定输出结果数据集中累计发生率的变量名称，“timelist”选项指定预测的时间点。“covariates=”指定数据集，该SAS数据集需包含所有解释变量的集合，数据集中的所有变量都会复制到“out=”数据集中。

2.3.2.2 Fine-Gray模型

处理竞争风险的模型更为常用的模型是 1999年Fine和Gray提出的次分布风险（subdistribution hazard）模型[16]，即Fine-Gray模型。不同于原因别cox模型将发生竞争风险事件的研究对象视为删失，Fine-Gray模型将发生竞争风险事件的研究对象仍纳入风险集，相当于为竞争事件定义了无限长的删失时间，或删失时间晚于最后观测失效的时间，保留了更多存在竞争风险的研究对象[17]。例如某研究在开发老年乳腺癌患者10年复发风险预测模型时，便分别使用Cox回归模型与Fine-Gray回归模型做对比，结果显示Cox回归模型高估了复发的发生率，在竞争风险较高的人群中，应用Fine-Gray模型更为适合[18]。又如，某研究在对重症监护室中患者患压力损伤（压疮、褥疮）的风险开发预测模型时，考虑到死亡为获得性压力损伤事件的竞争风险事件，研究者便采用了Fine-Gray模型[19]。

在SAS中，Fine-Gray模型可使用PROC PHREG过程步实现，以“heart”数据集为例，通过性别（sex），身高（height），体重（weight），舒张压（Diastolic）等来预测生存结局（s_c），并将预测结果保存到数据集"pred",结局累计发生率变量名为“phat”，SAS的实现代码如何下：

proc phreg data=heart;
  class sex;
  model AgeAtDeath*s_c(**0**)= sex height weight Diastolic Smoking Cholesterol/eventcode=1; 
  baseline covariates=heart out=pred cif=phat timelist=70; 
run;

在“model”语句中，使用eventcode=选项指定结局变量中结局的编码；在“baseline”语句中可使用out=选项指定输出结果的数据集名称，cif=选项指定输出结果数据集中累计发生率的变量名称，timelist选项指定预测的时间点。covariates=指定数据集，该SAS数据集需包含所有解释变量的集合，数据集中的所有变量都会复制到out=数据集中。

3 SAS实现总结

SAS作为专业的统计分析软件，支持各种统计模型，如一般线性回归模型（PROC GLM），广义线性回归模型（PROC GENMOD），logistic回归模型（PROC LOGISTIC），线性混合模型（PROC MIXED），广义线性混合模型（PROC GLIMMIX）等。我们可以使用这些PROC步建立预测模型，并通过指定的相应的选项（option）导出预测值或预测概率。各模型的SAS建模及预测结果导出代码总结见下表（表2）：

表2 SAS中不同结局类型的预测模型建模方法总结

4 总结

如何选择模型是建立临床预测模型的关键，预测结局的变量类型是决定模型类型的关键。若结局为连续变量, 则常用线性回归模型，对于偏态的连续变量，则可考虑广义线性模型，并选择合适的分布和连接函数；若结局为二分类变量，常选择Logistic回归模型；若结局为多分类变量，则依据是否有序，可以选择无序logistic回归，或者有序logistic回归；若结局为时间-事件变量，通常选择Cox回归模型。同时还应提前考虑是否存在竞争风险，若存在竞争风险，可选择原因别Cox模型或者Fine-Gray模型。

在SAS中可以通过PROC REG，PROC GLM，ROC GENMOD，PROC MIXED，以及PROC GLIMMIX等过程步，配合dist=，以及link=等选项，拟合常见的各种统计模型，并通过predicted=，predprobs=，outp=，survival=，cif=等选项将预测值输出到out=指定的数据集。

此外，除了传统的统计模型,一些基于机器学习的算法,如决策树、随机森林、神经网络等也在预测模型研究中得到广泛的应用，鉴于在数据维度有限情况下，机器学习的模型预测能力并不优于传统预测模型方法，因此本文未做展开描述，读者可查询相关资料。