餓死抬頭,這些年來很多人一直找我做過許多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。有興趣的人可以點下面這個連結看線上版本的第十二章:

SAS for Mixed Models

不過書裡面的程式只能計算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 是「一組」的數據,所以通常要乘以二才是全部樣本的大小。

創作者介紹

ToTo 奇妙の冒險

cchien 發表在 痞客邦 PIXNET 留言(0) 人氣()