樣本量估計及其在nQuery和SAS軟件上的實現*
——相關分析(二)

南方醫科大學公共衛生學院生物統計學系(510515) 錢晨堅 吳研鵬 段重陽 陳平雁

本刊之前已介紹了單樣本的kappa系數檢驗(二分類變量)、相關系數檢驗(連續變量)和Lin和諧系數檢驗(連續變量)等有關相關分析的樣本量估計方法[1-2]，本文將進一步介紹相關分析中單樣本Cronbachα系數檢驗、單樣本組內相關系數檢驗和兩獨立樣本Pearson相關系數檢驗的樣本量估計方法，前兩種方法主要用于量表的評價。

相關分析

1.單樣本相關性分析

(1)差異性檢驗

①Cronbachα系數檢驗

方法：Feldt等[3]提出的單樣本Cronbachα系數檢驗樣本量的估計方法，是建立在自由度為v1和v2的F分布上的，檢驗效能的計算公式如下：

(1)

式中，v1=n-1，v2=(n-1)(k-1)，n為樣本量，k為量表的條目數；α為檢驗水準，s取1表示單側檢驗，取2表示雙側檢驗，1-β為檢驗效能；F1-α(v1，v2)表示

自由度為v1和v2的F分布上側α分位數，ProbF表示F分布的累積分布函數；ρα0和ρα1分別表示已知總體的Cronbachα系數和預期總體的Cronbachα系數。

計算樣本量時，設定樣本量n的初始值為2，通過不斷增加樣本量直至檢驗效能滿足設定條件為止，最終得到的所需樣本量。

【例1】欲評價某生活質量量表(SF-36)的信度，已知該量表共36個條目，即k=36。預期該量表的Cronbachα系數為0.89，以Cronbachα系數等于0.8為目標值(即原假設ρα0=0.8)。試估計單側檢驗水準α=0.025，檢驗效能為90%，該量表的Cronbachα系數不低于0.8所需樣本量。

nQuery Advanced 8.6.0.0實現：設置單側檢驗水準α=0.025，檢驗效能1-β=90%，其他數據相應代入，在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：Test for One Coefficient(Cronbach)Alpha

在彈出的樣本量計算窗口將各參數值鍵入，如圖1所示，結果n=60。即本研究所需樣本量為60例。

圖1 nQueryAdvanced8.6.0.0關于例1樣本量估計的參數設置與計算結果

SAS 9.4軟件實現：

/*Cronbach alpha系數檢驗*/

%macro cronbach(

alpha=/*檢驗水準*/

，side=/*單雙側檢驗，1表示單側檢驗，2表示雙側檢驗*/

，r0=/*已知總體的Cronbachα系數*/

，r1=/*預期總體的Cronbachα系數*/

，k=/*量表條目數*/

，power=/*檢驗效能(%)*/

)；

data a；

r0=&r0；r1=&r1；alpha=&alpha；k=&k；side=&side；power=&power/100；

/*輸出錯誤信息*/

if(alpha>0.2 | alpha<0)then do；

error=1；

put “error：Test significance level′s range：0-0.2”；

end；

if(side^=1 & side^=2)then do；

error=1；

put “error：Side′s range：1 OR 2”；

end；

if(r0<=0 | r0>=1)then do；

error=1；

put “error：r0′s range：0

end；

if(r1<=0 | r1>=1)then do；

error=1；

put “error：r1′s range：0

end；

if(k<2)then do；

error=1；

put “error：k must be greater than 1”；

end；

if(&power>100 | &power<0)then do；

error=1；

put “error：Power′s range：0-100”；

end；

/*如有錯誤，跳出循環*/

if(error=1)then stop；

power=0；n=2；

do while(power<&power)；

f1=FINV(1-alpha/side，n-1，(k-1)*(n-1))；

f2=FINV(alpha/side，n-1，(k-1)*(n-1))；

c=(1-r1)/(1-r0)；

/*r0

if(r0 < r1)then power=100*(1-PROBF(c*f1，n-1，(k-1)*(n-1)))；

/*r0>r1*/

if(r0 > r1)then power=100*(PROBF(c*f2，n-1，(k-1)*(n-1)))；

if(power>=&power)then leave；

else n=n+1；

end；

power=round(power，0.01)；

run；

/*結果輸出*/

proc print data=a label；

var alpha side r0 r1 k power n；

label

alpha=“Test Significance Level”

side=“1 or 2 sided”

r0=“Null Coefficient Alpha”

r1=“Alternative Coefficient Alpha”

k=“Number of Raters”

power=“Power(%)”

n=“Sample Size”；

quit；

%mend cronbach；

%cronbach(alpha=0.025，side=1，r0=0.8，r1=0.89，k=36，power=90)

SAS運行結果：

圖2 SAS 9.4 關于例1樣本量估計的參數設置與計算結果

②組內相關系數檢驗

方法：Donner和Eliasziw[4]提出的單樣本組內相關系數(interclass correlation coefficient，ICC)單側檢驗的樣本量估計方法，是建立在自由度為v1和v2的F分布上的，檢驗效能的計算公式如下：

(2)

式中，v1=n-1，v2=n(m-1)，n為樣本量，m為處理組數或重復因素的水平數；α為檢驗水準，1-β為檢驗效能；F1-α(v1，v2)表示自由度為v1和v2的F分布上側α分位數，ProbF為F分布的累積分布函數；ρI0和ρI1為已知總體的組內相關系數和預期總體的組內相關系數。

計算樣本量時，設定樣本量n的初始值為2，通過不斷增加樣本量直至檢驗效能滿足設定條件為止，最終得到的n表示所需樣本量。

【例2】某研究使用葡萄糖氧化酶法對一批血液樣本進行血糖濃度測定，欲評價不同實驗者使用該法的測定結果一致性，即組內相關系數的大小。該研究下，參與一致性測定的實驗者人數為3，即m=3。預期該研究的組內相關系數為0.85，以組內相關系數大于0.75(信度良好)為目標值(即原假設ρI0=0.75)。試估計單側檢驗水準α=0.025，檢驗效能為90%，欲驗證該研究的組內相關系數不低于0.75所需血液樣本的樣本量。

nQuery Advanced 8.6.0.0實現：設置單側檢驗水準α=0.025，檢驗效能1-β=90%，其他數據相應代入，在nQuery Advanced 8.6.0.0主菜單選擇：

方法框中選擇：Test for Intraclass(Intracluster)Correlation

在彈出的樣本量計算窗口將各參數值鍵入，如圖3所示，結果n=92。即本研究所需樣本量為92例。

圖3 nQueryAdvanced8.6.0.0關于例2樣本量估計的參數設置與計算結果

SAS 9.4軟件實現：

/*組內相關系數檢驗*/

%macro ICC(

alpha=/*檢驗水準*/

，r0=/*已知總體的組內相關系數*/

，r1=/*預期總體的組內相關系數*/

，m=/*處理組數或重復因素的水平數*/

，power=/*檢驗效能(%)*/

)；

data a；

r0=&r0；r1=&r1；alpha=&alpha；m=&m；power=&power/100；

/*輸出錯誤信息*/

if(alpha>0.2 | alpha<0)then do；

error=1；

put “error：Test significance level′s range：0-0.2”；

end；

if(r0<0 | r0>=1)then do；

error=1；

put “error：r0′s range：0<=r0<1”；

end；

if(r1<0 | r1>=1)then do；

error=1；

put “error：r1′s range：0<=r1<1”；

end；

if(r0>=r1)then do；

error=1；

put “error：r1 must be greater than r0”；

end；

if(m<2)then do；

error=1；

put “error：m must be greater than 1”；

end；

if(&power>100 | &power<0)then do；

error=1；

put “error：Power′s range：0-100”；

end；

/*如有錯誤，跳出循環*/

if(error=1)then stop；

power=0；n=2；

do while(power<&power)；

f=FINV(1-alpha，n-1，n*(m-1))；

c=(1+m*r0/(1-r0))/(1+m*r1/(1-r1))；

power=100*(1-PROBF(c*f，n-1，n*(m-1)))；

if(power>=&power)then leave；

else n=n+1；

end；

power=round(power，0.01)；

run；/*結果輸出*/

proc print data=a label；

var alpha r0 r1 m power n；

label

alpha=“Test Significance Level”

r0=“Null Intracluster Correlation”

r1=“Alternative Intracluster Correlation”

m=“Number of Measurements/Raters”

power=“Power(%)”

n=“Sample Size”；

quit；

%mend ICC；

%ICC(alpha=0.025，r0=0.75，r1=0.85，m=3，power=90)

SAS運行結果：

圖4 SAS 9.4 關于例2樣本量估計的參數設置與計算結果

2.兩獨立樣本相關性分析

(1)差異性檢驗

①兩獨立樣本Pearson相關系數檢驗

方法：Zar[5]給出了兩獨立樣本Pearson相關系數檢驗的樣本量估計方法。基于大樣本近似正態分布理論，檢驗效能的計算公式如下：

(3)

計算樣本量時，設定樣本量n1的初始值為4，n2=n1/R，R為兩組樣本量比值，通過不斷增加樣本量直至檢驗效能滿足設定條件為止，最終得到的n1和n2分別表示兩組所需樣本量。

【例3】某研究欲比較不同組織細胞的端粒DNA長度與年齡的相關性。據既往研究，外周血白細胞中端粒DNA長度與年齡的相關性系數為0.79，心肌細胞中端粒DNA長度與年齡的相關性系數為0.87。兩相關系數分別從兩組人群獲得，試估計雙側檢驗水準α=0.05，檢驗效能為80%，兩組樣本量比例R=1的情況下，能夠發現這兩種不同細胞的端粒DNA長度與年齡的相關系數存在差異所需的樣本量。

nQuery Advanced 8.6.0.0實現：設置檢驗水準α=0.05，雙側檢驗，檢驗效能1-β=80%，其他數據相應代入。在nQueryAdvanced 8.6.0.0主菜單選擇：

方法框中選擇：Two Correlations

在彈出的樣本量計算窗口將各參數值鍵入，如圖5所示，結果n1=233，n2=233。即本研究兩組所需樣本量各為233例，共需466例。

圖5 nQueryAdvanced8.6.0.0關于例3樣本量估計的參數設置與計算結果

SAS 9.4軟件實現：

/*兩獨立樣本相關系數的檢驗*/

%macro twocor(

alpha=/*檢驗水準*/

，side=/*單雙側檢驗，1表示單側檢驗，2表示雙側檢驗*/

，p1=/*樣本1的相關系數*/

，p2=/*樣本2的相關系數*/

，R=/*兩組樣本量比例*/

，power=/*檢驗效能(%)*/

)；

data a；

p1=&p1；p2=&p2；alpha=&alpha；R=&R；side=&side；

/*輸出錯誤信息*/

if(alpha>0.2 | alpha<0)then do；

error=1；

put “error：Test significance level′s range：0-0.2”；

end；

if(side^=1 & side^=2)then do；

error=1；

put “error：Side′s range：1 OR 2”；

end；

if(p1<=-1 | p1>=1)then do；

error=1；

put “error：P′s range：-1

end；

if(p2<=-1 | p2>=1)then do；

error=1；

put “error：P′s range：-1

end；

if(&power>100 | &power<0)then do；

error=1；

put “error：Power′s range：0-100”；

end；

/*如有錯誤，跳出循環*/

if(error=1)then stop；

z1=0.5*log((1+p1)/(1-p1))；

z2=0.5*log((1+p2)/(1-p2))；

power=0；

n1=4；

/*直到效能達到要求，跳出循環*/

do while(power<&power)；

sigma=sqrt(1/(n1-3)+1/(n1/R-3))；

/*雙側*/

if side=2 then

power=100*(probnorm((z1-z2)/sigma-probit(1-alpha/2))+

probnorm(-((z1-z2)/sigma)-probit(1-alpha/2)))；

/*p1

if side=1 and p1

power=100*(probnorm((z1-z2)/sigma-probit(1-alpha)))；

/*p1>p2，單側*/

if side=1 and p1>p2 then

power=100*(probnorm(-((z1-z2)/sigma)-

probit(1-alpha)))；

if(power >=&power)then leave；

else n1=n1+1；

n2=n1/R；

end；

n2=ceil(n1/R)；R=n1/n2；

power=round(power，0.01)；

run；

/*結果輸出*/

proc print data=a label；

var alpha side p1 p2 R n1 n2 power；

label

alpha=“Test significance level”

side=“1 or 2 sided test”

p1=“Control Correlation”

p2=“Treatment Correlation”

R=“Sample Size Ratio”

n1=“Control Sample Size”

n2=“Treatment Sample Size”

power=“Power(%)”；

quit；

%mend twocor；

%twocor(alpha=0.05，side=2，p1=0.79，p2=0.87，R=1，power=80)

SAS運行結果：

圖6 SAS 9.4關于例3樣本量估計的參數設置和計算結果