陶麗新，劉一松，胡良平

(1.首都醫(yī)科大學(xué)公共衛(wèi)生學(xué)院，北京 100069；2.世界中醫(yī)藥聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029；3.北京醫(yī)普科諾科技有限公司，北京 100190；4.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心，北京 100850*通信作者：胡良平：E-mail: lphu812@sina.com)

?

基于SAS軟件實(shí)現(xiàn)樣本含量估計及應(yīng)用

陶麗新1,2，劉一松3，胡良平2,4*

(1.首都醫(yī)科大學(xué)公共衛(wèi)生學(xué)院，北京 100069；2.世界中醫(yī)藥聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京 100029；3.北京醫(yī)普科諾科技有限公司，北京 100190；4.軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心，北京 100850*通信作者：胡良平：E-mail: lphu812@sina.com)

本文的目的是使讀者能方便快捷地運(yùn)用SAS軟件中的POWER過程和GLMPOWER過程實(shí)現(xiàn)樣本含量估計。在不同的場合下估計樣本含量需要提供不同的前提條件，即使僅限于假設(shè)檢驗時估計樣本含量，也必須進(jìn)一步弄清對應(yīng)的設(shè)計類型、結(jié)果變量的性質(zhì)、某些先驗知識和對結(jié)果精確度的要求。本文通過一些實(shí)例，介紹了估計樣本含量與檢驗效能的SAS實(shí)現(xiàn)方法。讀者只需要修改本文中所呈現(xiàn)的SAS程序中的少量參數(shù)，就可方便地使用SAS軟件實(shí)現(xiàn)樣本含量與檢驗效能估計。事實(shí)說明，盡管SAS軟件非常難學(xué)難用，但借助現(xiàn)成的SAS程序，讀者可以輕松自如地解決很多與統(tǒng)計分析有關(guān)的具體問題。

SAS軟件；樣本含量；檢驗效能；均值；標(biāo)準(zhǔn)差；率

1 估計樣本含量的必要性與前提條件

在開展試驗性科學(xué)研究中，往往需要在相同試驗條件下做多次重復(fù)試驗，以便使隨機(jī)變量的變化規(guī)律能充分顯露出來。文獻(xiàn)[1]就“重復(fù)原則與樣本含量”、“重復(fù)的三層含義”、“估計樣本含量和檢驗效能的意義”、“估計樣本含量需要的前提條件”等內(nèi)容，作了較為詳細(xì)的介紹，本文將直接運(yùn)用SAS軟件中的POWER過程和GLMPOWER過程估計樣本含量和檢驗效能。

2 定量資料假設(shè)檢驗時樣本含量估計[2-5]

【例1】某研究者欲比較兩種物理療法增加肌肉彈性的效果。希望采用t檢驗比較兩總體均值之間差別并希望得出具有統(tǒng)計學(xué)意義的檢驗結(jié)果，設(shè)檢驗水準(zhǔn)α為0.05，檢驗效能1-β為0.9，假定的均值和標(biāo)準(zhǔn)差的數(shù)值(在實(shí)踐中，這些數(shù)據(jù)應(yīng)基于文獻(xiàn)資料或預(yù)實(shí)驗的結(jié)果)將列在下面。擬回答的問題是：至少需要檢測多少例受試對象？

標(biāo)準(zhǔn)治療方法的肌肉彈性均值為13，如果新療法的肌肉彈性被定位14.0、14.5和15.0三種情況，且假定分別有兩種可能的標(biāo)準(zhǔn)差，即σ=1.2與σ=1.7；假定研究者還希望兩組的樣本含量具有下列三種比例關(guān)系，即：①每個組的樣本含量相同(1:1)；②新療法的樣本含量是標(biāo)準(zhǔn)療法的兩倍(1:2)；③新療法的樣本含量是標(biāo)準(zhǔn)療法的三倍(1:3)。

【分析與解答】根據(jù)前面給出的條件可知，均值有3種可能的取值、標(biāo)準(zhǔn)差有兩種可能的取值、兩組樣本含量的比例有三種可能的取值，故這樣一共有3×2×3=18種情況。也就是說，要在18種條件下估計所需要的樣本含量。

使用下面的SAS程序，可以獲得與上面18種條件對應(yīng)的樣本含量估計結(jié)果：

proc power;

twosamplemeans

groupmeans = (13 14) (13 14.5) (13 15)

stddev = 1.2 1.7

groupweights = 1 | 1 2 3

power = 0.9

ntotal = .;

run;

【程序說明】在power過程中使用twosamplemeans語句來決定18種情況下為了達(dá)到90%的檢驗效能所需要的樣本含量。將選項“ntotal=.”指定樣本含量為結(jié)果變量；用“groupmenas=”選項指定假設(shè)的均值，對每種情況的均值使用匹配標(biāo)記，即用圓括號隔開不同的均值組合條件；使用“stddev=”選項為每種情況指定標(biāo)準(zhǔn)差；使用“groupweights=”選項指定權(quán)重比例，可以再使用匹配標(biāo)記，即使用“|”符號隔開不同的取值。

【輸出結(jié)果】

Computed N Total

【結(jié)果解讀】最好的設(shè)計類型是(均值差最大為2，小的標(biāo)準(zhǔn)差為1.2和平衡設(shè)計)，達(dá)到至少0.9的檢驗效能需要的樣本含量為18(n1=n2=9)；最差的設(shè)計為(小的均值差為1，大的標(biāo)準(zhǔn)差為1.7，1:3的非平衡設(shè)計)，達(dá)到至少0.9的檢驗效能需要樣本含量為164(n1=41、n2=123)。在“固定設(shè)計中”指定的檢驗效能為0.9，在“計算的樣本含量”表格中輸出根據(jù)樣本比例調(diào)整的實(shí)際的檢驗效能(見上面輸出結(jié)果中倒數(shù)第2列)。

值得一提的是，在上面的SAS程序中，若給定總樣本含量，而將power設(shè)置為“缺省值”，就是在各種不同條件下，估計檢驗效能。如將上面的SAS程序修改如下：

procpower; twosamplemeans groupmeans=(1314)(1314.5)(1315) stddev=1.21.7 groupweights=1|123 power=. ntotal=100;run;

讀者可以自己運(yùn)行上面這段SAS程序，會得到各種條件下的檢驗效能估計值，因篇幅所限，此輸出結(jié)果從略。在下面的其他SAS程序中，也可進(jìn)行同樣的修改。

【分析與解答】這是一個結(jié)果變量為定量變量的單因素三水平設(shè)計試驗研究的樣本含量估計問題。所需的SAS程序如下。

odshtml; procpower; onewayanovagroupmeans=290|658|763stddev=174155127 groupweights=(111)alpha=0.05Ntotal=.power=0.8; run;quit; odshtmlclose;

【程序說明】選項“groupmeans =290 | 658 | 763”定義各水平組指標(biāo)總體均值的估計值；選項“stddev = 174 155 127”定義幾種總體標(biāo)準(zhǔn)差的估計值；選項“groupweights = (1 1 1)”定義各水平組樣本含量的分配比例。對于類似的問題，上述程序中的數(shù)字需要根據(jù)用戶的實(shí)際情況進(jìn)行修改。

【輸出結(jié)果】

ComputedNTotalIndexStdDevActualPowerNTotal11740.8661221550.9331231270.9119

【結(jié)果解讀】分析結(jié)果表明：本研究各類人群只需隨機(jī)抽取4例，3類人群共12例即可滿足要求。

【例3】用三種方法治療腦卒中抑郁患者，觀察神經(jīng)功能康復(fù)狀況。估計治療后三種方法斯堪地那維亞卒中量表(Scemelinavian Stroke Scale,SSS)評分均數(shù)分別為11.0、23.0、9.0，標(biāo)準(zhǔn)差分別為3、3、2。希望各組樣本含量的比例或權(quán)重(本例取1:1:1)、各組的平均標(biāo)準(zhǔn)差(保守的做法，可采用所有組的最大標(biāo)準(zhǔn)差)、允許犯第Ⅰ類(或假陽性)錯誤的概率α(通常取其值為0.05)、給定檢驗效能1-β的值[相當(dāng)于給定允許犯第Ⅱ類(或假陰性)錯誤的概率β(通常取其值為0.2或0.1)，本例取1-β=0.9]。如果要得到三組之間的差別以及兩兩比較的結(jié)果都有統(tǒng)計學(xué)意義的結(jié)論，問每組各需要多少例患者？

【分析與解答】 這是做單因素三水平設(shè)計一元定量資料均值比較及兩兩比較之前進(jìn)行樣本含量估計的問題。所需的SAS程序如下。

data d11_20;

input method $ score CellWgt;

datalines;

A 11.0 1

B 23.0 1

C 9.0 1

;

run;

ods html;

PROC GLMPOWER data=d11_20;

class method;

model score=method;

weight CellWgt;

contrast "A vs. B" method 1 -1 0;

contrast "A vs. C" method 1 0 -1;

contrast "B vs. C" method 0 1 -1;

POWER stddev = 3.0 alpha = 0.05 power = 0.90 ntotal = .；

run; quit;

ods html close;

【程序說明】 data步定義影響因素各水平組指標(biāo)的總體均值的估計值及各水平組的樣本含量分配比例；GLMPOWER過程中的選項“stddev=3.0”中的3.0為各水平組指標(biāo)的總體標(biāo)準(zhǔn)差的估計值，本例取三個水平組總體標(biāo)準(zhǔn)差估計值的最大值帶入計算，也可以同時考慮多個總體標(biāo)準(zhǔn)差的估計值，將選項寫成“stddev=2.0 2.5 3.0”或“stddev=2.0 to 3.0 by 0.5”的形式；程序中的數(shù)字需要根據(jù)實(shí)際情況進(jìn)行修改。

具體地說，在數(shù)據(jù)步中：method實(shí)際上指“試驗因素”，若用factor更直觀，它是一個字符型變量，其具體取值或水平可用一個字母或一個單詞、甚至可用一句話(當(dāng)字符個數(shù)超過8時，不能僅用一個美元符號，如最長為20個字符，可寫成：method $20.)來表示；“score”指各組的算術(shù)平均值，若用“mean”更直觀；“CellWgt”指各組的樣本含量比例或權(quán)重，本例取1:1:1，根據(jù)實(shí)際需要，也可取其它比例關(guān)系，如1:2:3；

在過程步中，有三個Contrast語句，相當(dāng)于單因素三水平之間的兩兩比較，分別表示A與B比較、A與C比較、B與C比較；POWER語句中寫了4個內(nèi)容，可以寫在同一行上，寫成4行是為了看起來更清楚一些，前3行代表提供的3個已知條件，它們分別為平均的或最大的樣本標(biāo)準(zhǔn)差、α值和1-β的值；而“ntotal = .”表示需要估計總樣本含量，即三組樣本含量之總和。

【輸出結(jié)果】

ComputedNTotalIndexTypeSourceTestDFErrorDFActualPowerNTotal1Effectmethod260.98992ContrastAvs.B160.98193ContrastAvs.C11410.9001444ContrastBvs.C160.9979

【結(jié)果解讀】由計算結(jié)果可知：若要求通過單因素三水平設(shè)計一元定量資料的方差分析及兩兩比較有90%的把握發(fā)現(xiàn)三個組SSS評分的總體均值不完全相等，則每組只需3例患者，三組共需9例患者；但若要求通過單因素三水平設(shè)計一元定量資料的方差分析及兩兩比較有90%的把握發(fā)現(xiàn)三個組中任何兩組的SSS評分的總體均值不相等，則每組需要48例患者，三組共需144例患者。

為什么會出現(xiàn)這樣大的差異呢？問題的癥結(jié)在于A組與C組的均值過于接近，它們之間的數(shù)值之差量為11-9=2，要想把總體上存在如此小的“差距”都能發(fā)現(xiàn)出來，必須具有非常大的樣本含量。

試想：若A組的均值不是11，而是16，可以計算得出：每組只需要5例，三組共需要15例即可。程序修改之后的輸出結(jié)果如下：

ComputedNTotalIndexTypeSourceTestDFErrorDFActualPowerNTotal1Effectmethod260.97692ContrastAvs.B1120.922153ContrastAvs.C1120.922154ContrastBvs.C160.9979

【分析與解答】 這是做具有一個重復(fù)測量因素的兩因素設(shè)計一元定量資料均值比較之前進(jìn)行樣本含量估計的問題。所需的SAS程序如下。

%MACROcanshu; alpha=0.05;beta=0.10;k=5;pc=0.75;delta=10;deltae=190.641;deltau=457.333;%MENDcanshu;%MACROrepeatmesure(a,b,c,d);FILEPRINT;u005=put((-probit(alpha/2)),8.3);u02=put((-probit(beta)),8.3);n=int((2*(deltau+(1+(k-1)*pc)*deltae/k)*(u005+u02)*(u005+u02)/(delta*delta)))+1;PUT#2@15‘兩組各需’n‘例?！?RUN;%MENDrepeatmesure;DATA;%canshu;%repeatmesure(a,b,c,d);RUN;

【程序說明】第2行中“alpha=0.05”中的0.05為可能犯第Ⅰ類錯誤的概率，“beta=0.10”中的0.10為可能犯第Ⅱ類錯誤的概率，“k=5”中的5為重復(fù)測量次數(shù)，“pc=0.75”中的0.75為條件相關(guān)系數(shù)，“delta=10”中的10為兩組間可以識別的最小差異，“deltae=190.641”中的190.641為重復(fù)測量誤差，“deltau=457.333”中的457.333為個體間差異的方差。

【輸出結(jié)果】每組需要129例。

3 生存分析時樣本含量估計

【例5】某臨床試驗擬比較TIPS手術(shù)與脊髓改道手術(shù)延長出血性食管曲張患者的生存時間有無差別，試驗組為TIPS手術(shù)組，對照組為脊髓改道手術(shù)組，隨訪期至少1年。根據(jù)既往研究可知對照組總體1年生存率是45%，預(yù)期試驗組總體1年生存率將達(dá)到65%，采用Log-Rank檢驗比較兩組生存率，試估計樣本含量(采用雙側(cè)檢驗，檢驗水準(zhǔn)α=0.05，檢驗效能1-β=0.85)。

【分析與解答】 對照組生存率為45%、試驗組生存率為65%，對照組和試驗組的終點(diǎn)事件發(fā)生率分別為：π1=1-0.65=0.35，πc=1-0.45=0.55。這是對生存資料分析之前的樣本含量估計問題。所需的SAS程序如下。

%letpt=0.35;%letpc=0.55;%letalpha=0.05;%letbeta=0.15; dataa6_1; theta=log10(1-&pt)/log10(1-&pc); theta1=round(theta,0.001); e=((theta+1)/(theta-1))**2*(probit(1-&alpha/2)+probit(1-&beta))**2; e1=int(e)+1; n=2/(&pt+&pc)*e; n1=int(n)+1; fileprint; PUT#3@15‘風(fēng)險比為‘theta1’;兩組終點(diǎn)事件發(fā)生總例數(shù)不小于‘e1’;兩組所需總的樣本含量為‘n1’。’; run;

【程序說明】程序當(dāng)中pt、pc、alpha、beta分別代表試驗組終點(diǎn)事件發(fā)生率πt、對照組終點(diǎn)事件發(fā)生率πc、檢驗水準(zhǔn)α、犯假陰性錯誤的概率β。

【輸出結(jié)果】風(fēng)險比為0.539；兩組終點(diǎn)事件發(fā)生總例數(shù)不小于101；兩組所需總的樣本含量為223。

說明：樣本含量與檢驗效能估計的內(nèi)容非常豐富，本文的目的僅僅希望起到拋磚引玉的作用。這方面的統(tǒng)計軟件和書籍[6-7]很多，讀者需要時，還可參閱有關(guān)文獻(xiàn)。

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

[2] SAS Institute Inc. SAS/STAT 9.3 User’s Guide[M]. Cary, NC: SAS Institute Inc, 2011: 3361-3400,5729-5962.

[3] 胡良平. 統(tǒng)計學(xué)三型理論在實(shí)驗設(shè)計中的應(yīng)用[M]. 北京: 人民軍醫(yī)出版社, 2006: 215-265.

[4] 胡良平. 科研設(shè)計與統(tǒng)計分析[M]. 北京: 軍事醫(yī)學(xué)科學(xué)出版社, 2012: 136-227.

[5] 胡良平. SAS常用統(tǒng)計分析教程[M]. 北京: 電子工業(yè)出版社, 2015: 235-252.

[6] Shein-Chung Chow, Jun Shao, Hansheng Wang. Sample Size Calculations in Clinical Research[M]. 2nd Edition. Chapman & Hall/CRC, 2008: 49-408.

[7] Thomas P. Ryan. Sample Size Determination and Power[M]. Wiley, 2013: 17-352.

(本文編輯：吳俊林)

Sample size estimation and its application by using SAS software

TAOLi-Xin1,2,LIUYi-song3,HULiang-ping2,4*

(1.PublicHealthAcademy,CapitalUniversityofMedicalSciences,Beijing100069,China;2.SpecialtyCommitteeofClinicalScientificResearchStatisticsofWorldFederationofChineseMedicineSocieties,Beijing100029,China;3.BeijingImprove-QualityTechCo.,Ltd.,Beijing100190,China;4.ConsultingCenterofBiomedicalStatistics,AcademyofMilitaryMedicalSciences,Beijing100850,China*Correspondingauthor:HULiang-ping,E-mail:lphu812@sina.com)

The paper aims to help the readers to estimate the sample size by using the POWER and GLMPOWER procedures in SAS software conveniently and easily. It is necessary for a user to supply the prerequisites of sample size estimation under the different situations, even if limited under the case of the hypothesis testing. Usually the premise conditions of sample size estimation, which should be given, are as follows: the design type, the characters of outcomes, the prior knowledge, and the accuracy of the results. The methods of estimating the sample size and power by SAS were introduced through several real examples in this article. It is convenient for the readers to realize the estimation of sample size and power by using SAS. The only thing for them to do is to modify a few parameters in the given SAS programs. Although SAS software is very difficult for people to learn and use, as a matter of fact, the users can solve the concrete problems concerned with statistical analysis by means of revisions of ready-made SAS programs easily.

SAS software; Sample size; Power; Mean value; Standard deviation; Rate

R195.1

A doi:10.11886/j.issn.1007-3256.2016.05.004

2016-10-11)