診斷試驗樣本量估計及其在nQuery和SAS軟件上的實現*

南方醫科大學公共衛生學院生物統計學系(510515) 吳研鵬 錢晨堅 陳平雁 段重陽

1.單樣本靈敏度/特異度的檢驗

方法：Li和Fine(2004)[1]提出了篩檢試驗基于單樣本靈敏度/特異度檢驗的樣本量估計方法，其檢驗效能基于二項分布的確切概率法，采用兩步迭代計算。以靈敏度為例計算公式如下：

第一步：

(1)

(2)

第二步：

(3)

將公式(1)至(3)中的Se0和Se1分別替換為1-Sp0和1-Sp1，N+替換為N-，就是特異度的單樣本樣本量估計公式，其中，Sp0和Sp1分別為原假設和備擇假設下的特異度，N-為陰性組樣本量，最終總樣本量為n=N-/(1-p)。

需要特別指出，上述公式用于全人群篩檢時p是人群患病率，如果用于醫院就診人群的診斷，p應該代表就診人群的陽性病例的比例，即試驗樣本是隨機就診人群而不是選擇人群，否則所估計的樣本量會偏高或偏低，其中若選擇人群p偏高，則所估計樣本量偏高，反之偏低[2-4]。

【例1】 某研究旨在評價高危型人乳頭瘤病毒(HR-HPV-DNA)檢測較超薄液基細胞學(TCT)檢測在宮頸癌篩查中的診斷價值，篩查對象為35～64歲已婚婦女，以病理活檢結果為金標準。根據既往研究結果，該人群的宮頸癌患病率為3.32/萬。預期HR-HPV-DNA檢測的靈敏度為0.850，以既往研究報告TCT診斷宮頸癌的靈敏度0.765為目標值(即原假設Se0=0.765)。設定單側檢驗水準α=0.025，檢驗效能為80%，試估計HR-HPV-DNA檢測技術的靈敏度不低于0.765所需樣本量。

nQuery Advanced 8.6.0.0實現：設置單側檢驗水準α=0.025，檢驗效能1-β=80%，其他數據相應代入，在nQuery Advanced 8.6.0.0 主菜單選擇：

方法框中選擇：One Sensitivity

在彈出的樣本量計算窗口將各參數值鍵入，如圖1第1列數據所示，結果為n=502857，即本試驗的最少樣本量為502857例。

若該檢測方法應用于醫院內就診人群，且假設醫院既往接受該檢測方法的患者中病例活檢陽性的比例為20%，如圖1第2列所示，則HR-HPV-DNA檢測技術不低于0.765所需樣本量為880例。由此可見，基于篩查和診斷的目的不同，樣本量有很大差別，能否準確估計樣本量，關鍵之一在于篩查人群患病率或就診人群陽性率的正確估計。

圖1 nQuery Advanced 8.6.0.0 關于例1樣本量估計的參數設置與計算結果

SAS 9.4軟件實現：

%macro sensitivity(alpha= /* 檢驗水準 */

，side= /* 單雙側檢驗，1表示單側檢驗，2表示雙側檢驗 */

，se0= /* 原假設下的靈敏度，即

靈敏度目標值 */

，se1= /* 備擇假設下的靈敏度，即預期的檢驗靈敏度 */

，p= /* 陽性率(人群患病率)*/

，power=/* 檢驗效能 */

)；

/*判斷輸入的參數是否符合要求，并輸出日志提示 */

%if(&alpha>1 or &alpha<0)%then %do；

%put ERR%STR(OR)：Test significance level must be in 0-1；

%goto exit；

%end；

%if(&side^=1 and &side^=2)%then %do；

%put ERR%STR(OR)：s=1 or 2；

%goto exit；

%end；

%if(&p>1 or & p<0)% then %do；

%put ERR%STR(OR)：Overall Event Rate(P)must be in 0-1；

%goto exit；

%end；

%if(&power>100 or & power<0)% then %do；

%put ERR%STR(OR)：Power must be in 0-100；

%goto exit；

%end；

/*定義函數ca，用于計算原假設下X分布的上、下α布的分位數 */

proc datasets nolist nodetails NoWarn lib=work；

delete funcsol/memtype=data；

quit；

proc fcmp outlib=work.funcsol.conversion；

function ca(q，n，s)；

n1=1；

x=0；

%**從陽性樣本量初始值1開始迭代 *；

do until(x gt q)；

x=cdf('binomial'，n1，s，n)；

n1=n1+1；

end；

n1=n1-1；

return(n1)；

endsub；

run；

options cmplib=(work.funcsol)；

/*樣本量估計 */

data a；

%**參數輸入 *；

side=&side；

se0=&se0；

se1=&se1；

alpha=&alpha；

p=&p；

power1=&power；

%**兩步迭代計算樣本量

第一步：調用函數ca，迭代計算原假設下X分布的上、下α布的分位數

第二步：不斷增加陽性樣本量，迭代計算檢驗效能，直至滿足設定條件*；

n1=2；

power=0；

do until(power>=&power)；

x1=ca(1-alpha/side，n1，se0)；

x2=ca(alpha/side，n1，se0)；

if side=2 then power=(1-probbnml(se1，n1，x1)+probbnml(se1，n1，x2-1))*100；

if side=1 then do；

if se0

else power=(probbnml(se1，n1，x2-1))*100；

end；

n1=n1+1；

end；

n1=n1-1；

power=round(power，0.01)；

n=floor(n1/p)；%**計算最終總樣本量*；

run；

/*結果輸出 */

proc print data=a label；

var alpha side se0 se1 p n power；

label alpha="Test Significance Level"

side="1 or 2 sided"

se0="Null Hypothesis Sensitivity"

se1="Alternative Hypothesis Sensitivity"

p='Prevalence'

n="Sample Size"

power="Power(%)"；

quit；

%exit：

%mend sensitivity；

%sensitivity(alpha=0.025，side=1，se0=0.765，se1=0.850，p=0.00035，power=80)

%sensitivity(alpha=0.025，side=1，se0=0.765，se1=0.850，p=0.20，power=80)

SAS運行結果：

圖2 SAS 9.4 關于例1樣本量估計的參數設置和計算結果

【例2】以例1為例，預期HR-HPV-DNA合并TCT檢測法診斷的特異度為0.807，以既往研究報告TCT診斷宮頸癌的特異度0.732為目標值(即原假設Sp0=0.732)，設定單側檢驗水準α=0.025，檢驗效能為80%下，試估計HR-HPV-DNA合并TCT檢測技術的特異度不低于0.732所需樣本量。

nQuery Advanced 8.6.0.0實現：設置單側檢驗水準α=0.025，檢驗效能1-β=80%，其他數據相應代入，在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：One Specificity

彈出的樣本量計算窗口將各參數值鍵入，如圖3第1列數據所示，結果為n=251。即本試驗的最少樣本量為251例。另外，考慮該檢測技術在醫院就診人群的應用，且假設醫院既往接受該檢測方法的患者中病例活檢陽性的比例為20%，如圖3第2列數據所示，則HR-HPV-DNA檢測技術特異度不低于0.732所需樣本量為313例。

圖3 nQuery Advanced 8.6.0.0關于例2樣本量估計的參數設置與計算結果

SAS軟件實現：

在單樣本靈敏度的SAS樣本量估算程序中，將輸入參數Se0和Se1直接替換為Sp0和Sp1，迭代公式中，Se0和Se1分別替換為1-Sp0和1-Sp1，最后計算出總樣量n=N-/(1-p)。

SAS運行結果：

圖4 SAS 9.4 關于例2樣本量估計的參數設置和結果

2.單樣本AUC檢驗

方法：Obuchowski和McClish[5]以及Hanley和McNeil[6-7]給出了觀測變量為等級類型和連續類型下的單樣本AUC(ROC曲線下的面積)檢驗的樣本量和檢驗效能的估計方法，其計算公式如下：

(4)

(5)

等級變量類型的資料常見于放射學研究，例如用CT掃描識別組織病變，觀測變量記錄為“正常”，“可能正常”，“不確定”，“可能不正常”，“不正常”。Obuchowski和McClish給出了等級變量下V(θ)估計公式：

(6)

(7)

(8)

式中，θ表示AUC；f和g分別表示由A和B構成的函數式。A和B構造基于如下雙正態法：

(9)

(10)

等級型觀測變量下，B值一般無法直接得到，可通過與觀測變量密切相關的連續測量指標，間接估計陰性組和陽性組中對應的潛在連續型觀測的標準差比值，也可直接根據兩組觀測的離散程度估計其比值；不確定情況下，通常設B為1，得到樣本量的保守估計[5]。

連續型觀測變量下，可基于AUC值直接估計方差，Hanley和McNeil給出了連續型變量下的V(θ)估計公式：

(11)

樣本量計算中，根據不同觀測數據類型，直接將θ=θ0和θ=θ1帶入相應V(θ)計算公式中即可得到V(θ0)與V(θ1)的估計值。

【例3】 某研究欲評估磁共振成像(MRI)在診斷疑似膝關節損傷患者是否發生關節軟骨損傷的臨床價值，關節鏡檢查作為金標準，影像科醫生根據每個病人MRI顯像做如下評分：1=肯定沒有軟骨異常；2=可能沒有軟骨異常；3=懷疑軟骨異常；4=可能有軟骨異常；5=肯定軟骨異常。研究顯示疑似膝關節損傷患者中60%存在軟骨異常，因此正常組與異常組樣本量比值設為0.67(即0.4/0.6)，假設正常組和異常組的觀測標準差相同，即B=1。預期MRI診斷AUC為0.850，AUC的目標值為0.80，設定雙側檢驗水準α=0.05，檢驗效能為80%，試估計MRI診斷的AUC高于目標值所需樣本量。

nQuery Advanced 8.6.0.0實現：設置雙側檢驗水準α=0.05(相當于單側0.025)，檢驗效能1-β=80%，其他數據相應代入，在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：One Roc Curve

在彈出的樣本量計算窗口將各參數值鍵入，如圖5所示，結果為422和283，即本試驗所需樣本量為陽性組422例和陰性組283例，共需705例。

圖5 nQuery Advanced 8.6.0.0關于例3樣本量估計的參數設置與計算結果

SAS 9.4軟件實現：

%macro oneroc(

alpha=/* 檢驗水準*/

，side=/* 單雙側檢驗，1表示單側檢驗，2表示雙側檢驗 */

，type=/* 觀測變量數據類型，1表示等級變量類型，2表示連續變量類型 */

，theta0=/* 原假設下AUC，即AUC目標值 */

，theta1=/* 備擇假設下AUC，即預期AUC */

，B=/* 等級變量類型中，陰性組和陽性組潛在連續觀測的標準差比值*/

，R=/* 總樣本量中陰性與陽性例數的比值*/

，power=/* 檢驗效能*/

)；

/* 判斷輸入的參數是否符合要求，并輸出日志提示*/

%if(&alpha>1 or &alpha<0)%then%do；

%put ERR%STR(OR)：Test significance level must be in 0-1；

%goto exit；

%end；

%if(&side^=1 and &side^=2)%then%do；

%put ERR%STR(OR)：s=1 or 2；

%goto exit；

%end；

%if(&type^=1 and &type^=2)%then%do；

%put ERR%STR(OR)：type=1 or 2；

%goto exit；

%end；

%if(&power>100 or &power<0)%then%do；

%put ERR%STR(OR)：Power must be in 0-100；

%goto exit；

%end；

/* 定義函數V，用于計算原假設和備擇假設下的AUC方差*/

proc datasets nolist nodetails NoWarn lib=work；

delete funcsol/memtype=data；

quit；

proc fcmp outlib=work.funcsol.conversion；

function V(theta，B，R，type)；

pi=constant(′pi′)；

if type=2 then do；

v=theta/(R*(2-theta))+2*theta*theta/(1+theta)-theta*theta*(1+R)/R；

end；

if type=1 then do；

A=PROBIT(theta)*sqrt(1+B**2)；

E=exp(-A*A/(2+2*B**2))；

g=A*B*E/(sqrt(2*pi*((1+B**2)**3)))；

f=E/sqrt(2*pi*(1+B**2))；

v=f**2*(1+B**2/R+A**2/2)+g**2*(B**2)*(1+R)/(2*R)；

end；

return(v)；

endsub；

run；

options cmplib=(work.funcsol)；

/* 樣本量估計 */

data oneroc；

%**參數輸入 *；

theta0=&theta0；

theta1=&theta1；

side=&side；

alpha=&alpha；

B=&B；

R=&R；

type=&type；

power=&power；

%**樣本量計算 *；

N_p1=(probit(1-alpha/side)*(sqrt(V(theta0，B，R，type)))+probit(power/100)*(sqrt(V(theta1，B，R，type))))**2/((theta1-theta0)**2)；

N_p=ceil(N_p1)；

N_n=N_p*R；

N_n=ceil(N_n)；

power=probnorm((sqrt(N_p)*ABS(theta1-theta0)-probit(1-alpha/side)*

sqrt(V(theta0，B，R，type)))/sqrt(V(theta1，B，R，type)))*100；power=round(power，0.01)；

run；

data oneroc1；

set oneroc；

if type=2 then q=′Continuous′；

if type=1 then q=′Discrete′；

run；

/* 結果輸出 */

proc print data=oneroc1 label；

var alpha side q theta0 theta1 B R N_p N_n power；

label alpha=“Test significance level”

side=“1 or 2 sided test”

q=“Data Type”

theta0=“Null Hypothesis AUC”

theta1=“Alternative Hypothesis AUC”

B=“Standard Deviation Ratio”

R=“Sample Size Ratio”

N_p=“Positive Group Sample Size”

N_n=“Negative Group Sample Size”

power=“Power(%)”；

quit；

%exit：

%mend oneroc；

%oneroc(alpha=0.05，side=2，type=1，theta0=0.800，theta1=0.850，B=1，R=0.67，power=80)

SAS運行結果：

圖6 SAS 9.4 關于例3樣本量估計的參數設置與計算結果

3.配對樣本AUC檢驗

方法：與單樣本AUC檢驗不同，將目標AUC值θ0替代成同人群下的另一種診斷AUC值θ2，即進行配對樣本AUC檢驗，其樣本量和檢驗效能的估計方法與單樣本AUC檢驗同時提出[5，7]，其計算公式如下：

(12)

(13)

式中，Δ表示兩組診斷下的AUC差值，即θ2-θ1；V0(Δ)和V1(Δ)分別表示Δ在原假設和備擇假設下的方差函數，其計算公式如下：

V0(Δ)=2V(θ1)-2C(θ1，θ1)

(14)

V1(Δ)=V(θ1)+V(θ2)-2C(θ1，θ2)

(15)

V(θi)的意義和計算方法同上述單樣本AUC檢驗，計算方法見公式(6)至(11)。C(θi，θj)表示θi和θj協方差或近似協方差，其計算方法也取決于觀測變量是等級型還是連續型。

對于等級變量，C(θi，θj)的估計公式：

(16)

式中，r-和r+分別表示雙正態法下(見公式(9))，陰性組中兩診斷觀測的相關系數和陽性組中兩診斷觀測的相關系數，等級型觀測下通常不能直接得出，計算樣本量時可參考兩組等級觀測間秩相關系數，考慮一系列合理取值。下標i，j標注的參數分別是診斷試驗i，j中對應的參數，上述公式中其余參數意義和計算方法與在單樣本AUC檢驗中一致，這里不再贅述。

對于連續型變量，Hanley和McNeil給出了V(θi)(見公式(11))和C(θi，θj)的估計公式：

(17)

式中r表示θ1和θ2的相關系數，利用查表法，根據r+和r-的平均值及θ1和θ2平均值，可在 Hanley和McNeil提供的表中查到對應r值，不在表中的值可通過一定合理外推得出[7](r值表見附錄表1)。與等級類型中的雙正態法估計方法不同，這里r+和r-可基于陽性組和陰性組中的連續觀測，使用Pearson相關系數或秩相關系數直接估計得出。

表1 r值表*

【例4】某研究欲比較血氧飽和度下降指數(ODI)和呼吸紊亂指數(RDI)對是否患阻塞性睡眠呼吸暫停綜合征(OSA)的診斷價值。對同一組病人的OSA診斷預試驗中，ODI指數診斷法的AUC為0.85，RDI指數診斷法的AUC為0.90；兩診斷指數在OSA陽性組中的相關系數r+為0.80，陰性組中的相關系數r-為0.76；兩診斷指數的陰性組與陽性組標準差比值B1和B2分別為1。試估計雙側檢驗水準α=0.05，檢驗效能為80%，陰性組和陽性組的樣本比例R=4，能夠發現兩種指數診斷的AUC差異所需的樣本量。

nQuery Advanced 8.6.0.0實現：設置檢驗水準α=0.05，雙側檢驗，檢驗效能1-β=80%，其他數據相應帶入。在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：Two Roc Curve

在彈出的樣本量計算窗口將各參數值鍵入，如圖7所示，結果為104和416。即本試驗所需樣本量為陽性組104例，陰性組416例，共需520例。

圖7 nQuery Advanced 8.6.0.0 關于例4 樣本量估計的參數設置與計算結果

SAS 9.4軟件實現：

%macrotworoc(

alpha= /*Test significance level*/

，side= /*1：one sided test，2：two sided test*/

，type= /*1：Discrete，2：Continuous*/

，theta1= /*Reference Test AUC*/

，theta2= /*New Test AUC*/

，r1= /*Correlation for Positive Group*/

，r2= /*Correlation for Negative Group*/

，B1= /*Test 1 St.Dev.Ratio*/

，B2= /*Test 2 St.Dev.Ratio*/

，R= /*Sample Size Ratio*/

，power=/*Power(%)*/

，rvalue=/*根據r+和r-的平均值、AUC的平均值查表進行賦值，不能為空。比如此例為0.72*/

)；

%let error=；

/*判斷輸入的參數是否符合要求，并輸出日志提示 */

%if %length(&rvalue.)=0 %then %do；

%let error=1；

%put ERROR：宏參數(rvalue)不能為空。；

%end；

%if(&alpha.>1 | &alpha.<0)%then %do；

%let error=1；

%put %str(ERROR：Test significance level%'s range：0-1)；

%end；

%if(&side.^=1 & &side.^=2)%then %do；

%let error=1；

%put %str(ERROR：Side%'s range：1 OR 2)；

%end；

%if(&type.^=1 & &type.^=2)%then %do；

%let error=1；

%put %str(ERROR：type%'s range：1 OR 2)；

%end；

%if(%sysevalf(&theta1.+&theta2.)<1.4 | %sysevalf(&theta1.+&theta2.)>1.95)%then %do；

%let error=1；

%put %str(ERROR：sum%(theta1，theta2%)%'s range：1.4<=P=<1.95)；

%end；

%*r取值范圍超過r值表；

%if(%sysevalf(&r1.+&r2.)<0.04 | %sysevalf(&r1.+&r2.)>1.8)%then %do；

%let error=1；

%put %str(ERROR：sum(r1，r2)%'s range：0.04

%end；

%if(&power.>100 | &power.<0)%then %do；

%let error=1；

%put %str(ERROR：Power%'s range：0-100)；

%end；

%if &error.=1 %then %goto exit；

proc datasets nolist nodetails NoWarn lib=work；

delete fuction/memtype=data；

quit；

/*定義函數vv、c，用于計算AUC的方差、協方差 */

proc fcmp outlib=work.fuction.conversion；

function vv(theta，B，R，type)；

pi=constant('pi')；

%*對于觀測變量為連續型，求AUC的方差；

if type=2 then v=theta/(R*(2-theta))

+2*theta*theta/(1+theta)- theta*theta*(1+R)/R；

%*對于觀測變量為等級型，求AUC的方差；

if type=1 then do；

A=probit(theta)*sqrt(1+B**2)；

E=exp(-A*A/(2+2* B**2))；

g=-A*B*E/sqrt(2*pi*((1+B**2)**3))；

f=E/sqrt(2*pi*(1+B**2))；

v=f**2*(1+B**2/R+A**2/2)+g**2*(B**2)*(1+R)/(2*R)；

end；

return(v)；

endsub；

run；

proc fcmp outlib=work.fuction.conversion；

function C(theta1，theta2，B1，B2，r1，r2，R，type)；

A1=probit(theta1)* sqrt(1+B1*B1)；

E1=exp(-A1*A1/(2+2*B1*B1))；

g1=-A1*B1*E1/sqrt(2*constant('pi')*((1+B1**2)**3))；

f1=E1/sqrt(2*constant('pi')*(1+B1**2))；

v1=f1**2*(1+B1**2/R+A1**2/2)+g1**2*(B1**2)*(1+R)/(2*R)；

A2=probit(theta2)* sqrt(1+B2*B2)；

E2=exp(-A2*A2/(2+2*B2*B2))；

g2=-A2*B2*E2/sqrt(2*constant('pi')*((1+B2**2)**3))；

f2=E2/sqrt(2*constant('pi')*(1+B2**2))；

v2 =f2**2*(1+B2**2/R+A2**2/2)+g2**2*(B2**2)*(1+R)/(2*R)；

%*對于觀測變量為等級型，求AUC的協方差；

if type=1 then c=f1*f2*(r1+B1*B2*r2/R+A1*A2*r1*r1/2)

+g1*g2*B1*B2*(r2**2 +R*(r1**2))/2/R

+f1*g2*A1*B2*r1*r1/2

+f2*g1*A2*B1*r1*r1/2；

%*對于觀測變量為連續型，求AUC的協方差；

if type=2 then do；

v1=theta1/(R*(2-theta1))

+2*theta1*theta1/(1+theta1)- theta1*theta1*(1+R)/R；

v2=theta2/(R*(2-theta2))+2*theta2*theta2/(1+theta2)

- theta2*theta2*(1+R)/R；

/*&rvalue為查表所得，比如此例為rvalue=0.72 */

c=&rvalue.*sqrt(v1*v2)；

end；

return(c)；

endsub；

run；

optionscmplib=(work.fuction)； %*調用定義的函數；

data tworoc；

%**---參數輸入----**；

zero1=2*vv(&theta1，&B1，&R，&type)； %*求得V0(Δ)中的：2V(θ1)；

one1=vv(&theta1，&B1，&R，&type)

+vv(&theta2，&B2，&R，&type)；%*求得V1(Δ)中的：V(θ1)+V(θ2)；

%*求得V0(Δ)中的：2C(θ1，θ1)；

zero2=2*C(&theta1，&theta1，&B1，&B1，&r1，&r2，&R，&type)；

%*求得V1(Δ)中的：2C(θ1，θ2)；

one2=2*C(&theta1，&theta2，&B1，&B2，&r1，&r2，&R，&type)；

alpha=&alpha；

side=&side；

theta1=&theta1；

theta2=&theta2；

r1=&r1；

r2=&r2；

B1=&B1；

B2=&B2；

R=&R；

type=&type；

power=&power；

power=power/100；

V0=zero1-zero2； %*求得V0(Δ)；

V1=one1-one2； %*求得V1(Δ)；

N_p=ceil((probit(1-alpha/side)*sqrt(V0)+probit(power)*sqrt(V1))**2

/((&theta1-&theta2)**2))；%*陽性個體數；

N_n=ceil(N_p*R)； %*陰性個體數；

R=N_n/ N_p；

power=round(100*probnorm((sqrt(N_p)*ABS(&theta1-&theta2)

-probit(1-alpha/side)*sqrt(V0))/sqrt(V1))，0.01)；

run；

data tworoc1；

set tworoc；

if type=2 then q='Continuous'；

if type=1 then q='Discrete'；

run；

/*結果輸出*/

proc print data=tworoc1 label；

var alpha side q theta1 theta2 r1 r2 B1 B2 R N_p N_n power；

label alpha="Test significance level"

side="1 or 2 sided test"

q="Data Type"

theta1="Reference Test AUC"

theta2="New Test AUC"

r1="Correlation for Positive Group"

r2="Correlation for Negative Group"

B1="Test 1St.Dev.Ratio"

B2="Test 2St.Dev.Ratio"

R="Sample Size Ratio"

N_p="Positive Group Sample Size"

N_n="Negative Group Sample Size"

power="Power(%)"；

quit；

%exit：

%mend；

%tworoc(alpha=0.05，side=2，type=2，theta1=0.85，theta2=0.90，

r1=0.80，r2=0.76，B1=1，B2=1，R=4，power=80，rvalue=0.72)

SAS運行結果：

SAS 9.4陽性組樣本量與nQuery Advanced的對應結果相差3，是因為nQuery Advanced中在連續型變量計算中，對r的取值與基于查表法提供的r值稍有不同(經驗證前者r的取值約為0.727，而表格1中r的取值為0.72)。

圖8 SAS 9.4 關于例4樣本量估計的參數設置與計算結果