R语言实战__第10章 功效分析

第10章 功效分析

主要内容：

• 判断所需样本量
• 计算效应值
• 评价统计功效

功效分析可以帮助在给定置信度的情况下，判断达到给定效应值时所需的样本量。也可在给定置信度水平情况下，计算在某样本量内能检测到给定效应值的概率。

10.1 假设检验速览

统计假设检验中，首先对总体分布参数设定一个假设（零假设$H_0$），然后从总体分布中抽样，通过样本计算所得统计量来对总体参数进行推断。假定零假设为真，如果计算获得观测样本的统计量的概率非常小，便可以拒绝厡假设，接收他的对立面（称为备择假设或者研究假设$H_1$）。

由于使用取自总体的样本数据来对总体做推断。判断可能有一下四种结果。

假设检验结果.png

在研究过程中，关注四个量：

• 样本大小：实验设计中每种条件/组中观测的数目；
• 显著性水平：由Ⅰ型错误的概率来定义。可以看作发现效应不发生的概率；
• 功效：1-P（Ⅱ型错误）。可以看作真实效应发生的概率；
• 效应值：在备择假设或研究假设下效应的量，依赖于假设检验中使用的统计方法。

可以直接控制样本大小和显著性水平，但对功效和效应值的影响是间接的。

研究目标是维持一个可接受的显著性水平，尽量使用较少的样本，然后最大化统计检验的功效。也就是说，最大化发现真是效应的几率，并最小化发现错误效应的几率，同时把研究成本控制在合理的范围内。

四个量（样本大小、显著性水平、功效和效应值）紧密相关，给定其中任意三个量，便可推算第四个量。

10.2 用pwr包做功效分析

pwr包可实现功效分析。表中列出一些重要函数。对于每个函数，用户可设定四个量（样本大小、显著形式水平、功效和效应值）中的三个量，第四个量由软件计算。

pwr包中的函数.png

四个量中，效应值最难规定。

10.2.1 t检验

对于t检验，pwr.t.test()函数提供了许多有用的功效分析选项，格式为：
pwr.t.test(n = , d = , sig.level = , power = , alternative = )
其中：

• n为样本大小
• d为效应值，即标准化的均值之差$$d = \frac {\mu_1 - \mu_2} {\sigma}$$ $\mu1=$组1均值， $\mu2=$组2均值，$\sigma^2=$误差方差
• sig.level为显著性水平（默认0.05）
• power为功效水平
• type为检验类型：双样本t检验（two.sample）、单样本t检验（one.sample）、相依样本t检验（paired）。默认双样本t检验
• alternative指统计检验是双侧检验（two.sided）还是单侧检验（less或greater）。默认双侧检验。

> library(pwr)
> pwr.t.test(d= 0.8, sig.level=0.5, power=0.9, type="two.sample", alternative="two.sided")

     Two-sample t test power calculation 

              n = 11.76996
              d = 0.8
      sig.level = 0.5
          power = 0.9
    alternative = two.sided

NOTE: n is number in *each* group

> pwr.t.test(n=20, d=0.5, sig.level=0.01, type="two.sample", alternative="two.sided")

     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
> pwr.t2n.test(

10.2.2 方差分析

pwr.anova.test()函数可以对平衡单因素分析进行功效分析。格式为：
pwr.anova.test(k= , n= , f= , sig.level= , power= )
其中，k是组的个数，n是各组中的样本大小。
对于单因素方差分析，效应值可通过f来衡量：
$$ f = \sqrt { \frac {\sum_{i-1}^{k} {p_i × (\mu_i - \mu)^2}} {\sigma^2} } $$
其中：
$ p_i = n_i / N$
$ n_i = 组i的观测数目 $
$ N = 总观测数目 $
$ \mu_i = 组i均值 $
$ \mu = 总体均值 $
$ \sigma^2 = 组内误差方差 $

例：5个组做单因素方差分析，要达到功效0.8，效应值0.25，显著性水平0.05，计算各组所需样本大小：

> #方差分析
> pwr.anova.test(k=5, f=0.25, sig.level=.05, power=.8)

     Balanced one-way analysis of variance power calculation 

              k = 5
              n = 39.1534
              f = 0.25
      sig.level = 0.05
          power = 0.8

NOTE: n is number in each group

> k：分组个数，n：组内样本大小，f：效应值（根据统计方法计算），sig.level：显著性水平，power：功效

本例中需要估计在同方差时五个组的均值。以计算$f$如果不知道，参见10.2.7节。

10.2.3 相关性

pwr.r.test()函数可对相关性分析进行功效分析。格式如下：
pwr.r.test(n= , r= , sig.level= , power= , alternative= )
其中，n是观测数目，r是效应值（通过相关系数衡量），alternative指定显著性检验是双边检验（two.sided）还是单边检验（less或greater）。

> pwr.r.test(r=0.25, sig.level=0.05, power=0.90, alternative="greater")

     approximate correlation power calculation (arctangh transformation) 

              n = 133.2803
              r = 0.25
      sig.level = 0.05
          power = 0.9
    alternative = greater

10.2.4 线性模型

pwr.f2.test()函数可进行线性模型的功效分析。格式为：
pwr.f2.test(u= , v= , f2= , sig.level= , power= )
其中，$u$和$v$分别是分子自由度和分母自由度，f2是效应值。
$$f^2 = \frac {R^2} {1-R^2}$$
$$ f^2 = \frac {{R}{AB}^{2} - {R} {A}^{2}} {1-{R}{AB}^{2}} $$
其中：
$$R^2=多重相关性的总体平方值$$
$$R_A^2 = 集合A中变量对总体方差的贡献率$$
$$R{AB}^2= 集合A和B中变量对总体方差的贡献率$$

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

> 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=193。

10.2.5 比例检验

pwr.2p.test()函数对比例进行功效分析。格式为：
pwr.2p.test(h= , n= ,sig.level= , power= )
其中：h是效应值，n是各组相同的样本量。效应值h定义如下：
$$h= 2 arcsin(\sqrt{P_1})-2arcsin({\sqrt{P_2}})$$
可用h=ES.h(p1, p2)函数进行计算。
当个组中n不同时，使用函数：
pwr.2p2n.test(h= , n1= , n2= , sig.level= , power= )
alternative=选项可设定检验是双尾检验（two.sided）还是单尾检验（less或greater）。默认双尾检验。

> pwr.2p.test(h=ES.h(0.65,0.6), sig.level=0.05, power=0.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

10.2.6 卡方检验

卡方检验常用来评价两个类别型变量的关系。典型的零假设是变量之间独立，备择假设是不独立。pwr.chisq.test()函数可以评估卡方检验的功效、效应值和所需的样本大小。格式为：
pwr.chisq.test(w= , N= , df= , sig.level= , power= )
其中，w是效应值，N是总样本大小，df是自由度。此外，效应值w如下定义：
$$ w = \sqrt { \frac {\sum_{i-1}^{m} {(p 0_i - p1_i)^2}} {p 0_i} }$$
其中：
$p0_i = H_0 时第$i$单元格中的概率$
$p1_i = H_1 时第$i$单元格中的概率$
此处从1到m进行求和，累加符号上m指列联表中单元格的数目。函数ES.w2(P)可以计算双因素列联表中备择假设的效应值，p是一个假设的双因素概率表。

> prob <- matrix(c(.42, .28, .03, .07, .10, .10), byrow=TRUE, nrow=3)
> pwr.chisq.test(w = ES.w2(prob), df=2, sig.level=.05, power=.9)

     Chi squared power calculation 

              w = 0.1853198
              N = 368.4528
             df = 2
      sig.level = 0.05
          power = 0.9

NOTE: N is the number of observations

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

功效分析中，预期效应值较难确定。需要对主题有了解并结合经验。

行为科学领域，Cohen提出一个基准，可为各种统计检验划分“小”、“中”、“大”三种效应值。

Cohen效应值基准.png

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

> library(pwr)
> es <- seq(.1, .5, .01)
> nes <- length(es)
> nes
[1] 41
> #seq生成向量，开始，结束，间距
> 
> samize <- 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")

效应值与样本数量.png

图形可以看出，随着样本量超过200，再增加样本效应值的减小已经不明显了。

10.3 绘制功效分析图形

假设对于相关系数统计显著性的检验，要计算一系列效应值和功效水平下所需的样本量，此时可用pwr.r.test()函数和for循环来完成任务。
代码清单10-2 检验各种效应值下的相关性所需的样本量曲线

> library(pwr)
> r <- seq(.1, .5, .01)
> nr <- length(r)
> 
> p <- seq(.4, .9, .1)
> np <- length(p)
> 
> samsize <- array(numeric(nr*np), dim=c(nr, np))
> for (i in 1:np) {
+  for (j in 1:nr) {
+    result <- pwr.r.test(n=NULL, r=r[j],
+    sig.level=.05, power=p[i],
+    alternative="two.sided")
+    samsize[j,i] <- ceiling(result$n)
+  }
+ }
> xrange <- range(r)
> yrange <- round(range(samsize))
> colors <- rainbow(length(p))
> plot(xrange, yrange, type="n", 
+      xlab="Correlation Coefficient (r)",
+      ylab="Sample Size (n)")
> 
> for (i in 1:np) {
+   lines(r, samsize[ , i], type="l", lwd=2, col=colors[i])
+ }
> 
> abline(v=0, h=seq(0, yrange[2], 50), lty=2, col="grey89")
> abline(h = 0, v=seq(xrange[1], xrange[2], .02), lty=2, col="gray89")
> title("Sample Size Estimation for Corelation Studies\n
+       Sig = .05 (Two-tailed)")
> legend("topright", title="Power", as.character(p), fill=colors)

不同效应值下各相关性所需样本量.png

10.4 其他软件包

piface包提供了一个与R交互的Java图形用户截面（GUI），包含各种计算样本量的方法。该GUI允许用户交互性地修改研究参数，观察对其他参数的影响。

ps：此包不可用，暂时没找到替代下载方式。

专业化的功效分析软件包.png

10.5 小结

功效分析不仅可以帮助判断在给定置信度和效应值的前提下所需的样本量，也能说明在给定样本量时检测到要求效应值的概率。

功效分析的一个重要附加效益是引起方向性的转变，它鼓励不要仅仅关注于二值型（即效应存在与否）的假设检验，而应该思考效应值增加的意义。研究结果要既包含p值又包含效应值。它们不仅帮助判断研究的实际意义，还提供用于未来研究的信息。

---

——2017.3.9
——于510

ps：工作室网络抽风两天，整个人都不好了。