南方醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計(jì)學(xué)系(510515) 吳研鵬 錢(qián)晨堅(jiān) 陳平雁 段重陽(yáng)

1.單樣本靈敏度/特異度的檢驗(yàn)

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

第一步：

(1)

(2)

第二步：

(3)

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

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

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

nQuery Advanced 8.6.0.0實(shí)現(xiàn)：設(shè)置單側(cè)檢驗(yàn)水準(zhǔn)α=0.025，檢驗(yàn)效能1-β=80%，其他數(shù)據(jù)相應(yīng)代入，在nQuery Advanced 8.6.0.0 主菜單選擇：

方法框中選擇：One Sensitivity

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

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

圖1 nQuery Advanced 8.6.0.0 關(guān)于例1樣本量估計(jì)的參數(shù)設(shè)置與計(jì)算結(jié)果

SAS 9.4軟件實(shí)現(xiàn)：

%macro sensitivity(alpha= /* 檢驗(yàn)水準(zhǔn) */

，side= /* 單雙側(cè)檢驗(yàn)，1表示單側(cè)檢驗(yàn)，2表示雙側(cè)檢驗(yàn) */

，se0= /* 原假設(shè)下的靈敏度，即

靈敏度目標(biāo)值 */

，se1= /* 備擇假設(shè)下的靈敏度，即預(yù)期的檢驗(yàn)靈敏度 */

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

，power=/* 檢驗(yàn)效能 */

)；

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

%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；

/*定義函數(shù)ca，用于計(jì)算原假設(shè)下X分布的上、下α布的分位數(shù) */

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；

%**從陽(yáng)性樣本量初始值1開(kāi)始迭代 *；

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)；

/*樣本量估計(jì) */

data a；

%**參數(shù)輸入 *；

side=&side；

se0=&se0；

se1=&se1；

alpha=&alpha；

p=&p；

power1=&power；

%**兩步迭代計(jì)算樣本量

第一步：調(diào)用函數(shù)ca，迭代計(jì)算原假設(shè)下X分布的上、下α布的分位數(shù)

第二步：不斷增加陽(yáng)性樣本量，迭代計(jì)算檢驗(yàn)效能，直至滿足設(shè)定條件*；

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)；%**計(jì)算最終總樣本量*；

run；

/*結(jié)果輸出 */

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運(yùn)行結(jié)果：

圖2 SAS 9.4 關(guān)于例1樣本量估計(jì)的參數(shù)設(shè)置和計(jì)算結(jié)果

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

nQuery Advanced 8.6.0.0實(shí)現(xiàn)：設(shè)置單側(cè)檢驗(yàn)水準(zhǔn)α=0.025，檢驗(yàn)效能1-β=80%，其他數(shù)據(jù)相應(yīng)代入，在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：One Specificity

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

圖3 nQuery Advanced 8.6.0.0關(guān)于例2樣本量估計(jì)的參數(shù)設(shè)置與計(jì)算結(jié)果

SAS軟件實(shí)現(xiàn)：

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

SAS運(yùn)行結(jié)果：

圖4 SAS 9.4 關(guān)于例2樣本量估計(jì)的參數(shù)設(shè)置和結(jié)果

2.單樣本AUC檢驗(yàn)

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

(4)

(5)

等級(jí)變量類(lèi)型的資料常見(jiàn)于放射學(xué)研究，例如用CT掃描識(shí)別組織病變，觀測(cè)變量記錄為“正常”，“可能正?！?，“不確定”，“可能不正?！?，“不正常”。Obuchowski和McClish給出了等級(jí)變量下V(θ)估計(jì)公式：

(6)

(7)

(8)

式中，θ表示AUC；f和g分別表示由A和B構(gòu)成的函數(shù)式。A和B構(gòu)造基于如下雙正態(tài)法：

(9)

(10)

等級(jí)型觀測(cè)變量下，B值一般無(wú)法直接得到，可通過(guò)與觀測(cè)變量密切相關(guān)的連續(xù)測(cè)量指標(biāo)，間接估計(jì)陰性組和陽(yáng)性組中對(duì)應(yīng)的潛在連續(xù)型觀測(cè)的標(biāo)準(zhǔn)差比值，也可直接根據(jù)兩組觀測(cè)的離散程度估計(jì)其比值；不確定情況下，通常設(shè)B為1，得到樣本量的保守估計(jì)[5]。

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

(11)

樣本量計(jì)算中，根據(jù)不同觀測(cè)數(shù)據(jù)類(lèi)型，直接將θ=θ0和θ=θ1帶入相應(yīng)V(θ)計(jì)算公式中即可得到V(θ0)與V(θ1)的估計(jì)值。

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

nQuery Advanced 8.6.0.0實(shí)現(xiàn)：設(shè)置雙側(cè)檢驗(yàn)水準(zhǔn)α=0.05(相當(dāng)于單側(cè)0.025)，檢驗(yàn)效能1-β=80%，其他數(shù)據(jù)相應(yīng)代入，在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：One Roc Curve

在彈出的樣本量計(jì)算窗口將各參數(shù)值鍵入，如圖5所示，結(jié)果為422和283，即本試驗(yàn)所需樣本量為陽(yáng)性組422例和陰性組283例，共需705例。

圖5 nQuery Advanced 8.6.0.0關(guān)于例3樣本量估計(jì)的參數(shù)設(shè)置與計(jì)算結(jié)果

SAS 9.4軟件實(shí)現(xiàn)：

%macro oneroc(

alpha=/* 檢驗(yàn)水準(zhǔn)*/

，side=/* 單雙側(cè)檢驗(yàn)，1表示單側(cè)檢驗(yàn)，2表示雙側(cè)檢驗(yàn) */

，type=/* 觀測(cè)變量數(shù)據(jù)類(lèi)型，1表示等級(jí)變量類(lèi)型，2表示連續(xù)變量類(lèi)型 */

，theta0=/* 原假設(shè)下AUC，即AUC目標(biāo)值 */

，theta1=/* 備擇假設(shè)下AUC，即預(yù)期AUC */

，B=/* 等級(jí)變量類(lèi)型中，陰性組和陽(yáng)性組潛在連續(xù)觀測(cè)的標(biāo)準(zhǔn)差比值*/

，R=/* 總樣本量中陰性與陽(yáng)性例數(shù)的比值*/

，power=/* 檢驗(yàn)效能*/

)；

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

%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；

/* 定義函數(shù)V，用于計(jì)算原假設(shè)和備擇假設(shè)下的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)；

/* 樣本量估計(jì) */

data oneroc；

%**參數(shù)輸入 *；

theta0=&theta0；

theta1=&theta1；

side=&side；

alpha=&alpha；

B=&B；

R=&R；

type=&type；

power=&power；

%**樣本量計(jì)算 *；

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；

/* 結(jié)果輸出 */

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運(yùn)行結(jié)果：

圖6 SAS 9.4 關(guān)于例3樣本量估計(jì)的參數(shù)設(shè)置與計(jì)算結(jié)果

3.配對(duì)樣本AUC檢驗(yàn)

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

(12)

(13)

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

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

(14)

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

(15)

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

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

(16)

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

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

(17)

式中r表示θ1和θ2的相關(guān)系數(shù)，利用查表法，根據(jù)r+和r-的平均值及θ1和θ2平均值，可在 Hanley和McNeil提供的表中查到對(duì)應(yīng)r值，不在表中的值可通過(guò)一定合理外推得出[7](r值表見(jiàn)附錄表1)。與等級(jí)類(lèi)型中的雙正態(tài)法估計(jì)方法不同，這里r+和r-可基于陽(yáng)性組和陰性組中的連續(xù)觀測(cè)，使用Pearson相關(guān)系數(shù)或秩相關(guān)系數(shù)直接估計(jì)得出。

表1 r值表*

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

nQuery Advanced 8.6.0.0實(shí)現(xiàn)：設(shè)置檢驗(yàn)水準(zhǔn)α=0.05，雙側(cè)檢驗(yàn)，檢驗(yàn)效能1-β=80%，其他數(shù)據(jù)相應(yīng)帶入。在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：Two Roc Curve

在彈出的樣本量計(jì)算窗口將各參數(shù)值鍵入，如圖7所示，結(jié)果為104和416。即本試驗(yàn)所需樣本量為陽(yáng)性組104例，陰性組416例，共需520例。

圖7 nQuery Advanced 8.6.0.0 關(guān)于例4 樣本量估計(jì)的參數(shù)設(shè)置與計(jì)算結(jié)果

SAS 9.4軟件實(shí)現(xiàn)：

%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=/*根據(jù)r+和r-的平均值、AUC的平均值查表進(jìn)行賦值，不能為空。比如此例為0.72*/

)；

%let error=；

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

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

%let error=1；

%put ERROR：宏參數(shù)(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取值范圍超過(guò)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；

/*定義函數(shù)vv、c，用于計(jì)算AUC的方差、協(xié)方差 */

proc fcmp outlib=work.fuction.conversion；

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

pi=constant('pi')；

%*對(duì)于觀測(cè)變量為連續(xù)型，求AUC的方差；

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

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

%*對(duì)于觀測(cè)變量為等級(jí)型，求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)；

%*對(duì)于觀測(cè)變量為等級(jí)型，求AUC的協(xié)方差；

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；

%*對(duì)于觀測(cè)變量為連續(xù)型，求AUC的協(xié)方差；

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)； %*調(diào)用定義的函數(shù)；

data tworoc；

%**---參數(shù)輸入----**；

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))；%*陽(yáng)性個(gè)體數(shù)；

N_n=ceil(N_p*R)； %*陰性個(gè)體數(shù)；

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；

/*結(jié)果輸出*/

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運(yùn)行結(jié)果：

SAS 9.4陽(yáng)性組樣本量與nQuery Advanced的對(duì)應(yīng)結(jié)果相差3，是因?yàn)閚Query Advanced中在連續(xù)型變量計(jì)算中，對(duì)r的取值與基于查表法提供的r值稍有不同(經(jīng)驗(yàn)證前者r的取值約為0.727，而表格1中r的取值為0.72)。

圖8 SAS 9.4 關(guān)于例4樣本量估計(jì)的參數(shù)設(shè)置與計(jì)算結(jié)果