定性指標單臂多中心臨床試驗中心效應的調整及SAS程序實現*

北京協和醫學院，中國醫學科學院，國家心血管病中心，阜外醫院，醫學統計部(102300)王闖世 范肖雪 趙延延 白銀曉 王 楊

【提 要】 目的 探討定性指標單臂多中心臨床試驗中調整中心效應的策略及其應用。方法 介紹使用倒方差和D-L法計算中心加權效應值在定性指標單臂多中心臨床試驗中調整中心效應的基本原理和相關考慮，如權重的計算，率的轉換等，編寫SAS程序，并以“人工心臟瓣膜”試驗為實例說明。結果 計算中心加權效應值處理單臂多中心臨床試驗中心效應時，若各中心效應同質，各中心權重等于其效應值方差的倒數；可根據Q值判斷是否需要對各中心權重進行調整。當分中心結局發生率趨近0或1時，建議對原始率進行數據轉換。實例中，使用雙反正弦轉換法的結果顯示，中心加權復合終點事件率的95%置信區間上限低于目標值，拒絕原假設，提示研究產品的療效優于設定的目標療效。結論 定性指標單臂多中心臨床試驗調整中心效應可通過計算中心加權效應值，若分中心結局發生率接近0或1，建議先對原始效應值進行數據轉換。

在多中心臨床試驗(multicenter clinical trial)中，因各中心在受試者基線特征、臨床實踐等方面可能存在差異而產生的中心效應(center effects)[1]是統計分析階段不可忽視的問題，《藥物臨床試驗的生物統計學指導原則》和《ICH E9臨床試驗的統計學指導原則》中均指出需考慮中心效應，選擇恰當的中心效應調整策略是準確評價藥物或醫療器械有效性和安全性的前提。多中心隨機對照試驗(randomized controlled trial，RCT)中，根據主要結局指標的類型(定性或定量)，可采用相應中心效應調整的方法，如協方差分析、Cochran-Mantel-Haenszel(CMH)檢驗和多水平模型(multilevel model)等[2-3]，其中CMH檢驗的本質是計算加權效應值，各中心權重與其樣本量或效應值方差相關。但在單臂多中心臨床試驗中，目前缺少相關方法介紹。本文借鑒RCT中CMH方法調整中心效應的思路，在定性指標的單臂臨床試驗中通過“加權”的方法以調整中心效應，并編寫SAS程序。

材料與方法

1.基本原理

(1)中心加權效應值

CMH方法是主要結局指標為二分類定性變量(如手術是否成功)的多中心RCT中調整中心效應的常用策略之一，其基本原理是把各中心作為“層”，基于各層的四格表數據分別計算效應值并通過一定的算法合并，達到調整中心效應的目的。根據所采用的算法，SAS中提供兩個效應估計值：Mantel-Haenszel(MH)估計值和logit估計值，均為“中心水平”的加權，MH法把各中心樣本量的倒數作為權重[4]，logit法把各中心效應值自然對數的方差取倒數作為權重[5]。類似地，在定性指標的單臂臨床試驗中，可借鑒上述中心加權的策略，通過使用倒方差等算法計算中心加權效應值調整中心效應。

(1)

θi～N(θ，τ2)

(2)

各中心間效應估計值的差異包括隨機抽樣誤差(si)和中心間變異(τ2)；若各中心的效應同質，即各個中心的效應相同或一致，不同中心來自同一總體，則各中心效應估計值的差異為隨機抽樣誤差(si)，沒有中心間變異(τ2=0)。根據倒方差法[6]和D-L法[7]，兩種情形下，中心加權效應值θw及其標準誤Se(θw)的計算分別如公式(4)、(5)和公式(11)、(12)所示。

若各中心效應同質，各中心權重wi取其效應估計值方差的倒數，計算如下：

(3)

(4)

(5)

相應的95%置信區間為：

θw±1.96Se(θw)

(6)

在單臂臨床試驗中，通常會基于主要結局指標預先設定有臨床意義的目標值Δ(objective performance criteria)，在單側0.025的置信水平下，根據公式(6)計算其95%CI是否包含Δ判定是否拒絕原假設：對于高優指標(如手術成功)，若95%CI下限大于Δ，則拒絕原假設；對于低優指標(如死亡)，若95%CI上限小于Δ，則拒絕原假設。亦可通過標準化變換(或z變換)[8]：

(7)

若z>1.96，則在單側0.025的置信水平下拒絕原假設。

各中心效應是否不同(即是否有中心效應)可使用下面公式判斷[7]：

Q=∑wi(θi-θw)2

(8)

(9)

其中，

(10)

(11)

(12)

類似地，把根據公式(11)和(12)得到的加權效應值95%CI與目標值Δ比較，判斷是否拒絕原假設。

(2)效應值(率)的轉換

對于二分類定性結局指標，結局發生數服從二項分布B(n，π)，π為結局發生概率，可以證明當n足夠大，且nπ和nπ(1-π)都大于5時，二項分布B(n，π)近似正態分布N[nπ，nπ(1-π)][8]，前文中心加權效應值的區間估計也是在該前提下完成。但當結局發生率趨近0或1時，可能會出現：(1)按照近似正態估計得到的置信區間下限小于0或者上限大于1；(2)結局發生率趨近于0或者1的中心因方差趨于0而被賦予過大的權重[6]。為解決上述問題，可考慮對各中心效應值(即率)進行數據轉換，計算轉換后的中心加權效應值，再反算。本文主要介紹兩種轉換：logit轉換和雙反正弦轉換(double arcsine transformation)。

logit轉換是醫學研究中常用的一種數據轉換，轉換后的變量θt及其方差var(θt)計算如下：

(13)

(14)

其中，n表示總人數。通過公式(4)～(6)或公式(8)～(12)得到θt的中心加權估計值θtw后，反算未轉換的效應估計值θw：

(15)

需要特別注意的是，當θ等于0或者1時，無法直接進行logit轉換，可以先進行“連續校正”[9]，如對該中心的結局和非結局發生數各加0.5。

雙反正弦轉換[10-11]，又稱為Freeman-Tukey轉換，轉換后的變量θt及其方差var(θt)計算公式如下：

(16)

(17)

其中，y表示結局發生數，n表示總人數。類似地，計算得到θtw后，反算θw的公式如下：

(18)

或使用下面簡易公式反算：

(19)

2.實例資料來源

本文以一項評價置入人工主動脈瓣膜安全性和有效性的多中心、單組目標值臨床試驗為例。該臨床試驗共有5個中心參與，其主要終點指標是術后12個月內的全因死亡或嚴重卒中的發生，評價術后12個月復合事件發生率是否低于目標值(35%)，全分析集FAS(full analysis set)和符合方案集PPS(per protocol set)分別有101例和97例受試者，各中心入組情況見表1。結果將分別展示不考慮中心效應與使用上述權重法調整中心效應的總事件發生率及其95%CI，并提供調整中心效應的SAS程序；同時將展示使用原始率與轉換率計算的加權率及95%CI，以比較上述數據轉換方法。

結 果

1.不考慮中心效應

若不考慮中心效應，對各中心數據進行直接累加估計總體事件發生率及其95%CI，無論是漸進正態法或者精確概率法，FAS和PPS對應的事件發生率的95%CI上限均遠低于目標值35%(表1)，拒絕原假設，可以認為研究產品的療效優于設定的目標療效。

2.加權法調整中心效應

FAS集中，使用原始率，雙反正弦轉換率或logit轉換率均未提示異質性(P>0.1)，使用各中心率的方差倒數作為其相應權重，得到率的加權均值95%CI上限低于目標值，拒絕原假設；考慮中心間變異調整中心權重后結果仍穩定。PPS集中，未調整中心權重計算的加權率95%CI上限也均低于目標值，但雙反正弦轉換和logit轉換法提示存在異質性，調整中心權重后，logit法加權率的95%置信區間上限為35.3%，略高于目標值，但考慮到01中心事件發生數為0，使用logit轉換時進行了“連續校正”，其結果解讀應慎重。

本實例中各中心事件發生數均較少(<5，見表1)，當使用漸進正態估計各中心事件發生率時，95%CI下限“截斷”為0，這種情況下使用原始率方差的倒數作為各中心權重的處理不太合理，提示應對率進行數據轉換。雙反正弦轉換法中各中心權重與中心樣本量成正比，而logit轉換法中權重與中心樣本量和事件發生率均相關，但考慮到存在中心事件發生數為0的情況，logit法需要進行“連續校正”，本實例中更推薦使用雙反正弦轉換法對率轉換后計算加權效應值。

表1 各中心術后12個月復合終點事件發生率

表2 術后12個月復合終點事件率的加權均值(調整中心效應)

3.程序實現

基于第一部分介紹的中心效應調整的策略，我們編寫了相應的SAS宏程序(%S_CG_CenAdj)以便讀者參考。并以上述試驗為例，提供調用該宏程序的代碼，詳見文后附錄。

討 論

以單組目標值設計為代表的單臂臨床試驗在創新醫療器械產品評價中比較常見，在該研究設計下準確評估器械產品的有效性和安全性非常重要，但目前尚缺少單臂臨床試驗中心效應調整策略的研究。借鑒RCT中CMH法調整中心效應的思路，本文介紹了定性指標單臂多中心臨床試驗中調整中心效應的策略及考慮，并以某“人工心臟瓣膜”的多中心單臂試驗為例介紹其應用。

根據ICH E9及中國食品藥品監督管理局發布的《藥物臨床試驗的生物統計學指導原則》，在試驗數據庫鎖定前，統計師需要制定統計分析計劃書并預先指明計劃使用的統計分析方法，尤其是針對主要終點指標的分析；對于多中心臨床試驗，還需明確中心效應調整的策略。目前在單臂多中心臨床試驗中，往往沒有考慮中心效應問題，僅對各中心的結局發生數和總例數進行簡單地累加來計算總體結局發生率。而實踐中各中心受試者招募能力、結局發生率等可能存在較大差異，如本研究實例中，02和04中心結局發生率分別是33.3%(2/6)和5.7%(2/35)，忽視中心效應可能會得到有偏性的效應估計值。本文提出通過計算中心加權效應值以調整中心效應，其思路與多中心RCT中CMH法調整中心效應類似；該策略的倒方差加權法和D-L法與單組率的薈萃分析算法一致，如把各中心看作薈萃分析中每個納入的原始研究，既往亦有研究介紹薈萃分析在二分類定性指標的多中心RCT中處理中心效應的應用[12-13]。

若分中心結局發生(低優指標)或未發生(高優指標)例數較少(如≤5)，建議對效應值(即率)轉換后計算加權效應值；使用logit轉換時，當分中心結局發生率為0或1時，需進行“連續校正”，其結果解讀應慎重。目前關于轉換方法的優劣仍存在爭議，有研究顯示雙反正弦轉換優于logit轉換[6]，但其他研究發現雙反正弦轉換可能會產生令人誤導的結果[14]，主要因為該法在反算率時僅需要一個樣本量[10]，并建議可通過開展一系列敏感性分析(如使用一定范圍的多個樣本量)來檢查結果的穩定性。關于不同轉換方法的比較有待通過模擬研究進一步探討。實際研究中，可開展系列分析(如權重調整前后，率的轉換與否及不同轉換方法等)，提供多種場景下的分析結果以綜合判斷結果的穩定性。

本研究也有其局限性，主要考慮了結局指標為二分類定性指標的單臂試驗，未涉及多分類定性指標或定量指標。實際研究中，二分類定性指標在單臂試驗中更常見，而多分類定性指標和定量指標可根據臨床專業知識轉化為二分類定性指標后開展相關分析。

附 錄

宏程序%S_CG_CenAdj：

/*宏參數使用說明：*/

/*dataset：原始數據集，包含各中心結局事件發生例數及總受試者例數；*/

/*type：取值p，darcsin或logit，表示使用原始率，雙反正弦轉換率或logit轉換率；*/

/*event：各中心結局事件例數；*/

/*total：各中心受試者例數；*/

/*siteno：中心號；*/

/*name：標識率的轉換方法，取值“原始率”、“雙反正弦轉換率”或“logit轉換率”，與type取值相對應，如type取值為p，name取值為“原始率”；*/

/*output：生成新數據集，用于輸出結果*/

%macro S_CG_CenAdj(dataset，type，event，total，siteno，name，output)；

/*數據預處理：包括原始率的轉換及各中心權重(未調整前)計算*/

data weight_&type.；

set &dataset.；

noevent=&total.-&event.；/*未發生事件例數=總例數-事件發生例數*/

total_in=1/&total.；/* 各中心受試者例數的倒數 */

/*原始率*/

%if &type.=p %then%do；

if &event.=0 or &event.=&total.then do；/*當事件發生率為0或1時，*/

event1=&event.+0.5；/*事件數+0.5*/

noevent1=noevent+0.5；/*非事件數+0.5*/

total1=event1+noevent1；/*總例數+1*/

end；

else do；

event1=&event.；

noevent1=noevent；

total1=&total.；

end；

&type._p=&event./&total.；/*各中心事件率，無需上述+0.5校正*/

var=((event1/total1)*(noevent1/total1))/total1；/*各中心事件率的方差*/

%end；

/*率的logit轉換*/

%if &type.=logit %then%do；

if &event.=0 or &event.=&total.then do；/*當事件發生率為0或1時，*/

&event.=&event.+0.5；/*事件數+0.5*/

&total.=&total.+1；/*總例數+1*/

end；

&type._p=log((&event./&total.)/(noevent/&total.))；/*logit轉換：logit(p)=log(p/(1-p))*/

var=(1/(&event.))+(1/(noevent))；/*各中心logit轉換率的方差*/

%end；

/*率的雙反正弦轉換，即Freeman-Tukey轉換*/

%if &type.=darcsin %then%do；

&type._p=arsin(sqrt(&event./(&total.+1)))+arsin(sqrt((&event.+1)/(&total.+1)))；/*FT轉換*/

var=1/(&total.+0.5)；/*各中心FT轉換率的方差*/

%end；

weight=1/var；/*各中心權重*/

p_w=weight*&type._p；/* 權重×事件發生率 */

weight2=weight*weight；/* 權重×權重 */

run；

/*生成宏變量：計算Q統計量、加權率(未調整權重)等*/

proc sql noprint；

select sum(p_w)，/*各中心權重×事件發生率之和*/

sum(weight)，/*各中心權重之和*/

sum(p_w)/sum(weight)，/*率的加權平均(未調整權重)即：加權率*/

1/sqrt(sum(weight))，/*加權率的標準誤*/

sum(weight2)，/*各中心權重×權重之和*/

sum(&total.)，/*總樣本量*/

count(distinct &siteno.)，/*中心數量*/

count(distinct &siteno.)/sum(total_in)/*諧波均數harmonic mean*/

into：p_w_sum，：weight_sum，：effect_p0，：se，：weight2_sum，：total_sum，：C，：h_mean

from weight_&type.；

select sum(weight*(&effect_p0.-&type._p)**2)into：Q /*Q值*/

from weight_&type.；

quit；

/*異質性檢驗，并計算中心間變異(D值)*/

data Q_&type.；

p_Q=1-probchi(&Q.，&C.-1)；/*異質性檢驗*/

if p_Q<0.0001 then pvalue_Q='<0.0001'；else pvalue_Q=put(p_Q，6.4)；

D=max(0，(&Q.-(&C.-1))/(&weight_sum.-(&weight2_sum./&weight_sum.)))；

call symput(′D′，D)；

run；

/*生成宏變量：用于計算加權率(調整權重后)及其95%置信區間的相關指標*/

proc sql noprint；

create table weight_&type.

as select *，1/(&D.+1/weight)as adjust_weight，/*調整后的各中心權重*/

&type._p*calculated adjust_weight as p_adw /*調整后權重×事件發生率*/

from weight_&type.；

select sum(p_adw)，/*各中心調整后權重×事件發生率之和*/

sum(adjust_weight)，/*各中心調整后權重之和*/

sum(p_adw)/sum(adjust_weight)，/*率的加權平均(調整權重后)即：加權率*/

1/sqrt(sum(adjust_weight))/*加權率的標準誤*/

into：p_adw_sum，：adw_sum，：ad_effect_p0，：ad_se

from weight_&type.；

quit；

/*計算事件率及其95%置信區間：原始率，logit轉換率和雙反正弦轉換率；率的反算；調整權重與否*/

data effect_&type.；

set Q_&type.；

Q=put(&Q.，6.2)；

/*原始率*/

% if &type.=p %then%do；

effect_p=&effect_p0.；/*不調整權重effect size*/

lower=max(effect_p-probit(0.975)*&se.，0)；/*95%置信區間下限*/

upper=min(effect_p+probit(0.975)*&se.，1)；/*95%置信區間上限*/

ad_effect_p=&ad_effect_p0.；/*調整權重effect size*/

ad_lower=max(ad_effect_p-probit(0.975)*&ad_se.，0)；/*95%置信區間下限*/

ad_upper=min(ad_effect_p+probit(0.975)*&ad_se.，1)；/*95%置信區間上限*/

%end；

/*logit轉換*/

% if &type.=logit %then%do；

effect_p=exp(&effect_p0.)/(1+exp(&effect_p0.))；/*不調整權重effect size*/

lower=exp(&effect_p0.-probit(0.975)*&se.)/(1+exp(&effect_p0.-probit(0.975)*&se.))；/*95%置信區間下限*/

upper=exp(&effect_p0.+probit(0.975)*&se.)/(1+exp(&effect_p0.+probit(0.975)*&se.))；/*95%置信區間上限*/

ad_effect_p=exp(&ad_effect_p0.)/(1+exp(&ad_effect_p0.))；/*調整權重effect size*/

ad_lower=exp(&ad_effect_p0.-probit(0.975)*&ad_se.)/(1+exp(&ad_effect_p0.-probit(0.975)*&ad_se.))；/*95%置信區間下限*/

ad_upper=exp(&ad_effect_p0.+probit(0.975)*&ad_se.)/(1+exp(&ad_effect_p0.+probit(0.975)*&ad_se.))；/*95%置信區間上限*/

%end；

/*雙重反正弦轉換*/

%if &type.=darcsin %then%do；

effect_p=0.5*(1-sign(cos(&effect_p0.))*sqrt((1-(sin(&effect_p0.)+(sin(&effect_p0.)-1/sin(&effect_p0.))/(&h_mean.))**2)))；/*不調整權重effect size*/

lower0=&effect_p0.-probit(0.975)*&se.；

upper0=&effect_p0.+probit(0.975)*&se.；

lower=0.5*(1-sign(cos(lower0))*sqrt((1-(sin(lower0)+(sin(lower0)-1/sin(lower0))/(&h_mean.))**2)))；/*95%置信區間下限*/

upper=0.5*(1-sign(cos(upper0))*sqrt((1-(sin(upper0)+(sin(upper0)-1/sin(upper0))/(&h_mean.))**2)))；/*95%置信區間上限*/

ad_effect_p=0.5*(1-sign(cos(&ad_effect_p0.))*sqrt((1-(sin(&ad_effect_p0.)+(sin(&ad_effect_p0.)-1/sin(&ad_effect_p0.))/(&h_mean.))**2)))；/*調整權重effect size*/

ad_lower0=&ad_effect_p0.-probit(0.975)*&ad_se.；

ad_upper0=&ad_effect_p0.+probit(0.975)*&ad_se.；

ad_lower=0.5*(1-sign(cos(ad_lower0))*sqrt((1-(sin(ad_lower0)+(sin(ad_lower0)-1/sin(ad_lower0))/(&h_mean.))**2)))；/*95%置信區間下限*/

ad_upper=0.5*(1-sign(cos(ad_upper0))*sqrt((1-(sin(ad_upper0)+(sin(ad_upper0)-1/sin(ad_upper0))/(&h_mean.))**2)))；/*95%置信區間上限*/

/*如果95%置信區間下限或上限小于minimum，則“截斷”為0；如果95%置信區間下限或上限大于maximum，則“截斷”為1*/

minimum=(arsin(sqrt(0/(&h_mean.+1)))+arsin(sqrt((0+1)/(&h_mean.+1))))；

maximum=(arsin(sqrt(&h_mean./(&h_mean.+1)))+arsin(sqrt((&h_mean.+1)/(&h_mean.+1))))；

if lower0

if &effect_p0.

if upper0

if lower0>maximum then lower=1；

if &effect_p0.>maximum then effect_p=1；

if upper0>maximum then upper=1；

if ad_lower0

if &ad_effect_p0.

if ad_upper0

if ad_lower0>maximum then ad_lower=1；

if &ad_effect_p0.>maximum then ad_effect_p=1；

if ad_upper0>maximum then ad_upper=1；

%end；

run；

/*生成輸出數據集，保留關鍵指標*/

data &output.；

length label $50；

retain label effect CI ad_effect ad_CI Q pvalue_Q；

set effect_&type.；

label=“&name.”；

effect=compress(put(effect_p*100，6.2)||'%')；

lower=put(lower*100，6.2)；

upper=put(upper*100，6.2)；

CI=compress(′(′||lower||′%′||′，′||upper||′%′||′)′)；

ad_effect=compress(put(ad_effect_p*100，6.2)||′%′)；

ad_lower=put(ad_lower*100，6.2)；

ad_upper=put(ad_upper*100，6.2)；

ad_CI=compress(′(′||ad_lower||′%′||′，′||ad_upper||′%′||′)′)；

label label=“數據轉換方法” effect=“中心權重法：事件發生率” CI=“中心權重法：95%置信區間” ad_effect=“調整中心權重：事件發生率”

ad_CI=“調整中心權重：95%置信區間” Q=“Q統計量” pvalue_Q=“Q檢驗P值”；

keep label effect CI ad_effect ad_CI Q pvalue_Q；

run；

%mend S_CG_CenAdj；

宏程序調用示例：

%S_CG_CenAdj(example，p，event，total，siteno，原始率，p1)；

%S_CG_CenAdj(example，logit，event，total，siteno，logit轉換率，p2)；

%S_CG_CenAdj(example，darcsin，event，total，siteno，雙反正弦轉換率，p3)；