R in action读书笔记(13)第十章 功效分析

时间:2024-07-09 18:07:14

功效分析

功效分析可以帮助在给定置信度的情况下,判断检测到给定效应值时所需的样本量。反过来,它也可以帮助你在给定置信度水平情况下,计算在某样本量内能检测到给定效应值的概率。如果概率低得难以接受,修改或者放弃这个实验将是一个明智的选择。

10.1假设检验速览

在研究过程时,研究者通常关注四个量:样本大小、显著性水平、功效和效应值。样本大小指的是实验设计中每种条件/组中观测的数目。显著性水平(也称为alpha)由I型错误的概率来定义。也可以把它看做是发现效应不发生的概率。功效通过1减去II型错误的概率来定义。我们可以把它看做是真实效应发生的概率。效应值指的是在备择或研究假设下效应的量。效应值的表达式依赖于假设检验中使用的统计方法。

R in action读书笔记(13)第十章 功效分析

10.2 用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.1 t 检验

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

function (n = NULL, d = NULL, sig.level = 0.05, power = NULL,

type = c("two.sample","one.sample", "paired"), alternative =c("two.sided", "less", "greater"))

n为样本大小。

d为效应值,即标准化的均值之差。d =(m1-m2)/s

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=.8,sig.level=.05,power=.9,type="two.sample",alternative="two.sided")

Two-samplet test power calculation

n =33.82554

d =0.8

sig.level =0.05

power=0.9

alternative = two.sided

NOTE: n is number in*each* group

10.2.2 方差分析

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

function (k = NULL, n = NULL, f = NULL, sig.level =0.05, power = NULL)其中,k是组的个数,n是各组中的样本大小

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

R in action读书笔记(13)第十章 功效分析

其中,pi = ni/N,

ni = 组i的观测数目

N = 总观测数目

μi

= 组i均值

μ= 总体均值

σ2 = 组内误差方差

> pwr.anova.test(k=5,f=.25,sig.level=.05,power=.8)

Balancedone-way analysis of variance power calculation

k= 5 #结果表明,总样本大小为5 × 39,即195

n= 39.1534

f= 0.25

sig.level= 0.05

power= 0.8

NOTE: n is number in each group

10.2.3 相关性

pwr.r.test()函数可以对相关性分析进行功效分析

function (n = NULL, r = NULL, sig.level = 0.05, power = NULL, 
alternative = c("two.sided", "less", "greater")) 

n是观测数目,r是效应值(通过线性相关系数衡量),sig.level是显著性水平,power是功

效水平,alternative指定显著性检验是双边检验(tow.sided)还是单边检验(less或greater)。

>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

10.2.4 线性模型

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

function (u = NULL, v = NULL, f2 = NULL, sig.level = 0.05, power = NULL) ,u和v分别是分子*度和分母*度,f2是效应值。

R in action读书笔记(13)第十章 功效分析

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

> pwr.f2.test(u=3,f2=.0769,sig.level=.05,power=.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()函数进行功效分析。格式为:

function (h = NULL, n = NULL, sig.level = 0.05, power = NULL, 
    alternative = c("two.sided", "less", "greater")) 

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

R in action读书笔记(13)第十章 功效分析

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

当各组中n不相同时,则使用函数:

function (h = NULL, n1 = NULL, n2 = NULL, sig.level = 0.05, power = NULL, alternative = c("two.sided", "less", "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

10.2.6 卡方检验

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

function (w = NULL, N = NULL, df = NULL, sig.level = 0.05, power = NULL) 

其中,w是效应值,N是总样本大小,df是*度。此处,效应值w如下定义:

R in action读书笔记(13)第十章 功效分析

p0i = H0时第i单元格中的概率

p1i = H1时第i单元格中的概率

> prob<-matrix(c(.42,.28,.03,.07,.10,.10),byrow=TRUE,nrow=3)
> ES.w2(prob)
[1] 0.1853198
> pwr.chisq.test(w=.1853198,df=2,sig.level=.05,power=.9)
     Chi squared power calculation 
              w = 0.1853198
              N = 368.4529
             df = 2
      sig.level = 0.05
          power = 0.9
NOTE: N is the number of observations

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

功效分析中,预期效应值是最难决定的参数。

R in action读书笔记(13)第十章 功效分析

单因素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")

R in action读书笔记(13)第十章 功效分析

10.3 绘制功效分析图形

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

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="correlationcoefficient (r)",

ylab="sample size(n)")

for(i in 1:np){

lines(r,samsize[,1],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="grey89")

title("sanple size ")

legend("topright",title="power",as.character(p),fill=colors)

R in action读书笔记(13)第十章 功效分析

10.4 其他软件包

asypow 通过渐进似然比方法计算功效

PwrGSD 组序列设计的功效分析

pamm 混合模型中随机效应的功效分析

powerSurvEpi 流行病研究的生存分析中功效和样本量的计算

powerpkg 患病同胞配对法和TDT(TransmissionDisequilibrium Test,传送不均衡检验)设

计的功效分析

powerGWASinteractionGWAS交互作用的功效计算

pedantics 一些有助于种群基因研究功效分析的函数

gap 一些病例队列研究设计中计算功效和样本量的函数

ssize.fdr 微阵列实验中样本量的计算

R in action读书笔记(13)第十章 功效分析