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

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

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

1.2.1.9 兩組泊松分布均數比較

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

(1-34)

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

(1-35)

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

(1-36)

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

(1-37)

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

(1-38)

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

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

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

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

在nQuery+nTerim 4.0主菜單選擇：

Goal：⊙Means

Number of Groups：⊙Two

Analysis Method：⊙Test

方法框中選擇：Two Poisson Means

在彈出的樣本量估計窗口，將各參數值鍵入，選擇統(tǒng)計量，結果如圖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 檢驗水準;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關于例1-25樣本量估計的參數設置與計算結果

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運行結果見圖1-62。

圖1-62 SAS 9.4關于例1-25樣本量 估計的參數設置與計算結果

1.2.1.10 兩組負二項分布均數比較

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

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

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

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

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

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

在nQuery+nTerim 4.0主菜單選擇：

Goal：⊙Means

Number of Groups：⊙Two

Analysis Method：⊙Test

方法框中選擇：Two Negative Binomial Rates

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

圖1-63 nQuery+nTerim 4.0關于例1-26樣本量估計的參數設置與計算結果

SAS程序：

OPTIONS LINESIZE=MAX;proc IML;

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

/*a 檢驗水準； r2 對照組平均事件率；rr=r1/r2；miut個體平均暴露時間；k 離散參數；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運行結果：

圖1-64 SAS 9.4關于例1-26樣本量估計的 參數設置與計算結果

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

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

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

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

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

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

[7]張斌，呂朵，陳平雁等.樣本量估計及其在nQuery和SAS軟件上的實現(xiàn)—均數比較(七)，中國衛(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.

(責任編輯：郭海強)

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

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

△通信作者：陳平雁