SAS軟件實(shí)現(xiàn)Serfling回歸模型

寧波市疾病預(yù)防控制中心(315010)

李 寧 張 良 紀(jì) 威 韓航濤 樊明鎖 俞延峰

【提 要】 目的 為方便研究人員使用Serfling回歸模型做疾病負(fù)擔(dān)研究，開發(fā)SAS宏程序用于相關(guān)數(shù)據(jù)整理、模型應(yīng)用和結(jié)果展示。方法 以寧波市2015年1月1日到2018年12月31日0～3歲兒童流感和肺炎(ICD-10：J09～J18)門診數(shù)據(jù)為實(shí)例，介紹和驗(yàn)證SAS宏程序%SERFLING(IN_DSN=,AGE=,SEX=,CLIP=,WEEKS_EPI=)。結(jié)果 按照SAS宏程序%SERFLING的參數(shù)說明，在指定輸入數(shù)據(jù)集(IN_DSN=)基礎(chǔ)上按需選擇年齡范圍(AGE=)、性別(SEX =)、是否將不足一周的頭尾數(shù)據(jù)刪除(CLIP =)以及是否使用傳統(tǒng)Serfling回歸模型(WEEKS_EPI =)，就可以自動(dòng)完成對(duì)原始數(shù)據(jù)的整理、Serfling回歸模型應(yīng)用和相關(guān)結(jié)果展示。結(jié)論 SAS宏程序%SERFLING方便了Serfling回歸模型的應(yīng)用，具有一定的實(shí)用性。

1963年Serfling博士首次提出Serfling回歸模型，并且用于肺炎和流感導(dǎo)致的超額死亡研究[1]，之后該模型被廣泛應(yīng)用于各類季節(jié)性傳染病導(dǎo)致的超額發(fā)病率、住院率和死亡率等研究[2-4]。2015年Xiaoli Wang等提出調(diào)整Serfling回歸模型，該模型通過迭代方式自動(dòng)排除大于閾值的按周統(tǒng)計(jì)數(shù)，解決了人為判斷疾病流行期的困難，而且使擬合效果更好[4]。

Serfling回歸模型的基本思想簡(jiǎn)單明了，但在文獻(xiàn)中尚未查到相應(yīng)SAS程序。本研究設(shè)計(jì)了名為%SERFLING的SAS宏程序，用于Serfling回歸模型的數(shù)據(jù)整理、模型應(yīng)用和結(jié)果展示。該SAS宏程序只需指定輸入數(shù)據(jù)集(IN_DSN=)以及可選的幾個(gè)參數(shù)如年齡范圍(AGE=)、性別(SEX=)、流行期(WEEKS_EPI=)等，就可自動(dòng)完成對(duì)傳統(tǒng)或調(diào)整Serfling回歸模型的應(yīng)用，提高了工作效率。

資料與方法

1.資料

通過寧波市區(qū)域衛(wèi)生信息平臺(tái)，收集57家醫(yī)療機(jī)構(gòu)2015年1月1日至2018年12月31日0～3歲兒童的流感和肺炎門診數(shù)據(jù)(ICD-10：J09～J18)共計(jì)243202例，收集字段包括門診時(shí)間(DATE)、性別(SEX)和年齡(AGE)，如表1所示。

表1 數(shù)據(jù)集TEST部分觀測(cè)

2.方法

Serfling回歸方程如下：

Yt=a+bt+ct2+dt3+fsin(2πt/52)+gcos(2πt/52)+hsin(4πt/52)+icos(4πt/52)+et

其中Yt代表第t周統(tǒng)計(jì)數(shù)，t代表連續(xù)的周次(t=1,2,3,…)，a代表截距；b、c、d、f、g、h、i為回歸系數(shù)，其中a+bt+ct2+dt3反映長(zhǎng)期趨勢(shì)，fsin(2πt/52)+gcos(2πt/52)+hsin(4πt/52)+icos(4πt/52)反映季節(jié)周期性；et是誤差項(xiàng)。模型的擬合優(yōu)度通過決定系數(shù)R2判斷。

傳統(tǒng)Serfling回歸模型需要人為刪除流行期對(duì)應(yīng)統(tǒng)計(jì)數(shù)來擬合基線[1-2]。調(diào)整Serfling回歸模型通過計(jì)算機(jī)迭代來排除流行期數(shù)據(jù)，第一次迭代先用全部按周統(tǒng)計(jì)數(shù)擬合Serfling回歸模型，從全部按周統(tǒng)計(jì)數(shù)中刪除大于擬合值的數(shù)據(jù)，第二次迭代用上一次迭代保留的按周統(tǒng)計(jì)數(shù)擬合Serfling回歸模型，重新從全部按周統(tǒng)計(jì)數(shù)中刪除大于擬合值95%可信區(qū)間上限的按周統(tǒng)計(jì)數(shù)，第三次迭代用上一次迭代保留的按周統(tǒng)計(jì)數(shù)擬合Serfling回歸模型，重新從全部按周統(tǒng)計(jì)數(shù)中刪除大于擬合值95%可信區(qū)間上限的按周統(tǒng)計(jì)數(shù)，不斷重復(fù)直到R2下降，取R2最大的Serfling回歸模型來擬合基線[5-6]。

3.程序

/*SAS宏程序%SERFLING，參數(shù)說明見表2*/

%MACRO SERFLING(IN_DSN=,AGE=,SEX=,CLIP=,WEEKS_EPI=);

/*數(shù)據(jù)集IN_DSN按DATE排序*/

PROC SORT DATA=&IN_DSN OUT=TEMP_1SELECT;

BY DATE;RUN;

/*如果指定AGE=，則篩選相應(yīng)年齡范圍*/

/*如果未指定AGE=，則直接跳轉(zhuǎn)%SKIP_AGE:*/

%IF"&AGE"=""%THEN%GOTO SKIP_AGE;

DATA _NULL_;/*提取參數(shù)AGE=相關(guān)信息*/

IF FIND("&AGE",′-′)THEN DO;

AGE_LOWER=INPUT(SCAN("&AGE",1,′-′),BEST.);

AGE_UPPER=INPUT(SCAN("&AGE",2,′-′),BEST.);

END;ELSE DO;

AGE_LOWER=&AGE;AGE_UPPER=&AGE;END;

CALL SYMPUT(′AGE_LOWER′,AGE_LOWER);

CALL SYMPUT(′AGE_UPPER′,AGE_UPPER);RUN;

DATA TEMP_1SELECT;/*篩選相應(yīng)年齡范圍*/

SET TEMP_1SELECT;

WHERE &AGE_LOWER<=AGE<=&AGE_UPPER;RUN;

%SKIP_AGE:

/*如果指定SEX=，則篩選相應(yīng)性別*/

/*如果未指定SEX=，則直接跳轉(zhuǎn)%SKIP_SEX:*/

%IF "&SEX"=""%THEN %GOTO SKIP_SEX;

DATA _NULL_;/*提取參數(shù)SEX=相關(guān)信息*/

IF 0 THEN SET TEMP_1SELECT;

SEX_VTYPE=VTYPE(SEX);

CALL SYMPUT(′SEX_VTYPE′,SEX_VTYPE);RUN;

DATA TEMP_1SELECT;/*篩選相應(yīng)性別*/

SET TEMP_1SELECT;WHERE

%IF &SEX_VTYPE=N %THEN SEX=&SEX;

%ELSE SEX="&SEX";RUN;

%SKIP_SEX:

/*對(duì)數(shù)據(jù)按周匯總整理*/

PROC TIMESERIES DATA=TEMP_1SELECT OUT=TEMP_2WEEK;

ID DATE INTERVAL=WEEK ACCUMULATE=TOTAL;

VAR COUNT;RUN;

/*如果指定CLIP=YES，則根據(jù)實(shí)際修剪頭尾*/

%IF"&CLIP"=""%THEN %GOTO SKIP_CLIP;

PROC SQL NOPRINT;

SELECT MAX(DATE),MIN(DATE)

INTO:DATE_MAX,:DATE_MIN

FROM &IN_DSN;QUIT;

DATA TEMP_2WEEK;/*根據(jù)實(shí)際修剪頭尾*/

SET TEMP_2WEEK END=LAST;

IF _N_=1 AND DATE^=&DATE_MIN THEN DELETE;

IF LAST AND DATE^=&DATE_MAX-6 THEN DELETE;RUN;

%SKIP_CLIP:

/*為SERFLING回歸方程生成新變量*/

DATA TEMP_3SERFLING;SET TEMP_2WEEK;

WEEK=WEEK(DATE);/*用于傳統(tǒng)SERFLING回歸模型*/

COUNT_EXCLUDE=COUNT;

T=_N_;/*生成連續(xù)周次*/

T_B=T;T_C=T*T;T_D=T*T*T;

T_F=SIN(2*CONSTANT(′PI′)*T/52);

T_G=COS(2*CONSTANT(′PI′)*T/52);

T_H=SIN(4*CONSTANT(′PI′)*T/52);

T_I=COS(4*CONSTANT(′PI′)*T/52);RUN;

/*如果未指定WEEKS_EPI=，則使用調(diào)整SERFLING回歸模型*/

/*如果指定WEEKS_EPI=，則直接跳轉(zhuǎn)%SKIP_ADJUST:，使用傳統(tǒng)SERFLING回歸模型*/

%IF"&WEEKS_EPI"^=""%THEN %GOTO SKIP_ADJUST;

/*第一次迭代*/

PROC REG DATA=TEMP_3SERFLING PLOTS=NONE OUTEST=TEMP_4FITSTAT RSQUARE;/*輸出R2*/

MODEL COUNT_EXCLUDE=T_B T_C T_D T_F T_G T_H T_I / SELECTION=MAXR;

OUTPUT OUT=TEMP_6PREDICTED PREDICTED=PREDICTED;/*輸出擬合值*/

RUN;QUIT;

DATA _NULL_;SET TEMP_4FITSTAT;

CALL SYMPUT(′R2′,_RSQ_);RUN;

DATA TEMP_6PREDICTED;SET TEMP_6PREDICTED;

THRESHOLD=PREDICTED;RUN;/*閾值*/

/*第N次迭代*/

%NEXT_SERFLING:

DATA TEMP_6PREDICTED_NEW(DROP=PREDICTED:UCL:LCL:);

SET TEMP_6PREDICTED;COUNT_EXCLUDE=COUNT;

IF COUNT>THRESHOLD THEN CALL MISSING(COUNT_EXCLUDE);RUN;

ODS OUTPUT SELPARMEST=TEMP_5PARM_NEW;

PROC REG DATA=TEMP_6PREDICTED_NEW PLOTS=NONE OUTEST=TEMP_4FITSTAT_NEW RSQUARE;

MODEL COUNT_EXCLUDE=T_B T_C T_D T_F T_G T_H T_I / SELECTION=MAXR;

OUTPUT OUT=TEMP_6PREDICTED_NEW PREDICTED=PREDICTED UCL=UCL LCL=LCL;/*輸出擬合值95%CI上限*/

RUN;QUIT;

ODS OUTPUT CLOSE;

DATA _NULL_;SET TEMP_4FITSTAT_NEW;

CALL SYMPUT(′R2_NEW′,_RSQ_);RUN;

%IF &R2_NEW<=&R2 %THEN %GOTO NEXT_PLOT;/*不斷重復(fù)直到R2下降*/

%LET R2=&R2_NEW;

DATA TEMP_4FITSTAT;SET TEMP_4FITSTAT_NEW;RUN;

DATA TEMP_5PARM;SET TEMP_5PARM_NEW;RUN;

DATA TEMP_6PREDICTED;SET TEMP_6PREDICTED_NEW;

THRESHOLD=UCL;RUN;/*閾值*/

%GOTO NEXT_SERFLING;/*繼續(xù)迭代*/

%SKIP_ADJUST:

/*使用傳統(tǒng)SERFLING回歸模型*/

DATA _NULL_;/*提取參數(shù)WEEKS_EPI=相關(guān)信息*/

FORMAT WEEKS $1000.;WEEKS=0;

DO I=1 TO COUNT("&WEEKS_EPI",′ ′)+1;

WEEK=SCAN("&WEEKS_EPI",I,′ ′);

IF FIND(WEEK,′-′)THEN DO;

WEEK_LOWER=SCAN(WEEK,1,′-′);

WEEK_UPPER=SCAN(WEEK,2,′-′);

WEEKS=CATX(′ ′,WEEKS,′OR′,WEEK_LOWER,′<=WEEK<=′,WEEK_UPPER);

END;ELSE WEEKS=CATX(′ ′,WEEKS,′OR′,′WEEK=′,WEEK);END;

CALL SYMPUT(′WEEKS′,WEEKS);RUN;

DATA TEMP_3SERFLING;SET TEMP_3SERFLING;

IF &WEEKS THEN CALL MISSING(COUNT_EXCLUDE);RUN;

ODS OUTPUT SELPARMEST=TEMP_5PARM;

PROC REG DATA=TEMP_3SERFLING PLOTS=NONE OUTEST=TEMP_4FITSTAT RSQUARE;

MODEL COUNT_EXCLUDE=T_B T_C T_D T_F T_G T_H T_I / SELECTION=MAXR;

OUTPUT OUT=TEMP_6PREDICTED PREDICTED=PREDICTED UCL=UCL LCL=LCL;

RUN;QUIT;ODS OUTPUT CLOSE;

%NEXT_PLOT:

/*刪除無統(tǒng)計(jì)學(xué)意義的項(xiàng)后，重新擬合模型*/

PROC SORT DATA=TEMP_5PARM;BY DESCENDING STEP;RUN;

DATA TEMP_5PARM;SET TEMP_5PARM;RETAIN MAX_STEP;

IF _N_=1 THEN MAX_STEP=STEP;IF STEP=MAX_STEP;RUN;

DATA _NULL_;LENGTH REGRESSORS $1000;

SET TEMP_5PARM END=LAST;RETAIN REGRESSORS ′′;

IF FIND(VARIABLE,′Intercept′)THEN DO;

IF PROBF<0.05 THEN CALL SYMPUT(′INTERCEPT′,′′);

ELSE CALL SYMPUT(′INTERCEPT′,′NOINT′);END;

ELSE IF PROBF<0.05 THEN REGRESSORS=CATX(′ ′,REGRESSORS,VARIABLE);

IF LAST THEN CALL SYMPUT(′REGRESSORS′,REGRESSORS);RUN;

ODS OUTPUT FITSTATISTICS=TEMP_4FITSTAT PARAMETERESTIMATES=TEMP_5PARM;

PROC REG DATA=TEMP_6PREDICTED(DROP=PREDICTED:UCL:LCL:)PLOTS=NONE;

MODEL COUNT_EXCLUDE=®RESSORS / &INTERCEPT;

OUTPUT OUT=TEMP_6PREDICTED PREDICTED=PREDICTED UCL=UCL LCL=LCL;

RUN;QUIT;ODS OUTPUT CLOSE;

/*結(jié)果展示*/

ODS HTML IMAGE_DPI=300;ODS GRAPHICS /NOBORDER;

PROC SGPLOT DATA=TEMP_6PREDICTED;

SERIES X=DATE Y=COUNT;

SERIES X=DATE Y=PREDICTED / LINEATTRS=(PATTERN=DOT);

SERIES X=DATE Y=UCL / LINEATTRS=(PATTERN=SHORTDASH);

XAXIS GRID;YAXIS GRID;RUN;

/*刪除無關(guān)數(shù)據(jù)集*/

PROC DATASETS LIB=WORK NOLIST;

DELETE TEMP_4FITSTAT_NEW TEMP_5PARM_NEW TEMP_6PREDICTED_NEW;RUN;QUIT;

/*超額病例數(shù)*/

DATA TEMP_6PREDICTED;SET TEMP_6PREDICTED;

IF COUNT>UCL THEN EPI=1;RUN;

DATA TEMP_6PREDICTED;SET TEMP_6PREDICTED;

BY EPI NOTSORTED;

IF FIRST.EPI AND LAST.EPI THEN CALL MISSING(EPI);

IF EPI=1 THEN DO;EXCESS=COUNT-PREDICTED;

LCL_EXCESS=COUNT-UCL;UCL_EXCESS=COUNT-LCL;END;RUN;

/*匯總SERFLING回歸模型擬合結(jié)果*/

DATA TEMP_6PREDICTED;SET TEMP_6PREDICTED;

IF EPI=1 THEN COUNT_EPI=COUNT;

IF MDY(7,1,2015)<=DATE

ELSE IF DATE

ELSE IF DATE

PROC SQL;CREATE TABLE TEMP_7SUMMARY AS

SELECT YEAR,SUM(EPI)AS EPI_WEEKS,

SUM(EXCESS)AS EXCESS,SUM(LCL_EXCESS)AS LCL,SUM(UCL_EXCESS)AS UCL,

SUM(EXCESS)/SUM(COUNT_EPI)AS EXCESS_EPI,

SUM(LCL_EXCESS)/SUM(COUNT_EPI)AS LCL_EPI,

SUM(UCL_EXCESS)/SUM(COUNT_EPI)AS UCL_EPI,

SUM(EXCESS)/SUM(COUNT)AS EXCESS_ALL,

SUM(LCL_EXCESS)/SUM(COUNT)AS LCL_ALL,

SUM(UCL_EXCESS)/SUM(COUNT)AS UCL_ALL

FROM TEMP_6PREDICTED

GROUP BY YEAR;RUN;QUIT;

%MEND SERFLING;

表2 宏程序%SERFLING參數(shù)說明

結(jié) 果

對(duì)以上代碼進(jìn)行編譯后就可以根據(jù)表2參數(shù)說明調(diào)用SAS宏程序%SERFLING，例如：%SERFLING(IN_DSN=TEST,AGE=0-3,CLIP=YES)。WORK邏輯庫(kù)依次生成7個(gè)臨時(shí)數(shù)據(jù)集：TEMP_1SELECT、TEMP_2WEEK、TEMP_3SERFLING、TEMP_4FITSTAT、TEMP_5PARM和TEMP_6PREDICTED、TEMP_7SUMMARY。

原始數(shù)據(jù)集TEST(表1)包含4個(gè)變量：日期變量DATE、性別變量SEX、年齡變量AGE、計(jì)數(shù)變量COUNT，由于是個(gè)案形式，所有COUNT=1；數(shù)據(jù)集TEMP_1SELECT根據(jù)參數(shù)AGE=0-3，完成對(duì)原始數(shù)據(jù)集TEST的篩選，保留0～3歲兒童觀測(cè)；數(shù)據(jù)集TEMP_2WEEK在數(shù)據(jù)集TEMP_1SELECT基礎(chǔ)上對(duì)COUNT變量按周匯總，并且刪除監(jiān)測(cè)不足一周的頭尾數(shù)據(jù)(CLIP=YES)；數(shù)據(jù)集TEMP_3SERFLING根據(jù)Serfling回歸方程Yt=a+bt+ct2+dt3+fsin(2πt/52)+gcos(2πt/52)+hsin(4πt/52)+icos(4πt/52)+et構(gòu)建第一次迭代所需數(shù)據(jù)，變量COUNT_EXCLUDE即方程中的Yt，變量T為連續(xù)周次即方程中的t，T_B、T_C、T_D、T_F、T_G、T_H、T_I分別表示方程中的t、t2、t3、sin(2πt/52)、cos(2πt/52)、sin(4πt/52)、cos(4πt/52)，變量WEEK表示根據(jù)變量DATE計(jì)算的當(dāng)年周次(用于傳統(tǒng)Serfling回歸模型)；數(shù)據(jù)集TEMP_4FITSTAT記錄了最終的R2值；數(shù)據(jù)集TEMP_5PARM(表3)記錄了Serfling回歸模型最終的參數(shù)估計(jì)；數(shù)據(jù)集TEMP_6PREDICTED(表4)根據(jù)最終的Serfling回歸模型計(jì)算擬合值(PREDICTED)、擬合值的95%可信區(qū)間上限(UCL)以及在此基礎(chǔ)上判斷的流行周(EPI，定義為連續(xù)2周以上實(shí)際按周統(tǒng)計(jì)數(shù)超過擬合值的95%可信區(qū)間上限，即COUNT>UCL)、超額數(shù)(EXCESS，定義為流行周期間實(shí)際數(shù)與擬合值的差值，即EXCESS=COUNT-PREDICTED)；圖1以線圖方式直觀展示實(shí)際觀察值、基線和流行閾值，其中實(shí)線表示實(shí)際按周統(tǒng)計(jì)數(shù)(COUNT)，點(diǎn)線表示擬合值(PREDICTED)，短虛線表示擬合值的95%可信區(qū)間上限(UCL)；數(shù)據(jù)集TEMP_7SUMMARY(表5)在數(shù)據(jù)集TEMP_6PREDICTED基礎(chǔ)上，根據(jù)自定義的流行年度(YEAR，例如將每年7月至次年6月作為一個(gè)監(jiān)測(cè)年度)匯總流行周數(shù)(EPI_WEEKS)、超額數(shù)(EXCESS)及其95%可信區(qū)間(LCL、UCL)、超額數(shù)占流行周觀察總數(shù)比例(EXCESS_EPI)及其95%可信區(qū)間(LCL_EPI、UCL_EPI)、超額數(shù)占全年度觀察總數(shù)比例(EXCESS_ALL)及其95%可信區(qū)間(LCL_ALL、UCL_ALL)，經(jīng)過簡(jiǎn)單整理后可得表6。

表3 % SERFLING輸出結(jié)果：數(shù)據(jù)集TEMP_5PARM部分變量和觀測(cè)

結(jié)合圖1和表6可以發(fā)現(xiàn)，2017/2018年度流感和肺炎流行周中，0～3歲兒童超額就診病例數(shù)及其占全年度觀察病例數(shù)比例與往年相比有明顯上升，且流行持續(xù)17周，高于前兩個(gè)監(jiān)測(cè)年度，提示可能出現(xiàn)病原變化，需高度關(guān)注。

表4 %SERFLING輸出結(jié)果：數(shù)據(jù)集TEMP_6PREDICTED部分變量和觀測(cè)

*：變量COUNT_EXCLUDE中的“-”表明該周病例數(shù)(COUNT)未用于最后的Serfling回歸方程；變量EPI中的“-”表明該周最終被判斷為非流行周。

圖1 %SERFLING輸出結(jié)果：實(shí)際觀察值、基線和流行閾值

年度EPI_WEEKSEXCESSLCLUCLEXCESS_EPILCL_EPIUCL_EPIEXCESS_ALLLCL_ALLUCL_ALL55767 4248 7285 0.4062 0.2992 0.5132 0.1810 0.1333 0.2287 2015/201663677 1855 5498 0.2700 0.1362 0.4038 0.0662 0.0334 0.0989 2016/2017136415 2452 10378 0.3571 0.1365 0.5778 0.0744 0.0284 0.1203 2017/20181713297 8136 18458 0.3649 0.2233 0.5065 0.1955 0.1196 0.2713

表6 2014-2015至2017-2018年度寧波市0～3歲兒童流感和肺炎超額病例

討 論

Serfling回歸模型常用于流感的季節(jié)性和疾病負(fù)擔(dān)研究[2-3,7]，評(píng)估流感的季節(jié)性流行或大流行對(duì)人群健康的影響，對(duì)流感防控具有重要指導(dǎo)意義。近年來該模型也用于其他病毒性疾病的研究[6]，可為病毒研究和疾病防控提供依據(jù)。鑒于目前無Serfling回歸模型程序?qū)崿F(xiàn)的相關(guān)文獻(xiàn)，本文設(shè)計(jì)了名為%SERFLING的SAS宏程序，集數(shù)據(jù)整理、模型應(yīng)用和結(jié)果展示于一體，方便了Serfling回歸模型的使用，也提高了工作效率，具有一定的實(shí)用性。