王雨萌　孫瑞華　黃　傲　徐　凱　李　歡

?

·計算機(jī)應(yīng)用·

單因素多水平臨床試驗(yàn)定量指標(biāo)統(tǒng)計分析報表的SAS宏實(shí)現(xiàn)

王雨萌1孫瑞華2△黃傲1徐凱1李歡1

SAS軟件輸出的統(tǒng)計結(jié)果多而復(fù)雜，常需要人為將統(tǒng)計結(jié)果復(fù)制粘貼到統(tǒng)計報告相應(yīng)位置上，這一過程耗時耗力且容易出錯。為此，可根據(jù)需要編制相應(yīng)的SAS宏程序，直接生成針對各種類型資料的統(tǒng)計分析報表。以往對于成組計量或計數(shù)資料的統(tǒng)計分析報表的SAS宏程序已經(jīng)有了一定的探討[1-2]，本文主要介紹單因素多水平臨床試驗(yàn)定量指標(biāo)統(tǒng)計分析報表自動實(shí)現(xiàn)的SAS宏程序。

資料與統(tǒng)計方法

根據(jù)統(tǒng)計學(xué)原理，對于符合正態(tài)分布的定量指標(biāo)，通常采用均數(shù)、標(biāo)準(zhǔn)差、中位數(shù)、最小值、最大值以及95%置信區(qū)間來描述其集中和離散趨勢，而非正態(tài)分布的定量指標(biāo)，則采用均數(shù)、中位數(shù)、最小值、最大值以及四分位間距來進(jìn)行統(tǒng)計描述。

對于多水平定量指標(biāo)的組間比較，若滿足方差分析的前提條件，則采用方差分析處理資料；若指標(biāo)滿足正態(tài)性但不滿足方差齊性要求，則選用Welch近似方差分析，此時Welch近似方差分析的結(jié)果較方差分析的結(jié)果更穩(wěn)定；當(dāng)指標(biāo)既不滿足正態(tài)性也不滿足方差齊性要求時，則選用Kruskal-Wallis秩和檢驗(yàn)。方差分析與Welch近似方差分析的統(tǒng)計量用F_value表示，Kruskal-Wallis秩和檢驗(yàn)的統(tǒng)計量用H表示。為直接得到如表1所示的統(tǒng)計報表，編寫如下三組定量指標(biāo)統(tǒng)計分析的SAS宏程序。

表1　三組基線一般資料的組間比較(FAS)

定義宏及宏變量

在該宏程序中設(shè)置四個宏變量&database、&var、&index和&datasty，分別表示數(shù)據(jù)庫、統(tǒng)計量、指標(biāo)項(xiàng)和統(tǒng)計分析數(shù)據(jù)集，并假定分組變量為group。

%macro sanzujl(database=,var=,index=,datasty=);

生成描述性統(tǒng)計量及正態(tài)檢驗(yàn)結(jié)果

利用output語句將描述性統(tǒng)計量和正態(tài)檢驗(yàn)的結(jié)果輸出到SAS指定的數(shù)據(jù)集中，對該數(shù)據(jù)集進(jìn)行拆分后再橫向合并，使三組的描述性統(tǒng)計量和正態(tài)檢驗(yàn)的結(jié)果顯示在一行，方便后面的調(diào)用和輸出。

/*利用univariate程序輸出部分描述性結(jié)果到result1數(shù)據(jù)集*/

proc univariate noprint;

var &var;

by group;

where &datasty=1;

output out=result1 n=n mean=mean median=med std=std min=min max=max nmiss=nmiss qrange=qrange PROBn=Pnor;

/*利用means程序輸出95%置信區(qū)間到result2數(shù)據(jù)集*/

proc means clm noprint data=&database;

var &var;

by group;

where &datasty=1;

output out=result2(drop=_TYPE_ _FREQ_)uclm=uclm lclm=lclm;

run;

/*合并result1和result2為數(shù)據(jù)集result*/

data result;

merge result1 result2;

by group;

run;

/*將三組描述性結(jié)果分別輸出到獨(dú)立的數(shù)據(jù)集中，并定義變量的長度*/

data a(where=(group=“A”))b(where=(group=“B”))c(where=(group=“C”));

set result;

mean=round(mean,0.01);

med=round(med,0.01);

std=round(std,0.01);

Pnor=round(Pnor,0.01);

lclm=round(lclm,0.01);

uclm=round(uclm,0.01);

qrange=round(qrange,0.01);

min=round(min,0.01);

max=round(max,0.01);

run;

/*將三個獨(dú)立的數(shù)據(jù)集分別重新命變量名，再合并到result數(shù)據(jù)集中*/

data result(drop=group);

merge a(rename=(n=NA nmiss=MA mean=MeanA med=medA std=StdA qrange=qrangeA lclm=lclmA uclm=uclmA Min=MinA Max=MaxA Pnor=PnorA))

b(rename=(n=NB nmiss=MB mean=MeanB med=medB std=StdB qrange=qrangeB lclm=lclmB uclm=uclmB Min=MinB Max=MaxB Pnor=PnorB))

c(rename=(n=NC nmiss=MC mean=MeanC med=medC std=StdC qrange=qrangeC lclm=lclmC uclm=uclmC Min=MinC Max=MaxC Pnor=PnorC));

run;

生成方差分析及近似方差分析結(jié)果

利用SAS的ODS功能將方差分析、Welch近似方差分析以及方差齊性檢驗(yàn)的結(jié)果分別輸出到指定的數(shù)據(jù)集中，將所有有用變量橫向合并到一個數(shù)據(jù)集。

/*利用anova程序輸出方差分析結(jié)果到ModelANOVA數(shù)據(jù)集、輸出方差齊性檢驗(yàn)結(jié)果到HOVFTest1數(shù)據(jù)集、輸出welch近似方差分析結(jié)果到Welch1數(shù)據(jù)集*/

ods listing close;

ods output ModelANOVA=ModelANOVA HOVFTest=HOVFTest1 Welch=Welch1;

proc anova data=&database;

where &datasty=1;

class group;

model &var=group;

means group/hovtest welch;

run;

ods listing;

/*方差齊性檢驗(yàn)保留有用變量*/

data HOVFTest(keep=FValue ProbF);

set HOVFTest1;

if _n_ ^=1 then delete;

run;

/*welch近似方差分析保留有用變量*/

data Welch(keep=FValue ProbF);

set Welch1;

if _n_ ^=1 then delete;

run;

/*將方差分析、Welch近似方差分析、方差齊性檢驗(yàn)的結(jié)果合并到數(shù)據(jù)集test中,保留有用變量，添加方差分析與Welch近似方差分析統(tǒng)計量標(biāo)簽*/

data test(keep=FValue ProbF stat FValueH ProbFH FValueW ProbFW);

merge ModelANOVA(keep=FValue ProbF)HOVFTest(rename=(FValue=FValueH ProbF=ProbFH))Welch(rename=(FValue=FValueW ProbF=ProbFW));

stat=“F_value”;

run;

生成Kruskal-Wallis秩和檢驗(yàn)結(jié)果

利用SAS的ODS功能將秩和檢驗(yàn)的結(jié)果輸出到指定的數(shù)據(jù)集，將統(tǒng)計量和P值拆分后橫向合并到一個數(shù)據(jù)集。

/*利用npar1way wilcoxon程序輸出秩和結(jié)果到KruskalWallisTest數(shù)據(jù)集*/

ods listing close;

Ods Output KruskalWallisTest=KruskalWallisTest;

proc npar1way wilcoxon data=&database;

where &datasty=1;

var &var;

class group;

run;

ods listing;

/*只保留統(tǒng)計量和P值*/

data KruskalWallisTest1;

set KruskalWallisTest;

if Label1=“DF” then delete;

run;

/*將統(tǒng)計量和P值分別獨(dú)立保存*/

data d(where=(Name1=“_KW_”))e(where=(Name1=“P_KW”));

set KruskalWallisTest1;

run;

/*將秩和檢驗(yàn)的統(tǒng)計量和P值合并到kwttest數(shù)據(jù)集，并定義變量的長度和統(tǒng)計量的標(biāo)簽*/

data kwttest(keep=cValue1d cValue1e statkw);

merge d(rename=(Name1=Name1d cValue1=cValue1d nValue1=nValue1d))e(rename=(Name1=Name1e cValue1=cValue1e nValue1=nValue1e));

cValue1d=round(cValue1d,0.01);

label cValue1d=“Chi-Square”;

label cValue1e=“Pr > Chi-Square”;

statkw=“H”;

run;

輸出統(tǒng)計結(jié)果

將描述性統(tǒng)計量、正態(tài)檢驗(yàn)結(jié)果、方差分析、Welch近似方差分析以及秩和檢驗(yàn)結(jié)果的數(shù)據(jù)集橫向合并，生成打印的數(shù)據(jù)集，從該數(shù)據(jù)集中輸出特定條件下的描述性統(tǒng)計量和P值。

/*將數(shù)據(jù)集result、test、kwttest合并*/

data _null_;

merge result test kwttest;

/*如果三組均正態(tài)且方差齊，用方差分析*/

if PnorA>0.05 and PnorB>0.05 and PnorC>0.05 and ProbFH>0.05 then do;

file print notitle;

put #1 @3 “&index”

#2 @5 “例數(shù)(缺失)” @22 NA‘(‘MA’)’ @38 NB‘(‘MB’)’ @55 NC‘(‘MC’)’ @70 stat’=‘FValue @88 ProbF

#3 @5 “均數(shù)±標(biāo)準(zhǔn)差” @22 MeanA‘±’StdA @38 MeanB‘±’StdB @55 MeanC‘±’StdC

#4 @5 “中位數(shù)” @22 MedA @38 MedB @55 MedC

#5 @5 “最小值-最大值” @22 MinA‘-’MaxA @38 MinB‘-’MaxB @55 MinC‘-’MaxC

#6 @5 ‘95%CI’ @22 lclmA‘-’uclmA @38 lclmB‘-’uclmB @55 lclmC‘-’uclmC;

end;

/*如果三組均正態(tài)但方差不齊，用welch方差分析*/

else if PnorA>0.05 and PnorB>0.05 and PnorC>0.05 and ProbFH<0.05 then do;

file print notitle;

put #1 @3 “&index”

#2 @5 “例數(shù)(缺失)”@22 NA‘(‘MA’)’ @38 NB‘(‘MB’)’ @55 NC‘(‘MC’)’ @70 stat’=‘FValueW @88 ProbFW

#3 @5 “均數(shù)±標(biāo)準(zhǔn)差” @22 MeanA‘±’StdA @37 MeanB‘±’StdB @55 MeanC‘±’StdC

#4 @5 “中位數(shù)” @22 MedA @38 MedB @55 MedC

#5 @5 “最小值-最大值” @22 MinA‘-’MaxA @38 MinB‘-’MaxB @55 MinC‘-’MaxC

#6 @5 ‘95%CI’ @22 lclmA‘-’uclmA @38 lclmB‘-’uclmB @55 lclmC‘-’uclmC;

end;

/*如果三組有一組非正態(tài)，用秩和檢驗(yàn)*/

else if PnorA<0.05 or PnorB<0.05 or PnorC<0.05 then do;

file print notitle;

put #1 @3 “&index”

#2 @5 “例數(shù)(缺失)”@22 NA‘(‘MA’)’ @38 NB‘(‘MB’)’ @55 NC‘(‘MC’)’ @70 statkw’=‘cValue1d @88 cValue1e

#3 @5 “均數(shù)” @22 MeanA @38 MeanB @55 MeanC

#4 @5 “中位數(shù)” @22 MedA @38 MedB @55 MedC

#5 @5 “最小值-最大值” @22 MinA‘-’MaxA @38 MinB‘-’MaxB @55 MinC‘-’MaxC

#6 @5 “四分位間距” @22 qrangeA @38 qrangeB @55 qrangeC;

end;

run;

%mend sanzujl;

打印表頭及表格線

最后，打印出表頭及表格線。

%macro ctit(tit);

data _null_;

file print n=ps notitles;put #1 @5 “&tit”

#2 @2 94*‘_’

#3 @5 ‘指標(biāo)項(xiàng)’ @22 ‘A組’ @38 ‘B組’ @55 ‘C組’ @70 ‘統(tǒng)計量’ @88 ‘P值’

#4 @2 94*‘_’;

run;

%mend ctit;

%macro cleg;

data _null_;

file print n=ps notitles;

put #1 @2 94*‘_’;

run;

%mend cleg;

運(yùn)行宏程序

將以上所有SAS宏程序提交SAS系統(tǒng)運(yùn)行，即可完成如表1所示的統(tǒng)計分析表格。

/*其中，knfps為數(shù)據(jù)庫，設(shè)定數(shù)據(jù)集為FAS數(shù)據(jù)集，要分析比較的變量是weigh和high，分別表示體重和身高。*/

%ctit(表1.三組基線一般資料組間比較(FAS));

%sanzujl(database=knfps,var=weigh,datasty=FAS,index=體重);

%sanzujl(database=knfps,var=high,datasty=FAS,index=身高);

%cleg;

討　　論

目前無論美國FDA和國內(nèi)SFDA都要求采用世界公認(rèn)的統(tǒng)計分析軟件SAS來進(jìn)行統(tǒng)計分析，但是SAS軟件自身的統(tǒng)計分析結(jié)果比較豐富，通常只需要將其中一小部分關(guān)鍵的結(jié)果按照一定格式制作成表格[3]。盡管SAS提供了如proc report、proc template等制表的功能，但定義過程較為復(fù)雜，需要對SAS軟件有一定的掌握程度。本文介紹了一段簡單實(shí)用且可操作性非常強(qiáng)的SAS宏程序供廣大讀者調(diào)用，其運(yùn)行結(jié)果經(jīng)驗(yàn)證準(zhǔn)確而可靠，適合應(yīng)用于單因素多水平定量指標(biāo)的組間比較。

[1]鄒建東，熊寧寧，卜擎燕，等.正態(tài)分布定量指標(biāo)統(tǒng)計分析報表的SAS宏程序.中國臨床藥理學(xué)與治療學(xué)，2004，9(7)：838-840.

[2]鄒建東，熊寧寧，卜擎燕，等.四格表指標(biāo)統(tǒng)計分析報表的SAS宏程序.中國臨床藥理學(xué)與治療學(xué)，2005，10(3)：357-360.

[3]殷紅.臨床試驗(yàn)中統(tǒng)計分析報告自動化生成的研究與應(yīng)用.復(fù)旦大學(xué)，2009.

(責(zé)任編輯：劉壯)

孫瑞華，E-mail：sunruihua@263.net

1.北京中醫(yī)藥大學(xué)管理學(xué)院(100029)

2.中日友好醫(yī)院科研處