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

張 斌 呂 朵 段重陽 陳平雁△

1.3 多樣本的均數(shù)比較

1.3.1 差異性檢驗

1.3.1.1 One-way ANOVA

方法:O'Brien和 Muller(1993)〔2〕給出的 One-way ANOVA樣本量估計是建立在自由度為G-1，N-G，非中心參數(shù)為N·V/σ2的非中心F分布上。其檢驗效能的計算公式為:

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

【例1－22】一項有關(guān)降血壓藥的臨床試驗，設(shè)置3個處理組，即安慰劑、陽性對照藥和新藥，以舒張壓下降值為主要療效評價指標(biāo)。由以往研究結(jié)果獲知，安慰劑可使舒張壓平均下降5mmHg，陽性對照藥可下降12mmHg，公共標(biāo)準(zhǔn)差為6mmHg。我們預(yù)期新藥的效果與陽性對照藥相當(dāng)，即可使舒張壓下降12mmHg。若設(shè)定檢驗效能為90%，試估計樣本量。

nQuery Advisor 7.0實現(xiàn):設(shè)定檢驗水準(zhǔn) α=0.05;檢驗效能取1－β=90% 。依據(jù)上述基礎(chǔ)數(shù)據(jù)可知，μ1=5，μ2=12，μ3=12，σ =6。在 nQuery Advisor 7.0主菜單選擇:

Goal:Make Conclusion Using:⊙Means

Number of Groups:⊙ ＞ Two

Analysis Method:⊙Test

方法框中選擇:One-way analysis of variance。

Assistants:⊙Compute Effect Size

在彈出的標(biāo)準(zhǔn)差計算窗口將各參數(shù)鍵入，如圖1－52所示，結(jié)果為V=10.899。

圖1－52 nQuery Advisor 7.0關(guān)于例1-22樣本量估計的參數(shù)計算結(jié)果

將計算結(jié)果V和其他參數(shù)鍵入樣本量計算窗口，如圖1－53所示，結(jié)果為n=15。

圖1－53 nQuery Advisor 7.0關(guān)于例1-22樣本量估計的參數(shù)設(shè)置與計算結(jié)果

SAS9.2軟件實現(xiàn):

%let u={5 12 12};

%let r={1 1 1};PROC IML;

start MGT0(a，G，sd，power);error=0;

if(a＞1|a＜0)then do;error=1;print“error”“Test significance level must be in 0-1”;end;

if(sd ＜0)then do;error=1;print“error”“standard deviation must be ＞ =0”;end;

if(G＜0|ceil(G)^=G)then do;error=1;print“error”“The Number of groups must be positive integer”;end;

if(power＞100|power＜1)then do;error=1;print“error”“Power(%)must be in 1-100”;end;

if(error=1)then stop;

if(error=0)then do;

V=sum((＆u-sum(＆u#＆r/sum(＆r)))##2#＆r/sum(＆r));es=V/sd##2;total=G+1;

do until(pw＞=power/100);ncp=total*es;df1=G-1;df2=total-G;f=FINV(1-a，df1，df2);

pw=1-PROBF(f，df1，df2，ncp);total=total+0.01;end;

total=ceil(total-0.01);n=ceil(total/G);

print a［label=“Test significance level”］

G［label=“Number of groups”］

V［label=“Variance of means”］

sd［label=“Common standard deviation”］

es［label=“Effect size”］

power［label=“Power(%)”］

n［label=“n per group”］;end;

finish MGT0;

run MGT0(0.05，3，6，90);quit;

SAS運行結(jié)果:

圖1－54 SAS9.2關(guān)于例1-22樣本量估計的參數(shù)設(shè)置與計算結(jié)果

式中，C代表每組的樣本均數(shù)和對應(yīng)對比系數(shù)的乘積的和，即∑μici;D=∑c2i;G為組數(shù);σ為樣本標(biāo)準(zhǔn)差。

當(dāng)各組的樣本量相等時，檢驗效能的計算公式與樣本量不相等的情況一樣，但每一組的對比系數(shù)應(yīng)都為1。

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

【例1－23】接例22實驗設(shè)計將安慰劑組、低劑量新藥組、高劑量新藥組與標(biāo)準(zhǔn)降血壓藥進行比較。我們估計安慰劑組的血壓會有5mmHg的下降，標(biāo)準(zhǔn)藥組下降 12mmHg，低劑量和高劑量分別下降10.5mmHg和13.5mmHg。我們再估計血壓的標(biāo)準(zhǔn)差為6mmHg。為了保證足夠的樣本量，我們假設(shè)每組的標(biāo)準(zhǔn)差都為6mmHg。本研究若在90%的檢驗效能條件下，試估計樣本量。

仍然是有關(guān)降血壓藥的臨床試驗，設(shè)置4個處理組，即安慰劑、陽性對照藥、低劑量新藥組和高劑量新藥組，以舒張壓下降值為主要療效評價指標(biāo)。由以往研究結(jié)果獲知，安慰劑可使舒張壓平均下降5mmHg，陽性對照藥可下降12mmHg。由預(yù)試驗得到的數(shù)據(jù)顯示，低劑量新藥和高劑量新藥分別使舒張壓下降10.5mmHg和 13.5mmHg。假定公共標(biāo)準(zhǔn)差為6mmHg，若設(shè)定檢驗效能為90%，試估計樣本量。

nQuery Advisor 7.0實現(xiàn):設(shè)定檢驗水準(zhǔn) α=0.05;雙側(cè)檢驗，即 s=2;檢驗效能取1－β=90% 。依據(jù)上述基礎(chǔ)數(shù)據(jù)估計，μ1=5，μ2=10.5，μ3=13.5，μ4=12，σ=6，。在nQuery Advisor 7.0主菜單選擇:

Goal:Make Conclusion Using:⊙Means

Number of Groups:⊙ ＞ Two

Analysis Method:⊙Test

方法框中選擇:Single One-way between means contrast。

Assistants:⊙Compute Effect Size

在彈出的標(biāo)準(zhǔn)差計算窗口將各參數(shù)鍵入，如圖1－55所示，結(jié)果為C=3，D=1.414。將計算結(jié)果傳輸至主對話框，如圖1－56所示，結(jié)果為n=85。

圖1－55 nQuery Advisor 7.0關(guān)于例1-23樣本量估計的參數(shù)計算結(jié)果

圖1－56 nQuery Advisor 7.0關(guān)于例1-23樣本量估計的參數(shù)設(shè)置與計算結(jié)果

SAS9.2軟件實現(xiàn):

%let u={5 10.5 13.5 12};

%let coe={0-1 1 0};

%let r={1 1 1 1};

proc IML;

start MGT1(a，s，G，sd，power);error=0;

if(a＞1|a＜0)then do;error=1;print“error”“Test significance level must be in 0-1”;end;

if(s^=1 ＆ s^=2)then do;error=1;print“error”“s=1 or 2”;end;

if(sd ＜0)then do;error=1;print“error”“Standard deviation must be ＞ =0”;end;

if(G＜0|ceil(G)^=G)then do;error=1;print“error”“The Number of groups must be positive integer”;end;

if(power＞100|power＜1)then do;error=1;print“error”“Power(%)must be in 1-100”;end;

if(sum(＆coe)^=0)then do;error=1;print“error”“The sum of coe must be 1”;end;

if(error=1)then stop;

if(error=0)then do;

C=sum(＆u#＆coe);D=sqrt(sum(＆coe##2/＆r));

es=abs(C)/(sd#D);n=2;*n of group 1;

do until(pw ＞ =power/100);total=n#sum(＆r);*the total N;ncp=n#es##2;df1=1;df2=total-G;if(s=2)then do;f1=finv(a，df1，df2);f2=finv(1-a，df1，df2);

pw=probf(f1，df1，df2，ncp)+1-probf(f2，df1，df2，ncp);

end;if(s=1)then do;f=finv(1-2*a，df1，df2);

pw=1-probf(f，df1，df2，ncp);

end;n=n+0.01;end;

n=ceil(n-0.01);

print a［label=“Test significance level”］

s［label=“1 or 2 sided test”］

G［label=“Number of groups”］

sd［label=“Common standard deviation”］

es［label=“Effect size”］

power［label=“Power(%)”］

n［label=“n per group”］;end;

finish MGT1;

run MGT1(0.05，2，4，6，90);quit;

SAS運行結(jié)果:

圖1－57 SAS9.2關(guān)于例1-23樣本量估計的參數(shù)設(shè)置與計算結(jié)果

1.3.1.3 雙因素方差分析

方法:根據(jù)O'Brien和 Muller(1993)〔2〕提出的方法，雙因素方差分析要考慮兩個因素及其交互項對樣本量的影響。而知道任何一個因素的檢驗效能可以計算另一個因素及交互因素的檢驗效能。各個因素檢驗效能的估計是建立在各自的自由度及非中心參數(shù)的F分布。其各自檢驗效能的計算公式為:

此時的μi是A因素某個水平下B因素的均數(shù);μj是B因素某個水平下A因素的均數(shù);μij是某個水平某個因素的均數(shù);ˉμ是所有數(shù)據(jù)的均數(shù)。

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

【例1－24】一項臨床試驗旨在評價兩種作用于心臟收縮的新藥控制收縮壓的療效，采用2×3析因設(shè)計，考慮兩個因素:一個是性別(因素 A)，有男、女2個水平;另一個是藥物(因素B)，有3個水平，分別是新藥C、新藥D和陽性對照藥，評價指標(biāo)為收縮壓水平，詳細(xì)背景見文獻〔17〕。預(yù)實驗數(shù)據(jù)顯示:男性服用新藥C、新藥D和陽性對照藥的平均收縮壓分別為130mmHg、128mmHg、125mmHg;而女性服用新藥 C、新藥D和陽性對照藥的平均收縮壓分別為125mmHg、121mmHg、118mmHg;估計標(biāo)準(zhǔn)差為6 mm-Hg。采用平衡設(shè)計，設(shè)定藥物間的檢驗效能為90%，試估計樣本量。

nQuery Advisor 7.0實現(xiàn):設(shè)定檢驗水準(zhǔn) α=0.05;檢驗效能取1-βB=90% 。由題意可知，μ11=130，μ12=128，μ13=125，μ21=125，μ22=121，μ23=118，σ=6;

在nQuery Advisor 7.0主菜單選擇:

Goal:Make Conclusion Using:⊙Means

Number of Groups:⊙ ＞Two

Analysis Method:⊙Test

方法框中選擇:Two-way analysis of variance。

Assistants:⊙Compute Effect Size

在彈出的計算窗口將各參數(shù)鍵入，如圖1－58所示，結(jié)果為 VA=10.028、VB=6.000、VAB=0.222。將計算結(jié)果V傳輸至主對話框，鍵入其他參數(shù)后求得n=14(圖1－59)。

圖1－58 nQuery Advisor 7.0關(guān)于例1-24樣本量估計的參數(shù)計算結(jié)果

圖1－59 nQuery Advisor 7.0關(guān)于例1-24樣本量估計的參數(shù)設(shè)置與計算結(jié)果

SAS9.2軟件實現(xiàn):

%let means={130 128 125，125 121 118};

PROC IML;

start MGT2(a，lA，lB，sd，Apower，Bpower，AB-power);error=0;

if(a＞1|a＜0)then do;error=1;print“error”“Test significance level must be in 0-1”;end;

if(sd ＜0)then do;error=1;print“error”“Standard deviation must be ＞ =0”;end;

if(lA＜0|ceil(lA)^=lA)then do;error=1;print“error”“The number of factor A levels must be positive integer”;end;

if(lB＜0|ceil(lB)^=lB)then do;error=1;print“error”“The number of factor A levels must be positive integer”;end;

if(error=1)then stop;

if(error=0)then do;

*compute V;Va=sum((＆means［，:］-＆means［:］)##2)/lA;Vb=sum((＆means［:，］-＆means［:］)##2)/lB;Vab=sum((＆means-＆means［:］)##2)/(lA#lB)-Va-Vb;*interaction;

Aes=Va/(sd##2);Bes=Vb/(sd##2);ABes=Vab/(sd##2);

n=2;*n of every group;

do until(Apw＞=Apower/100＆Bpw＞=Bpower/100＆ABpw＞=ABpower/100);total=n*lA*lB;*the total N;df=lA*lB*(n-1);Ancp=total*Aes;dfa=lA-1;fa=finv(1-a，dfa，df);

Apw=1-probf(fa，dfa，df，Ancp);Bncp=total*Bes;dfb=lB-1;fb=finv(1-a，dfb，df);

Bpw=1-probf(fb，dfb，df，Bncp);ABncp=total*ABes;dfab=(lA-1)*(lB-1);fab=finv(1-a，dfab，df);

ABpw=1-probf(fab，dfab，df，ABncp);n=n+0.01;end;

n=ceil(n-0.01);total=ceil(total);

Apw=100*Apw;Bpw=100*Bpw;ABpw=100*ABpw;

print a［label=“Test significance level”］

lA［label=“Number of factor A levels”］

lB［label=“Number of factor B levels”］

Va［label=“Variance in means，VA”］

Vb［label=“Variance in means，VB”］

Vab［label=“Variance in means，VAB”］

sd［label=“Common standard deviation”］

Apw［label=“Power for A(%)”］

Bpw［label=“Power for B(%)”］

ABpw［label=“Power for AB(%)”］

n［label=“n per group”］;end;finish MGT2;

run MGT2(0.05，2，3，6，0，90，0);quit;

SAS運行結(jié)果:

圖1－60 nQuery Advisor 7.0關(guān)于例1-24樣本量估計的參數(shù)設(shè)置與計算結(jié)果

致謝 我們在此感謝nQuery Advisor軟件全球供銷商“愛爾蘭Statistical Solutions有限公司”為本研究提供的支持和幫助。

1．Elashoff JD．nQuery Advisor User's Guide．Ireland:Statistical Solutions Ltd．，2007．

2．O'Brien RG，Muller KE．Applied analysis of variance in behavioral science．New York:Marcel Dekker，1993:297-344．

3．Dixon WJ，Massey FJ．Introduction to Statistical Analysis．4th Edition．New York:McGraw-Hill，1983．

4．Overall JE，Doyle SR．Estimating sample sizes for repeated measures designs．Controlled Clinical Trials，1994，15:100-123．

5．Muller KE，Barton CN．Approximate power for repeated-measures ANOVA lacking sphericity．Journal of the American Statistical Association，1989，84:549-555．

6．Machin D，Campbell MJ．Statistical Tables for Design of Clinical Tri-als．Oxford:Blackwell Scientific Publications，1987．

7．Moser BK，Stevens GR，Watts CL．The two-sample t test versus satterthwaite approximate F test．Commun．Statist．-Theory Meth，1989，18:1963-3975．

8．Diletti E，Hauschke D，Steinijans VW．Sample size determination for bioequivalence assessment by means of confidence intervals．Int．Journal of Clinical Pharmacology，1991:29．

9．Noether GE．Sample size determination for some common nonparametric statistics．Journal of the American Statistical Association，1987，82:645-647．

10．Kolassa J．A comparison of size and power calculations for the Wilcoxon statistic for ordered categorical data．Statistics in Medicine，1995，14:1577-1581．

11．Senn S．Cross-over Trials in Clinical Research(2nd Edition)．New York:Wiley，2002．

12．Schuirmann DJ．A comparison of the two one-sided tests procedure and the power approach for assessing the equivalence of average bioavailability．Journal of Pharmacokinetics and Biopharmaceutics，1987，15:657-680．

13．Phillips KE．Power of the two one-sided tests procedure in bioequivalence．Journal of Pharmacokinetics and Biopharmaceutics，1990，18:137-143．

14．Owen DB．A special case of a bivariate non-central t-distribution．Biometrika 1965，52:437-446．

15．Chow SC，Liu JP．Design and Analysis of Bioavailability and Bioequivalence Studies．New York:Marcel Dekker，Inc．，1992．

16．Hauschke D，Kieser M，Diletti E，Burke M．Sample size determination for proving equivalence based on the ratio of two means for normally distributed data．Statistics in Medicine，1999，18:93-105．

17．SASInstitute，Inc．Getting Started with the SAS Power and Sample Size Application．North Carolina:SAS Institute，Inc．，2004．