进行简单的模拟以确认统计检验的功效？

机器算法验证 r 可能性 t检验 模拟 统计能力
2022-04-12 04:57:19

我试图做一个简单的模拟来验证我在 R 中使用包 pwr 找到的一些功率计算可能看起来很奇怪（https://cran.r-project.org/web/packages/pwr/pwr. .pdf )。

我所做的是使用 N 为 25 和样本为 0.5 的单样本 t 检验的最简单示例：

          pwr.t.test(n = 25,d = .5,type= "one.sample",sig.level = .05)

结果表明，在这种情况下应该有大约 67% 的功效： 单样本 t 检验功效计算

          n = 25
          d = 0.5
  sig.level = 0.05
      power = 0.6697077
alternative = two.sided

然后，我做了一个小型模拟，从平均值为 107.5、标准差为 15（得到 ad=.5）的分布中抽样 100,000 次：

                  g <- matrix(NA,nrow = 100000,ncol=1)
                  for(i in 1:100000){
                  x <- rnorm(n = 25,mean = 107.5,sd=15)
                  g[i] = mean(x)}

最后，我计算了这些样本超过抽样分布的 95% CI 的次数：

        upper = 100 + (1.96*(15/sqrt(25)))
        lower = 100 - (1.96*(15/sqrt(25)))
        Outcomes <- ifelse((g > upper | g < lower),"Reject Null      Hypothesis","Fail to Reject Null Hypothesis")
        table(Outcomes)/100000

        Outcomes
        Fail to Reject Null Hypothesis         Reject Null Hypothesis 
                   0.29726                        0.70274

如您所见，pwr 包的结果与我的模拟结果之间存在微小差异（~3%）。它看起来很小但不平凡。我的逻辑中是否有错误？

1个回答

您在模拟中使用的检验不是 t 检验，因为您使用的是 1.96 而不是 t 值，并且您使用的是真实标准差而不是样本标准差。您的模拟近似于相应“z 检验”的功效：

> pwr.norm.test(d = 0.5, n = 25, power = NULL)

     Mean power calculation for normal distribution with known variance 

              d = 0.5
              n = 25
      sig.level = 0.05
          power = 0.705418
    alternative = two.sided

以下代码显示了如何验证 给出的功率值pwr.t.test。

n <- 25      # sample size
mu <- 107.5  # true mean
sigma <- 15  # true SD
mu0 <- 100   # mean under the null hypothesis

reps <- 100000  # number of simulations

## p-value approach:

pvalues <- numeric(reps)

set.seed(1)

for (i in 1:reps) {
  x <- rnorm(n, mu, sigma)
  t.stat <- (mean(x) - mu0)/(sd(x)/sqrt(n))
  pvalues[i] <- 2*(1 - pt(abs(t.stat), n-1))
  # alternatively: pvalues[i] <- t.test(x, mu = mu0)$p.value
}

> mean(pvalues < 0.05)
[1] 0.66907

## Confidence interval approach:

outsideCI <- numeric(reps) # 1 if mu0 not in 95% CI, otherwise 0

set.seed(2)

for (i in 1:reps) {
  x <- rnorm(n, mu, sigma)
  CI.lower <- mean(x) - qt(0.975, n-1)*sd(x)/sqrt(n)
  CI.upper <- mean(x) + qt(0.975, n-1)*sd(x)/sqrt(n)
  outsideCI[i] <- ifelse(mu0 < CI.lower | mu0 > CI.upper, 1, 0)
}

> mean(outsideCI)
[1] 0.66893
其它你可能感兴趣的问题
上一篇为什么神经网络中不经常使用高斯激活函数？ 下一篇Gamma 形状和尺度参数的先验