​R语言实战(10)——功效分析

引言:在前两章的学习中我们了解了如何用R来实现回归和方差分析这两种常用的统计方法,在实验设计中,如何确定合适的样本量呢?本章我们将一起来学习用R进行功效分析,来解决这个问题。  

功效分析

  • 可以在给定置信度和效应值的情况下,判断所需的样本量。
  • 可以在给定置信度水平和样本量的情况下,计算检测到要求效应值的概率。

10.1 用pwr包做功效分析

表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不相等的两样本)

对于每个函数,用户可以设定四个量(样本大小、显著性水平、功效和效应值)中的三个量,计算出第四个量。

1. t 检验

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

pwr.t.test(n=, d=, sig.level=, power=, type=, alternative= )
  1. n为样本大小
  2. d为效应值,即标准化的均值之差
    d=(μ1-μ2)/σ,μ1=组1均值,μ2=组2均值,σ^2 =误差方差
  3. sig.level表示显著性水平(默认为0.05)
  4. power为功效水平
  5. type指检验类型:双样本t检验(“two.sample”)、单样本t检验(“one.sample”)或相依样本t检验(“paired”),默认为双样本t检验
  6. “alternative”指统计检验是双侧检验(“two.sided”)还是单侧检验(“less”或”greater”),默认为双侧检验

例:使用双尾独立样本t检验来比较有无使用手机的驾驶员的反应时间均值。
根据经验:反应时间有1.25s的标准偏差,认定反应时间1s的差值是巨大的差异。

问题

  • 设定要检测的效应值为d=1/1.25=0.8或者更大

  • 如果差异存在,希望有90%的把握检测到它:功效水平90%

  • 由于随机变异性的存在,希望有95%的把握不会误报差异显著:显著性水平0.05

计算所需样本量n:

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 = 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个标准偏差的差异:效应值d=0.5

  • 将误报差异的几率限制在1%内:显著性水平0.01

  • 能获得的受试者只有40人:n=20

计算检测到这么大总体均值差异的概率:功效水平:

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.1439551    alternative = two.sided
NOTE: n is number in *each* group

因变量的标准差为1.25s,有低于14%的可能性断言差值为0.625s或者不显著(d=0.5=0.625/1.25)。即,有86%的可能性错过要寻找的效应值。因此,可能需要慎重考虑要投入到该研究中的时间和精力。

若两组中样本大小不同,可用:

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

此时可通过设定效应值d计算功效水平power,或设定功效水平power计算效应值d。

示例:

pwr.t2n.test(n1=18, n2=22, power=.9, alternative="two.sided")     test power calculation                                     n1 = 18                          n2 = 22                            d = 1.057295            sig.level = 0.05       #默认值为0.05                    power = 0.9        alternative = two.sided

2. 方差分析

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

pwr.anova.test(k=, n=, f=, sig.level=, power= )
  1. k是组的个数
  2. n是各组中的样本大小

单因素方差分析,效应值可通过f来衡量:

pi=ni/N
ni = 组i的观测数目
N = 总观测数目
μi = 组i均值
μ = 总体均值
σ^2 = 组内误差方差

例:对五个组做单因素方差分析,要达到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.1534            f = 0.25     sig.level = 0.05        power = 0.8
   NOTE: n is number in each group

总样本大小为5×39,即195。

3. 相关性

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

pwr.r.test(n=, r=, sig.level=, power=, alternative= )
  1. n是观测数目
  2. r是效应值(通过线性相关系数衡量)
  3. sig.level是显著性水平
  4. power是功效水平
  5. alternative指定显著性检验是双边检验(“tow.sided”)还是单边检验(“less” 或”greater”)

例:研究抑郁与孤独的关系
H0: ρ≤0.25 和 H1: ρ > 0.25;ρ是两个心理变量的总体相关性大小。

  • 设定显著性水平为0.05

  • 如果H0是错误的,想有90%的信心拒绝H0:power=0.9

计算所需样本量:

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

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

4. 线性模型

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

pwr.f2.test(u=, v=, f2=, sig.level=, power= )
  1. u是分子自由度
  2. v是分母自由度
  3. f2是效应值

R^2=多重相关性的总体平方值
RA^2=集合A中变量对总体方差的解释率
RAB^2=集合A和B中变量对总体方差的解释率

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

例:研究老板的领导风格对员工满意度的影响,是否超过薪水和工作小费对员工满意度的影响
领导风格可用四个变量来评估,薪水和小费与三个变量有关;经验表明,薪水和小费能够解释约30%的员工满意度的方差;从现实出发,领导风格至少能解释35%的方差。

  • 假定显著性水平为0.05:sig.level=0.05

  • 在90%的置信度情况:power=0.90

  • u=3(总预测变量数减去集合B中的预测变量数)

  • 效应值为f2=(0.35–0.30)/(1–0.35) = 0.076

计算需要多少受试者才能得到这样的方差贡献率:

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。

5. 比例检验

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

pwr.2p.test(h=, n=, sig.level=, power=)
  1. h是效应值
  2. n是各组相同的样本量
  3. alternative=选项可以设定检验是双尾检验(“two.sided”)还是单尾检验(“less”或 “greater”),默认是双尾检验

h=2arcsin(√p1)-2arcsin(√p2),可用ES.h(p1, p2)函数进行计算。
各组中n不相同时,则使用函数:

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

例:假定对某流行药物能缓解60%使用者的症状感到怀疑,一种更贵的新药如果能缓解65%使用者的症状,就会被投放到市场中。

  • 假设想有90%的把握得出新药更有效的结论:power=0.9

  • 希望有95%的把握不会误得结论:显著性水平0.05

  • 只对评价新药是否比标准药物更好感兴趣,因此只需用单边检验:alternative=”greater”

计算能够检测到两种药物存在这一特定的差异所需的样本量:

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个人试用已有药物。

6. 卡方检验

卡方检验常常用来评价两个类别型变量的关系。

pwr.chisq.test()函数可以评估卡方检验的功效、效应值和所需的样本大小,格式为:

pwr.chisq.test(w=, N=, df=, sig.level=, power= )
  1. w是效应值
  2. N是总样本大小
  3. df是自由度

p0i=H0时第i单元格中的概率
p1i=H1时第i单元格中的概率 此处从1到m进行求和,m是列联表中单元格的数目

函数ES.w2(P)可以计算双因素列联表中备择假设的效应值,P是一个假设的双因素概率表。

例:研究人种与工作晋升的关系

预期样本中70%是白种人,10%是美国黑人,20%是西班牙裔人;相比30%的美国黑人和50%的西班牙裔人,60%的白种人更容易晋升。

表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)                                              [10.1853198

计算所需的样本大小:

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

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

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

功效分析中,预期效应值是最难决定的参数,通常需要对主题有一定的了解,并有相应的测量经验。
Cohen效应值基准可为各种统计检验划分“小”“中”“大”三种效应值。
表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

Cohen的基准值根据许多社科类研究得出的一般性建议,对于特殊的研究领域可能并不适用。

其他可选择的方法是改变研究参数,记录其对诸如样本大小和功效等方面的影响。

例:五个分组的单因素方差分析(显著性水平为0.05)

计算单因素ANOVA中检测显著效应所需的样本大小:

library(pwr)es <- seq(.1, .5, .01)nes <- length(es)samsize <- NULLfor (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-1 五分组的单因素ANOVA中检测显著效应所需的样本大小

图形有助于估计不同条件时的影响值。各组样本量高于200个观测时,再增加样本效果不大

10.2 绘制功效分析图形

对于相关系数统计显著性的检验, 可用pwr.r.test()函数和for循环来完成任务,计算一系列效应值和功效水平下所需的样本量。

检验各种效应值下的相关性所需的样本量曲线

library(pwr)#seq函数生成一系列效应值r(H1时的相关系数)和功效水平pr <- seq(.1,.5,.01)nr <- length(r)p <- seq(.4,.9,.1)np <- length(p)#获取样本大小,将其存储在数组samsize中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)" )#使用曲线形式(lines)而不是点形式(points)来添加功效曲线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 Correlation Studiesn      Sig=0.05 (Two-tailed)")legend("topright", title="Power", as.character(p),       fill=colors)

图10-2 在不同功效水平下检测到显著的相关性所需的样本量

在40%的置信度下,要检测到0.20的相关性,需要约75的样本量;
在90%的置信度下,要检测到相同的相关性,需要大约185个额外的观测(n=260)。

10.3 其他软件包

表10-4 专业化的功效分析软件包

软件包 目的
asypow 通过渐进似然比方法计算功效
longpower 纵向数据中的样本量的计算
PwrGSD 组序列设计的功效分析
pamm 混合模型中随机效应的功效分析
powerSurvEpi 流行病研究的生存分析中功效和样本量的计算
powerMediation 线性、Logistic、泊松和 Cox 回归的中介效应中功效和样本量的计算
powerpkg 患病同胞配对法和TDT(Transmission Disequilibrium Test,传送不均衡检验)设计的功效分析
powerGWASinteraction GWAS 交互作用的功效计算
pedantics 一些有助于种群基因研究功效分析的函数
gap 一些病例队列研究设计中计算功效和样本量的函数
ssize.fdr 微阵列实验中样本量的计算

10.4 小结

本章我们主要学习了使用pwr包中函数对常见的统计方法(包括t检验、卡方检验、比例检验、ANOVA和回归)进行功效和样本量的计算。在之前的学习中我们了解了一些基本图形的可视化,下一章我们将学习如何将多元关系可视化的方法。敬请期待~

生物信息学

R语言实战(9)——方差分析

2019-12-17 21:42:49

生物信息学

R语言实战(11)——中级绘图

2019-12-17 21:46:03

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索