基于R軟件估計樣本含量與檢驗效能及其應用

郭春雪，胡良平,2*

(1.軍事醫學科學院生物醫學統計學咨詢中心，北京 100850；2.世界中醫藥學會聯合會臨床科研統計學專業委員會，北京 100029

基于R軟件估計樣本含量與檢驗效能及其應用

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

(1.軍事醫學科學院生物醫學統計學咨詢中心，北京 100850；2.世界中醫藥學會聯合會臨床科研統計學專業委員會，北京 100029

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

R軟件；樣本含量；檢驗效能；假設檢驗；均值；率

1 估計樣本含量與檢驗效能的概述

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

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

1.1.1 與結果精確度有關的前提條件

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

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

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

1.1.2 與評價指標有關的前提條件

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

1.1.3 與設計類型和比較類型有關的前提條件

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

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

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

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

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

以上提及的6種場合分別為：成組設計一元定量資料均值的雙側差異性檢驗和單側差異性檢驗、單組設計一元定量資料均值的假設檢驗；成組設計一元定性資料比例或率的雙側差異性檢驗和單側差異性檢驗、單組設計一元定性資料比例或率的假設檢驗。

值得注意的是，R軟件中的內容十分豐富，其“程序包、小插件和函數”多如牛毛，而且每項內容放置在何處，可能的確無人精準知曉。只能是發現什么，調用什么，很難窮盡！

特別提示：當用戶加載了程序包“samplesize”后，就可在R軟件環境中，使用如下命令：>help(samplesize)，就可進入有關樣本大小估計的幫助窗口，此窗口內列示出了幾十個以“samplesize”或“sample.size”開頭的程序包或函數；以“power”開頭的用于估計樣本大小和檢驗效能的函數被放置在“stats”程序包中，使用命令“>help(stats)”后可進入此程序包的幫助信息查詢窗口，選擇其中的“index”，就可按26個字母順序去查看相應字母開頭的函數(例如，選擇字母P，就可迅速顯示此程序包中以字母P開頭的全部函數名及其功能的解說信息)。

2 幾種場合下估計樣本含量與檢驗效能的實例[2-6]

2.1 評價指標為定量變量的場合

【例1】某研究者觀察氯沙坦與伊貝沙坦治療對伴高尿酸血癥的原按發性高血壓患者血清尿酸水平的影響并評價其降壓療效。采用多中心、隨機、雙盲、平行對照設計。預試驗的結果表明，收縮壓改變值的情況見表1。使用雙側差異性檢驗評價兩種藥物的降壓效果的差別是否具有統計學意義。取α=0.05，β=0.20，試估計該試驗所需的樣本含量。

表1 兩組患者治療6周后收縮壓下降幅度(mmHg)

解答：需要獲得兩樣本均值之差量delta的數值和兩樣本標準差之均值sd的數值。然后，調用power.t.test()函數。

> delta=14.87-13.29;delta

[1] 1.58

以上是求出兩樣本均值之差量delta的數值為1.58。

> sd=(6.10+5.84)/2;sd

[1] 5.97

以上是求出兩樣本標準差之均值sd的數值為5.97。

> power.t.test(power=0.80,sig.level=0.05,delta=1.58,sd=5.97)

以上語句的目的是調用power.t.test()函數，其中的四個參數分別給定了具體的數值。事實上，還有三個參數取默認值，第一個為設計類型：type=c(“two.sample”,“one.sample”,“paired”)，默認值為“two.sample”，即成組設計；第二個為備擇假設：alternative=c(“two.sided”,“one.sided”)，默認值為“two.sided”；第三個為“strict=T or F”，等號后面只能選定一個，其默認值為“strict=F”或“strict=FALSE”，其含義是：指定在雙側檢驗時是否使用嚴格解釋。還剩下一個參數(即n)的值未給定，需要R軟件計算。

以下是輸出結果：

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

由輸出結果可知：每組應選取約226例。

【例2】某研究者觀察氯沙坦與伊貝沙坦治療對伴高尿酸血癥的原發性高血壓患者血清尿酸水平的影響并評價其降壓療效。采用多中心、隨機、雙盲、平行對照設計。隨機選取320例受試者，治療6周后收縮壓改變值的情況見表2。使用雙側差異性檢驗評價兩種藥物的降壓效果的差別是否具有統計學意義。取α=0.05，β=0.20，試估計該試驗的檢驗效能。

表2 兩組患者治療6周后收縮壓下降幅度(mmHg)

解答：需要獲得兩樣本均值之差量delta的數值和兩樣本標準差之均值sd的數值。然后，調用power.t.test()函數。

> delta=14.87-13.29;delta

[1] 1.58

以上是求出兩樣本均值之差量delta的數值為1.58。

> sd=(6.10+5.84)/2;sd

[1] 5.97

以上是求出兩樣本標準差之均值sd的數值為5.97。

> power.t.test(n=160,sig.level=0.05,delta=1.58,sd=5.97)

以上語句的目的是調用power.t.test()函數，其中的四個參數分別給定了具體的數值。三個默認參數前已述及，不再贅述。還剩下一個參數(即power)的值未給定，需要R軟件計算。

以下是輸出結果：

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

以上結果表明：每組用160例，其檢驗效能僅為65.54%<80.0%(常規的要求)，犯假陰性錯誤的概率(34.46%)過大。

【例3】為觀察神經功能康復情況，使用三種方法分別治療腦卒中抑郁患者，估計治療后三種方法的SSS評分均值分別為11.0、10.0、9.0，組間方差相等且都為9，組內方差分別為4、5、6、9四種取值條件下，取α=0.05，β=0.10，要求得到三組間差別有統計學意義的結論，每組各需要患者多少例(三組所需患者人數相等)？

解答：情形一，在組內方差為4的條件下；

> power.anova.test(group=3,between.var=9,within.var=4,sig.level=0.05,power=0.90)

此條件下輸出的結果如下：

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

以上結果表明：每組只需要5例。

情形二，在組內方差為5的條件下；

> power.anova.test(group=3,between.var=9,within.var=5,sig.level=0.05,power=0.90)

此條件下輸出的結果如下：

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

以上結果表明：每組只需要5例。

情形三，在組內方差為6的條件下；

> power.anova.test(group=3,between.var=9,within.var=6,sig.level=0.05,power=0.90)

此條件下輸出的結果如下：

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

以上結果表明：每組只需要6例。

情形四，在組內方差為9的條件下；

> power.anova.test(group=3,between.var=9,within.var=9,sig.level=0.05,power=0.90)

此條件下輸出的結果如下：

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

以上結果表明：每組只需要8例。

2.2 評價指標為定性變量的場合

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

解答：需要獲得估計成組設計兩比例或率差異性檢驗時樣本含量所需要的基本信息： p1=0.58、p2=0.46、sig.level=0.05、power=0.80，將各組樣本含量n留作待估計的參數。然后，調用power.prop.test()函數。

> power.prop.test(p1=0.58,p2=0.46,sig.level=0.05,power=0.80)

以上語句的目的是調用power.prop.test()函數，其中的四個參數分別給定了具體的數值。兩個默認參數如下。第一個為備擇假設：alternative=c(“two.sided”,“one.sided”)，默認值為“two.sided”；第二個為“strict=T or F”，等號后面只能選定一個，其默認值為“strict=F”或“strict=FALSE”，其含義是：指定在雙側檢驗時是否使用嚴格解釋；還剩下一個參數(即power)的值未給定，需要R軟件計算。

以下是輸出結果：

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

以上結果表明：各組需要約271例。

【例5】在一個Ⅱ期臨床試驗中，已知對照組有效率p1=30%、試驗組有效率p2=60%、三個限制條件分別為：伽瑪C=15%、伽瑪E=15%、伽瑪Delta=15%、對照組與試驗組樣本含量的比例分別為1:3與1:1兩種條件下，試估計各組的樣本含量。

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

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

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

以上語句的目的是調用Sample.Size()函數，其中七個參數的數值均被給定，最后兩個參數若不出現，就取默認值，即1:1。

以下是輸出的結果：

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

以上內容實際上是給定的前提條件。

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]”給出了優化設計下的兩組各需要20例；第二部分“見上面的[2]”給出了1:1條件下的兩組各需要20例；第三部分“見上面的[3]”給出了1:3條件下的對照組需要13例、試驗組需要39例。

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

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

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

此語句與前面的語句相比，最后兩個參數取默認值，即對照組與試驗組樣本含量之比為1:1。

以下是輸出的結果：

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

以上內容實際上是給定的前提條件。

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]”給出了優化設計下的兩組各需要20例；第二部分“見上面的[2]”給出了1:1條件下的兩組各需要20例。

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

筆者認為：這種設計要慎用！

[1] 張效嘉, 胡良平. 精神衛生科研如何嚴格遵守試驗設計四原則之重復原則[J]. 四川精神衛生, 2016, 29(4): 303-306.

[2] 黃文, 王正林. 數據挖掘: R語言實戰[M]. 北京: 電子工業出版社, 2015: 34-39.

[3] 李詩羽, 張飛, 王正林. 數據分析：R語言實戰[M]. 北京: 電子工業出版社, 2015: 88-156.

[4] 方匡南, 朱建平, 姜葉飛. R 數據分析：方法與案例詳解[M]. 北京: 電子工業出版社, 2015: 54-168.

[5] Joseph Adler. R語言核心技術手冊[M]. 2版. 劉思喆, 李艦, 陳鋼, 等譯. 北京: 電子工業出版社, 2015: 417-421.

[6] 胡良平, 陶麗新. 臨床試驗設計與統計分析[M]. 北京: 軍事醫學科學出版社, 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

國家高技術研究發展計劃課題資助(2015AA020102)

2016-12-06)