四川大學華西公共衛生學院四川大學華西第四醫院(610041) 張 露 鄭倩雯 李亞文 李宓兒 王 菊 朱彩蓉
【提 要】 目的 探索在SAS和stata中如何實現生物等效性檢驗,并且比較兩個軟件實現方法的異同點和結果的一致性。方法 利用現有加拿大指南中的2×2交叉設計生物等效性試驗的數據集,結合生物等效性檢驗的原理,給出SAS實現的宏程序以及stata程序,并且比較二者結果之間的異同。結果 通過SAS宏語言能夠實現生物等效性檢驗結果的快速輸出。由于stata中“pkequiv”默認的計算原理在最后需要取反對數的過程與SAS有所不同,導致最后結果出現差異,通過添加對最終置信區間取反對數的命令后與SAS結果一致。結論 本研究可作為利用常用軟件進行生物等效性檢驗時參考,同時提示應注意根據自身數據集分布正確使用“pkequiv”命令。
2016年以來,中國高度重視仿制藥和原研藥的一致性評價工作[1],出臺了很多仿制藥與原研藥一致性評價的統計指導原則類文件。然而生物等效性評價作為仿制藥與原研藥一致性評價中不可缺少的一環,目前尚未出臺如何利用軟件實現生物等效性檢驗的具體操作指南。在目前能夠實現生物等效性檢驗的眾多軟件中,國外更推薦在SAS中進行生物等效性檢驗[2]。相對于SAS程序而言,stata在進行生物等效性檢驗時,有專門的命令“pkequiv”,僅需要一個命令即可進行檢驗,十分簡便。故本文將對比stata和SAS進行生物等效性檢驗的異同,為相關工作者利用常用軟件進行生物等效性檢驗提供參考。
生物等效性試驗最常用也是各國指導性文件中最推薦的試驗設計是2×2交叉設計[3],下面將簡單介紹在該試驗設計下的生物等效性檢驗的原理。
生物等效性檢驗最常用的方法之一是置信區間法[4]。需要進行生物等效性檢驗的藥代動力學參數包括峰濃度(cmax)、血藥濃度0到t時間下面積(area under the curve from zero to t time point,AUCT)及血藥濃度0到無窮大時間下面積(area under the curve from zero to infinity,AUCI),這三個參數通常情況下呈正偏態分布,在經過對數轉換后呈正態分布。故首先需要對這三個參數進行對數轉換。進行對數轉換后的參數通過方差分解,獲得殘差均方MSresidual。
利用上述的殘差均方計算估計的標準誤,并計算試驗制劑與參比制劑差值的90%置信區間[5]。
(1)
(2)
根據殘差均方MSresidual計算個體內變異性[6]。
(3)
對置信區間取反對數后,獲得試驗制劑與參比制劑幾何均值比的90%置信區間。在個體內變異性小于30%的情況下,如果這三個參數幾何均值比的90%置信區間在(80%,125%)范圍內,則說明具有生物等效性[7]。
現以加拿大生物等效性試驗指南中的原始數據為例[6],來展示生物等效應檢驗在SAS與stata兩個軟件中的實現。

表1 原始數據格式
seq:給藥順序,1代表TR,2代表RT;period:1代表第一階段,2代表第二階段;trt:藥物,1表示R,2表示T。
1.SAS程序 由于在進行生物等效性檢驗時,三個藥代動力學參數Cmax、AUCT及AUCI在SAS中計算原理完全一致,故我們利用SAS宏和DO循環輸出其生物等效性結果。其中產生殘差的關鍵步驟是“proc mixed”過程。以下為SAS程序:
%macro be(data=);
/*分別建立三個參數的數據集并對參數進行對數轉換*/
data a1;set&data(keep=id trt seq period cmax rename=(cmax=var));logvar=log(var);run;
data a2;set&data(keep=id trt seq period auct rename=(auct=var));logvar=log(var);run;
data a3;set&data(keep=id trt seq period auci rename=(auci=var));logvar=log(var);run;
/*對數據進行混合線性模型分析,獲得方差分析差值估計、殘差*/
%do i = 1 %to 3;
ods output estimates=b&i;odsoutput covparms=c&i;
proc mixed data=a&i;
class trt period seq id;
model logvar=trt period seq/cl alpha=0.1;
random id(seq);
estimate ′T-R′ trt-1 1/cl alpha=0.1;
run;
/*對試驗制劑與參比制劑差值及其90%置信區間取反對數并調整幾何均值比和置信區間的格式*/
data b&I;setb&I;
pe=exp(estimate)*100;lowerCI=exp(lower)*100;upperCI=exp(upper)*100;
CI=cat(′(′,strip(put(lowerCI,8.2)),′%′,′,′,strip(put(upperCI,8.2)),′%′,′)′);
GMR=cat(strip(put(pe,8.2)),′%′);
keep CI GMR;run;
/*計算個體內變異*/
data c&i;set c&i(firstobs=2);
_cv=100*sqrt(exp(estimate)-1);
cv=cat(strip(put(_cv,8.2)),′%′);
keep cv;run;
/*對其中一個參數的各個生物等效性檢驗結果進行橫向合并*/
data d&i;merge b&i c&i;run;
%end;
/*將三個藥代動力學參數的結果進行縱向合并,增加新變量z指示參數名稱*/
data d;setd1 d2 d3;length z $20.;
if _n_=1 then z="Cmax";if _n_=2 then z="AUCT";if _n_=3 then z="AUCI";
run;
%mend;
/*調用宏*/
%be(data=be);
2.stata程序 在stata中,進行生物等效性檢驗有“pkcross”和“pkequiv” 兩個關鍵命令,這兩個命令包含的算法是在stata中預先設定好的。其中“pkcross”將產生方差分析表,獲得殘差均方,后者用來計算個體變異。而“pkequiv”過程則計算比值的90%置信區間。
generate outcome1=log(cmax)
pkcross outcome1,id(id)period(period)sequence(seq)treatment(trt)carryover(none)
pkequiv outcome1 trt period seq id
generate outcome2=log(auct)
pkcross outcome2,id(id)period(period)sequence(seq)treatment(trt)carryover(none)
pkequiv outcome2 trt period seq id
generate outcome3=log(auci)
pkcross outcome3,id(id)period(period)sequence(seq)treatment(trt)carryover(none)
pkequiv outcome3 trt period seq id
SAS生物等效性檢驗與stata的結果對比見表2。

表2 SAS與stata的生物等效性檢驗結果
根據結果對比,stata給出了試驗制劑和參比制劑之間的比值及其置信區間,與SAS給出的幾何均值比和置信區間之間差異巨大,但是其個體內變異卻保持一致。
表3僅以Cmax為例,對比SAS、stata以及指南中給出的結果,SAS的結果與指南基本一致,出現細微差別的原因是指南在進行差值90%置信區間計算時采用的是保留四位小數的估計值。stata的結果中差值及其90%置信區間與SAS結果保持一致,而其比值及其90%置信區間差別明顯。

表3 SAS與stata對Cmax進行生物等效性的檢驗結果對比
通過stata的幫助文獻,“pkequiv”作為一個高度合成的命令,其算法在軟件中是默認的,stata在方差分析、計算差值及其90%置信區間時,其原理與SAS一致,但是SAS在計算幾何均值比及其90%置信區間時是將差值的結果取反對數回去后獲得,而在stata中并沒有計算其幾何均值比,而是通過如下的方法計算試驗制劑和參比制劑之間的比值:
(4)
(5)
(6)

return list
di exp(r(uci))
di exp(r(lci))
通過上述的討論,我們可以看出,默認的stata中的生物等效性檢驗過程實際上是適用于原始藥代動力學參數呈正態分布的數據。盡管stata的原理與SAS在最后一步略有不同而導致結果不一致,但是通過命令的補充,同樣可對藥代動力學參數經對數轉化后呈正態分布的數據進行生物等效性檢驗,且命令更簡單。
在本研究中,SAS利用宏可以多次調用處理同類數據,所有操作都不需要人為的復制粘貼數據;SAS在進行臨床試驗的統計分析時更具有權威性。stata命令相對來說更簡潔,但stata的內置“pkequiv”很容易導致初次使用者由于不了解使用條件而產生誤用。因此本研究更推薦SAS作為生物等效性檢驗的軟件工具。