3

我设计了 3000 个实验，因此在一个实验中有 4 组（治疗），每组有 50 个人（受试者）。对于每个实验，我都会做一个标准的单向方差分析，并证明它们的 p.value 在零假设下是否具有单一概率函数，但是 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")
4

1 回答 1

9

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

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

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

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

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
于 2013-06-29T16:03:14.240 回答