餓死抬頭，這些年來很多人一直找我做過許多power analysis來求power或樣本數，而SAS內建的PROC POWER和PROC GLMPOWER這兩個程序已經能夠把最基本的各種假設檢定、ANOVA和GLM下的power analysis給搞定，順便還可以畫出相當精美的圖形。
不過，如果遇到使用GLMM時，早期他的power analysis是無解的，也很少看有人去做過，不然就是退回到用GLM的power analysis結果來替代。但至從「SAS for Mixed Models」這本書升級到第二版時，裡面就附上了一個完整的章節教大家怎樣去做GLMM的power analysis。有興趣的人可以點下面這個連結看線上版本的第十二章：
不過書裡面的程式只能計算power，沒有辦法去推算達到某個程度的power時需要多少樣本。不過有人在網路上提供了下列這段程式碼，專門來計算Mixed model with repeated measurement的所需樣本：
n = function(m,s2b,rho,delta,alpha,power) {
floor(1+2*(qnorm(1-alpha/2)+qnorm(power))**2*s2b*(1+(m-1)*rho)/(m*delta**2))
}
n # number of patients required per group
m <- c(4) # number of repeated measures
s2b <- c(76) # between-patient variance
rho <- c(0.53) # intraclass correlation
delta <- c(5) # difference to be detected
alpha <- c(0.05) # significance level
power <- c(0.80) # power
n(m,s2b,rho,delta,alpha,power)
這程式可以在R和S+底下跑，裡面有六個參數需要使用者自行設定：
m：重複測量次數（如果你說你每個subject的重複測量次數不同，那就取最大的那個）。
s2b: between-subject的變異數，也就是Var(yij)。如果不知道這個數據，就參考過去文獻去抓一個近似值出來。
rho: 這就是ICC，亦是重複觀測值之間的相關係數。這個參數可自訂一個合理的數據。
delta: 組間可能的差距。此數據可自訂或計算(z_(1-alpha/2)+z_beta)*SE(ti-tj)。
alpha: 顯著水準。
power: 預期要達到的檢定力。
注意上面這個程式只能用在 normal data，如果不是 normal data，就得改用下面這個程式：
n = function(m,a,b,phi,rho,delta,alpha,power) {
floor(1+2*(qnorm(1-alpha/2)+qnorm(power))**2*(phi*a/b)*(1+(m-1)*rho)/(m*delta**2))
}
n # number of patients required per group
m <- c(4) # number of repeated measures
a <- 1 # denominator term for binomial data or offset term for count data
b <- 0.24 # expected variance (u(1-u) for binary data and u for count data
phi <- 1 # dispersion parameter
rho <- c(0.5) # intraclass correlation
delta <- c(0.693) # difference to be detected
alpha <- c(0.05) # significance level
power <- c(0.80) # power
n(m,a,b,phi,rho,delta,alpha,power)
這邊多了幾個新參數：
a: binomial data的分母或是count data的offset。
b: binary data和count data的期望變異數。
phi: 這個參數的中文我不會翻譯。他是count data裡面E(X)/V(X)的值。
如果沒有或不會R和S+，那我將他們改寫成下面這個SAS macro，請安心服用：
%macro n_mixed_repeat(normal= ,m= ,s2b= ,rho= ,delta= ,alpha= ,power= , a=0, b=0, phi=1);
proc iml;
%if &normal=1 %then %do;
n=floor(1+2*(quantile('NORMAL',(1-&alpha/2))+quantile('NORMAL',&power))**2*&s2b*(1+(&m-1)*&rho)/(&m*&delta**2));
%end;
%else %if &normal=0 %then %do;
n=floor(1+2*(quantile('NORMAL',(1-&alpha/2))+quantile('NORMAL',&power))**2*(&phi*&a/&b)*(1+(&m-1)*&rho)/(&m*&delta**2));
%end;
%else %do;
print 'Please input normal=1 or 0';
%end;
print n;
quit;
%mend;
另外需注意的是，算出來的 n 是「一組」的數據，所以通常要乘以二才是全部樣本的大小。
留言列表
