用pwr包做功效分析

R語言中的pwr包可以實現很好的功效分析，包中包含了一系列非常重要的函數。對於每個函數，用戶可以設定四個量（樣本大小、顯著性水平、功效和效應值）中的三個量，第四個量將由軟件計算出來。

四個量中，效應值是最難規定的，計算效應值通常需要一些相關估計的經驗和對過去研究知識的理解。但是如果在一個特定的研究中，你對需要的效應值一無所知，又該怎麼辦呢？接下來將介紹pwr包在常見統計檢驗（t檢驗、方差分析、相關性、線性模型、比例檢驗、卡方檢驗）中的應用。在調用以上函數時，請確保已經安裝並載入pwr包。

t檢驗的功效分析

對於t檢驗，pwr.t.test()函數提供了許多有用的功效分析選項，格式為：

pwr.t.test(n=, d=, sig.level=, power=, alternative=)

其中，

1、n為樣本大小；

2、d為效應值，即標準化的均值之差；

3、sig.level表示顯著性水平（默認為0.05）；

4、power為功效水平；

5、type表示檢驗類型：雙樣本t檢驗（two.sample）、單樣本t檢驗（one.sample）或相依樣本t檢驗（paired）。默認為雙樣本t檢驗。

6、alternative表示統計檢驗時雙側檢驗（two.sided）還是單側檢驗（less或greater）。默認為雙側檢驗。

例子：我們想評價使用手機對駕駛員反應時間的影響，零為假設H0：u1-u2=0，u1是駕駛員使用手機時的反應時間均值，u2是駕駛員不使用手機時的反應時間均值。假設我們拒絕該零假設，備擇假設就是H1：u1-u2 != 0，即兩種條件下反應時間的均值不相等。現挑選一個由不同個體構成的樣本，將他們隨機分配到任意一種情況下，參與者邊打手機，邊在一個模擬器中應對一系列駕駛挑戰；第二種情況，參與者在一個模擬器中完成一系列相同的駕駛挑戰，但不打手機。然後評估每個個體的總體反應時間。

假定將使用雙尾獨立樣本t檢驗來比較兩種情況下駕駛員的反應時間均值。如果根據過去的經驗知道反應時間有1.25s的標準偏差，並認定反應時間1s的差值是巨大的差異，那麼在這個研究中，可設定要檢測的效應值為d=1/1.25=0.8或者更大。另外，如果差異存在，我們希望有90%的把握檢測到它，由於隨機變異性的存在，我們也希望有95%的把握不會誤報差異顯著。這時，對於該研究需要多少受試者呢？

> library(pwr)

> pwr.t.test(d=.8, sig.level = .05, power = .9, type = "two.sample", alternative = "two.side")

Two-sample t test power calculation

n = 33.82555

d = 0.8

sig.level = 0.05

power = 0.9

alternative = two.sided

NOTE: n is number in *each* group

結果表明，每組中需要34個受試者（總共68人），這樣才能保證有90%的把握檢測到0.8的效應值，並且最多5%的可能性會誤報差異存在。

現在變化一下問題。假定在比較這兩種情況時，我們想檢測到總體均值0.5個標準誤差的差異，並且將誤報差異的機率限制在1%內。此外，受試者只有40人，那麼在該研究中能檢測到這麼大總體均值差異的概率是多少呢？

假定每種情況下受試者數目相同，可以進行如下操作：

> pwr.t.test(n=20, d=.5, sig.level = .01, type = "two.sample", alternative = "two.side")

Two-sample t test power calculation

n = 20

d = 0.5

sig.level = 0.01

power = 0.1439551

alternative = two.sided

NOTE: n is number in *each* group

結果表明：在0.01的先驗顯著性水平下，每組20個受試者，因變量的標準差為1.25s，有低於14%的可能性斷言差值為0.625s或者不顯著（d=0.5=0.625/1.25）。換句話說，我們有86%的可能錯過我們要尋找的效應值。因此，可能需要慎重考慮要投入到該研究中的時間和經歷。

上面的例子都是假定兩組樣本大小相等，如果兩組中樣本大小不同，可用函數：

pwr.t2n.test(n1=, n2=, d=, sig.level=, power=, alternative=)

其中n1和n2是兩組樣本大小，其他參數含義與pwr.t.test()相同，可以嘗試改變pwr.t2n.test()函數的參數值，看看輸出的效應值如何變化。