By

文章目录
  1. 功效分析
    1. 1、假设检验速览
    2. 2、用 pwr 包做功效分析
      1. 2.1、t检验
      2. 2.2、方差分析
      3. 2.3、相关性
      4. 2.4、线性模型
      5. 2.5、比例检验
      6. 2.6、卡方检验
    3. 3、绘制功效图形

功效分析

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

1、假设检验速览

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

功效分析就是给定三个值得情况下，求得另外一个值

2、用 pwr 包做功效分析

2.1、t检验

1
2
3
4
5
6
7
pwr.t.test(n=,d=,sig.level=,power=,alternative=)
#n: 样本大小
#d: 效应值,d=(u1-u2)/o, u1与u2为两组均值，o为标准差
#sig.level: 显著水平
#power: 功效水平
#type: 检验类型，双样本t检验（two.sample）、单样本t检验（one.sample）或相依样本t检验（paired）。默认为双样本t检验。
#alternative: 指统计检验是双侧检验（two.sided）还是单侧检验（less或greater）。默认为双侧检验
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
> ##t检验
> library(pwr)
> 
> #实例1、手机与驾驶反应时间实验，根据过去经验知道反应时间的标准差为1.25s,当前认为相差1s表示巨大差异，此时可
> #设定d=1/1.25=0.8或更大，如果差异存在，你希望90%把握检测到，由于随机变异性存在，希望95%把握不会误报差异显著
> pwr.t.test(d=0.8,sig.level = 0.05,power = 0.9,type = "two.sample",alternative = "two.side")
     Two-sample t 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
> 
> #结果显示需要每组34个受试者，一共2组68人
> 
> #实例2、如果在比较实例1的2种情况时候，想检测到总体均值0.5个标准差的差异，将误报的几率限制在%1，只有40个受试者，
> #那么能检测到这种差异的概率是多大？
> pwr.t.test(n=20,d=0.5,sig.level = 0.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
> 
> #结果表明只有14%的把握能检测到这种差异，即差值为0.5*1.25=0.625s,86%的可能性会错过，此时考虑是否应该进行该测试
> 
> #如果两组的样本大小不同，使用：
> pwr.t2n.test(n1=20,n2=30,d=0.5,sig.level = 0.01,alternative = "two.side")
     t test power calculation 
             n1 = 20
             n2 = 30
              d = 0.5
      sig.level = 0.01
          power = 0.183246
    alternative = two.sided

2.2、方差分析

效应值f计算公式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
> ##方差分析
> #pwr.anova.test(k=组数,n=各组样本大小,f=效应值,sig.level=显著水平,power=功效)
> pwr.anova.test(k=5,f=0.25,sig.level=0.05,power=0.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个组的均值

2.3、相关性

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> ##相关性
> #pwr.r.test(n=观测数,r=效应值（线性相关系数）,sig.level = 显著水平,power = 功效,alternative =c("two.sided", "less","greater"))
> #抑郁与孤独的关系，设定相关系数为0.25,大于0.25表示有关系拒绝零假设，显著水平为0.05,90%的信心拒绝零假设，
> pwr.r.test(r=0.25,sig.level = 0.05,power = 0.9,alternative ="greater")
     approximate correlation power calculation (arctangh transformation) 
              n = 133.2803
              r = 0.25
      sig.level = 0.05
          power = 0.9
    alternative = greater
> #需要134个受试者

2.4、线性模型

效应值f2计算公式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
> ##线性模型
> #pwr.f2.test(u=分子自由度,v=分母自由度,f2=效应值,sig.level = ,power = )
> #当要评价一组预测变量对结果的影响程度时，适宜用第一个公式来计算f2；当要评价一组预
> #测变量对结果的影响超过第二组变量（协变量）多少时，适宜用第二个公式。
> #现假设你想研究老板的领导风格对员工满意度的影响，是否超过薪水和工作小费对员工满意度的影响。
> #领导风格可用四个变量来评估，薪水和小费与三个变量有关。过去的经验表明，薪水
> #和小费能够解释约30%的员工满意度的方差。而从现实出发，领导风格至少能解释35%的方差。
> #假定显著性水平为0.05，那么在90%的置信度情况下，你需要多少受试者才能得到这样的方差贡
> #献率呢？
> 
> #此时u=3(总预测变量-集合B预测变量)，f2=(0.35-0.3)/(1-0.35)=0.0769
> pwr.f2.test(u=3,f2=0.0769,sig.level =0.05 ,power = 0.9)
     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

2.5、比例检验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
> ##比例检验
> #pwr.2p.test(h=效应值,n=各组相同的样本量,sig.level = ,power = )
> #当n不同时候，使用
> #pwr.2p.test(h=效应值,n1=,n2=,sig.level = ,power = )
> #h使用ES.h(p1,p2)计算
> #假定你对某流行药物能缓解60%使用者的症状感到怀疑。而一种更贵的新药如果能缓解65%
> #使用者的症状，就会被投放到市场中。此时，在研究中你需要多少受试者才能够检测到两种药物
> #存在这一特定的差异？
> 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
> 
> #每组1605个人

2.6、卡方检验

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
##卡方检验
#常常用来评价两个类别型变量的关系。典型的零假设是变量之间独立，备择假设是不独立。 
pwr.chisq.test(w=效应值,N=总样本大小,df=自由度,sig.level = ,power = )
#ES.w2(P)可以用来计算双因素列联表备择假设的效应值，P为一个假设的双因素概率表
#假设你想研究人种与工作晋升的关系。你预期样本中70%是白种人， 10%是美国黑人， 20%西班牙裔人。而且，你认为相比30%的美国黑人和50%的西班牙裔人， 60%的
#白种人更容易晋升。晋升概率如表所示。
#|人种|晋升比例|未晋升|
#|白  |0.42    |  0.28|
#|黑  |0.03    |  0.07|
#|西裔|0.10    |  0.10|
prob<-matrix(c(0.42,0.28,0.03,0.07,0.10,0.10),byrow = TRUE,nrow = 3)
ES.w2(prob)
pwr.chisq.test(w=ES.w2(prob),df=2,sig.level = 0.05,power = 0.9)
#需要369个测试者

3、绘制功效图形

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
##绘制功效分析图形
library(pwr)
#生成一系列相关系数和功效值
r<-seq(0.1,0.5,0.01)
nr<-length(r)
p<-seq(0.4,0.9,0.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 = 0.05,power=p[i],alternative = "two.side")
    samsize[j,i]<-ceiling(result$n)
  }
}
#返回最小值与最大值
xrange<-range(r)
#返回最小值与最大值，向上取整
yrange<-round(range(samsize))
#产生一组彩虹色
colors<-rainbow(length(p))
#绘制坐标系
plot(xrange,yrange,type="n",xlab="相关系数",ylab="样本大小")
#每一种功效绘制曲线
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],0.02),lty=2,col="grey89")
title("相关性分析的样本大小估计，显著水平0.05，双尾")
legend("topright",title="功效",as.character(p),fill = colors)

文章目录
  1. 功效分析
    1. 1、假设检验速览
    2. 2、用 pwr 包做功效分析
      1. 2.1、t检验
      2. 2.2、方差分析
      3. 2.3、相关性
      4. 2.4、线性模型
      5. 2.5、比例检验
      6. 2.6、卡方检验
    3. 3、绘制功效图形