我更喜欢通过模拟进行超越基础的功率分析。使用预制包装，我永远不太确定正在做出哪些假设。

使用 R 模拟功率非常简单（而且价格合理）。

1. 决定你认为你的数据应该是什么样子以及你将如何分析它
2. 编写一个函数或一组表达式来模拟给定关系和样本大小的数据并进行分析（最好使用函数，因为您可以将样本大小和参数转换为参数，以便更容易尝试不同的值）。函数或代码应返回 p 值或其他测试统计信息。
3. 使用该replicate函数从上面运行代码多次（我通常从大约 100 次开始，以了解需要多长时间并获得正确的总体区域，然后将其增加到 1,000，有时为 10,000 或 100,000我将使用的最终值）。您拒绝原假设的次数比例就是功效。
4. 对另一组条件重做上述操作。

这是一个带有序数回归的简单示例：

library(rms)

tmpfun <- function(n, beta0, beta1, beta2) {
    x <- runif(n, 0, 10)
    eta1 <- beta0 + beta1*x
    eta2 <- eta1 + beta2
    p1 <- exp(eta1)/(1+exp(eta1))
    p2 <- exp(eta2)/(1+exp(eta2))
    tmp <- runif(n)
    y <- (tmp < p1) + (tmp < p2)
    fit <- lrm(y~x)
    fit$stats[5]
}

out <- replicate(1000, tmpfun(100, -1/2, 1/4, 1/4))
mean( out < 0.05 )

除了 Snow 的出色示例之外，我相信您还可以通过从具有您效果的现有数据集重新采样来进行功率模拟。不完全是一个引导程序，因为您不是用相同的n进行采样，而是相同的想法。

所以这里有一个例子：我做了一个小的自我实验，结果是一个正的点估计，但因为它很小，在序数逻辑回归中几乎没有统计学意义。有了这个点估计，我需要多大的n ？对于各种可能的n，我多次生成数据集并运行序数逻辑回归并看到p值有多小：

library(boot)
library(rms)
npt <- read.csv("http://www.gwern.net/docs/nootropics/2013-gwern-noopept.csv")
newNoopeptPower <- function(dt, indices) {
    d <- dt[sample(nrow(dt), n, replace=TRUE), ] # new dataset, possibly larger than the original
    lmodel <- lrm(MP ~ Noopept + Magtein, data = d)
    return(anova(lmodel)[7])
}
alpha <- 0.05
for (n in seq(from = 300, to = 600, by = 30)) {
   bs <- boot(data=npt, statistic=newNoopeptPower, R=10000, parallel="multicore", ncpus=4)
   print(c(n, sum(bs$t<=alpha)/length(bs$t)))
}

输出（对我来说）：

[1] 300.0000   0.1823
[1] 330.0000   0.1925
[1] 360.0000   0.2083
[1] 390.0000   0.2143
[1] 420.0000   0.2318
[1] 450.0000   0.2462
[1] 480.000   0.258
[1] 510.0000   0.2825
[1] 540.0000   0.2855
[1] 570.0000   0.3184
[1] 600.0000   0.3175

在这种情况下，在n = 600 时，功率为 32%。不是很鼓舞人心。

（如果我的模拟方法是错误的，请有人告诉我。我正在写几篇医学论文，讨论计划临床试验的功率模拟，但我完全不确定我的精确实施。）

我会在 Snow 的答案中添加另一件事（这适用于通过模拟进行的任何功率分析） - 请注意您是在寻找 1 尾测试还是 2 尾测试。像 G*Power 这样的流行程序默认为 1 尾测试，如果您试图查看您的模拟是否与它们匹配（当您学习如何执行此操作时总是一个好主意），您需要先检查一下。

为了让 Snow 运行 1-tailed 测试，我会在函数输入中添加一个名为“tail”的参数，并在函数本身中添加如下内容：

 #two-tail test
  if (tail==2) fit$stats[5]

  #one-tail test
  if (tail==1){
    if (fit$coefficients[5]>0) {
          fit$stats[5]/2
    } else 1

1 尾版本基本上检查系数是否为正，然后将 p 值减半。