王崇旭 王詩(shī)遠(yuǎn) 劉 沛

東南大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(210009)

疫苗臨床試驗(yàn)中統(tǒng)計(jì)圖形的自動(dòng)化輸出

王崇旭 王詩(shī)遠(yuǎn) 劉 沛△

東南大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(210009)

SAS作為當(dāng)前主流的統(tǒng)計(jì)分析軟件而廣泛應(yīng)用于臨床試驗(yàn),其統(tǒng)計(jì)分析功能全面、可靠，可滿足臨床試驗(yàn)統(tǒng)計(jì)工作各方面的要求,并且被世界各國(guó)藥物評(píng)審機(jī)構(gòu)所認(rèn)可。然而SAS的繪圖功能相對(duì)薄弱，編制通用和具有特點(diǎn)的SAS繪圖宏程序，通常需要花費(fèi)較多的時(shí)間[1]。對(duì)于臨床試驗(yàn)中一些具有特點(diǎn)的圖，如試驗(yàn)流程圖、不良事件火山圖，在SAS中難以實(shí)現(xiàn)或者需要編寫(xiě)非常復(fù)雜的繪圖語(yǔ)句才能實(shí)現(xiàn)。R軟件的ggplot2包繪圖功能十分強(qiáng)大，和SAS相比它的繪圖更加簡(jiǎn)單靈活。此外R軟件中有豐富的繪圖包，可以方便快捷地實(shí)現(xiàn)各種繪圖功能。因此，兼顧SAS強(qiáng)大的數(shù)據(jù)處理功能和R軟件靈活的繪圖功能，如果能在疫苗臨床試驗(yàn)中聯(lián)合SAS與R，利用SAS繪圖宏程序，能夠方便快捷地實(shí)現(xiàn)美觀、高質(zhì)素圖形的批量自動(dòng)化輸出，從而提高臨床試驗(yàn)統(tǒng)計(jì)分析的質(zhì)量。

R軟件ggplot2包簡(jiǎn)介

ggplot2是一個(gè)用來(lái)繪制統(tǒng)計(jì)圖形的R軟件包。和其他繪圖軟件不同，ggplot2采用了圖層的設(shè)計(jì)方式，以一種結(jié)構(gòu)化的方式來(lái)設(shè)計(jì)和構(gòu)建圖形。和R的基礎(chǔ)圖形系統(tǒng)相比，ggplot2的功能更加強(qiáng)大，這源于其背后特定的圖形語(yǔ)法。ggplot2的圖形語(yǔ)法由一系列獨(dú)立的圖形部件組成，并能以許多種不同的方式組合起來(lái)，獨(dú)立的圖形部件使得圖形的可擴(kuò)展性很強(qiáng)，可以創(chuàng)建有創(chuàng)新性的、新穎的圖形類(lèi)型。ggplot2可以繪制出美觀的圖形，同時(shí)能減少繪圖需要的參數(shù)，避免諸多繁瑣的細(xì)節(jié)，例如添加圖例等[2]。用ggplot2繪圖時(shí)，圖形的每個(gè)部分可以依次進(jìn)行構(gòu)建，之后還可以進(jìn)行編輯。ggplot2精心挑選了一系列預(yù)設(shè)圖形，因此大部分情形下可以快速地繪制出許多高質(zhì)量的圖形，無(wú)需花費(fèi)太多時(shí)間調(diào)整圖形外觀。如果在格式上還有額外的要求，也可以利用ggplot2中的主題系統(tǒng)進(jìn)行定制。此外，使用ggplot2中特有的函數(shù)ggsave()可以將圖形存儲(chǔ)為高分辨率的多種圖形格式。

圖形自動(dòng)化輸出實(shí)現(xiàn)方法

1.SAS和R的連接

由于SAS和R各自有很多的版本，不同版本的SAS支持與相應(yīng)版本的R的連接，因此在連接前要保證SAS版本和R軟件版本的可連接性，本文使用SAS9.4版本和R3.2.5版本進(jìn)行連接。

首先在SAS軟件的安裝路徑中找到sasv9.cfg文件，用記事本打開(kāi)后，在最后一行添加“-rlang”選項(xiàng)，這個(gè)選項(xiàng)用于控制在SAS里面運(yùn)行R代碼的權(quán)限。然后利用SAS進(jìn)行數(shù)據(jù)處理和統(tǒng)計(jì)分析，將處理好的數(shù)據(jù)和統(tǒng)計(jì)分析結(jié)果輸出為csv文件。之后通過(guò)file語(yǔ)句和put語(yǔ)句，將R軟件的繪圖腳本程序輸出為“.r”后綴的外部文件。最后在SAS中用IML模塊提供R的接口[3]，用%include語(yǔ)句調(diào)用R軟件的程序腳本，即可實(shí)現(xiàn)SAS與R的連接。

2.高分辨率圖形的批量自動(dòng)化輸出

在疫苗臨床試驗(yàn)統(tǒng)計(jì)分析過(guò)程中，常常需要批量生成某種統(tǒng)計(jì)圖形，如抗體滴度逆向累積分布圖和不良事件火山圖。利用SAS生成疫苗臨床試驗(yàn)中幾種常用圖形的宏程序，將原始數(shù)據(jù)集名、變量名列表、SAS數(shù)據(jù)步中的篩選條件等設(shè)置成宏參數(shù)，將SAS數(shù)據(jù)處理和統(tǒng)計(jì)分析過(guò)程產(chǎn)生的一些變量設(shè)置成宏變量，在SAS數(shù)據(jù)處理及統(tǒng)計(jì)分析宏文本和R繪圖宏文本中以宏變量的形式引用已定義的宏參數(shù)。通過(guò)重復(fù)調(diào)用宏程序，每次調(diào)用時(shí)定義不同的宏參數(shù)，即可實(shí)現(xiàn)多種篩選條件下某種統(tǒng)計(jì)圖形的批量自動(dòng)化輸出。

為了滿足統(tǒng)計(jì)分析報(bào)告或雜志出版圖形印刷的要求，通常需要輸出高質(zhì)素的圖形。ggplot2包中特有的ggsave()函數(shù)通過(guò)dpi參數(shù)控制圖形的分辨率，ggsave()函數(shù)還可以將圖形保存成多種格式，如png、pdf。在SAS中每次調(diào)用圖形宏程序時(shí)都會(huì)運(yùn)行一次ggsave()函數(shù)，通過(guò)多次調(diào)用宏程序，即可實(shí)現(xiàn)高質(zhì)素圖形在外部文件中的批量保存。

疫苗臨床試驗(yàn)中圖形自動(dòng)化輸出實(shí)例

1.試驗(yàn)流程圖

在臨床試驗(yàn)統(tǒng)計(jì)分析報(bào)告中，除文字、表格描述外，應(yīng)采用流程圖的方式描述受試者的分布情況，包括篩選例數(shù)、篩選失敗例數(shù)及原因、參與隨機(jī)化的例數(shù)、各組脫落或剔除受試者的例數(shù)、百分比等，以及方案偏離情況、各分析數(shù)據(jù)集的分布等。

通常情況下，人們習(xí)慣于在菜單式操作軟件中采用插入形狀的方式繪制流程圖的框架，然后將SAS輸出的例數(shù)匯總結(jié)果手動(dòng)填寫(xiě)到流程圖的文本中。這種做法一方面比較繁瑣，繪制不同項(xiàng)目的流程圖時(shí)不僅需要多次編寫(xiě)匯總例數(shù)的SAS程序，而且繪制流程圖時(shí)需要在軟件之間反復(fù)切換。另一方面在手動(dòng)填寫(xiě)例數(shù)時(shí)很容易出現(xiàn)人為錯(cuò)誤，影響統(tǒng)計(jì)分析質(zhì)量。如果能夠在SAS中生成例數(shù)匯總的宏程序，將各種分布情況的例數(shù)計(jì)算出來(lái)并定義為宏變量，在R中用程序性的語(yǔ)言繪制流程圖并在繪圖程序腳本中引用定義的宏變量，通過(guò)連接 SAS和R，即可實(shí)現(xiàn)流程圖的自動(dòng)化輸出，這樣不僅能夠避免人為錯(cuò)誤，而且大大提高流程圖的制作效率[4]。

R中有專(zhuān)門(mén)繪制流程圖的diagram包，繪圖的流程如下：(1)載入繪制流程圖所需的軟件包：shape、diagram和Cairo，其中shape包和diagram包用于繪制流程圖，Cairo包用于高質(zhì)量圖形渲染。(2)用CairoPNG()設(shè)定圖形的保存路徑及尺寸。(3)用openplotmat(sub="",family=" ",cex.sub= )語(yǔ)句設(shè)定圖形的標(biāo)題及標(biāo)題字體和大小。(4)用elpos<-coordinates (c(1,2,1,2,2,1))語(yǔ)句設(shè)定流程圖元素框的布局。(5)用treearrow(from=elpos[#,],to=elpos[#,],arr.pos =,arr.length=)語(yǔ)句設(shè)定不同元素框間連接線的位置和相關(guān)屬性。(6)用textrect(elpos[#,],lab=′脫落&n.人′,family=" ",cex=)語(yǔ)句設(shè)定每個(gè)元素框中的內(nèi)容及屬性，將引用宏變量的文本填入元素框。圖1為通過(guò)上述方法自動(dòng)化輸出的某個(gè)疫苗臨床試驗(yàn)的流程圖，試驗(yàn)流程圖描述了臨床試驗(yàn)中受試者的分布情況。

圖1 試驗(yàn)流程圖

2.逆向累積分布圖

逆向累積分布圖是一種用來(lái)顯示免疫原型數(shù)據(jù)分布的圖形工具，特別適合于疫苗組間分布的可視性圖形比較[5]。該圖的橫坐標(biāo)表示抗體滴度或濃度的對(duì)數(shù)轉(zhuǎn)換值，縱坐標(biāo)表示至少有此抗體水平的個(gè)體所占的百分比。以23價(jià)肺炎球菌多糖疫苗Ⅲ期臨床試驗(yàn)為例，批量生成試驗(yàn)組和對(duì)照組免前免后各個(gè)血清型抗體滴度的逆向累積分布圖。

23價(jià)肺炎球菌多糖疫苗由23種血清型肺炎球菌多糖抗原組成，能夠?qū)τ?3種最常見(jiàn)血清型引起的肺炎球菌感染性疾病產(chǎn)生保護(hù)。為了直觀地對(duì)試驗(yàn)組和對(duì)照組免前免后每種血清型抗體滴度進(jìn)行比較，通過(guò)連接SAS與R，利用SAS逆向累積分布圖宏程序，實(shí)現(xiàn)試驗(yàn)組和對(duì)照組免前免后每種血清型抗體滴度的批量自動(dòng)化輸出。圖2為FAS集血清型1免前抗體逆向累積分布圖，從圖中可以直觀地看出FAS集試驗(yàn)組和對(duì)照組血清型1免前抗體分布情況。實(shí)現(xiàn)步驟如下：

(1)用SAS進(jìn)行數(shù)據(jù)處理，將原始數(shù)據(jù)處理成繪圖需要的數(shù)據(jù)形式。

(2)用“proc export”將處理好的數(shù)據(jù)導(dǎo)出為“.csv”文件。

(3)將R的繪圖程序腳本保存在“.r”的文件中，以便在SAS調(diào)用。具體程序如下：

data _null_;

file ′C:submit.r′;/* 保存繪圖程序腳本 */

put ′submit /R;′;

put ′csv=read.csv("F:\csv.csv",header=T)′;/* R讀入SAS導(dǎo)出的數(shù)據(jù) */

put ′library(ggplot2)′;/* 載入ggplot2包和grid包 */

put ′library(grid)′;

put "p <- ggplot(data=csv,mapping=aes(x=&AntiBody.,y=comulative))";

put "p + geom_point(aes(colour =&Grplab.,size=0.5 ))+";

put "scale_x_continuous(trans=′log′,limits=c(1,&length.),breaks=c(1,4,16,64,256,512,1024,2048,4096))+";

put ′labs(x="抗體滴度",y="累積頻率(%)",title="′ "&title."""′)+ ′;/* 繪制逆向累積分布圖,并設(shè)定繪圖參數(shù) */

put "theme( axis.title.x=element_text(size=16,colour=′black′,face=′bold′),axis.title.y=element_text(size=16,colour=′black′,face=′bold′))+";/* 主題系統(tǒng)精雕細(xì)琢 */

put "ggsave(′C:\&title..png′,width=5,height=7,dpi=700)";/* 保存高質(zhì)素的圖形 */

put "endsubmit;";

run;

(4)在SAS IML模塊中用%include語(yǔ)句調(diào)用R的程序腳本文件：

proc iml;

%include "C:submit.r";

quit;

(5)將上述數(shù)據(jù)處理和圖形自動(dòng)化輸出全過(guò)程封裝成帶參數(shù)的逆向累積分布圖宏程序。

(6)調(diào)用帶參數(shù)的逆向累積分布圖的宏程序，實(shí)現(xiàn)各血清型免前免后抗體逆向累積分布圖的批量自動(dòng)化輸出。下面是調(diào)用宏程序繪制血清型1免前抗體逆向累積分布圖的SAS程序：

%stone(Datain=crf,Group=grp,AntiBody=M1_MQ,Condition=%str(fas=1),title=FAS集血清型1免前抗體逆向累積分布曲線);/* 輸出FAS集免前血清型1 抗體滴度的逆向累積分布圖 */

圖2 FAS集血清型1免前抗體逆向累積分布曲線

3.不良事件火山圖

篩選組間差異有統(tǒng)計(jì)學(xué)意義的不良事件，明確其與研究藥物或疫苗的關(guān)系是臨床試驗(yàn)安全性評(píng)價(jià)的主要目的之一。在統(tǒng)計(jì)學(xué)中，火山圖是通過(guò)散點(diǎn)圖的形式快速識(shí)別大型數(shù)據(jù)集中同類(lèi)型數(shù)據(jù)的組間差異[6]。和統(tǒng)計(jì)分析表格和文字描述相比，不良事件火山圖能夠幫助研究者快速直觀地發(fā)現(xiàn)組間差異有統(tǒng)計(jì)學(xué)意義的不良事件。圖3是用R軟件ggplot2包繪制的疫苗臨床試驗(yàn)中的總體不良事件火山圖。圖中x軸為不良事件組間發(fā)生率差(試驗(yàn)組-對(duì)照組)，y軸為-log10(P)。對(duì)P值進(jìn)行負(fù)對(duì)數(shù)變換,一方面是為了讓P值變化更加平穩(wěn)，另一方面是為了讓組間率差小而P值大的點(diǎn)和組間率差大而P值小的點(diǎn)呈現(xiàn)中間低兩側(cè)高類(lèi)似于火山噴發(fā)的分布。P值是根據(jù)不同的不良事件資料由SAS軟件進(jìn)行卡方檢驗(yàn)、校正卡方檢驗(yàn)、Fisher確切概率法得到的，虛線為-log10(0.05)，即在虛線上方的不良事件組間差異有統(tǒng)計(jì)學(xué)意義。點(diǎn)的大小代表不良事件的發(fā)生率，組間發(fā)生率差異越大，點(diǎn)的顏色越深。從圖3中可以看出疼痛和腫這兩種不良事件發(fā)生率組間差異有統(tǒng)計(jì)學(xué)意義，其中腫的發(fā)生率試驗(yàn)組小于對(duì)照組(P<0.05)，疼痛的發(fā)生率試驗(yàn)組大于對(duì)照組(P<0.05)。

在疫苗臨床試驗(yàn)中，不僅需要分析總體不良事件發(fā)生率，還要分析不同針次、不同部位、不同等級(jí)的不良事件發(fā)生率，如第1針接種后不良反應(yīng)發(fā)生率、局部不良反應(yīng)發(fā)生率、全身不良反應(yīng)發(fā)生率、1級(jí)不良反應(yīng)發(fā)生率、第1針接種后局部1級(jí)不良反應(yīng)發(fā)生率等等。通過(guò)連接SAS與R，利用SAS不良事件火山圖宏程序，實(shí)現(xiàn)各種分析條件下的不良事件火山圖的批量自動(dòng)化輸出，實(shí)現(xiàn)步驟如下：

圖3 不良事件火山圖

(1)用SAS進(jìn)行數(shù)據(jù)處理和統(tǒng)計(jì)分析，將統(tǒng)計(jì)分析結(jié)果處理成繪圖需要的數(shù)據(jù)形式。

(2)用“proc export”將處理好的數(shù)據(jù)導(dǎo)出為“.csv”文件。

(3)將R的繪圖程序腳本保存在“.r”的文件中，以便在SAS調(diào)用。程序如下：

data _null_;

file "C:submitae.r";/* 保存繪圖程序腳本 */

put "submit /R;";

put ′csv=read.csv("F:\csv2.csv",header=T)′;/* R讀入SAS導(dǎo)出的數(shù)據(jù) */

put ′library(ggplot2)′;/* 載入ggplot2包和grid包 */

put ′library(grid)′;

put "gred<-scale_color_gradient2(low=′yellow′,high=′purple′,limits =c(&min.,&max.),midpoint=0)";/* 設(shè)置點(diǎn)的顏色 */

put ′q <- ggplot(data=csv,mapping=aes(x=d,y=-log10(p)))′;

put ′q +geom_point(aes(size=p3,color=d,alpha=abs(d)))+′;

put "scale_x_continuous(limits=c(&min.,&max.)) +";

put ′scale_y_continuous(limits=c(-0.5,2))+′;/* 繪制散點(diǎn)圖并設(shè)置x軸和y軸的范圍 */

put "geom_hline(aes(yintercept=-log10(0.05)),colour=′red′,linetype=′longdash′)+";/* 繪制p=0.05處的水平參照線 */

put ′geom_text(aes(x=d,y=-log10(p),label=a,size=p3/10,fontface=2))+′;/* 在點(diǎn)內(nèi)填入相應(yīng)的不良事件描述 */

put ′labs(x="RD(發(fā)生率差%)",y=expression( -log[10](P.Value)),title="′ " & title." ′",color="RD")+′;/* 設(shè)置x軸、y軸的標(biāo)題和圖形標(biāo)題 */

put "theme(axis.text.x=element_text(size=12),axis.text.y=element_text(size=12))+";/*主題系統(tǒng)精雕細(xì)琢*/

put "ggsave(′C:\&title..png′,width=7,height=5,dpi=700) ";/* 保存高質(zhì)素的圖形 */

put "endsubmit;";

(4)在SAS IML模塊中用%include語(yǔ)句調(diào)用R的程序腳本文件：

proc iml;

%include "C:submitae.r";

quit;

(5)將上述數(shù)據(jù)處理、統(tǒng)計(jì)分析和圖形自動(dòng)化輸出全過(guò)程封裝成帶參數(shù)的不良事件火山圖宏程序。

(6)調(diào)用帶參數(shù)的不良事件火山圖的宏程序，實(shí)現(xiàn)不同針次、不同部位、不同等級(jí)的不良事件火山圖的批量自動(dòng)化輸出。下面是調(diào)用宏程序繪制疫苗接種后1級(jí)不良反應(yīng)火山圖的SAS程序：

%stone(Grp1=′試驗(yàn)組′,Grp2=′對(duì)照組′,Condition=%str(serious=1),title=疫苗接種后1級(jí)不良反應(yīng)火山圖 );/* 輸出疫苗接種后嚴(yán)重程度為1級(jí)的不良事件火山圖 */

討 論

SAS軟件擅于數(shù)據(jù)處理和統(tǒng)計(jì)分析，R軟件擅長(zhǎng)繪制多種圖形且ggplot2包繪圖簡(jiǎn)單靈活、功能強(qiáng)大，在疫苗臨床試驗(yàn)中聯(lián)合SAS和R，可以實(shí)現(xiàn)統(tǒng)計(jì)軟件之間的互補(bǔ)，而且連接過(guò)程簡(jiǎn)單易操作。雖然SAS 9.2開(kāi)始，ODS Graphic System的出現(xiàn)，使得SAS在統(tǒng)計(jì)繪圖美觀性和便捷性方面比SAS/Graph模塊有所改觀，而且SAS提供GTL(graphic template language)用于定制繪圖模板，但是定制具有特色的統(tǒng)計(jì)圖模板需要比較復(fù)雜的SAS語(yǔ)法且對(duì)SAS編程能力要求很高。和SAS相比，R軟件ggplot2包繪圖簡(jiǎn)單靈活，采用圖層語(yǔ)法簡(jiǎn)單易懂，而且R中有一些獨(dú)特的包和函數(shù)，可以方便快捷地實(shí)現(xiàn)某些功能。在醫(yī)學(xué)研究和藥物臨床試驗(yàn)中聯(lián)合SAS與R分析數(shù)據(jù)在國(guó)內(nèi)外已有討論[7]。和之前的研究相比，本研究做了如下幾方面的改進(jìn)和探索：(1)之前的研究通過(guò)聯(lián)合SAS與R實(shí)現(xiàn)單個(gè)統(tǒng)計(jì)圖形的輸出，本研究從實(shí)際應(yīng)用角度出發(fā)，利用SAS宏程序?qū)崿F(xiàn)批量統(tǒng)計(jì)圖形的自動(dòng)化輸出。(2)ggplot2包作為R軟件繪圖的一大特色，在前人的研究中并未提及，本文對(duì)ggplot2包做了介紹并在實(shí)際應(yīng)用中使用ggplot2繪制統(tǒng)計(jì)圖形。(3)簡(jiǎn)化了SAS與R連接和圖形輸出的方法，直接在SAS的IML模塊中調(diào)用R的繪圖腳本，并用ggsave( )函數(shù)保存高分辨率圖形。(4)介紹了疫苗臨床試驗(yàn)中具有特點(diǎn)的統(tǒng)計(jì)圖形，給出了繪制這些統(tǒng)計(jì)圖形的方法。本文通過(guò)連接SAS與R，在SAS中生成并調(diào)用R的繪圖腳本，借助于SAS繪圖宏程序，最終實(shí)現(xiàn)疫苗臨床試驗(yàn)中常用的試驗(yàn)流程圖、逆向累積分布圖、不良事件火山圖的批量自動(dòng)化輸出，大大提高了臨床試驗(yàn)統(tǒng)計(jì)分析的效率。

[1] 張雷,吳瑩,侯艷,等.臨床試驗(yàn)統(tǒng)計(jì)分析報(bào)告SAS與R語(yǔ)言圖形表達(dá)宏程序的實(shí)現(xiàn).中國(guó)衛(wèi)生統(tǒng)計(jì),2013,30(2):294-299.

[2] Hadley Wickham.ggplot2:數(shù)據(jù)分析與圖形藝術(shù).統(tǒng)計(jì)之都譯.第1版.西安：西安交通大學(xué)出版社,2013,1-5.

[3] 金哲植,劉光華.R與SAS的集成在統(tǒng)計(jì)學(xué)教學(xué)中的應(yīng)用研究.科教導(dǎo)刊,2013，(24):103-104.

[4] 才忠勝,唐立波,陳志輝,等.流程圖的自動(dòng)生成系統(tǒng).計(jì)算機(jī)系統(tǒng)應(yīng)用,2013，(12):43-46.

[5] 陳召青,王詩(shī)遠(yuǎn),楊舒靜,等.逆向累積分布圖及其在疫苗免疫原性數(shù)據(jù)分析中的應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2014，31(2):337-339.

[6] 袁延楠,姚晨,閻小妍,等.臨床試驗(yàn)中不良事件數(shù)據(jù)的可視化評(píng)價(jià).中國(guó)新藥雜志,2012，(6):647-665.

[7] 郭貔,李克.聯(lián)合SAS和R分析醫(yī)學(xué)數(shù)據(jù)舉例.中國(guó)衛(wèi)生統(tǒng)計(jì),2012，29(2):288-290.

△通信作者：劉沛

郭海強(qiáng))