A small simulation study on tests for normality in R
I'm doing a little simulation research to judge how good the two tests of normality are. My plan is to generate many common distribution patterns of not too many cases and determine how often each test rejects the null hypothesis of normality.
The (incomplete) code I have so far is
library(nortest)
y<-replicate(10000,{
x<-rnorm(50)
ad.test(x)$p.value
ks.test(x,y=pnorm)$p.value
}
)
Now I would like to calculate the proportion of these p-values that are less than 0.05 for each test. Could you please tell me how can I do this? I apologize if this is a newbie question, but I'm new to R. myself.
Thank.
source to share
library(nortest)
nsim <- 10000
nx <- 50
set.seed(101)
y <- replicate(nsim,{
x<-rnorm(nx)
c(ad=ad.test(x)$p.value,
ks=ks.test(x,y=pnorm)$p.value)
}
)
apply(y<0.05,MARGIN=1,mean)
## ad ks
## 0.0534 0.0480
Usage MARGIN=1
says apply
to average across rows rather than columns - which is reasonable given the ordering that the inline simplification creates replicate()
.
For examples of this type, the I error rates of any standard tests will be extremely close to their nominal value (0.05 in this example).
source to share
Your code doesn't output p values. You can do something like this:
rep_test <- function(reps=10000) {
p_ks <- rep(NA, reps)
p_ad <- rep(NA, reps)
for (i in 1:reps) {
x <- rnorm(50)
p_ks[i] <- ks.test(x, y=pnorm)$p.value
p_ad[i] <- ad.test(x)$p.value
}
return(data.frame(cbind(p_ks, p_ad)))
}
sum(test$p_ks<.05)/10000
sum(test$p_ad<.05)/10000
source to share