3

私は 3000 の実験を設計したので、1 つの実験には 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")
4

1 に答える 1

9

コード内のランダム シードを変更する行をコメント アウトしたところ、P 値は .34 になりました。それは未知の種だったので、再現性のために、set.seed(1)もう一度実行しました。今回は P 値 0.98 を得ました。

なぜこれが違いを生むのかについては、私は PRNG の専門家ではありませんが、適切なジェネレーターは、すべての実用的な目的で、連続する描画が統計的に独立していることを保証します。最良のものは、より大きなラグに対しても同じことを保証します。たとえば、R のデフォルト PRNG である Mersenne Twister は、最大 623 (IIRC) のラグに対してそれを保証します。実際、シードに干渉すると、抽選の統計的特性が損なわれる可能性があります。

あなたのコードはまた、非常に非効率的な方法で物事を行っています。実験のリストを作成し、実験ごとに 1 つの項目を追加しています。各実験内で、マトリックスも作成し、観測ごとに行を追加します。次に、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 に答える