两个样本 t 检验的功效

机器算法验证 统计能力
2022-03-24 23:11:44

我试图了解两个独立样本 t 检验的功率计算（不假设方差相等，所以我使用了 Satterthwaite）。

这是我发现有助于理解该过程的图表：

在此处输入图像描述

所以我假设给定以下关于两个总体的信息并给出样本量：

mu1<-5
mu2<-6
sd1<-3
sd2<-2
n1<-20
n2<-20

我可以计算与具有 0.05 上尾概率相关的空值下的临界值：

df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )
CV<- qt(0.95,df) #equals 1.730018

然后计算备择假设（在这种情况下，我学到的是“非中心 t 分布”）。我使用非中心分布和上面找到的临界值计算了上图中的 beta。这是R中的完整脚本：

#under alternative
mu1<-5
mu2<-6
sd1<-3
sd2<-2
n1<-20
n2<-20


#Under null
Sp<-sqrt(((n1-1)*sd1^2+(n2-1)*sd2^2)/(n1+n2-2))
df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )
CV<- qt(0.95,df)


#under alternative
diff<-mu1-mu2
t<-(diff)/sqrt((sd1^2/n1)+ (sd2^2/n2))
ncp<-(diff/sqrt((sd1^2/n1)+(sd2^2/n2)))


#power
1-pt(t, df, ncp)

这给出了 0.4935132 的功率值。

这是正确的方法吗？我发现如果我使用其他功率计算软件（比如 SAS，我认为我已经设置了与下面的问题等效的软件）我会得到另一个答案（来自 SAS 是 0.33）。

SAS 代码：

proc power;
      twosamplemeans test=diff_satt
         meandiff = 1
         groupstddevs = 3 | 2
         groupweights = (1 1)
         ntotal = 40
         power = .
        sides=1;
   run;

最终，我希望获得一种理解，使我能够查看更复杂程序的模拟。

编辑：我发现了我的错误。本来应该

1-pt(CV, df, ncp) 不是 1-pt(t, df, ncp)

2个回答

您已经接近了，但需要进行一些小的更改：

  • 均值的真实差异通常被视为μ2−μ1，而不是相反。
  • G*电源使用n1+n2−2作为自由度t- 在这种情况下的分布（不同的方差，相同的组大小），按照 Cohen 的建议，如此处所述
  • SAS 可能会使用韦尔奇公式或萨特斯韦特公式来计算 df 给定的不等方差（在您引用的此 pdf中找到） - 结果中只有 2 个有效数字无法分辨（见下文）

如n1, n2, mu1, mu2, sd1, sd2您的问题中所定义：

> alpha   <- 0.05
> dfGP    <- n1+n2 - 2                     # degrees of freedom (used by G*Power)
> cvGP    <- qt(1-alpha, dfGP)             # crit. value for one-sided test (under the null)
> muDiff  <- mu2-mu1                       # true difference in means
> sigDiff <- sqrt((sd1^2/n1) + (sd2^2/n2)) # true SD for difference in empirical means
> ncp     <- muDiff / sigDiff              # noncentrality parameter (under alternative)
> 1-pt(cvGP, dfGP, ncp)                    # power
[1] 0.3348385

这与G*Power的结果相匹配，这是一个很好的解决这些问题的程序。它还显示 df、临界值和 ncp，因此您可以单独检查所有这些计算。

在此处输入图像描述

编辑：使用 Satterthwaite 的公式或 Welch 的公式变化不大（仍然是 0.33*）：

# Satterthwaite's formula
> var1  <- sd1^2
> var2  <- sd2^2
> num   <- (var1/n1 + var2/n2)^2
> denST <- var1^2/((n1-1)*n1^2) + var2^2/((n2-1)*n2^2)
> (dfST <- num/denST)
[1] 33.10309

> cvST <- qt(1-alpha, dfST)
> 1-pt(cvST, dfST, ncp)
[1] 0.3336495

# Welch's formula
> denW <- var1^2/((n1+1)*n1^2) + var2^2/((n2+1)*n2^2)
> (dfW <- (num/denW) - 2)
[1] 34.58763

> cvW   <- qt(1-alpha, dfW)
> 1-pt(cvW, dfW, ncp)
[1] 0.3340453

（注意我把一些变量名稍微改成了t, df,diff也是内置函数的名字，还要注意你的代码的分子df是错误的，它有一个错位^2，一个^2太多了，应该是((sd1^2/n1) + (sd2^2/n2))^2）

如果您主要对计算功率感兴趣（而不是通过手工学习）并且您已经在使用 R，那么请查看pwr包和pwr.t.testorpwr.t2n.test函数。（即使您手动学习，这些也可以很好地验证您的结果）。

其它你可能感兴趣的问题
上一篇超越 Fisher 内核 下一篇在二项式数据上使用 R 的良好拟合时的 NaN p 值