哈佛的一位統計分析師分享了一個名為 %mcpower 的巨集程式，能夠輕鬆的完成在線性迴歸、羅吉斯迴歸以及卜瓦松迴歸的蒙地卡羅檢定力模擬與估計。
%mcpower 主要包含 %glmpower 和 %manypower 兩個巨集程式。程式可從此處下載：
其中，%glmpower 是無法直接被執行的，所有的參數都是定義在 %manypower 裡面，格式如下：
%macro manypower(alt1=., alt2=., alt3=.,
n1=.,n2=.,n3=.,n4=.,n5=.,n6=.,n7=.,n8=.,n9=.,n10=.,
outcomedist = "NORMAL", null=0, scaley=1,
xdist="NORMAL", xloc=., xscale =.,
nreps=,plotrows=2, descending = desc, intercept = 0,
plotstyle=graphic, outputpower=YES, outputnsubs=YES, test=FALSE);不意外地，巨集變數很多，但因為模擬本來就是一件即為複雜的工作，所以參數多也是很正常的。這些參數主要可以分類為三大群：
假設巨集變數：
- null：所要檢定之虛無假設的值。預設值為 0。
- alt1, alt2, alt3：所要檢定之對立假設(Ha)的數值。若是線性迴歸，則檢定的項目是斜率，若是羅吉斯迴歸，則檢定的項目是odds ratio。若是卜瓦松迴歸，則檢定的項目是risk ratio。無預設值。
- n1, n2, ..., n10：每個模擬資料集裡面的樣本個數。最多可設定十個。無預設值。
- outcomedist：依變數的分配。可填入 "NORMAL"、"BERNOULLI" 和 "POISSON"。預設值為 "NORMAL"。(注意，雙引號要寫上）
- scaley：殘差變異數。若使用線性迴歸(即outcomedist = "NORMAL")，則可以省略這個參數的設定。預設值為 1。
- intercept：模擬的模型裡面的截距項。預設值為 0。
- xdist：應變數的分配。預設值是 "NORMAL"。
- xloc：應變數分配的第一個參數值。若 xdist = "NORMAL" 時，此值表示平均數。若 xdist = "POISSON" 時，則此值表示 lambda。若 xdist = "BERNOULLI" 時，則此值表示 p。
- xscale：應變數分配的第二個參數值。若 xdist 為 "POISSON" 或 "BERNOULLI" 時，則此值要省略，因為這兩個分配只有一個參數值。若 xdist = "NORMAL" 時，則可設定變異數。
- outputpower：輸出檢定力表格。預設值為 YES。
- outputnsubs：輸出樣本數表格。預設值為 YES。
- plotstyle：輸出表格的格式。如果設定為 sg，則圖形會用 PROC SGPLOT 程序來畫。如果設定為 graphic，則圖形會用 PROC GPLOT 來畫。如果沒有設定，則圖形就不會輸出。預設值為 graphic。
- plotrows：若 plotstyle = sg 時，則可設定圖形要有幾行。預設值為 2。
範例：
%manypower(null = 0, scaley=1,xdist = "bernoulli", xloc=.5, xscale=.,
nreps = 100, n1=1000,n2=1100,n3=1200,n4=1400,n5=1700,n6=2000,
alt1 = 1.15, alt2 = 1.20, alt3 = 1.25,
plotrows=2,plotstyle=graphic, outcomedist="poisson", intercept = 1);結果：
CONTACT INFORMATION
Your comments and questions would be welcome and valued. Contact the author at:
Ken Kleinman
Harvard Medical School and Harvard Pilgrim Health Care
133 Brookline Ave.
Boston, MA 02215
ken.kleinman@gmail.com
沒有留言:
張貼留言
要問問題的人請在文章下方的intensedebate欄位留言，請勿使用blogger預設的意見表單。今後用blogger意見表單留言的人我就不回應了。