#功效分析
#可以帮助在给定置信度的情况下，判断检测到给定效应值是所需的样本数
#反过来它也可以帮助你在给定置信度水平情况下，计算某样本量内能检测到给定效应值的概率
#研究过程的关注点
#样本大小，实验设计中每种条件/组中观测的数目
#显著性水平，由I型错误的概率来定义。也可以把它看做效应不发生的概率
#功效，多大的把握检测到它
####1-显著性水平(效应不发生的概率)=置信区间
#效应值，在备择或研究假设下效应的量，效应值的表达式依赖于假设检验中使用的统计方法
#研究目标是维持一个可接受的显著性水平，尽量使用较少的样本，然后最大化统计检验的功效
#最大化发现真实效应的几率，并最小化发现错误效应的几率，同时把研究成本控制在合理的范围内
#四个量紧密联系，给定其中任意三个量，并可以推算出第四个量
#结论：
#功效分析不仅可以帮助你判断在给定置信度和效应值的前提下所需的样本量，也能说明在给定样本量
#时检测到要求效应值的概率。对于限定误报效应显著性的可能性（I型错误）和正确检测真实效应（功效）的可能性的平衡
#典型的功效分析是一个交互性的过程，研究者会通过改变样本量，效应值，预期显著性水平和预期功效水平等参数
#来观测它们对于其他参数的影响
#功效分析的一个重要附加效应是引起方向性的转变，它鼓励不要仅仅关注于二值型（即效应值存在还是不存在）的假设检验
#而应该仔细思考效应值增加的意义，期刊编辑越来越多地要求作者在报告研究结果的时候即包含p值又包含效应值
#p值，针对零假设，显著性检验
#许多研究的特点是“零假设检验”（null hypothesis significance testing，NHST），
#这种情况下，研究人员主要都专注于将P值控制在低于统计显著性的假设值，这通常视为结果“真实性”的指标，
#也就是说研究发现不是侥幸得来
#效应量，能反馈误差有多大，相关性有多强
#例如，当我们说“吸烟和一个人的性生活满意度显著相关。（p < .001）”，
#但我们不知道它们的关联性有多强，如果我们发现相关系数是.1，那么实际上它们相关性是挺弱的；
#又如果相关系数是.7，那么实际上它们彼此间的关系是挺有趣的
#当同时提供效果量与置信区间（可能的效应值概率）时，
#相较于只有P值（不论P值大小）更能让我们更有效评估数据间的相关性
#效果量应该是怎么样呢？
#APA报告提供一些有用的建议：「如果度量单位于实际面来看是有意义的（比如一天抽几根烟），
#那么一般非标准单位（回归系数、平均数）会比标准单位（r或d）来得好，这也等同于说明效应值符合实践和理论。」
#效果量可以用原单位呈现，像是解决检验问题的平均数，这是最容易让人理解的方式。
#不过，像是Cohen's d这样的效应量标准（如units-free）单位有时候也是很有用的，
#这种效果量根据统计检验的不同而改变，在方差分析（ANOVA）中你可以用埃塔平方（部分埃塔平方也可）、ω平方或是F检验，
#您也可以利用以下几个连接找到一些有用的论文：
#范例一
install.packages("pwr")
library(pwr)
#pwr.t.test(n=,d=,sig.level=,power=,type=,alternative=)
#type指定检验类型，双样本t检验(two.sample),单样本t检验(one.sample),或者相依样本t检验paried,默认双样本检验
#如果两组中样本大小不同，使用pwr.t2n.test(n1=,n2=,d=,sig.level=,power=,alternative=)
pwr.t.test(d=0.8,sig.level=0.05,power=0.9,type="two.sample",
alternative="two.sided")
#n = 33.82555 样本大小
#d = 0.8 效应值，即标准化的均值之差，即统计量
#sig.level = 0.05 显著性水平，即alpha
#power = 0.9 功效水平
#统计检验是双侧检验two.sided，单侧检验less或greater,默认是双侧检验
#alternative = two.sided
#结果表明：每组中你需要34个受试者（双样本，所以总共68人），这样才能保证有90%的把握检测到
#0.8的效应值，并且最多5%的可能性会五宝差异存在
#范例二
pwr.t.test(n=20,d=0.5,sig.level=0.01,type="two.sample",
alternative="two.sided")
#n = 20 已知，总样本为40,40/2=20
#d = 0.5 已知，总体均值0.5个标准偏差的差异
#sig.level = 0.01 已知，误报差异的概率限制在1%以内
#power = 0.1439551 未知，待求
#alternative = two.sided 已知，双样本检验
#结论：在0.01的先验显著性水平下，每组20个受试者，已知因变量的标准差为1.25s(分母)，差异
#为0.5个标准偏差（结果），则差值为0.5*1.25=0.625(分子)
#有低于14%的可能性断言差值为0.625或者不显著
#换句话说，有86%（1-0.1439551=0.86）的可能性错过你要寻找的效应值
#方差分析
pwr.anova.test(k=,n=,f=,sig.level=,power=)
#可以对平衡单因素方差分析进行功效分析
#范例一
pwr.anova.test(k=5,f=0.25,sig.level=0.05,power=0.8)
#k = 5 分组的个数
#n = 39.1534 各组中的样本大小
#f = 0.25 F统计量
#sig.level = 0.05 显著性水平
#power = 0.8 功效值
#结论：用五个组做单因素方差分析，要达到0.8的功效，效应值为0.25，并选择0.05的显著性水平，
#计算各组需要的样本大小
#总样本大小为5*39=195
#相关性
#范例
pwr.r.test(r=0.25,sig.level=0.05,power=0.90,alternative="greater")
#n = 133.2803 观测数目
#r = 0.25 效应值(通过线性相关系数衡量)
#sig.level = 0.05 显著性水平
#power = 0.9 功效水平
#alternative = greater 双边检测two.sided,单边检测less or greater
#结论：需要134个受试者来评价抑郁与孤独的关系
#线性模型
#范例
pwr.f2.test(u=3,f2=0.0769,sig.level=0.05,power=0.90)
#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=193
#比例检验
#范例
#ES.h(p1,p2) 计算响应的效应大小 2*asin(sqrt(p1))-2*asin(sqrt(p2))
#p1 假设某一种更贵的新药如果能缓解65%使用者的症状，就会投放到市场中
#p2 假设某一种流行药物能缓解60%使用者的症状
#当各组中n不相同时
#pwr.2p2n.test(h=,n1=,n2=,sig.level=,power=)
pwr.2p.test(h=ES.h(0.65,0.6),sig.level=0.05,power=0.9,
alternative="greater")
#h = 0.1033347 效应值
#n = 1604.007 各组相同的样本量
#sig.level = 0.05 显著性水平，95%的把握不会误得结论
#power = 0.9 90%的把握得出新药更有效的结论
#alternative = greater two.sided双尾检验，单尾检验less 或greater
#结论：本研究需要1605个人试用新药，1605个人试用已有药物
#卡方检验
#范例
#人种和晋升之间的关系
#70% 白种人 10% 美国黑人 20% 西班牙裔人
#60% 白种人 30% 美国黑人 50% 西班牙裔人 容易晋升
#白种人 晋升0.7*0.6=0.42 不晋升0.7*（1-0.6）=0.28
#美国黑人 晋升0.1*0.3=0.03 不晋升0.1*（1-0.3）=0.07
#西班牙裔人 晋升0.2*0.5=0.1 不晋升0.2*（1-0.5）=0.1
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=0.1853,df=2,sig.level=0.05,power=0.9)
#w = 0.1853 效应值
#N = 368.5317 样本数
#df = 2 自由度，双因素列联表的自由度=(r-1)(c-1)
#sig.level = 0.05 显著性水平
#power = 0.9
#结论：结果表明，在既定的效应值，功效水平和显著性水平下，该研究需要369个受试者才能检验人种和晋升的关系
#在新情况中选择合适的效应值
#Cohen效应值基准
#社科类研究得出的一般性建议：
#统计方法 效应值测量 建议的效应值基准(小，中，大)
#方差分析 f值 0.10<0.25<0.40
#t检验 d值 0.20<0.50<0.80
#线性模型 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个受试者
#方差分析 f值 0.10<0.25<0.40
pwr.anova.test(k=5,n=25,sig.level=0.05,f=0.1) #反馈power = 0.1180955
pwr.anova.test(k=5,n=25,sig.level=0.05,f=0.25) #反馈power = 0.5738
pwr.anova.test(k=5,n=25,sig.level=0.05,f=0.40) #反馈power = 0.9569163
x <- pwr.anova.test(k=5,n=25,sig.level=0.05,f=0.40)
class(x) #返回"power.htest"
x$n #k,n,f,sig.leve,power 调用参数都可以访问
class(x$n) #返回 "numeric"
x$method #[1] "Balanced one-way analysis of variance power calculation"，抬头行信息
x$note #[1] "n is number in each group" 末尾行信息
#不小于当前小数的最小整数
ceiling(0.2) #返回1
#范例二
library(pwr)
es <- seq(0.1,0.5,0.01) #生成等差序列，从0.10开始，步长为0.01，到0.50结束
nes <- length(es) #序列长度为41
samsize <- NULL #生成空对象
for(i in 1:nes){
result <- pwr.anova.test(k=5,f=es[i],sig.level=0.05,power=0.9)
samsize[i] <- ceiling(result$n) #返回不比当前n小的最小整数
}#以效率值数列为依据，根据不同效率值，求不同样本N
plot(samsize,es,type="l",lwd=2,col="red",
ylab="Effect Size",
xlab="Sample Size (per cell)",
main="One Way ANOVA with Power=0.9 and Alpha=0.05") #基于样本N,效应值画散点图
#结论：各组样本量高于200个观测时，再增加样本已经效果不大了
#绘制功效分析图形
library(pwr)
r <- seq(0.1,0.5,0.01) #起始值0.1，结束值0.5，步长为0.01，生成相关系数向量
nr <- length(r) #统计相关系数向量长度
p <- seq(0.4,0.9,0.1)#其实质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.sided")
samsize[j,i] <- ceiling(result$n)
}
}
xrange <- range(r) #查看相关系数向量的值域，最小值0.1~最大值0.5，
yrange <- round(range(samsize)) #取结果矩阵的最小值~最大值，若存在小数则对它们取整
colors <- rainbow(length(p)) #基于效应值的不同取值生成彩虹颜色向量
plot(xrange,yrange,type="n",xlab="Correlation Coefficient (r)",
ylab="Sample Size (n)")#画空图，仅画出x轴，y轴
for(i in 1:np){
lines(r,samsize[,i],type="l",lwd=2,col=colors[i]
)
}
#基于所有效应值(效应值个数=结果矩阵列数)；
#取相关系数向量，结果集中每一列，画成一条直线；
#横向，基于纵轴刻度，画底格线
#seq(0,yrange[2],50) #起始值0，结束值yrang的第二值(最大值)，步长为50，生成序列
abline(h=seq(0,yrange[2],50),lty=2,col="grey89")
#纵向，基于横轴刻度，画底格线
#seq(xrange[1],xrange[2],0.02) #起始值xrang的第一值(最小值)，结束值xrang的第二值(最大值)，步长为0.02，生成序列
abline(v=seq(xrange[1],xrange[2],0.02),lty=2,col="gray89")
#加标题
title("Sample Size Estimation for Correlation Studies\n Sig=0.05 (Two-tailed)")
#加图例
legend("topright",title="Power",as.character(p),fill=colors)
扫一扫
494
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
