R语言中的功效分析

一、假设检验速览

在统计假设检验中，首先要对总体分布参数设定一个假设（零假设H0），然后从总体分布中抽样，通过样本计算所得的统计量来对总体参数进行推断。假定零假设为真，如果计算获得观测样本的统计量的概率非常小，便可以拒绝原假设，接受它的对立面（称作备择假设或者研究假设H1）。

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

二、用pwr 包做功效分析

对于每个函数，用户可以设定四个量（样本大小、显著性水平、功效和效应值）中的三个量，第四个量将由软件计算出来。

1、t检验

如果你根据过去的经验知道反应时间有1.25 s的标准偏差，并认定反应时间1 s的差值是巨大的差异，那么在这个研究中，可设定要检测的效应值为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.sided")

 Two-sample t test power calculation 
              n = 33.82554
              d = 0.8
      sig.level = 0.05
          power = 0.9
    alternative = two.sided
n=34*2=68

每组中你需要34个受试者（总共68人），这样才能保证有90%的把握检测到0.8的
效应值，并且最多5%的可能性会误报差异存在。

如果二个样本不同怎么办？

2、方差分析

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
总样本大小为5 × 39，即195。注意，本例中需要估计在同方差时五个组的均值。

3、相关性

pwr.r.test(r = 0.25, sig.level = 0.05, power = 0.9, 
    alternative = "greater")

        n = 133.2803
              r = 0.25
      sig.level = 0.05
          power = 0.9
    alternative = greater

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

4、线性模型

pwr.f2.test(u = 3, f2 = 0.0769, sig.level = 0.05, 
    power = 0.9)

5、卡方检验

prob <- matrix(c(0.42, 0.28, 0.03, 0.07, 0.1, 0.1), 
    byrow = TRUE, nrow = 3)
ES.w2(prob)  计算效应值
pwr.chisq.test(w = 0.1853, df = 3, sig.level = 0.05, 
    power = 0.9)

Chi squared power calculation 


              w = 0.1853
              N = 412.7286
             df = 3
      sig.level = 0.05
          power = 0.9

效应值是最难计算的

查看样本与功效之间的关系

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")

实验设计中，这样的图形有助于估计不同条件时的影响值。例如，从图形可以看到各组样本量高于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)