基于R軟件估計樣本含量與檢驗(yàn)效能及其應(yīng)用

郭春雪，胡良平,2*

(1.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心，北京 100850；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029

基于R軟件估計樣本含量與檢驗(yàn)效能及其應(yīng)用

郭春雪1，胡良平1,2*

(1.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心，北京 100850；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029

本文目的是使讀者快速掌握用R軟件估計樣本含量和檢驗(yàn)效能的方法。通過R軟件中的stats包中的三個函數(shù)，即power.t.test()、power.prop.test()和power.anova.test()，可以很方便地估計若干種場合下的樣本含量或檢驗(yàn)效能。事實(shí)表明：R軟件易于獲取、易學(xué)易用、功能強(qiáng)大、適用面寬，可以方便快捷地解決試驗(yàn)設(shè)計中的樣本含量與檢驗(yàn)效能估計問題。

R軟件；樣本含量；檢驗(yàn)效能；假設(shè)檢驗(yàn)；均值；率

1 估計樣本含量與檢驗(yàn)效能的概述

1.1 估計樣本含量與檢驗(yàn)效能的前提條件[1]

關(guān)于估計樣本含量與檢驗(yàn)效能的概念和前提條件，盡管在文獻(xiàn)[1]中已作了介紹，為了便于讀者閱讀本文，此處仍以概要的形式總結(jié)如下。在試驗(yàn)設(shè)計中，擬對定量指標(biāo)的平均值或定性指標(biāo)的率進(jìn)行假設(shè)檢驗(yàn)時，常需提供與結(jié)果精確度、評價指標(biāo)、設(shè)計類型和比較類型有關(guān)的前提條件。

1.1.1 與結(jié)果精確度有關(guān)的前提條件

①定出檢驗(yàn)水準(zhǔn)：即事先規(guī)定本批試驗(yàn)允許犯Ⅰ型(或假陽性)錯誤的概率α，通常規(guī)定α=0.05，同時應(yīng)明確單雙側(cè)檢驗(yàn)。α定得越小，研究所需樣本含量就越大。

②提出期望的檢驗(yàn)效能(或稱把握度)1-β：即在指定的α水準(zhǔn)下，若比較的總體之間確實(shí)存在著差別，該試驗(yàn)可以發(fā)現(xiàn)差別的概率。檢驗(yàn)效能越大，所需樣本含量越多。在科研設(shè)計時，檢驗(yàn)效能一般取0.8或以上比較適宜。

③估計實(shí)驗(yàn)過程中的樣本損耗。假設(shè)研究者估計本批實(shí)驗(yàn)過程中將有10%的受試者脫落而無法完成實(shí)驗(yàn)，則應(yīng)將通過計算得到的樣本量除以0.9，將此時得到的結(jié)果作為該實(shí)驗(yàn)最終需要的樣本量。

1.1.2 與評價指標(biāo)有關(guān)的前提條件

必須知道由樣本推斷總體的一些信息。在比較兩總體均數(shù)或概率之間的差別是否具有統(tǒng)計學(xué)意義時，需要知道總體參數(shù)間差值δ的信息。如兩總體均數(shù)間的差值δ=μ1-μ2的信息(或有關(guān)于μ1和μ2的估計值)，兩總體概率間的差值δ=π1-π2的信息(或有關(guān)于π1和π2的估計值)。此外，確定兩均數(shù)比較的樣本含量時，還需要有關(guān)總體標(biāo)準(zhǔn)差σ的信息(或有關(guān)于總體標(biāo)準(zhǔn)差σ的估計值)。若希望進(jìn)行非劣效性檢驗(yàn)、等效性檢驗(yàn)或優(yōu)效性檢驗(yàn)時，需要提供在臨床上有意義的界值δ(此界值一般應(yīng)由多位不同地區(qū)且學(xué)術(shù)權(quán)威性高、經(jīng)驗(yàn)豐富的臨床和統(tǒng)計學(xué)專家共同討論來商定)。這些信息可以通過查閱資料、借鑒前人的經(jīng)驗(yàn)或進(jìn)行預(yù)試驗(yàn)尋找參考值。

1.1.3 與設(shè)計類型和比較類型有關(guān)的前提條件

前面提到“兩總體”，其真實(shí)含義是指所采用的是“單因素兩水平設(shè)計(常簡稱為成組設(shè)計)”。換句話說，擬采用什么試驗(yàn)設(shè)計類型(除了單因素兩水平設(shè)計之外，還有單組設(shè)計、配對設(shè)計、單因素多水平設(shè)計、某種特定的多因素設(shè)計等設(shè)計類型)是估計樣本含量的重要前提條件之一；而擬采用的比較類型(包括差異性檢驗(yàn)、非劣效性檢驗(yàn)、等效性檢驗(yàn)或優(yōu)效性檢驗(yàn))也是估計樣本含量的重要前提條件之一。

1.2 R軟件中可用于估計樣本含量與檢驗(yàn)效能的程序包及函數(shù)[2-5]

在R軟件的stats包中，有三個函數(shù)，即power.t.test()、power.prop.test()和power.anova.test()，可用于估計樣本含量或檢驗(yàn)效能。

在R軟件的sample.size包中，n.ttest( )、samplesize-package( )這兩個函數(shù)可用于估計樣本含量。其中，n.ttest()函數(shù)可用于配對設(shè)計和非配對設(shè)計(即成組設(shè)計)一元定量資料t檢驗(yàn)時，估計樣本含量；而samplesize-package( )函數(shù)可用于多種成組設(shè)計場合下假設(shè)檢驗(yàn)時估計樣本含量(注：前述提及的“多種成組設(shè)計場合”指“成組設(shè)計與配對設(shè)計一元定量資料t檢驗(yàn)時”、“Welch近似t檢驗(yàn)時”、“有序資料中帶有或不帶有結(jié)的Wilcoxon-Mann-Whitney檢驗(yàn)時”，對最后的場合，R軟件的samplesize包中，還有n.wilcox.ord()函數(shù)可用于估計樣本含量)。

值得一提的是，在R軟件的samplesize4surveys包中，有12個函數(shù)，即b4ddm()、b4ddp()、b4dm()、b4dp()、b4m()、b4p()；ss4ddmH()、ss4ddpH()、ss4dmH()、ss4dpH()、ss4dH()、ss4pH()，其中，前6個函數(shù)用于6種場合下估計檢驗(yàn)效能，后6個函數(shù)用于前述6種場合下估計樣本含量。

以上提及的6種場合分別為：成組設(shè)計一元定量資料均值的雙側(cè)差異性檢驗(yàn)和單側(cè)差異性檢驗(yàn)、單組設(shè)計一元定量資料均值的假設(shè)檢驗(yàn)；成組設(shè)計一元定性資料比例或率的雙側(cè)差異性檢驗(yàn)和單側(cè)差異性檢驗(yàn)、單組設(shè)計一元定性資料比例或率的假設(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)樣本大小估計的幫助窗口，此窗口內(nèi)列示出了幾十個以“samplesize”或“sample.size”開頭的程序包或函數(shù)；以“power”開頭的用于估計樣本大小和檢驗(yàn)效能的函數(shù)被放置在“stats”程序包中，使用命令“>help(stats)”后可進(jìn)入此程序包的幫助信息查詢窗口，選擇其中的“index”，就可按26個字母順序去查看相應(yīng)字母開頭的函數(shù)(例如，選擇字母P，就可迅速顯示此程序包中以字母P開頭的全部函數(shù)名及其功能的解說信息)。

2 幾種場合下估計樣本含量與檢驗(yàn)效能的實(shí)例[2-6]

2.1 評價指標(biāo)為定量變量的場合

【例1】某研究者觀察氯沙坦與伊貝沙坦治療對伴高尿酸血癥的原按發(fā)性高血壓患者血清尿酸水平的影響并評價其降壓療效。采用多中心、隨機(jī)、雙盲、平行對照設(shè)計。預(yù)試驗(yàn)的結(jié)果表明，收縮壓改變值的情況見表1。使用雙側(cè)差異性檢驗(yàn)評價兩種藥物的降壓效果的差別是否具有統(tǒng)計學(xué)意義。取α=0.05，β=0.20，試估計該試驗(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ù)，其中的四個參數(shù)分別給定了具體的數(shù)值。事實(shí)上，還有三個參數(shù)取默認(rèn)值，第一個為設(shè)計類型：type=c(“two.sample”,“one.sample”,“paired”)，默認(rèn)值為“two.sample”，即成組設(shè)計；第二個為備擇假設(shè)：alternative=c(“two.sided”,“one.sided”)，默認(rèn)值為“two.sided”；第三個為“strict=T or F”，等號后面只能選定一個，其默認(rèn)值為“strict=F”或“strict=FALSE”，其含義是：指定在雙側(cè)檢驗(yàn)時是否使用嚴(yán)格解釋。還剩下一個參數(shù)(即n)的值未給定，需要R軟件計算。

以下是輸出結(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】某研究者觀察氯沙坦與伊貝沙坦治療對伴高尿酸血癥的原發(fā)性高血壓患者血清尿酸水平的影響并評價其降壓療效。采用多中心、隨機(jī)、雙盲、平行對照設(shè)計。隨機(jī)選取320例受試者，治療6周后收縮壓改變值的情況見表2。使用雙側(cè)差異性檢驗(yàn)評價兩種藥物的降壓效果的差別是否具有統(tǒng)計學(xué)意義。取α=0.05，β=0.20，試估計該試驗(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ù)，其中的四個參數(shù)分別給定了具體的數(shù)值。三個默認(rèn)參數(shù)前已述及，不再贅述。還剩下一個參數(shù)(即power)的值未給定，需要R軟件計算。

以下是輸出結(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ī)的要求)，犯假陰性錯誤的概率(34.46%)過大。

【例3】為觀察神經(jīng)功能康復(fù)情況，使用三種方法分別治療腦卒中抑郁患者，估計治療后三種方法的SSS評分均值分別為11.0、10.0、9.0，組間方差相等且都為9，組內(nèi)方差分別為4、5、6、9四種取值條件下，取α=0.05，β=0.10，要求得到三組間差別有統(tǒng)計學(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 評價指標(biāo)為定性變量的場合

【例4】一個新的抗腫瘤藥物A與臨床有效藥物B對照進(jìn)行臨床試驗(yàn)，選取一定數(shù)目且符合要求的患者隨機(jī)均分成兩組，分別接受A藥和B藥治療。預(yù)試驗(yàn)結(jié)果為A藥的有效率是58.0%，B藥的有效率是46.0%。欲使用雙側(cè)差異性檢驗(yàn)評價兩種藥物的降壓效果的差別且希望得出具有統(tǒng)計學(xué)意義的結(jié)果，取α=0.05，β=0.20，試估計該試驗(yàn)中各組至少需要多大的樣本含量？

解答：需要獲得估計成組設(shè)計兩比例或率差異性檢驗(yàn)時樣本含量所需要的基本信息： p1=0.58、p2=0.46、sig.level=0.05、power=0.80，將各組樣本含量n留作待估計的參數(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ù)，其中的四個參數(shù)分別給定了具體的數(shù)值。兩個默認(rèn)參數(shù)如下。第一個為備擇假設(shè)：alternative=c(“two.sided”,“one.sided”)，默認(rèn)值為“two.sided”；第二個為“strict=T or F”，等號后面只能選定一個，其默認(rèn)值為“strict=F”或“strict=FALSE”，其含義是：指定在雙側(cè)檢驗(yàn)時是否使用嚴(yán)格解釋；還剩下一個參數(shù)(即power)的值未給定，需要R軟件計算。

以下是輸出結(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】在一個Ⅱ期臨床試驗(yàn)中，已知對照組有效率p1=30%、試驗(yàn)組有效率p2=60%、三個限制條件分別為：伽瑪C=15%、伽瑪E=15%、伽瑪Delta=15%、對照組與試驗(yàn)組樣本含量的比例分別為1:3與1:1兩種條件下，試估計各組的樣本含量。

【解答】在臨床試驗(yàn)中，有所謂的“雙臂試驗(yàn)”，即一個試驗(yàn)組與一個對照組比較。在R軟件的sample.size包中，有一個Sample.Size()函數(shù)可用于Ⅱ期臨床試驗(yàn)且采用雙臂優(yōu)化設(shè)計中比較兩個比例或率(此法由Mayo等于2010年提出，可采用固定或靈活的分配比例，可以基于多種限制條件下進(jìn)行優(yōu)化設(shè)計)時，估計樣本含量或檢驗(yàn)效能。

情形一，對照組與試驗(yàn)組樣本含量的比例為1:3條件下估計各組的樣本含量

> Sample.Size(0.3, 0.6, 0.15, 0.15, 0.15, Allratio_c = 1, Allratio_e = 3)

以上語句的目的是調(diào)用Sample.Size()函數(shù)，其中七個參數(shù)的數(shù)值均被給定，最后兩個參數(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è)計下的兩組各需要20例；第二部分“見上面的[2]”給出了1:1條件下的兩組各需要20例；第三部分“見上面的[3]”給出了1:3條件下的對照組需要13例、試驗(yàn)組需要39例。

說明：由以上給定條件和輸出結(jié)果可知，即便給定的前提條件是兩組樣本含量之比為1:3，但也將默認(rèn)的前提條件1:1多對應(yīng)的樣本含量估計出來了。

情形二，對照組與試驗(yàn)組樣本含量的比例為1:1條件下估計各組的樣本含量

> Sample.Size(0.3, 0.6, 0.15, 0.15, 0.15)

此語句與前面的語句相比，最后兩個參數(shù)取默認(rèn)值，即對照組與試驗(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è)計下的兩組各需要20例；第二部分“見上面的[2]”給出了1:1條件下的兩組各需要20例。

值得注意的是：此設(shè)計并沒有交代清楚：sig.level=？power=？

筆者認(rèn)為：這種設(shè)計要慎用！

[1] 張效嘉, 胡良平. 精神衛(wèi)生科研如何嚴(yán)格遵守試驗(yàn)設(shè)計四原則之重復(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ù)手冊[M]. 2版. 劉思喆, 李艦, 陳鋼, 等譯. 北京: 電子工業(yè)出版社, 2015: 417-421.

[6] 胡良平, 陶麗新. 臨床試驗(yàn)設(shè)計與統(tǒng)計分析[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ā)展計劃課題資助(2015AA020102)

2016-12-06)