机器算法验证 - 计算 R 中 Fisher 精确检验的功效 - 吾爱随笔录

计算 R 中 Fisher 精确检验的功效

机器算法验证 r 假设检验统计能力列联表渔民精确测试

2022-03-23 16:26:34

假设我有以下 3 x 3 列联表：

      |    T1   |    T2   |    T3   |
------+---------+---------+---------+---
  G1  |   18    |   15    |   65    | 
------+---------+---------+---------+---
  G2  |   20    |   10    |   70    |
------+---------+---------+---------+---
  G3  |   15    |   55    |   30    |

然后我可以在 R 中运行 Fisher 精确检验（使用 Monte Carlo 模拟选项），如下所示：

table = matrix(c(18,20,15,15,10,55,65,70,30), 3, 3)
fisher.test(table, simulate.p.value=TRUE)

，产生以下结果（p 值）：

Fisher's Exact Test for Count Data with simulated p-value (based on
    2000 replicates)

data:  table
p-value = 0.0004998
alternative hypothesis: two.sided

在我的假设检验中，我还将有兴趣观察 II 型错误的概率（假阴性率，beta）以及当原假设为假时正确拒绝原假设 (H0) 的概率（统计功效，1 - 测试版）。有人知道 R 和 Fisher 对 mxn（m > 2 和 n > 2）列联表的精确检验是否可行？你可能知道我可以用其他一些工具来计算这个吗？谢谢！

PS：请注意，上面的列联表只是一个示例，实际获得的 p 值可能 > alpha（在我的情况下为 0.05），这是更有趣的情况。

1个回答

您在这里要求的是事后功率分析。（更具体地说，“正确拒绝原假设的概率”是幂，1 次幂是 beta，“II 型错误的概率”。你要两者，但我们只需要一个知道另一个.) 我们将您现有的数据集作为真实数据生成过程的替代假设/模型。我不知道有一个专门的、预先存在的函数（例如，在pwr包中）来执行此操作，但是，是的，这可以在R. 你只需要模拟它。有关（相当多）有关功率分析的更多信息，并在中对其进行模拟R，您应该在此处阅读我的答案：逻辑回归功率分析的模拟 - 设计的实验. 在这种情况下，我将给出一个快速的、经过调整的版本来处理 Fisher 的精确检验。（我通常编写代码尽可能接近伪代码，以便更广泛地理解它，但因为这可能需要很长时间才能运行，所以我尝试尽可能多地移出for循环，并使用一些R的独特能力。）


table = matrix(c(18,20,15,15,10,55,65,70,30), 3, 3)
table
#      [,1] [,2] [,3]
# [1,]   18   15   65
# [2,]   20   10   70
# [3,]   15   55   30
N = sum(table)  # this is the total number of observations
N
# [1] 298
probs = prop.table(table)       
      # these are the probabilities of an observation
probs                           #  being in any given cell
#            [,1]       [,2]      [,3]
# [1,] 0.06040268 0.05033557 0.2181208
# [2,] 0.06711409 0.03355705 0.2348993
# [3,] 0.05033557 0.18456376 0.1006711
probs.v = as.vector(probs)      
      # notice that the probabilities read column-wise
probs.v
# [1] 0.06040268 0.06711409 0.05033557 0.05033557 0.03355705 
    0.18456376 0.21812081
# [8] 0.23489933 0.10067114
cuts = c(0, cumsum(probs.v))    
        # notice that I add a 0 on the front
cuts
# [1] 0.00000000 0.06040268 0.12751678 0.17785235 0.22818792 
      0.26174497
# [7] 0.44630872 0.66442953 0.89932886 1.00000000

set.seed(4941)      # this makes it exactly reproducible
B      = 10000      # number of iterations in simulation
vals   = runif(N*B) # generate random values / probabilities
cats   = cut(vals, breaks=cuts, labels=c("11", "21", "31", 
             "12", "22", "32", "13", "23", "33"))
cats   = matrix(cats, nrow=N, ncol=B, byrow=F)
counts = apply(cats, 2, function(x){ as.vector(table(x)) })

rm(table, N, vals, probs, probs.v, cuts, cats) 
p.vals = vector(length=B) # this will store the outputs
ptm = proc.time()         # this lets me time the simulation
for(i in 1:B){
  mat       = matrix(counts[,i], nrow=3, ncol=3, byrow=T)
  p.vals[i] = fisher.test(mat, simulate.p.value=T)$p.value
}
proc.time() - ptm               # not too bad, really
#  user  system elapsed 
# 28.66    0.32   29.08 
#
mean(p.vals>=.05)               
     # the estimated probability of type II errors is 0
# [1] 0
c(0, 3/B)                       
     # using the rule of 3 to estimate the 95% CI
# [1] 0e+00 3e-04

考虑到您的数据与 Fisher 精确检验中的原假设的偏离程度以及您拥有的数据量，此模拟在 10,000 次迭代中不会出现单一的 II 类错误。因为每次迭代都可以理解为从具有概率的二项分布中抽取 $p$ （我们将其估计为观察到的 II 型错误的比例），该模拟实际上是具有一些随机变异性的估计。我们可以形成一个 95% 的置信区间来界定 II 类错误的真实概率。为了绕过我们实际上没有发现任何类型 II 错误的事实，我们将使用规则 3 ( $3/N$ ) 来估计 CI 的上限。因此，真正的 II 型错误率的 95% CI 是 $[0,\ 0.0003]$ .

另一方面，@rvl 在评论中指出“[p] ost hoc 权力是一种愚蠢的做法”。这在很大程度上是正确的。我看到人们提出这样的论点，实际上，“我的结果并不显着，但我没有任何权力，所以没有理由相信我的理论是错误的”，这在任何层面上都是相当奇怪的. 另一方面，由于您的结果很重要，因此也不清楚知道您的研究的事后能力有什么不同。我发现理解事后权力可以在教学上帮助人们开始理解这个话题。我们也可以以此为起点进行先验功效分析，以规划未来的研究。

其它你可能感兴趣的问题

上一篇如何测量不使用 drop out 的多层感知器神经网络中发生的协同适应？下一篇将中位数扩展到大数据整数分布？