R语言实战笔记--第十章 功效分析&样本量

R语言实战笔记–第十章 功效分析&样本量

标签（空格分隔）： R语言 功效分析 样本量

功效及两类错误

　　功效及样本量在概率论与数理统计的假设检验部分里面有说明，可以返回查看原理，这里以R语言实战的描述简单复习一下。

　　以书中例子说明，两组玩手机开车试验，它们的统计量应该是双总体，标准差未知的配对t检验，所以，若它们没有什么区别，则零假设为它们的反应时间之差为0，所以实际统计量应该是(x¯1−x¯2)/（Sw1n1+1n2−−−−−−−√），其中Sw=(n1−1)S21+(n2−1)S22n1+n2−2，公式需求方差齐性。

　　而书中简化的公式它需求样本量相等，在此条件下，才服从2n−2的t分布，否则服从n1+n2−2的t分布。

　　回归功效分析，首先简单了解一下两类错误，用一个简单例子说明：

　　背景or需求：某个盒子里面装着一个人，他/她用变声器跟我们说话，我们想确认“它”是男的还是女的。

　　假设：我们假设他是男的，那么，备择假设自然是它不是男的。

　　方法：检验方法是我们要用什么检验方法，统计量是什么，这就涉及了检验方法问题，假如我们问“你早餐吃什么？”那么很明显，方法错误，因为这个问题根本不可能有效检验性别，假如我们问“你是站着撒尿的吗？”这个问题很明显，它具有极大的概率（置信水平，比如这个问题可以检验99%的男人是站着撒尿，对应的显著水平是1%，即只有1%的男性不是站着撒尿）可以检验出“它”到底是不是男的。

　　检验结果：如果他回答“是”，那么我们有极大的概率判断，他是男的，如果回答“否”，那么我们判断他不是男的。

　　两类错误：先说正确的情况，他是男的，他回答了“是”，我们判断“他是男的”与她是女的，她回答了“否”，我们判断“他不是男的”，这两个情况，都是我们判断正确了。说说第一类错误，他是男的，但是，他回答了“否”，我们由此判断他不是男的，此时，为第一类错误，即相对我们的原假设来说，注意，是相对原假设来说，他是男的，但他不站着撒尿，由于他的奇葩行为（即原假设中的异常值）导致我们的检验结果出错，这是第一类错误，这类错误的概率由我们的显著水平给出。第二类错误，就是相对备择假设来说，注意，是相对备择假设来说，她是女的，但她站着撒尿，让我们错误的判断“她是男的”，这类错误是由于备择假设的异常值导致我们检验出错，一般来说，这个概率我们是不知道的，因为我们在设定检验问题的时候，通常是可以通过显著水平去确定这个检验问题能保证多少概率能判断出对还是不不对，但这个检验问题没办法确定不是男人而站着撒尿的概率。

　　功效的定义为（1-第二类错误的概率），而其它在个量分别为样本大小、显著水平以及效应值（为总体的标准化均值差，它依赖于历史经验或估计经验，通常，在检验两者是否相等的时候，效应值为0）

　　

功效分析

　　书中给出的是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不相等的两样本）

t检验

　　书中的描述太过专业，我试着转化一下更好理解。首先看函数，与书中一样：

pwr.t.test(n=,d=,sig.level=,power=,type=,alternative=)

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

n(n1,n2)为样本大小，d为效应值（标准化的均值差，d=μ1−μ2σ），sig.level为显著水平，power为功效，type为检验类型（单样本one.sample，双样本two.sample，配对样本paired，默认为双样本），alternative为检验方式（单侧less或greater，双侧two.sided，默认为双侧）

　　四个变量任意输入三个得到第四个，检验的类型和方式可以使用默认值也可以输入，注意区分样本数相同及样本数不同时使用的函数是不一样的。然后再解释一下书中例子：开车玩手机与开车不玩手机，在偶到突发情况的反应速度差异。若我们需要得到样本大小，则需要设定其它三个参数（即效应值、显著水平以及功效值），下面分别介绍这些值的概念以及怎么得到：

　　样本大小n：本次试验所进行的样本大小，通常大小是由其它三个参数确定，但有时候，我们只能够获得这么大的样本，需要确认功效情况也比较常见。

　　效应值d：由历史水平给出标准差（原假设的标准差，在例子里，是不玩手机时的标准差），假如我们只是检验它们是否不一样，那么功效就是0，假如我们不单单检测他们是否不一样，还规定相差1s就是巨大的差异，且我们要检测出这种巨大的差异，那么功效就是1/1.25=0.8，通俗的来说，效应就是一个这样的数：样本与预期的均值差的标准化参数。

　　显著水平α：一般为0.05，就是确定检验统计量的误报率，换个说法，就是有（1-显著水平）%的概率不会误报，以上面说明两类错误的例子来说，就是你相信问题“站着撒尿”能判断出他是男性的把握是多少，比如说99%，就是你相信99%的男性都是站着撒尿的，比如60%，就是你相信60%的男性是站着撒尿的，剩下的40%你不相信，从实际上看，显然应该选择前者（99%）。

　　功效power：与显著水平几乎对立，也继续上面的例子，99%的男性站着撒尿，那么如果她是女性，那么我们有多大把握拒绝原假设（即女性回答不是站着撒尿的概率）呢，这个就是功效值，也就是1-II型错误的概率。

　　用一句语总结这四个参数的关系就是：对双样本检验，以（得到）两个样本均为n大小（或分别是n1，n2大小）的样本，以（得到）α的误报率，以（得到）d的标准化差异（可接受差异），假如出现异常，我 也能以（得到）power的概率去拒绝异常。可以任意把其中一个换成得到。

　　

方差功效分析

　　此处只能对平衡单因素方差分析进行，函数如下：

pwr.anova.test(k=,n=,f=,sig.level=,power=)

其中k是组数，n是样本大小，f是效应值，其它同上。对单因素方差分析，效应f=∑ki=1pi(μi−μ)2σ2−−−−−−−−−−√，pi=niN，ni组i的观测数目，N总观测数目，μi组i均值，μ总体均值，σ2组内误差方差

相关功效分析

　　函数如下：

pwr.r.test(n=,r=,sig.level=,power=,alternative=)r为相关系数，其它同上

线性模型功效分析

　　对线性模型，如多元回归分析的功效分析，函数如下：

pwr.f2.test(u=,v=,f2=,sig.level=,power=)

其中u和v分别是分子自由度和分母自由度，分子自由度就是预测变量的个数（含常量），也是回归平方和的自由度（p），在计算两组变量的时候，它还需要送去第二组变量（待估组）的变量数; 分母自由度就是残差平方和的自由度(n-p-1)；f2是效应值它有两个公式可选，当需求是评价一组变量对结果的影响程度时，使用f2=R21−R2,R2是多重相关性的总体平方值（不知道是什么。。。初步理解就是线性拟合后的R平方，不知道对不对）；当需求是评价一组变量(A)超过另一组变量(B)多少时，适用公式f2=R2AB−R2A1−R2AB，R2A为集合A中变量对总体方差的解释率，R2AB为集合A与B共同的解释率。

　　书中的翻译有坑，例子没问题，翻译过来就出问题了，原翻译是，领导风格至少有5%以上的影响，所以R2AB就是薪水+小费（为A集合），领导风格为B集中的和即35%，而书中翻译是直接翻译是领导风格至少解释35%的方差，显然，如果按翻译的说法，R2AB就应该为65%才对！

比例检验功效分析

　　先看函数：

pwr.2p.test(h=,n=,sig.level=,power=，alternative=),效应值h=2arcsinp1−−√−2arcsinp2−−√，可以使用ES.h(p1,p2)计算得到

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

卡方检验功效分析

　　函数如下：

pwr.chisq.test(w=,n=,df=.sig.level=,power=),效应值w=∑mi−1(p0i−p1i)2p0i−−−−−−−−−−−−√，其中p0为第i单元格在零假设时的概率p1为备择假设时的概率。可以由函数ES.w2(p)得到，其中p=matrix(c(p1,p2),byrow=T,nrow=r)，自由度df=(r-1)(c-1)

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

　　在功效分析中，预期的效应值是最难决定的参数，它通常需要对主题有一定的了解并有相应的测量经验，例如使用过去的的研究数据得到的参考值。但当进行全新的研究时，并无任何参考值时，可使用一些建议的基准值进行研究，如下表：

统计方法	效应值测量	小效应值基准	中效应值基准	大效应值基准
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

　　注意的是，本表中的建议效应值是在社科类研究给出的一些基准，并不适用所有研究领域！

　　另外，可以由绘图帮助确认效应值，代码如下：

es <- seq(0.1, 0.5, 0.01)
nes <- length(es)
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)
}
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")

　　从图中可以看出，当n>200时，再增加样本，效应值的变体不多了，意味着，200个样本，几乎确定可以检测到最小的效应值了（即两样本的最小差异，但实际上需不需要到这么小，或者需要更小，就需要结合实际云考虑了，但在这张图的背景下，只能达到这么多，如果不需要那么小的效应值，那么可以减少一些样本量）。

　　

绘制功效分析图形

　　其实绘制功效分析图形的原理就是列举法，把每个点按一定的等差列出来，把它们绘制成一条条的直线，即可以得到各种情况下的图形。以下是书中代码：

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.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], 0.02), lty = 2, col = "grey89")
title("Sample Size Estimation for Correlation Studies\nSig=0.05 (Two-tailed)")
legend("topright", title = "Power", as.character(p), fill = colors)

　　

图形如下，显示的是在各功效值下，不同样本量与相关系数的关系：