10.2　用pwr包做功效分析

10.2　用pwr包做功效分析

Stéphane Champely开发的pwr包可以实现Cohen（1988）描述的功效分析。表10-1列出了一些非常重要的函数。对于每个函数，用户可以设定四个量（样本大小、显著性水平、功效和效应值）中的三个量，第四个量将由软件计算出来。 表10-1　pwr包中的函数

函　　数	功效计算的对象
`pwr.2p.test()`	两比例（n相等）
`pwr.2p2n.test()`	两比例（n不相等）
`pwr.anova.test()`	平衡的单因素ANOVA
`pwr.chisq.test()`	卡方检验
`pwr.f2.test()`	广义线性模型
`pwr.p.test()`	比例（单样本）
`pwr.r.test()`	相关系数
`pwr.t.test()`	t检验（单样本、两样本、配对）
`pwr.t2n.test()`	t检验（n不相等的两样本）

四个量中，效应值是最难规定的。计算效应值通常需要一些相关估计的经验和对过去研究知识的理解。但是如果在一个特定的研究中，你对需要的效应值一无所知，该怎么做呢？10.2.7节将会讨论这个难题。本节接下来介绍pwr包在常见统计检验中的应用。在调用以上函数时，请确定已经安装并载入pwr包。

10.2.1　t检验

对于t检验，pwr.t.test()函数提供了许多有用的功效分析选项，格式为：

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

其中元素解释如下。

n为样本大小。
d为效应值，即标准化的均值之差。

10.2　用pwr包做功效分析 - 图1

其中

μ1= 组1均值

μ2= 组2均值

σ2 = 误差方差

sig.level表示显著性水平（默认为0.05）。
power为功效水平。
type指检验类型：双样本t检验（two.sample）、单样本t检验（one.sample）或相依样本t检验（paired）。默认为双样本t检验。
alternative指统计检验是双侧检验（two.sided）还是单侧检验（less或greater）。默认为双侧检验。

让我们举例说明函数的用法。仍继续10.1节使用手机与驾驶反应时间的实验，假定将使用双尾独立样本t检验来比较两种情况下驾驶员的反应时间均值。

如果你根据过去的经验知道反应时间有1.25 s的标准偏差，并认定反应时间1 s的差值是巨大的差异，那么在这个研究中，可设定要检测的效应值为d=1/1.25=0.8或者更大。另外，如果差异存在，你希望有90%的把握检测到它，由于随机变异性的存在，你也希望有95%的把握不会误报差异显著。这时，对于该研究需要多少受试者呢？

将这些信息输入到pwr.t.test()函数中，形式如下：

> library(pwr)
> pwr.t.test(d=.8, sig.level=.05, power=.9, type="two.sample",
  alternative="two.sided")
    Two-sample t test power calculation
             n = 34
             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.sided")
     Two-sample t test power calculation
              n = 20
              d = 0.5
      sig.level = 0.01
          power = 0.14
    alternative = two.sided
 NOTE: n is number in *each* group

结果表明，在0.01的先验显著性水平下，每组20个受试者，因变量的标准差为1.25 s，有低于14%的可能性断言差值为0.625 s或者不显著（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()¹函数中的参数值，看看输出的效应值如何变化。

^{1 R中函数名称后面最好加上()。——译者注}

10.2.2　方差分析

pwr.anova.test()函数可以对平衡单因素方差分析进行功效分析。格式为：

pwr.anova.test(k=, n=, f=, sig.level=, power=)

其中，k是组的个数，n是各组中的样本大小。

对于单因素方差分析，效应值可通过f来衡量：

10.2　用pwr包做功效分析 - 图2

其中，p_i = n_i/ N，

n_i = 组i的观测数目

N = 总观测数目

μ_i = 组i均值

μ= 总体均值

σ² = 组内误差方差

让我们举例说明函数用法。现对五个组做单因素方差分析，要达到0.8的功效，效应值为0.25，并选择0.05的显著性水平，计算各组需要的样本大小。代码如下：

> pwr.anova.test(k=5, f=.25, sig.level=.05, power=.8)
     Balanced one-way analysis of variance power calculation
              k = 5
              n = 39
              f = 0.25
      sig.level = 0.05
          power = 0.8
 NOTE: n is number in each group

结果表明，总样本大小为5 × 39，即195。注意，本例中需要估计在同方差时五个组的均值。如果你对上述情况都一无所知，10.2.7节提供的方法可能会有所帮助。

10.2.3　相关性

pwr.r.test()函数可以对相关性分析进行功效分析。格式如下：

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

其中，n是观测数目，r是效应值（通过线性相关系数衡量），sig.level是显著性水平，power是功效水平，alternative指定显著性检验是双边检验（tow.sided）还是单边检验（less或greater）。

假定正在研究抑郁与孤独的关系。你的零假设和研究假设为：

H₀:ρ ≤ 0.25 和 H₁:ρ > 0.25

其中，ρ是两个心理变量的总体相关性大小。你设定显著性水平为0.05，而且如果H₀是错误的，你想有90%的信心拒绝H₀，那么研究需要多少观测呢？下面的代码给出了答案：

> pwr.r.test(r=.25, sig.level=.05, power=.90, alternative="greater")
     approximate correlation power calculation (arctangh transformation)
              n = 134
              r = 0.25
      sig.level = 0.05
          power = 0.9
    alternative = greater

因此，要满足以上要求，你需要134个受试者来评价抑郁与孤独的关系，以便在零假设为假的情况下有90%的信心拒绝它。

10.2.4　线性模型

对于线性模型（比如多元回归），pwr.f2.test()函数可以完成相应的功效分析，格式为：

pwr.f2.test(u=, v=, f2=, sig.level=, power=)

其中，u和v分别是分子自由度和分母自由度，f2是效应值。

10.2　用pwr包做功效分析 - 图3

其中

R² = 多重相关性的总体平方值

10.2　用pwr包做功效分析 - 图4

其中

R_A² = 集合A中变量对总体方差的解释率²

^{2 也常称作方差贡献率。——译者注}

R_AB² = 集合A和 B中变量对总体方差的解释率

当要评价一组预测变量对结果的影响程度时，适宜用第一个公式来计算f2；当要评价一组预测变量对结果的影响超过第二组变量（协变量）多少时，适宜用第二个公式。

现假设你想研究老板的领导风格对员工满意度的影响，是否超过薪水和工作小费对员工满意度的影响。领导风格可用四个变量来评估，薪水和小费与三个变量有关。过去的经验表明，薪水和小费能够解释约30%的员工满意度的方差。而从现实出发，领导风格至少能解释35%的方差。假定显著性水平为0.05，那么在90%的置信度情况下，你需要多少受试者才能得到这样的方差贡献率呢？

此处，sig.level = 0.05，power = 0.90，u = 3（总预测变量数减去集合B中的预测变量数），效应值为f2 = (0.35 - 0.30)/(1 - 0.35) = 0.076 9。将这些信息输入到函数中：

> pwr.f2.test(u=3, f2=0.0769, sig.level=0.05, power=0.90)
     Multiple regression power calculation
              u = 3
              v = 184.2426
             f2 = 0.0769
      sig.level = 0.05
          power = 0.9

在多元回归中，分母的自由度等于N - k - 1，N是总观测数，k是预测变量数。本例中，N - 7 - 1 = 185，即需要样本大小N = 185 + 7 + 1 = 193。

10.2.5　比例检验

当比较两个比例时，可使用pwr.2p.test()函数进行功效分析。格式为：

pwr.2p.test(h=, n=, sig.level=, power=)

其中，h是效应值，n是各组相同的样本量。效应值h定义如下：

10.2　用pwr包做功效分析 - 图5

可用ES.h(p1, p2)函数进行计算。

当各组中n不相同时，则使用函数：

pwr.2p2n.test(h =, n1 =, n2 =, sig.level=, power=).

alternative =选项可以设定检验是双尾检验（two.sided）还是单尾检验（less或greater）。默认是双尾检验。

假定你对某流行药物能缓解60%使用者的症状感到怀疑。而一种更贵的新药如果能缓解65%使用者的症状，就会被投放到市场中。此时，在研究中你需要多少受试者才能够检测到两种药物存在这一特定的差异？

假设你想有90%的把握得出新药更有效的结论，并且希望有95%的把握不会误得结论。另外，你只对评价新药是否比标准药物更好感兴趣，因此只需用单边检验，代码如下：

> pwr.2p.test(h=ES.h(.65, .6), sig.level=.05, power=.9,
              alternative="greater")
     Difference of proportion power calculation for binomial
     distribution (arcsine transformation)
              h = 0.1033347
              n = 1604.007
      sig.level = 0.05
          power = 0.9
    alternative = greater
 NOTE: same sample sizes

根据结果可知，为满足以上要求，在本研究中需要1605个人试用新药，1605个人试用已有药物。

10.2.6　卡方检验

卡方检验常常用来评价两个类别型变量的关系。典型的零假设是变量之间独立，备择假设是不独立。pwr.chisq.test()函数可以评估卡方检验的功效、效应值和所需的样本大小。格式为：

pwr.chisq.test(w =, N = , df = , sig.level =, power = )

其中，w是效应值，N是总样本大小，df是自由度。此处，效应值w如下定义：

10.2　用pwr包做功效分析 - 图6

其中

p 0_i = H₀时第i单元格中的概率

p 1_i = H₁时第i单元格中的概率

此处从1到m进行求和，连加号上的m指的是列联表中单元格的数目。函数ES.w2(P)可以计算双因素列联表中备择假设的效应值，P是一个假设的双因素概率表。

举一个简单的例子，假设你想研究人种与工作晋升的关系。你预期样本中70%是白种人，10%是美国黑人，20%西班牙裔人。而且，你认为相比30%的美国黑人和50%的西班牙裔人，60%的白种人更容易晋升。研究假设的晋升概率如表10-2所示。 表10-2　研究假设下预期晋升的人群比例

人　　种	晋升比例	未晋升者比例
白种人	0.42	0.28
美国黑人	0.03	0.07
西班牙裔	0.10	0.10

从表中看到，你预期总人数的42%是晋升的白种人（0.42 = 0.70 × 0.60），总人数的7%是未晋升的美国黑人（0.07 = 0.10 × 0.70）。让我们取0.05的显著性水平和0.90的预期功效水平。双因素列联表的自由度为(r-1)(c-1)，r是行数，c是列数。编写如下代码，你可以计算假设的效应值：

> prob <- matrix(c(.42, .28, .03, .07, .10, .10), byrow=TRUE, nrow=3)
> ES.w2(prob)
[1] 0.1853198

使用该信息，你又可以计算所需的样本大小：

> pwr.chisq.test(w=.1853, df=2, sig.level=.05, power=.9)
     Chi squared power calculation
              w = 0.1853
              N = 368.5317
             df = 2
      sig.level = 0.05
          power = 0.9
 NOTE: N is the number of observations

结果表明，在既定的效应值、功效水平和显著性水平下，该研究需要369个受试者才能检验人种与工作晋升的关系。

10.2.7　在新情况中选择合适的效应值

功效分析中，预期效应值是最难决定的参数。它通常需要你对主题有一定的了解，并有相应的测量经验。例如，过去研究中的数据可以用来计算效应值，这能为后面深层次的研究提供一些参考。

但是当面对全新的研究情况，没有任何过去的经验可借鉴时，你能做些什么呢？在行为科学领域，Cohen（1988）曾尝试提出一个基准，可为各种统计检验划分“小”、“中”、“大”三种效应值。表10-3列出了这些基准值。 表10-3　Cohen效应值基准

统计方法	效应值测量	建议的效应值基准
		小	中	大
t检验	d	0.20	0.50	0.80
方差分析	f	0.10	0.25	0.40
线性模型	f2	0.02	0.15	0.35
比例检验	h	0.20	0.50	0.80
卡方检验	w	0.10	0.30	0.50

当你对研究的效应值一无所知时，这个表可给你提供一些指引。例如，假如你想在0.05的显著性水平下，对5个组、每组25个受试者的设计进行单因素方差分析，那么拒绝错误零假设（也就是发现真实的效应值）的概率是多大呢？

使用pwr.anova.test()函数和表10-3中f的建议值，得到对于小效应值功效水平为0.118，中等效应值的为0.574，大效应值的为0.957。给定样本大小的限制，在大效应值时你才可能发现要研究的效应。

另外，你一定要牢记Cohen的基准值仅仅是根据许多社科类研究得出的一般性建议，对于特殊的研究领域可能并不适用。其他可选择的方法是改变研究参数，记录其对诸如样本大小和功效等方面的影响。仍以五个分组的单因素方差分析（显著性水平为0.05）为例，代码清单10-1计算了为检测一系列效应值所需的样本大小，结果见图10-2。

代码清单10-1　单因素ANOVA中检测显著效应所需的样本大小

library(pwr)
es <- seq(.1, .5, .01)
nes <- length(es)
samsize <- NULL
for (i in 1:nes){
    result <- pwr.anova.test(k=5, f=es[i], sig.level=.05, power=.9)
    samsize[i] <- ceiling(result$n)
}
plot(samsize,es, type="l", lwd=2, col="red",
     ylab="Effect Size",
     xlab="Sample Size (per cell)",
     main="One Way ANOVA with Power=.90 and Alpha=.05")

10.2　用pwr包做功效分析 - 图7

图10-2　五分组的单因素ANOVA中检测显著效应所需的样本大小（假定0.90的功效和0.05的显著性水平）

实验设计中，这样的图形有助于估计不同条件时的影响值。例如，从图形可以看到各组样本量高于200个观测时，再增加样本已经效果不大了。下一节我们将看看其他图形示例。

10.2 用pwr包做功效分析