郭春雪,胡良平,2*
(1.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢中心,北京 100850;2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會(huì),北京 100029
基于R軟件估計(jì)樣本含量與檢驗(yàn)效能及其應(yīng)用
郭春雪1,胡良平1,2*
(1.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢中心,北京 100850;2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會(huì),北京 100029
本文目的是使讀者快速掌握用R軟件估計(jì)樣本含量和檢驗(yàn)效能的方法。通過R軟件中的stats包中的三個(gè)函數(shù),即power.t.test()、power.prop.test()和power.anova.test(),可以很方便地估計(jì)若干種場(chǎng)合下的樣本含量或檢驗(yàn)效能。事實(shí)表明:R軟件易于獲取、易學(xué)易用、功能強(qiáng)大、適用面寬,可以方便快捷地解決試驗(yàn)設(shè)計(jì)中的樣本含量與檢驗(yàn)效能估計(jì)問題。
R軟件;樣本含量;檢驗(yàn)效能;假設(shè)檢驗(yàn);均值;率
1.1 估計(jì)樣本含量與檢驗(yàn)效能的前提條件[1]
關(guān)于估計(jì)樣本含量與檢驗(yàn)效能的概念和前提條件,盡管在文獻(xiàn)[1]中已作了介紹,為了便于讀者閱讀本文,此處仍以概要的形式總結(jié)如下。在試驗(yàn)設(shè)計(jì)中,擬對(duì)定量指標(biāo)的平均值或定性指標(biāo)的率進(jìn)行假設(shè)檢驗(yàn)時(shí),常需提供與結(jié)果精確度、評(píng)價(jià)指標(biāo)、設(shè)計(jì)類型和比較類型有關(guān)的前提條件。
1.1.1 與結(jié)果精確度有關(guān)的前提條件
①定出檢驗(yàn)水準(zhǔn):即事先規(guī)定本批試驗(yàn)允許犯Ⅰ型(或假陽性)錯(cuò)誤的概率α,通常規(guī)定α=0.05,同時(shí)應(yīng)明確單雙側(cè)檢驗(yàn)。α定得越小,研究所需樣本含量就越大。
②提出期望的檢驗(yàn)效能(或稱把握度)1-β:即在指定的α水準(zhǔn)下,若比較的總體之間確實(shí)存在著差別,該試驗(yàn)可以發(fā)現(xiàn)差別的概率。檢驗(yàn)效能越大,所需樣本含量越多。在科研設(shè)計(jì)時(shí),檢驗(yàn)效能一般取0.8或以上比較適宜。
③估計(jì)實(shí)驗(yàn)過程中的樣本損耗。假設(shè)研究者估計(jì)本批實(shí)驗(yàn)過程中將有10%的受試者脫落而無法完成實(shí)驗(yàn),則應(yīng)將通過計(jì)算得到的樣本量除以0.9,將此時(shí)得到的結(jié)果作為該實(shí)驗(yàn)最終需要的樣本量。
1.1.2 與評(píng)價(jià)指標(biāo)有關(guān)的前提條件
必須知道由樣本推斷總體的一些信息。在比較兩總體均數(shù)或概率之間的差別是否具有統(tǒng)計(jì)學(xué)意義時(shí),需要知道總體參數(shù)間差值δ的信息。如兩總體均數(shù)間的差值δ=μ1-μ2的信息(或有關(guān)于μ1和μ2的估計(jì)值),兩總體概率間的差值δ=π1-π2的信息(或有關(guān)于π1和π2的估計(jì)值)。此外,確定兩均數(shù)比較的樣本含量時(shí),還需要有關(guān)總體標(biāo)準(zhǔn)差σ的信息(或有關(guān)于總體標(biāo)準(zhǔn)差σ的估計(jì)值)。若希望進(jìn)行非劣效性檢驗(yàn)、等效性檢驗(yàn)或優(yōu)效性檢驗(yàn)時(shí),需要提供在臨床上有意義的界值δ(此界值一般應(yīng)由多位不同地區(qū)且學(xué)術(shù)權(quán)威性高、經(jīng)驗(yàn)豐富的臨床和統(tǒng)計(jì)學(xué)專家共同討論來商定)。這些信息可以通過查閱資料、借鑒前人的經(jīng)驗(yàn)或進(jìn)行預(yù)試驗(yàn)尋找參考值。
1.1.3 與設(shè)計(jì)類型和比較類型有關(guān)的前提條件
前面提到“兩總體”,其真實(shí)含義是指所采用的是“單因素兩水平設(shè)計(jì)(常簡稱為成組設(shè)計(jì))”。換句話說,擬采用什么試驗(yàn)設(shè)計(jì)類型(除了單因素兩水平設(shè)計(jì)之外,還有單組設(shè)計(jì)、配對(duì)設(shè)計(jì)、單因素多水平設(shè)計(jì)、某種特定的多因素設(shè)計(jì)等設(shè)計(jì)類型)是估計(jì)樣本含量的重要前提條件之一;而擬采用的比較類型(包括差異性檢驗(yàn)、非劣效性檢驗(yàn)、等效性檢驗(yàn)或優(yōu)效性檢驗(yàn))也是估計(jì)樣本含量的重要前提條件之一。
1.2 R軟件中可用于估計(jì)樣本含量與檢驗(yàn)效能的程序包及函數(shù)[2-5]
在R軟件的stats包中,有三個(gè)函數(shù),即power.t.test()、power.prop.test()和power.anova.test(),可用于估計(jì)樣本含量或檢驗(yàn)效能。
在R軟件的sample.size包中,n.ttest( )、samplesize-package( )這兩個(gè)函數(shù)可用于估計(jì)樣本含量。其中,n.ttest()函數(shù)可用于配對(duì)設(shè)計(jì)和非配對(duì)設(shè)計(jì)(即成組設(shè)計(jì))一元定量資料t檢驗(yàn)時(shí),估計(jì)樣本含量;而samplesize-package( )函數(shù)可用于多種成組設(shè)計(jì)場(chǎng)合下假設(shè)檢驗(yàn)時(shí)估計(jì)樣本含量(注:前述提及的“多種成組設(shè)計(jì)場(chǎng)合”指“成組設(shè)計(jì)與配對(duì)設(shè)計(jì)一元定量資料t檢驗(yàn)時(shí)”、“Welch近似t檢驗(yàn)時(shí)”、“有序資料中帶有或不帶有結(jié)的Wilcoxon-Mann-Whitney檢驗(yàn)時(shí)”,對(duì)最后的場(chǎng)合,R軟件的samplesize包中,還有n.wilcox.ord()函數(shù)可用于估計(jì)樣本含量)。
值得一提的是,在R軟件的samplesize4surveys包中,有12個(gè)函數(shù),即b4ddm()、b4ddp()、b4dm()、b4dp()、b4m()、b4p();ss4ddmH()、ss4ddpH()、ss4dmH()、ss4dpH()、ss4dH()、ss4pH(),其中,前6個(gè)函數(shù)用于6種場(chǎng)合下估計(jì)檢驗(yàn)效能,后6個(gè)函數(shù)用于前述6種場(chǎng)合下估計(jì)樣本含量。
以上提及的6種場(chǎng)合分別為:成組設(shè)計(jì)一元定量資料均值的雙側(cè)差異性檢驗(yàn)和單側(cè)差異性檢驗(yàn)、單組設(shè)計(jì)一元定量資料均值的假設(shè)檢驗(yàn);成組設(shè)計(jì)一元定性資料比例或率的雙側(cè)差異性檢驗(yàn)和單側(cè)差異性檢驗(yàn)、單組設(shè)計(jì)一元定性資料比例或率的假設(shè)檢驗(yàn)。
值得注意的是,R軟件中的內(nèi)容十分豐富,其“程序包、小插件和函數(shù)”多如牛毛,而且每項(xiàng)內(nèi)容放置在何處,可能的確無人精準(zhǔn)知曉。只能是發(fā)現(xiàn)什么,調(diào)用什么,很難窮盡!
特別提示:當(dāng)用戶加載了程序包“samplesize”后,就可在R軟件環(huán)境中,使用如下命令:>help(samplesize),就可進(jìn)入有關(guān)樣本大小估計(jì)的幫助窗口,此窗口內(nèi)列示出了幾十個(gè)以“samplesize”或“sample.size”開頭的程序包或函數(shù);以“power”開頭的用于估計(jì)樣本大小和檢驗(yàn)效能的函數(shù)被放置在“stats”程序包中,使用命令“>help(stats)”后可進(jìn)入此程序包的幫助信息查詢窗口,選擇其中的“index”,就可按26個(gè)字母順序去查看相應(yīng)字母開頭的函數(shù)(例如,選擇字母P,就可迅速顯示此程序包中以字母P開頭的全部函數(shù)名及其功能的解說信息)。
2.1 評(píng)價(jià)指標(biāo)為定量變量的場(chǎng)合
【例1】某研究者觀察氯沙坦與伊貝沙坦治療對(duì)伴高尿酸血癥的原按發(fā)性高血壓患者血清尿酸水平的影響并評(píng)價(jià)其降壓療效。采用多中心、隨機(jī)、雙盲、平行對(duì)照設(shè)計(jì)。預(yù)試驗(yàn)的結(jié)果表明,收縮壓改變值的情況見表1。使用雙側(cè)差異性檢驗(yàn)評(píng)價(jià)兩種藥物的降壓效果的差別是否具有統(tǒng)計(jì)學(xué)意義。取α=0.05,β=0.20,試估計(jì)該試驗(yàn)所需的樣本含量。
表1 兩組患者治療6周后收縮壓下降幅度(mmHg)
解答:需要獲得兩樣本均值之差量delta的數(shù)值和兩樣本標(biāo)準(zhǔn)差之均值sd的數(shù)值。然后,調(diào)用power.t.test()函數(shù)。
> delta=14.87-13.29;delta
[1] 1.58
以上是求出兩樣本均值之差量delta的數(shù)值為1.58。
> sd=(6.10+5.84)/2;sd
[1] 5.97
以上是求出兩樣本標(biāo)準(zhǔn)差之均值sd的數(shù)值為5.97。
> power.t.test(power=0.80,sig.level=0.05,delta=1.58,sd=5.97)
以上語句的目的是調(diào)用power.t.test()函數(shù),其中的四個(gè)參數(shù)分別給定了具體的數(shù)值。事實(shí)上,還有三個(gè)參數(shù)取默認(rèn)值,第一個(gè)為設(shè)計(jì)類型:type=c(“two.sample”,“one.sample”,“paired”),默認(rèn)值為“two.sample”,即成組設(shè)計(jì);第二個(gè)為備擇假設(shè):alternative=c(“two.sided”,“one.sided”),默認(rèn)值為“two.sided”;第三個(gè)為“strict=T or F”,等號(hào)后面只能選定一個(gè),其默認(rèn)值為“strict=F”或“strict=FALSE”,其含義是:指定在雙側(cè)檢驗(yàn)時(shí)是否使用嚴(yán)格解釋。還剩下一個(gè)參數(shù)(即n)的值未給定,需要R軟件計(jì)算。
以下是輸出結(jié)果:
Two-sample t test power calculation
n = 225.08
delta = 1.58
sd = 5.97
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in*each*group
由輸出結(jié)果可知:每組應(yīng)選取約226例。
【例2】某研究者觀察氯沙坦與伊貝沙坦治療對(duì)伴高尿酸血癥的原發(fā)性高血壓患者血清尿酸水平的影響并評(píng)價(jià)其降壓療效。采用多中心、隨機(jī)、雙盲、平行對(duì)照設(shè)計(jì)。隨機(jī)選取320例受試者,治療6周后收縮壓改變值的情況見表2。使用雙側(cè)差異性檢驗(yàn)評(píng)價(jià)兩種藥物的降壓效果的差別是否具有統(tǒng)計(jì)學(xué)意義。取α=0.05,β=0.20,試估計(jì)該試驗(yàn)的檢驗(yàn)效能。
表2 兩組患者治療6周后收縮壓下降幅度(mmHg)
解答:需要獲得兩樣本均值之差量delta的數(shù)值和兩樣本標(biāo)準(zhǔn)差之均值sd的數(shù)值。然后,調(diào)用power.t.test()函數(shù)。
> delta=14.87-13.29;delta
[1] 1.58
以上是求出兩樣本均值之差量delta的數(shù)值為1.58。
> sd=(6.10+5.84)/2;sd
[1] 5.97
以上是求出兩樣本標(biāo)準(zhǔn)差之均值sd的數(shù)值為5.97。
> power.t.test(n=160,sig.level=0.05,delta=1.58,sd=5.97)
以上語句的目的是調(diào)用power.t.test()函數(shù),其中的四個(gè)參數(shù)分別給定了具體的數(shù)值。三個(gè)默認(rèn)參數(shù)前已述及,不再贅述。還剩下一個(gè)參數(shù)(即power)的值未給定,需要R軟件計(jì)算。
以下是輸出結(jié)果:
Two-sample t test power calculation
n = 160
delta = 1.58
sd = 5.97
sig.level = 0.05
power = 0.6554376
alternative = two.sided
NOTE: n is number in*each*group
以上結(jié)果表明:每組用160例,其檢驗(yàn)效能僅為65.54%<80.0%(常規(guī)的要求),犯假陰性錯(cuò)誤的概率(34.46%)過大。
【例3】為觀察神經(jīng)功能康復(fù)情況,使用三種方法分別治療腦卒中抑郁患者,估計(jì)治療后三種方法的SSS評(píng)分均值分別為11.0、10.0、9.0,組間方差相等且都為9,組內(nèi)方差分別為4、5、6、9四種取值條件下,取α=0.05,β=0.10,要求得到三組間差別有統(tǒng)計(jì)學(xué)意義的結(jié)論,每組各需要患者多少例(三組所需患者人數(shù)相等)?
解答:情形一,在組內(nèi)方差為4的條件下;
> power.anova.test(group=3,between.var=9,within.var=4,sig.level=0.05,power=0.90)
此條件下輸出的結(jié)果如下:
Balanced one-way analysis of variance power calculation
groups = 3
n = 4.017349
between.var = 9
within.var = 4
sig.level = 0.05
power = 0.9
NOTE: n is number in each group
以上結(jié)果表明:每組只需要5例。
情形二,在組內(nèi)方差為5的條件下;
> power.anova.test(group=3,between.var=9,within.var=5,sig.level=0.05,power=0.90)
此條件下輸出的結(jié)果如下:
Balanced one-way analysis of variance power calculation
groups = 3
n = 4.688307
between.var = 9
within.var = 5
sig.level = 0.05
power = 0.9
NOTE: n is number in each group
以上結(jié)果表明:每組只需要5例。
情形三,在組內(nèi)方差為6的條件下;
> power.anova.test(group=3,between.var=9,within.var=6,sig.level=0.05,power=0.90)
此條件下輸出的結(jié)果如下:
Balanced one-way analysis of variance power calculation
groups = 3
n = 5.36743
between.var = 9
within.var = 6
sig.level = 0.05
power = 0.9
NOTE: n is number in each group
以上結(jié)果表明:每組只需要6例。
情形四,在組內(nèi)方差為9的條件下;
> power.anova.test(group=3,between.var=9,within.var=9,sig.level=0.05,power=0.90)
此條件下輸出的結(jié)果如下:
Balanced one-way analysis of variance power calculation
groups = 3
n = 7.431865
between.var = 9
within.var = 9
sig.level = 0.05
power = 0.9
NOTE: n is number in each group
以上結(jié)果表明:每組只需要8例。
2.2 評(píng)價(jià)指標(biāo)為定性變量的場(chǎng)合
【例4】一個(gè)新的抗腫瘤藥物A與臨床有效藥物B對(duì)照進(jìn)行臨床試驗(yàn),選取一定數(shù)目且符合要求的患者隨機(jī)均分成兩組,分別接受A藥和B藥治療。預(yù)試驗(yàn)結(jié)果為A藥的有效率是58.0%,B藥的有效率是46.0%。欲使用雙側(cè)差異性檢驗(yàn)評(píng)價(jià)兩種藥物的降壓效果的差別且希望得出具有統(tǒng)計(jì)學(xué)意義的結(jié)果,取α=0.05,β=0.20,試估計(jì)該試驗(yàn)中各組至少需要多大的樣本含量?
解答:需要獲得估計(jì)成組設(shè)計(jì)兩比例或率差異性檢驗(yàn)時(shí)樣本含量所需要的基本信息: p1=0.58、p2=0.46、sig.level=0.05、power=0.80,將各組樣本含量n留作待估計(jì)的參數(shù)。然后,調(diào)用power.prop.test()函數(shù)。
> power.prop.test(p1=0.58,p2=0.46,sig.level=0.05,power=0.80)
以上語句的目的是調(diào)用power.prop.test()函數(shù),其中的四個(gè)參數(shù)分別給定了具體的數(shù)值。兩個(gè)默認(rèn)參數(shù)如下。第一個(gè)為備擇假設(shè):alternative=c(“two.sided”,“one.sided”),默認(rèn)值為“two.sided”;第二個(gè)為“strict=T or F”,等號(hào)后面只能選定一個(gè),其默認(rèn)值為“strict=F”或“strict=FALSE”,其含義是:指定在雙側(cè)檢驗(yàn)時(shí)是否使用嚴(yán)格解釋;還剩下一個(gè)參數(shù)(即power)的值未給定,需要R軟件計(jì)算。
以下是輸出結(jié)果:
Two-sample comparison of proportions power calculation
n = 270.9126
p1= 0.58
p2= 0.46
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in*each*group
以上結(jié)果表明:各組需要約271例。
【例5】在一個(gè)Ⅱ期臨床試驗(yàn)中,已知對(duì)照組有效率p1=30%、試驗(yàn)組有效率p2=60%、三個(gè)限制條件分別為:伽瑪C=15%、伽瑪E=15%、伽瑪Delta=15%、對(duì)照組與試驗(yàn)組樣本含量的比例分別為1:3與1:1兩種條件下,試估計(jì)各組的樣本含量。
【解答】在臨床試驗(yàn)中,有所謂的“雙臂試驗(yàn)”,即一個(gè)試驗(yàn)組與一個(gè)對(duì)照組比較。在R軟件的sample.size包中,有一個(gè)Sample.Size()函數(shù)可用于Ⅱ期臨床試驗(yàn)且采用雙臂優(yōu)化設(shè)計(jì)中比較兩個(gè)比例或率(此法由Mayo等于2010年提出,可采用固定或靈活的分配比例,可以基于多種限制條件下進(jìn)行優(yōu)化設(shè)計(jì))時(shí),估計(jì)樣本含量或檢驗(yàn)效能。
情形一,對(duì)照組與試驗(yàn)組樣本含量的比例為1:3條件下估計(jì)各組的樣本含量
> Sample.Size(0.3, 0.6, 0.15, 0.15, 0.15, Allratio_c = 1, Allratio_e = 3)
以上語句的目的是調(diào)用Sample.Size()函數(shù),其中七個(gè)參數(shù)的數(shù)值均被給定,最后兩個(gè)參數(shù)若不出現(xiàn),就取默認(rèn)值,即1:1。
以下是輸出的結(jié)果:
Specified values for parameters:
Response rates:
control = 0.3 experiment = 0.6
Upper bounds for constriants:
gammaC = 0.15 gammaE = 0.15 gammaDelta = 0.15
以上內(nèi)容實(shí)際上是給定的前提條件。
Required sample sizes:
[1] Optimal Design:
nc = 20 ne = 20 n = 40
[2] 1 to 1 Allocation Design:
nc = 20 ne = 20 n = 40
[3] 1 to 3 Allocation Design:
nc = 13 ne = 39 n = 52
第一部分“見上面的[1]”給出了優(yōu)化設(shè)計(jì)下的兩組各需要20例;第二部分“見上面的[2]”給出了1:1條件下的兩組各需要20例;第三部分“見上面的[3]”給出了1:3條件下的對(duì)照組需要13例、試驗(yàn)組需要39例。
說明:由以上給定條件和輸出結(jié)果可知,即便給定的前提條件是兩組樣本含量之比為1:3,但也將默認(rèn)的前提條件1:1多對(duì)應(yīng)的樣本含量估計(jì)出來了。
情形二,對(duì)照組與試驗(yàn)組樣本含量的比例為1:1條件下估計(jì)各組的樣本含量
> Sample.Size(0.3, 0.6, 0.15, 0.15, 0.15)
此語句與前面的語句相比,最后兩個(gè)參數(shù)取默認(rèn)值,即對(duì)照組與試驗(yàn)組樣本含量之比為1:1。
以下是輸出的結(jié)果:
Specified values for parameters:
Response rates:
control = 0.3 experiment = 0.6
Upper bounds for constriants:
gammaC = 0.15 gammaE = 0.15 gammaDelta = 0.15
以上內(nèi)容實(shí)際上是給定的前提條件。
Required sample sizes:
[1] Optimal Design:
nc = 20 ne = 20 n = 40
[2] 1 to 1 Allocation Design:
nc = 20 ne = 20 n = 40
第一部分“見上面的[1]”給出了優(yōu)化設(shè)計(jì)下的兩組各需要20例;第二部分“見上面的[2]”給出了1:1條件下的兩組各需要20例。
值得注意的是:此設(shè)計(jì)并沒有交代清楚:sig.level=?power=?
筆者認(rèn)為:這種設(shè)計(jì)要慎用!
[1] 張效嘉, 胡良平. 精神衛(wèi)生科研如何嚴(yán)格遵守試驗(yàn)設(shè)計(jì)四原則之重復(fù)原則[J]. 四川精神衛(wèi)生, 2016, 29(4): 303-306.
[2] 黃文, 王正林. 數(shù)據(jù)挖掘: R語言實(shí)戰(zhàn)[M]. 北京: 電子工業(yè)出版社, 2015: 34-39.
[3] 李詩羽, 張飛, 王正林. 數(shù)據(jù)分析:R語言實(shí)戰(zhàn)[M]. 北京: 電子工業(yè)出版社, 2015: 88-156.
[4] 方匡南, 朱建平, 姜葉飛. R 數(shù)據(jù)分析:方法與案例詳解[M]. 北京: 電子工業(yè)出版社, 2015: 54-168.
[5] Joseph Adler. R語言核心技術(shù)手冊(cè)[M]. 2版. 劉思喆, 李艦, 陳鋼, 等譯. 北京: 電子工業(yè)出版社, 2015: 417-421.
[6] 胡良平, 陶麗新. 臨床試驗(yàn)設(shè)計(jì)與統(tǒng)計(jì)分析[M]. 北京: 軍事醫(yī)學(xué)科學(xué)出版社, 2013: 101-134.
(本文編輯:吳俊林)
The estimation of sample size and power and its application based on R software
GuoChunxue1,HuLiangping1,2*
(1.ConsultingCenterofBiomedicalStatistics,AcademyofMilitaryMedicalSciences,Beijing100850,China; 2.SpecialtyCommitteeofClinicalScientificResearchStatisticsofWorldFederationofChineseMedicineSocieties,Beijing100029,China
*Correspondingauthor:HuLiangping,E-mail:lphu812@sina.com)
The paper aims to help the readers to grasp the method of estimating the sample size and power with R software. By using the three functions [power.t.test(), power.prop.test() and power.anova.test()] of stat in R software, it is convenient for readers to realize the estimation of sample size and power by using R software under the different situations. The methods of estimating the sample size and power by R software were introduced through several real examples in this article. Since that R is very easy for people to learn and use, and has the advantages of powerful functions and wide application, the users can solve the concrete problems concerned with the estimation of sample size and power in experimental designs conveniently and easily.
R software; Sample size; Power; Hypothesis testing; Mean value; Rate
*通信作者:胡良平,E-mail:lphu812@sina.com)
R195.1
A
10.11886/j.issn.1007-3256.2016.06.004
國家高技術(shù)研究發(fā)展計(jì)劃課題資助(2015AA020102)
2016-12-06)