R 中的 I 类错误模拟

Question

我正在尝试使用 Monte Carlo 模拟计算双变量正态数据相关性检验的第 i 类错误率和功效。

但是我得到了 I 类错误和功率的意外值。（I 类错误为 0.864）

我想知道我是否做错了什么。谁能帮帮我？

set.seed(160230)
library("mvtnorm", lib.loc="~/R/win-library/3.4")
sigma= matrix(c(1,0.8,0.8,1),2,2) 
mu <- c(0,0)
#bivariate normal data
sim=replicate(n=1000 , rmvnorm(10,mean=mu , sigma = sigma))

pval1=c()

for(i in 1:1000)
{
 pval1[i]=cor.test(sim[,1,i],sim[,2,i],method = c("pearson"))$p.value

  }
#type1 error rate
mean(pval1<0.05)

#power
mean(pval3>0.05)

Answer 1

您的代码没问题，但您的模拟设置有误。

在您的代码中，您

模拟具有强相关性的双变量数据，rho=0.8。
检验假设 H0: rho=0。

因此，您是在备择假设下模拟数据，这就是您得到 0.864 结果的原因。这本质上是你对那个特定选择的权力。您可以改为执行以下操作：

先模拟原假设下的数据

sigma <- matrix(c(1,0,0,1),2,2) 
mu <- c(0,0)
#bivariate normal data under H0
sim <- replicate(n=1000, rmvnorm(10, mean=mu, sigma = sigma))

# Test the actual level under H0

result <- sapply(1:1000, function(i) { 
    cor.test(sim[,1,i],sim[,2,i],method = c("pearson"))$p.value})

mean(result < 0.05)

给出的值约为 0.05。在替代方案下，您可以使用相关系数为 0.8（或其他一些数字）的代码。您可以使用以下代码对此进行概括，以轻松获得多个相关性的功效。

rho <- seq(0, .9, .1)
pwr <- sapply(rho, function(r) {
    sigma <- matrix(c(1,r,r,1),2,2) 
    mu <- c(0,0)
    #bivariate normal data 
    sim <- replicate(n=1000, rmvnorm(10, mean=mu, sigma = sigma))

    # Test the actual level
    result <- sapply(1:1000, function(i) { 
        cor.test(sim[,1,i],sim[,2,i],method = c("pearson"))$p.value})

    mean(result < 0.05)
})

然后您可以通过绘制关系来查看相关性对功效的影响

plot(rho, pwr, type="l", xlab=expression(rho), ylab="Power")

R 中的 I 类错误模拟

type I error simulation in R

r

simulation

montecarlo