寻求计算单样本比例假设检验功效的R内置函数

确实，R内置的power.prop.test()只支持双样本比例检验，不过我们可以自己动手实现单样本的计算，或者用模拟、扩展包的方式来解决你的硬币抛投问题。

先明确你的问题参数：

• 原假设 ( H_0: p = 0.5 )（硬币无偏）
• 备择假设 ( H_1: p = 0.49 ) 或 ( p = 0.51 )（硬币有偏）
• 显著性水平 ( \alpha = 0.05 )（默认）
• 目标功效（比如80%或90%）：即当硬币确实有偏时，我们能正确拒绝原假设的概率

下面提供几种可行的方法：

方法1：基于正态近似手动实现函数（纯内置函数，无需扩展包） #

当样本量较大时，二项分布可以用正态分布近似，我们可以推导样本量和功效的计算公式，编写一个自定义函数：

# 单样本比例检验的功效/样本量计算函数（正态近似）
power.prop.one.sample <- function(p0, p1, alpha = 0.05, power = 0.8, alternative = "two.sided") {
  # 处理检验方向，计算临界z值
  if (alternative == "two.sided") {
    z_alpha <- qnorm(1 - alpha/2)
  } else if (alternative == "greater") {
    z_alpha <- qnorm(1 - alpha)
  } else if (alternative == "less") {
    z_alpha <- qnorm(alpha)
  } else {
    stop("alternative must be 'two.sided', 'greater', or 'less'")
  }
  
  # 对应功效的z值
  z_beta <- qnorm(power)
  
  # 计算所需样本量（正态近似公式）
  numerator <- (z_alpha * sqrt(p0*(1-p0)) + z_beta * sqrt(p1*(1-p1)))^2
  denominator <- (p1 - p0)^2
  n <- ceiling(numerator / denominator)
  
  # 计算实际达到的功效（验证）
  se0 <- sqrt(p0*(1-p0)/n)
  se1 <- sqrt(p1*(1-p1)/n)
  
  if (alternative == "two.sided") {
    # 双侧检验的临界值
    crit_low <- p0 - z_alpha * se0
    crit_high <- p0 + z_alpha * se0
    # 备择假设下，样本比例落在拒绝域的概率
    actual_power <- pnorm(crit_low, mean = p1, sd = se1) + (1 - pnorm(crit_high, mean = p1, sd = se1))
  } else if (alternative == "greater") {
    crit <- p0 + z_alpha * se0
    actual_power <- 1 - pnorm(crit, mean = p1, sd = se1)
  } else {
    crit <- p0 - z_alpha * se0
    actual_power <- pnorm(crit, mean = p1, sd = se1)
  }
  
  # 返回结果列表
  return(list(
    required_sample_size = n,
    actual_power = round(actual_power, 4),
    null_proportion = p0,
    alternative_proportion = p1,
    alpha = alpha
  ))
}

代入你的硬币问题参数计算 #

# 原假设概率
p0 <- 0.5
# 备择假设概率（以0.51为例，0.49结果对称）
p1 <- 0.51

# 计算双侧检验、功效80%、α=0.05所需的样本量
result <- power.prop.one.sample(p0 = p0, p1 = p1, alpha = 0.05, power = 0.8, alternative = "two.sided")

cat("需要抛硬币的次数：", result$required_sample_size, "\n")
cat("实际能达到的功效：", result$actual_power, "\n")

运行后你会得到需要的抛投次数（大概是~19200次左右，因为p0和p1差异很小，需要大样本）。

方法2：模拟法计算功效（更准确，无近似假设） #

如果担心正态近似的偏差，可以用模拟的方式，通过多次重复实验来估算功效：

# 模拟单样本比例检验的功效
simulate_power <- function(p0, p1, n, alpha = 0.05, alternative = "two.sided", n_sim = 10000) {
  # 重复n_sim次实验
  p_values <- replicate(n_sim, {
    # 生成n次抛投的结果，真实概率为p1
    heads <- rbinom(1, size = n, prob = p1)
    # 做单样本二项检验
    test_result <- binom.test(heads, n = n, p = p0, alternative = alternative)
    test_result$p.value
  })
  
  # 功效是拒绝原假设的比例（p值<α）
  power <- mean(p_values < alpha)
  return(round(power, 4))
}

# 用方法1得到的样本量验证功效
simulated_power <- simulate_power(p0 = 0.5, p1 = 0.51, n = result$required_sample_size, n_sim = 10000)
cat("模拟得到的功效：", simulated_power, "\n")

模拟法的结果会和正态近似的结果非常接近，因为这里样本量很大，二项分布已经很接近正态了。

方法3：用扩展包pwr简化计算（推荐，若允许安装包） #

如果你不介意安装扩展包，pwr包提供了专门的单样本比例检验功效函数pwr.p.test()：

# 先安装包（仅第一次需要）
# install.packages("pwr")
library(pwr)

# 计算效应量（反正弦转换的差异）
effect_size <- ES.h(p1, p0)

# 计算样本量
pwr_result <- pwr.p.test(h = effect_size, sig.level = 0.05, power = 0.8, alternative = "two.sided")
print(pwr_result)

这个函数的结果和我们手动实现的正态近似方法一致，因为它底层也是用的正态近似。

内容的提问来源于stack exchange，提问作者Andrew Davidson