在R中进行正态性检验的小型模拟研究

Question

我正在进行一项小型模拟研究，以评估两个正态性检验的好坏。我的计划是生成许多不太多的正态分布样本，并确定每个测试有多少次拒绝正态性的零假设。

我目前有的（不完整）代码如下：

  library(nortest)
  y<-replicate(10000,{
     x<-rnorm(50)
     ad.test(x)$p.value
     ks.test(x,y=pnorm)$p.value
   }
   )

现在我想计算每个测试中小于0.05的p值的比例。请问我该如何做？如果这是一个新手问题，我很抱歉，因为我自己对R也很陌生。

谢谢。

- JohnK

Answer 1

如果你分别运行每个测试，那么可以简单地计算存储在y中小于0.05的值的数量。

y<-replicate(1000,{
     x<-rnorm(50)
     ks.test(x,y=pnorm)$p.value})
length(which(y<0.05))

- darwin

---

谢谢，但我希望每个测试都在同一个样本上进行。 - JohnK

Answer 2

 library(nortest)
 nsim <- 10000
 nx <- 50

 set.seed(101)
 y <- replicate(nsim,{
    x<-rnorm(nx)
    c(ad=ad.test(x)$p.value,
      ks=ks.test(x,y=pnorm)$p.value)
  }
 )
 apply(y<0.05,MARGIN=1,mean)
 ##     ad     ks 
 ## 0.0534 0.0480

使用MARGIN=1告诉apply沿行取平均值，而不是列--这是合理的，考虑到replicate()内置简化产生的排序。

对于这种类型的示例，任何标准测试的类型I错误率将非常接近其名义值（在此示例中为0.05）。

- Ben Bolker

---

你可以通过一次性将所有的正常值选取并放入矩阵中来加快速度... - Ben Bolker

谢谢，这就可以了。但我想问一下，为什么你在 apply 函数的第二个参数中放了 1？ - JohnK

Answer 3

你的代码没有输出p值。你可以像这样做：

rep_test <- function(reps=10000) {

  p_ks <- rep(NA, reps)
  p_ad <- rep(NA, reps)

  for (i in 1:reps) {
    x <- rnorm(50)
    p_ks[i] <- ks.test(x, y=pnorm)$p.value
    p_ad[i] <- ad.test(x)$p.value
  }

  return(data.frame(cbind(p_ks, p_ad)))
}

sum(test$p_ks<.05)/10000
sum(test$p_ad<.05)/10000

- eamcvey