预赛：

• 正如G*Power 手册中所讨论的，有几种不同类型的功率分析，具体取决于您要解决的问题。（也就是说，、效应大小、和幂是相互关联的；指定其中任意三个可以让您求解第四个。） $N$$ES$$\alpha$

  • 在您的描述中，您想知道适当的来捕获您使用和 power = 80% 指定的响应率。这是先验的力量。 $N$$\alpha=.05$
  • 我们可以从事后功率（确定功率给定、响应率和 alpha）开始，因为这在概念上更简单，然后向上移动$N$

• 除了@GregSnow 的出色帖子之外，还可以在此处找到另一个非常棒的基于仿真的 CV 功率分析指南：计算统计功率。总结一下基本思路：

  1. 找出您希望能够检测到的效果
  2. 从那个可能世界生成 N 个数据
  3. 运行您打算对这些虚假数据进行的分析
  4. 根据您选择的 alpha 存储结果是否“显着”
  5. 重复多次（处（事后）功率的估计$B$$N$
  6. 要确定先验功率，请搜索可能的以找到产生所需功率的值 $N$

• 您是否会在特定迭代中发现显着性可以理解为概率（其中是幂）的伯努利试验的结果。次迭代中找到的比例使我们能够逼近真实的。为了获得更好的近似值，我们可以增加，尽管这也会使模拟花费更长的时间。 $p$$p$$B$$p$$B$

• 在 R 中，生成具有给定“成功”概率的二进制数据的主要方法是?rbinom

  • 例如，要以概率 p 获得 10 次伯努利试验的成功次数，代码为rbinom(n=10, size=1, prob=p),（您可能希望将结果分配给一个变量进行存储）
  • 您还可以使用?runif不太优雅地生成此类数据，例如，ifelse(runif(1)<=p, 1, 0)
  • 如果您认为结果是由潜在高斯变量介导的，您可以使用?rnorm将潜在变量生成为协变量的函数，然后将它们转换为概率pnorm()并在rbinom()代码中使用这些概率。

• 您声明您将“包括多项式项 Var1*Var1) 以解释任何曲率”。这里有一个混乱；多项式项可以帮助我们解释曲率，但这是一个交互项——它不会以这种方式帮助我们。尽管如此，您的回复率要求我们在模型中同时包含平方项和交互项。具体来说，您的模型将需要包括：、和，超出基本条款。 $var1^2$$var1*var2$$var1^2*var2$

• 虽然是在不同问题的背景下写的，但我在这里的回答：logit 和 probit 模型之间的差异有很多关于这些类型模型的基本信息。

• 正如当有多个假设时存在不同类型的 I 类错误率（例如，每个对比错误率、家庭错误率和每个家庭错误率），是否存在不同类型的幂*（例如，对于单个预先指定的效果，适用于任何效果，&适用于所有效果）。您还可以寻求检测特定效果组合的能力，或者寻求同时测试整个模型的能力。根据您对 SAS 代码的描述，我猜测它正在寻找后者。但是，根据您对情况的描述，我假设您希望至少检测交互效果。

  • *参考：Maxwell, SE (2004)。心理学研究中动力不足的研究持续存在：原因、后果和补救措施。 心理学方法, 9 , 2 , pp. 147-163。
  • 您的效果非常小（不要与低响应率混淆），因此我们会发现很难达到良好的效果。
  • 请注意，尽管这些听起来都非常相似，但它们非常不一样（例如，很可能得到一个没有显着影响的显着模型——在这里讨论：回归如何显着但所有预测变量都不是显着？，或显着影响，但模型不显着 - 在这里讨论：线性回归中系数的显着性：显着 t 检验与非显着 F 统计量），这将在下面说明。

• 要以不同的方式思考与权力相关的问题，请参阅我的回答：如何在证明样本量合理的情况下报告估计相关性的一般精度。

R中逻辑回归的简单事后功效：

假设您假设的回复率代表了世界上的真实情况，并且您已经发出了 10,000 封信。检测这些影响的能力是什么？（请注意，我以编写“可笑的低效”代码而闻名，以下内容旨在易于理解，而不是针对效率进行优化；实际上，它很慢。）

set.seed(1)

repetitions = 1000
N = 10000
n = N/8
var1  = c(   .03,    .03,    .03,    .03,    .06,    .06,    .09,   .09)
var2  = c(     0,      0,      0,      1,      0,      1,      0,     1)
rates = c(0.0025, 0.0025, 0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002)

var1    = rep(var1, times=n)
var2    = rep(var2, times=n)
var12   = var1**2
var1x2  = var1 *var2
var12x2 = var12*var2

significant = matrix(nrow=repetitions, ncol=7)

startT = proc.time()[3]
for(i in 1:repetitions){
  responses          = rbinom(n=N, size=1, prob=rates)
  model              = glm(responses~var1+var2+var12+var1x2+var12x2, 
                           family=binomial(link="logit"))
  significant[i,1:5] = (summary(model)$coefficients[2:6,4]<.05)
  significant[i,6]   = sum(significant[i,1:5])
  modelDev           = model$null.deviance-model$deviance
  significant[i,7]   = (1-pchisq(modelDev, 5))<.05
}
endT = proc.time()[3]
endT-startT

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.042
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.017
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.035
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.019
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.022
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 0.168
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.001
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.065
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.017

所以我们看到 10,000 个字母并没有真正达到 80% 的功率（任何种类的）来检测这些响应率。（我不太确定 SAS 代码在做什么，无法解释这些方法之间的明显差异，但这段代码在概念上很简单——如果速度很慢的话——我花了一些时间检查它，我认为这些结果是合理的。）

逻辑回归的基于模拟的先验功效：

从这里开始，我们的想法是简单地搜索可能的，直到我们找到一个产生您感兴趣的功率类型所需水平的值。您可以编写任何搜索策略来处理它都可以（在理论）。考虑到捕捉如此小的效果所需的，值得考虑如何更有效地做到这一点。我的典型方法只是蛮力，即评估我可能合理考虑（但是请注意，我通常只考虑一个很小的范围，而且我通常使用非常小 ——至少与此相比。） $N$$N$$N$$N$

相反，我在这里的策略是将可能的括起来，以了解权力的范围。因此，我选择了一个并重新运行代码（启动相同的种子，nb 这需要一个半小时才能运行）。结果如下： $N$$N$

sum(significant[,1])/repetitions      # pre-specified effect power for var1
[1] 0.115
sum(significant[,2])/repetitions      # pre-specified effect power for var2
[1] 0.091
sum(significant[,3])/repetitions      # pre-specified effect power for var12
[1] 0.059
sum(significant[,4])/repetitions      # pre-specified effect power for var1X2
[1] 0.606
sum(significant[,5])/repetitions      # pre-specified effect power for var12X2
[1] 0.913
sum(significant[,7])/repetitions      # power for likelihood ratio test of model
[1] 1
sum(significant[,6]==5)/repetitions   # all effects power
[1] 0.005
sum(significant[,6]>0)/repetitions    # any effect power
[1] 0.96
sum(significant[,4]&significant[,5])/repetitions   # power for interaction terms
[1] 0.606

从这里我们可以看出，你的影响的大小变化很大，因此你检测它们的能力也不同。例如，的影响特别难以检测，即使有 50 万个字母，也只有 6% 的时间显着。另一方面，作为一个整体的模型总是明显优于空模型。其他可能性介于两者之间。尽管每次迭代都会丢弃大部分“数据”，但仍有可能进行大量探索。例如，我们可以使用矩阵来评估不同变量显着概率之间的相关性。 $var1^2$significant

最后我应该指出，由于您的情况复杂且很大，这并不像我在最初的评论中所怀疑/声称的那么简单。但是，您当然可以从我在这里提出的内容中了解总体上如何做到这一点，以及功率分析中涉及的问题。HTH。 $N$

@Gung 的回答非常有助于理解。这是我将使用的方法：

mydat <- data.frame( v1 = rep( c(3,6,9), each=2 ),
    v2 = rep( 0:1, 3 ), 
    resp=c(0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002) )

fit0 <- glm( resp ~ poly(v1, 2, raw=TRUE)*v2, data=mydat,
    weight=rep(100000,6), family=binomial)
b0 <- coef(fit0)


simfunc <- function( beta=b0, n=10000 ) {
    w <- sample(1:6, n, replace=TRUE, prob=c(3, rep(1,5)))
    mydat2 <- mydat[w, 1:2]
    eta <- with(mydat2,  cbind( 1, v1, 
                v1^2, v2,
                v1*v2,
                v1^2*v2 ) %*% beta )
    p <- exp(eta)/(1+exp(eta))
    mydat2$resp <- rbinom(n, 1, p)

    fit1 <- glm( resp ~ poly(v1, 2)*v2, data=mydat2,
        family=binomial)
    fit2 <- update(fit1, .~ poly(v1,2) )
    anova(fit1,fit2, test='Chisq')[2,5]
}

out <- replicate(100, simfunc(b0, 10000))
mean( out <= 0.05 )
hist(out)
abline(v=0.05, col='lightgrey')

这个功能测试v2的整体效果，可以换模型看其他类型的测试。我喜欢把它写成一个函数，这样当我想测试不同的东西时，我可以改变函数参数。您还可以添加进度条，或使用并行包来加快速度。

在这里我只做了 100 次重复，我通常从那个水平开始找到近似的样本量，然后当我在正确的球场上时增加迭代（当你有 20% 的功率时，不需要浪费时间在 10,000 次迭代上）。