Kolmogorov Smirnov Test in R

Question

한 실험에서 4 그룹 (치료)이되도록 각 그룹에 50 명의 개체 (피험자)가 있도록 3000 회 실험을 설계했습니다. 각각의 실험에서 표준 일원 분산 분석을 수행하고 그들의 p.values가 null 가설 하에서 단일 확률 함수를 갖고 있는지 증명하지만, 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 값이 .34 인 줄을 주석 처리했습니다. 그것은 알려지지 않은 씨앗 이었기 때문에 재현성을 위해 set.seed(1)을 실행하고 다시 실행했습니다. 이번에는 0.98의 P 값을 얻었다.

이것이 차이를 만드는 이유는 PRNGs의 전문가는 아니지만 괜찮은 발전기는 연속적인 추첨이 모든 실제적인 목적을 위해 통계적으로 독립적임을 보장합니다. 가장 좋은 것들은 더 큰 시차를 위해 동일하게 할 것입니다. 예를 들어 메르 센 트위스터는 R의 기본 PRNG가 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