公告

[公告]
2014/01/17
由於已經是faculty的關係,不太有足夠時間寫部落格。因此更新的速度會相當緩慢。再加上近幾年來SAS GLOBAL FORUM沒有出現讓我覺得驚艷的技術文件,所以能分享的文章相對也減少許多。若有人推薦值得分享的SAS技術文件,請利用『問題討論區』告知。

2013/07/19
臉書留言板的功能因為有不明原因故障,因此特此移除。而intensedebate的留言板因管理不易,也一併移除。目前已經開啟內建的 G+ 留言系統,所以請有需要留言的朋友,可直接至『問題討論區』裡面留言。


2015年9月7日 星期一

A Generalized Approach to Estimating Sample Sizes

原文載點:http://support.sas.com/resources/papers/proceedings12/336-2012.pdf

幾年前曾經寫過幾篇關於介紹PROC POWER和PROC GLMPOWER功能的技術文件。本篇特地再引述一篇計算樣本的技術文件的用意是這一篇特別介紹了如何用模擬的方式來做檢定力分析(Power Analysis),尤其當模型中涉及到交互作用項時。想看原文的可以直接跳到原文的第四頁來看。

若我們要用下面這個模型來做減肥藥的分析:


其中Yweightloss是服藥後減掉的體重,Xdrug是一個二項變數(1=減肥藥, 0=安慰劑),Xcalorie是一個連續型卡路里變數。這個模型還包含四個未知參數:截距項是α,兩個獨立線性項(β1, β2)還有一個交互作用項β3。假設根據之前的研究這四個參數的估計值依序是是10, 5, -0.004, 0.0025,而卡路里變數服從一個常態分配N(2500, 1000^2),而最後面的殘差項則服從常態分配N(0, 5),最後樣本裡面約只有30%的人有服用真正的減肥藥,而另外70%的人則服用安慰劑,則可以用下面的程式來模擬一百萬筆數據:


其中前面九個%let statement是來宣告模擬樣本數,估計參數值,以及卡路里變數和殘差項的分配參數值。下面的data step則是跑了一個一百萬次的迴圈,其中用ranuni(seed)來從均勻分布U(0,1)裡面模擬一個隨機亂數,讓這個隨機亂數<0.3時則該模擬值屬於那30%的實驗組(服用減肥藥,drug=1),而其餘的則變成對照組(服用安慰劑,drug=0)。至於卡路里變數則是用其平均數加上標準差乘上一個隨機亂數來生成,但我個人認為比較快的方法是直接套用rand("Normal", &mean_cal, &sd_cal)來生成。同理,殘差項也可以用rand('Normal', 0, &sd_error)來生成。最後剩下依變數Y,則把剛剛生成的那些變數套進模型公式裡面即可生成。

接著,用一個PROC GLM來分析這一百萬筆模擬數據:


估計參數如下所示:


這個結果被一個ODS OUPUT指令另存在一個叫做pe的新檔中。最後,為了計算要多少樣本才能檢測到那個交互作用項的存在,則需要用到下面這個公式:

(PS)原文誤植。期望值裡面應該是T值而非X值。

在顯著水準95%及檢定力80%的情況下,α和β分別是0.05以及0.8,所以在雙尾檢定下,這個樣本估計程式如下所示:


其中z_alpha和z_beta是拿來計算Zα和Zβ這兩個臨界值。expected_t則是計算T值的期望值:

E(T) = |T|/N^0.5

最後便可得到


由此可知需要至少152個樣本才能檢定出那個交互作用項出來。不過根據我個人的經驗,只模擬一份資料不夠,大概模擬個一千次再來求平均所得到的樣本會比較穩健。

CONTACT INFORMATION 
Wuchen Zhao
University of Southern California Department of Preventive Medicine
Division of Biostatistics
1540 Alcazar Street
CHP 222 Los Angeles, CA 90089-9010
E-mail: zhaowuchen@gmail.com
CODE { display: block; /* fixes a strange ie margin bug */ font-family: Courier New; font-size: 8pt; overflow:auto; background: #f0f0f0 url(http://klcintw.images.googlepages.com/Code_BG.gif) left top repeat-y; border: 1px solid #ccc; padding: 10px 10px 10px 21px; max-height:200px; height:200px; // for IE6 line-height: 1.2em; }