r - 在 R 中使用 Kolmogorov Smirnov 检验

Question

我设计了 3000 个实验，因此在一个实验中有 4 组(治疗组)，每组有 50 个人(受试者)。对于每个实验，我都会做一个标准的单向方差分析，并证明它们的 p.values 在零假设下是否具有单一概率函数，但 ks.test 拒绝了这个假设，我不明白为什么？

subject<-50
treatment<-4
experiment<-list()
R<-3000
seed<-split(1:(R*subject),1:R)
for(i in 1:R){
  e<-c()
  for(j in 1:subject){
    set.seed(seed[[i]][j]) 
    e<-c(e,rmvnorm(mean=rep(0,treatment),sigma=diag(3,4),n=1,method="chol"))
   }
  experiment<-c(experiment,list(matrix(e,subject,treatment,byrow=T)))
 }

 p.values<-c()
for(e in experiment){
  d<-data.frame(response=c(e),treatment=factor(rep(1:treatment,each=subject)))
  p.values<-c(p.values,anova(lm(response~treatment,d))[1,"Pr(>F)"])
 }

 ks.test(p.values, punif,alternative = "two.sided")

Answer 1

我注释掉了代码中更改随机种子的行，得到的 P 值为 0.34。那是一个未知的种子，所以为了重现性，我做了 set.seed(1) 并再次运行它。这一次，我的 P 值为 0.98。

至于为什么这会产生影响，我不是 PRNG 方面的专家，但任何体面的生成器都可以确保连续抽签在所有实际用途中都是统计独立的。最好的将确保更大的延迟相同，例如，R 的默认 PRNG 的 Mersenne Twister 保证它的延迟高达 623 (IIRC)。事实上，干预种子很可能会损害抽签的统计特性。

您的代码也在以一种非常低效的方式做事。您正在为实验创建一个列表，并为每个实验添加一个项目。 在每个实验中，您还可以创建一个矩阵，并为每个观察添加一行。然后你对 P 值做一些非常相似的事情。我看看能不能解决。

这就是我要替换您的代码的方式。严格来说，我可以通过避免使用公式、创建裸模型矩阵并直接调用 lm.fit 来使其更加紧凑。但这意味着必须手动编写 ANOVA 测试代码，而不是简单地调用 anova，这比它的值(value)更麻烦。

set.seed(1) # or any other number you like

x <- factor(rep(seq_len(treatment), each=subject))
p.values <- sapply(seq_len(R), function(r) {
    y <- rnorm(subject * treatment, s=3)
    anova(lm(y ~ x))[1,"Pr(>F)"]
})
ks.test(p.values, punif,alternative = "two.sided")


        One-sample Kolmogorov-Smirnov test

data:  p.values
D = 0.0121, p-value = 0.772
alternative hypothesis: two-sided