R中最小效应假设(非中心f分布)的功率计算

Question

在工作中，我被要求计算Wald测试统计量的统计能力。我们正在针对1.54%的解释方差的零假设进行测试，这就是为什么我们比较的F分布是非中心的。

现在，我被告知首先在显著性水平的零假设下计算临界F值(我们使用.05)。我在R中这样做：

`groups <- 2
N <- 4440
PV <- 0.0154

min.eff.hyp <- function(N, groups=2, PV=.01, alpha=.05){
  df1 <- groups - 1 # numerator degrees of freedom
  df2 <- N-groups # denominator degrees of freedom
  ncp <- ((N-1)*PV)/(1-PV) # non-centrality parameter
  g <- (df1+ncp)^2/(df1+2*ncp) # adjusted df1
  k <- (df1+ncp)/df1 # to adjust for non-centrality
  crit.central <- qf(1-alpha,g,df2)
  crit.val <- crit.central * k
  return(paste('Kritischer F-Wert (PV=',PV,'): ',round(crit.val,3), sep = ''))
}

min.eff.hyp(N,groups,PV)`

然后我应该使用这个关键的F值来确定在另一种假设下的alpha‘。然后，这个alpha‘可以用来计算统计功率。

不幸的是，我不知道如何获得alpha‘。我再次尝试使用R：pf(crit.val,g,df2,ncp)

我认为，这应该计算我的数据的概率，假设另一种假设为真。但说实话，我现在很迷茫。我不知道如何在功率计算中实现非中心性，而且不知何故，我找不到以前遇到过同样问题的人，并且实际上找到了解决方案。

当我针对最小效应假设进行测试，从而与非中心分布进行比较时，如何计算测试的统计能力？

非常感谢您抽出时间来帮助我！

你好，玛丽亚

Answer 1

我们可以计算F-检验统计量的统计功率，只需：# of groups，# of percentage N，显著性水平alpha和解释方差PV的百分比。

一种方法是使用带有估计的非中心性参数NCP的非中心F分布。另一种方法，如您的示例所示，是使用调整后的中心F分布(Murphy et al. 2009)。

仿真结果表明，第一种方法对power的估计精度更高。

>        theoretical  simulation  method_1   method_2
> Power    0.4631185      0.4647  0.460971  0.5521764
> F mean    4.518145    4.529506  4.497807          -
> PV       0.0089153   0.0089509         -          -
> NCP            3.5           -  3.479744          -

非中心F分布的模拟：

## Simulation to get PV and power
groups <- 2
N <- 500
alpha <- 0.05
ncp <- 3.5 #Non-Centrality Parameter

#simulated Sum of Squares for Model and Error
set.seed(1)
df1 <- groups - 1 #main effect degrees of freedom 
df2 <- N - groups #error degrees of freedom
SSM = SSE = rep(NA,10000)
for (i in 1:10000)
  {
    SSM[i] <- sum(rnorm(df1, mean=sqrt(ncp/df1))^2)
    SSE[i] <- sum(rnorm(df2)^2)
  }

#mean F
F <- SSM/df1 / (SSE/df2)
mean(F) #estimated
> 4.529506
(df1+ncp) / df1 / ((df2-2)/df2) #theoretical
> 4.518145

#PV: Percentage of explained Variance
PV <- SSM / (SSM + SSE)
mean(PV)
> 0.008950899
(ncp + df1) / (ncp + df1 + df2) #theoretical
> 0.008915272

#statistical Power
F_crit <- qf(1-alpha, df1, df2) #critical F value
sum(F > F_crit) / length(F) #estimated
> 0.4647
1 - pf(F_crit, df1, df2, ncp) #theoretical
> 0.4631185

方法1(非中心F分布)：

groups <- 2
N <- 500
alpha <- 0.05
PV <- 0.008950899
ncp <- NA

df1 <- groups - 1
df2 <- N - groups

#mean F
PV/df1 / ((1-PV) / df2)
> 4.497807

#estimated Non-Centrality Parameter
ncp <- (df2 - 2) * PV / (1-PV) - df1
> 3.479744

#estimated statistical power
F_crit <- qf(1-alpha, df1, df2) #critical F value
1 - pf(F_crit, df1, df2, ncp=ncp)
> 0.460971

方法2(调整后的中心F分布)：

lambda <- df2 * PV / (1 - PV) # another definition of non-centrality parameter
g <- (df1 + lambda)^2 / (df1 + 2 * lambda) # adjusted df1
k <- (df1 + lambda) / df1 # to adjust for non-centrality

#estimated statistical power
F_crit <- qf(1-alpha, df1, df2)
1 - pf(F_crit/k, g, df2)
> 0.5521764

方法1背后的代数(估计非中心性参数)：

​

。

​

​

其中最后一个等价物来自wikipedia。

​

。