在R中使用Kolmogorov Smirnov測試

我設計了3000個實驗，因此在一個實驗中有4組（治療），每組有50個個體（受試者）。對於每個實驗，我做一個標準的單向方差分析，並證明它們的p.values在無效假設下是否具有單一概率函數，但ks.test否定了這個假設，我不明白爲什麼？在R中使用Kolmogorov Smirnov測試

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")

來源

2013-06-29 Klaus

短語「多重比較校正」對您來說意味着什麼？ – zwol

不，不是真的，我讀了它在維基百科，不能看到我的模擬的相關性。我設計了獨立實驗，並且只對一個假設進行了每個實驗的測試，而不是更多。 – Klaus

@Zack這是一個模擬研究。 OP正在計算每個實驗的單個P值，但重複該過程多次以檢查P值統計的屬性。 –

我註釋掉了代碼中改變隨機種子的行，並得到了0.34的P值。這是一個未知的種子，所以爲了再現性，我做了set.seed(1)並再次運行它。這一次，我得到了0.98的P值。至於爲什麼這會有所作爲，我並不是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-29 16:03:14

如果你只使用set.seed（1）爲每個rn代，這改變了實驗的設計。關於你對代碼效率的評論，我知道代碼在這種情況下效率不高，但我使用這個框架來處理複雜的事情，並通過複製粘貼上面的特定設計來嘗試。我的問題只針對r.n.代。 – Klaus

我在循環中不使用'set.seed（1）'。我在開始時設置種子ONCE，然後運行你的代碼（在註釋掉重置種子的行之後）。 –

好吧，我誤解了set.seed的含義。我想我必須在每次電話會議之前設置它。生成函數以獲得可重複的數據。第二件事是，如果我使用rnorm來設置種子，但是如果我使用rmvnorm，那麼我得到了「良好的價值」。所以我認爲設置一顆種子就能解決問題。 – Klaus

在R中使用Kolmogorov Smirnov測試

回答

相關問題