樣本量估計及其在nQuery+nTerim和SAS軟件上的實現(xiàn)
——均數(shù)比較(八)*

南方醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系 詹志穎 曹穎姝 段重陽 陳平雁

樣本量估計及其在nQuery+nTerim和SAS軟件上的實現(xiàn)
——均數(shù)比較(八)*

南方醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系 詹志穎 曹穎姝 段重陽 陳平雁△

1.2.1.9 兩組泊松分布均數(shù)比較

1.非限制性最大似然估計(MLE，簡記W1)

(1-34)

2.限制性最大似然估計(CMLE，簡記W2)

(1-35)

3.對數(shù)轉(zhuǎn)換的非限制性最大似然估計(ln(MLE)，簡記W3)

(1-36)

4.對數(shù)轉(zhuǎn)換的限制性最大似然估計(ln(CMLE)，簡記W4)

(1-37)

5.方差穩(wěn)健估計(VS，簡記W5)

(1-38)

式中，n1與n2、γ1與γ2以及t1與t2分別為試驗組與對照組的樣本量、單位時間內(nèi)事件發(fā)生率以及個體平均觀察時間；c=R0/R1，其中R0和R1分別為零假設(shè)和備擇假設(shè)下兩組發(fā)生率的比值；ρ=R0t2n2/t1n1。

在計算樣本量時，一般先設(shè)定樣本量初始值，然后迭代樣本量直到所得的檢驗效能滿足條件。此時的樣本量，即研究所需的最小樣本量。

【例1-25】某研究欲探索使用荷爾蒙人群(試驗組)與未使用荷爾蒙人群(對照組)單位時間內(nèi)冠心病發(fā)病率是否相同。要求試驗組與對照組樣本量比例為2：1，觀察時間為2年。根據(jù)以往研究，試驗組冠心病發(fā)病率為0.0005，對照組發(fā)病率是試驗組的4倍。設(shè)定檢驗水準(zhǔn)為0.05，試估計當(dāng)檢驗效能為90%時所需的樣本量[8]。

nQuery+nTerim 4.0實現(xiàn)：設(shè)定檢驗水準(zhǔn)α=0.05，檢驗效能1-β=90%。原假設(shè)為試驗組與對照組發(fā)病率相同，即總體R0=1；備擇假設(shè)為對照組發(fā)病率是試驗組的4倍，即總體R1=4。

在nQuery+nTerim 4.0主菜單選擇：

Goal：⊙Means

Number of Groups：⊙Two

Analysis Method：⊙Test

方法框中選擇：Two Poisson Means

在彈出的樣本量估計窗口，將各參數(shù)值鍵入，選擇統(tǒng)計量，結(jié)果如圖1-61所示，即W1(MLE):試驗組需要8564人，對照組需要4282人，總樣本量為12846人；W2(CMLE):試驗組需要6889人，對照組需要3445人，總樣本量為10334人；W3(ln(MLE))和W4(ln(CMLE)):試驗組需要6685人，對照組需要3343人，總樣本量為10028人；W5(VS)):試驗組需要8590人，對照組需要4295人，總樣本量為12885人。W1與W5方法估計的樣本量偏大，其余三種方法相對來說較接近且偏少。

SAS程序：

OPTIONS LINESIZE=MAX;proc IML;

start PSMS(a,R0,R1,stat,t1,t2,rate1,nr,power);

/*a 檢驗水準(zhǔn);R0、R1分別為H0、H1下兩組發(fā)生率比值;stat檢驗統(tǒng)計量(1:W1,2:W2,3:W3,4:W4,5:W5);t1、t2分別為試驗組、對照組個體平均觀察時間;rate1試驗組發(fā)生率;nr=n2/n1;power檢驗效能*/

error=0;if(a>1|a<0)|(R0<0)|(R1<0)|(stat^=1 & stat^=2 & stat^=3 & stat^=4 & stat^=5)|

(t1<0)|(t2<0)|(rate1<0)|(nr<0)|(power>100|power< 1) then error=1;

圖1-61 nQuery+nTerim 4.0關(guān)于例1-25樣本量估計的參數(shù)設(shè)置與計算結(jié)果

if(error=1) then stop;if(error=0) then do;n1=1;rou=R0*t2*nr/t1;c=R0/R1;

do until(pwer>=power/100);n1=n1+1;n2=ceil(n1*nr);if(stat=1) then do;ts=“W1”;

pwer=probnorm((sqrt(n1*rate1*t1)*abs(1-c))/sqrt(c/rou+c**2)-probit(1-a));end;

else if(stat=2) then do;ts=“W2”;pwer=probnorm((sqrt(n1*rate1*t1)*abs(1-c))/sqrt(c/rou+c**2)-probit(1-a)*sqrt((c+rou)/(1+c*rou)));end;else if(stat=3) then do;ts=“W3”;pwer=probnorm(sqrt(n1*rate1*t1)*abs(log(c))/sqrt(c/rou+1)-probit(1-a));end;

else if(stat=4) then do;ts=“W4”;pwer=probnorm(sqrt(n1*rate1*t1)*abs(log(c))/sqrt(c/rou+1)-probit(1-a)*(sqrt(c*(rou**2+2*rou+1))/(c+rou)));end;else do;ts=“W5”;pwer=probnorm((2*abs(1-sqrt(c))*sqrt(n1*rate1*t1+3/8)-probit(1-a)*sqrt(c/rou+c))/sqrt(c/rou+1));end;end;pwer=round(pwer*100,0.000000000001);print a R0 R1 ts t1 t2 rate1 n1 n2 nr pwer[label=“power(%)”];end;finish PSMS;

run PSMS(0.05,1,4,1,2,2,0.0005,0.5,90);

run PSMS(0.05,1,4,2,2,2,0.0005,0.5,90);

run PSMS(0.05,1,4,3,2,2,0.0005,0.5,90);

run PSMS(0.05,1,4,4,2,2,0.0005,0.5,90);

run PSMS(0.05,1,4,5,2,2,0.0005,0.5,90);quit;

SAS運行結(jié)果見圖1-62。

圖1-62 SAS 9.4關(guān)于例1-25樣本量 估計的參數(shù)設(shè)置與計算結(jié)果

1.2.1.10 兩組負二項分布均數(shù)比較

方法：Zhu和Lakkis(2014)[9]提出的兩組負二項分布均數(shù)差異性檢驗的樣本量估計方法是建立在均值為ln(r1/r2)，方差為σ2=(1/n2)[(1/μt)[1/r2+1/(θr1)]+(1+θ)k/θ]的正態(tài)分布基礎(chǔ)上，其檢驗效能的計算公式為：

(1-39)

(1)Reference Rate(RR):

V0(RR)=(1+θ)/(θμtr2)+(1+θ)k/θ

(1-40)

(2)True Rates(TR):

V0(TR)=(1/μt)[1/r2+1/(θr1)]+(1+θ)k/θ

(1-41)

(3)Maximum Likelihood(ML)

V0(ML)=(1+θ)2/[θμt(r2+θr1)]+(1+θ)k/θ

(1-42)

在計算樣本量時，一般先設(shè)定樣本量初始值，然后迭代樣本量直到所得的檢驗效能滿足條件。此時的樣本量，即研究所需的最小樣本量。

【例1-26】慢性阻塞性肺病的惡化事件研究中，對照組事件率為0.8人年，試驗組較對照組有15%的減少(r1/r2=0.85)，平均暴露時間為0.75年，離散參數(shù)為0.7，且為均衡設(shè)計(θ=1)。設(shè)定檢驗水準(zhǔn)為0.05，試估計當(dāng)檢驗效能為80%時研究所需樣本量[9]。

nQuery+nTerim 4.0實現(xiàn)：

設(shè)定檢驗水準(zhǔn)α=0.05，檢驗效能1-β=80%。原假設(shè)為試驗人群與對照人群的事件率相同，即總體r1=r2；備擇假設(shè)為對照人群事件率與試驗人群不同，即總體r1≠r2。

在nQuery+nTerim 4.0主菜單選擇：

Goal：⊙Means

Number of Groups：⊙Two

Analysis Method：⊙Test

方法框中選擇：Two Negative Binomial Rates

在彈出的樣本量估計窗口，將各參數(shù)值鍵入，結(jié)果如圖3-63所示，即Reference Rate方法，兩組各需1433人，總樣本量為2866人；True Rates方法，兩組各需1494人，總樣本量為2988人；Maximum Likelihood方法，兩組各需1490人，總樣本量為2980人。

圖1-63 nQuery+nTerim 4.0關(guān)于例1-26樣本量估計的參數(shù)設(shè)置與計算結(jié)果

SAS程序：

OPTIONS LINESIZE=MAX;proc IML;

start NBR(a,r2,rr,miut,k,rv,nr,power);

/*a 檢驗水準(zhǔn)； r2 對照組平均事件率；rr=r1/r2；miut個體平均暴露時間；k 離散參數(shù)；rv方差(V0)的計算方法；1：RR；2：TR；3：ML；nr=n1/n2，n1與n2分別為試驗組與對照組的樣本量；power 檢驗效能*/

error=0;if(a>1|a<0)|(r2<0)|(rr<0)|(miut<0)|(k<0)|(rv^=1 &rv^=2 &rv^=3)|(nr<0)|(power>100|power< 1) then error=1;if(error=1) then stop;if(error=0) then do;n2=2;r1=r2*rr;

v1=(1/r2+1/(nr*r1))/miut+(1+nr)*k/nr;if(rv=1) then do;rvlbl=“RR”;

v0=(1+nr)/(nr*miut*r2)+(1+nr)*k/nr;end;else if(rv=2) then do;rvlbl=“TR”;

v0=(1/r2+1/(nr*r1))/miut+(1+nr)*k/nr ;end;else do;rvlbl=“ML”;

v0=(1+nr)**2/(nr*miut*(r2+nr*r1))+(1+nr)*k/nr ;end;do until(pwer>=power);n2=n2+1;pwer=100*probnorm((sqrt(n2)*abs(log(rr))-probit(1-a/2)*sqrt(v0))/sqrt(v1));end;n1=nr*n2;pwer=round(pwer,0.0000001);print a r2 rr miut k rvlbl[lable=“rv”] n1 n2 nr pwer[label=“power(%)”];end;finish NBR;

run NBR(0.05,0.8,0.85,0.75,0.7,1,1,80);

run NBR(0.05,0.8,0.85,0.75,0.7,2,1,80);

run NBR(0.05,0.8,0.85,0.75,0.7,3,1,80);

quit;

SAS運行結(jié)果：

圖1-64 SAS 9.4關(guān)于例1-26樣本量估計的 參數(shù)設(shè)置與計算結(jié)果

[1]呂朵，段重陽，陳平雁.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數(shù)比較(一)，中國衛(wèi)生統(tǒng)計，2012，29(1)：127-131.

[2]鄧居敏，呂朵，陳平雁.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數(shù)比較(二)，中國衛(wèi)生統(tǒng)計，2012，29(1)：132-134.

[3]段重陽，呂朵，陳平雁.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數(shù)比較(三)，中國衛(wèi)生統(tǒng)計，2012，29(2)：275-278.

[4]段重陽，呂朵，陳平雁.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數(shù)比較(四)，中國衛(wèi)生統(tǒng)計，2012，29(2)：279-283.

[5]呂朵，段重陽，陳平雁.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數(shù)比較(五)，中國衛(wèi)生統(tǒng)計，2012，29(3)：451-455.

[6]段重陽，呂朵，陳平雁.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數(shù)比較(六)，中國衛(wèi)生統(tǒng)計，2012，29(3)：456-459.

[7]張斌，呂朵，陳平雁等.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數(shù)比較(七)，中國衛(wèi)生統(tǒng)計，2012，29(4)：598-602.

[8]Gu K,Ng HKT,Tang ML,et al.Testing the ratio of two poisson rates.Biometrical Journal,2008,50(2):283-298.

[9]Zhu H,Lakkis H.Sample size calculation for comparing two negative binomial rates.Statistics in medicine,2014,33(3):376-387.

(責(zé)任編輯：郭海強)

*本研究得到國家自然科學(xué)基金資助(編號：81673270)

編者按：本文為南方醫(yī)科大學(xué)陳平雁教授團隊2012年發(fā)表于本刊的《樣本量估計及其在nQuery+nTerim和SAS軟件上的實現(xiàn)—均數(shù)比較》系列文章[1-7]的后續(xù)部分。前期主要考慮了連續(xù)變量和等級變量的情形，本文將介紹離散變量，即兩組泊松分布及負二項分布均數(shù)比較的樣本量估計方法。文中的公式和實例序號均依照前期的系列文章順序編排，以保持原有結(jié)構(gòu)。

△通信作者：陳平雁