疫苗臨床試驗中統(tǒng)計圖形的自動化輸出

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

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

疫苗臨床試驗中統(tǒng)計圖形的自動化輸出

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

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

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

R軟件ggplot2包簡介

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

圖形自動化輸出實現(xiàn)方法

1.SAS和R的連接

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

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

2.高分辨率圖形的批量自動化輸出

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

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

疫苗臨床試驗中圖形自動化輸出實例

1.試驗流程圖

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

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

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

圖1 試驗流程圖

2.逆向累積分布圖

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

23價肺炎球菌多糖疫苗由23種血清型肺炎球菌多糖抗原組成，能夠?qū)τ?3種最常見血清型引起的肺炎球菌感染性疾病產(chǎn)生保護(hù)。為了直觀地對試驗組和對照組免前免后每種血清型抗體滴度進(jìn)行比較，通過連接SAS與R，利用SAS逆向累積分布圖宏程序，實現(xiàn)試驗組和對照組免前免后每種血清型抗體滴度的批量自動化輸出。圖2為FAS集血清型1免前抗體逆向累積分布圖，從圖中可以直觀地看出FAS集試驗組和對照組血清型1免前抗體分布情況。實現(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語句調(diào)用R的程序腳本文件：

proc iml;

%include "C:submit.r";

quit;

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

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

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

圖3 不良事件火山圖

(1)用SAS進(jìn)行數(shù)據(jù)處理和統(tǒng)計分析，將統(tǒng)計分析結(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è)置點的顏色 */

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))+′;/* 繪制散點圖并設(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))+′;/* 在點內(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語句調(diào)用R的程序腳本文件：

proc iml;

%include "C:submitae.r";

quit;

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

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

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

討 論

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

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

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

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

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

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

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

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

△通信作者：劉沛

郭海強)