幾年前曾經寫過幾篇關於介紹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
沒有留言:
張貼留言
要問問題的人請在文章下方的intensedebate欄位留言，請勿使用blogger預設的意見表單。今後用blogger意見表單留言的人我就不回應了。